JGonzalez/vlaplex
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Reorganizing calculation of db_dphi

Browse source 
I was not assuming that the first iteration in the N-R method for the
Poisson equation was fully quasi-neutral becuase db_dphi was calculated
as if we had a distribution of electrons different from ions.

It should be fixed now
This commit is contained in:
Jorge Gonzalez 2024-12-14 11:42:24 +01:00
commit 968c6ee787
1 changed files with 9 additions and 6 deletions
Download patch file Download diff file Expand all files Collapse all files

15
vlaplex.f90
View file

@ -36,6 +36,8 @@ program VlaPlEx
real(dp), parameter:: gamma_i = 1.0_dp ! Adiabatic coefficient for ions
real(dp), parameter:: m_e = 9.1093837e-31_dp ! Electron mass in kg
real(dp), parameter:: gamma_e = 4.0_dp / 3.0_dp ! Adiabatic coefficient for electrons
real(dp), parameter:: gamma_e_exp = 1.0_dp /(gamma_e - 1.0_dp) ! Exponent for polytropic electrons
real(dp), parameter:: gamma_e_dexp = (2.0_dp - gamma_e)/(gamma_e - 1.0_dp) ! Exponent for polytropic db_dphi
real(dp):: r0, rf
real(dp), allocatable, dimension(:):: r
@ -229,7 +231,7 @@ program VlaPlEx
time = t * dt + t0
call boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_bc)
call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc)
! u_bc = sqrt(Zave_bc * Temp_bc)
u_bc = sqrt(Zave_bc * Temp_bc)
! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
f0 = f0 * n_bc / (sum(f0)*dv)
@ -279,16 +281,13 @@ program VlaPlEx
! Assume quasi-neutrality to start iterating
n_e = Zave * n_i
db_dphi = 0.0_dp
! Solve Poission (maximum number of iterations, break if convergence is reached before)
do k = 1, 2000
! Store previous value
phi_old = phi
! Diagonal matrix for Newton integration scheme
! db_dphi = n_e / T_e ! Isotropic
db_dphi = Zave(1) * n_i(1) / (gamma_e * T_e) * &
(1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi_old-phi0)/T_e)**((2.0_dp - gamma_e)/(gamma_e - 1.0_dp)) ! Polytropic
diag = -2.0_dp / dr**2 - db_dphi
diag_low = 1.0_dp / dr**2 - 1.0_dp / (r(2:nr) * dr)
diag_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * dr)
@ -316,7 +315,11 @@ program VlaPlEx
! Calculate distribution of electrons
! 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)) ! Polytropic
n_e = Zave(1) * n_i(1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic
! Diagonal matrix for Newton integration scheme
! db_dphi = n_e / T_e ! Isotropic
db_dphi = Zave(1) * n_i(1) / (gamma_e * T_e) * &
(1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_dexp !Polytropic
! Check if the solution has converged
phiConv = maxval(abs(Res),1)