From af5a4fae27b90411446a20d425f7f5e463a5a615 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Thu, 26 Jun 2025 12:56:35 +0200 Subject: [PATCH] Remove geometry coefficients --- src/vlaplex.f90 | 25 ++++++++----------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/src/vlaplex.f90 b/src/vlaplex.f90 index 7fc9ebe..4ff8f24 100644 --- a/src/vlaplex.f90 +++ b/src/vlaplex.f90 @@ -163,8 +163,8 @@ program VlaPlEx b = 0.0_dp db_dphi = 0.0_dp diag = -2.0_dp / dr**2 - 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) + 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) diag(1) = 1.0_dp ! Dirichlet diag_high(1) = 0.0_dp ! Dirichlet ! diag_high(1) = 2.0_dp / dr**2 ! Neumann @@ -247,13 +247,13 @@ program VlaPlEx do i = 1, nr ! Advect negative velocity if (i < nr) then - f_i(iz,i,1:j0-1) = f_i_old(iz,i,1:j0-1) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(iz,i+1,1:j0-1) - & - r(i )**2*f_i_old(iz,i ,1:j0-1)) + f_i(iz,i,1:j0-1) = f_i_old(iz,i,1:j0-1) - v(1:j0-1)*dt/dr*(f_i_old(iz,i+1,1:j0-1) - & + f_i_old(iz,i ,1:j0-1)) end if ! Advect positive velocity if (i > 1) then - f_i(iz,i,j0:nv) = f_i_old(iz,i, j0:nv) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(iz,i , j0:nv) - & - r(i-1)**2*f_i_old(iz,i-1, j0:nv)) + f_i(iz,i,j0:nv) = f_i_old(iz,i, j0:nv) - v( j0:nv)*dt/dr*(f_i_old(iz,i , j0:nv) - & + f_i_old(iz,i-1, j0:nv)) end if n_i(iz,i) = sum(f_i(iz,i,:))*dv @@ -283,8 +283,8 @@ program VlaPlEx phi_old = phi 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) + 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) diag(1) = 1.0_dp ! Dirichlet diag_high(1) = 0.0_dp ! Dirichlet ! diag(nr) = 1.0_dp ! Dirichlet @@ -322,15 +322,6 @@ program VlaPlEx end if - ! ! Calculate new potential to ensure 0 current at the edge - ! if (n_i(nr) > n_epsilon) then - ! phiF = phi0 + T_e * log((2.0_dp*sqrt(pi)*Zave(nr)*n_i(nr)*u_i(nr)) / (Zave(1)*n_i(1)*sqrt(m_i*T_e/m_e))) - ! - ! else - ! phiF = phi(nr-5) - ! - ! end if - end do ! Calculate electric field