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

Cleaning

Browse source 
Trying to do polytropic electrons.
This commit is contained in:
Jorge Gonzalez 2024-10-18 09:51:15 +02:00
commit 49edcd0df9
1 changed files with 8 additions and 5 deletions
Download patch file Download diff file Expand all files Collapse all files

13
plasmaExpansion.f90
View file

@ -70,7 +70,7 @@ program plasmaExpansion
real(dp), allocatable, dimension(:):: b real(dp), allocatable, dimension(:):: b
integer:: info integer:: info
real(dp), allocatable, dimension(:):: phi, phi_old, E real(dp), allocatable, dimension(:):: phi, phi_old, E, db_dphi
real(dp):: phiConv real(dp):: phiConv
real(dp):: phi0 real(dp):: phi0
real(dp):: T_e real(dp):: T_e
@ -131,7 +131,7 @@ program plasmaExpansion
end if end if
t0 = 0.0_dp t0 = 0.0_dp
tf = 1.0e-6_dp / t_ref tf = 2.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/c_s
nt = nint((tf - t0) / dt) nt = nint((tf - t0) / dt)
@ -176,10 +176,12 @@ program plasmaExpansion
! Allocate matrix for Poisson equation ! Allocate matrix for Poisson equation
allocate(diag(1:nr), diag_low(1:nr-1), diag_high(1:nr-1)) allocate(diag(1:nr), diag_low(1:nr-1), diag_high(1:nr-1))
allocate(b(1:nr)) allocate(b(1:nr))
allocate(db_dphi(1:nr))
diag = 0.0_dp diag = 0.0_dp
diag_low = 0.0_dp diag_low = 0.0_dp
diag_high = 0.0_dp diag_high = 0.0_dp
b = 0.0_dp b = 0.0_dp
db_dphi = 0.0_dp
diag = -2.0_dp / dr**2 diag = -2.0_dp / dr**2
diag_low = 1.0_dp / dr**2 - 1.0_dp / (r(2:nr) * 2.0_dp * dr) 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_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * 2.0_dp * dr)
@ -278,7 +280,8 @@ program plasmaExpansion
phi_old = phi phi_old = phi
! Diagonal matrix for Newton integration scheme ! Diagonal matrix for Newton integration scheme
diag = -2.0_dp / dr**2 - n_e / T_e db_dphi = n_e / T_e ! Isotropic
diag = -2.0_dp / dr**2 - db_dphi
diag_low = 1.0_dp / dr**2 - 1.0_dp / (r(2:nr) * 2.0_dp * dr) 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_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * 2.0_dp * dr)
diag(1) = 1.0_dp ! Dirichlet diag(1) = 1.0_dp ! Dirichlet
@ -286,7 +289,7 @@ program plasmaExpansion
! 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) / T_e ! Neumann diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann
! Calculate charge density ! Calculate charge density
b = -(Zave*n_i - n_e) b = -(Zave*n_i - n_e)
@ -307,7 +310,7 @@ program plasmaExpansion
! Calculate distribution of electrons ! Calculate distribution of electrons
T_e = 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) * 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 ! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working
! Check if the solution has converged ! Check if the solution has converged
phiConv = maxval(abs(Res),1) phiConv = maxval(abs(Res),1)