From 968c6ee7879b4a89ced2f813917e2dabcdb11ef8 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 14 Dec 2024 11:42:24 +0100 Subject: [PATCH] Reorganizing calculation of db_dphi I was not assuming that the first iteration in the N-R method for the Poisson equation was fully quasi-neutral becuase db_dphi was calculated as if we had a distribution of electrons different from ions. It should be fixed now --- vlaplex.f90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/vlaplex.f90 b/vlaplex.f90 index 61717d3..869ce81 100644 --- a/vlaplex.f90 +++ b/vlaplex.f90 @@ -36,6 +36,8 @@ program VlaPlEx 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 = 4.0_dp / 3.0_dp ! Adiabatic coefficient for electrons + real(dp), parameter:: gamma_e_exp = 1.0_dp /(gamma_e - 1.0_dp) ! Exponent for polytropic electrons + real(dp), parameter:: gamma_e_dexp = (2.0_dp - gamma_e)/(gamma_e - 1.0_dp) ! Exponent for polytropic db_dphi real(dp):: r0, rf real(dp), allocatable, dimension(:):: r @@ -229,7 +231,7 @@ program VlaPlEx time = t * dt + t0 call boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_bc) call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc) - ! u_bc = sqrt(Zave_bc * Temp_bc) + u_bc = sqrt(Zave_bc * Temp_bc) ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) f0 = f0 * n_bc / (sum(f0)*dv) @@ -279,16 +281,13 @@ program VlaPlEx ! Assume quasi-neutrality to start iterating n_e = Zave * n_i + db_dphi = 0.0_dp ! Solve Poission (maximum number of iterations, break if convergence is reached before) do k = 1, 2000 ! Store previous value 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 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) @@ -316,7 +315,11 @@ program VlaPlEx ! 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) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic + ! 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-phi0)/T_e)**gamma_e_dexp !Polytropic ! Check if the solution has converged phiConv = maxval(abs(Res),1)