From f3b2c71df5c7d374670b06d5c66177a76c68c0d7 Mon Sep 17 00:00:00 2001 From: JHendrikx Date: Tue, 4 Feb 2025 17:25:14 +0100 Subject: [PATCH] Fix calculations with Z list --- vlaplex.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/vlaplex.f90 b/vlaplex.f90 index 453b667..736c5de 100644 --- a/vlaplex.f90 +++ b/vlaplex.f90 @@ -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