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

Cartesian coordinates working as expected

Browse source
This commit is contained in:
Jorge Gonzalez 2024-10-01 09:07:23 +02:00
commit 82935b0f40
1 changed files with 47 additions and 51 deletions
Download patch file Download diff file Expand all files Collapse all files

98
plasmaExpansion.f90
View file

@ -354,13 +354,13 @@ program plasmaExpansion
! 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 = 5.0_dp * eV_to_K / Temp_ref TempF = 60.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 = c_s!sqrt(Temp0) u_bc0 = c_s!sqrt(Temp0)
u_bcF = sqrt(TempF) u_bcF = u_bc0!sqrt(TempF)
n_bc0 = n_ecr / T_to_Z(Temp0) n_bc0 = n_ecr / T_to_Z(Temp0)
n_bcF = n_ecr*1.0e-1 / T_to_Z(Temp0) n_bcF = n_bc0!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
@ -395,7 +395,7 @@ program plasmaExpansion
t0 = 0.0_dp t0 = 0.0_dp
tf = 1.0e-6_dp / t_ref tf = 1.0e-6_dp / t_ref
! tf = 1.0e1_dp * (rf - r0) / c_s ! tf = 1.0e1_dp * (rf - r0) / c_s
dt = 1.0e-2_dp*dr/c_s dt = 1.0e-2_dp*dr/(10.0_dp*c_s)
nt = nint((tf - t0) / dt) nt = nint((tf - t0) / dt)
dt = (tf - t0) / float(nt) dt = (tf - t0) / float(nt)
@ -430,7 +430,7 @@ program plasmaExpansion
u_i = 0.0_dp u_i = 0.0_dp
E_i = 0.0_dp E_i = 0.0_dp
T_i = 0.0_dp T_i = 0.0_dp
n_e = 0.0_dp n_e = 1.0_dp
Zave = 0.0_dp Zave = 0.0_dp
phi = 0.0_dp phi = 0.0_dp
phi_old = 0.0_dp phi_old = 0.0_dp
@ -445,14 +445,14 @@ program plasmaExpansion
diag_high = 0.0_dp diag_high = 0.0_dp
b = 0.0_dp b = 0.0_dp
diag = -2.0_dp / dr**2 diag = -2.0_dp / dr**2
diag_low = 1.0_dp / dr**2 * r(1:nr-1) / r(2:nr) diag_low = 1.0_dp / dr**2! * r(1:nr-1) / r(2:nr)
diag_high = 1.0_dp / dr**2 * r(2:nr) / r(1:nr-1) diag_high = 1.0_dp / dr**2! * r(2:nr) / r(1:nr-1)
! diag(1) = 1.0_dp ! Dirichlet diag(1) = 1.0_dp ! Dirichlet
! diag_high(1) = 0.0_dp ! Dirichlet diag_high(1) = 0.0_dp ! Dirichlet
diag_high(1) = 2.0_dp / dr**2 ! Neumann ! diag_high(1) = 2.0_dp / dr**2 ! Neumann
diag(nr) = 1.0_dp ! Dirichlet ! diag(nr) = 1.0_dp ! Dirichlet
diag_low(nr-1) = 0.0_dp ! Dirichlet ! diag_low(nr-1) = 0.0_dp ! Dirichlet
! diag_low(nr-1) = 2.0_dp / dr**2 ! Neumann diag_low(nr-1) = 2.0_dp / dr**2 ! Neumann
allocate(A(1:nr,1:nr)) allocate(A(1:nr,1:nr))
A = 0.0_dp A = 0.0_dp
@ -470,10 +470,9 @@ program plasmaExpansion
Res = 0.0_dp Res = 0.0_dp
! Set boundary values ! Set boundary values
! phi0 = 0.0_dp / phi_ref ! Dirichlet phi0 = 0.0_dp / phi_ref ! Dirichlet
! phi(1) = phi0 ! Dirichlet ! phi0 = phi(1) ! Neumann
phi0 = phi(1) ! Neumann ! phiF = 0.0_dp ! Dirichlet
phiF = 0.0_dp ! Dirichlet
allocate(f0(j0:nv)) allocate(f0(j0:nv))
f0 = 0.0_dp f0 = 0.0_dp
@ -520,13 +519,13 @@ program plasmaExpansion
do i = 1, nr do i = 1, nr
! Advect negative velocity ! Advect negative velocity
if (i < nr) then if (i < nr) then
f_i(i,1:j0-1) = f_i_old(i,1:j0-1) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(i+1,1:j0-1) - & f_i(i,1:j0-1) = f_i_old(i,1:j0-1) - v(1:j0-1)*dt/dr*(f_i_old(i+1,1:j0-1) - &
r(i )**2*f_i_old(i ,1:j0-1)) f_i_old(i ,1:j0-1))
end if end if
! Advect positive velocity ! Advect positive velocity
if (i > 1) then if (i > 1) then
f_i(i,j0:nv) = f_i_old(i, j0:nv) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(i , j0:nv) - & f_i(i,j0:nv) = f_i_old(i, j0:nv) - v( j0:nv)*dt/dr*(f_i_old(i , j0:nv) - &
r(i-1)**2*f_i_old(i-1, j0:nv)) f_i_old(i-1, j0:nv))
end if end if
n_i(i) = sum(f_i(i,:)*dv) n_i(i) = sum(f_i(i,:)*dv)
@ -546,29 +545,26 @@ 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, 500 do k = 1, 1000
! Store previous value ! Store previous value
phi_old = phi phi_old = phi
! Calculate distribution of electrons
n_e = Zave(1) * n_i(1) * exp((phi_old - phi0) / T_i(1))
! Diagonal matrix for Newton integration scheme ! Diagonal matrix for Newton integration scheme
diag = -2.0_dp / dr**2 - n_e diag = -2.0_dp / dr**2 - n_e
diag_low = 1.0_dp / dr**2 * r(1:nr-1) / r(2:nr) diag_low = 1.0_dp / dr**2! * r(1:nr-1) / r(2:nr)
diag_high = 1.0_dp / dr**2 * r(2:nr) / r(1:nr-1) diag_high = 1.0_dp / dr**2! * r(2:nr) / r(1:nr-1)
! diag(1) = 1.0_dp ! Dirichlet diag(1) = 1.0_dp ! Dirichlet
! diag_high(1) = 0.0_dp ! Dirichlet diag_high(1) = 0.0_dp ! Dirichlet
diag_high(1) = 2.0_dp / dr**2 - n_e(1) ! Neumann ! diag_high(1) = 2.0_dp / dr**2 - n_e(1) ! Neumann
diag(nr) = 1.0_dp ! Dirichlet ! diag(nr) = 1.0_dp ! Dirichlet
diag_low(nr-1) = 0.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) ! Neumann
! Calculate charge density ! Calculate charge density
b = -(Zave*n_i - n_e) b = -(Zave*n_i - n_e)
! Apply boundary conditions ! Apply boundary conditions
! b(1) = phi0 ! Dirichlet b(1) = phi0 ! Dirichlet
b(nr) = phiF ! Dirichlet ! b(nr) = phiF ! Dirichlet
! Calculate residual ! Calculate residual
!$omp parallel workshare !$omp parallel workshare
@ -578,39 +574,39 @@ 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) ! Ensures smooth transition in Neumann boundary
phi0 = phi(1) ! Neumann ! Calculate distribution of electrons
n_e = Zave(1) * n_i(1) * exp((phi - phi0) / T_i(1))
! 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-1_dp) then if (phiConv < 1.0e-5_dp) then
exit exit
end if end if
! Calculate new potential to ensure 0 current at the edge ! ! Calculate new potential to ensure 0 current at the edge
if (n_i(nr) > 1.0e-10_dp) then ! 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))) ! 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 ! else
phiF = phi(nr-5) ! phiF = phi(nr-5)
!
end if ! end if
end do end do
! Calculate electric field ! Calculate electric field
! E(1) = - (phi(2) - phi(1)) / dr ! Dirichlet E(1) = - (phi(2) - phi(1)) / dr ! Dirichlet
E(1) = 0.0_dp ! Neumann ! E(1) = 0.0_dp ! Neumann
!$omp parallel do !$omp parallel do
do i = 2, nr-1 do i = 2, nr-1
E(i) = - 0.5_dp*(phi(i+1) - phi(i-1)) / dr E(i) = - 0.5_dp*(phi(i+1) - phi(i-1)) / dr
end do end do
!$omp end parallel do !$omp end parallel do
E(nr) = - (phi(nr) - phi(nr-1)) / dr ! Dirichlet ! E(nr) = - (phi(nr) - phi(nr-1)) / dr ! Dirichlet
! E(nr) = 0.0_dp ! Neumann E(nr) = 0.0_dp ! Neumann
! Trick to avoid problems at the sheath
! Update intermediate f ! Update intermediate f
f_i_old = f_i f_i_old = f_i