Some updates

This commit is contained in:
Jorge Gonzalez 2024-10-01 18:01:27 +02:00
commit 7ffdf8e65c

View file

@ -349,16 +349,13 @@ program plasmaExpansion
L_ref = u_ref * t_ref L_ref = u_ref * t_ref
phi_ref = kb * Temp_ref / qe phi_ref = kb * Temp_ref / qe
! Set position to calculate cumulative sum of f (non-dimensional units)
rCum = 5.0e-3 / L_ref
! Set input parameters (remember these have to be in non-dimensional units) ! Set input parameters (remember these have to be in non-dimensional units)
Temp0 = 60.0_dp * eV_to_K / Temp_ref Temp0 = 60.0_dp * eV_to_K / Temp_ref
TempF = 60.0_dp * eV_to_K / Temp_ref TempF = 10.0_dp * eV_to_K / Temp_ref
n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref n_ecr = 1.0e19_dp * cm3_to_m3 / n_ref
c_s = sqrt(T_to_Z(Temp0) * gam * Temp0) c_s = sqrt(T_to_Z(Temp0) * gam * Temp0)
u_bc0 = c_s!sqrt(Temp0) u_bc0 = c_s!sqrt(Temp0)
u_bcF = u_bc0!sqrt(TempF) u_bcF = sqrt(TempF)
n_bc0 = n_ecr / T_to_Z(Temp0) n_bc0 = n_ecr / T_to_Z(Temp0)
n_bcF = n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0) n_bcF = n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0)
@ -373,6 +370,9 @@ program plasmaExpansion
r(i) = dr * float(i-1) + r0 r(i) = dr * float(i-1) + r0
end do end do
! Set position to calculate cumulative sum of f (non-dimensional units)
rCum = 5.0e-3 / L_ref
! Index for cumulative sum ! Index for cumulative sum
rCum_index = minloc(abs(r - rCum), 1) rCum_index = minloc(abs(r - rCum), 1)
@ -395,7 +395,7 @@ program plasmaExpansion
t0 = 0.0_dp t0 = 0.0_dp
tf = 1.0e-6_dp / t_ref tf = 1.0e-6_dp / t_ref
! tf = 1.0e1_dp * (rf - r0) / c_s ! tf = 1.0e1_dp * (rf - r0) / c_s
dt = 1.0e-1_dp*dr/(10.0_dp*c_s) dt = 1.0e-2_dp*dr/(10.0_dp*c_s)
nt = nint((tf - t0) / dt) nt = nint((tf - t0) / dt)
dt = (tf - t0) / float(nt) dt = (tf - t0) / float(nt)
@ -500,7 +500,8 @@ program plasmaExpansion
n_bc = n_bcF n_bc = n_bcF
end if 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, T_to_Z(Temp_bc))
f0(j0:nv) = n_bc / sqrt(PI*Temp_bc) * 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 = f0 * n_bc / (sum(f0)*dv)
! Boundary conditions ! Boundary conditions
! r = r0, v>0 ! r = r0, v>0
@ -528,10 +529,10 @@ program plasmaExpansion
r(i-1)**2*f_i_old(i-1, j0:nv)) r(i-1)**2*f_i_old(i-1, j0:nv))
end if end if
n_i(i) = sum(f_i(i,:)*dv) n_i(i) = sum(f_i(i,:))*dv
if (n_i(i) > 1.0e-10_dp) then if (n_i(i) > 1.0e-10_dp) then
u_i(i) = sum(v(:) *f_i(i,:)*dv) / n_i(i) 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) 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 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) = T_to_Z(T_i(i))