From c39e9ab3fc4753ac3adaf271f3b34125f10f0366 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 27 Sep 2024 17:56:03 +0200 Subject: [PATCH] Trying to implement floating potential These things are never easy. --- .gitignore | 1 + plasmaExpansion.f90 | 47 ++++++++++++++++++++++---------------- scripts_python/plotCumF.py | 3 ++- scripts_python/plotMom.py | 11 +++++---- scripts_python/readPhi.py | 3 ++- 5 files changed, 38 insertions(+), 27 deletions(-) diff --git a/.gitignore b/.gitignore index 9d36d42..eb6b4d9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ plasmaExpansion +*.pyc *.csv *.mod diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90 index b58eaa5..4e906d5 100644 --- a/plasmaExpansion.f90 +++ b/plasmaExpansion.f90 @@ -40,7 +40,7 @@ module output integer, parameter:: dataF_id = 30 integer, parameter:: dataPhi_id = 40 integer, parameter:: dataCum_id = 50 - character(len=7), parameter:: formatFloat = 'ES0.4e3' + character(len=7), parameter:: formatFloat = 'ES0.6e3' character(len=3), parameter:: formatSep = '","' character(len=7):: formatTime @@ -103,7 +103,7 @@ module output end subroutine writeOutputF - subroutine writeOutputPhi(t, dt, nr, r, phi, n_e) + subroutine writeOutputPhi(t, dt, nr, r, phi, E, n_e) use constantParameters, only: eV_to_K use referenceValues, only: L_ref, phi_ref, t_ref, n_ref @@ -112,6 +112,7 @@ module output real(dp), intent(in):: dt real(dp), intent(in):: r(1:nr) real(dp), intent(in):: phi(1:nr) + real(dp), intent(in):: E(1:nr) real(dp), intent(in):: n_e(1:nr) character(:), allocatable:: filename integer:: i @@ -124,11 +125,12 @@ module output open(unit=dataPhi_id, file=pathOutput//filename) write(dataPhi_id, '(A)') "t (s)" write(dataPhi_id, '('//formatFloat//')') t*dt*t_ref - write(dataPhi_id, '(A,2('//formatSep//',A))') "r (m)","phi (V)","n_e (m^-3)" + write(dataPhi_id, '(A,3('//formatSep//',A))') "r (m)","phi (V)","E (V m^-1)","n_e (m^-3)" do i = 1, nr - write(dataPhi_id, '('//formatFloat//',2('//formatSep //','//formatFloat//'))') & + write(dataPhi_id, '('//formatFloat//',3('//formatSep //','//formatFloat//'))') & r(i)*L_ref, & phi(i)*phi_ref, & + E(i)*phi_ref/L_ref, & n_e(i)*n_ref end do @@ -328,7 +330,7 @@ program plasmaExpansion real(dp), allocatable, dimension(:):: phi, phi_old, E real(dp):: phiConv - real(dp):: phi0 + real(dp):: phi0, phiF, JF integer:: k real(dp), allocatable, dimension(:):: fCum_i @@ -352,17 +354,17 @@ program plasmaExpansion ! Set input parameters (remember these have to be in non-dimensional units) Temp0 = 60.0_dp * eV_to_K / Temp_ref - TempF = 10.0_dp * eV_to_K / Temp_ref + TempF = 60.0_dp * eV_to_K / Temp_ref n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref c_s = sqrt(T_to_Z(Temp0) * gam * Temp0) u_bc0 = sqrt(Temp0) u_bcF = sqrt(TempF) n_bc0 = n_ecr / T_to_Z(Temp0) - n_bcF = n_ecr*1.0e-1 / T_to_Z(Temp0) + n_bcF = n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0) ! Set domain boundaries (non-dimensional units) r0 = 200.0e-6_dp / L_ref - rf = 1.0e-2_dp / L_ref + rf = 6.0e-3_dp / L_ref dr = 1.0e3_dp nr = nint((rf - r0) / dr) + 1 dr = (rf - r0) / float(nr-1) @@ -470,8 +472,9 @@ program plasmaExpansion ! Set boundary values ! phi0 = 0.0_dp / phi_ref ! Dirichlet ! phi(1) = phi0 ! Dirichlet - phi0 = phi(1) ! Neumann - phi(nr) = 0.0_dp ! Dirichlet + phi0 = phi(1) ! Neumann + phiF = 0.0_dp ! Dirichlet + JF = 0.0_dp ! Current at the edge for floating potential allocate(f0(j0:nv)) f0 = 0.0_dp @@ -481,7 +484,7 @@ program plasmaExpansion t = 0 call writeOutputRef() call writeOutputF(t, dt, nr, r, nv, v, f_i_old) - call writeOutputPhi(t, dt, nr, r, phi, n_e) + call writeOutputPhi(t, dt, nr, r, phi, E, n_e) call writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Zave) ! Main loop @@ -528,7 +531,7 @@ program plasmaExpansion end if n_i(i) = sum(f_i(i,:)*dv) - if (n_i(i) > 0.0_dp) then + if (n_i(i) > 1.0e-10_dp) then u_i(i) = sum(v(:) *f_i(i,:)*dv) / n_i(i) E_i(i) = sum(v(:)**2*f_i(i,:)*dv) / n_i(i) T_i(i) = 2.0_dp*E_i(i) - 2.0_dp*u_i(i)**2 @@ -544,7 +547,7 @@ program plasmaExpansion !$omp end parallel do ! Solve Poission (maximum number of iterations, break if convergence is reached before) - do k = 1, 200 + do k = 1, 100 ! Store previous value phi_old = phi @@ -566,20 +569,25 @@ program plasmaExpansion b = -(Zave*n_i - n_e) ! Apply boundary conditions ! b(1) = phi0 ! Dirichlet - b(nr) = 0.0_dp ! Dirichlet + b(nr) = phiF ! Dirichlet ! Calculate residual !$omp parallel workshare Res = -(MATMUL(A, phi_old) - b) !$omp end parallel workshare - ! Res(1) = 0.0_dp ! Dirichlet - Res(nr) = 0.0_dp ! Dirichlet ! Iterate system call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info) phi = phi_old + Res - phi0 = phi(1) ! Neumann - phi(nr-10:nr) = phi(nr-11) + phi(1) = phi(2) ! Neumann + phi0 = phi(1) ! Neumann + ! phi(nr) = phi(nr-5) ! Dirichlet quasi-floating + + ! Calculate updated current + JF = Zave(nr)*n_i(nr)*u_i(nr) - n_e(nr)*sqrt(qe*T_i(1)*Temp_ref/9.1e-31_dp)/u_ref + + ! Iterate floating potential + phiF = phi(nr-5)! + 1.0e-2_dp*JF ! Check if the solution has converged phiConv = maxval(abs(Res),1) @@ -602,7 +610,6 @@ program plasmaExpansion E(nr) = - (phi(nr) - phi(nr-1)) / dr ! Dirichlet ! E(nr) = 0.0_dp ! Neumann ! Trick to avoid problems at the sheath - E(nr-5:nr) = 0.0_dp ! Update intermediate f f_i_old = f_i @@ -644,7 +651,7 @@ program plasmaExpansion ! Write output if (mod(t,everyOutput) == 0 .or. t == nt) then call writeOutputF(t, dt, nr, r, nv, v, f_i_old) - call writeOutputPhi(t, dt, nr, r, phi, n_e) + call writeOutputPhi(t, dt, nr, r, phi, E, n_e) call writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Zave) call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i) diff --git a/scripts_python/plotCumF.py b/scripts_python/plotCumF.py index 4bab15b..cc8a0a3 100644 --- a/scripts_python/plotCumF.py +++ b/scripts_python/plotCumF.py @@ -8,8 +8,9 @@ from scipy.constants import e, k m_i = 1.9712e-25 -paths = ['../quasiNeutral_fullAblation/','../Poisson_fullAblation/'] +# paths = ['../quasiNeutral_fullAblation/','../Poisson_fullAblation/'] # paths = ['../quasiNeutral_partialAblation/','../Poisson_partialAblation/'] +paths = ['../2024-09-27_17.41.21/'] for path in paths: filesCum_i = sorted(glob.glob(path+'time_*_fCum_i.csv')) start = 0 diff --git a/scripts_python/plotMom.py b/scripts_python/plotMom.py index 5816516..ce5a45c 100644 --- a/scripts_python/plotMom.py +++ b/scripts_python/plotMom.py @@ -8,23 +8,24 @@ from scipy.constants import e, k # paths = ['../quasiNeutral_fullAblation/','../2024-09-26_11.48.04/'] # paths = ['../2024-09-26_12.47.11/'] -paths = ['../quasiNeutral_partialAblation/','../2024-09-26_13.58.24/'] +# paths = ['../quasiNeutral_partialAblation/','../2024-09-26_13.58.24/'] # path = '../quasiNeutral_fullAblation/' # path = '../quasiNeutral_partialAblatio/' - +paths = ['../2024-09-27_17.41.21/'] for path in paths: filesPhi = sorted(glob.glob(path+'time_*_phi.csv')) filesMom_i = sorted(glob.glob(path+'time_*_mom_i.csv')) start = 0 end = len(filesMom_i) - every = 100 + every = 50 fig, ax = plt.subplots(4, sharex='all') for fileMom_i, filePhi in zip(filesMom_i[start:end+1:every], filesPhi[start:end+1:every]): - time, r, phi, n_e = readPhi.read(filePhi) + time, r, phi, E, n_e = readPhi.read(filePhi) time, r, n_i, u_i, T_i, Zave = readMom.read(fileMom_i) - ax[0].plot(r, phi, label='t = {:.3f} ns'.format(time*1e9)) + # ax[0].plot(r, phi, label='t = {:.3f} ns'.format(time*1e9)) + ax[0].plot(r, E, label='t = {:.3f} ns'.format(time*1e9)) ax[1].set_yscale('log') ax[1].set_ylim([1e14,2e25]) ax[1].plot(r, Zave*n_i) diff --git a/scripts_python/readPhi.py b/scripts_python/readPhi.py index 9be2ced..e0b0fd5 100644 --- a/scripts_python/readPhi.py +++ b/scripts_python/readPhi.py @@ -8,8 +8,9 @@ def read(filename): df = pandas.read_csv(filename,skiprows=2) x = df['r (m)'].to_numpy() phi = df['phi (V)'].to_numpy() + E = df['E (V m^-1)'].to_numpy() n_e = df['n_e (m^-3)'].to_numpy() - return time, x, phi, n_e + return time, x, phi, E, n_e