2195 lines
62 KiB
Fortran
2195 lines
62 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-2022 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: MESH - handling adaptive mesh structure
|
|
!!
|
|
!!******************************************************************************
|
|
!
|
|
module mesh
|
|
|
|
implicit none
|
|
|
|
! interface for the domain setup
|
|
!
|
|
abstract interface
|
|
subroutine setup_domain_iface(status)
|
|
integer, intent(out) :: status
|
|
end subroutine
|
|
end interface
|
|
|
|
! interfaces for the problem setup
|
|
!
|
|
abstract interface
|
|
subroutine setup_problem_iface(pdata)
|
|
use blocks, only : block_data
|
|
type(block_data), pointer, intent(inout) :: pdata
|
|
end subroutine
|
|
end interface
|
|
|
|
! pointer to the problem domain setup subroutine
|
|
!
|
|
procedure(setup_domain_iface) , pointer, save :: setup_domain => null()
|
|
|
|
! pointer to the problem initial setup subroutine
|
|
!
|
|
procedure(setup_problem_iface), pointer, save :: setup_problem => null()
|
|
|
|
! the flag indicating that the block structure or distribution has changed
|
|
!
|
|
logical, save :: changed = .true.
|
|
|
|
! the dimensions of the temporary array used in the block prolongation
|
|
!
|
|
integer, dimension(3), save :: pm = 1
|
|
|
|
private
|
|
|
|
public :: initialize_mesh, finalize_mesh, print_mesh
|
|
public :: generate_mesh, update_mesh
|
|
public :: redistribute_blocks
|
|
public :: setup_domain, setup_problem
|
|
public :: changed
|
|
|
|
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
!
|
|
contains
|
|
!
|
|
!===============================================================================
|
|
!!
|
|
!!*** PUBLIC SUBROUTINES *****************************************************
|
|
!!
|
|
!===============================================================================
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine INITIALIZE_MESH:
|
|
! --------------------------
|
|
!
|
|
! Subroutine initializes mesh module and its variables.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! verbose - the verbose flag;
|
|
! status - the subroutine call status;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine initialize_mesh(verbose, status)
|
|
|
|
use coordinates, only : ncells, nghosts
|
|
use refinement , only : initialize_refinement
|
|
|
|
implicit none
|
|
|
|
logical, intent(in) :: verbose
|
|
integer, intent(out) :: status
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
setup_domain => setup_domain_default
|
|
|
|
pm(1:NDIMS) = 2 * (ncells + nghosts)
|
|
|
|
call initialize_refinement(verbose, status)
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine initialize_mesh
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine FINALIZE_MESH:
|
|
! ------------------------
|
|
!
|
|
! Subroutine releases memory used by the module variables.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine finalize_mesh(status)
|
|
|
|
use refinement, only : finalize_refinement
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
nullify(setup_problem)
|
|
|
|
call finalize_refinement(status)
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine finalize_mesh
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine PRINT_MESH:
|
|
! ---------------------
|
|
!
|
|
! Subroutine prints module parameters.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! verbose - flag determining if the subroutine should be verbose;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine print_mesh(verbose)
|
|
|
|
use refinement, only : print_refinement
|
|
|
|
implicit none
|
|
|
|
logical, intent(in) :: verbose
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
call print_refinement(verbose)
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine print_mesh
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine GENERATE_MESH:
|
|
! ------------------------
|
|
!
|
|
! Subroutine generates the iinitial block structure by creating blocks
|
|
! according to the geometry, initial problems and refinement criterion.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine generate_mesh(status)
|
|
|
|
use blocks , only : block_meta, block_data, list_meta
|
|
use blocks , only : allocate_datablock, deallocate_datablock
|
|
use blocks , only : append_datablock
|
|
use blocks , only : link_blocks, unlink_blocks, refine_block
|
|
use blocks , only : get_mblocks, get_nleafs
|
|
use blocks , only : set_neighbors_refine
|
|
use blocks , only : build_leaf_list, build_datablock_list
|
|
#ifdef DEBUG
|
|
use blocks , only : check_neighbors
|
|
#endif /* DEBUG */
|
|
use coordinates, only : minlev, maxlev
|
|
use helpers , only : print_section, print_message
|
|
use mpitools , only : master, nproc, npmax, nodes, lprocs
|
|
use refinement , only : check_refinement_criterion
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
type(block_meta), pointer :: pmeta
|
|
type(block_data), pointer :: pdata
|
|
|
|
logical :: refine, flag
|
|
integer(kind=4) :: level, lev
|
|
integer(kind=4) :: nc, nr, np, nl
|
|
|
|
integer(kind=4), dimension(nodes) :: nb
|
|
integer(kind=4), dimension(lprocs,nodes) :: lb
|
|
|
|
character(len=*), parameter :: loc = 'MESH::generate_mesh()'
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
call setup_domain(status)
|
|
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot setup the computational domain!")
|
|
go to 100
|
|
end if
|
|
|
|
call allocate_datablock(pdata, status)
|
|
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot allocate a temporary data block!")
|
|
go to 100
|
|
end if
|
|
|
|
! at this point we assume, that the initial structure of blocks
|
|
! according to the defined geometry is already created; no refinement
|
|
! is done yet; we fill out the coarse blocks with the initial condition
|
|
!
|
|
call print_section(master, "Generating the initial mesh")
|
|
if (master) write(*,"(4x,a16,11x,'=')",advance="no") "generating level"
|
|
|
|
refine = .true.
|
|
level = 1
|
|
do while (level < maxlev .and. refine)
|
|
|
|
!! DETERMINE THE REFINEMENT OF ALL BLOCKS AT THE CURRENT LEVEL
|
|
!!
|
|
! iterate over all meta blocks at the current level, and initialize the problem
|
|
! for them using temporary data block
|
|
!
|
|
refine = .false.
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
|
|
if (pmeta%level == level) then
|
|
call link_blocks(pmeta, pdata)
|
|
call setup_problem(pdata)
|
|
|
|
! check the refinement criterion for the current block, and do not allow for
|
|
! the derefinement; if the level is lower then the minimum level set it to be
|
|
! refined;
|
|
!
|
|
pmeta%refine = max(0, check_refinement_criterion(pdata))
|
|
if (level < minlev) pmeta%refine = 1
|
|
|
|
! if there is only one block, refine it anyway because the resolution for
|
|
! the problem initiation may be too small
|
|
!
|
|
if (get_mblocks() == 1 .and. level == 1) pmeta%refine = 1
|
|
|
|
call unlink_blocks(pmeta, pdata)
|
|
|
|
! set the refine flag, if there is any block selected for refinement or
|
|
! derefinement
|
|
!
|
|
if (pmeta%refine /= 0) refine = .true.
|
|
end if ! pmeta%level == level
|
|
pmeta => pmeta%next
|
|
end do ! pmeta
|
|
|
|
! refine or derefine only if it is needed
|
|
!
|
|
if (refine) then
|
|
|
|
! print the currently processed level
|
|
!
|
|
if (master) write(*, "(1x,i0)", advance = "no") level
|
|
|
|
!! STEP DOWN AND SELECT BLOCKS WHICH NEED TO BE REFINED
|
|
!!
|
|
! walk through all levels down from the current level and check if neighbors
|
|
! of the refined blocks need to be refined as well; there is no need for
|
|
! checking the blocks at the lowest level;
|
|
!
|
|
do lev = level, 2, -1
|
|
|
|
! iterate over all meta blocks at the level lev and if the current block is
|
|
! selected for the refinement and its neighbors are at lower levels select them
|
|
! for refinement too;
|
|
!
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
|
|
if (pmeta%leaf .and. pmeta%level == lev .and. pmeta%refine == 1) &
|
|
call set_neighbors_refine(pmeta)
|
|
|
|
pmeta => pmeta%next
|
|
end do ! pmeta
|
|
|
|
end do ! lev = level, 2, -1
|
|
|
|
!! REFINE ALL BLOCKS FROM THE LOWEST LEVEL UP
|
|
!!
|
|
! walk through the levels starting from the lowest to the current level
|
|
!
|
|
do lev = 1, level
|
|
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
|
|
! check if the current block is at the level lev and refine if if
|
|
! it is selected for refinement (refine without creating new data blocks)
|
|
!
|
|
if (pmeta%level == lev .and. pmeta%refine == 1) then
|
|
|
|
call refine_block(pmeta, .false., status)
|
|
|
|
! quit if the block couldn't be refined
|
|
!
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot refine the meta block!")
|
|
call deallocate_datablock(pdata, status)
|
|
status = 1
|
|
go to 100
|
|
end if
|
|
|
|
end if
|
|
|
|
pmeta => pmeta%next
|
|
end do ! pmeta
|
|
|
|
end do ! lev = 1, level
|
|
|
|
! increase the level number
|
|
!
|
|
level = level + 1
|
|
|
|
end if ! refine
|
|
|
|
end do ! level = 1, maxlev
|
|
|
|
! print the currently processed level
|
|
!
|
|
if (master) write(*, "(1x,i0)") level
|
|
|
|
! deallocate temporary data block
|
|
!
|
|
call deallocate_datablock(pdata, status)
|
|
|
|
! check if the data block deallocation succeeded
|
|
!
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot deallocate the temporary data block!")
|
|
go to 100
|
|
end if
|
|
|
|
!! ASSOCIATE THE META BLOCKS WITH PROCESSES AND CREATE INITIAL DATA BLOCKS
|
|
!!
|
|
! divide leaf blocks between all processes, use the number of data blocks to do
|
|
! this, but keep at the same process blocks with the same parent
|
|
!
|
|
nl = mod(get_nleafs(), nodes)
|
|
nb(:) = get_nleafs() / nodes
|
|
nb(1:nl) = nb(1:nl) + 1
|
|
do nc = 1, nodes
|
|
nl = mod(nb(nc), lprocs)
|
|
lb( : ,nc) = nb(nc) / lprocs
|
|
lb(1:nl,nc) = lb(1:nl,nc) + 1
|
|
end do
|
|
|
|
! reset the rank, processor and block numbers
|
|
!
|
|
nc = 1
|
|
nr = 1
|
|
np = 0
|
|
nl = 0
|
|
|
|
! iterate over all meta blocks
|
|
!
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
|
|
! assign the process number to the current block
|
|
!
|
|
pmeta%process = np
|
|
|
|
! check if the current block is the leaf
|
|
!
|
|
if (pmeta%leaf) then
|
|
|
|
! increase the number of leafs for the current process
|
|
!
|
|
nl = nl + 1
|
|
|
|
! if the block belongs to the current process, append a new data block, link it
|
|
! to the current meta block and initialize the problem
|
|
!
|
|
if (pmeta%process == nproc) then
|
|
|
|
call append_datablock(pdata, status)
|
|
|
|
if (status == 0) then
|
|
call link_blocks(pmeta, pdata)
|
|
call setup_problem(pdata)
|
|
else
|
|
call print_message(loc, "Cannot append new a data block!")
|
|
go to 100
|
|
end if
|
|
|
|
end if ! leaf and belongs to the current process
|
|
|
|
! if the number of leafs for the current process exceeds the number of assigned
|
|
! blocks, reset the counter and increase the process number
|
|
!
|
|
if (nl >= lb(nr,nc)) then
|
|
|
|
np = min(npmax, np + 1)
|
|
nl = 0
|
|
nr = nr + 1
|
|
if (nr > lprocs) then
|
|
nr = 1
|
|
nc = nc + 1
|
|
end if
|
|
flag = nc <= nodes
|
|
if (flag) flag = lb(nr,nc) == 0
|
|
do while(flag)
|
|
np = min(npmax, np + 1)
|
|
nr = nr + 1
|
|
if (nr > lprocs) then
|
|
nr = 1
|
|
nc = nc + 1
|
|
end if
|
|
flag = nc <= nodes
|
|
if (flag) flag = lb(nr,nc) == 0
|
|
end do
|
|
|
|
end if ! nl >= lb(nr,np)
|
|
|
|
end if ! leaf
|
|
pmeta => pmeta%next
|
|
|
|
end do ! pmeta
|
|
|
|
! update the list of leafs
|
|
!
|
|
call build_leaf_list(status)
|
|
|
|
! quit if the data block allocation failed
|
|
!
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot build the list of leafs!")
|
|
go to 100
|
|
end if
|
|
|
|
#ifdef DEBUG
|
|
! check if neighbors are consistent after mesh generation
|
|
!
|
|
call check_neighbors()
|
|
#endif /* DEBUG */
|
|
|
|
! error entry point
|
|
!
|
|
100 continue
|
|
|
|
! update the vector of data blocks
|
|
!
|
|
call build_datablock_list(status)
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine generate_mesh
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine UPDATE_MESH:
|
|
! ----------------------
|
|
!
|
|
! Subroutine checks the refinement criterion for each block, and refines
|
|
! or restricts it if necessary by prolongating or restricting its data.
|
|
! In the MPI version the data blocks are redistributed among all processes
|
|
! after the mesh update.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine update_mesh(status)
|
|
|
|
use blocks, only : block_data
|
|
use blocks, only : build_leaf_list, build_datablock_list
|
|
#ifdef DEBUG
|
|
use blocks, only : check_neighbors
|
|
#endif /* DEBUG */
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! check the refinement criterion of all data blocks at the current process
|
|
!
|
|
call check_block_refinement(status)
|
|
if (status /= 0) go to 100
|
|
|
|
! update neighbor refinement flags, if they need to be refined as well
|
|
!
|
|
call update_neighbor_refinement()
|
|
|
|
! prepare siblings of blocks marked for restriction
|
|
!
|
|
call prepare_sibling_derefinement(status)
|
|
if (status /= 0) go to 100
|
|
|
|
! restrict selected blocks
|
|
!
|
|
call derefine_selected_blocks(status)
|
|
if (status /= 0) go to 100
|
|
|
|
! prolong selected blocks
|
|
!
|
|
call refine_selected_blocks(status)
|
|
if (status /= 0) go to 100
|
|
|
|
! update the list of leafs
|
|
!
|
|
call build_leaf_list(status)
|
|
if (status /= 0) go to 100
|
|
|
|
#ifdef MPI
|
|
call redistribute_blocks(status)
|
|
if (status /= 0) go to 100
|
|
#endif /* MPI */
|
|
|
|
#ifdef DEBUG
|
|
! check if neighbors are consistent after mesh refinement
|
|
!
|
|
call check_neighbors()
|
|
#endif /* DEBUG */
|
|
|
|
100 continue
|
|
|
|
! update the vector of data blocks
|
|
!
|
|
call build_datablock_list(status)
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine update_mesh
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine REDISTRIBUTE_BLOCKS:
|
|
! ------------------------------
|
|
!
|
|
! Subroutine redistributes data blocks between processes.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the call status;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine redistribute_blocks(status)
|
|
|
|
#ifdef MPI
|
|
use blocks , only : block_meta, block_data, list_meta
|
|
use blocks , only : get_nleafs, get_last_id
|
|
use blocks , only : append_datablock, remove_datablock, link_blocks
|
|
use helpers , only : print_message
|
|
use mpitools, only : check_status
|
|
use mpitools, only : npmax, nproc, nodes, lprocs
|
|
use mpitools, only : send_array, receive_array
|
|
#endif /* MPI */
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
#ifdef MPI
|
|
logical :: flag
|
|
integer(kind=4) :: np, nl, nc, nr
|
|
|
|
type(block_meta), pointer :: pmeta
|
|
type(block_data), pointer :: pdata
|
|
|
|
integer(kind=4) :: tag1, tag2
|
|
|
|
integer(kind=4), dimension(nodes) :: nb
|
|
integer(kind=4), dimension(lprocs,nodes) :: lb
|
|
|
|
character(len=*), parameter :: loc = 'MESH::redistribute_blocks()'
|
|
#endif /* MPI */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
#ifdef MPI
|
|
! calculate the new data block distribution between processes
|
|
!
|
|
nl = mod(get_nleafs(), nodes)
|
|
nb(:) = get_nleafs() / nodes
|
|
nb(1:nl) = nb(1:nl) + 1
|
|
do nc = 1, nodes
|
|
nl = mod(nb(nc), lprocs)
|
|
lb( : ,nc) = nb(nc) / lprocs
|
|
lb(1:nl,nc) = lb(1:nl,nc) + 1
|
|
end do
|
|
|
|
nc = 1
|
|
nr = 1
|
|
np = 0
|
|
nl = 0
|
|
|
|
tag2 = 0
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
|
|
! exchange data blocks between processes
|
|
!
|
|
if (pmeta%process /= np) then
|
|
if (pmeta%leaf) then
|
|
changed = .true.
|
|
|
|
tag1 = tag2 + 1
|
|
tag2 = tag1 + 1
|
|
|
|
if (nproc == pmeta%process) then
|
|
|
|
call send_array(np, tag1, pmeta%data%uu)
|
|
call send_array(np, tag2, pmeta%data%q)
|
|
|
|
call remove_datablock(pmeta%data, status)
|
|
if (status /= 0) then
|
|
call print_message(loc, "Could not remove data block!")
|
|
go to 1000
|
|
end if
|
|
|
|
end if
|
|
|
|
if (nproc == np) then
|
|
|
|
call append_datablock(pdata, status)
|
|
if (status /= 0) then
|
|
call print_message(loc, "Could not append new data block!")
|
|
go to 1000
|
|
end if
|
|
call link_blocks(pmeta, pdata)
|
|
|
|
call receive_array(pmeta%process, tag1, pmeta%data%uu)
|
|
call receive_array(pmeta%process, tag2, pmeta%data%q)
|
|
|
|
end if
|
|
end if
|
|
|
|
pmeta%process = np
|
|
|
|
end if
|
|
|
|
! determine the new block distribution
|
|
!
|
|
if (pmeta%leaf) then
|
|
nl = nl + 1
|
|
|
|
if (nl >= lb(nr,nc)) then
|
|
|
|
np = min(npmax, np + 1)
|
|
nl = 0
|
|
nr = nr + 1
|
|
if (nr > lprocs) then
|
|
nr = 1
|
|
nc = nc + 1
|
|
end if
|
|
flag = nc <= nodes
|
|
if (flag) flag = lb(nr,nc) == 0
|
|
do while(flag)
|
|
np = min(npmax, np + 1)
|
|
nr = nr + 1
|
|
if (nr > lprocs) then
|
|
nr = 1
|
|
nc = nc + 1
|
|
end if
|
|
flag = nc <= nodes
|
|
if (flag) flag = lb(nr,nc) == 0
|
|
end do
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
pmeta => pmeta%next
|
|
end do
|
|
|
|
1000 continue
|
|
if (check_status(status /= 0)) &
|
|
call print_message(loc, "Could not redistribute blocks!")
|
|
#endif /* MPI */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine redistribute_blocks
|
|
!
|
|
!===============================================================================
|
|
!!
|
|
!!*** PRIVATE SUBROUTINES ****************************************************
|
|
!!
|
|
!===============================================================================
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine SETUP_DOMAIN_DEFAULT:
|
|
! -------------------------------
|
|
!
|
|
! Subroutine sets the default domain of N₁xN₂xN₃ blocks in the right
|
|
! configuration.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine setup_domain_default(status)
|
|
|
|
use blocks , only : pointer_meta, block_meta, block_data
|
|
use blocks , only : append_metablock, append_datablock, link_blocks
|
|
use blocks , only : metablock_set_leaf, metablock_set_level
|
|
use blocks , only : metablock_set_configuration
|
|
use blocks , only : metablock_set_coordinates, metablock_set_bounds
|
|
use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen
|
|
use coordinates, only : ddims => domain_base_dims
|
|
use coordinates, only : periodic
|
|
use helpers , only : print_message
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
integer :: i, j, k, n, p, ic, jc, kc
|
|
real(kind=8) :: xl, xmn, xmx, yl, ymn, ymx, zl, zmn, zmx
|
|
|
|
integer, dimension(3) :: pos, del
|
|
|
|
type(block_meta), pointer :: pmeta
|
|
|
|
integer, dimension(:,:,:), allocatable :: cfg
|
|
integer, dimension(:) , allocatable :: im, ip
|
|
integer, dimension(:) , allocatable :: jm, jp
|
|
#if NDIMS == 3
|
|
integer, dimension(:) , allocatable :: km, kp
|
|
#endif /* NDIMS == 3 */
|
|
|
|
type(pointer_meta), dimension(:,:,:), allocatable :: block_array
|
|
|
|
character(len=*), parameter :: loc = 'DOMAINS::setup_domain_default()'
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
! obtain the number of blocks
|
|
!
|
|
n = product(ddims(1:NDIMS))
|
|
|
|
!! PREPARE BLOCK CONFIGURATION ARRAY
|
|
!!
|
|
! allocate the configuration array
|
|
!
|
|
allocate(cfg(ddims(1),ddims(2),ddims(3)), stat = status)
|
|
|
|
! check if the allocation succeeded
|
|
!
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot allocate the configuration array!")
|
|
return
|
|
end if
|
|
|
|
! set the block configurations
|
|
!
|
|
cfg(1:ddims(1),1:ddims(2):2,1:ddims(3):2) = 12
|
|
|
|
if (ddims(2) > 1) then
|
|
cfg(1:ddims(1),2:ddims(2):2,1:ddims(3):2) = 43
|
|
cfg( ddims(1),1:ddims(2) ,1:ddims(3):2) = 13
|
|
end if
|
|
|
|
if (ddims(3) > 1) then
|
|
cfg(1:ddims(1),1:ddims(2):2,2:ddims(3):2) = 65
|
|
if (ddims(2) > 1) then
|
|
cfg(1:ddims(1),2:ddims(2):2,2:ddims(3):2) = 78
|
|
cfg( ddims(1),1:ddims(2) ,2:ddims(3):2) = 75
|
|
end if
|
|
if (ddims(1) == 1 .or. mod(ddims(2),2) == 1) then
|
|
cfg( ddims(1), ddims(2) ,1:ddims(3) ) = 15
|
|
else
|
|
cfg( 1 , ddims(2) ,1:ddims(3) ) = 48
|
|
end if
|
|
end if
|
|
|
|
!! ALLOCATE AND GENERATE META BLOCK CHAIN AND SET BLOCK CONFIGURATIONS
|
|
!!
|
|
! allocate the block pointer array
|
|
!
|
|
allocate(block_array(ddims(1),ddims(2),ddims(3)), stat = status)
|
|
|
|
! check if the allocation succeeded
|
|
!
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot allocate a pointer array for domain base!")
|
|
if (allocated(cfg)) deallocate(cfg, stat = status)
|
|
status = 1
|
|
return
|
|
end if ! allocation
|
|
|
|
! generate the gray code for a given configuration and link the block in
|
|
! the proper order
|
|
!
|
|
pos(:) = (/ 0, 0, 0 /)
|
|
del(:) = (/ 1, 1, 1 /)
|
|
|
|
p = 1
|
|
do k = 1, ddims(3)
|
|
if (del(3) == 1) pos(3) = pos(3) + del(3)
|
|
do j = 1, ddims(2)
|
|
if (del(2) == 1) pos(2) = pos(2) + del(2)
|
|
do i = 1, ddims(1)
|
|
if (del(1) == 1) pos(1) = pos(1) + del(1)
|
|
|
|
! append a new metablock
|
|
!
|
|
call append_metablock(block_array(pos(1),pos(2),pos(3))%ptr, status)
|
|
|
|
! set the configuration type, if new block was appended successfully
|
|
!
|
|
if (status == 0) then
|
|
call metablock_set_configuration( &
|
|
block_array(pos(1),pos(2),pos(3))%ptr &
|
|
, cfg(pos(1),pos(2),pos(3)))
|
|
else ! status
|
|
call print_message(loc, &
|
|
"Cannot append a new metablock to the domain base!")
|
|
if (allocated(block_array)) deallocate(block_array, stat = status)
|
|
if (allocated(cfg)) deallocate(cfg , stat = status)
|
|
status = 1
|
|
return
|
|
end if ! status
|
|
|
|
p = p + 1
|
|
|
|
if (del(1) == -1) pos(1) = pos(1) + del(1)
|
|
end do
|
|
if (del(2) == -1) pos(2) = pos(2) + del(2)
|
|
del(1) = - del(1)
|
|
end do
|
|
if (del(3) == -1) pos(3) = pos(3) + del(3)
|
|
del(2) = - del(2)
|
|
end do
|
|
|
|
! deallocate the configuration array
|
|
!
|
|
deallocate(cfg, stat = status)
|
|
|
|
! check if the deallocation succeeded
|
|
!
|
|
if (status /= 0) &
|
|
call print_message(loc, "Cannot deallocate the configuration array!")
|
|
|
|
!! FILL OUT THE REMAINING FIELDS AND ALLOCATE AND ASSOCIATE DATA BLOCKS
|
|
!!
|
|
! calculate block sizes
|
|
!
|
|
xl = xlen / ddims(1)
|
|
yl = ylen / ddims(2)
|
|
zl = zlen / ddims(3)
|
|
|
|
! fill out block structure fields
|
|
!
|
|
do k = 1, ddims(3)
|
|
|
|
! claculate the block position along Z
|
|
!
|
|
kc = k - 1
|
|
|
|
! calculate the Z bounds
|
|
!
|
|
zmn = zmin + kc * zl
|
|
zmx = zmin + k * zl
|
|
|
|
do j = 1, ddims(2)
|
|
|
|
! claculate the block position along Y
|
|
!
|
|
jc = j - 1
|
|
|
|
! calculate the Y bounds
|
|
!
|
|
ymn = ymin + jc * yl
|
|
ymx = ymin + j * yl
|
|
|
|
do i = 1, ddims(1)
|
|
|
|
! claculate the block position along Y
|
|
!
|
|
ic = i - 1
|
|
|
|
! calculate the Z bounds
|
|
!
|
|
xmn = xmin + ic * xl
|
|
xmx = xmin + i * xl
|
|
|
|
! assign a pointer
|
|
!
|
|
pmeta => block_array(i,j,k)%ptr
|
|
|
|
! mark it as the leaf
|
|
!
|
|
call metablock_set_leaf(pmeta)
|
|
|
|
! set the level
|
|
!
|
|
call metablock_set_level(pmeta, 1)
|
|
|
|
! set block coordinates
|
|
!
|
|
call metablock_set_coordinates(pmeta, (/ ic, jc, kc /))
|
|
|
|
! set the bounds
|
|
!
|
|
call metablock_set_bounds(pmeta, xmn, xmx, ymn, ymx, zmn, zmx)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
!! ASSIGN THE BLOCK NEIGHBORS
|
|
!!
|
|
! allocate indices
|
|
!
|
|
#if NDIMS == 3
|
|
allocate(im(ddims(1)), ip(ddims(1)), jm(ddims(2)), jp(ddims(2)), &
|
|
km(ddims(3)), kp(ddims(3)), stat = status)
|
|
#else /* NDIMS == 3 */
|
|
allocate(im(ddims(1)), ip(ddims(1)), jm(ddims(2)), jp(ddims(2)), &
|
|
stat = status)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! check if the allocation succeeded
|
|
!
|
|
if (status == 0) then
|
|
|
|
! generate indices
|
|
!
|
|
im(:) = cshift((/(i, i = 1, ddims(1))/),-1)
|
|
ip(:) = cshift((/(i, i = 1, ddims(1))/), 1)
|
|
jm(:) = cshift((/(j, j = 1, ddims(2))/),-1)
|
|
jp(:) = cshift((/(j, j = 1, ddims(2))/), 1)
|
|
#if NDIMS == 3
|
|
km(:) = cshift((/(k, k = 1, ddims(3))/),-1)
|
|
kp(:) = cshift((/(k, k = 1, ddims(3))/), 1)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! check periodicity and reset the edge indices if box is not periodic
|
|
!
|
|
if (.not. periodic(1)) then
|
|
im( 1) = 0
|
|
ip(ddims(1)) = 0
|
|
end if
|
|
if (.not. periodic(2)) then
|
|
jm( 1) = 0
|
|
jp(ddims(2)) = 0
|
|
end if
|
|
#if NDIMS == 3
|
|
if (.not. periodic(3)) then
|
|
km( 1) = 0
|
|
kp(ddims(3)) = 0
|
|
end if
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! iterate over all initial blocks
|
|
!
|
|
do k = 1, ddims(3)
|
|
do j = 1, ddims(2)
|
|
do i = 1, ddims(1)
|
|
|
|
! assign pmeta with the current block
|
|
!
|
|
pmeta => block_array(i,j,k)%ptr
|
|
|
|
#if NDIMS == 3
|
|
! assign face neighbor pointers
|
|
!
|
|
if (im(i) > 0) then
|
|
pmeta%faces(1,1,1,1)%ptr => block_array(im(i),j,k)%ptr
|
|
pmeta%faces(1,2,1,1)%ptr => block_array(im(i),j,k)%ptr
|
|
pmeta%faces(1,1,2,1)%ptr => block_array(im(i),j,k)%ptr
|
|
pmeta%faces(1,2,2,1)%ptr => block_array(im(i),j,k)%ptr
|
|
end if
|
|
if (ip(i) > 0) then
|
|
pmeta%faces(2,1,1,1)%ptr => block_array(ip(i),j,k)%ptr
|
|
pmeta%faces(2,2,1,1)%ptr => block_array(ip(i),j,k)%ptr
|
|
pmeta%faces(2,1,2,1)%ptr => block_array(ip(i),j,k)%ptr
|
|
pmeta%faces(2,2,2,1)%ptr => block_array(ip(i),j,k)%ptr
|
|
end if
|
|
|
|
if (jm(j) > 0) then
|
|
pmeta%faces(1,1,1,2)%ptr => block_array(i,jm(j),k)%ptr
|
|
pmeta%faces(2,1,1,2)%ptr => block_array(i,jm(j),k)%ptr
|
|
pmeta%faces(1,1,2,2)%ptr => block_array(i,jm(j),k)%ptr
|
|
pmeta%faces(2,1,2,2)%ptr => block_array(i,jm(j),k)%ptr
|
|
end if
|
|
if (jp(j) > 0) then
|
|
pmeta%faces(1,2,1,2)%ptr => block_array(i,jp(j),k)%ptr
|
|
pmeta%faces(2,2,1,2)%ptr => block_array(i,jp(j),k)%ptr
|
|
pmeta%faces(1,2,2,2)%ptr => block_array(i,jp(j),k)%ptr
|
|
pmeta%faces(2,2,2,2)%ptr => block_array(i,jp(j),k)%ptr
|
|
end if
|
|
|
|
if (km(k) > 0) then
|
|
pmeta%faces(1,1,1,3)%ptr => block_array(i,j,km(k))%ptr
|
|
pmeta%faces(2,1,1,3)%ptr => block_array(i,j,km(k))%ptr
|
|
pmeta%faces(1,2,1,3)%ptr => block_array(i,j,km(k))%ptr
|
|
pmeta%faces(2,2,1,3)%ptr => block_array(i,j,km(k))%ptr
|
|
end if
|
|
if (kp(k) > 0) then
|
|
pmeta%faces(1,1,2,3)%ptr => block_array(i,j,kp(k))%ptr
|
|
pmeta%faces(2,1,2,3)%ptr => block_array(i,j,kp(k))%ptr
|
|
pmeta%faces(1,2,2,3)%ptr => block_array(i,j,kp(k))%ptr
|
|
pmeta%faces(2,2,2,3)%ptr => block_array(i,j,kp(k))%ptr
|
|
end if
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! assign edge neighbor pointers
|
|
!
|
|
#if NDIMS == 2
|
|
if (im(i) > 0) then
|
|
pmeta%edges(1,1,2)%ptr => block_array(im(i),j,k)%ptr
|
|
pmeta%edges(1,2,2)%ptr => block_array(im(i),j,k)%ptr
|
|
end if
|
|
if (ip(i) > 0) then
|
|
pmeta%edges(2,1,2)%ptr => block_array(ip(i),j,k)%ptr
|
|
pmeta%edges(2,2,2)%ptr => block_array(ip(i),j,k)%ptr
|
|
end if
|
|
if (jm(j) > 0) then
|
|
pmeta%edges(1,1,1)%ptr => block_array(i,jm(j),k)%ptr
|
|
pmeta%edges(2,1,1)%ptr => block_array(i,jm(j),k)%ptr
|
|
end if
|
|
if (jp(j) > 0) then
|
|
pmeta%edges(1,2,1)%ptr => block_array(i,jp(j),k)%ptr
|
|
pmeta%edges(2,2,1)%ptr => block_array(i,jp(j),k)%ptr
|
|
end if
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
if (jm(j) > 0 .and. km(k) > 0) then
|
|
pmeta%edges(1,1,1,1)%ptr => block_array(i,jm(j),km(k))%ptr
|
|
pmeta%edges(2,1,1,1)%ptr => block_array(i,jm(j),km(k))%ptr
|
|
end if
|
|
if (jp(j) > 0 .and. km(k) > 0) then
|
|
pmeta%edges(1,2,1,1)%ptr => block_array(i,jp(j),km(k))%ptr
|
|
pmeta%edges(2,2,1,1)%ptr => block_array(i,jp(j),km(k))%ptr
|
|
end if
|
|
if (jm(j) > 0 .and. kp(k) > 0) then
|
|
pmeta%edges(1,1,2,1)%ptr => block_array(i,jm(j),kp(k))%ptr
|
|
pmeta%edges(2,1,2,1)%ptr => block_array(i,jm(j),kp(k))%ptr
|
|
end if
|
|
if (jp(j) > 0 .and. kp(k) > 0) then
|
|
pmeta%edges(1,2,2,1)%ptr => block_array(i,jp(j),kp(k))%ptr
|
|
pmeta%edges(2,2,2,1)%ptr => block_array(i,jp(j),kp(k))%ptr
|
|
end if
|
|
|
|
if (im(i) > 0 .and. km(k) > 0) then
|
|
pmeta%edges(1,1,1,2)%ptr => block_array(im(i),j,km(k))%ptr
|
|
pmeta%edges(1,2,1,2)%ptr => block_array(im(i),j,km(k))%ptr
|
|
end if
|
|
if (ip(i) > 0 .and. km(k) > 0) then
|
|
pmeta%edges(2,1,1,2)%ptr => block_array(ip(i),j,km(k))%ptr
|
|
pmeta%edges(2,2,1,2)%ptr => block_array(ip(i),j,km(k))%ptr
|
|
end if
|
|
if (im(i) > 0 .and. kp(k) > 0) then
|
|
pmeta%edges(1,1,2,2)%ptr => block_array(im(i),j,kp(k))%ptr
|
|
pmeta%edges(1,2,2,2)%ptr => block_array(im(i),j,kp(k))%ptr
|
|
end if
|
|
if (ip(i) > 0 .and. kp(k) > 0) then
|
|
pmeta%edges(2,1,2,2)%ptr => block_array(ip(i),j,kp(k))%ptr
|
|
pmeta%edges(2,2,2,2)%ptr => block_array(ip(i),j,kp(k))%ptr
|
|
end if
|
|
|
|
if (im(i) > 0 .and. jm(j) > 0) then
|
|
pmeta%edges(1,1,1,3)%ptr => block_array(im(i),jm(j),k)%ptr
|
|
pmeta%edges(1,1,2,3)%ptr => block_array(im(i),jm(j),k)%ptr
|
|
end if
|
|
if (ip(i) > 0 .and. jm(j) > 0) then
|
|
pmeta%edges(2,1,1,3)%ptr => block_array(ip(i),jm(j),k)%ptr
|
|
pmeta%edges(2,1,2,3)%ptr => block_array(ip(i),jm(j),k)%ptr
|
|
end if
|
|
if (im(i) > 0 .and. jp(j) > 0) then
|
|
pmeta%edges(1,2,1,3)%ptr => block_array(im(i),jp(j),k)%ptr
|
|
pmeta%edges(1,2,2,3)%ptr => block_array(im(i),jp(j),k)%ptr
|
|
end if
|
|
if (ip(i) > 0 .and. jp(j) > 0) then
|
|
pmeta%edges(2,2,1,3)%ptr => block_array(ip(i),jp(j),k)%ptr
|
|
pmeta%edges(2,2,2,3)%ptr => block_array(ip(i),jp(j),k)%ptr
|
|
end if
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! assign corner neighbor pointers
|
|
!
|
|
#if NDIMS == 2
|
|
if (im(i) > 0 .and. jm(j) > 0) &
|
|
pmeta%corners(1,1)%ptr => block_array(im(i),jm(j),k)%ptr
|
|
if (ip(i) > 0 .and. jm(j) > 0) &
|
|
pmeta%corners(2,1)%ptr => block_array(ip(i),jm(j),k)%ptr
|
|
if (im(i) > 0 .and. jp(j) > 0) &
|
|
pmeta%corners(1,2)%ptr => block_array(im(i),jp(j),k)%ptr
|
|
if (ip(i) > 0 .and. jp(j) > 0) &
|
|
pmeta%corners(2,2)%ptr => block_array(ip(i),jp(j),k)%ptr
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
if (im(i) > 0 .and. jm(j) > 0 .and. km(k) > 0) &
|
|
pmeta%corners(1,1,1)%ptr => block_array(im(i),jm(j),km(k))%ptr
|
|
if (ip(i) > 0 .and. jm(j) > 0 .and. km(k) > 0) &
|
|
pmeta%corners(2,1,1)%ptr => block_array(ip(i),jm(j),km(k))%ptr
|
|
if (im(i) > 0 .and. jp(j) > 0 .and. km(k) > 0) &
|
|
pmeta%corners(1,2,1)%ptr => block_array(im(i),jp(j),km(k))%ptr
|
|
if (ip(i) > 0 .and. jp(j) > 0 .and. km(k) > 0) &
|
|
pmeta%corners(2,2,1)%ptr => block_array(ip(i),jp(j),km(k))%ptr
|
|
if (im(i) > 0 .and. jm(j) > 0 .and. kp(k) > 0) &
|
|
pmeta%corners(1,1,2)%ptr => block_array(im(i),jm(j),kp(k))%ptr
|
|
if (ip(i) > 0 .and. jm(j) > 0 .and. kp(k) > 0) &
|
|
pmeta%corners(2,1,2)%ptr => block_array(ip(i),jm(j),kp(k))%ptr
|
|
if (im(i) > 0 .and. jp(j) > 0 .and. kp(k) > 0) &
|
|
pmeta%corners(1,2,2)%ptr => block_array(im(i),jp(j),kp(k))%ptr
|
|
if (ip(i) > 0 .and. jp(j) > 0 .and. kp(k) > 0) &
|
|
pmeta%corners(2,2,2)%ptr => block_array(ip(i),jp(j),kp(k))%ptr
|
|
#endif /* NDIMS == 3 */
|
|
end do ! over i
|
|
end do ! over j
|
|
end do ! over k
|
|
|
|
! deallocate indices
|
|
!
|
|
#if NDIMS == 3
|
|
deallocate(im, ip, jm, jp, km, kp, stat = status)
|
|
#else /* NDIMS == 3 */
|
|
deallocate(im, ip, jm, jp, stat = status)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! check if the deallocation succeeded
|
|
!
|
|
if (status /= 0) &
|
|
call print_message(loc, "Cannot deallocate the configuration indices!")
|
|
|
|
else ! allocation
|
|
call print_message(loc, "Cannot allocate the configuration indices!")
|
|
status = 1
|
|
end if ! allocation
|
|
|
|
! deallocate the block pointer array
|
|
!
|
|
deallocate(block_array, stat = status)
|
|
|
|
! check if the deallocation succeeded
|
|
!
|
|
if (status /= 0) &
|
|
call print_message(loc, &
|
|
"Cannot deallocate the pointer array for domain base!")
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine setup_domain_default
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine CHECK_BLOCK_REFINEMENT:
|
|
! ---------------------------------
|
|
!
|
|
! Subroutine scans over all data blocks, gets and corrects their refinement
|
|
! flags. If the MPI is used, the refinement flags are syncronized among all
|
|
! processes.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine check_block_refinement(status)
|
|
|
|
use blocks , only : block_meta, list_meta
|
|
use blocks , only : block_data, data_blocks, get_dblocks
|
|
#ifdef MPI
|
|
use blocks , only : get_nleafs
|
|
#endif /* MPI */
|
|
use coordinates, only : minlev, maxlev
|
|
use helpers , only : print_message
|
|
#ifdef MPI
|
|
use mpitools , only : reduce_sum
|
|
#endif /* MPI */
|
|
use refinement , only : check_refinement_criterion
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
type(block_meta), pointer :: pmeta
|
|
type(block_data), pointer :: pdata
|
|
|
|
integer :: l, n
|
|
#ifdef MPI
|
|
integer(kind=4) :: nl
|
|
|
|
integer(kind=4), dimension(:), allocatable :: ibuf
|
|
|
|
character(len=*), parameter :: loc = 'MESH::check_block_refinement()'
|
|
#endif /* MPI */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
n = get_dblocks()
|
|
|
|
! 1) reset the refinement flag for all meta blocks
|
|
!
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
pmeta%refine = 0
|
|
pmeta => pmeta%next
|
|
end do ! iterate over meta blocks
|
|
|
|
! 2) determine the refinement of data block from the current process
|
|
!
|
|
!$omp parallel do default(shared) private(pdata,pmeta)
|
|
do l = 1, n
|
|
pdata => data_blocks(l)%ptr
|
|
pmeta => pdata%meta
|
|
|
|
! check the refinement criterion for the current data block
|
|
!
|
|
pmeta%refine = check_refinement_criterion(pdata)
|
|
|
|
! correct the refinement flag for the minimum and maximum levels
|
|
!
|
|
if (pmeta%level < minlev) pmeta%refine = 1
|
|
if (pmeta%level == minlev) pmeta%refine = max(0, pmeta%refine)
|
|
if (pmeta%level == maxlev) pmeta%refine = min(0, pmeta%refine)
|
|
if (pmeta%level > maxlev) pmeta%refine = -1
|
|
|
|
end do
|
|
!$omp end parallel do
|
|
|
|
#ifdef MPI
|
|
! 3) synchronize the refinement flags between all processes
|
|
!
|
|
nl = get_nleafs()
|
|
|
|
allocate(ibuf(nl), stat = status)
|
|
if (status == 0) then
|
|
ibuf(1:nl) = 0
|
|
l = 0
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
if (pmeta%leaf) then
|
|
l = l + 1
|
|
|
|
ibuf(l) = pmeta%refine
|
|
end if
|
|
|
|
pmeta => pmeta%next
|
|
end do
|
|
|
|
call reduce_sum(ibuf(1:nl))
|
|
|
|
l = 0
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
if (pmeta%leaf) then
|
|
l = l + 1
|
|
|
|
pmeta%refine = ibuf(l)
|
|
|
|
end if
|
|
|
|
pmeta => pmeta%next
|
|
end do
|
|
|
|
deallocate(ibuf, stat = status)
|
|
|
|
if (status /= 0) &
|
|
call print_message(loc, &
|
|
"Refinement flag buffer could not be deallocated!")
|
|
|
|
else
|
|
call print_message(loc, "Refinement flag buffer could not be allocated!")
|
|
end if
|
|
#endif /* MPI */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine check_block_refinement
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine UPDATE_NEIGHBOR_REFINEMENT:
|
|
! -------------------------------------
|
|
!
|
|
! Subroutine scans over all neighbors of blocks selected for refinement or
|
|
! derefinement, and if necessary selects them to be refined as well, or
|
|
! cancels their derefinement.
|
|
!
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine update_neighbor_refinement()
|
|
|
|
! import external procedures and variables
|
|
!
|
|
use blocks , only : block_meta, list_meta
|
|
use blocks , only : set_neighbors_refine
|
|
use coordinates , only : toplev
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! local pointers
|
|
!
|
|
type(block_meta), pointer :: pmeta
|
|
|
|
! local variables
|
|
!
|
|
integer(kind=4) :: l
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! iterate down over all levels and correct the refinement of neighbor blocks
|
|
!
|
|
do l = toplev, 1, -1
|
|
|
|
! assign pmeta to the first meta block on the list
|
|
!
|
|
pmeta => list_meta
|
|
|
|
! iterate over all meta blocks
|
|
!
|
|
do while (associated(pmeta))
|
|
|
|
! check only leafs at the current level
|
|
!
|
|
if (pmeta%leaf .and. pmeta%level == l) then
|
|
|
|
! correct neighbor refinement flags
|
|
!
|
|
call set_neighbors_refine(pmeta)
|
|
|
|
end if ! the leaf at level l
|
|
|
|
! assign pmeta to the next meta block
|
|
!
|
|
pmeta => pmeta%next
|
|
|
|
end do ! iterate over meta blocks
|
|
|
|
end do ! levels
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine update_neighbor_refinement
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine PREPARE_SIBLING_DEREFINEMENT:
|
|
! ---------------------------------------
|
|
!
|
|
! Subroutine scans over all blocks selected for derefinement and checks if
|
|
! their siblings can be derefined as well. If any of them cannot be
|
|
! derefined, the derefinement of all siblings is canceled. Then, if MPI is
|
|
! used, the subroutine brings back all siblings together to lay on
|
|
! the same process.
|
|
!
|
|
! Note: This subroutine sets %refine flag of the parent to -1 to let the next
|
|
! executed subroutine derefine_selected_blocks() which parent block
|
|
! has to be derefined. That subroutine resets %refine flag of the
|
|
! parent after performing full restriction.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine prepare_sibling_derefinement(status)
|
|
|
|
use blocks , only : block_meta, list_meta
|
|
#ifdef MPI
|
|
use blocks , only : block_data
|
|
#endif /* MPI */
|
|
use blocks , only : nchildren
|
|
use blocks , only : set_neighbors_refine
|
|
#ifdef MPI
|
|
use blocks , only : append_datablock, remove_datablock, link_blocks
|
|
#endif /* MPI */
|
|
use coordinates, only : toplev
|
|
use helpers , only : print_message
|
|
#ifdef MPI
|
|
use mpitools , only : nproc
|
|
use mpitools , only : send_array, receive_array
|
|
#endif /* MPI */
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
type(block_meta), pointer :: pmeta, pparent, pchild
|
|
#ifdef MPI
|
|
type(block_data), pointer :: pdata
|
|
#endif /* MPI */
|
|
|
|
integer(kind=4) :: l, p
|
|
#ifdef MPI
|
|
integer(kind=4) :: tag
|
|
#endif /* MPI */
|
|
|
|
logical, dimension(nchildren) :: flag
|
|
|
|
#ifdef MPI
|
|
character(len=*), parameter :: loc = 'MESH::prepare_sibling_derefinement()'
|
|
#endif /* MPI */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
! 1) the block can be derefined only if all its children are selected for
|
|
! the derefined; otherwise reset the derefinement both for the parent
|
|
! and children;
|
|
!
|
|
do l = 2, toplev
|
|
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
|
|
if (pmeta%leaf .and. pmeta%level == l) then
|
|
if (pmeta%refine == -1) then
|
|
pparent => pmeta%parent
|
|
|
|
do p = 1, nchildren
|
|
pchild => pparent%child(p)%ptr
|
|
|
|
flag(p) = pchild%leaf .and. pchild%refine == -1
|
|
end do
|
|
|
|
if (all(flag)) then
|
|
pparent%refine = -1
|
|
else
|
|
do p = 1, nchildren
|
|
pchild => pparent%child(p)%ptr
|
|
|
|
pchild%refine = max(0, pchild%refine)
|
|
end do
|
|
end if
|
|
end if
|
|
end if
|
|
|
|
pmeta => pmeta%next
|
|
end do
|
|
|
|
end do
|
|
|
|
#ifdef MPI
|
|
! 2) bring all siblings to the same process, so we can merge them into a newly
|
|
! created data block of the parent;
|
|
!
|
|
tag = 0
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
|
|
if (.not. pmeta%leaf) then
|
|
|
|
if (pmeta%refine == -1) then
|
|
|
|
pchild => pmeta%child(1)%ptr
|
|
|
|
pmeta%process = pchild%process
|
|
|
|
do p = 2, nchildren
|
|
pchild => pmeta%child(p)%ptr
|
|
|
|
if (pchild%process /= pmeta%process) then
|
|
|
|
tag = tag + 1
|
|
|
|
if (pchild%process == nproc) then
|
|
call send_array(pmeta%process, tag, pchild%data%uu)
|
|
|
|
call remove_datablock(pchild%data, status)
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot remove the data block!")
|
|
go to 100
|
|
end if
|
|
end if
|
|
|
|
if (pmeta%process == nproc) then
|
|
call append_datablock(pdata, status)
|
|
if (status == 0) then
|
|
call link_blocks(pchild, pdata)
|
|
else
|
|
call print_message(loc, "Cannot append a new data block!")
|
|
go to 100
|
|
end if
|
|
|
|
call receive_array(pchild%process, tag, pdata%uu)
|
|
end if
|
|
|
|
pchild%process = pmeta%process
|
|
|
|
end if
|
|
end do
|
|
|
|
end if
|
|
|
|
end if
|
|
|
|
pmeta => pmeta%next
|
|
end do
|
|
|
|
100 continue
|
|
#endif /* MPI */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine prepare_sibling_derefinement
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine DEREFINE_SELECTED_BLOCKS:
|
|
! -----------------------------------
|
|
!
|
|
! Subroutine scans over all blocks and restrict those selected.
|
|
!
|
|
! Note: This subroutine resets the flag %refine set in subroutine
|
|
! prepare_sibling_derefinement().
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine derefine_selected_blocks(status)
|
|
|
|
use blocks , only : block_meta, block_data, list_meta
|
|
use blocks , only : append_datablock, link_blocks, derefine_block
|
|
use coordinates, only : toplev
|
|
use helpers , only : print_message
|
|
#ifdef MPI
|
|
use mpitools , only : nproc
|
|
#endif /* MPI */
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
type(block_meta), pointer :: pmeta
|
|
type(block_data), pointer :: pdata
|
|
|
|
integer(kind=4) :: l
|
|
|
|
character(len=*), parameter :: loc = 'MESH::derefine_selected_blocks()'
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! reset the status flag
|
|
!
|
|
status = 0
|
|
|
|
! iterate over levels and restrict the blocks selected for restriction
|
|
!
|
|
do l = toplev - 1, 1, -1
|
|
|
|
! iterate over all meta blocks
|
|
!
|
|
pmeta => list_meta
|
|
do while (associated(pmeta))
|
|
|
|
! process non-leafs at the current level selected for restriction
|
|
!
|
|
if (.not. pmeta%leaf .and. pmeta%level == l &
|
|
.and. pmeta%refine == -1) then
|
|
|
|
! indicate that the block structure will change
|
|
!
|
|
changed = .true.
|
|
|
|
#ifdef MPI
|
|
! check if pmeta belongs to the current process
|
|
!
|
|
if (pmeta%process == nproc) then
|
|
#endif /* MPI */
|
|
|
|
! check if a data block is associated with pmeta, if not create one
|
|
!
|
|
if (.not. associated(pmeta%data)) then
|
|
|
|
! append new data block and link it with the current pmeta
|
|
!
|
|
call append_datablock(pdata, status)
|
|
|
|
if (status == 0) then
|
|
call link_blocks(pmeta, pdata)
|
|
else
|
|
call print_message(loc, "Cannot append a new data block!")
|
|
go to 100
|
|
end if
|
|
|
|
end if ! no data block associated
|
|
|
|
! perform the block restriction
|
|
!
|
|
call restrict_block(pmeta)
|
|
|
|
#ifdef MPI
|
|
end if ! pmeta belongs to the current process
|
|
#endif /* MPI */
|
|
|
|
! perform the mesh derefinement, and reset the refinement flag of
|
|
! the current block
|
|
!
|
|
call derefine_block(pmeta, status)
|
|
|
|
if (status == 0) then
|
|
pmeta%refine = 0
|
|
else
|
|
call print_message(loc, "Cannot derefine the current meta block!")
|
|
go to 100
|
|
end if
|
|
|
|
end if ! non-leaf at current level selected for derefinement
|
|
|
|
pmeta => pmeta%next
|
|
end do ! iterate over meta blocks
|
|
|
|
end do ! levels
|
|
|
|
! error entry level
|
|
!
|
|
100 continue
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine derefine_selected_blocks
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine REFINE_SELECTED_BLOCKS:
|
|
! ---------------------------------
|
|
!
|
|
! Subroutine scans over all blocks and prolongates those selected.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine refine_selected_blocks(status)
|
|
|
|
use blocks , only : block_meta, list_meta
|
|
use blocks , only : refine_block, remove_datablock
|
|
use coordinates, only : toplev
|
|
use helpers , only : print_message
|
|
#ifdef MPI
|
|
use mpitools , only : nproc
|
|
#endif /* MPI */
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
type(block_meta), pointer :: pmeta, pparent
|
|
|
|
integer(kind=4) :: l
|
|
|
|
character(len=*), parameter :: loc = 'MESH::refine_selected_blocks()'
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! reset the status flag
|
|
!
|
|
status = 0
|
|
|
|
! iterate over all levels and prolong those selected for prolongation
|
|
!
|
|
do l = 1, toplev - 1
|
|
|
|
! assign pmeta to the first meta block on the list
|
|
!
|
|
pmeta => list_meta
|
|
|
|
! iterate over all meta blocks
|
|
!
|
|
do while (associated(pmeta))
|
|
|
|
! process leafs at the current level selected for prolongation
|
|
!
|
|
if (pmeta%leaf .and. pmeta%level == l .and. pmeta%refine == 1) then
|
|
|
|
! indicate that the block structure will change
|
|
!
|
|
changed = .true.
|
|
|
|
! assign pparent with the new parent block
|
|
!
|
|
pparent => pmeta
|
|
|
|
#ifdef MPI
|
|
! check if pmeta belongs to the current process
|
|
!
|
|
if (pmeta%process == nproc) then
|
|
#endif /* MPI */
|
|
|
|
! prepare child blocks with allocating the data blocks
|
|
!
|
|
call refine_block(pmeta, .true., status)
|
|
|
|
! perform the data prolongation
|
|
!
|
|
call prolong_block(pparent, status)
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot prolong the current data block!")
|
|
go to 100
|
|
end if
|
|
|
|
! remove the data block associated with the new parent
|
|
!
|
|
call remove_datablock(pparent%data, status)
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot remove the current data block!")
|
|
go to 100
|
|
end if
|
|
|
|
#ifdef MPI
|
|
else ! pmeta belongs to the current process
|
|
|
|
! prepare child blocks without actually allocating the data blocks
|
|
!
|
|
call refine_block(pmeta, .false., status)
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot refine the current meta block!")
|
|
go to 100
|
|
end if
|
|
|
|
end if ! pmeta belongs to the current process
|
|
|
|
! redistribute blocks equally among all processors
|
|
!
|
|
call redistribute_blocks(status)
|
|
if (status /= 0) then
|
|
call print_message(loc, "Cannot redistribute blocks!")
|
|
go to 100
|
|
end if
|
|
#endif /* MPI */
|
|
|
|
end if ! leaf at current level selected for refinement
|
|
|
|
! assign pmeta to the next meta block
|
|
!
|
|
pmeta => pmeta%next
|
|
|
|
end do ! iterate over meta blocks
|
|
|
|
end do ! levels
|
|
|
|
! error entry level
|
|
!
|
|
100 continue
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine refine_selected_blocks
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine PROLONG_BLOCK:
|
|
! ------------------------
|
|
!
|
|
! Subroutine prolongs variables from a data blocks linked to the input
|
|
! meta block and copy the resulting array of variables to
|
|
! the newly created children data block. The process of data restriction
|
|
! conserves stored variables.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pmeta - the input meta block;
|
|
! status - the subroutine call status: 0 for success, otherwise failure;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine prolong_block(pmeta, status)
|
|
|
|
use blocks , only : block_meta, block_data, nchildren
|
|
use coordinates , only : ni => ncells, nn => bcells
|
|
use coordinates , only : nh => nghosts_half
|
|
use coordinates , only : nb, ne
|
|
use equations , only : nv, positive
|
|
use helpers , only : print_message
|
|
use interpolations, only : limiter_prol
|
|
use workspace , only : resize_workspace, work, work_in_use
|
|
|
|
implicit none
|
|
|
|
type(block_meta), pointer, intent(inout) :: pmeta
|
|
integer , intent(out) :: status
|
|
|
|
type(block_meta), pointer :: pchild
|
|
type(block_data), pointer :: pdata
|
|
|
|
logical, save :: first = .true.
|
|
|
|
integer :: n, p, nl, nu
|
|
integer :: i, im, ip
|
|
integer :: j, jm, jp
|
|
#if NDIMS == 3
|
|
integer :: k, km, kp
|
|
#else /* NDIMS == 3 */
|
|
integer :: k
|
|
#endif /* NDIMS == 3 */
|
|
real(kind=8) :: dul, dur
|
|
#if NDIMS == 3
|
|
real(kind=8) :: du1, du2, du3, du4
|
|
#else /* NDIMS == 3 */
|
|
real(kind=8) :: du1, du2
|
|
#endif /* NDIMS == 3 */
|
|
|
|
character(len=80) :: msg
|
|
|
|
integer , dimension(NDIMS) :: l, u
|
|
real(kind=8), dimension(NDIMS) :: du
|
|
|
|
real(kind=8), dimension(:,:,:), pointer, save :: tmp
|
|
|
|
integer, save :: nt
|
|
!$ integer :: omp_get_thread_num
|
|
!$omp threadprivate(first, nt, tmp)
|
|
|
|
character(len=*), parameter :: loc = 'MESH::prolong_block()'
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
nt = 0
|
|
!$ nt = omp_get_thread_num()
|
|
status = 0
|
|
|
|
if (first) then
|
|
n = product(pm(:))
|
|
|
|
call resize_workspace(n, status)
|
|
if (status /= 0) then
|
|
call print_message(loc, "Could not resize the workspace!")
|
|
go to 100
|
|
end if
|
|
|
|
tmp(1:pm(1),1:pm(2),1:pm(3)) => work(1:n,nt)
|
|
|
|
first = .false.
|
|
end if
|
|
|
|
if (work_in_use(nt)) &
|
|
call print_message(loc, "Workspace is being used right now! " // &
|
|
"Corruptions can occur!")
|
|
|
|
work_in_use(nt) = .true.
|
|
|
|
#if NDIMS == 2
|
|
k = 1
|
|
#endif /* NDIMS == 2 */
|
|
|
|
nl = nb - nh
|
|
nu = ne + nh
|
|
|
|
pdata => pmeta%data
|
|
|
|
do n = 1, nv
|
|
|
|
#if NDIMS == 3
|
|
do k = nl, nu
|
|
km = k - 1
|
|
kp = k + 1
|
|
l(3) = 2 * (k - nl) + 1
|
|
u(3) = l(3) + 1
|
|
#endif /* NDIMS == 3 */
|
|
do j = nl, nu
|
|
jm = j - 1
|
|
jp = j + 1
|
|
l(2) = 2 * (j - nl) + 1
|
|
u(2) = l(2) + 1
|
|
do i = nl, nu
|
|
im = i - 1
|
|
ip = i + 1
|
|
l(1) = 2 * (i - nl) + 1
|
|
u(1) = l(1) + 1
|
|
|
|
dul = pdata%u(n,i ,j,k) - pdata%u(n,im,j,k)
|
|
dur = pdata%u(n,ip,j,k) - pdata%u(n,i ,j,k)
|
|
du(1) = limiter_prol(2.5d-01, dul, dur)
|
|
|
|
dul = pdata%u(n,i,j ,k) - pdata%u(n,i,jm,k)
|
|
dur = pdata%u(n,i,jp,k) - pdata%u(n,i,j ,k)
|
|
du(2) = limiter_prol(2.5d-01, dul, dur)
|
|
|
|
#if NDIMS == 3
|
|
dul = pdata%u(n,i,j,k ) - pdata%u(n,i,j,km)
|
|
dur = pdata%u(n,i,j,kp) - pdata%u(n,i,j,k )
|
|
du(3) = limiter_prol(2.5d-01, dul, dur)
|
|
#endif /* NDIMS == 3 */
|
|
|
|
if (positive(n) .and. pdata%u(n,i,j,k) < sum(abs(du(1:NDIMS)))) then
|
|
if (pdata%u(n,i,j,k) > 0.0d+00) then
|
|
do while (pdata%u(n,i,j,k) <= sum(abs(du(1:NDIMS))))
|
|
du(:) = 0.5d+00 * du(:)
|
|
end do
|
|
else
|
|
write(msg,"(a,3i4,a)") &
|
|
"Positive variable is not positive at (", i, j, k, " )"
|
|
call print_message(loc, msg)
|
|
du(:) = 0.0d+00
|
|
end if
|
|
end if
|
|
|
|
#if NDIMS == 2
|
|
du1 = du(1) + du(2)
|
|
du2 = du(1) - du(2)
|
|
tmp(l(1),l(2),k) = pdata%u(n,i,j,k) - du1
|
|
tmp(u(1),l(2),k) = pdata%u(n,i,j,k) + du2
|
|
tmp(l(1),u(2),k) = pdata%u(n,i,j,k) - du2
|
|
tmp(u(1),u(2),k) = pdata%u(n,i,j,k) + du1
|
|
#endif /* NDIMS == 2 */
|
|
|
|
#if NDIMS == 3
|
|
du1 = du(1) + du(2) + du(3)
|
|
du2 = du(1) - du(2) - du(3)
|
|
du3 = du(1) - du(2) + du(3)
|
|
du4 = du(1) + du(2) - du(3)
|
|
tmp(l(1),l(2),l(3)) = pdata%u(n,i,j,k) - du1
|
|
tmp(u(1),l(2),l(3)) = pdata%u(n,i,j,k) + du2
|
|
tmp(l(1),u(2),l(3)) = pdata%u(n,i,j,k) - du3
|
|
tmp(u(1),u(2),l(3)) = pdata%u(n,i,j,k) + du4
|
|
tmp(l(1),l(2),u(3)) = pdata%u(n,i,j,k) - du4
|
|
tmp(u(1),l(2),u(3)) = pdata%u(n,i,j,k) + du3
|
|
tmp(l(1),u(2),u(3)) = pdata%u(n,i,j,k) - du2
|
|
tmp(u(1),u(2),u(3)) = pdata%u(n,i,j,k) + du1
|
|
#endif /* NDIMS == 3 */
|
|
end do
|
|
end do
|
|
#if NDIMS == 3
|
|
end do
|
|
#endif /* NDIMS == 3 */
|
|
|
|
do p = 1, nchildren
|
|
|
|
pchild => pmeta%child(p)%ptr
|
|
|
|
l(1:NDIMS) = pchild%pos(1:NDIMS) * ni + 1
|
|
u(1:NDIMS) = l(1:NDIMS) + nn - 1
|
|
|
|
#if NDIMS == 2
|
|
pchild%data%u(n,:,:,1) = tmp(l(1):u(1),l(2):u(2),k)
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
pchild%data%u(n,:,:,:) = tmp(l(1):u(1),l(2):u(2),l(3):u(3))
|
|
#endif /* NDIMS == 3 */
|
|
|
|
end do ! nchildren
|
|
end do ! n = 1, nv
|
|
|
|
work_in_use(nt) = .false.
|
|
|
|
100 continue
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine prolong_block
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine RESTRICT_BLOCK:
|
|
! -------------------------
|
|
!
|
|
! Subroutine restricts variables from children data blocks linked to
|
|
! the input meta block and copy the resulting array of variables to
|
|
! the associated data block. The process of data restriction conserves
|
|
! stored variables.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! pmeta - the input meta block;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine restrict_block(pmeta)
|
|
|
|
! import external procedures and variables
|
|
!
|
|
use blocks , only : block_meta, block_data, nchildren
|
|
use coordinates , only : nn => bcells, nh => bcells_half
|
|
use coordinates , only : ng => nghosts_half
|
|
use coordinates , only : nb, ne
|
|
use equations , only : nv
|
|
|
|
! local variables are not implicit by default
|
|
!
|
|
implicit none
|
|
|
|
! subroutine arguments
|
|
!
|
|
type(block_meta), pointer, intent(inout) :: pmeta
|
|
|
|
! local variables
|
|
!
|
|
integer :: p
|
|
integer :: il, iu, ip, is, it
|
|
integer :: jl, ju, jp, js, jt
|
|
#if NDIMS == 3
|
|
integer :: kl, ku, kp, ks, kt
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! local arrays
|
|
!
|
|
integer, dimension(NDIMS) :: pos
|
|
|
|
! local pointers
|
|
!
|
|
type(block_data), pointer :: pparent, pchild
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
! assign the parent data pointer
|
|
!
|
|
pparent => pmeta%data
|
|
|
|
! iterate over all children
|
|
!
|
|
do p = 1, nchildren
|
|
|
|
! assign a pointer to the current child
|
|
!
|
|
pchild => pmeta%child(p)%ptr%data
|
|
|
|
! obtain the child position in the parent block
|
|
!
|
|
pos(1:NDIMS) = pchild%meta%pos(1:NDIMS)
|
|
|
|
! calculate the bound indices of the source nad destination arrays
|
|
!
|
|
if (pos(1) == 0) then
|
|
il = 1
|
|
iu = ne
|
|
is = nb - ng
|
|
it = nh
|
|
else
|
|
il = nb
|
|
iu = nn
|
|
is = nh + 1
|
|
it = ne + ng
|
|
end if
|
|
ip = il + 1
|
|
if (pos(2) == 0) then
|
|
jl = 1
|
|
ju = ne
|
|
js = nb - ng
|
|
jt = nh
|
|
else
|
|
jl = nb
|
|
ju = nn
|
|
js = nh + 1
|
|
jt = ne + ng
|
|
end if
|
|
jp = jl + 1
|
|
#if NDIMS == 3
|
|
if (pos(3) == 0) then
|
|
kl = 1
|
|
ku = ne
|
|
ks = nb - ng
|
|
kt = nh
|
|
else
|
|
kl = nb
|
|
ku = nn
|
|
ks = nh + 1
|
|
kt = ne + ng
|
|
end if
|
|
kp = kl + 1
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! restrict conserved variables from the current child and copy the resulting
|
|
! array to the proper location of the parent data block
|
|
!
|
|
#if NDIMS == 2
|
|
pparent%u(1:nv,is:it,js:jt, : ) = &
|
|
2.50d-01 * ((pchild%u(1:nv,il:iu:2,jl:ju:2, : ) &
|
|
+ pchild%u(1:nv,ip:iu:2,jp:ju:2, : )) &
|
|
+ (pchild%u(1:nv,il:iu:2,jp:ju:2, : ) &
|
|
+ pchild%u(1:nv,ip:iu:2,jl:ju:2, : )))
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
pparent%u(1:nv,is:it,js:jt,ks:kt) = &
|
|
1.25d-01 * (((pchild%u(1:nv,il:iu:2,jl:ju:2,kl:ku:2) &
|
|
+ pchild%u(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) &
|
|
+ (pchild%u(1:nv,il:iu:2,jl:ju:2,kp:ku:2) &
|
|
+ pchild%u(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) &
|
|
+ ((pchild%u(1:nv,il:iu:2,jp:ju:2,kp:ku:2) &
|
|
+ pchild%u(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) &
|
|
+ (pchild%u(1:nv,il:iu:2,jp:ju:2,kl:ku:2) &
|
|
+ pchild%u(1:nv,ip:iu:2,jl:ju:2,kp:ku:2))))
|
|
#endif /* NDIMS == 3 */
|
|
|
|
end do ! p = 1, nchildren
|
|
|
|
! prepare indices for the border filling
|
|
!
|
|
is = nb - ng
|
|
it = ne + ng
|
|
js = nb - ng
|
|
jt = ne + ng
|
|
#if NDIMS == 3
|
|
ks = nb - ng
|
|
kt = ne + ng
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! fill out the borders
|
|
!
|
|
#if NDIMS == 2
|
|
do il = is - 1, 1, -1
|
|
pparent%u(1:nv,il,js:jt, : ) = pparent%u(1:nv,is,js:jt, : )
|
|
end do
|
|
do iu = it + 1, nn
|
|
pparent%u(1:nv,iu,js:jt, : ) = pparent%u(1:nv,it,js:jt, : )
|
|
end do
|
|
do jl = js - 1, 1, -1
|
|
pparent%u(1:nv, : ,jl, : ) = pparent%u(1:nv, : ,js, : )
|
|
end do
|
|
do ju = jt + 1, nn
|
|
pparent%u(1:nv, : ,ju, : ) = pparent%u(1:nv, : ,jt, : )
|
|
end do
|
|
#endif /* NDIMS == 2 */
|
|
#if NDIMS == 3
|
|
do il = is - 1, 1, -1
|
|
pparent%u(1:nv,il,js:jt,ks:kt) = pparent%u(1:nv,is,js:jt,ks:kt)
|
|
end do
|
|
do iu = it + 1, nn
|
|
pparent%u(1:nv,iu,js:jt,ks:kt) = pparent%u(1:nv,it,js:jt,ks:kt)
|
|
end do
|
|
do jl = js - 1, 1, -1
|
|
pparent%u(1:nv, : ,jl,ks:kt) = pparent%u(1:nv, : ,js,ks:kt)
|
|
end do
|
|
do ju = jt + 1, nn
|
|
pparent%u(1:nv, : ,ju,ks:kt) = pparent%u(1:nv, : ,jt,ks:kt)
|
|
end do
|
|
do kl = ks - 1, 1, -1
|
|
pparent%u(1:nv, : , : ,kl) = pparent%u(1:nv, : , : ,ks)
|
|
end do
|
|
do ku = kt + 1, nn
|
|
pparent%u(1:nv, : , : ,ku) = pparent%u(1:nv, : , : ,kt)
|
|
end do
|
|
#endif /* NDIMS == 3 */
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine restrict_block
|
|
|
|
!===============================================================================
|
|
!
|
|
end module
|