From 47fe2fb55dde1e1c209c02e7af7cb01835f734a6 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 16 Oct 2024 17:56:36 +0200 Subject: [PATCH] Preparing for polytropic electrons --- plasmaExpansion.f90 | 37 +++++++++---------------------------- 1 file changed, 9 insertions(+), 28 deletions(-) diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90 index 45354ed..ee875a5 100644 --- a/plasmaExpansion.f90 +++ b/plasmaExpansion.f90 @@ -52,7 +52,6 @@ program plasmaExpansion real(dp):: Zave_bc ! Average charge state real(dp):: u_bc ! Injection velocity real(dp):: n_bc ! Injection density - ! real(dp):: n_ecr ! Electron critical density for the laser real(dp):: c_s ! Ion sound speed type(tableBC):: boundaryConditions character(:), allocatable:: bc_file @@ -74,6 +73,7 @@ program plasmaExpansion real(dp), allocatable, dimension(:):: phi, phi_old, E real(dp):: phiConv real(dp):: phi0 + real(dp):: T_e ! real(dp):: phiF integer:: k @@ -93,16 +93,7 @@ program plasmaExpansion phi_ref = kb * Temp_ref / qe ! Set input parameters (remember these have to be in non-dimensional units) - ! Temp0 = 30.0_dp * eV_to_K / Temp_ref - ! TempF = 5.0_dp * eV_to_K / Temp_ref - ! Zave0 = 9.0_dp - ! ZaveF = 9.0_dp - ! n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref c_s = sqrt(9.0_dp * gamma_i * 1.0_dp) - ! u_bc0 = sqrt(Temp0) - ! u_bcF = sqrt(TempF) - ! n_bc0 = n_ecr / Zave0 - ! n_bcF = n_bc0*1.0e-2_dp!n_ecr*1.0e-1 / T_to_Z(Temp0) bc_file = 'bc.csv' call boundaryConditions%init(bc_file) @@ -146,9 +137,6 @@ program plasmaExpansion nt = nint((tf - t0) / dt) dt = (tf - t0) / float(nt) - ! t_bc0 = nint(100.0e-9_dp / t_ref / dt) - ! t_bcF = nint(150.0e-9_dp / t_ref / dt) - everyOutput = nint(1.0e-9_dp/t_ref/dt) if (everyOutput == 0) then everyOutput = 1 @@ -178,6 +166,7 @@ program plasmaExpansion E_i = 0.0_dp T_i = 0.0_dp n_e = 1.0_dp + T_e = 1.0_dp Zave = 0.0_dp phi = 0.0_dp phi_old = 0.0_dp @@ -235,20 +224,10 @@ program plasmaExpansion ! Main loop do t = 1, nt time = t * dt + t0 - ! 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 boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_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(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) ! Boundary conditions @@ -299,7 +278,7 @@ program plasmaExpansion phi_old = phi ! Diagonal matrix for Newton integration scheme - diag = -2.0_dp / dr**2 - n_e + diag = -2.0_dp / dr**2 - n_e / T_e 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(1) = 1.0_dp ! Dirichlet @@ -307,7 +286,7 @@ program plasmaExpansion ! diag_high(1) = 2.0_dp / dr**2 - n_e(1) ! Neumann ! diag(nr) = 1.0_dp ! Dirichlet ! diag_low(nr-1) = 0.0_dp ! Dirichlet - diag_low(nr-1) = 2.0_dp / dr**2 - n_e(nr) ! Neumann + diag_low(nr-1) = 2.0_dp / dr**2 - n_e(nr) / T_e ! Neumann ! Calculate charge density b = -(Zave*n_i - n_e) @@ -326,7 +305,9 @@ program plasmaExpansion ! phi0=phi(1) ! Calculate distribution of electrons - n_e = Zave(1) * n_i(1) * exp((phi - phi0) / 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) * (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 ! Check if the solution has converged phiConv = maxval(abs(Res),1)