JGonzalez
052a4dc05e
Small set of corrections to the tag v2.5. Includes changes to python scripts to plot data. Includes new 'fast' setup conditions that allow to output cases in an hour or so with still a good CFL condition and grid resolution.
466 lines
14 KiB
Fortran
466 lines
14 KiB
Fortran
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.0e1_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 = 2.0e-7_dp / t_ref
|
|
! tf = 1.0e1_dp * (rf - r0) / c_s
|
|
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 = 14
|
|
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) = 1.0_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
|
|
|
|
! Find new \bar{Z}_i based on T and density
|
|
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)
|
|
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)
|
|
z_inj = minloc(abs(Z_list - Zave_bc),1)
|
|
Zave_bc = Z_list(z_inj)
|
|
end do
|
|
u_bc = sqrt(Zave_bc * Temp_bc)
|
|
call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc)
|
|
|
|
! 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 ! Isotropic
|
|
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
|
|
|