From eb09a04baf61f2a2195d9f5c7925fa97c46f7c0d Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 12 Nov 2024 10:05:44 +0100 Subject: [PATCH] Trying to fix E at r0 --- plasmaExpansion.f90 | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90 index ab18b34..a4f492b 100644 --- a/plasmaExpansion.f90 +++ b/plasmaExpansion.f90 @@ -98,9 +98,9 @@ program plasmaExpansion call boundaryConditions%init(bc_file) ! Set domain boundaries (non-dimensional units) - r0 = 200.0e-6_dp / L_ref - rf = 6.0e-3_dp / L_ref - dr = 1.0e3_dp + r0 = 10.0e-6_dp / L_ref + rf = 50.0e-6_dp / L_ref + dr = 2.0e1_dp nr = nint((rf - r0) / dr) + 1 dr = (rf - r0) / float(nr-1) allocate(r(1:nr)) @@ -109,7 +109,7 @@ program plasmaExpansion end do ! Set position to calculate cumulative sum of f (non-dimensional units) - rCum = 5.0e-3 / L_ref + rCum = 40.0e-6 / L_ref ! Index for cumulative sum rCum_index = minloc(abs(r - rCum), 1) @@ -131,7 +131,7 @@ program plasmaExpansion end if t0 = 0.0_dp - tf = 2.0e-6_dp / t_ref + tf = 3.0e-8_dp / t_ref ! tf = 1.0e1_dp * (rf - r0) / c_s dt = 1.0e-2_dp*dr/c_s nt = nint((tf - t0) / dt) @@ -183,8 +183,8 @@ program plasmaExpansion b = 0.0_dp db_dphi = 0.0_dp diag = -2.0_dp / dr**2 - 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_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) diag(1) = 1.0_dp ! Dirichlet diag_high(1) = 0.0_dp ! Dirichlet ! diag_high(1) = 2.0_dp / dr**2 ! Neumann @@ -228,6 +228,7 @@ program plasmaExpansion 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) ! 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) @@ -282,12 +283,12 @@ program plasmaExpansion ! Diagonal matrix for Newton integration scheme db_dphi = n_e / T_e ! Isotropic diag = -2.0_dp / dr**2 - db_dphi - 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_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) diag(1) = 1.0_dp ! Dirichlet diag_high(1) = 0.0_dp ! Dirichlet - ! diag_high(1) = 2.0_dp / dr**2 - n_e(1) ! Neumann - ! diag(nr) = 1.0_dp ! Dirichlet + ! diag_high(1) = 2.0_dp / dr**2 - db_dphi(1) ! Neumann + ! diag(nr) = 2.0_dp ! Dirichlet ! diag_low(nr-1) = 0.0_dp ! Dirichlet diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann @@ -295,7 +296,7 @@ program plasmaExpansion b = -(Zave*n_i - n_e) ! Apply boundary conditions b(1) = phi0 ! Dirichlet - ! b(nr) = phiF ! Dirichlet + ! b(nr) = 0.0_dp ! Dirichlet ! Calculate residual !$omp parallel workshare @@ -305,10 +306,10 @@ program plasmaExpansion ! Iterate system call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info) phi = phi_old + 1.0e-1_dp*Res - ! phi0=phi(1) + phi0=phi(1) ! Neumann ! Calculate distribution of electrons - T_e = 1.0_dp + T_e = Temp_bc 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_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working