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

Input file

Browse source 
This should not have been a single commit.

Now we have input files.

Also, I've restructured the code and renamed modules.
This commit is contained in:
Jorge Gonzalez 2025-04-10 11:56:05 +02:00
commit 50b4258c8f
22 changed files with 446 additions and 251 deletions
Download patch file Download diff file Expand all files Collapse all files

459
vlaplex.f90
View file

@ -1,459 +0,0 @@
module eos
use constantParameters, only: dp
implicit none
private
public:: T_to_Z
contains
pure function T_to_Z(T) result(Z)
use constantParameters, only: eV_to_K
use referenceValues, only: Temp_ref
implicit none
real(dp), intent(in):: T
real(dp):: Z
Z = 22.5 * (Temp_ref * T / eV_to_K / 100.0)**0.6
end function T_to_Z
end module eos
program VlaPlEx
use constantParameters, only: dp, kb, qe, eps_0, ev_to_K, cm3_to_m3, PI
use output
use referenceValues
use eos, only: T_to_Z
use moduleTableBC
use moduleTableTtoZ
use moduleTableTtoZne
use omp_lib
implicit none
real(dp), parameter:: m_i = 1.9712258e-25_dp ! Tin atom mass in kg
real(dp), parameter:: gamma_i = 1.0_dp ! Adiabatic coefficient for ions
real(dp), parameter:: m_e = 9.1093837e-31_dp ! Electron mass in kg
real(dp), parameter:: gamma_e = 4.0_dp / 3.0_dp ! Adiabatic coefficient for electrons
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):: r0, rf
real(dp), allocatable, dimension(:):: r
real(dp):: v0, vf
real(dp), allocatable, dimension(:):: v
real(dp):: t0, tf
real(dp):: time
real(dp):: dr, dv, dt
integer:: nr, nv, nt, nz
integer:: nzMin, nzMax
integer:: i, iz, j, t, z_inj
integer:: j0 ! First integer of positive velocity
real(dp):: Temp_bc ! Temperature
real(dp):: Zave_bc, Zave_bc_old ! Average charge state
real(dp):: u_bc ! Injection velocity
real(dp):: n_bc ! Injection density
real(dp):: c_s ! Ion sound speed
type(tableBC):: boundaryConditions
character(:), allocatable:: bc_file
type(tableTtoZ):: TtoZ
character(:), allocatable:: TtoZ_file
type(tableTtoZne):: TtoZne
character(:), allocatable:: TtoZne_file
real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old
real(dp), allocatable, dimension(:):: f0 ! Boundary at r = x_0
real(dp), allocatable, dimension(:,:):: n_i
real(dp), allocatable, dimension(:):: sum_ni
real(dp), allocatable, dimension(:,:):: u_i
real(dp), allocatable, dimension(:):: E_i
real(dp), allocatable, dimension(:,:):: T_i
real(dp), allocatable, dimension(:):: n_e
real(dp), allocatable, dimension(:):: Zave
real(dp), allocatable, dimension(:):: Z_list
real(dp), allocatable, dimension(:):: diag, diag_low, diag_high
real(dp), allocatable, dimension(:,:):: A
real(dp), allocatable, dimension(:):: Res
real(dp), allocatable, dimension(:):: b
integer:: info
real(dp), allocatable, dimension(:):: phi, phi_old, E, db_dphi
real(dp):: phiConv
real(dp):: phi0
real(dp):: T_e
! real(dp):: phiF
integer:: k
real(dp), allocatable, dimension(:,:):: fCum_i
real(dp):: rCum
integer:: rCum_index
! Set number of threads
call omp_set_num_threads(16)
! Set reference numbers (in SI units)
Temp_ref = 30.0_dp * eV_to_K
n_ref = 1.0e20_dp * cm3_to_m3
t_ref = sqrt(eps_0 * m_i / (n_ref * 1.0_dp * qe**2)) ! 1.0_dp represents Z = 1 for reference values
u_ref = sqrt(kb * Temp_ref / m_i)
L_ref = u_ref * t_ref
phi_ref = kb * Temp_ref / qe
! Set input parameters (remember these have to be in non-dimensional units)
c_s = sqrt(11.0_dp * gamma_i * 1.0_dp)
bc_file = 'bc.csv'
call boundaryConditions%init(bc_file)
TtoZ_file = 'average_TtoZ_curve.csv'
call TtoZ%init(TtoZ_file)
TtoZne_file = 'Zave_values_per_Ne.csv'
call TtoZne%init(TtoZne_file)
! Set domain boundaries (non-dimensional units)
r0 = 10.0e-6_dp / L_ref
rf = 2.0e-3_dp / L_ref
dr = 1.0e-6_dp / L_ref
nr = nint((rf - r0) / dr) + 1
dr = (rf - r0) / float(nr-1)
allocate(r(1:nr))
do i = 1, nr
r(i) = dr * float(i-1) + r0
end do
! Set position to calculate cumulative sum of f (non-dimensional units)
rCum = 1.0e-3 / L_ref
! Index for cumulative sum
rCum_index = minloc(abs(r - rCum), 1)
v0 =-0.5e1_dp*c_s
vf = 1.5e1_dp*c_s
dv = 2.0e-1_dp
nv = nint((vf - v0) / dv) + 1
dv = (vf - v0) / float(nv-1)
allocate(v(1:nv))
do j = 1, nv
v(j) = dv * float(j-1) + v0
end do
! Shift v mesh so it passes by 0
v = v - (minval(abs(v)))
j0 = minloc(abs(v), 1)
if (v(j0) < 0.0_dp) then
j0 = j0 + 1
end if
t0 = 0.0_dp
tf = 3.0e-7_dp / t_ref
dt = 5.0e-2_dp*dr/c_s
nt = nint((tf - t0) / dt)
dt = (tf - t0) / float(nt)
everyOutput = nint(1.0e-9_dp/t_ref/dt)
if (everyOutput == 0) then
everyOutput = 1
end if
everyWrite = everyOutput/10
if (everyWrite == 0) then
everyWrite = 1
end if
write(*, '(A,ES0.4e3)') 'CFL: ', dt*vf/dr
nzMin = 3
nzMax = 13
nz = nzMax - nzMin + 1
nz = nz + 1 ! Add bin for low Z plasma
! Allocate vectors
allocate(f_i(1:nz,1:nr,1:nv), f_i_old(1:nz,1:nr,1:nv))
allocate(n_i(1:nz,1:nr))
allocate(sum_ni(1:nr))
allocate(u_i(1:nz,1:nr), E_i(1:nr), T_i(1:nz,1:nr))
allocate(Zave(1:nr))
allocate(Z_list(1:nz))
allocate(n_e(1:nr))
allocate(phi(1:nr), phi_old(1:nr), E(1:nr))
allocate(fCum_i(1:nz,1:nv))
f_i = 0.0_dp
f_i_old = 0.0_dp
n_i = 0.0_dp
sum_ni = 0.0_dp
u_i = 0.0_dp
E_i = 0.0_dp
T_i = 0.0_dp
n_e = 0.0_dp
T_e = 0.0_dp
Zave = 0.0_dp
Z_list(1) = 0.1_dp ! Low Z bin
do iz = nzMin, nzMax
Z_list(iz-nzMin+1+1) = float(iz)
end do
Zave_bc_old = 0.0_dp
phi = 0.0_dp
phi_old = 0.0_dp
E = 0.0_dp
fCum_i = 0.0_dp
! Allocate matrix for Poisson equation
allocate(diag(1:nr), diag_low(1:nr-1), diag_high(1:nr-1))
allocate(b(1:nr))
allocate(db_dphi(1:nr))
diag = 0.0_dp
diag_low = 0.0_dp
diag_high = 0.0_dp
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(1) = 1.0_dp ! Dirichlet
diag_high(1) = 0.0_dp ! Dirichlet
! diag_high(1) = 2.0_dp / dr**2 ! Neumann
! diag(nr) = 1.0_dp ! Dirichlet
! diag_low(nr-1) = 0.0_dp ! Dirichlet
diag_low(nr-1) = 2.0_dp / dr**2 ! Neumann
allocate(A(1:nr,1:nr))
A = 0.0_dp
A(1,1) = diag(1)
A(1,2) = diag_high(1)
do i = 2, nr - 1
A(i, i-1) = diag_low(i-1)
A(i, i) = diag(i)
A(i, i+1) = diag_high(i)
end do
A(nr,nr-1) = diag_low(nr-1)
A(nr,nr) = diag(nr)
allocate(Res(1:nr))
Res = 0.0_dp
! Set boundary values
phi0 = 1.0e2_dp / phi_ref ! Dirichlet
phi(1) = phi0 ! Dirichlet
! phi0 = phi(1) ! Neumann
allocate(f0(j0:nv))
f0 = 0.0_dp
! Output initial values
call createPath()
call setTimeFormat(nt)
t = 0
call writeOutputRef()
! call writeOutputF(t, dt, nr, r, nv, v, f_i_old)
call writeOutputFCum(t, dt, nz, r(rCum_index), nv, v, fCum_i, Z_list)
call writeOutputPhi(t, dt, nr, r, phi, E, n_e)
call writeOutputMom(t, dt, nz, nr, r, n_i, u_i, T_i, Z_list)
call writeOutputZList(nz, Z_list)
! Main loop
do t = 1, nt
time = t * dt + t0
! 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
call TtoZne%get(Temp_bc, n_bc, Zave_bc)
! Assign Z(T,n) to bin
z_inj = minloc(abs(Z_list - Zave_bc),1)
! 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)
f0 = f0 * n_bc / (sum(f0)*dv)
f_i_old(:,1,j0:nv) = 0.0_dp
f_i_old(z_inj,1,j0:nv) = f0
f_i(:,1,j0:nv) = f_i_old(:,1,j0:nv)
T_e = Temp_bc
! r = rf, v<0
f_i_old(:,nr,1:j0-1) = 0.0_dp
f_i(:,nr,1:j0-1) = f_i_old(:,nr,1:j0-1)
! set edge velocities to 0
f_i_old(:,:,1) = 0.0_dp
f_i_old(:,:,nv) = 0.0_dp
sum_ni = 0.0_dp
! Advect in the r direction
do iz = 1, nz
if (all(n_i(iz,:) < n_epsilon) .and. iz .ne. z_inj) then
cycle
end if
!$omp parallel do
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))
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))
end if
n_i(iz,i) = sum(f_i(iz,i,:))*dv
if (n_i(iz,i) > n_epsilon) 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
else
f_i(iz,i,:) = 0.0_dp
n_i(iz,i) = 0.0_dp
u_i(iz,i) = 0.0_dp
T_i(iz,i) = 0.0_dp
end if
end do
!$omp end parallel do
sum_ni = sum_ni + Z_list(iz) * n_i(iz,:)
end do
! Assume quasi-neutrality to start iterating
n_e = sum_ni
db_dphi = 0.0_dp
! Solve Poission (maximum number of iterations, break if convergence is reached before)
do k = 1, 2000
! Store previous value
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(1) = 1.0_dp ! Dirichlet
diag_high(1) = 0.0_dp ! Dirichlet
! diag(nr) = 1.0_dp ! Dirichlet
! diag_low(nr-1) = 0.0_dp ! Dirichlet
diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann
! Calculate charge density
b = - (sum_ni - n_e)
! Apply boundary conditions
b(1) = phi0 ! Dirichlet
! b(nr) = 0.0_dp ! Dirichlet
! Calculate residual
!$omp parallel workshare
Res = -(MATMUL(A, phi_old) - b)
!$omp end parallel workshare
! Iterate system
call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info)
phi = phi_old + Res
! 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
! 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
! Check if the solution has converged
phiConv = maxval(abs(Res),1)
if (phiConv < 1.0e-6_dp) then
exit
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
E(1) = - (phi(2) - phi(1)) / dr ! Dirichlet
! E(1) = 0.0_dp ! Neumann
!$omp parallel do
do i = 2, nr-1
E(i) = - 0.5_dp*(phi(i+1) - phi(i-1)) / dr
end do
!$omp end parallel do
! E(nr) = - (phi(nr) - phi(nr-1)) / dr ! Dirichlet
E(nr) = 0.0_dp ! Neumann
! Update intermediate f
f_i_old = f_i
do iz = 1, nz
if (all(n_i(iz,:) < n_epsilon) .and. iz .ne. z_inj) then
cycle
end if
! Advect in the v direction
! 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) - 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) - 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) - 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) - Z_list(iz)*E(i)*dt/dv*(f_i_old(iz,i,3:nv) - f_i_old(iz,i,2:nv-1))
end if
end do
!$omp end parallel do
! 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) - 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) - 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
! Reset values for next iteration
f_i_old = f_i
do iz = 1, nz
if (all(n_i(iz,:) < n_epsilon) .and. iz .ne. z_inj) then
cycle
end if
fCum_i(iz,:) = fCum_i(iz,:) + f_i_old(iz,rCum_index,:)
end do
! Write output
if (mod(t,everyOutput) == 0 .or. t == nt) then
! call writeOutputF(t, dt, nz, nr, r, nv, v, f_i_old, Z_list)
call writeOutputPhi(t, dt, nr, r, phi, E, n_e)
call writeOutputMom(t, dt, nz, nr, r, n_i, u_i, T_i, Z_list)
call writeOutputFCum(t, dt, nz, r(rCum_index), nv, v, fCum_i, Z_list)
end if
! Write progress
if (mod(t,everyWrite) == 0) then
write (*, '(I10, A, I10)' ) t, '/', nt
write (*, '(A, ES0.4e3,","ES0.4e3)') 'phi max,min: ', maxval(phi)*phi_ref, minval(phi)*phi_ref
end if
end do
end program VlaPlEx