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

Try to simulate expansion by reducing f at each z

Browse source
This commit is contained in:
Jorge Gonzalez 2026-01-11 18:33:38 +01:00
commit a551da69ca
1 changed files with 14 additions and 13 deletions
Download patch file Download diff file Expand all files Collapse all files

27
src/vlaplex.f90
View file

@ -14,6 +14,7 @@ program VlaPlEx
real(dp), parameter:: gamma_e_exp = 1.0_dp /(gamma_e - 1.0_dp) ! Exponent for polytropic electrons
real(dp), parameter:: gamma_e_dexp = (2.0_dp - gamma_e)/(gamma_e - 1.0_dp) ! Exponent for polytropic db_dphi
real(dp), parameter:: n_epsilon = 1.0e-16_dp
real(dp), parameter:: cosTheta = 0.995_dp
real(dp):: r0, rf
real(dp), allocatable, dimension(:):: r
@ -53,7 +54,7 @@ program VlaPlEx
real(dp), allocatable, dimension(:):: phi, phi_old, E, db_dphi
real(dp):: phiConv
real(dp):: phi0
real(dp):: T_e
real(dp):: T_e0
! real(dp):: phiF
integer:: k
@ -145,7 +146,7 @@ program VlaPlEx
E_i = 0.0_dp
T_i = 0.0_dp
n_e = 0.0_dp
T_e = 0.0_dp
T_e0 = 0.0_dp
Zave = 0.0_dp
Zave_bc_old = 0.0_dp
phi = 0.0_dp
@ -211,7 +212,7 @@ program VlaPlEx
! Get boundary conditions for specific time
call boundaryConditions%get(time, n_bc, u_bc, Temp_bc)
! Find new \bar{Z}_i based on T_e = Temp_bc and n_e = n_bc
! Find new \bar{Z}_i based on T_e0 = Temp_bc and n_e = n_bc
call Tene_to_Z%get(Temp_bc, n_bc, Zave_bc)
! Assign Z(T,n) to bin
z_inj = minloc(abs(Zlist - Zave_bc),1)
@ -229,7 +230,7 @@ program VlaPlEx
f_i_old(z_inj,1,j0:nv) = f0
f_i(:,1,j0:nv) = f_i_old(:,1,j0:nv)
T_e = Temp_bc
T_e0 = Temp_bc
! r = rf, v<0
f_i_old(:,nr,1:j0-1) = 0.0_dp
@ -247,13 +248,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*(f_i_old(iz,i+1,1:j0-1) - &
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)*cosTheta*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*(f_i_old(iz,i , j0:nv) - &
f_i_old(iz,i-1, j0:nv))
f_i(iz,i,j0:nv) = f_i_old(iz,i, j0:nv) - v( j0:nv)*cosTheta*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
@ -308,12 +309,12 @@ program VlaPlEx
! phi0=phi(1) ! Neumann
! Calculate distribution of electrons
! 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
! n_e = sum_ni(1) * exp((phi- phi0) / T_e0) ! Isothermal (Boltzmann)
n_e = sum_ni(1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e0)**gamma_e_exp !Polytropic
! Diagonal matrix for Newton integration scheme
! db_dphi = n_e / T_e ! Isothermal (Boltzmann)
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
! db_dphi = n_e / T_e0 ! Isothermal (Boltzmann)
db_dphi = sum_ni(1) / (gamma_e * T_e0) * &
(1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e0)**gamma_e_dexp !Polytropic
! Check if the solution has converged
phiConv = maxval(abs(Res),1)