810 lines
25 KiB
Fortran
810 lines
25 KiB
Fortran
!!******************************************************************************
|
|
!!
|
|
!! 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-2019 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 domain block dimensions
|
|
!
|
|
integer , save :: nc = 8, in = 8, jn = 8, kn = 1
|
|
|
|
! the number of ghost zones
|
|
!
|
|
integer , save :: ng = 2, nh = 1, nd = 4
|
|
|
|
! the domain block dimensions including the ghost zones
|
|
!
|
|
integer , save :: im = 12, jm = 12, km = 1
|
|
|
|
! the domain block dimensions including the ghost zones
|
|
!
|
|
integer , save :: it = 11, jt = 11, kt = 1
|
|
|
|
! the domain division
|
|
!
|
|
integer , save :: ir = 1, jr = 1, kr = 1
|
|
|
|
! the limits of refinement level
|
|
!
|
|
integer , save :: minlev = 1, maxlev = 1, toplev = 1
|
|
|
|
! block indices
|
|
!
|
|
integer , save :: ih = 6, jh = 6, kh = 1
|
|
integer , save :: ib = 3, jb = 3, kb = 1
|
|
integer , save :: ie = 10, je = 10, ke = 1
|
|
integer , save :: ibl = 2, jbl = 2, kbl = 1
|
|
integer , save :: ibu = 4, jbu = 4, kbu = 1
|
|
integer , save :: iel = 9, jel = 9, kel = 1
|
|
integer , save :: ieu = 11, jeu = 11, keu = 1
|
|
|
|
! 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:
|
|
!
|
|
! verbose - flag determining if the subroutine should be verbose;
|
|
! iret - return flag of the procedure execution status;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine initialize_coordinates(tlev, verbose, iret)
|
|
|
|
! include external procedures
|
|
!
|
|
use parameters, only : get_parameter
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! subroutine arguments
|
|
!
|
|
logical, intent(in) :: verbose
|
|
integer, intent(in) :: tlev
|
|
integer, intent(inout) :: iret
|
|
|
|
! local variables
|
|
!
|
|
integer :: i, j, k, l, p, q, r, ff
|
|
integer :: fi, fj, fk
|
|
integer :: ni, nj, nk, nm, np, nr, nq, ns, nt, nu
|
|
logical :: info
|
|
|
|
! local arrays
|
|
!
|
|
integer(kind=4), dimension(3) :: cm, rm, dm
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! obtain the number of cells along each block dimension
|
|
!
|
|
call get_parameter("ncells" , nc)
|
|
|
|
! obtain the number of ghost zones
|
|
!
|
|
call get_parameter("nghosts", ng)
|
|
|
|
! set the block dimensions
|
|
!
|
|
in = nc
|
|
jn = nc
|
|
#if NDIMS == 3
|
|
kn = nc
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! calculate half and double of the number of ghose zones
|
|
!
|
|
nh = ng / 2
|
|
nd = ng * 2
|
|
|
|
! calculate the block dimensions including ghost cells
|
|
!
|
|
im = in + 2 * ng
|
|
jm = jn + 2 * ng
|
|
#if NDIMS == 3
|
|
km = kn + 2 * ng
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! prepare indices
|
|
!
|
|
it = im - nh + 1
|
|
jt = jm - nh + 1
|
|
#if NDIMS == 3
|
|
kt = km - nh + 1
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! calculate block indices
|
|
!
|
|
ih = im / 2
|
|
ib = ng + 1
|
|
ie = ng + in
|
|
ibl = ib - 1
|
|
ibu = ib + ng - 1
|
|
iel = ie - ng + 1
|
|
ieu = ie + 1
|
|
jh = jm / 2
|
|
jb = ng + 1
|
|
je = ng + jn
|
|
jbl = jb - 1
|
|
jbu = jb + ng - 1
|
|
jel = je - ng + 1
|
|
jeu = je + 1
|
|
#if NDIMS == 3
|
|
kh = km / 2
|
|
kb = ng + 1
|
|
ke = ng + kn
|
|
kbl = kb - 1
|
|
kbu = kb + ng - 1
|
|
kel = ke - ng + 1
|
|
keu = ke + 1
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! obtain the refinement level bounds
|
|
!
|
|
call get_parameter("minlev", minlev)
|
|
call get_parameter("maxlev", maxlev)
|
|
|
|
! set the top level
|
|
!
|
|
toplev = max(tlev, maxlev)
|
|
|
|
! obtain the domain base division
|
|
!
|
|
call get_parameter("xblocks", ir )
|
|
call get_parameter("yblocks", jr )
|
|
#if NDIMS == 3
|
|
call get_parameter("zblocks", kr )
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! obtain the domain bounds
|
|
!
|
|
call get_parameter("xmin" , xmin )
|
|
call get_parameter("xmax" , xmax )
|
|
call get_parameter("ymin" , ymin )
|
|
call get_parameter("ymax" , ymax )
|
|
#if NDIMS == 3
|
|
call get_parameter("zmin" , zmin )
|
|
call get_parameter("zmax" , zmax )
|
|
#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
|
|
!
|
|
allocate(ax (toplev, im))
|
|
allocate(ay (toplev, jm))
|
|
allocate(az (toplev, km))
|
|
allocate(adx (toplev))
|
|
allocate(ady (toplev))
|
|
allocate(adz (toplev))
|
|
allocate(adr (toplev))
|
|
allocate(adxi (toplev))
|
|
allocate(adyi (toplev))
|
|
allocate(adzi (toplev))
|
|
allocate(advol(toplev))
|
|
|
|
! 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
|
|
!
|
|
ff = 2**(l - 1)
|
|
ni = in * ff
|
|
nj = jn * ff
|
|
nk = kn * ff
|
|
|
|
! calculate the cell sizes for each level
|
|
!
|
|
adx (l) = xlen / (ir * ni)
|
|
ady (l) = ylen / (jr * nj)
|
|
#if NDIMS == 3
|
|
adz (l) = zlen / (kr * nk)
|
|
#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, im)/) - ng - 5.0d-01) * adx(l)
|
|
ay(l,:) = ((/(j, j = 1, jm)/) - ng - 5.0d-01) * ady(l)
|
|
#if NDIMS == 3
|
|
az(l,:) = ((/(k, k = 1, km)/) - 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 = nc + ng
|
|
nm = nc - ng
|
|
nr = nc - nd
|
|
nq = nc - nh
|
|
ns = nc / 2
|
|
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) = ib + 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) = jb + fj * ns
|
|
|
|
edges_dc(i,q,1)%l(1) = ib + fi * ns
|
|
edges_dc(i,q,1)%l(2) = jb + fj * nm
|
|
edges_dc(p,j,2)%l(1) = ib + fi * nm
|
|
edges_dc(p,j,2)%l(2) = jb + 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) = ib + 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) = jb + fj * ns
|
|
|
|
edges_dr(i,q,1)%l(1) = ib
|
|
edges_dr(i,q,1)%l(2) = jb + fj * nr
|
|
edges_dr(p,j,2)%l(1) = ib + fi * nr
|
|
edges_dr(p,j,2)%l(2) = jb
|
|
|
|
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(:) + (/ nc, nd /) - 1
|
|
edges_dr(p,j,2)%u(:) = edges_dr(p,j,2)%l(:) + (/ nd, nc /) - 1
|
|
|
|
! for edge prolongation
|
|
!
|
|
edges_gp(i,j,1)%l(1) = ib - 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) = jb - fj * ng
|
|
|
|
edges_dp(i,q,1)%l(1) = ib + fi * nu
|
|
edges_dp(i,q,1)%l(2) = jb + fj * nq
|
|
edges_dp(p,j,2)%l(1) = ib + fi * nq
|
|
edges_dp(p,j,2)%l(2) = jb + 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) = ib + fi * nm
|
|
corners_dc(p,q)%l(2) = jb + 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) = ib + fi * nr
|
|
corners_dr(p,q)%l(2) = jb + 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) = ib + fi * nq
|
|
corners_dp(p,q)%l(2) = jb + 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) = jb + fj * ns
|
|
faces_gc(i,j,k,1)%l(3) = kb + fk * ns
|
|
faces_gc(i,j,k,2)%l(1) = ib + fi * ns
|
|
faces_gc(i,j,k,2)%l(2) = 1 + fj * np
|
|
faces_gc(i,j,k,2)%l(3) = kb + fk * ns
|
|
faces_gc(i,j,k,3)%l(1) = ib + fi * ns
|
|
faces_gc(i,j,k,3)%l(2) = jb + fj * ns
|
|
faces_gc(i,j,k,3)%l(3) = 1 + fk * np
|
|
|
|
faces_dc(p,j,k,1)%l(1) = ib + fi * nm
|
|
faces_dc(p,j,k,1)%l(2) = jb + fj * ns
|
|
faces_dc(p,j,k,1)%l(3) = kb + fk * ns
|
|
faces_dc(i,q,k,2)%l(1) = ib + fi * ns
|
|
faces_dc(i,q,k,2)%l(2) = jb + fj * nm
|
|
faces_dc(i,q,k,2)%l(3) = kb + fk * ns
|
|
faces_dc(i,j,r,3)%l(1) = ib + fi * ns
|
|
faces_dc(i,j,r,3)%l(2) = jb + fj * ns
|
|
faces_dc(i,j,r,3)%l(3) = kb + 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) = jb + fj * ns
|
|
faces_gr(i,j,k,1)%l(3) = kb + fk * ns
|
|
faces_gr(i,j,k,2)%l(1) = ib + fi * ns
|
|
faces_gr(i,j,k,2)%l(2) = 1 + fj * np
|
|
faces_gr(i,j,k,2)%l(3) = kb + fk * ns
|
|
faces_gr(i,j,k,3)%l(1) = ib + fi * ns
|
|
faces_gr(i,j,k,3)%l(2) = jb + fj * ns
|
|
faces_gr(i,j,k,3)%l(3) = 1 + fk * np
|
|
|
|
faces_dr(p,j,k,1)%l(1) = ib + fi * nr
|
|
faces_dr(p,j,k,1)%l(2) = jb
|
|
faces_dr(p,j,k,1)%l(3) = kb
|
|
faces_dr(i,q,k,2)%l(1) = ib
|
|
faces_dr(i,q,k,2)%l(2) = jb + fj * nr
|
|
faces_dr(i,q,k,2)%l(3) = kb
|
|
faces_dr(i,j,r,3)%l(1) = ib
|
|
faces_dr(i,j,r,3)%l(2) = jb
|
|
faces_dr(i,j,r,3)%l(3) = kb + 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, nc, nc /) - 1
|
|
faces_dr(i,q,k,2)%u(:) = faces_dr(i,q,k,2)%l(:) + (/ nc, nd, nc /) - 1
|
|
faces_dr(i,j,r,3)%u(:) = faces_dr(i,j,r,3)%l(:) + (/ nc, nc, nd /) - 1
|
|
|
|
! for face prolongation
|
|
!
|
|
faces_gp(i,j,k,1)%l(1) = 1 + fi * np
|
|
faces_gp(i,j,k,1)%l(2) = jb - fj * ng
|
|
faces_gp(i,j,k,1)%l(3) = kb - fk * ng
|
|
faces_gp(i,j,k,2)%l(1) = ib - fi * ng
|
|
faces_gp(i,j,k,2)%l(2) = 1 + fj * np
|
|
faces_gp(i,j,k,2)%l(3) = kb - fk * ng
|
|
faces_gp(i,j,k,3)%l(1) = ib - fi * ng
|
|
faces_gp(i,j,k,3)%l(2) = jb - fj * ng
|
|
faces_gp(i,j,k,3)%l(3) = 1 + fk * np
|
|
|
|
faces_dp(p,j,k,1)%l(1) = ib + fi * nq
|
|
faces_dp(p,j,k,1)%l(2) = jb + fj * nu
|
|
faces_dp(p,j,k,1)%l(3) = kb + fk * nu
|
|
faces_dp(i,q,k,2)%l(1) = ib + fi * nu
|
|
faces_dp(i,q,k,2)%l(2) = jb + fj * nq
|
|
faces_dp(i,q,k,2)%l(3) = kb + fk * nu
|
|
faces_dp(i,j,r,3)%l(1) = ib + fi * nu
|
|
faces_dp(i,j,r,3)%l(2) = jb + fj * nu
|
|
faces_dp(i,j,r,3)%l(3) = kb + 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) = ib + 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) = jb + 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) = kb + fk * ns
|
|
|
|
edges_dc(i,q,r,1)%l(1) = ib + fi * ns
|
|
edges_dc(i,q,r,1)%l(2) = jb + fj * nm
|
|
edges_dc(i,q,r,1)%l(3) = kb + fk * nm
|
|
edges_dc(p,j,r,2)%l(1) = jb + fi * nm
|
|
edges_dc(p,j,r,2)%l(2) = jb + fj * ns
|
|
edges_dc(p,j,r,2)%l(3) = kb + fk * nm
|
|
edges_dc(p,q,k,3)%l(1) = ib + fi * nm
|
|
edges_dc(p,q,k,3)%l(2) = jb + fj * nm
|
|
edges_dc(p,q,k,3)%l(3) = kb + 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) = ib + 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) = jb + 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) = kb + fk * ns
|
|
|
|
edges_dr(i,q,r,1)%l(1) = ib
|
|
edges_dr(i,q,r,1)%l(2) = jb + fj * nr
|
|
edges_dr(i,q,r,1)%l(3) = kb + fk * nr
|
|
edges_dr(p,j,r,2)%l(1) = ib + fi * nr
|
|
edges_dr(p,j,r,2)%l(2) = jb
|
|
edges_dr(p,j,r,2)%l(3) = kb + fk * nr
|
|
edges_dr(p,q,k,3)%l(1) = ib + fi * nr
|
|
edges_dr(p,q,k,3)%l(2) = jb + fj * nr
|
|
edges_dr(p,q,k,3)%l(3) = kb
|
|
|
|
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(:) + (/ nc, nd, nd /) - 1
|
|
edges_dr(p,j,r,2)%u(:) = edges_dr(p,j,r,2)%l(:) + (/ nd, nc, nd /) - 1
|
|
edges_dr(p,q,k,3)%u(:) = edges_dr(p,q,k,3)%l(:) + (/ nd, nd, nc /) - 1
|
|
|
|
! for edge prolongation
|
|
!
|
|
edges_gp(i,j,k,1)%l(1) = ib - 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) = jb - 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) = kb - fk * ng
|
|
|
|
edges_dp(i,q,r,1)%l(1) = ib + fi * nu
|
|
edges_dp(i,q,r,1)%l(2) = jb + fj * nq
|
|
edges_dp(i,q,r,1)%l(3) = kb + fk * nq
|
|
edges_dp(p,j,r,2)%l(1) = ib + fi * nq
|
|
edges_dp(p,j,r,2)%l(2) = jb + fj * nu
|
|
edges_dp(p,j,r,2)%l(3) = kb + fk * nq
|
|
edges_dp(p,q,k,3)%l(1) = ib + fi * nq
|
|
edges_dp(p,q,k,3)%l(2) = jb + fj * nq
|
|
edges_dp(p,q,k,3)%l(3) = kb + 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) = ib + fi * nm
|
|
corners_dc(p,q,r)%l(2) = jb + fj * nm
|
|
corners_dc(p,q,r)%l(3) = kb + 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) = ib + fi * nr
|
|
corners_dr(p,q,r)%l(2) = jb + fj * nr
|
|
corners_dr(p,q,r)%l(3) = kb + 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) = ib + fi * nq
|
|
corners_dp(p,q,r)%l(2) = jb + fj * nq
|
|
corners_dp(p,q,r)%l(3) = kb + 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 */
|
|
|
|
! print general information about the level resolutions
|
|
!
|
|
if (verbose) then
|
|
|
|
! the base resolution
|
|
!
|
|
cm(1) = ir * in
|
|
cm(2) = jr * jn
|
|
cm(3) = kr * kn
|
|
|
|
! the effective resolution
|
|
!
|
|
ff = 2**(maxlev - 1)
|
|
rm(1) = cm(1) * ff
|
|
rm(2) = cm(2) * ff
|
|
rm(3) = cm(3) * ff
|
|
|
|
! the top level block division
|
|
!
|
|
dm(1) = rm(1) / in
|
|
dm(2) = rm(2) / jn
|
|
dm(3) = rm(3) / kn
|
|
|
|
! obtain the maximum number of block
|
|
!
|
|
ff = product(dm(1:NDIMS))
|
|
|
|
! print info
|
|
!
|
|
write(*,"(4x,a, 1x,i6 )" ) "refinement to level =", toplev
|
|
write(*,"(4x,a,3(1x,i6 ))") "base configuration =", ir, jr, kr
|
|
write(*,"(4x,a,3(1x,i6 ))") "top level blocks =", dm(1:NDIMS)
|
|
write(*,"(4x,a, 3x,i18)" ) "maximum cover blocks =", ff
|
|
write(*,"(4x,a,3(1x,i6 ))") "base resolution =", cm(1:NDIMS)
|
|
write(*,"(4x,a,3(1x,i6 ))") "effective resolution =", rm(1:NDIMS)
|
|
|
|
end if ! verbose
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine initialize_coordinates
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine FINALIZE_COORDINATES:
|
|
! -------------------------------
|
|
!
|
|
! Subroutine deallocates mesh coordinates.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! iret - return flag of the procedure execution status;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine finalize_coordinates(iret)
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! subroutine arguments
|
|
!
|
|
integer, intent(inout) :: iret
|
|
!
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! deallocating coordinate variables
|
|
!
|
|
if (allocated(ax) ) deallocate(ax)
|
|
if (allocated(ay) ) deallocate(ay)
|
|
if (allocated(az) ) deallocate(az)
|
|
if (allocated(adx) ) deallocate(adx)
|
|
if (allocated(ady) ) deallocate(ady)
|
|
if (allocated(adz) ) deallocate(adz)
|
|
if (allocated(adr) ) deallocate(adr)
|
|
if (allocated(adxi) ) deallocate(adxi)
|
|
if (allocated(adyi) ) deallocate(adyi)
|
|
if (allocated(adzi) ) deallocate(adzi)
|
|
if (allocated(advol)) deallocate(advol)
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine finalize_coordinates
|
|
|
|
!===============================================================================
|
|
!
|
|
end module
|