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

Partial case with Poisson equation

Browse source 
It reproduces Diko's peak!
This commit is contained in:
Jorge Gonzalez 2024-09-28 20:55:03 +02:00
commit fea4f21e17
1 changed files with 21 additions and 20 deletions
Download patch file Download diff file Expand all files Collapse all files

41
plasmaExpansion.f90
View file

@ -21,7 +21,7 @@ module referenceValues
implicit none implicit none
real(dp):: L_ref, t_ref, n_ref, u_ref, Temp_ref ! Reference values real(dp):: L_ref, t_ref, n_ref, u_ref, Temp_ref ! Reference values
real(dp):: phi_ref, p_ref ! Reference values real(dp):: phi_ref ! Reference values
end module referenceValues end module referenceValues
@ -291,7 +291,8 @@ program plasmaExpansion
use omp_lib use omp_lib
implicit none implicit none
real(dp), parameter:: m_i = 1.9712258e-25_dp ! Tin aton mass in kg real(dp), parameter:: m_i = 1.9712258e-25_dp ! Tin atom mass in kg
real(dp), parameter:: m_e = 9.1093837e-31_dp ! Electron mass in kg
real(dp), parameter:: gam = 1.0_dp ! Adiabatic coefficient real(dp), parameter:: gam = 1.0_dp ! Adiabatic coefficient
real(dp):: r0, rf real(dp):: r0, rf
@ -330,7 +331,7 @@ program plasmaExpansion
real(dp), allocatable, dimension(:):: phi, phi_old, E real(dp), allocatable, dimension(:):: phi, phi_old, E
real(dp):: phiConv real(dp):: phiConv
real(dp):: phi0, phiF, JF real(dp):: phi0, phiF
integer:: k integer:: k
real(dp), allocatable, dimension(:):: fCum_i real(dp), allocatable, dimension(:):: fCum_i
@ -347,20 +348,19 @@ program plasmaExpansion
u_ref = sqrt(kb * Temp_ref / m_i) u_ref = sqrt(kb * Temp_ref / m_i)
L_ref = u_ref * t_ref L_ref = u_ref * t_ref
phi_ref = kb * Temp_ref / qe phi_ref = kb * Temp_ref / qe
p_ref = n_ref * kb * Temp_ref
! Set position to calculate cumulative sum of f (non-dimensional units) ! Set position to calculate cumulative sum of f (non-dimensional units)
rCum = 5.0e-3 / L_ref rCum = 5.0e-3 / L_ref
! Set input parameters (remember these have to be in non-dimensional units) ! Set input parameters (remember these have to be in non-dimensional units)
Temp0 = 60.0_dp * eV_to_K / Temp_ref Temp0 = 60.0_dp * eV_to_K / Temp_ref
TempF = 60.0_dp * eV_to_K / Temp_ref TempF = 5.0_dp * eV_to_K / Temp_ref
n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref
c_s = sqrt(T_to_Z(Temp0) * gam * Temp0) c_s = sqrt(T_to_Z(Temp0) * gam * Temp0)
u_bc0 = sqrt(Temp0) u_bc0 = c_s!sqrt(Temp0)
u_bcF = sqrt(TempF) u_bcF = sqrt(TempF)
n_bc0 = n_ecr / T_to_Z(Temp0) n_bc0 = n_ecr / T_to_Z(Temp0)
n_bcF = n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0) n_bcF = n_ecr*1.0e-1 / T_to_Z(Temp0)
! Set domain boundaries (non-dimensional units) ! Set domain boundaries (non-dimensional units)
r0 = 200.0e-6_dp / L_ref r0 = 200.0e-6_dp / L_ref
@ -399,8 +399,8 @@ program plasmaExpansion
nt = nint((tf - t0) / dt) nt = nint((tf - t0) / dt)
dt = (tf - t0) / float(nt) dt = (tf - t0) / float(nt)
t_bc0 = nint( 20.0e-9_dp / t_ref / dt) t_bc0 = nint(100.0e-9_dp / t_ref / dt)
t_bcF = nint( 25.0e-9_dp / t_ref / dt) t_bcF = nint(105.0e-9_dp / t_ref / dt)
everyOutput = nint(1.0e-9_dp/t_ref/dt) everyOutput = nint(1.0e-9_dp/t_ref/dt)
if (everyOutput == 0) then if (everyOutput == 0) then
@ -474,7 +474,6 @@ program plasmaExpansion
! phi(1) = phi0 ! Dirichlet ! phi(1) = phi0 ! Dirichlet
phi0 = phi(1) ! Neumann phi0 = phi(1) ! Neumann
phiF = 0.0_dp ! Dirichlet phiF = 0.0_dp ! Dirichlet
JF = 0.0_dp ! Current at the edge for floating potential
allocate(f0(j0:nv)) allocate(f0(j0:nv))
f0 = 0.0_dp f0 = 0.0_dp
@ -547,7 +546,7 @@ program plasmaExpansion
!$omp end parallel do !$omp end parallel do
! Solve Poission (maximum number of iterations, break if convergence is reached before) ! Solve Poission (maximum number of iterations, break if convergence is reached before)
do k = 1, 100 do k = 1, 500
! Store previous value ! Store previous value
phi_old = phi phi_old = phi
@ -579,23 +578,25 @@ program plasmaExpansion
! Iterate system ! Iterate system
call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info) call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info)
phi = phi_old + Res phi = phi_old + Res
phi(1) = phi(2) ! Neumann phi(1) = phi(2) ! Ensures smooth transition in Neumann boundary
phi0 = phi(1) ! Neumann phi0 = phi(1) ! Neumann
! phi(nr) = phi(nr-5) ! Dirichlet quasi-floating
! Calculate updated current
JF = Zave(nr)*n_i(nr)*u_i(nr) - n_e(nr)*sqrt(qe*T_i(1)*Temp_ref/9.1e-31_dp)/u_ref
! Iterate floating potential
phiF = phi(nr-5)! + 1.0e-2_dp*JF
! Check if the solution has converged ! Check if the solution has converged
phiConv = maxval(abs(Res),1) phiConv = maxval(abs(Res),1)
if (phiConv < 1.0e-2_dp) then if (phiConv < 1.0e-1_dp) then
exit exit
end if end if
! Calculate new potential to ensure 0 current at the edge
if (n_i(nr) > 1.0e-10_dp) then
phiF = phi0 + T_i(1) * log((2.0_dp*sqrt(pi)*Zave(nr)*n_i(nr)*u_i(nr)) / (Zave(1)*n_i(1)*sqrt(m_i*T_i(1)/m_e)))
else
phiF = phi(nr-5)
end if
end do end do
! Calculate electric field ! Calculate electric field