diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
index e2a0b94a9f31087bb2ffe8b1173fd72963bab935..b43fa68dd42e6783c41cae69e3396b0f145fdbb7 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/multigrid.f90
@@ -473,7 +473,7 @@ contains
type(scalar3d), intent(inout) :: phicoarse
! local var
integer :: ix,iy,iz
- integer :: ix_min, ix_max, iy_min, iy_max, n
+ integer :: ix_min, ix_max, iy_min, iy_max, n ,nz
real(kind=rl) :: xn,xs,xw,xe
real , dimension(:,:,:) , pointer , contiguous :: zphifine_st , zphicoarse_st
@@ -483,13 +483,14 @@ contains
ix_max = phicoarse%icompx_max
iy_min = phicoarse%icompy_min
iy_max = phicoarse%icompy_max
+ nz = phicoarse%grid_param%nz
! three dimensional cell average
if (mg_param%restriction == REST_CELLAVERAGE) then
! Do not coarsen in z-direction
if (LUseO) then
do ix=ix_min,ix_max
do iy=iy_min,iy_max
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
phicoarse%s(iz,iy,ix) = &
phifine%s(iz ,2*iy ,2*ix ) + &
phifine%s(iz ,2*iy-1,2*ix ) + &
@@ -503,7 +504,7 @@ contains
zphifine_st => phifine%st
zphicoarse_st => phicoarse%st
!$acc kernels loop independent collapse(3)
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
do iy=iy_min,iy_max
do ix=ix_min,ix_max
zphicoarse_st(ix,iy,iz) = &
@@ -529,7 +530,7 @@ contains
xn=1.0
if (iy==1) xs=0.0
if (iy==n) xn=0.0
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
phicoarse%s(iz,iy,ix) = 0.25_rl * ( &
phifine%s(iz,2*iy+1,2*ix-1) * xn + &
phifine%s(iz,2*iy+1,2*ix ) * xn + &
@@ -548,7 +549,7 @@ contains
end do
end if
if (LUseT) then
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
do iy=iy_min,iy_max
xs=1.0
xn=1.0
@@ -587,19 +588,21 @@ contains
type(scalar3d), intent(in) :: phifine
type(scalar3d), intent(inout) :: phicoarse
integer :: ix,iy,iz
- integer :: ix_min, ix_max, iy_min, iy_max
+ integer :: ix_min, ix_max, iy_min, iy_max , nz
ix_min = phicoarse%icompx_min
ix_max = phicoarse%icompx_max
iy_min = phicoarse%icompy_min
iy_max = phicoarse%icompy_max
+ nz = phicoarse%grid_param%nz
+
! three dimensional cell average
if (mg_param%restriction == REST_CELLAVERAGE) then
! Do not coarsen in z-direction
if (LUseO) then
do ix=ix_min,ix_max
do iy=iy_min,iy_max
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
phicoarse%s(iz,iy,ix) = &
phifine%s(iz ,2*iy ,2*ix ) + &
phifine%s(iz ,2*iy-1,2*ix ) + &
@@ -610,7 +613,7 @@ contains
end do
end if
if (LUseT) then
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
do iy=iy_min,iy_max
do ix=ix_min,ix_max
phicoarse%st(ix,iy,iz) = &
@@ -755,14 +758,16 @@ contains
subroutine loop_over_grid_constant_mnh(iblock)
implicit none
integer, intent(in) :: iblock
- integer :: ix,iy,iz
+ integer :: ix,iy,iz, nz
+
+ nz = phicoarse%grid_param%nz
if (LUseO) then
do ix=ixmin(iblock),ixmax(iblock)
do iy=iymin(iblock),iymax(iblock)
do dix = -1,0
do diy = -1,0
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
phifine%s(iz,2*iy+diy,2*ix+dix) = phicoarse%s(iz,iy,ix)
end do
end do
@@ -775,7 +780,7 @@ contains
do iy=iymin(iblock),iymax(iblock)
do dix = -1,0
do diy = -1,0
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
phifine%st(2*ix+dix,2*iy+diy,iz) = phicoarse%st(ix,iy,iz)
end do
end do
@@ -797,6 +802,10 @@ contains
real , dimension(:,:,:) , pointer , contiguous :: zphifine_st , zphicoarse_st
+ integer :: nz
+
+ nz = phicoarse%grid_param%nz
+
! optimisation for newman MNH case : all coef constant
rhox = 0.25_rl
rhoy = 0.25_rl
@@ -805,7 +814,7 @@ contains
do ix=ixmin(iblock),ixmax(iblock)
do iy=iymin(iblock),iymax(iblock)
! Piecewise linear interpolation
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
do dix = -1,0
do diy = -1,0
phifine%s(iz,2*iy+diy,2*ix+dix) = &
@@ -827,7 +836,7 @@ contains
zphicoarse_st => phicoarse%st
!$acc kernels loop independent collapse(5)
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
do diy = -1,0
do dix = -1,0
do iy=iymin(iblock),iymax(iblock)
@@ -981,7 +990,7 @@ contains
do iy=iymin(iblock),iymax(iblock)
do dix = -1,0
do diy = -1,0
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
phifine%s(iz,2*iy+diy,2*ix+dix) = phicoarse%s(iz,iy,ix)
end do
end do
@@ -1001,7 +1010,7 @@ contains
do iy=iymin(iblock),iymax(iblock)
#ifdef PIECEWISELINEAR
! Piecewise linear interpolation
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
do dix = -1,0
do diy = -1,0
if ( (ix+(2*dix+1)+phicoarse%ix_min-1 < 1 ) .or. &
@@ -1105,7 +1114,7 @@ contains
end do
! invert A
call invertA(A)
- do iz=1,phicoarse%grid_param%nz
+ do iz=1,nz
! Calculate gradient on each level
dxu(1:2) = dx(1,1:2)*phicoarse%s(iz,iy ,ix-1) &
+ dx(2,1:2)*phicoarse%s(iz,iy ,ix+1) &
@@ -1498,6 +1507,7 @@ contains
real(kind=rl) :: alpha, beta, pq, rz, rz_old
integer :: ierr
real(kind=rl) , pointer , contiguous , dimension(:,:,:) :: z_one_st
+ integer :: icompx_max, icompy_max, nz
solvertype = solver_param%solvertype
resreduction = solver_param%resreduction
@@ -1524,22 +1534,27 @@ contains
!
! Init 1 vector = z_one
call create_scalar3d(MPI_COMM_HORIZ,bRHS%grid_param,halo_size,z_one)
+ !
+ icompx_max = z_one%icompx_max
+ icompy_max = z_one%icompy_max
+ nz = z_one%grid_param%nz
+ !
if (LUseO) then
z_one%s(:,:,:) = 0.0_rl
- z_one%s(1:z_one%grid_param%nz,1:z_one%icompy_max,1:z_one%icompx_max) = 1.0_rl
+ z_one%s(1:nz,1:icompy_max,1:icompx_max) = 1.0_rl
end if
if (LUseT) then
z_one_st => z_one%st
!$acc kernels
z_one_st(:,:,:) = 0.0_rl
- z_one_st(1:z_one%icompx_max,1:z_one%icompy_max,1:z_one%grid_param%nz) = 1.0_rl
+ z_one_st(1:icompx_max,1:icompy_max,1:nz) = 1.0_rl
!$acc end kernels
end if
! Mean / Norm of B
call scalarprod_mnh(pproc,z_one,z_one, mean_initial )
call scalarprod_mnh(pproc,z_one,bRHS, mean_initial )
norm_initial = l2norm_mnh(bRHS,.true.)
- mean_initial = mean_initial / (( z_one%grid_param%nz ) * ( z_one%grid_param%n )**2)
+ mean_initial = mean_initial / (( nz ) * ( z_one%grid_param%n )**2)
norm_initial = mean_initial / norm_initial
if (LMean) then
! b <- b -mean_initial * z_one
diff --git a/src/ZSOLVER/tensorproductmultigrid_Source/profiles.f90 b/src/ZSOLVER/tensorproductmultigrid_Source/profiles.f90
index 72338e703e1ab7a938de88d82c72d8b43c2ac763..70e7508d52411b9bd9e9c28b5dd81054d9f7a86f 100644
--- a/src/ZSOLVER/tensorproductmultigrid_Source/profiles.f90
+++ b/src/ZSOLVER/tensorproductmultigrid_Source/profiles.f90
@@ -58,7 +58,7 @@ private
type(grid_parameters), intent(in) :: grid_param
type(model_parameters), intent(in) :: model_param
type(scalar3d), intent(inout) :: u
- integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max
+ integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max, nz
real(kind=rl) :: x, y, z
real(kind=rl) :: rho, sigma, theta, phi, r, b_low, b_up, pi
@@ -71,12 +71,14 @@ private
ix_max = u%ix_max
iy_min = u%iy_min
iy_max = u%iy_max
+ nz = u%grid_param%nz
+
IF (.NOT. PRESENT(KIB) ) THEN
! Initialise RHS
if (LUseO) then
do ix=ix_min, ix_max
do iy=iy_min, iy_max
- do iz=1,u%grid_param%nz
+ do iz=1,nz
u%s(iz,iy-iy_min+1,ix-ix_min+1) = 0.0_rl
@@ -87,7 +89,7 @@ private
if (LUseT) then
zu_st => u%st
!$acc kernels loop independent collapse(3)
- do iz=1,u%grid_param%nz
+ do iz=1,nz
do iy=iy_min, iy_max
do ix=ix_min, ix_max
@@ -103,7 +105,7 @@ private
if (LUseO) then
do ix=ix_min, ix_max
do iy=iy_min, iy_max
- do iz=1,u%grid_param%nz
+ do iz=1,nz
u%s(iz,iy-iy_min+1,ix-ix_min+1) = PU(IX-ix_min+KIB,IY-iy_min+KJB,IZ)
end do
end do
@@ -112,7 +114,7 @@ private
if (LUseT) then
zu_st => u%st
!$acc kernels loop independent collapse(3)
- do iz=1,u%grid_param%nz
+ do iz=1,nz
do iy=iy_min, iy_max
do ix=ix_min, ix_max
zu_st(ix-ix_min+1,iy-iy_min+1,iz) = PU(IX-ix_min+KIB,IY-iy_min+KJB,IZ)
@@ -131,7 +133,7 @@ private
type(grid_parameters), intent(in) :: grid_param
type(model_parameters), intent(in) :: model_param
type(scalar3d), intent(inout) :: u
- integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max
+ integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max, nz
real(kind=rl) :: x, y, z
real(kind=rl) :: rho, sigma, theta, phi, r, b_low, b_up, pi
@@ -144,6 +146,8 @@ private
ix_max = u%ix_max
iy_min = u%iy_min
iy_max = u%iy_max
+ nz = u%grid_param%nz
+
IF (.NOT. PRESENT(KIB) ) THEN
!
ELSE
@@ -151,7 +155,7 @@ private
if (LUseO) then
do ix=ix_min, ix_max
do iy=iy_min, iy_max
- do iz=1,u%grid_param%nz
+ do iz=1,nz
PU(IX-ix_min+KIB,IY-iy_min+KJB,IZ) = u%s(iz,iy-iy_min+1,ix-ix_min+1)
end do
end do
@@ -159,7 +163,7 @@ private
else
zu_st => u%st
!$acc kernels loop independent collapse(3)
- do iz=1,u%grid_param%nz
+ do iz=1,nz
do iy=iy_min, iy_max
do ix=ix_min, ix_max
PU(IX-ix_min+KIB,IY-iy_min+KJB,IZ) = zu_st(ix-ix_min+1,iy-iy_min+1,iz)
@@ -178,7 +182,7 @@ private
type(grid_parameters), intent(in) :: grid_param
type(model_parameters), intent(in) :: model_param
type(scalar3d), intent(inout) :: b
- integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max
+ integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max, nz, n
real(kind=rl) :: x, y, z
real(kind=rl) :: rho, sigma, theta, phi, r, b_low, b_up, pi
@@ -191,15 +195,18 @@ private
ix_max = b%ix_max
iy_min = b%iy_min
iy_max = b%iy_max
+ nz = b%grid_param%nz
+ n = b%grid_param%n
+ !
IF (.NOT. PRESENT(KIB) ) THEN
! Initialise RHS
if (LUseO) then
do ix=ix_min, ix_max
do iy=iy_min, iy_max
- do iz=1,b%grid_param%nz
- x = 1.0_rl*((ix-0.5_rl)/(1.0_rl*b%grid_param%n))
- y = 1.0_rl*((iy-0.5_rl)/(1.0_rl*b%grid_param%n))
- z = 1.0_rl*((iz-0.5_rl)/(1.0_rl*b%grid_param%nz))
+ do iz=1,nz
+ x = 1.0_rl*((ix-0.5_rl)/(1.0_rl*n))
+ y = 1.0_rl*((iy-0.5_rl)/(1.0_rl*n))
+ z = 1.0_rl*((iz-0.5_rl)/(1.0_rl*nz))
b%s(iz,iy-iy_min+1,ix-ix_min+1) = 0.0_rl
@@ -219,12 +226,12 @@ private
if (LUseT) then
zb_st => b%st
!$acc kernels loop independent collapse(3)
- do iz=1,b%grid_param%nz
+ do iz=1,nz
do iy=iy_min, iy_max
do ix=ix_min, ix_max
- x = 1.0_rl*((ix-0.5_rl)/(1.0_rl*b%grid_param%n))
- y = 1.0_rl*((iy-0.5_rl)/(1.0_rl*b%grid_param%n))
- z = 1.0_rl*((iz-0.5_rl)/(1.0_rl*b%grid_param%nz))
+ x = 1.0_rl*((ix-0.5_rl)/(1.0_rl*n))
+ y = 1.0_rl*((iy-0.5_rl)/(1.0_rl*n))
+ z = 1.0_rl*((iz-0.5_rl)/(1.0_rl*nz))
zb_st(ix-ix_min+1,iy-iy_min+1,iz) = 0.0_rl
@@ -247,7 +254,7 @@ private
if (LUseO) then
do ix=ix_min, ix_max
do iy=iy_min, iy_max
- do iz=1,b%grid_param%nz
+ do iz=1,nz
b%s(iz,iy-iy_min+1,ix-ix_min+1) = PY(IX-ix_min+KIB,IY-iy_min+KJB,IZ)
end do
end do
@@ -256,7 +263,7 @@ private
if (LUseT) then
zb_st => b%st
!$acc kernels loop independent collapse(3)
- do iz=1,b%grid_param%nz
+ do iz=1,nz
do iy=iy_min, iy_max
do ix=ix_min, ix_max
zb_st(ix-ix_min+1,iy-iy_min+1,iz) = PY(IX-ix_min+KIB,IY-iy_min+KJB,IZ)
@@ -275,7 +282,7 @@ private
type(grid_parameters), intent(in) :: grid_param
type(model_parameters), intent(in) :: model_param
type(scalar3d), intent(inout) :: b
- integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max
+ integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max, nz, n
real(kind=rl) :: x, y, z
real(kind=rl) :: rho, sigma, theta, phi, r, b_low, b_up, pi
@@ -287,16 +294,18 @@ private
ix_max = b%ix_max
iy_min = b%iy_min
iy_max = b%iy_max
+ nz = b%grid_param%nz
+ n = b%grid_param%n
b_low = 1.0_rl+0.25*b%grid_param%H
b_up = 1.0_rl+0.75*b%grid_param%H
pi = 4.0_rl*atan2(1.0_rl,1.0_rl)
! Initialise RHS
do ix=ix_min, ix_max
do iy=iy_min, iy_max
- do iz=1,b%grid_param%nz
- x = 1.0_rl*((ix-0.5_rl)/(1.0_rl*b%grid_param%n))
- y = 1.0_rl*((iy-0.5_rl)/(1.0_rl*b%grid_param%n))
- z = 1.0_rl*((iz-0.5_rl)/(1.0_rl*b%grid_param%nz))
+ do iz=1,nz
+ x = 1.0_rl*((ix-0.5_rl)/(1.0_rl*n))
+ y = 1.0_rl*((iy-0.5_rl)/(1.0_rl*n))
+ z = 1.0_rl*((iz-0.5_rl)/(1.0_rl*nz))
#ifdef TESTCONVERGENCE
! RHS for analytical solution x*(1-x)*y*(1-y)*z*(1-z)
if (grid_param%vertbc == VERTBC_DIRICHLET) then
@@ -359,21 +368,23 @@ private
implicit none
type(grid_parameters), intent(in) :: grid_param
type(scalar3d), intent(inout) :: u
- integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max
+ integer :: ix, iy, iz, ix_min, ix_max, iy_min, iy_max, nz, n
real(kind=rl) :: x, y, z
ix_min = u%ix_min
ix_max = u%ix_max
iy_min = u%iy_min
iy_max = u%iy_max
+ nz = u%grid_param%nz
+ n = u%grid_param%n
! Initialise RHS
do ix=ix_min, ix_max
do iy=iy_min, iy_max
- do iz=1,u%grid_param%nz
- x = u%grid_param%L*((ix-0.5_rl)/(1.0_rl*u%grid_param%n))
- y = u%grid_param%L*((iy-0.5_rl)/(1.0_rl*u%grid_param%n))
- z = u%grid_param%H*((iz-0.5_rl)/(1.0_rl*u%grid_param%nz))
+ do iz=1,nz
+ x = u%grid_param%L*((ix-0.5_rl)/(1.0_rl*n))
+ y = u%grid_param%L*((iy-0.5_rl)/(1.0_rl*n))
+ z = u%grid_param%H*((iz-0.5_rl)/(1.0_rl*nz))
if (grid_param%vertbc == VERTBC_DIRICHLET) then
u%s(iz,iy-iy_min+1,ix-ix_min+1) &
= x*(1.0_rl-x) &