Trying to implement floating potential

These things are never easy.

.gitignore

@@ -1,3 +1,4 @@
 plasmaExpansion
+*.pyc
 *.csv
 *.mod

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
+  phi0 = phi(1) ! Neumann
-  phi(nr) = 0.0_dp ! Dirichlet
+  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(1)  = phi(2)    ! Neumann
-      phi(nr-10:nr) = phi(nr-11)
+      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)
 

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

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)

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