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

Nothing important

Browse source
This commit is contained in:
Jorge Gonzalez 2024-11-13 12:52:54 +01:00
commit 4511aede8c
1 changed files with 9 additions and 11 deletions
Download patch file Download diff file Expand all files Collapse all files

20
plasmaExpansion.f90
View file

@ -99,7 +99,7 @@ program plasmaExpansion
! Set domain boundaries (non-dimensional units) ! Set domain boundaries (non-dimensional units)
r0 = 10.0e-6_dp / L_ref r0 = 10.0e-6_dp / L_ref
rf = 50.0e-6_dp / L_ref rf = 100.0e-6_dp / L_ref
dr = 2.0e1_dp dr = 2.0e1_dp
nr = nint((rf - r0) / dr) + 1 nr = nint((rf - r0) / dr) + 1
dr = (rf - r0) / float(nr-1) dr = (rf - r0) / float(nr-1)
@ -109,7 +109,7 @@ program plasmaExpansion
end do end do
! Set position to calculate cumulative sum of f (non-dimensional units) ! Set position to calculate cumulative sum of f (non-dimensional units)
rCum = 40.0e-6 / L_ref rCum = 80.0e-6 / L_ref
! Index for cumulative sum ! Index for cumulative sum
rCum_index = minloc(abs(r - rCum), 1) rCum_index = minloc(abs(r - rCum), 1)
@ -208,10 +208,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.0e0_dp / phi_ref ! Dirichlet
phi(1) = phi0 ! Dirichlet phi(1) = phi0 ! Dirichlet
! phi0 = phi(1) ! Neumann ! phi0 = phi(1) ! Neumann
! phiF = 0.0_dp ! Dirichlet
allocate(f0(j0:nv)) allocate(f0(j0:nv))
f0 = 0.0_dp f0 = 0.0_dp
@ -282,11 +281,6 @@ program plasmaExpansion
! Store previous value ! Store previous value
phi_old = phi phi_old = phi
! Calculate distribution of electrons
! n_e = Zave * n_i ! Quasi-neutral
n_e = Zave(1) * n_i(1) * exp((phi_old - phi0) / T_e) ! Isothermal (Boltzmann)
! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi_old-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working
! Diagonal matrix for Newton integration scheme ! Diagonal matrix for Newton integration scheme
db_dphi = n_e / T_e ! Isotropic db_dphi = n_e / T_e ! Isotropic
diag = -2.0_dp / dr**2 - db_dphi diag = -2.0_dp / dr**2 - db_dphi
@ -294,8 +288,7 @@ program plasmaExpansion
diag_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * dr) diag_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * dr)
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 - db_dphi(1) ! Neumann ! diag(nr) = 1.0_dp ! Dirichlet
! diag(nr) = 2.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 - db_dphi(nr) ! Neumann diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann
@ -315,6 +308,11 @@ program plasmaExpansion
phi = phi_old + 1.0e-1_dp*Res phi = phi_old + 1.0e-1_dp*Res
! phi0=phi(1) ! Neumann ! phi0=phi(1) ! Neumann
! Calculate distribution of electrons
! n_e = Zave * n_i ! Quasi-neutral
n_e = Zave(1) * n_i(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
! 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)
if (phiConv < 1.0e-3_dp) then if (phiConv < 1.0e-3_dp) then