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

Preparing for polytropic electrons

Browse source
This commit is contained in:
Jorge Gonzalez 2024-10-16 17:56:36 +02:00
commit 47fe2fb55d
1 changed files with 9 additions and 28 deletions
Download patch file Download diff file Expand all files Collapse all files

37
plasmaExpansion.f90
View file

@ -52,7 +52,6 @@ program plasmaExpansion
real(dp):: Zave_bc ! Average charge state
real(dp):: u_bc ! Injection velocity
real(dp):: n_bc ! Injection density
! real(dp):: n_ecr ! Electron critical density for the laser
real(dp):: c_s ! Ion sound speed
type(tableBC):: boundaryConditions
character(:), allocatable:: bc_file
@ -74,6 +73,7 @@ program plasmaExpansion
real(dp), allocatable, dimension(:):: phi, phi_old, E
real(dp):: phiConv
real(dp):: phi0
real(dp):: T_e
! real(dp):: phiF
integer:: k
@ -93,16 +93,7 @@ program plasmaExpansion
phi_ref = kb * Temp_ref / qe
! Set input parameters (remember these have to be in non-dimensional units)
! Temp0 = 30.0_dp * eV_to_K / Temp_ref
! TempF = 5.0_dp * eV_to_K / Temp_ref
! Zave0 = 9.0_dp
! ZaveF = 9.0_dp
! n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref
c_s = sqrt(9.0_dp * gamma_i * 1.0_dp)
! u_bc0 = sqrt(Temp0)
! u_bcF = sqrt(TempF)
! n_bc0 = n_ecr / Zave0
! n_bcF = n_bc0*1.0e-2_dp!n_ecr*1.0e-1 / T_to_Z(Temp0)
bc_file = 'bc.csv'
call boundaryConditions%init(bc_file)
@ -146,9 +137,6 @@ program plasmaExpansion
nt = nint((tf - t0) / dt)
dt = (tf - t0) / float(nt)
! t_bc0 = nint(100.0e-9_dp / t_ref / dt)
! t_bcF = nint(150.0e-9_dp / t_ref / dt)
everyOutput = nint(1.0e-9_dp/t_ref/dt)
if (everyOutput == 0) then
everyOutput = 1
@ -178,6 +166,7 @@ program plasmaExpansion
E_i = 0.0_dp
T_i = 0.0_dp
n_e = 1.0_dp
T_e = 1.0_dp
Zave = 0.0_dp
phi = 0.0_dp
phi_old = 0.0_dp
@ -235,20 +224,10 @@ program plasmaExpansion
! Main loop
do t = 1, nt
time = t * dt + t0
! if (t > t_bc0 .and. t <= t_bcF) then
! Temp_bc = (TempF - Temp0) / float(t_bcF - t_bc0)*(t - t_bc0) + Temp0
! Zave_bc = (ZaveF - Zave0) / float(t_bcF - t_bc0)*(t - t_bc0) + Zave0
! u_bc = (u_bcF - u_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + u_bc0
! n_bc = (n_bcF - n_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + n_bc0
! else if (t > t_bcF) then
! Temp_bc = TempF
! Zave_bc = ZaveF
! u_bc = u_bcF
! n_bc = n_bcF
! end if
call boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_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) = 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)
! Boundary conditions
@ -299,7 +278,7 @@ program plasmaExpansion
phi_old = phi
! Diagonal matrix for Newton integration scheme
diag = -2.0_dp / dr**2 - n_e
diag = -2.0_dp / dr**2 - n_e / T_e
diag_low = 1.0_dp / dr**2 - 1.0_dp / (r(2:nr) * 2.0_dp * dr)
diag_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * 2.0_dp * dr)
diag(1) = 1.0_dp ! Dirichlet
@ -307,7 +286,7 @@ program plasmaExpansion
! diag_high(1) = 2.0_dp / dr**2 - n_e(1) ! Neumann
! diag(nr) = 1.0_dp ! Dirichlet
! diag_low(nr-1) = 0.0_dp ! Dirichlet
diag_low(nr-1) = 2.0_dp / dr**2 - n_e(nr) ! Neumann
diag_low(nr-1) = 2.0_dp / dr**2 - n_e(nr) / T_e ! Neumann
! Calculate charge density
b = -(Zave*n_i - n_e)
@ -326,7 +305,9 @@ program plasmaExpansion
! phi0=phi(1)
! Calculate distribution of electrons
n_e = Zave(1) * n_i(1) * exp((phi - phi0) / 1.0_dp)
T_e = 1.0_dp
n_e = Zave(1) * n_i(1) * exp((phi - phi0) / T_e)
! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_i(1))**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working
! Check if the solution has converged
phiConv = maxval(abs(Res),1)