Potential at r0 is different from 0 now. Posibility of quasi-neutrality properly accounting for the electric field.

plasmaExpansion.f90

@@ -131,7 +131,7 @@ program plasmaExpansion
   end if
 
   t0 = 0.0_dp
-  tf = 3.0e-8_dp / t_ref
+  tf = 1.0e-7_dp / t_ref
   ! tf = 1.0e1_dp * (rf - r0) / c_s
   dt = 1.0e-2_dp*dr/c_s
   nt = nint((tf - t0) / dt) 
@@ -208,7 +208,7 @@ program plasmaExpansion
   Res = 0.0_dp
 
   ! Set boundary values
-  phi0 = 0.0e0_dp / phi_ref ! Dirichlet
+  phi0 = 1.0e2_dp / phi_ref ! Dirichlet
   phi(1) = phi0 ! Dirichlet
   ! phi0 = phi(1) ! Neumann
   allocate(f0(j0:nv))
@@ -276,6 +276,9 @@ program plasmaExpansion
     end do
     !$omp end parallel do
 
+    ! Assume quasi-neutrality to start iterating
+    n_e = Zave * n_i
+
     ! Solve Poission (maximum number of iterations, break if convergence is reached before)
     do k = 1, 1000
       ! Store previous value
@@ -309,7 +312,6 @@ program plasmaExpansion
       ! phi0=phi(1) ! Neumann
 
       ! Calculate distribution of electrons
-      ! n_e = Zave * n_i ! Quasi-neutral
       n_e = Zave(1) * n_i(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
       ! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working