JGonzalez/vlaplex
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Trying to fix E at r0

Browse source
This commit is contained in:
Jorge Gonzalez 2024-11-12 10:05:44 +01:00
commit eb09a04baf
1 changed files with 15 additions and 14 deletions
Download patch file Download diff file Expand all files Collapse all files

29
plasmaExpansion.f90
View file

@ -98,9 +98,9 @@ program plasmaExpansion
call boundaryConditions%init(bc_file) call boundaryConditions%init(bc_file)
! Set domain boundaries (non-dimensional units) ! Set domain boundaries (non-dimensional units)
r0 = 200.0e-6_dp / L_ref r0 = 10.0e-6_dp / L_ref
rf = 6.0e-3_dp / L_ref rf = 50.0e-6_dp / L_ref
dr = 1.0e3_dp dr = 2.0e1_dp
nr = nint((rf - r0) / dr) + 1 nr = nint((rf - r0) / dr) + 1
dr = (rf - r0) / float(nr-1) dr = (rf - r0) / float(nr-1)
allocate(r(1:nr)) allocate(r(1:nr))
@ -109,7 +109,7 @@ program plasmaExpansion
end do end do
! Set position to calculate cumulative sum of f (non-dimensional units) ! 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 ! Index for cumulative sum
rCum_index = minloc(abs(r - rCum), 1) rCum_index = minloc(abs(r - rCum), 1)
@ -131,7 +131,7 @@ program plasmaExpansion
end if end if
t0 = 0.0_dp 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 ! tf = 1.0e1_dp * (rf - r0) / c_s
dt = 1.0e-2_dp*dr/c_s dt = 1.0e-2_dp*dr/c_s
nt = nint((tf - t0) / dt) nt = nint((tf - t0) / dt)
@ -183,8 +183,8 @@ program plasmaExpansion
b = 0.0_dp b = 0.0_dp
db_dphi = 0.0_dp db_dphi = 0.0_dp
diag = -2.0_dp / dr**2 diag = -2.0_dp / dr**2
diag_low = 1.0_dp / dr**2 - 1.0_dp / (r(2:nr) * 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) * 2.0_dp * dr) diag_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * dr)
diag(1) = 1.0_dp ! Dirichlet diag(1) = 1.0_dp ! Dirichlet
diag_high(1) = 0.0_dp ! Dirichlet diag_high(1) = 0.0_dp ! Dirichlet
! diag_high(1) = 2.0_dp / dr**2 ! Neumann ! diag_high(1) = 2.0_dp / dr**2 ! Neumann
@ -228,6 +228,7 @@ program plasmaExpansion
time = t * dt + t0 time = t * dt + t0
call boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_bc) call boundaryConditions%get(time, n_bc, u_bc, Temp_bc, Zave_bc)
call writeOutputBoundary(t, dt, 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) = 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(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) f0 = f0 * n_bc / (sum(f0)*dv)
@ -282,12 +283,12 @@ program plasmaExpansion
! Diagonal matrix for Newton integration scheme ! Diagonal matrix for Newton integration scheme
db_dphi = n_e / T_e ! Isotropic db_dphi = n_e / T_e ! Isotropic
diag = -2.0_dp / dr**2 - db_dphi 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_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) * 2.0_dp * dr) diag_high = 1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * dr)
diag(1) = 1.0_dp ! Dirichlet diag(1) = 1.0_dp ! Dirichlet
diag_high(1) = 0.0_dp ! Dirichlet diag_high(1) = 0.0_dp ! Dirichlet
! diag_high(1) = 2.0_dp / dr**2 - n_e(1) ! Neumann ! diag_high(1) = 2.0_dp / dr**2 - db_dphi(1) ! Neumann
! diag(nr) = 1.0_dp ! Dirichlet ! diag(nr) = 2.0_dp ! Dirichlet
! diag_low(nr-1) = 0.0_dp ! Dirichlet ! diag_low(nr-1) = 0.0_dp ! Dirichlet
diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann 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) b = -(Zave*n_i - n_e)
! Apply boundary conditions ! Apply boundary conditions
b(1) = phi0 ! Dirichlet b(1) = phi0 ! Dirichlet
! b(nr) = phiF ! Dirichlet ! b(nr) = 0.0_dp ! Dirichlet
! Calculate residual ! Calculate residual
!$omp parallel workshare !$omp parallel workshare
@ -305,10 +306,10 @@ program plasmaExpansion
! Iterate system ! Iterate system
call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info) call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info)
phi = phi_old + 1.0e-1_dp*Res phi = phi_old + 1.0e-1_dp*Res
! phi0=phi(1) phi0=phi(1) ! Neumann
! Calculate distribution of electrons ! 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) * 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 ! 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