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

Remove geometry coefficients

Browse source
This commit is contained in:
Jorge Gonzalez 2025-06-26 12:56:35 +02:00
commit af5a4fae27
1 changed files with 8 additions and 17 deletions
Download patch file Download diff file Expand all files Collapse all files

25
src/vlaplex.f90
View file

@ -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