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

Fix calculations with Z list

Browse source
This commit is contained in:
JHendrikx 2025-02-04 17:25:14 +01:00
commit f3b2c71df5
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

@ -301,7 +301,7 @@ program VlaPlEx
end if
n_i(iz,i) = sum(f_i(iz,i,:))*dv
if (n_i(1,i) > 1.0e-10_dp) then
if (n_i(iz,i) > 1.0e-10_dp) then
u_i(iz,i) = sum(v(:) *f_i(iz,i,:))*dv / n_i(iz,i)
E_i(i) = sum(v(:)**2*f_i(iz,i,:))*dv / n_i(iz,i)
T_i(iz,i) = 2.0_dp*E_i(i) - 2.0_dp*u_i(iz,i)**2
@ -315,11 +315,11 @@ program VlaPlEx
end do
!$omp end parallel do
sum_ni = sum_ni + Zave * n_i(iz,:)
sum_ni = sum_ni + Z_list(iz) * n_i(iz,:)
end do
! Assume quasi-neutrality to start iterating
n_e = 1.0_dp/nz * sum_ni
n_e = sum_ni
db_dphi = 0.0_dp
! Solve Poission (maximum number of iterations, break if convergence is reached before)
@ -337,7 +337,7 @@ program VlaPlEx
diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann
! Calculate charge density
b = - 1.0_dp/nz * sum_ni + n_e
b = - (sum_ni - n_e)
! Apply boundary conditions
b(1) = phi0 ! Dirichlet
! b(nr) = 0.0_dp ! Dirichlet
@ -398,18 +398,18 @@ program VlaPlEx
! i = 1, v<0
i = 1
if (E(i) >= 0.0_dp) then
f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,2:j0-2) - f_i_old(iz,i,1:j0-3))
f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,2:j0-2) - f_i_old(iz,i,1:j0-3))
else
f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,3:j0-1) - f_i_old(iz,i,2:j0-2))
f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,3:j0-1) - f_i_old(iz,i,2:j0-2))
end if
! i = 2, nr-1; all v
!$omp parallel do
do i = 2, nr-1
if (E(i) >= 0.0_dp) then
f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,2:nv-1) - f_i_old(iz,i,1:nv-2))
f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,2:nv-1) - f_i_old(iz,i,1:nv-2))
else
f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,3:nv) - f_i_old(iz,i,2:nv-1))
f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,3:nv) - f_i_old(iz,i,2:nv-1))
end if
@ -418,9 +418,9 @@ program VlaPlEx
! i = nr, v>=0
i = nr
if (E(i) >= 0.0_dp) then
f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,j0+1:nv-1) - f_i_old(iz,i,j0:nv-2))
f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,j0+1:nv-1) - f_i_old(iz,i,j0:nv-2))
else
f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,j0+2:nv) - f_i_old(iz,i,j0+1:nv-1))
f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,j0+2:nv) - f_i_old(iz,i,j0+1:nv-1))
end if
end do