!!******************************************************************************
!!
!! 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-2020 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
#ifdef PROFILE
! import external subroutines
!
use timers, only : set_timer, start_timer, stop_timer
#endif /* PROFILE */
! module variables are not implicit by default
!
implicit none
#ifdef PROFILE
! timer indices
!
integer , save :: imi, ims, img, imu, ima, imp, imr
#endif /* PROFILE */
! file handler for the mesh statistics
!
integer, save :: funit = 11
! allocatable array for prolongation
!
real(kind=8), dimension(:,:,:,:), allocatable :: up
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_mesh, finalize_mesh
public :: generate_mesh, update_mesh
public :: redistribute_blocks, store_mesh_stats
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
contains
!
!===============================================================================
!!
!!*** PUBLIC SUBROUTINES *****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INITIALIZE_MESH:
! --------------------------
!
! Subroutine initializes mesh module and its variables.
!
! Arguments:
!
! nrun - the restarted execution count;
! status - the subroutine call status: 0 for success, otherwise failure;
!
!===============================================================================
!
subroutine initialize_mesh(nrun, status)
! import external procedures and variables
!
use coordinates , only : ni => ncells, ng => nghosts, toplev
use equations , only : nv
use iso_fortran_env, only : error_unit
use mpitools , only : master, nprocs
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(in) :: nrun
integer, intent(out) :: status
! local variables
!
character(len=64) :: fn
integer :: l, n
! local arrays
!
integer, dimension(3) :: pm = 1
! local parameters
!
character(len=*), parameter :: loc = 'MESH::initialize_mesh()'
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! set timer descriptions
!
call set_timer('mesh:: initialization' , imi)
call set_timer('mesh:: statistics' , ims)
call set_timer('mesh:: initial generation', img)
call set_timer('mesh:: adaptive update' , imu)
call set_timer('mesh:: autobalancing' , ima)
call set_timer('mesh:: block restriction' , imr)
call set_timer('mesh:: block prolongation', imp)
! start accounting time for module initialization/finalization
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! only master prepares the mesh statistics file
!
if (master) then
! generate the mesh statistics file name
!
write(fn, "('mesh_',i2.2,'.dat')") nrun
! create a new mesh statistics file
!
#ifdef __INTEL_COMPILER
open(unit = funit, file = fn, form = 'formatted', status = 'replace' &
, buffered = 'yes')
#else /* __INTEL_COMPILER */
open(unit = funit, file = fn, form = 'formatted', status = 'replace')
#endif /* __INTEL_COMPILER */
! write the mesh statistics header
!
write(funit, "('#',2(4x,a4),10x,a8,6x,a10,5x,a6,6x,a5,4x,a20)") &
'step', 'time', 'coverage', 'efficiency' &
, 'blocks', 'leafs', 'block distributions:'
! write the mesh distributions header
!
write(funit, "('#',76x,a8)", advance="no") 'level = '
do l = 1, toplev
write(funit, "(1x,i9)", advance="no") l
end do
#ifdef MPI
write(funit, "(4x,a6)", advance="no") 'cpu = '
do n = 1, nprocs
write(funit, "(1x,i9)", advance="no") n
end do
#endif /* MPI */
write(funit, "('' )")
write(funit, "('#')")
end if ! master
! prepare dimensions
!
pm(1:NDIMS) = 2 * (ni + ng)
! allocate array for the prolongation array
!
allocate(up(nv, pm(1), pm(2), pm(3)), stat = status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Prolongation array could not be allocated!"
end if
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
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)
! import external procedures and variables
!
use iso_fortran_env, only : error_unit
use mpitools , only : master
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
! local parameters
!
character(len=*), parameter :: loc = 'MESH::finalize_mesh()'
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for module initialization/finalization
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! only master posses a handler to the mesh statistics file
!
if (master) close(funit)
! deallocate prolongation array
!
if (allocated(up)) then
deallocate(up, stat = status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Prolongation array could not be deallocated!"
end if
end if
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine finalize_mesh
!
!===============================================================================
!
! subroutine STORE_MESH_STATS:
! ---------------------------
!
! Subroutine calculates and stores various mesh statistics.
!
! Arguments:
!
! step - the integration iteration step number;
! time - the physical time;
!
!===============================================================================
!
subroutine store_mesh_stats(step, time)
! import external procedures and variables
!
use blocks , only : block_meta, block_leaf
use blocks , only : list_meta, list_leaf
use blocks , only : ndims
use blocks , only : get_mblocks, get_nleafs
use coordinates , only : ddims => domain_base_dims
use coordinates , only : ni => ncells, nn => bcells
use coordinates , only : nd => nghosts_double
use coordinates , only : toplev
use mpitools , only : master, nprocs
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer , intent(in) :: step
real(kind=8), intent(in) :: time
! local variables
!
integer(kind=4) :: l, n
real(kind=8) :: cv, ef, ff
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_leaf), pointer :: pleaf
! local saved variables
!
logical , save :: first = .true.
integer(kind=4), save :: nm = 0, nl = 0
real(kind=8) , save :: fcv = 1.0d+00, fef = 1.0d+00
! local arrays
!
integer(kind=4), dimension(toplev) :: ldist
#ifdef MPI
integer(kind=4), dimension(nprocs) :: cdist
#endif /* MPI */
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the mesh statistics
!
call start_timer(ims)
#endif /* PROFILE */
! process and store the mesh statistics only on master node
!
if (master) then
! and only if the number of blocks changed
!
if (nm /= get_mblocks() .or. nl /= get_nleafs()) then
! get the maximum number of block and level, calculated only once
!
if (first) then
! determine the number of base blocks
!
n = 0
pmeta => list_meta
do while(associated(pmeta))
if (pmeta%level == 1) n = n + 1
pmeta => pmeta%next
end do ! pmeta
! calculate the maximum number of blocks
!
n = n * 2**(ndims*(toplev - 1))
! calculate the coverage and efficiency factors
!
ff = 2**(toplev - 1)
fcv = 1.0d+00 / n
fef = 1.0d+00 * (ddims(1) * ni * ff + nd) / nn
fef = fef * (ddims(2) * ni * ff + nd) / nn
#if NDIMS == 3
fef = fef * (ddims(3) * ni * ff + nd) / nn
#endif /* NDIMS == 3 */
! reset the first execution flag
!
first = .false.
end if ! first
! get the new numbers of meta blocks and leafs
!
nm = get_mblocks()
nl = get_nleafs()
! calculate the coverage (the number of leafs divided by the maximum number
! of blocks) and the efficiency (the total number of cells in a corresponding
! uniform mesh divided by the number of cells for adaptive mesh case)
!
cv = fcv * nl
ef = fef / nl
! get the block level and block process distributions
!
ldist(:) = 0
#ifdef MPI
cdist(:) = 0
#endif /* MPI */
pleaf => list_leaf
do while(associated(pleaf))
pmeta => pleaf%meta
ldist(pmeta%level) = ldist(pmeta%level) + 1
#ifdef MPI
cdist(pmeta%process+1) = cdist(pmeta%process+1) + 1
#endif /* MPI */
pleaf => pleaf%next
end do ! pleaf
! store the block statistics
!
write(funit, "(i9,3es14.6e3,2(2x,i9))", advance="no") &
step, time, cv, ef, nm, nl
! store the block level distribution
!
write(funit,"(12x)", advance="no")
do l = 1, toplev
write(funit,"(1x,i9)", advance="no") ldist(l)
end do ! l = 1, toplev
#ifdef MPI
! store the process level distribution
!
write(funit,"(10x)", advance="no")
do l = 1, nprocs
write(funit,"(1x,i9)", advance="no") cdist(l)
end do ! l = 1, nprocs
#endif /* MPI */
! append the new line character
!
write(funit,"('')")
end if ! number of blocks or leafs changed
end if ! master
#ifdef PROFILE
! stop accounting time for the mesh statistics
!
call stop_timer(ims)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine store_mesh_stats
!
!===============================================================================
!
! 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)
! import external procedures and variables
!
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
#ifdef DEBUG
use blocks , only : check_neighbors
#endif /* DEBUG */
use coordinates , only : minlev, maxlev
use domains , only : setup_domain
use helpers , only : print_section
use iso_fortran_env, only : error_unit
use mpitools , only : master, nproc, nprocs, npmax
use problems , only : setup_problem
use refinement , only : check_refinement_criterion
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
! local variables
!
logical :: refine
integer(kind=4) :: level, lev
integer(kind=4) :: np, nl
integer(kind=4), dimension(0:npmax) :: lb
! local parameters
!
character(len=*), parameter :: loc = 'MESH::generate_mesh()'
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the initial mesh generation
!
call start_timer(img)
#endif /* PROFILE */
! allocate the initial lowest level structure of blocks
!
call setup_domain(status)
! quit if the domain setup failed
!
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot setup domain!"
go to 100
end if
! allocate temporary data block to determine the initial block structure
!
call allocate_datablock(pdata, status)
! quit if the block allocation failed
!
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot allocate 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
!
! print the levels while generating them
!
call print_section(master, "Generating the initial mesh")
if (master) write(*,"(4x,a16,11x,'=')",advance="no") "generating level"
! iterate over all levels up to the maximum one
!
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
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot refine 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
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot deallocate 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(), nprocs) - 1
lb( : ) = get_nleafs() / nprocs
lb(0:nl) = lb(0:nl) + 1
! reset the processor and block numbers
!
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
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot append new 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(np)) then
! reset the leaf counter for the current process
!
nl = 0
! increase the process number
!
np = min(npmax, np + 1)
end if ! nl >= lb(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
write(error_unit,"('[',a,']: ',a)") trim(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
#ifdef PROFILE
! stop accounting time for the initial mesh generation
!
call stop_timer(img)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
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)
! import external procedures and variables
!
use blocks, only : build_leaf_list
#ifdef DEBUG
use blocks, only : check_neighbors
#endif /* DEBUG */
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the adaptive mesh refinement update
!
call start_timer(imu)
#endif /* PROFILE */
! check the refinement criterion of all data blocks at the current process
!
call check_data_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
! redistribute blocks equally among all processors
!
call redistribute_blocks()
#endif /* MPI */
#ifdef DEBUG
! check if neighbors are consistent after mesh refinement
!
call check_neighbors()
#endif /* DEBUG */
! error entry point
!
100 continue
#ifdef PROFILE
! stop accounting time for the adaptive mesh refinement update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_mesh
!
!===============================================================================
!
! subroutine REDISTRIBUTE_BLOCKS:
! ------------------------------
!
! Subroutine redistributes data blocks between processes.
!
!
!===============================================================================
!
subroutine redistribute_blocks()
#ifdef MPI
! import external procedures and variables
!
use blocks , only : block_meta, block_data, list_meta
use blocks , only : get_nleafs
use blocks , only : append_datablock, remove_datablock, link_blocks
use coordinates , only : nn => bcells
use equations , only : nv
use mpitools , only : nprocs, npmax, nproc
use mpitools , only : send_array, receive_array
#endif /* MPI */
! local variables are not implicit by default
!
implicit none
#ifdef MPI
! local variables
!
integer :: status
integer(kind=4) :: np, nl
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
! tag for the MPI data exchange
!
integer(kind=4) :: itag
! array for number of data block for autobalancing
!
integer(kind=4), dimension(0:npmax) :: lb
! local buffer for data block exchange
!
#if NDIMS == 3
real(kind=8) , dimension(3,nv,nn,nn,nn) :: rbuf
#else /* NDIMS == 3 */
real(kind=8) , dimension(3,nv,nn,nn, 1) :: rbuf
#endif /* NDIMS == 3 */
#endif /* MPI */
!-------------------------------------------------------------------------------
!
#ifdef MPI
#ifdef PROFILE
! start accounting time for the block redistribution
!
call start_timer(ima)
#endif /* PROFILE */
! calculate the new data block division between processes
!
nl = mod(get_nleafs(), nprocs) - 1
lb( : ) = get_nleafs() / nprocs
lb(0:nl) = lb(0:nl) + 1
! reset the processor and leaf numbers
!
np = 0
nl = 0
! set the pointer to the first block on the meta block list
!
pmeta => list_meta
! iterate over all meta blocks and reassign their process numbers
!
do while (associated(pmeta))
! consider only meta blocks which belong to active processes
!
if (pmeta%process < nprocs) then
! check if the block belongs to another process
!
if (pmeta%process /= np) then
! check if the block is the leaf
!
if (pmeta%leaf) then
! generate a tag for communication
!
itag = 16 * (np * nprocs + pmeta%process)
! sends the block to the right process
!
if (nproc == pmeta%process) then
! copy data to buffer
!
rbuf(1,:,:,:,:) = pmeta%data%q (:,:,:,:)
rbuf(2,:,:,:,:) = pmeta%data%u0(:,:,:,:)
rbuf(3,:,:,:,:) = pmeta%data%u1(:,:,:,:)
! send data
!
call send_array(np, itag, rbuf)
! remove data block from the current process
!
call remove_datablock(pmeta%data, status)
! send data block
!
end if ! nproc == pmeta%process
! receive the block from another process
!
if (nproc == np) then
! allocate a new data block and link it with the current meta block
!
call append_datablock(pdata, status)
call link_blocks(pmeta, pdata)
! receive the data
!
call receive_array(pmeta%process, itag, rbuf)
! coppy the buffer to data block
!
pmeta%data%q (:,:,:,:) = rbuf(1,:,:,:,:)
pmeta%data%u0(:,:,:,:) = rbuf(2,:,:,:,:)
pmeta%data%u1(:,:,:,:) = rbuf(3,:,:,:,:)
end if ! nproc == n
end if ! leaf
! set new processor number
!
pmeta%process = np
end if ! pmeta%process /= np
! increase the number of blocks on the current process; if it exceeds the
! allowed number reset the counter and increase the processor number
!
if (pmeta%leaf) then
! increase the number of leafs for the current process
!
nl = nl + 1
! 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(np)) then
! reset the leaf counter for the current process
!
nl = 0
! increase the process number
!
np = min(npmax, np + 1)
end if ! l >= lb(n)
end if ! leaf
end if ! pmeta%process < nprocs
! assign the pointer to the next meta block
!
pmeta => pmeta%next
end do ! pmeta
#ifdef PROFILE
! stop accounting time for the block redistribution
!
call stop_timer(ima)
#endif /* PROFILE */
#endif /* MPI */
!-------------------------------------------------------------------------------
!
end subroutine redistribute_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)
! import external procedures and variables
!
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 interpolations , only : limiter_prol
use iso_fortran_env, only : error_unit
! local variables are not implicit by default
!
implicit none
! input arguments
!
type(block_meta), pointer, intent(inout) :: pmeta
integer , intent(out) :: status
! local variables
!
integer :: p
integer :: i, il, iu, ic, ip
integer :: j, jl, ju, jc, jp
integer :: k, kc
#if NDIMS == 3
integer :: kl, ku, kp
#endif /* NDIMS == 3 */
real(kind=8) :: dul, dur, du1, du2
#if NDIMS == 3
real(kind=8) :: du3, du4
#endif /* NDIMS == 3 */
! local pointers
!
type(block_meta), pointer :: pchild
type(block_data), pointer :: pdata
! local arrays
!
real(kind=8), dimension(NDIMS) :: du
! local parameters
!
character(len=*), parameter :: loc = 'MESH::prolong_block()'
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the block prolongation
!
call start_timer(imp)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! assign the pdata pointer
!
pdata => pmeta%data
! prepare indices
!
il = nb - nh
iu = ne + nh
jl = nb - nh
ju = ne + nh
#if NDIMS == 3
kl = nb - nh
ku = ne + nh
#endif /* NDIMS == 3 */
! expand the block variables
!
#if NDIMS == 2
k = 1
kc = 1
#endif /* NDIMS == 2 */
#if NDIMS == 3
do k = kl, ku
kc = 2 * (k - kl) + 1
kp = kc + 1
#endif /* NDIMS == 3 */
do j = jl, ju
jc = 2 * (j - jl) + 1
jp = jc + 1
do i = il, iu
ic = 2 * (i - il) + 1
ip = ic + 1
do p = 1, nv
du(1) = 1.25d-01 * (pdata%u(p,i+1,j,k) - pdata%u(p,i-1,j,k))
du(2) = 1.25d-01 * (pdata%u(p,i,j+1,k) - pdata%u(p,i,j-1,k))
#if NDIMS == 3
du(3) = 1.25d-01 * (pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k-1))
#endif /* NDIMS == 3 */
if (positive(p) .and. pdata%u(p,i,j,k) < sum(abs(du(1:NDIMS)))) then
if (pdata%u(p,i,j,k) > 0.0d+00) then
dul = pdata%u(p,i ,j,k) - pdata%u(p,i-1,j,k)
dur = pdata%u(p,i+1,j,k) - pdata%u(p,i ,j,k)
du(1) = limiter_prol(2.5d-01, dul, dur)
dul = pdata%u(p,i,j ,k) - pdata%u(p,i,j-1,k)
dur = pdata%u(p,i,j+1,k) - pdata%u(p,i,j ,k)
du(2) = limiter_prol(2.5d-01, dul, dur)
#if NDIMS == 3
dul = pdata%u(p,i,j,k ) - pdata%u(p,i,j,k-1)
dur = pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k )
du(3) = limiter_prol(2.5d-01, dul, dur)
#endif /* NDIMS == 3 */
do while (pdata%u(p,i,j,k) <= sum(abs(du(1:NDIMS))))
du(:) = 0.5d+00 * du(:)
end do
else
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
, "Positive variable is not positive at (", i, j, k, " )"
du(:) = 0.0d+00
end if
end if
#if NDIMS == 2
du1 = du(1) + du(2)
du2 = du(1) - du(2)
up(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1
up(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2
up(p,ic,jp,kc) = pdata%u(p,i,j,k) - du2
up(p,ip,jp,kc) = pdata%u(p,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)
up(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1
up(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2
up(p,ic,jp,kc) = pdata%u(p,i,j,k) - du3
up(p,ip,jp,kc) = pdata%u(p,i,j,k) + du4
up(p,ic,jc,kp) = pdata%u(p,i,j,k) - du4
up(p,ip,jc,kp) = pdata%u(p,i,j,k) + du3
up(p,ic,jp,kp) = pdata%u(p,i,j,k) - du2
up(p,ip,jp,kp) = pdata%u(p,i,j,k) + du1
#endif /* NDIMS == 3 */
end do
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
! iterate over all children
!
do p = 1, nchildren
! assign pointer to the current child
!
pchild => pmeta%child(p)%ptr
! obtain the position of child in the parent block
!
ic = pchild%pos(1)
jc = pchild%pos(2)
#if NDIMS == 3
kc = pchild%pos(3)
#endif /* NDIMS == 3 */
! calculate indices of the current child subdomain
!
il = 1 + ic * ni
jl = 1 + jc * ni
#if NDIMS == 3
kl = 1 + kc * ni
#endif /* NDIMS == 3 */
iu = il + nn - 1
ju = jl + nn - 1
#if NDIMS == 3
ku = kl + nn - 1
#endif /* NDIMS == 3 */
! copy data to the current child
!
#if NDIMS == 2
pchild%data%u(1:nv,:,:,:) = up(1:nv,il:iu,jl:ju, : )
#endif /* NDIMS == 2 */
#if NDIMS == 3
pchild%data%u(1:nv,:,:,:) = up(1:nv,il:iu,jl:ju,kl:ku)
#endif /* NDIMS == 3 */
end do ! nchildren
#ifdef PROFILE
! stop accounting time for the block prolongation
!
call stop_timer(imp)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
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
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the block restriction
!
call start_timer(imr)
#endif /* PROFILE */
! 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 */
#ifdef PROFILE
! stop accounting time for the block restriction
!
call stop_timer(imr)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine restrict_block
!
!===============================================================================
!!
!!*** PRIVATE SUBROUTINES ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine CHECK_DATA_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_data_block_refinement(status)
! import external procedures and variables
!
use blocks , only : block_meta, block_data, list_meta, list_data
#ifdef MPI
use blocks , only : get_nleafs
#endif /* MPI */
use coordinates , only : minlev, maxlev
use iso_fortran_env, only : error_unit
#ifdef MPI
#ifdef DEBUG
use mpitools , only : nproc
#endif /* DEBUG */
use mpitools , only : reduce_sum
#endif /* MPI */
use refinement , only : check_refinement_criterion
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
#ifdef MPI
! local variables
!
integer(kind=4) :: nl, l
! array for storing the refinement flags
!
integer(kind=4), dimension(:), allocatable :: ibuf
#endif /* MPI */
! local parameters
!
character(len=*), parameter :: loc = 'MESH::check_data_block_refinement()'
!-------------------------------------------------------------------------------
!
! reset the status flag
!
status = 0
! 1) reset the refinement flag for all meta blocks
!
! iterate over 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
!
! iterate over all data blocks
!
pdata => list_data
do while (associated(pdata))
! assign pmeta to the meta block associated with pdata
!
pmeta => pdata%meta
#ifdef DEBUG
! check if pmeta is associated
!
if (associated(pmeta)) then
! check if pmeta is a leaf
!
if (pmeta%leaf) then
#endif /* DEBUG */
! 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
#ifdef DEBUG
else ! pmeta is a leaf
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Associated meta block is not a leaf!"
end if ! pmeta is a leaf
else ! pmeta associated
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "No meta block associated with the current data block!"
end if ! pmeta associated
#endif /* DEBUG */
pdata => pdata%next
end do ! iterate over data blocks
#ifdef MPI
! 3) synchronize the refinement flags between all processes
!
! get the number of leafs
!
nl = get_nleafs()
! allocate a buffer for the refinement flags
!
allocate(ibuf(nl), stat = status)
! check if the buffer was allocated successfully
!
if (status == 0) then
! reset the buffer
!
ibuf(1:nl) = 0
! reset the leaf block counter
!
l = 0
! iterate over all meta blocks
!
pmeta => list_meta
do while (associated(pmeta))
! process only leafs
!
if (pmeta%leaf) then
! increase the leaf block counter
!
l = l + 1
! store the refinement flag in the buffer
!
ibuf(l) = pmeta%refine
end if ! pmeta is a leaf
pmeta => pmeta%next
end do ! iterate over meta blocks
! update refinement flags across all processes
!
call reduce_sum(ibuf(1:nl))
! reset the leaf block counter
!
l = 0
! iterate over all meta blocks
!
pmeta => list_meta
do while (associated(pmeta))
! process only leafs
!
if (pmeta%leaf) then
! increase the leaf block counter
!
l = l + 1
#ifdef DEBUG
! check if the MPI update process does not change the local refinement flags
!
if (pmeta%process == nproc .and. pmeta%refine /= ibuf(l)) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Refinement flag does not match after MPI reduction!"
end if
#endif /* DEBUG */
! restore the refinement flags
!
pmeta%refine = ibuf(l)
end if ! pmeta is a leaf
pmeta => pmeta%next
end do ! iterate over meta blocks
! deallocate the refinement flag buffer
!
deallocate(ibuf, stat = status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Refinement flag buffer could not be deallocated!"
end if
else ! buffer couldn't be allocated
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Refinement flag buffer could not be allocated!"
end if ! buffer couldn't be allocated
#endif /* MPI */
!-------------------------------------------------------------------------------
!
end subroutine check_data_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)
! import external procedures and variables
!
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
#ifdef MPI
use coordinates , only : nn => bcells
use equations , only : nv
#endif /* MPI */
use iso_fortran_env, only : error_unit
#ifdef MPI
use mpitools , only : nprocs, nproc
use mpitools , only : send_array, receive_array
#endif /* MPI */
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
! local pointers
!
type(block_meta), pointer :: pmeta, pparent, pchild
#ifdef MPI
type(block_data), pointer :: pdata
#endif /* MPI */
! local variables
!
logical :: flag
integer(kind=4) :: l, p
#ifdef MPI
! tag for the MPI data exchange
!
integer(kind=4) :: itag
! local buffer for data block exchange
!
#if NDIMS == 3
real(kind=8), dimension(nv,nn,nn,nn) :: rbuf
#else /* NDIMS == 3 */
real(kind=8), dimension(nv,nn,nn, 1) :: rbuf
#endif /* NDIMS == 3 */
#endif /* MPI */
! local parameters
!
character(len=*), parameter :: loc = 'MESH::check_children_derefinement()'
!-------------------------------------------------------------------------------
!
! reset the status flag
!
status = 0
! 1) check if the siblings of the block selected for derefinement, can be
! derefined as well, if not cancel the derefinement of all siblings
!
! iterate over levels and check sibling derefinement
!
do l = 2, toplev
! iterate over all meta blocks
!
pmeta => list_meta
do while (associated(pmeta))
! check only leafs at the current level
!
if (pmeta%leaf .and. pmeta%level == l) then
! check if block is selected for derefinement
!
if (pmeta%refine == -1) then
! assign pparent to the parent block of pmeta
!
pparent => pmeta%parent
#ifdef DEBUG
! check if pparent is associated
!
if (associated(pparent)) then
#endif /* DEBUG */
! reset derefinement flag
!
flag = .true.
! iterate over all children
!
do p = 1, nchildren
! assign pchild to the current child
!
pchild => pparent%child(p)%ptr
#ifdef DEBUG
! check if pchild is associated
!
if (associated(pchild)) then
#endif /* DEBUG */
! check if the current child is a leaf and selected for derefinement as well
!
flag = flag .and. (pchild%leaf .and. pchild%refine == -1)
#ifdef DEBUG
else ! pchild is associated
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Child does not exist!"
status = 1
go to 100
end if ! pparent is associated
#endif /* DEBUG */
end do ! over all children
! if children can be derefined, set the refine flag of the parent to -1,
! otherwise, cancel the derefinement of all siblings
!
if (flag) then
pparent%refine = -1
else
! iterate over all children
!
do p = 1, nchildren
! assign pchild to the current child
!
pchild => pparent%child(p)%ptr
! reset its derefinement
!
pchild%refine = max(0, pchild%refine)
end do ! children
end if ! ~flag
#ifdef DEBUG
else ! pparent is associated
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Current meta block has no parent!"
status = 1
go to 100
end if ! pparent is associated
#endif /* DEBUG */
end if ! %refine = -1
end if ! only leafs at level l
pmeta => pmeta%next
end do ! iterate over meta blocks
end do ! levels
#ifdef MPI
! 2) bring all siblings together to the same process
!
! assign pmeta to the first meta block on the list
!
pmeta => list_meta
! iterate over all meta blocks
!
do while (associated(pmeta))
! process only parent blocks (not leafs)
!
if (.not. pmeta%leaf) then
! check if the first child is selected for derefinement
!
if (pmeta%refine == -1) then
! assign pchild with the first child
!
pchild => pmeta%child(1)%ptr
! set the parent process to be the same as the first child
!
pmeta%process = pchild%process
! iterate over remaining children and if they are not on the same process,
! bring them to the parent's one
!
do p = 2, nchildren
! assign pchild to the current child
!
pchild => pmeta%child(p)%ptr
! if pchild belongs to a different process move its data block to the process
! of its parent
!
if (pchild%process /= pmeta%process) then
! generate the tag for communication
!
itag = pchild%process * nprocs + pmeta%process + nprocs + p + 1
! send data block from the current child to the parent process and deallocate it
!
if (pchild%process == nproc) then
! assign pdata to the daba block of the current child
!
pdata => pchild%data
#ifdef DEBUG
! check if pdata is associated
!
if (associated(pdata)) then
#endif /* DEBUG */
! copy data to the local buffer
!
rbuf(:,:,:,:) = pdata%u(:,:,:,:)
#ifdef DEBUG
else ! pdata associated
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Child has no data block associated!"
status = 1
go to 100
end if ! pdata associated
#endif /* DEBUG */
! send data
!
call send_array(pmeta%process, itag, rbuf)
! deallocate the associated data block (it has to be pchild%data, and not pdata,
! otherwise, pchild%data won't be nullified)
!
call remove_datablock(pchild%data, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot remove data block!"
go to 100
end if
end if ! pchild%process == nproc
! allocate data block at the curent child, and receive its data from
! a different process
!
if (pmeta%process == nproc) then
! receive the data
!
call receive_array(pchild%process, itag, rbuf(:,:,:,:))
! allocate data block for the current child
!
call append_datablock(pdata, status)
if (status == 0) then
call link_blocks(pchild, pdata)
! copy buffer to data block
!
pdata%u(:,:,:,:) = rbuf(:,:,:,:)
else
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot append data block!"
go to 100
end if
end if ! pmeta%process == nproc
! set the current processor of the block
!
pchild%process = pmeta%process
end if ! pchild belongs to a different process
end do ! children
end if ! pmeta children are selected for derefinement
end if ! the block is parent
pmeta => pmeta%next
end do ! iterate over meta blocks
#endif /* MPI */
! error entry point
!
100 continue
!-------------------------------------------------------------------------------
!
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)
! import external procedures and variables
!
use blocks , only : block_meta, block_data, list_meta
use blocks , only : append_datablock, link_blocks, derefine_block
use coordinates , only : toplev
use iso_fortran_env, only : error_unit
#ifdef MPI
use mpitools , only : nproc
#endif /* MPI */
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
! local variables
!
integer(kind=4) :: l
! local parameters
!
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
#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
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot append 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
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot derefine 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)
! import external procedures and variables
!
use blocks , only : block_meta, list_meta
use blocks , only : refine_block, remove_datablock
use coordinates , only : toplev
use iso_fortran_env, only : error_unit
#ifdef MPI
use mpitools , only : nproc
#endif /* MPI */
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
! local pointers
!
type(block_meta), pointer :: pmeta, pparent
! local variables
!
integer(kind=4) :: l
! local parameters
!
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
! 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
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot prolong 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
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot remove 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
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot refine current meta block!"
go to 100
end if
end if ! pmeta belongs to the current process
#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
!===============================================================================
!
end module