gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/mesh.F90
Grzegorz Kowal 024fd4b82f MESH: Add a common workspace to be used in various modules. 
This workspace removes the necessity of allocating the local subroutine
arrays. It is allocated once in initialize_mesh() and can be used in any
place of the code. By default, the size of this workspace is controlled
by the block dimensions, the number of variables, and the simulation
dimension, but it can be controlled by a parameter 'workspace_size'.

Additionally, there is a logical flag 'work_in_use' exported from this
module to permit the verification if the workspace is currently in use.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2021-11-12 08:07:59 -03:00

2181 lines
58 KiB
Fortran
Raw Blame History

!!******************************************************************************
!!
!! This file is part of the AMUN source code, a program to perform
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!! adaptive mesh.
!!
!! Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: 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 */
! the handler of the mesh statistics file
!
integer, save :: funit = 11
! the size of the common workspace
!
integer, save :: nwork = 0
! the flag indicating that the workspace is in use
!
logical, save :: work_in_use = .false.
! 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
! the common workspace to use for local arrays
!
real(kind=8), dimension(:), allocatable, target :: work
private
public :: initialize_mesh, finalize_mesh
public :: generate_mesh, update_mesh
public :: redistribute_blocks, store_mesh_stats
public :: work, work_in_use
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
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)
use coordinates , only : nn => bcells, ni => ncells, ng => nghosts
use coordinates , only : toplev
use equations , only : nf
use iso_fortran_env, only : error_unit
use mpitools , only : master, nprocs
use parameters , only : get_parameter
implicit none
integer, intent(in) :: nrun
integer, intent(out) :: status
character(len=64) :: fn
integer :: l, n
character(len=*), parameter :: loc = 'MESH::initialize_mesh()'
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
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)
call start_timer(imi)
#endif /* PROFILE */
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, "(2x,a8)", advance="no") 'nproc = '
do n = 0, nprocs - 1
write(funit, "(1x,i9)", advance="no") n
end do
#endif /* MPI */
write(funit, "('' )")
write(funit, "('#')")
end if ! master
! prepare the dimensions of the prolongation array
!
pm(1:NDIMS) = 2 * (ni + ng)
! get the requested size of the workspace
!
call get_parameter('workspace_size', nwork)
nwork = max(nwork, (2 * NDIMS + 1) * nf * nn**NDIMS)
! allocate the workspace
!
allocate(work(nwork), stat=status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot allocate the common workspace!"
end if
#ifdef PROFILE
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)
use mpitools, only : master
implicit none
integer, intent(out) :: status
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
call start_timer(imi)
#endif /* PROFILE */
status = 0
if (master) close(funit)
if (allocated(work)) deallocate(work)
#ifdef PROFILE
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 helpers , only : flush_and_sync
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, nm, nl
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.
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 (changed) 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,"('')")
call flush_and_sync(funit)
end if ! number of blocks or leafs changed
end if ! master
! reset the flag indicating the block structure has changed
!
changed = .false.
#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)
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, nodes, lprocs
use problems , only : setup_problem
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()'
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
call start_timer(img)
#endif /* PROFILE */
call setup_domain(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot setup domain!"
go to 100
end if
call allocate_datablock(pdata, status)
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
!
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
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(), 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
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(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
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
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
use blocks , only : block_meta, block_data, list_meta
use blocks , only : get_nleafs, nregs
use blocks , only : append_datablock, remove_datablock, link_blocks
use coordinates , only : nn => bcells
use equations , only : nv
use mpitools , only : nprocs, npmax, nproc, nodes, lprocs
use mpitools , only : send_array, receive_array
#endif /* MPI */
implicit none
#ifdef MPI
logical :: flag
integer :: status
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
#endif /* MPI */
!-------------------------------------------------------------------------------
!
#ifdef MPI
#ifdef PROFILE
call start_timer(ima)
#endif /* PROFILE */
! calculate the new data block division 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
! reset the processor and leaf numbers
!
nc = 1
nr = 1
np = 0
nl = 0
pmeta => list_meta
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
if (pmeta%leaf) then
! indicate that the block structure will change
!
changed = .true.
! generate a tag for communication
!
tag1 = np * nprocs + pmeta%process
tag2 = tag1 + 1
! send the data block arrays to the destination process, and
! remove the data block from the current process
!
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)
end if ! nproc == pmeta%process
! on the destination process, allocate new data block, associate it with
! the corresponding meta block, and fill it with the received arrays
!
if (nproc == np) then
call append_datablock(pdata, status)
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 ! 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(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
end if ! pmeta%process < nprocs
pmeta => pmeta%next
end do ! pmeta
#ifdef PROFILE
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)
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
implicit none
type(block_meta), pointer, intent(inout) :: pmeta
integer , intent(out) :: status
type(block_meta), pointer :: pchild
type(block_data), pointer :: pdata
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 */
integer , dimension(NDIMS) :: l, u
real(kind=8), dimension(NDIMS) :: du
real(kind=8), dimension(:,:,:), allocatable :: work
character(len=*), parameter :: loc = 'MESH::prolong_block()'
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
call start_timer(imp)
#endif /* PROFILE */
status = 0
allocate(work(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
#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(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)
work(l(1),l(2),k) = pdata%u(n,i,j,k) - du1
work(u(1),l(2),k) = pdata%u(n,i,j,k) + du2
work(l(1),u(2),k) = pdata%u(n,i,j,k) - du2
work(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)
work(l(1),l(2),l(3)) = pdata%u(n,i,j,k) - du1
work(u(1),l(2),l(3)) = pdata%u(n,i,j,k) + du2
work(l(1),u(2),l(3)) = pdata%u(n,i,j,k) - du3
work(u(1),u(2),l(3)) = pdata%u(n,i,j,k) + du4
work(l(1),l(2),u(3)) = pdata%u(n,i,j,k) - du4
work(u(1),l(2),u(3)) = pdata%u(n,i,j,k) + du3
work(l(1),u(2),u(3)) = pdata%u(n,i,j,k) - du2
work(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) = work(l(1):u(1),l(2):u(2),k)
#endif /* NDIMS == 2 */
#if NDIMS == 3
pchild%data%u(n,:,:,:) = work(l(1):u(1),l(2):u(2),l(3):u(3))
#endif /* NDIMS == 3 */
end do ! nchildren
end do ! n = 1, nv
if (allocated(work)) then
deallocate(work, 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
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)
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 */
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) :: itag
#endif /* MPI */
logical, dimension(nchildren) :: flag
character(len=*), parameter :: loc = 'MESH::prepare_sibling_derefinement()'
!-------------------------------------------------------------------------------
!
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
pmeta => list_meta
do while (associated(pmeta))
if (pmeta%level == l) then
if (pmeta%leaf .and. 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 ! 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 (any(flag)) then
pparent%refine = -1
else
do p = 1, nchildren
pchild => pparent%child(p)%ptr
pchild%refine = max(0, pchild%refine)
end do ! children
end if ! ~flag
end if ! leafs selected for derefinement
end if ! only block 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
!
pmeta => list_meta
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
pchild => pmeta%child(1)%ptr
pmeta%process = pchild%process
! iterate over the remaining children and if any of them is not on
! the same process, bring it to the parent's one
!
do p = 2, nchildren
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 the child's data block from its process to the parent's process,
! and then deallocate it
!
if (pchild%process == nproc) then
call send_array(pmeta%process, itag, pchild%data%uu(:,:,:,:,:))
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
! on the parent's process, allocate a newdata block, and associate the data
! received from a different process
!
if (pmeta%process == nproc) then
call append_datablock(pdata, status)
if (status == 0) then
call link_blocks(pchild, pdata)
else
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot append data block!"
go to 100
end if
call receive_array(pchild%process, itag, pdata%uu(:,:,:,:,:))
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
! 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
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
! 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
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
! redistribute blocks equally among all processors
!
call redistribute_blocks()
#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
Reference in New Issue View Git Blame Copy Permalink