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

Polytropic electrons finally working!

Browse source 
There was a problem with one parenthesis in the wrong location. Now this
works and I am extremely happy, even if my face doesn't show it.
This commit is contained in:
Jorge Gonzalez 2024-12-11 08:50:54 +01:00
commit b093ca5298
1 changed files with 6 additions and 6 deletions
Download patch file Download diff file Expand all files Collapse all files

12
plasmaExpansion.f90
View file

@ -35,7 +35,7 @@ program plasmaExpansion
real(dp), parameter:: m_i = 1.9712258e-25_dp ! Tin atom mass in kg
real(dp), parameter:: gamma_i = 1.0_dp ! Adiabatic coefficient for ions
real(dp), parameter:: m_e = 9.1093837e-31_dp ! Electron mass in kg
real(dp), parameter:: gamma_e = 2.0_dp ! Adiabatic coefficient for electrons
real(dp), parameter:: gamma_e = 4.0_dp / 3.0_dp ! Adiabatic coefficient for electrons
real(dp):: r0, rf
real(dp), allocatable, dimension(:):: r
@ -285,9 +285,9 @@ program plasmaExpansion
phi_old = phi
! Diagonal matrix for Newton integration scheme
db_dphi = n_e / T_e ! Isotropic
! db_dphi = Zave(1) * n_i(1) / (gamma_e * T_e) * &
! (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi_old-phi0)/T_e)**((2.0_dp - gamma_e)/(gamma_e - 1.0_dp))) ! Polytropic
! db_dphi = n_e / T_e ! Isotropic
db_dphi = Zave(1) * n_i(1) / (gamma_e * T_e) * &
(1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi_old-phi0)/T_e)**((2.0_dp - gamma_e)/(gamma_e - 1.0_dp)) ! Polytropic
diag = -2.0_dp / dr**2 - db_dphi
diag_low = 1.0_dp / dr**2 - 1.0_dp / (r(2:nr) * dr)
diag_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * dr)
@ -314,8 +314,8 @@ program plasmaExpansion
! phi0=phi(1) ! Neumann
! Calculate distribution of electrons
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))) ! Polytropic
! 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)) ! Polytropic
! Check if the solution has converged
phiConv = maxval(abs(Res),1)