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

Change order of array dimensions

Browse source
This commit is contained in:
JHendrikx 2025-02-04 11:15:04 +01:00
commit 726b11d718
1 changed files with 38 additions and 38 deletions
Download patch file Download diff file Expand all files Collapse all files

76
vlaplex.f90
View file

@ -59,7 +59,7 @@ program VlaPlEx
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(:):: f0 ! Boundary at r = x_0
real(dp), allocatable, dimension(:,:):: n_i
real(dp), allocatable, dimension(:,:):: u_i
real(dp), allocatable, dimension(:):: E_i
@ -158,13 +158,13 @@ program VlaPlEx
nz = 2
! Allocate vectors
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(f_i(1:nz,1:nr,1:nv), f_i_old(1:nz,1:nr,1:nv))
allocate(n_i(1:nz,1:nr))
allocate(u_i(1:nz,1:nr), E_i(1:nr), T_i(1:nz,1:nr))
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,1:nz))
allocate(fCum_i(1:nz,1:nv))
f_i = 0.0_dp
f_i_old = 0.0_dp
n_i = 0.0_dp
@ -217,7 +217,7 @@ program VlaPlEx
phi0 = 1.0e2_dp / phi_ref ! Dirichlet
phi(1) = phi0 ! Dirichlet
! phi0 = phi(1) ! Neumann
allocate(f0(j0:nv,1:nz))
allocate(f0(j0:nv))
f0 = 0.0_dp
! Output initial values
@ -228,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(:,1), u_i(:,1), T_i(:,1), Zave)
call writeOutputMom(t, dt, nr, r, n_i(1,:), u_i(1,:), T_i(1,:), Zave)
! Main loop
do t = 1, nt
@ -238,14 +238,14 @@ program VlaPlEx
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,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)
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
! 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
f_i_old(iz,1,j0:nv) = f0
f_i(iz,1,j0:nv) = f_i_old(iz,1,j0:nv)
T_i(iz,1) = Temp_bc
end do
T_e = Temp_bc
print *, 'Time: ', time * t_ref
@ -258,36 +258,36 @@ program VlaPlEx
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,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))
f_i(iz,i,1:j0-1) = f_i_old(iz,i,1:j0-1) - v(1:j0-1)*dt/dr/r(i)**2*(r(i+1)**2*f_i_old(iz,i+1,1:j0-1) - &
r(i )**2*f_i_old(iz,i ,1:j0-1))
end if
! Advect positive velocity
if (i > 1) then
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))
f_i(iz,i,j0:nv) = f_i_old(iz,i, j0:nv) - v( j0:nv)*dt/dr/r(i)**2*(r(i )**2*f_i_old(iz,i , j0:nv) - &
r(i-1)**2*f_i_old(iz,i-1, j0:nv))
end if
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
n_i(iz,i) = sum(f_i(iz,i,:))*dv
if (n_i(1,i) > 1.0e-10_dp) then
u_i(iz,i) = sum(v(:) *f_i(iz,i,:))*dv / n_i(iz,i)
E_i(i) = sum(v(:)**2*f_i(iz,i,:))*dv / n_i(iz,i)
T_i(iz,i) = 2.0_dp*E_i(i) - 2.0_dp*u_i(iz,i)**2
Zave(i) = Zave_bc
else
u_i(i,iz) = 0.0_dp
T_i(i,iz) = 0.0_dp
u_i(iz,i) = 0.0_dp
T_i(iz,i) = 0.0_dp
Zave(i) = 0.0_dp
end if
@ -296,7 +296,7 @@ program VlaPlEx
!$omp end parallel do
! Assume quasi-neutrality to start iterating
n_e = Zave * n_i(:,1)
n_e = Zave * n_i(1,:)
!do iz = 1, nz
! n_e = n_e + Zave * n_i(:,iz)
!end do
@ -322,7 +322,7 @@ program VlaPlEx
!do iz = 1, nz
! b = b - (Zave * n_i(:,iz))
!end do
b = -(Zave * n_i(:,1) - n_e)
b = -(Zave * n_i(1,:) - n_e)
! Apply boundary conditions
b(1) = phi0 ! Dirichlet
! b(nr) = 0.0_dp ! Dirichlet
@ -383,18 +383,18 @@ program VlaPlEx
! i = 1, v<0
i = 1
if (E(i) >= 0.0_dp) then
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))
f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,2:j0-2) - f_i_old(iz,i,1:j0-3))
else
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))
f_i(iz,i,2:j0-2) = f_i_old(iz,i,2:j0-2) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,3:j0-1) - f_i_old(iz,i,2:j0-2))
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,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))
f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,2:nv-1) - f_i_old(iz,i,1:nv-2))
else
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))
f_i(iz,i,2:nv-1) = f_i_old(iz,i,2:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,3:nv) - f_i_old(iz,i,2:nv-1))
end if
@ -403,9 +403,9 @@ program VlaPlEx
! i = nr, v>=0
i = nr
if (E(i) >= 0.0_dp) then
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))
f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,j0+1:nv-1) - f_i_old(iz,i,j0:nv-2))
else
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))
f_i(iz,i,j0+1:nv-1) = f_i_old(iz,i,j0+1:nv-1) - Zave(i)*E(i)*dt/dv*(f_i_old(iz,i,j0+2:nv) - f_i_old(iz,i,j0+1:nv-1))
end if
end do
@ -413,14 +413,14 @@ program VlaPlEx
! Reset values for next iteration
f_i_old = f_i
do iz = 1, nz
fCum_i(:,iz) = fCum_i(:,iz) + f_i_old(rCum_index,:,iz)
fCum_i(iz,:) = fCum_i(iz,:) + f_i_old(iz,rCum_index,:)
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(:,1), u_i(:,1), T_i(:,1), Zave)
call writeOutputFCum(t, dt, r(rCum_index), nv, v, fCum_i(:,1))
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