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

Change b_i to sum_ni

Browse source
This commit is contained in:
JHendrikx 2025-02-04 12:10:58 +01:00
commit d7c23d5577
1 changed files with 10 additions and 10 deletions
Download patch file Download diff file Expand all files Collapse all files

20
vlaplex.f90
View file

@ -61,7 +61,7 @@ program VlaPlEx
real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old
real(dp), allocatable, dimension(:):: f0 ! Boundary at r = x_0
real(dp), allocatable, dimension(:,:):: n_i
real(dp), allocatable, dimension(:):: b_i
real(dp), allocatable, dimension(:):: sum_ni
real(dp), allocatable, dimension(:,:):: u_i
real(dp), allocatable, dimension(:):: E_i
real(dp), allocatable, dimension(:,:):: T_i
@ -161,7 +161,7 @@ program VlaPlEx
! Allocate vectors
allocate(f_i(1:nz,1:nr,1:nv), f_i_old(1:nz,1:nr,1:nv))
allocate(n_i(1:nz,1:nr))
allocate(b_i(1:nr))
allocate(sum_ni(1:nr))
allocate(u_i(1:nz,1:nr), E_i(1:nr), T_i(1:nz,1:nr))
allocate(Zave(1:nr))
allocate(n_e(1:nr))
@ -170,7 +170,7 @@ program VlaPlEx
f_i = 0.0_dp
f_i_old = 0.0_dp
n_i = 0.0_dp
b_i = 0.0_dp
sum_ni = 0.0_dp
u_i = 0.0_dp
E_i = 0.0_dp
T_i = 0.0_dp
@ -266,7 +266,7 @@ program VlaPlEx
! set edge velocities to 0
f_i_old(:,:,1) = 0.0_dp
f_i_old(:,:,nv) = 0.0_dp
b_i = 0.0_dp
sum_ni = 0.0_dp
! Advect in the r direction
!$omp parallel do
do iz = 1, nz
@ -296,12 +296,12 @@ program VlaPlEx
end if
end do
b_i = b_i + Zave * n_i(iz,:)
sum_ni = sum_ni + Zave * n_i(iz,:)
end do
!$omp end parallel do
! Assume quasi-neutrality to start iterating
n_e = 1.0_dp/nz * b_i
n_e = 1.0_dp/nz * sum_ni
db_dphi = 0.0_dp
! Solve Poission (maximum number of iterations, break if convergence is reached before)
@ -319,7 +319,7 @@ program VlaPlEx
diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann
! Calculate charge density
b = - 1.0_dp/nz * b_i + n_e
b = - 1.0_dp/nz * sum_ni + n_e
! Apply boundary conditions
b(1) = phi0 ! Dirichlet
! b(nr) = 0.0_dp ! Dirichlet
@ -335,11 +335,11 @@ program VlaPlEx
! phi0=phi(1) ! Neumann
! 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) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic
! n_e = sum_ni(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
n_e = sum_ni(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,1) / (gamma_e * T_e) * &
db_dphi = sum_ni(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