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

Add additional nZ rank to arrays

Browse source
This commit is contained in:
JHendrikx 2025-02-03 20:56:04 +01:00
commit 3754c0b910
1 changed files with 98 additions and 80 deletions
Download patch file Download diff file Expand all files Collapse all files

120
vlaplex.f90
View file

@ -46,8 +46,8 @@ program VlaPlEx
real(dp):: t0, tf
real(dp):: time
real(dp):: dr, dv, dt
integer:: nr, nv, nt
integer:: i, j, t
integer:: nr, nv, nt, nz
integer:: i, iz, j, t
integer:: j0 ! First integer of positive velocity
real(dp):: Temp_bc ! Temperature
@ -58,12 +58,12 @@ program VlaPlEx
type(tableBC):: boundaryConditions
character(:), allocatable:: bc_file
real(dp), allocatable, dimension(:,:):: f_i, f_i_old
real(dp), allocatable, dimension(:):: f0 ! Boundary at r = x_0
real(dp), allocatable, dimension(:):: n_i
real(dp), allocatable, dimension(:):: u_i
real(dp), allocatable, dimension(:,:,:):: f_i, f_i_old
real(dp), allocatable, dimension(:,:):: f0 ! Boundary at r = x_0
real(dp), allocatable, dimension(:,:):: n_i
real(dp), allocatable, dimension(:,:):: u_i
real(dp), allocatable, dimension(:):: E_i
real(dp), allocatable, dimension(:):: T_i
real(dp), allocatable, dimension(:,:):: T_i
real(dp), allocatable, dimension(:):: n_e
real(dp), allocatable, dimension(:):: Zave
real(dp), allocatable, dimension(:):: diag, diag_low, diag_high
@ -79,7 +79,7 @@ program VlaPlEx
! real(dp):: phiF
integer:: k
real(dp), allocatable, dimension(:):: fCum_i
real(dp), allocatable, dimension(:,:):: fCum_i
real(dp):: rCum
integer:: rCum_index
@ -96,7 +96,7 @@ program VlaPlEx
! Set input parameters (remember these have to be in non-dimensional units)
c_s = sqrt(11.0_dp * gamma_i * 1.0_dp)
bc_file = 'bc_80ns_T60Z16.csv'
bc_file = 'bc.csv'
call boundaryConditions%init(bc_file)
! Set domain boundaries (non-dimensional units)
@ -156,14 +156,15 @@ program VlaPlEx
write(*, '(A,ES0.4e3)') 'CFL: ', dt*vf/dr
nz = 2
! Allocate vectors
allocate(f_i(1:nr,1:nv), f_i_old(1:nr,1:nv))
allocate(n_i(1:nr))
allocate(u_i(1:nr), E_i(1:nr), T_i(1:nr))
allocate(f_i(1:nr,1:nv,1:nz), f_i_old(1:nr,1:nv,1:nz))
allocate(n_i(1:nr,1:nz))
allocate(u_i(1:nr,1:nz), E_i(1:nr), T_i(1:nr,1:nz))
allocate(Zave(1:nr))
allocate(n_e(1:nr))
allocate(phi(1:nr), phi_old(1:nr), E(1:nr))
allocate(fCum_i(1:nv))
allocate(fCum_i(1:nv,1:nz))
f_i = 0.0_dp
f_i_old = 0.0_dp
n_i = 0.0_dp
@ -216,7 +217,7 @@ program VlaPlEx
phi0 = 1.0e2_dp / phi_ref ! Dirichlet
phi(1) = phi0 ! Dirichlet
! phi0 = phi(1) ! Neumann
allocate(f0(j0:nv))
allocate(f0(j0:nv,1:nz))
f0 = 0.0_dp
! Output initial values
@ -227,7 +228,7 @@ program VlaPlEx
! call writeOutputF(t, dt, nr, r, nv, v, f_i_old)
call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i)
call writeOutputPhi(t, dt, nr, r, phi, E, n_e)
call writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Zave)
call writeOutputMom(t, dt, nr, r, n_i(:,1), u_i(:,1), T_i(:,1), Zave)
! Main loop
do t = 1, nt
@ -235,9 +236,17 @@ program VlaPlEx
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)
do iz = 1, nz
! 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)
f0(j0:nv,iz) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
f0(:,iz) = f0(:,iz) * n_bc / (sum(f0(:,iz))*dv)
! Boundary conditions
! r = r0, v>0
f_i_old(1,j0:nv,iz) = f0(:,iz)
f_i(1,j0:nv,iz) = f_i_old(1,j0:nv,iz)
T_i(1,iz) = Temp_bc
end do
T_e = Temp_bc
print *, 'Time: ', time * t_ref
print *, 'Temp_bc: ', Temp_bc
@ -246,50 +255,51 @@ program VlaPlEx
print *, '-------------------------'
! Boundary conditions
! r = r0, v>0
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) = 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)
f_i_old(nr,1:j0-1,:) = 0.0_dp
f_i(nr,1:j0-1,:) = f_i_old(nr,1:j0-1,:)
! set edge velocities to 0
f_i_old(:,1) = 0.0_dp
f_i_old(:,nv) = 0.0_dp
f_i_old(:,1,:) = 0.0_dp
f_i_old(:,nv,:) = 0.0_dp
! Advect in the r direction
!$omp parallel do
do iz = 1, nz
do i = 1, nr
! Advect negative velocity
if (i < nr) then
f_i(i,1:j0-1) = f_i_old(i,1:j0-1) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(i+1,1:j0-1) - &
r(i )**2*f_i_old(i ,1:j0-1))
f_i(i,1:j0-1,iz) = f_i_old(i,1:j0-1,iz) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(i+1,1:j0-1,iz) - &
r(i )**2*f_i_old(i ,1:j0-1,iz))
end if
! Advect positive velocity
if (i > 1) then
f_i(i,j0:nv) = f_i_old(i, j0:nv) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(i , j0:nv) - &
r(i-1)**2*f_i_old(i-1, j0:nv))
f_i(i,j0:nv,iz) = f_i_old(i, j0:nv,iz) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(i , j0:nv,iz) - &
r(i-1)**2*f_i_old(i-1, j0:nv,iz))
end if
n_i(i) = sum(f_i(i,:))*dv
if (n_i(i) > 1.0e-10_dp) then
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
n_i(i,iz) = sum(f_i(i,:,iz))*dv
if (n_i(i,1) > 1.0e-10_dp) then
u_i(i,iz) = sum(v(:) *f_i(i,:,iz))*dv / n_i(i,iz)
E_i(i) = sum(v(:)**2*f_i(i,:,iz))*dv / n_i(i,iz)
T_i(i,iz) = 2.0_dp*E_i(i) - 2.0_dp*u_i(i,iz)**2
Zave(i) = Zave_bc
else
u_i(i) = 0.0_dp
T_i(i) = 0.0_dp
u_i(i,iz) = 0.0_dp
T_i(i,iz) = 0.0_dp
Zave(i) = 0.0_dp
end if
end do
end do
!$omp end parallel do
! Assume quasi-neutrality to start iterating
n_e = Zave * n_i
n_e = Zave * n_i(:,1)
!do iz = 1, nz
! n_e = n_e + Zave * n_i(:,iz)
!end do
db_dphi = 0.0_dp
! Solve Poission (maximum number of iterations, break if convergence is reached before)
@ -307,7 +317,12 @@ program VlaPlEx
diag_low(nr-1) = 2.0_dp / dr**2 - db_dphi(nr) ! Neumann
! Calculate charge density
b = -(Zave*n_i - n_e)
!b = n_e
!b = b - (Zave * n_i(:,1))
!do iz = 1, nz
! b = b - (Zave * n_i(:,iz))
!end do
b = -(Zave * n_i(:,1) - n_e)
! Apply boundary conditions
b(1) = phi0 ! Dirichlet
! b(nr) = 0.0_dp ! Dirichlet
@ -324,10 +339,10 @@ program VlaPlEx
! Calculate distribution of electrons
! n_e = Zave(1) * n_i(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
n_e = Zave(1) * n_i(1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic
n_e = Zave(1) * n_i(1,1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_exp !Polytropic
! Diagonal matrix for Newton integration scheme
! db_dphi = n_e / T_e ! Isotropic
db_dphi = Zave(1) * n_i(1) / (gamma_e * T_e) * &
db_dphi = Zave(1) * n_i(1,1) / (gamma_e * T_e) * &
(1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**gamma_e_dexp !Polytropic
! Check if the solution has converged
@ -363,22 +378,23 @@ program VlaPlEx
! Update intermediate f
f_i_old = f_i
do iz = 1, nz
! Advect in the v direction
! i = 1, v<0
i = 1
if (E(i) >= 0.0_dp) then
f_i(i,2:j0-2) = f_i_old(i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:j0-2) - f_i_old(i,1:j0-3))
f_i(i,2:j0-2,iz) = f_i_old(i,2:j0-2,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:j0-2,iz) - f_i_old(i,1:j0-3,iz))
else
f_i(i,2:j0-2) = f_i_old(i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:j0-1) - f_i_old(i,2:j0-2))
f_i(i,2:j0-2,iz) = f_i_old(i,2:j0-2,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:j0-1,iz) - f_i_old(i,2:j0-2,iz))
end if
! i = 2, nr-1; all v
!$omp parallel do
do i = 2, nr-1
if (E(i) >= 0.0_dp) then
f_i(i,2:nv-1) = f_i_old(i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:nv-1) - f_i_old(i,1:nv-2))
f_i(i,2:nv-1,iz) = f_i_old(i,2:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,2:nv-1,iz) - f_i_old(i,1:nv-2,iz))
else
f_i(i,2:nv-1) = f_i_old(i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:nv) - f_i_old(i,2:nv-1))
f_i(i,2:nv-1,iz) = f_i_old(i,2:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,3:nv,iz) - f_i_old(i,2:nv-1,iz))
end if
@ -387,22 +403,24 @@ program VlaPlEx
! i = nr, v>=0
i = nr
if (E(i) >= 0.0_dp) then
f_i(i,j0+1:nv-1) = f_i_old(i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+1:nv-1) - f_i_old(i,j0:nv-2))
f_i(i,j0+1:nv-1,iz) = f_i_old(i,j0+1:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+1:nv-1,iz) - f_i_old(i,j0:nv-2,iz))
else
f_i(i,j0+1:nv-1) = f_i_old(i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+2:nv) - f_i_old(i,j0+1:nv-1))
f_i(i,j0+1:nv-1,iz) = f_i_old(i,j0+1:nv-1,iz) - Zave(i)*E(i)*dt/dv*(f_i_old(i,j0+2:nv,iz) - f_i_old(i,j0+1:nv-1,iz))
end if
end do
! Reset values for next iteration
f_i_old = f_i
fCum_i = fCum_i + f_i_old(rCum_index,:)
do iz = 1, nz
fCum_i(:,iz) = fCum_i(:,iz) + f_i_old(rCum_index,:,iz)
end do
! Write output
if (mod(t,everyOutput) == 0 .or. t == nt) then
! call writeOutputF(t, dt, nr, r, nv, v, f_i_old)
call writeOutputPhi(t, dt, nr, r, phi, E, n_e)
call writeOutputMom(t, dt, nr, r, n_i, u_i, T_i, Zave)
call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i)
call writeOutputMom(t, dt, nr, r, n_i(:,1), u_i(:,1), T_i(:,1), Zave)
call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i(:,1))
end if