gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/mesh.F90
Grzegorz Kowal 3263a08ec0 BOUNDARIES, MESH: Implement PPM block prolongation. 
Implement Colella and Woodward's Piecewise Parabolic Method (PPM) for
block prolongation. This method is also used for the prolongation of
boundaries from blocks at lower levels to higher levels.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2024-10-07 10:35:04 -03:00

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