!!******************************************************************************
!!
!! This file is part of the AMUN source code, a program to perform
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!! adaptive mesh.
!!
!! Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: COORDINATES
!!
!! This module provides variables and subroutines handling the coordinates
!! for all refinement levels.
!!
!!******************************************************************************
!
module coordinates
! module variables are not implicit by default
!
implicit none
! MODULE PARAMETERS:
! =================
!
! the number of interior cells along each block dimension (and half of it)
!
integer , save :: ncells = 8, ncells_half = 4
! the number of ghost zones at each side of the block (and half and double of it)
!
integer , save :: nghosts = 2, nghosts_half = 1, nghosts_double = 4
! the number of all block cells along each dimension (and half of it)
!
integer , save :: bcells = 12, bcells_half = 6
! the domain division - the number of blocks in each direction at the base level
!
integer, dimension(3), save :: domain_base_dims = 1
! the limits of refinement level
!
integer , save :: minlev = 1, maxlev = 1, toplev = 1
! block interior and ghost zone indices
!
integer , save :: nbl = 2, nb = 3, nbu = 4, nbm = 2, nbp = 4
integer , save :: nel = 9, ne = 10, neu = 11, nem = 9, nep = 11
! the domain bounds
!
real(kind=8), save :: xmin = 0.0d+00
real(kind=8), save :: xmax = 1.0d+00
real(kind=8), save :: xlen = 1.0d+00
real(kind=8), save :: ymin = 0.0d+00
real(kind=8), save :: ymax = 1.0d+00
real(kind=8), save :: ylen = 1.0d+00
real(kind=8), save :: zmin = 0.0d+00
real(kind=8), save :: zmax = 1.0d+00
real(kind=8), save :: zlen = 1.0d+00
! the domain volume and its inversion
!
real(kind=8), save :: vol = 1.0d+00
real(kind=8), save :: voli = 1.0d+00
! the domain boundary areas
!
real(kind=8), save :: xarea = 1.0d+00
real(kind=8), save :: yarea = 1.0d+00
real(kind=8), save :: zarea = 1.0d+00
! flags to indicate periodicity of boundaries
!
logical, dimension(3), save :: periodic = .true.
! the block coordinates for all levels of refinement
!
real(kind=8), dimension(:,:), allocatable, save :: ax , ay , az
real(kind=8), dimension(: ), allocatable, save :: adx , ady , adz, adr
real(kind=8), dimension(: ), allocatable, save :: adxi, adyi, adzi
real(kind=8), dimension(: ), allocatable, save :: advol
! define type for rectangular subarray description
!
type rectangular
integer, dimension(NDIMS) :: l ! indices of the lower corner
integer, dimension(NDIMS) :: u ! indices of the upper corner
end type rectangular
! the subarray indices to ghost and domain areas used for boundary exchange
! ('c' for copy, 'p' for prolongation, 'r' for restriction)
!
#if NDIMS == 2
type(rectangular), dimension(2,2,NDIMS) , save :: edges_dc , edges_gc
type(rectangular), dimension(2,2,NDIMS) , save :: edges_dp , edges_gp
type(rectangular), dimension(2,2,NDIMS) , save :: edges_dr , edges_gr
type(rectangular), dimension(2,2) , save :: corners_dc, corners_gc
type(rectangular), dimension(2,2) , save :: corners_dp, corners_gp
type(rectangular), dimension(2,2) , save :: corners_dr, corners_gr
#endif /* NDIMS == 2 */
#if NDIMS == 3
type(rectangular), dimension(2,2,2,NDIMS), save :: faces_dc , faces_gc
type(rectangular), dimension(2,2,2,NDIMS), save :: faces_dp , faces_gp
type(rectangular), dimension(2,2,2,NDIMS), save :: faces_dr , faces_gr
type(rectangular), dimension(2,2,2,NDIMS), save :: edges_dc , edges_gc
type(rectangular), dimension(2,2,2,NDIMS), save :: edges_dp , edges_gp
type(rectangular), dimension(2,2,2,NDIMS), save :: edges_dr , edges_gr
type(rectangular), dimension(2,2,2) , save :: corners_dc, corners_gc
type(rectangular), dimension(2,2,2) , save :: corners_dp, corners_gp
type(rectangular), dimension(2,2,2) , save :: corners_dr, corners_gr
#endif /* NDIMS == 3 */
! by default everything is private
!
public
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
contains
!
!===============================================================================
!
! subroutine INITIALIZE_COORDINATES:
! ---------------------------------
!
! Subroutine initializes mesh coordinates and other coordinate parameters.
!
! Arguments:
!
! ni - the number of interior cells per dimension in a block;
! ng - the number of ghosts cells at each side in a block;
! tlevel - the maximum refinement level for coordinate generation;
! bdims - the base level dimension in blocks;
! xmn, xmx, ymn, ymx, zmn, zmx -
! the domain extrema in each direction;
! verbose - flag determining if the subroutine should be verbose;
! status - return flag of the procedure execution status;
!
!===============================================================================
!
subroutine initialize_coordinates(ni, ng, tlevel, bdims, &
xmn, xmx, ymn, ymx, zmn, zmx, verbose, status)
use helpers , only : print_message
use parameters, only : get_parameter
implicit none
logical , intent(in) :: verbose
integer , intent(in) :: ni, ng, tlevel
integer, dimension(3), intent(in) :: bdims
real(kind=8) , intent(in) :: xmn, xmx, ymn, ymx, zmn, zmx
integer , intent(out) :: status
character(len=32) :: lbndry, ubndry
integer :: l, i, j, p, q
integer :: fi, fj
#if NDIMS == 3
integer :: k, r, fk
#endif /* NDIMS == 3 */
integer :: nh, nd, nn, nm, np, nr, nq, ns, nt, nu
character(len=80) :: msg
character(len=*), parameter :: loc = 'COORDINATES::initialize_coordinates()'
!
!-------------------------------------------------------------------------------
!
status = 0
! check if the number of cells and ghost zones are acceptable
!
if (ni < 4 .or. ng < 2 .or. ni < 2 * ng &
.or. mod(ni,2) == 1 .or. mod(ng,2) == 1) then
if (verbose) then
call print_message(loc, "Wrong number of block cells or ghost zones:")
if (ni < 4) then
write(msg,"(5x,a,i0,a)") "=> ncells = ", ni, " < 4."
call print_message(loc, msg)
end if
if (ng < 2) then
write(msg,"(5x,a,i0,a)") "=> nghosts = ", ng, " < 2."
call print_message(loc, msg)
end if
if (ni < 2 * ng) then
write(msg,"(5x,a,2(i0,a))") "=> ncells = ", ni, &
" < 2 * nghosts = 2 * ", ng, "."
call print_message(loc, msg)
end if
if (mod(ni,2) == 1) then
write(msg,"(5x,a,i0,a)") "=> ncells = ", ni, " is an odd number."
call print_message(loc, msg)
end if
if (mod(ng,2) == 1) then
write(msg,"(5x,a,i0,a)") "=> nghosts = ", ng, " is an odd number."
call print_message(loc, msg)
end if
end if
status = 1
return
end if
! check if the domain base dimensions are correct
!
if (minval(bdims) <= 0) then
if (verbose) then
call print_message(loc, "Wrong domain base dimensions:")
if (bdims(1) < 1) then
write(msg,"(5x,a,i0,a)") "=> xblocks = ", bdims(1), " < 1."
call print_message(loc, msg)
end if
if (bdims(2) < 1) then
write(msg,"(5x,a,i0,a)") "=> yblocks = ", bdims(2), " < 1."
call print_message(loc, msg)
end if
#if NDIMS == 3
if (bdims(3) < 1) then
write(msg,"(5x,a,i0,a)") "=> zblocks = ", bdims(3), " < 1."
call print_message(loc, msg)
end if
#endif /* NDIMS == 3 */
end if
status = 1
return
end if
! check if the domain limits are correct
!
if (xmn >= xmx .or. ymn >= ymx .or. zmn >= zmx) then
if (verbose) then
call print_message(loc, "Wrong domain limits:")
if (xmn >= xmx) then
write(msg,"(5x,a,2(es12.5,a))") "=> xmin = ", xmn, &
" >= xmax = ", xmx, "."
call print_message(loc, msg)
end if
if (ymn >= ymx) then
write(msg,"(5x,a,2(es12.5,a))") "=> ymin = ", ymn, &
" >= ymax = ", ymx, "."
call print_message(loc, msg)
end if
#if NDIMS == 3
if (zmn >= zmx) then
write(msg,"(5x,a,2(es12.5,a))") "=> zmin = ", zmn, &
" >= zmax = ", zmx, "."
call print_message(loc, msg)
end if
#endif /* NDIMS == 3 */
end if
status = 1
return
end if
! set the number of cells along each block dimension and
! the number of ghost zones
!
ncells = ni
ncells_half = ni / 2
nghosts = ng
nghosts_half = ng / 2
nghosts_double = ng * 2
bcells = ni + 2 * ng
bcells_half = ncells_half + ng
! local variables
!
nh = nghosts_half
nd = nghosts_double
! calculate block indices
!
nb = ng + 1
ne = ng + ni
nbl = nb - 1
nbu = nb + ng - 1
nel = ne - ng + 1
neu = ne + 1
nbm = nb - 1
nem = ne - 1
nbp = nb + 1
nep = ne + 1
! determine domain periodicity
!
lbndry = "periodic"
ubndry = "periodic"
call get_parameter("xlbndry", lbndry)
call get_parameter("xubndry", ubndry)
periodic(1) = lbndry == "periodic" .and. ubndry == "periodic"
lbndry = "periodic"
ubndry = "periodic"
call get_parameter("ylbndry", lbndry)
call get_parameter("yubndry", ubndry)
periodic(2) = lbndry == "periodic" .and. ubndry == "periodic"
#if NDIMS == 3
lbndry = "periodic"
ubndry = "periodic"
call get_parameter("ylbndry", lbndry)
call get_parameter("yubndry", ubndry)
periodic(3) = lbndry == "periodic" .and. ubndry == "periodic"
#endif /* NDIMS == 3 */
! obtain the refinement level bounds
!
call get_parameter("minlev", minlev)
call get_parameter("maxlev", maxlev)
! check if the number refinement levels are correct
!
if (minlev < 1 .or. maxlev < 1 .or. minlev > maxlev) then
if (verbose) then
call print_message(loc, "Wrong minimum or maximum refinement levels:")
if (minlev < 1) then
write(*,"(5x,a,i0,a)") "=> minlev = ", minlev, " < 1."
call print_message(loc, msg)
end if
if (maxlev < 1) then
write(*,"(5x,a,i0,a)") "=> maxlev = ", maxlev, " < 1."
call print_message(loc, msg)
end if
if (minlev > maxlev) then
write(*,"(5x,a,2(i0,a))") "=> minlev = ", minlev, &
" > maxlev = ", maxlev, "."
call print_message(loc, msg)
end if
end if
status = 1
return
end if
! set the top level
!
toplev = max(tlevel, maxlev)
! get the domain base dimensions
!
domain_base_dims(1:NDIMS) = bdims(1:NDIMS)
! set the domain bounds
!
xmin = xmn
xmax = xmx
ymin = ymn
ymax = ymx
#if NDIMS == 3
zmin = zmn
zmax = zmx
#endif /* NDIMS == 3 */
! calculate the domain sizes
!
xlen = xmax - xmin
ylen = ymax - ymin
#if NDIMS == 3
zlen = zmax - zmin
#endif /* NDIMS == 3 */
! calculate the domain volume
!
vol = xlen * ylen * zlen
voli = 1.0d+00 / vol
! calculate the boundary areas
!
xarea = ylen * zlen
yarea = xlen * zlen
zarea = xlen * ylen
! allocate space for coordinate variables
!
#if NDIMS == 3
allocate(ax(toplev, bcells), ay(toplev, bcells), az(toplev, bcells), &
#else /* NDIMS == 3 */
allocate(ax(toplev, bcells), ay(toplev, bcells), az(toplev, 1 ), &
#endif /* NDIMS == 3 */
adx(toplev), ady(toplev), adz(toplev), adr(toplev), &
adxi(toplev), adyi(toplev), adzi(toplev), advol(toplev), &
stat = status)
if (status > 0) then
call print_message(loc, "Cannot allocate memory for coordinates!")
return
end if
! reset all coordinate variables to initial values
!
ax (:,:) = 0.0d+00
ay (:,:) = 0.0d+00
az (:,:) = 0.0d+00
adx (:) = 1.0d+00
ady (:) = 1.0d+00
adz (:) = 1.0d+00
adr (:) = 1.0d+00
adxi (:) = 1.0d+00
adyi (:) = 1.0d+00
adzi (:) = 1.0d+00
advol(:) = 1.0d+00
! generate the coordinate variables for each level
!
do l = 1, toplev
! calculate the block resolution at each level
!
nn = ni * 2**(l - 1)
! calculate the cell sizes for each level
!
adx (l) = xlen / (domain_base_dims(1) * nn)
ady (l) = ylen / (domain_base_dims(2) * nn)
#if NDIMS == 3
adz (l) = zlen / (domain_base_dims(3) * nn)
#endif /* NDIMS == 3 */
#if NDIMS == 2
adr (l) = sqrt(adx(l)**2 + ady(l)**2)
#endif /* NDIMS == 2 */
#if NDIMS == 3
adr (l) = sqrt(adx(l)**2 + ady(l)**2 + adz(l)**2)
#endif /* NDIMS == 3 */
! calculate the inverse of cell size
!
adxi(l) = 1.0d+00 / adx(l)
adyi(l) = 1.0d+00 / ady(l)
#if NDIMS == 3
adzi(l) = 1.0d+00 / adz(l)
#endif /* NDIMS == 3 */
! calculate the block coordinates for each level
!
ax(l,:) = ((/(i, i = 1, bcells)/) - ng - 5.0d-01) * adx(l)
ay(l,:) = ((/(j, j = 1, bcells)/) - ng - 5.0d-01) * ady(l)
#if NDIMS == 3
az(l,:) = ((/(k, k = 1, bcells)/) - ng - 5.0d-01) * adz(l)
#endif /* NDIMS == 3 */
! calculate the cell volume at each level
!
advol(l) = adx(l) * ady(l) * adz(l)
end do ! l = 1, toplev
! initialize ghost subarray indices
!
np = ni + ng
nm = ni - ng
nr = ni - nd
nq = ni - nh
ns = ncells_half
nt = np / 2
nu = nm / 2
#if NDIMS == 2
do j = 1, 2
fj = j - 1
q = 3 - j
do i = 1, 2
fi = i - 1
p = 3 - i
! for edge copy
!
edges_gc(i,j,1)%l(1) = nb + fi * ns
edges_gc(i,j,1)%l(2) = 1 + fj * np
edges_gc(i,j,2)%l(1) = 1 + fi * np
edges_gc(i,j,2)%l(2) = nb + fj * ns
edges_dc(i,q,1)%l(1) = nb + fi * ns
edges_dc(i,q,1)%l(2) = nb + fj * nm
edges_dc(p,j,2)%l(1) = nb + fi * nm
edges_dc(p,j,2)%l(2) = nb + fj * ns
edges_gc(i,j,1)%u(:) = edges_gc(i,j,1)%l(:) + (/ ns, ng /) - 1
edges_gc(i,j,2)%u(:) = edges_gc(i,j,2)%l(:) + (/ ng, ns /) - 1
edges_dc(i,q,1)%u(:) = edges_dc(i,q,1)%l(:) + (/ ns, ng /) - 1
edges_dc(p,j,2)%u(:) = edges_dc(p,j,2)%l(:) + (/ ng, ns /) - 1
! for edge restriction
!
edges_gr(i,j,1)%l(1) = nb + fi * ns
edges_gr(i,j,1)%l(2) = 1 + fj * np
edges_gr(i,j,2)%l(1) = 1 + fi * np
edges_gr(i,j,2)%l(2) = nb + fj * ns
edges_dr(i,q,1)%l(1) = nb
edges_dr(i,q,1)%l(2) = nb + fj * nr
edges_dr(p,j,2)%l(1) = nb + fi * nr
edges_dr(p,j,2)%l(2) = nb
edges_gr(i,j,1)%u(:) = edges_gr(i,j,1)%l(:) + (/ ns, ng /) - 1
edges_gr(i,j,2)%u(:) = edges_gr(i,j,2)%l(:) + (/ ng, ns /) - 1
edges_dr(i,q,1)%u(:) = edges_dr(i,q,1)%l(:) + (/ ni, nd /) - 1
edges_dr(p,j,2)%u(:) = edges_dr(p,j,2)%l(:) + (/ nd, ni /) - 1
! for edge prolongation
!
edges_gp(i,j,1)%l(1) = nb - fi * ng
edges_gp(i,j,1)%l(2) = 1 + fj * np
edges_gp(i,j,2)%l(1) = 1 + fi * np
edges_gp(i,j,2)%l(2) = nb - fj * ng
edges_dp(i,q,1)%l(1) = nb + fi * nu
edges_dp(i,q,1)%l(2) = nb + fj * nq
edges_dp(p,j,2)%l(1) = nb + fi * nq
edges_dp(p,j,2)%l(2) = nb + fj * nu
edges_gp(i,j,1)%u(:) = edges_gp(i,j,1)%l(:) + (/ np, ng /) - 1
edges_gp(i,j,2)%u(:) = edges_gp(i,j,2)%l(:) + (/ ng, np /) - 1
edges_dp(i,q,1)%u(:) = edges_dp(i,q,1)%l(:) + (/ nt, nh /) - 1
edges_dp(p,j,2)%u(:) = edges_dp(p,j,2)%l(:) + (/ nh, nt /) - 1
! for corner copy
!
corners_gc(i,j)%l(1) = 1 + fi * np
corners_gc(i,j)%l(2) = 1 + fj * np
corners_dc(p,q)%l(1) = nb + fi * nm
corners_dc(p,q)%l(2) = nb + fj * nm
corners_gc(i,j)%u(:) = corners_gc(i,j)%l(:) + ng - 1
corners_dc(p,q)%u(:) = corners_dc(p,q)%l(:) + ng - 1
! for corner restriction
!
corners_gr(i,j)%l(1) = 1 + fi * np
corners_gr(i,j)%l(2) = 1 + fj * np
corners_dr(p,q)%l(1) = nb + fi * nr
corners_dr(p,q)%l(2) = nb + fj * nr
corners_gr(i,j)%u(:) = corners_gr(i,j)%l(:) + ng - 1
corners_dr(p,q)%u(:) = corners_dr(p,q)%l(:) + nd - 1
! for corner prolongation
!
corners_gp(i,j)%l(1) = 1 + fi * np
corners_gp(i,j)%l(2) = 1 + fj * np
corners_dp(p,q)%l(1) = nb + fi * nq
corners_dp(p,q)%l(2) = nb + fj * nq
corners_gp(i,j)%u(:) = corners_gp(i,j)%l(:) + ng - 1
corners_dp(p,q)%u(:) = corners_dp(p,q)%l(:) + nh - 1
end do ! i = 1, 2
end do ! j = 1, 2
#endif /* NDIMS == 2 */
#if NDIMS == 3
do k = 1, 2
fk = k - 1
r = 3 - k
do j = 1, 2
fj = j - 1
q = 3 - j
do i = 1, 2
fi = i - 1
p = 3 - i
! for face copy
!
faces_gc(i,j,k,1)%l(1) = 1 + fi * np
faces_gc(i,j,k,1)%l(2) = nb + fj * ns
faces_gc(i,j,k,1)%l(3) = nb + fk * ns
faces_gc(i,j,k,2)%l(1) = nb + fi * ns
faces_gc(i,j,k,2)%l(2) = 1 + fj * np
faces_gc(i,j,k,2)%l(3) = nb + fk * ns
faces_gc(i,j,k,3)%l(1) = nb + fi * ns
faces_gc(i,j,k,3)%l(2) = nb + fj * ns
faces_gc(i,j,k,3)%l(3) = 1 + fk * np
faces_dc(p,j,k,1)%l(1) = nb + fi * nm
faces_dc(p,j,k,1)%l(2) = nb + fj * ns
faces_dc(p,j,k,1)%l(3) = nb + fk * ns
faces_dc(i,q,k,2)%l(1) = nb + fi * ns
faces_dc(i,q,k,2)%l(2) = nb + fj * nm
faces_dc(i,q,k,2)%l(3) = nb + fk * ns
faces_dc(i,j,r,3)%l(1) = nb + fi * ns
faces_dc(i,j,r,3)%l(2) = nb + fj * ns
faces_dc(i,j,r,3)%l(3) = nb + fk * nm
faces_gc(i,j,k,1)%u(:) = faces_gc(i,j,k,1)%l(:) + (/ ng, ns, ns /) - 1
faces_gc(i,j,k,2)%u(:) = faces_gc(i,j,k,2)%l(:) + (/ ns, ng, ns /) - 1
faces_gc(i,j,k,3)%u(:) = faces_gc(i,j,k,3)%l(:) + (/ ns, ns, ng /) - 1
faces_dc(p,j,k,1)%u(:) = faces_dc(p,j,k,1)%l(:) + (/ ng, ns, ns /) - 1
faces_dc(i,q,k,2)%u(:) = faces_dc(i,q,k,2)%l(:) + (/ ns, ng, ns /) - 1
faces_dc(i,j,r,3)%u(:) = faces_dc(i,j,r,3)%l(:) + (/ ns, ns, ng /) - 1
! for face restriction
!
faces_gr(i,j,k,1)%l(1) = 1 + fi * np
faces_gr(i,j,k,1)%l(2) = nb + fj * ns
faces_gr(i,j,k,1)%l(3) = nb + fk * ns
faces_gr(i,j,k,2)%l(1) = nb + fi * ns
faces_gr(i,j,k,2)%l(2) = 1 + fj * np
faces_gr(i,j,k,2)%l(3) = nb + fk * ns
faces_gr(i,j,k,3)%l(1) = nb + fi * ns
faces_gr(i,j,k,3)%l(2) = nb + fj * ns
faces_gr(i,j,k,3)%l(3) = 1 + fk * np
faces_dr(p,j,k,1)%l(1) = nb + fi * nr
faces_dr(p,j,k,1)%l(2) = nb
faces_dr(p,j,k,1)%l(3) = nb
faces_dr(i,q,k,2)%l(1) = nb
faces_dr(i,q,k,2)%l(2) = nb + fj * nr
faces_dr(i,q,k,2)%l(3) = nb
faces_dr(i,j,r,3)%l(1) = nb
faces_dr(i,j,r,3)%l(2) = nb
faces_dr(i,j,r,3)%l(3) = nb + fk * nr
faces_gr(i,j,k,1)%u(:) = faces_gr(i,j,k,1)%l(:) + (/ ng, ns, ns /) - 1
faces_gr(i,j,k,2)%u(:) = faces_gr(i,j,k,2)%l(:) + (/ ns, ng, ns /) - 1
faces_gr(i,j,k,3)%u(:) = faces_gr(i,j,k,3)%l(:) + (/ ns, ns, ng /) - 1
faces_dr(p,j,k,1)%u(:) = faces_dr(p,j,k,1)%l(:) + (/ nd, ni, ni /) - 1
faces_dr(i,q,k,2)%u(:) = faces_dr(i,q,k,2)%l(:) + (/ ni, nd, ni /) - 1
faces_dr(i,j,r,3)%u(:) = faces_dr(i,j,r,3)%l(:) + (/ ni, ni, nd /) - 1
! for face prolongation
!
faces_gp(i,j,k,1)%l(1) = 1 + fi * np
faces_gp(i,j,k,1)%l(2) = nb - fj * ng
faces_gp(i,j,k,1)%l(3) = nb - fk * ng
faces_gp(i,j,k,2)%l(1) = nb - fi * ng
faces_gp(i,j,k,2)%l(2) = 1 + fj * np
faces_gp(i,j,k,2)%l(3) = nb - fk * ng
faces_gp(i,j,k,3)%l(1) = nb - fi * ng
faces_gp(i,j,k,3)%l(2) = nb - fj * ng
faces_gp(i,j,k,3)%l(3) = 1 + fk * np
faces_dp(p,j,k,1)%l(1) = nb + fi * nq
faces_dp(p,j,k,1)%l(2) = nb + fj * nu
faces_dp(p,j,k,1)%l(3) = nb + fk * nu
faces_dp(i,q,k,2)%l(1) = nb + fi * nu
faces_dp(i,q,k,2)%l(2) = nb + fj * nq
faces_dp(i,q,k,2)%l(3) = nb + fk * nu
faces_dp(i,j,r,3)%l(1) = nb + fi * nu
faces_dp(i,j,r,3)%l(2) = nb + fj * nu
faces_dp(i,j,r,3)%l(3) = nb + fk * nq
faces_gp(i,j,k,1)%u(:) = faces_gp(i,j,k,1)%l(:) + (/ ng, np, np /) - 1
faces_gp(i,j,k,2)%u(:) = faces_gp(i,j,k,2)%l(:) + (/ np, ng, np /) - 1
faces_gp(i,j,k,3)%u(:) = faces_gp(i,j,k,3)%l(:) + (/ np, np, ng /) - 1
faces_dp(p,j,k,1)%u(:) = faces_dp(p,j,k,1)%l(:) + (/ nh, nt, nt /) - 1
faces_dp(i,q,k,2)%u(:) = faces_dp(i,q,k,2)%l(:) + (/ nt, nh, nt /) - 1
faces_dp(i,j,r,3)%u(:) = faces_dp(i,j,r,3)%l(:) + (/ nt, nt, nh /) - 1
! for edge copy
!
edges_gc(i,j,k,1)%l(1) = nb + fi * ns
edges_gc(i,j,k,1)%l(2) = 1 + fj * np
edges_gc(i,j,k,1)%l(3) = 1 + fk * np
edges_gc(i,j,k,2)%l(1) = 1 + fi * np
edges_gc(i,j,k,2)%l(2) = nb + fj * ns
edges_gc(i,j,k,2)%l(3) = 1 + fk * np
edges_gc(i,j,k,3)%l(1) = 1 + fi * np
edges_gc(i,j,k,3)%l(2) = 1 + fj * np
edges_gc(i,j,k,3)%l(3) = nb + fk * ns
edges_dc(i,q,r,1)%l(1) = nb + fi * ns
edges_dc(i,q,r,1)%l(2) = nb + fj * nm
edges_dc(i,q,r,1)%l(3) = nb + fk * nm
edges_dc(p,j,r,2)%l(1) = nb + fi * nm
edges_dc(p,j,r,2)%l(2) = nb + fj * ns
edges_dc(p,j,r,2)%l(3) = nb + fk * nm
edges_dc(p,q,k,3)%l(1) = nb + fi * nm
edges_dc(p,q,k,3)%l(2) = nb + fj * nm
edges_dc(p,q,k,3)%l(3) = nb + fk * ns
edges_gc(i,j,k,1)%u(:) = edges_gc(i,j,k,1)%l(:) + (/ ns, ng, ng /) - 1
edges_gc(i,j,k,2)%u(:) = edges_gc(i,j,k,2)%l(:) + (/ ng, ns, ng /) - 1
edges_gc(i,j,k,3)%u(:) = edges_gc(i,j,k,3)%l(:) + (/ ng, ng, ns /) - 1
edges_dc(i,q,r,1)%u(:) = edges_dc(i,q,r,1)%l(:) + (/ ns, ng, ng /) - 1
edges_dc(p,j,r,2)%u(:) = edges_dc(p,j,r,2)%l(:) + (/ ng, ns, ng /) - 1
edges_dc(p,q,k,3)%u(:) = edges_dc(p,q,k,3)%l(:) + (/ ng, ng, ns /) - 1
! for edge restriction
!
edges_gr(i,j,k,1)%l(1) = nb + fi * ns
edges_gr(i,j,k,1)%l(2) = 1 + fj * np
edges_gr(i,j,k,1)%l(3) = 1 + fk * np
edges_gr(i,j,k,2)%l(1) = 1 + fi * np
edges_gr(i,j,k,2)%l(2) = nb + fj * ns
edges_gr(i,j,k,2)%l(3) = 1 + fk * np
edges_gr(i,j,k,3)%l(1) = 1 + fi * np
edges_gr(i,j,k,3)%l(2) = 1 + fj * np
edges_gr(i,j,k,3)%l(3) = nb + fk * ns
edges_dr(i,q,r,1)%l(1) = nb
edges_dr(i,q,r,1)%l(2) = nb + fj * nr
edges_dr(i,q,r,1)%l(3) = nb + fk * nr
edges_dr(p,j,r,2)%l(1) = nb + fi * nr
edges_dr(p,j,r,2)%l(2) = nb
edges_dr(p,j,r,2)%l(3) = nb + fk * nr
edges_dr(p,q,k,3)%l(1) = nb + fi * nr
edges_dr(p,q,k,3)%l(2) = nb + fj * nr
edges_dr(p,q,k,3)%l(3) = nb
edges_gr(i,j,k,1)%u(:) = edges_gr(i,j,k,1)%l(:) + (/ ns, ng, ng /) - 1
edges_gr(i,j,k,2)%u(:) = edges_gr(i,j,k,2)%l(:) + (/ ng, ns, ng /) - 1
edges_gr(i,j,k,3)%u(:) = edges_gr(i,j,k,3)%l(:) + (/ ng, ng, ns /) - 1
edges_dr(i,q,r,1)%u(:) = edges_dr(i,q,r,1)%l(:) + (/ ni, nd, nd /) - 1
edges_dr(p,j,r,2)%u(:) = edges_dr(p,j,r,2)%l(:) + (/ nd, ni, nd /) - 1
edges_dr(p,q,k,3)%u(:) = edges_dr(p,q,k,3)%l(:) + (/ nd, nd, ni /) - 1
! for edge prolongation
!
edges_gp(i,j,k,1)%l(1) = nb - fi * ng
edges_gp(i,j,k,1)%l(2) = 1 + fj * np
edges_gp(i,j,k,1)%l(3) = 1 + fk * np
edges_gp(i,j,k,2)%l(1) = 1 + fi * np
edges_gp(i,j,k,2)%l(2) = nb - fj * ng
edges_gp(i,j,k,2)%l(3) = 1 + fk * np
edges_gp(i,j,k,3)%l(1) = 1 + fi * np
edges_gp(i,j,k,3)%l(2) = 1 + fj * np
edges_gp(i,j,k,3)%l(3) = nb - fk * ng
edges_dp(i,q,r,1)%l(1) = nb + fi * nu
edges_dp(i,q,r,1)%l(2) = nb + fj * nq
edges_dp(i,q,r,1)%l(3) = nb + fk * nq
edges_dp(p,j,r,2)%l(1) = nb + fi * nq
edges_dp(p,j,r,2)%l(2) = nb + fj * nu
edges_dp(p,j,r,2)%l(3) = nb + fk * nq
edges_dp(p,q,k,3)%l(1) = nb + fi * nq
edges_dp(p,q,k,3)%l(2) = nb + fj * nq
edges_dp(p,q,k,3)%l(3) = nb + fk * nu
edges_gp(i,j,k,1)%u(:) = edges_gp(i,j,k,1)%l(:) + (/ np, ng, ng /) - 1
edges_gp(i,j,k,2)%u(:) = edges_gp(i,j,k,2)%l(:) + (/ ng, np, ng /) - 1
edges_gp(i,j,k,3)%u(:) = edges_gp(i,j,k,3)%l(:) + (/ ng, ng, np /) - 1
edges_dp(i,q,r,1)%u(:) = edges_dp(i,q,r,1)%l(:) + (/ nt, nh, nh /) - 1
edges_dp(p,j,r,2)%u(:) = edges_dp(p,j,r,2)%l(:) + (/ nh, nt, nh /) - 1
edges_dp(p,q,k,3)%u(:) = edges_dp(p,q,k,3)%l(:) + (/ nh, nh, nt /) - 1
! for corner copy
!
corners_gc(i,j,k)%l(1) = 1 + fi * np
corners_gc(i,j,k)%l(2) = 1 + fj * np
corners_gc(i,j,k)%l(3) = 1 + fk * np
corners_dc(p,q,r)%l(1) = nb + fi * nm
corners_dc(p,q,r)%l(2) = nb + fj * nm
corners_dc(p,q,r)%l(3) = nb + fk * nm
corners_gc(i,j,k)%u(:) = corners_gc(i,j,k)%l(:) + ng - 1
corners_dc(p,q,r)%u(:) = corners_dc(p,q,r)%l(:) + ng - 1
! for corner restriction
!
corners_gr(i,j,k)%l(1) = 1 + fi * np
corners_gr(i,j,k)%l(2) = 1 + fj * np
corners_gr(i,j,k)%l(3) = 1 + fk * np
corners_dr(p,q,r)%l(1) = nb + fi * nr
corners_dr(p,q,r)%l(2) = nb + fj * nr
corners_dr(p,q,r)%l(3) = nb + fk * nr
corners_gr(i,j,k)%u(:) = corners_gr(i,j,k)%l(:) + ng - 1
corners_dr(p,q,r)%u(:) = corners_dr(p,q,r)%l(:) + nd - 1
! for corner prolongation
!
corners_gp(i,j,k)%l(1) = 1 + fi * np
corners_gp(i,j,k)%l(2) = 1 + fj * np
corners_gp(i,j,k)%l(3) = 1 + fk * np
corners_dp(p,q,r)%l(1) = nb + fi * nq
corners_dp(p,q,r)%l(2) = nb + fj * nq
corners_dp(p,q,r)%l(3) = nb + fk * nq
corners_gp(i,j,k)%u(:) = corners_gp(i,j,k)%l(:) + ng - 1
corners_dp(p,q,r)%u(:) = corners_dp(p,q,r)%l(:) + nh - 1
end do ! i = 1, 2
end do ! j = 1, 2
end do ! k = 1, 2
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine initialize_coordinates
!
!===============================================================================
!
! subroutine FINALIZE_COORDINATES:
! -------------------------------
!
! Subroutine deallocates mesh coordinates.
!
! Arguments:
!
! status - return flag of the procedure execution status;
!
!===============================================================================
!
subroutine finalize_coordinates(status)
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
!
!-------------------------------------------------------------------------------
!
! reset the status flag
!
status = 0
! deallocating coordinate variables
!
if (allocated(ax) ) deallocate(ax , stat = status)
if (allocated(ay) ) deallocate(ay , stat = status)
if (allocated(az) ) deallocate(az , stat = status)
if (allocated(adx) ) deallocate(adx , stat = status)
if (allocated(ady) ) deallocate(ady , stat = status)
if (allocated(adz) ) deallocate(adz , stat = status)
if (allocated(adr) ) deallocate(adr , stat = status)
if (allocated(adxi) ) deallocate(adxi , stat = status)
if (allocated(adyi) ) deallocate(adyi , stat = status)
if (allocated(adzi) ) deallocate(adzi , stat = status)
if (allocated(advol)) deallocate(advol, stat = status)
!-------------------------------------------------------------------------------
!
end subroutine finalize_coordinates
!
!===============================================================================
!
! subroutine PRINT_COORDINATES:
! ----------------------------
!
! Subroutine print module parameters.
!
! Arguments:
!
! verbose - flag determining if the subroutine should be verbose;
!
!===============================================================================
!
subroutine print_coordinates(verbose)
use helpers , only : print_section, print_parameter
implicit none
logical, intent(in) :: verbose
character(len= 80) :: sfmt
character(len=255) :: msg
integer :: p, q
integer(kind=4), dimension(3) :: bm, rm, cm, fm
!-------------------------------------------------------------------------------
!
if (verbose) then
bm(:) = domain_base_dims(:)
rm(:) = bm(:) * 2**(maxlev - 1)
cm(:) = bm(:) * ncells
fm(:) = rm(:) * ncells
call print_section(verbose, "Geometry")
call print_parameter(verbose, "refinement to level" , maxlev )
call print_parameter(verbose, "number of block cells", ncells )
call print_parameter(verbose, "number of ghost zones", nghosts)
write(msg,"(i0)") maxval(fm(1:NDIMS))
p = len(trim(adjustl(msg)))
write(msg,"(i0)") maxval(rm(1:NDIMS))
q = len(trim(adjustl(msg)))
#if NDIMS == 3
write(sfmt,"(6(a,i0),a)") "('[',i", p, ",' x ',i", p, ",' x ',i", p, &
",'] ([',i", q, ",' x ',i", q, ",' x ',i", q, ",'] blocks)')"
#else /* NDIMS == 3 */
write(sfmt,"(4(a,i0),a)") "('[',i", p, ",' x ',i", p, &
",'] ([',i", q, ",' x ',i", q, ",'] blocks)')"
#endif /* NDIMS == 3 */
write(msg,sfmt) cm(1:NDIMS), bm(1:NDIMS)
call print_parameter(verbose, "base resolution" , msg)
write(msg,sfmt) fm(1:NDIMS), rm(1:NDIMS)
call print_parameter(verbose, "effective resolution" , msg)
call print_parameter(verbose, "X-bounds", xmin, xmax)
call print_parameter(verbose, "Y-bounds", ymin, ymax)
#if NDIMS == 3
call print_parameter(verbose, "Z-bounds", zmin, zmax)
#endif /* NDIMS */
end if
!-------------------------------------------------------------------------------
!
end subroutine print_coordinates
!===============================================================================
!
end module