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

n_e as boundary condition

Browse source 
Now n_e is given as the density at the boundary and n_i at the boundary
is calculated once Z is known.

This aims to eliminate the iterative process.
This commit is contained in:
Jorge Gonzalez 2025-04-10 08:49:01 +02:00
commit 0c27b98e2e
6 changed files with 38 additions and 38 deletions
Download patch file Download diff file Expand all files Collapse all files

28
vlaplex.f90
View file

@ -166,7 +166,7 @@ program VlaPlEx
write(*, '(A,ES0.4e3)') 'CFL: ', dt*vf/dr
nzMin = 3
nzMax = 14
nzMax = 13
nz = nzMax - nzMin + 1
nz = nz + 1 ! Add bin for low Z plasma
! Allocate vectors
@ -189,7 +189,7 @@ program VlaPlEx
n_e = 0.0_dp
T_e = 0.0_dp
Zave = 0.0_dp
Z_list(1) = 1.0_dp ! Low Z bin
Z_list(1) = 0.1_dp ! Low Z bin
do iz = nzMin, nzMax
Z_list(iz-nzMin+1+1) = float(iz)
end do
@ -255,23 +255,17 @@ program VlaPlEx
do t = 1, nt
time = t * dt + t0
! Find new \bar{Z}_i based on T and density
! Get boundary conditions for specific time
call boundaryConditions%get(time, n_bc, u_bc, Temp_bc)
! Reset previous value
Zave_bc_old = 0.0_dp
! Initial guess based on average table
call TtoZ%get(Temp_bc, Zave_bc)
! Find new \bar{Z}_i based on T_e = Temp_bc and n_e = n_bc
call TtoZne%get(Temp_bc, n_bc, Zave_bc)
! Assign Z(T,n) to bin
z_inj = minloc(abs(Z_list - Zave_bc),1)
Zave_bc = Z_list(z_inj)
! Start iterative process based on T, n_e table
do while (Zave_bc - Zave_bc_old > 0.1_dp)
Zave_bc_old = Zave_bc
call TtoZne%get(Temp_bc, Zave_bc * n_bc, Zave_bc)
end do
z_inj = minloc(abs(Z_list - Zave_bc),1)
Zave_bc = Z_list(z_inj)
u_bc = sqrt(Zave_bc * Temp_bc)
call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc)
! Calculate inject (sonic) speed
u_bc = sqrt(Z_list(z_inj)* Temp_bc)
! Calculate ion density to inject
n_bc = n_bc / Z_list(z_inj)
call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc, Z_list(z_inj))
! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)