From 0cfbdd2d0785580265edf4b9089d0da45918b786 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 1 Oct 2024 21:00:44 +0200 Subject: [PATCH] Case with Diko's peak So I'd to make some changes to the Newton iterative method, but it's working now. It's not giving noise, it converges, and with these conditions the case reproduces the Diko's peak. --- plasmaExpansion.f90 | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90 index d146046..2a4a331 100644 --- a/plasmaExpansion.f90 +++ b/plasmaExpansion.f90 @@ -307,6 +307,8 @@ program plasmaExpansion real(dp):: Temp_bc ! Temperature real(dp):: Temp0, TempF + real(dp):: Zave_bc ! Average charge state + real(dp):: Zave0, ZaveF real(dp):: n_ecr ! Electron critical density for the laser real(dp):: c_s ! Ion sound speed real(dp):: u_bc ! Injection velocity @@ -352,12 +354,14 @@ program plasmaExpansion ! Set input parameters (remember these have to be in non-dimensional units) Temp0 = 60.0_dp * eV_to_K / Temp_ref TempF = 10.0_dp * eV_to_K / Temp_ref + Zave0 = 12.0_dp + ZaveF = 3.0_dp n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref - c_s = sqrt(T_to_Z(Temp0) * gam * Temp0) + c_s = sqrt(Zave0 * gam * Temp0) u_bc0 = c_s!sqrt(Temp0) u_bcF = sqrt(TempF) - n_bc0 = n_ecr / T_to_Z(Temp0) - n_bcF = n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0) + n_bc0 = n_ecr / Zave0 + n_bcF = 1.0e-1*n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0) ! Set domain boundaries (non-dimensional units) r0 = 200.0e-6_dp / L_ref @@ -487,19 +491,22 @@ program plasmaExpansion ! Main loop Temp_bc = Temp0 + Zave_bc = Zave0 u_bc = u_bc0 n_bc = n_bc0 do t = 1, nt if (t > t_bc0 .and. t <= t_bcF) then Temp_bc = (TempF - Temp0) / float(t_bcF - t_bc0)*(t - t_bc0) + Temp0 + Zave_bc = (ZaveF - Zave0) / float(t_bcF - t_bc0)*(t - t_bc0) + Zave0 u_bc = (u_bcF - u_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + u_bc0 n_bc = (n_bcF - n_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + n_bc0 else if (t > t_bcF) then Temp_bc = TempF + Zave_bc = ZaveF u_bc = u_bcF n_bc = n_bcF end if - call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, T_to_Z(Temp_bc)) + call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc) f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) f0 = f0 * n_bc / (sum(f0)*dv) @@ -508,7 +515,7 @@ program plasmaExpansion f_i_old(1,j0:nv) = f0 f_i(1,j0:nv) = f_i_old(1,j0:nv) T_i(1) = Temp_bc - Zave(1) = T_to_Z(Temp_bc) + Zave(1) = Zave_bc ! r = rf, v<0 f_i_old(nr,1:j0-1) = 0.0_dp f_i(nr,1:j0-1) = f_i_old(nr,1:j0-1) @@ -534,7 +541,7 @@ program plasmaExpansion u_i(i) = sum(v(:) *f_i(i,:))*dv / n_i(i) E_i(i) = sum(v(:)**2*f_i(i,:))*dv / n_i(i) T_i(i) = 2.0_dp*E_i(i) - 2.0_dp*u_i(i)**2 - Zave(i) = T_to_Z(T_i(i)) + Zave(i) = Zave_bc else u_i(i) = 0.0_dp @@ -574,7 +581,7 @@ program plasmaExpansion ! Iterate system call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info) - phi = phi_old + Res + phi = phi_old + 1.0e-1_dp*Res ! phi0=phi(1) ! Calculate distribution of electrons