!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2020 Grzegorz Kowal !! !! 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 . !! !!****************************************************************************** !! !! module: MESH - handling adaptive mesh structure !! !!****************************************************************************** ! module mesh #ifdef PROFILE ! import external subroutines ! use timers, only : set_timer, start_timer, stop_timer #endif /* PROFILE */ ! module variables are not implicit by default ! implicit none #ifdef PROFILE ! timer indices ! integer , save :: imi, ims, img, imu, ima, imp, imr #endif /* PROFILE */ ! file handler for the mesh statistics ! integer, save :: funit = 11 ! allocatable array for prolongation ! real(kind=8), dimension(:,:,:,:), allocatable :: up ! by default everything is private ! private ! declare public subroutines ! public :: initialize_mesh, finalize_mesh public :: generate_mesh, update_mesh public :: redistribute_blocks, store_mesh_stats !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_MESH: ! -------------------------- ! ! Subroutine initializes mesh module and its variables. ! ! Arguments: ! ! nrun - the restarted execution count; ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine initialize_mesh(nrun, status) ! import external procedures and variables ! use coordinates , only : ni => ncells, ng => nghosts, toplev use equations , only : nv use iso_fortran_env, only : error_unit use mpitools , only : master, nprocs ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(in) :: nrun integer, intent(out) :: status ! local variables ! character(len=64) :: fn integer :: l, n ! local arrays ! integer, dimension(3) :: pm = 1 ! local parameters ! character(len=*), parameter :: loc = 'MESH::initialize_mesh()' ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! set timer descriptions ! call set_timer('mesh:: initialization' , imi) call set_timer('mesh:: statistics' , ims) call set_timer('mesh:: initial generation', img) call set_timer('mesh:: adaptive update' , imu) call set_timer('mesh:: autobalancing' , ima) call set_timer('mesh:: block restriction' , imr) call set_timer('mesh:: block prolongation', imp) ! start accounting time for module initialization/finalization ! call start_timer(imi) #endif /* PROFILE */ ! reset the status flag ! status = 0 ! only master prepares the mesh statistics file ! if (master) then ! generate the mesh statistics file name ! write(fn, "('mesh_',i2.2,'.dat')") nrun ! create a new mesh statistics file ! #ifdef __INTEL_COMPILER open(unit = funit, file = fn, form = 'formatted', status = 'replace' & , buffered = 'yes') #else /* __INTEL_COMPILER */ open(unit = funit, file = fn, form = 'formatted', status = 'replace') #endif /* __INTEL_COMPILER */ ! write the mesh statistics header ! write(funit, "('#',2(4x,a4),10x,a8,6x,a10,5x,a6,6x,a5,4x,a20)") & 'step', 'time', 'coverage', 'efficiency' & , 'blocks', 'leafs', 'block distributions:' ! write the mesh distributions header ! write(funit, "('#',76x,a8)", advance="no") 'level = ' do l = 1, toplev write(funit, "(1x,i9)", advance="no") l end do #ifdef MPI write(funit, "(4x,a6)", advance="no") 'cpu = ' do n = 1, nprocs write(funit, "(1x,i9)", advance="no") n end do #endif /* MPI */ write(funit, "('' )") write(funit, "('#')") end if ! master ! prepare dimensions ! pm(1:NDIMS) = 2 * (ni + ng) ! allocate array for the prolongation array ! allocate(up(nv, pm(1), pm(2), pm(3)), stat = status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Prolongation array could not be allocated!" end if #ifdef PROFILE ! stop accounting time for module initialization/finalization ! call stop_timer(imi) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine initialize_mesh ! !=============================================================================== ! ! subroutine FINALIZE_MESH: ! ------------------------ ! ! Subroutine releases memory used by the module variables. ! ! Arguments: ! ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine finalize_mesh(status) ! import external procedures and variables ! use iso_fortran_env, only : error_unit use mpitools , only : master ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status ! local parameters ! character(len=*), parameter :: loc = 'MESH::finalize_mesh()' ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for module initialization/finalization ! call start_timer(imi) #endif /* PROFILE */ ! reset the status flag ! status = 0 ! only master posses a handler to the mesh statistics file ! if (master) close(funit) ! deallocate prolongation array ! if (allocated(up)) then deallocate(up, stat = status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Prolongation array could not be deallocated!" end if end if #ifdef PROFILE ! stop accounting time for module initialization/finalization ! call stop_timer(imi) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine finalize_mesh ! !=============================================================================== ! ! subroutine STORE_MESH_STATS: ! --------------------------- ! ! Subroutine calculates and stores various mesh statistics. ! ! Arguments: ! ! step - the integration iteration step number; ! time - the physical time; ! !=============================================================================== ! subroutine store_mesh_stats(step, time) ! import external procedures and variables ! use blocks , only : block_meta, block_leaf use blocks , only : list_meta, list_leaf use blocks , only : ndims use blocks , only : get_mblocks, get_nleafs use coordinates , only : ddims => domain_base_dims use coordinates , only : ni => ncells, nn => bcells use coordinates , only : nd => nghosts_double use coordinates , only : toplev use mpitools , only : master, nprocs ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer , intent(in) :: step real(kind=8), intent(in) :: time ! local variables ! integer(kind=4) :: l, n real(kind=8) :: cv, ef, ff ! local pointers ! type(block_meta), pointer :: pmeta type(block_leaf), pointer :: pleaf ! local saved variables ! logical , save :: first = .true. integer(kind=4), save :: nm = 0, nl = 0 real(kind=8) , save :: fcv = 1.0d+00, fef = 1.0d+00 ! local arrays ! integer(kind=4), dimension(toplev) :: ldist #ifdef MPI integer(kind=4), dimension(nprocs) :: cdist #endif /* MPI */ !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the mesh statistics ! call start_timer(ims) #endif /* PROFILE */ ! process and store the mesh statistics only on master node ! if (master) then ! and only if the number of blocks changed ! if (nm /= get_mblocks() .or. nl /= get_nleafs()) then ! get the maximum number of block and level, calculated only once ! if (first) then ! determine the number of base blocks ! n = 0 pmeta => list_meta do while(associated(pmeta)) if (pmeta%level == 1) n = n + 1 pmeta => pmeta%next end do ! pmeta ! calculate the maximum number of blocks ! n = n * 2**(ndims*(toplev - 1)) ! calculate the coverage and efficiency factors ! ff = 2**(toplev - 1) fcv = 1.0d+00 / n fef = 1.0d+00 * (ddims(1) * ni * ff + nd) / nn fef = fef * (ddims(2) * ni * ff + nd) / nn #if NDIMS == 3 fef = fef * (ddims(3) * ni * ff + nd) / nn #endif /* NDIMS == 3 */ ! reset the first execution flag ! first = .false. end if ! first ! get the new numbers of meta blocks and leafs ! nm = get_mblocks() nl = get_nleafs() ! calculate the coverage (the number of leafs divided by the maximum number ! of blocks) and the efficiency (the total number of cells in a corresponding ! uniform mesh divided by the number of cells for adaptive mesh case) ! cv = fcv * nl ef = fef / nl ! get the block level and block process distributions ! ldist(:) = 0 #ifdef MPI cdist(:) = 0 #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta ldist(pmeta%level) = ldist(pmeta%level) + 1 #ifdef MPI cdist(pmeta%process+1) = cdist(pmeta%process+1) + 1 #endif /* MPI */ pleaf => pleaf%next end do ! pleaf ! store the block statistics ! write(funit, "(i9,3es14.6e3,2(2x,i9))", advance="no") & step, time, cv, ef, nm, nl ! store the block level distribution ! write(funit,"(12x)", advance="no") do l = 1, toplev write(funit,"(1x,i9)", advance="no") ldist(l) end do ! l = 1, toplev #ifdef MPI ! store the process level distribution ! write(funit,"(10x)", advance="no") do l = 1, nprocs write(funit,"(1x,i9)", advance="no") cdist(l) end do ! l = 1, nprocs #endif /* MPI */ ! append the new line character ! write(funit,"('')") end if ! number of blocks or leafs changed end if ! master #ifdef PROFILE ! stop accounting time for the mesh statistics ! call stop_timer(ims) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine store_mesh_stats ! !=============================================================================== ! ! subroutine GENERATE_MESH: ! ------------------------ ! ! Subroutine generates the iinitial block structure by creating blocks ! according to the geometry, initial problems and refinement criterion. ! ! Arguments: ! ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine generate_mesh(status) ! import external procedures and variables ! use blocks , only : block_meta, block_data, list_meta use blocks , only : allocate_datablock, deallocate_datablock use blocks , only : append_datablock use blocks , only : link_blocks, unlink_blocks, refine_block use blocks , only : get_mblocks, get_nleafs use blocks , only : set_neighbors_refine use blocks , only : build_leaf_list #ifdef DEBUG use blocks , only : check_neighbors #endif /* DEBUG */ use coordinates , only : minlev, maxlev use domains , only : setup_domain use helpers , only : print_section use iso_fortran_env, only : error_unit use mpitools , only : master, nproc, nprocs, npmax use problems , only : setup_problem use refinement , only : check_refinement_criterion ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! local variables ! logical :: refine integer(kind=4) :: level, lev integer(kind=4) :: np, nl integer(kind=4), dimension(0:npmax) :: lb ! local parameters ! character(len=*), parameter :: loc = 'MESH::generate_mesh()' !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the initial mesh generation ! call start_timer(img) #endif /* PROFILE */ ! allocate the initial lowest level structure of blocks ! call setup_domain(status) ! quit if the domain setup failed ! if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot setup domain!" go to 100 end if ! allocate temporary data block to determine the initial block structure ! call allocate_datablock(pdata, status) ! quit if the block allocation failed ! if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot allocate temporary data block!" go to 100 end if ! at this point we assume, that the initial structure of blocks ! according to the defined geometry is already created; no refinement ! is done yet; we fill out the coarse blocks with the initial condition ! ! print the levels while generating them ! call print_section(master, "Generating the initial mesh") if (master) write(*,"(4x,a16,11x,'=')",advance="no") "generating level" ! iterate over all levels up to the maximum one ! refine = .true. level = 1 do while (level < maxlev .and. refine) !! DETERMINE THE REFINEMENT OF ALL BLOCKS AT THE CURRENT LEVEL !! ! iterate over all meta blocks at the current level, and initialize the problem ! for them using temporary data block ! refine = .false. pmeta => list_meta do while (associated(pmeta)) if (pmeta%level == level) then call link_blocks(pmeta, pdata) call setup_problem(pdata) ! check the refinement criterion for the current block, and do not allow for ! the derefinement; if the level is lower then the minimum level set it to be ! refined; ! pmeta%refine = max(0, check_refinement_criterion(pdata)) if (level < minlev) pmeta%refine = 1 ! if there is only one block, refine it anyway because the resolution for ! the problem initiation may be too small ! if (get_mblocks() == 1 .and. level == 1) pmeta%refine = 1 call unlink_blocks(pmeta, pdata) ! set the refine flag, if there is any block selected for refinement or ! derefinement ! if (pmeta%refine /= 0) refine = .true. end if ! pmeta%level == level pmeta => pmeta%next end do ! pmeta ! refine or derefine only if it is needed ! if (refine) then ! print the currently processed level ! if (master) write(*, "(1x,i0)", advance = "no") level !! STEP DOWN AND SELECT BLOCKS WHICH NEED TO BE REFINED !! ! walk through all levels down from the current level and check if neighbors ! of the refined blocks need to be refined as well; there is no need for ! checking the blocks at the lowest level; ! do lev = level, 2, -1 ! iterate over all meta blocks at the level lev and if the current block is ! selected for the refinement and its neighbors are at lower levels select them ! for refinement too; ! pmeta => list_meta do while (associated(pmeta)) if (pmeta%leaf .and. pmeta%level == lev .and. pmeta%refine == 1) & call set_neighbors_refine(pmeta) pmeta => pmeta%next end do ! pmeta end do ! lev = level, 2, -1 !! REFINE ALL BLOCKS FROM THE LOWEST LEVEL UP !! ! walk through the levels starting from the lowest to the current level ! do lev = 1, level pmeta => list_meta do while (associated(pmeta)) ! check if the current block is at the level lev and refine if if ! it is selected for refinement (refine without creating new data blocks) ! if (pmeta%level == lev .and. pmeta%refine == 1) then call refine_block(pmeta, .false., status) ! quit if the block couldn't be refined ! if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot refine meta block!" call deallocate_datablock(pdata, status) status = 1 go to 100 end if end if pmeta => pmeta%next end do ! pmeta end do ! lev = 1, level ! increase the level number ! level = level + 1 end if ! refine end do ! level = 1, maxlev ! print the currently processed level ! if (master) write(*, "(1x,i0)") level ! deallocate temporary data block ! call deallocate_datablock(pdata, status) ! check if the data block deallocation succeeded ! if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot deallocate temporary data block!" go to 100 end if !! ASSOCIATE THE META BLOCKS WITH PROCESSES AND CREATE INITIAL DATA BLOCKS !! ! divide leaf blocks between all processes, use the number of data blocks to do ! this, but keep at the same process blocks with the same parent ! nl = mod(get_nleafs(), nprocs) - 1 lb( : ) = get_nleafs() / nprocs lb(0:nl) = lb(0:nl) + 1 ! reset the processor and block numbers ! np = 0 nl = 0 ! iterate over all meta blocks ! pmeta => list_meta do while (associated(pmeta)) ! assign the process number to the current block ! pmeta%process = np ! check if the current block is the leaf ! if (pmeta%leaf) then ! increase the number of leafs for the current process ! nl = nl + 1 ! if the block belongs to the current process, append a new data block, link it ! to the current meta block and initialize the problem ! if (pmeta%process == nproc) then call append_datablock(pdata, status) if (status == 0) then call link_blocks(pmeta, pdata) call setup_problem(pdata) else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot append new data block!" go to 100 end if end if ! leaf and belongs to the current process ! if the number of leafs for the current process exceeds the number of assigned ! blocks, reset the counter and increase the process number ! if (nl >= lb(np)) then ! reset the leaf counter for the current process ! nl = 0 ! increase the process number ! np = min(npmax, np + 1) end if ! nl >= lb(np) end if ! leaf pmeta => pmeta%next end do ! pmeta ! update the list of leafs ! call build_leaf_list(status) ! quit if the data block allocation failed ! if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot build the list of leafs!" go to 100 end if #ifdef DEBUG ! check if neighbors are consistent after mesh generation ! call check_neighbors() #endif /* DEBUG */ ! error entry point ! 100 continue #ifdef PROFILE ! stop accounting time for the initial mesh generation ! call stop_timer(img) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine generate_mesh ! !=============================================================================== ! ! subroutine UPDATE_MESH: ! ---------------------- ! ! Subroutine checks the refinement criterion for each block, and refines ! or restricts it if necessary by prolongating or restricting its data. ! In the MPI version the data blocks are redistributed among all processes ! after the mesh update. ! ! Arguments: ! ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine update_mesh(status) ! import external procedures and variables ! use blocks, only : build_leaf_list #ifdef DEBUG use blocks, only : check_neighbors #endif /* DEBUG */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the adaptive mesh refinement update ! call start_timer(imu) #endif /* PROFILE */ ! check the refinement criterion of all data blocks at the current process ! call check_data_block_refinement(status) if (status /= 0) go to 100 ! update neighbor refinement flags, if they need to be refined as well ! call update_neighbor_refinement() ! prepare siblings of blocks marked for restriction ! call prepare_sibling_derefinement(status) if (status /= 0) go to 100 ! restrict selected blocks ! call derefine_selected_blocks(status) if (status /= 0) go to 100 ! prolong selected blocks ! call refine_selected_blocks(status) if (status /= 0) go to 100 ! update the list of leafs ! call build_leaf_list(status) if (status /= 0) go to 100 #ifdef MPI ! redistribute blocks equally among all processors ! call redistribute_blocks() #endif /* MPI */ #ifdef DEBUG ! check if neighbors are consistent after mesh refinement ! call check_neighbors() #endif /* DEBUG */ ! error entry point ! 100 continue #ifdef PROFILE ! stop accounting time for the adaptive mesh refinement update ! call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine update_mesh ! !=============================================================================== ! ! subroutine REDISTRIBUTE_BLOCKS: ! ------------------------------ ! ! Subroutine redistributes data blocks between processes. ! ! !=============================================================================== ! subroutine redistribute_blocks() #ifdef MPI ! import external procedures and variables ! use blocks , only : block_meta, block_data, list_meta use blocks , only : get_nleafs use blocks , only : append_datablock, remove_datablock, link_blocks use coordinates , only : nn => bcells use equations , only : nv use mpitools , only : nprocs, npmax, nproc use mpitools , only : send_array, receive_array #endif /* MPI */ ! local variables are not implicit by default ! implicit none #ifdef MPI ! local variables ! integer :: status integer(kind=4) :: np, nl ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! tag for the MPI data exchange ! integer(kind=4) :: itag ! array for number of data block for autobalancing ! integer(kind=4), dimension(0:npmax) :: lb ! local buffer for data block exchange ! #if NDIMS == 3 real(kind=8) , dimension(3,nv,nn,nn,nn) :: rbuf #else /* NDIMS == 3 */ real(kind=8) , dimension(3,nv,nn,nn, 1) :: rbuf #endif /* NDIMS == 3 */ #endif /* MPI */ !------------------------------------------------------------------------------- ! #ifdef MPI #ifdef PROFILE ! start accounting time for the block redistribution ! call start_timer(ima) #endif /* PROFILE */ ! calculate the new data block division between processes ! nl = mod(get_nleafs(), nprocs) - 1 lb( : ) = get_nleafs() / nprocs lb(0:nl) = lb(0:nl) + 1 ! reset the processor and leaf numbers ! np = 0 nl = 0 ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks and reassign their process numbers ! do while (associated(pmeta)) ! consider only meta blocks which belong to active processes ! if (pmeta%process < nprocs) then ! check if the block belongs to another process ! if (pmeta%process /= np) then ! check if the block is the leaf ! if (pmeta%leaf) then ! generate a tag for communication ! itag = 16 * (np * nprocs + pmeta%process) ! sends the block to the right process ! if (nproc == pmeta%process) then ! copy data to buffer ! rbuf(1,:,:,:,:) = pmeta%data%q (:,:,:,:) rbuf(2,:,:,:,:) = pmeta%data%u0(:,:,:,:) rbuf(3,:,:,:,:) = pmeta%data%u1(:,:,:,:) ! send data ! call send_array(np, itag, rbuf) ! remove data block from the current process ! call remove_datablock(pmeta%data, status) ! send data block ! end if ! nproc == pmeta%process ! receive the block from another process ! if (nproc == np) then ! allocate a new data block and link it with the current meta block ! call append_datablock(pdata, status) call link_blocks(pmeta, pdata) ! receive the data ! call receive_array(pmeta%process, itag, rbuf) ! coppy the buffer to data block ! pmeta%data%q (:,:,:,:) = rbuf(1,:,:,:,:) pmeta%data%u0(:,:,:,:) = rbuf(2,:,:,:,:) pmeta%data%u1(:,:,:,:) = rbuf(3,:,:,:,:) end if ! nproc == n end if ! leaf ! set new processor number ! pmeta%process = np end if ! pmeta%process /= np ! increase the number of blocks on the current process; if it exceeds the ! allowed number reset the counter and increase the processor number ! if (pmeta%leaf) then ! increase the number of leafs for the current process ! nl = nl + 1 ! if the number of leafs for the current process exceeds the number of assigned ! blocks, reset the counter and increase the process number ! if (nl >= lb(np)) then ! reset the leaf counter for the current process ! nl = 0 ! increase the process number ! np = min(npmax, np + 1) end if ! l >= lb(n) end if ! leaf end if ! pmeta%process < nprocs ! assign the pointer to the next meta block ! pmeta => pmeta%next end do ! pmeta #ifdef PROFILE ! stop accounting time for the block redistribution ! call stop_timer(ima) #endif /* PROFILE */ #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine redistribute_blocks ! !=============================================================================== ! ! subroutine PROLONG_BLOCK: ! ------------------------ ! ! Subroutine prolongs variables from a data blocks linked to the input ! meta block and copy the resulting array of variables to ! the newly created children data block. The process of data restriction ! conserves stored variables. ! ! Arguments: ! ! pmeta - the input meta block; ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine prolong_block(pmeta, status) ! import external procedures and variables ! use blocks , only : block_meta, block_data, nchildren use coordinates , only : ni => ncells, nn => bcells use coordinates , only : nh => nghosts_half use coordinates , only : nb, ne use equations , only : nv, positive use interpolations , only : limiter_prol use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_meta), pointer, intent(inout) :: pmeta integer , intent(out) :: status ! local variables ! integer :: p integer :: i, il, iu, ic, ip integer :: j, jl, ju, jc, jp integer :: k, kc #if NDIMS == 3 integer :: kl, ku, kp #endif /* NDIMS == 3 */ real(kind=8) :: dul, dur, du1, du2 #if NDIMS == 3 real(kind=8) :: du3, du4 #endif /* NDIMS == 3 */ ! local pointers ! type(block_meta), pointer :: pchild type(block_data), pointer :: pdata ! local arrays ! real(kind=8), dimension(NDIMS) :: du ! local parameters ! character(len=*), parameter :: loc = 'MESH::prolong_block()' !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the block prolongation ! call start_timer(imp) #endif /* PROFILE */ ! reset the status flag ! status = 0 ! assign the pdata pointer ! pdata => pmeta%data ! prepare indices ! il = nb - nh iu = ne + nh jl = nb - nh ju = ne + nh #if NDIMS == 3 kl = nb - nh ku = ne + nh #endif /* NDIMS == 3 */ ! expand the block variables ! #if NDIMS == 2 k = 1 kc = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = kl, ku kc = 2 * (k - kl) + 1 kp = kc + 1 #endif /* NDIMS == 3 */ do j = jl, ju jc = 2 * (j - jl) + 1 jp = jc + 1 do i = il, iu ic = 2 * (i - il) + 1 ip = ic + 1 do p = 1, nv dul = pdata%u(p,i ,j,k) - pdata%u(p,i-1,j,k) dur = pdata%u(p,i+1,j,k) - pdata%u(p,i ,j,k) du(1) = limiter_prol(0.5d+00, dul, dur) dul = pdata%u(p,i,j ,k) - pdata%u(p,i,j-1,k) dur = pdata%u(p,i,j+1,k) - pdata%u(p,i,j ,k) du(2) = limiter_prol(0.5d+00, dul, dur) #if NDIMS == 3 dul = pdata%u(p,i,j,k ) - pdata%u(p,i,j,k-1) dur = pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k ) du(3) = limiter_prol(0.5d+00, dul, dur) #endif /* NDIMS == 3 */ if (positive(p)) then if (pdata%u(p,i,j,k) > 0.0d+00) then do while (pdata%u(p,i,j,k) <= sum(abs(du(1:NDIMS)))) du(:) = 0.5d+00 * du(:) end do else write(error_unit,"('[',a,']: ',a,i9,a,3i4,a)") trim(loc) & , "Positive variable is not positive for block ID:" & , pmeta%id, " at ( ", i, j, k, " )" du(:) = 0.0d+00 status = 1 end if end if du(:) = 0.5d+00 * du(:) #if NDIMS == 2 du1 = du(1) + du(2) du2 = du(1) - du(2) up(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1 up(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2 up(p,ic,jp,kc) = pdata%u(p,i,j,k) - du2 up(p,ip,jp,kc) = pdata%u(p,i,j,k) + du1 #endif /* NDIMS == 2 */ #if NDIMS == 3 du1 = du(1) + du(2) + du(3) du2 = du(1) - du(2) - du(3) du3 = du(1) - du(2) + du(3) du4 = du(1) + du(2) - du(3) up(p,ic,jc,kc) = pdata%u(p,i,j,k) - du1 up(p,ip,jc,kc) = pdata%u(p,i,j,k) + du2 up(p,ic,jp,kc) = pdata%u(p,i,j,k) - du3 up(p,ip,jp,kc) = pdata%u(p,i,j,k) + du4 up(p,ic,jc,kp) = pdata%u(p,i,j,k) - du4 up(p,ip,jc,kp) = pdata%u(p,i,j,k) + du3 up(p,ic,jp,kp) = pdata%u(p,i,j,k) - du2 up(p,ip,jp,kp) = pdata%u(p,i,j,k) + du1 #endif /* NDIMS == 3 */ end do end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ ! iterate over all children ! do p = 1, nchildren ! assign pointer to the current child ! pchild => pmeta%child(p)%ptr ! obtain the position of child in the parent block ! ic = pchild%pos(1) jc = pchild%pos(2) #if NDIMS == 3 kc = pchild%pos(3) #endif /* NDIMS == 3 */ ! calculate indices of the current child subdomain ! il = 1 + ic * ni jl = 1 + jc * ni #if NDIMS == 3 kl = 1 + kc * ni #endif /* NDIMS == 3 */ iu = il + nn - 1 ju = jl + nn - 1 #if NDIMS == 3 ku = kl + nn - 1 #endif /* NDIMS == 3 */ ! copy data to the current child ! #if NDIMS == 2 pchild%data%u(1:nv,:,:,:) = up(1:nv,il:iu,jl:ju, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pchild%data%u(1:nv,:,:,:) = up(1:nv,il:iu,jl:ju,kl:ku) #endif /* NDIMS == 3 */ end do ! nchildren #ifdef PROFILE ! stop accounting time for the block prolongation ! call stop_timer(imp) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine prolong_block ! !=============================================================================== ! ! subroutine RESTRICT_BLOCK: ! ------------------------- ! ! Subroutine restricts variables from children data blocks linked to ! the input meta block and copy the resulting array of variables to ! the associated data block. The process of data restriction conserves ! stored variables. ! ! Arguments: ! ! pmeta - the input meta block; ! !=============================================================================== ! subroutine restrict_block(pmeta) ! import external procedures and variables ! use blocks , only : block_meta, block_data, nchildren use coordinates , only : nn => bcells, nh => bcells_half use coordinates , only : ng => nghosts_half use coordinates , only : nb, ne use equations , only : nv ! local variables are not implicit by default ! implicit none ! subroutine arguments ! type(block_meta), pointer, intent(inout) :: pmeta ! local variables ! integer :: p integer :: il, iu, ip, is, it integer :: jl, ju, jp, js, jt #if NDIMS == 3 integer :: kl, ku, kp, ks, kt #endif /* NDIMS == 3 */ ! local arrays ! integer, dimension(NDIMS) :: pos ! local pointers ! type(block_data), pointer :: pparent, pchild !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the block restriction ! call start_timer(imr) #endif /* PROFILE */ ! assign the parent data pointer ! pparent => pmeta%data ! iterate over all children ! do p = 1, nchildren ! assign a pointer to the current child ! pchild => pmeta%child(p)%ptr%data ! obtain the child position in the parent block ! pos(1:NDIMS) = pchild%meta%pos(1:NDIMS) ! calculate the bound indices of the source nad destination arrays ! if (pos(1) == 0) then il = 1 iu = ne is = nb - ng it = nh else il = nb iu = nn is = nh + 1 it = ne + ng end if ip = il + 1 if (pos(2) == 0) then jl = 1 ju = ne js = nb - ng jt = nh else jl = nb ju = nn js = nh + 1 jt = ne + ng end if jp = jl + 1 #if NDIMS == 3 if (pos(3) == 0) then kl = 1 ku = ne ks = nb - ng kt = nh else kl = nb ku = nn ks = nh + 1 kt = ne + ng end if kp = kl + 1 #endif /* NDIMS == 3 */ ! restrict conserved variables from the current child and copy the resulting ! array to the proper location of the parent data block ! #if NDIMS == 2 pparent%u(1:nv,is:it,js:jt, : ) = & 2.50d-01 * ((pchild%u(1:nv,il:iu:2,jl:ju:2, : ) & + pchild%u(1:nv,ip:iu:2,jp:ju:2, : )) & + (pchild%u(1:nv,il:iu:2,jp:ju:2, : ) & + pchild%u(1:nv,ip:iu:2,jl:ju:2, : ))) #endif /* NDIMS == 2 */ #if NDIMS == 3 pparent%u(1:nv,is:it,js:jt,ks:kt) = & 1.25d-01 * (((pchild%u(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + pchild%u(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + (pchild%u(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + pchild%u(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + ((pchild%u(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + pchild%u(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + (pchild%u(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + pchild%u(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) #endif /* NDIMS == 3 */ end do ! p = 1, nchildren ! prepare indices for the border filling ! is = nb - ng it = ne + ng js = nb - ng jt = ne + ng #if NDIMS == 3 ks = nb - ng kt = ne + ng #endif /* NDIMS == 3 */ ! fill out the borders ! #if NDIMS == 2 do il = is - 1, 1, -1 pparent%u(1:nv,il,js:jt, : ) = pparent%u(1:nv,is,js:jt, : ) end do do iu = it + 1, nn pparent%u(1:nv,iu,js:jt, : ) = pparent%u(1:nv,it,js:jt, : ) end do do jl = js - 1, 1, -1 pparent%u(1:nv, : ,jl, : ) = pparent%u(1:nv, : ,js, : ) end do do ju = jt + 1, nn pparent%u(1:nv, : ,ju, : ) = pparent%u(1:nv, : ,jt, : ) end do #endif /* NDIMS == 2 */ #if NDIMS == 3 do il = is - 1, 1, -1 pparent%u(1:nv,il,js:jt,ks:kt) = pparent%u(1:nv,is,js:jt,ks:kt) end do do iu = it + 1, nn pparent%u(1:nv,iu,js:jt,ks:kt) = pparent%u(1:nv,it,js:jt,ks:kt) end do do jl = js - 1, 1, -1 pparent%u(1:nv, : ,jl,ks:kt) = pparent%u(1:nv, : ,js,ks:kt) end do do ju = jt + 1, nn pparent%u(1:nv, : ,ju,ks:kt) = pparent%u(1:nv, : ,jt,ks:kt) end do do kl = ks - 1, 1, -1 pparent%u(1:nv, : , : ,kl) = pparent%u(1:nv, : , : ,ks) end do do ku = kt + 1, nn pparent%u(1:nv, : , : ,ku) = pparent%u(1:nv, : , : ,kt) end do #endif /* NDIMS == 3 */ #ifdef PROFILE ! stop accounting time for the block restriction ! call stop_timer(imr) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine restrict_block ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine CHECK_DATA_BLOCK_REFINEMENT: ! -------------------------------------- ! ! Subroutine scans over all data blocks, gets and corrects their refinement ! flags. If the MPI is used, the refinement flags are syncronized among all ! processes. ! ! Arguments: ! ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine check_data_block_refinement(status) ! import external procedures and variables ! use blocks , only : block_meta, block_data, list_meta, list_data #ifdef MPI use blocks , only : get_nleafs #endif /* MPI */ use coordinates , only : minlev, maxlev use iso_fortran_env, only : error_unit #ifdef MPI #ifdef DEBUG use mpitools , only : nproc #endif /* DEBUG */ use mpitools , only : reduce_sum #endif /* MPI */ use refinement , only : check_refinement_criterion ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata #ifdef MPI ! local variables ! integer(kind=4) :: nl, l ! array for storing the refinement flags ! integer(kind=4), dimension(:), allocatable :: ibuf #endif /* MPI */ ! local parameters ! character(len=*), parameter :: loc = 'MESH::check_data_block_refinement()' !------------------------------------------------------------------------------- ! ! reset the status flag ! status = 0 ! 1) reset the refinement flag for all meta blocks ! ! iterate over all meta blocks ! pmeta => list_meta do while (associated(pmeta)) pmeta%refine = 0 pmeta => pmeta%next end do ! iterate over meta blocks ! 2) determine the refinement of data block from the current process ! ! iterate over all data blocks ! pdata => list_data do while (associated(pdata)) ! assign pmeta to the meta block associated with pdata ! pmeta => pdata%meta #ifdef DEBUG ! check if pmeta is associated ! if (associated(pmeta)) then ! check if pmeta is a leaf ! if (pmeta%leaf) then #endif /* DEBUG */ ! check the refinement criterion for the current data block ! pmeta%refine = check_refinement_criterion(pdata) ! correct the refinement flag for the minimum and maximum levels ! if (pmeta%level < minlev) pmeta%refine = 1 if (pmeta%level == minlev) pmeta%refine = max(0, pmeta%refine) if (pmeta%level == maxlev) pmeta%refine = min(0, pmeta%refine) if (pmeta%level > maxlev) pmeta%refine = -1 #ifdef DEBUG else ! pmeta is a leaf write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Associated meta block is not a leaf!" end if ! pmeta is a leaf else ! pmeta associated write(error_unit,"('[',a,']: ',a)") trim(loc) & , "No meta block associated with the current data block!" end if ! pmeta associated #endif /* DEBUG */ pdata => pdata%next end do ! iterate over data blocks #ifdef MPI ! 3) synchronize the refinement flags between all processes ! ! get the number of leafs ! nl = get_nleafs() ! allocate a buffer for the refinement flags ! allocate(ibuf(nl), stat = status) ! check if the buffer was allocated successfully ! if (status == 0) then ! reset the buffer ! ibuf(1:nl) = 0 ! reset the leaf block counter ! l = 0 ! iterate over all meta blocks ! pmeta => list_meta do while (associated(pmeta)) ! process only leafs ! if (pmeta%leaf) then ! increase the leaf block counter ! l = l + 1 ! store the refinement flag in the buffer ! ibuf(l) = pmeta%refine end if ! pmeta is a leaf pmeta => pmeta%next end do ! iterate over meta blocks ! update refinement flags across all processes ! call reduce_sum(ibuf(1:nl)) ! reset the leaf block counter ! l = 0 ! iterate over all meta blocks ! pmeta => list_meta do while (associated(pmeta)) ! process only leafs ! if (pmeta%leaf) then ! increase the leaf block counter ! l = l + 1 #ifdef DEBUG ! check if the MPI update process does not change the local refinement flags ! if (pmeta%process == nproc .and. pmeta%refine /= ibuf(l)) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Refinement flag does not match after MPI reduction!" end if #endif /* DEBUG */ ! restore the refinement flags ! pmeta%refine = ibuf(l) end if ! pmeta is a leaf pmeta => pmeta%next end do ! iterate over meta blocks ! deallocate the refinement flag buffer ! deallocate(ibuf, stat = status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Refinement flag buffer could not be deallocated!" end if else ! buffer couldn't be allocated write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Refinement flag buffer could not be allocated!" end if ! buffer couldn't be allocated #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine check_data_block_refinement ! !=============================================================================== ! ! subroutine UPDATE_NEIGHBOR_REFINEMENT: ! ------------------------------------- ! ! Subroutine scans over all neighbors of blocks selected for refinement or ! derefinement, and if necessary selects them to be refined as well, or ! cancels their derefinement. ! ! !=============================================================================== ! subroutine update_neighbor_refinement() ! import external procedures and variables ! use blocks , only : block_meta, list_meta use blocks , only : set_neighbors_refine use coordinates , only : toplev ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_meta), pointer :: pmeta ! local variables ! integer(kind=4) :: l !------------------------------------------------------------------------------- ! ! iterate down over all levels and correct the refinement of neighbor blocks ! do l = toplev, 1, -1 ! assign pmeta to the first meta block on the list ! pmeta => list_meta ! iterate over all meta blocks ! do while (associated(pmeta)) ! check only leafs at the current level ! if (pmeta%leaf .and. pmeta%level == l) then ! correct neighbor refinement flags ! call set_neighbors_refine(pmeta) end if ! the leaf at level l ! assign pmeta to the next meta block ! pmeta => pmeta%next end do ! iterate over meta blocks end do ! levels !------------------------------------------------------------------------------- ! end subroutine update_neighbor_refinement ! !=============================================================================== ! ! subroutine PREPARE_SIBLING_DEREFINEMENT: ! --------------------------------------- ! ! Subroutine scans over all blocks selected for derefinement and checks if ! their siblings can be derefined as well. If any of them cannot be ! derefined, the derefinement of all siblings is canceled. Then, if MPI is ! used, the subroutine brings back all siblings together to lay on ! the same process. ! ! Note: This subroutine sets %refine flag of the parent to -1 to let the next ! executed subroutine derefine_selected_blocks() which parent block ! has to be derefined. That subroutine resets %refine flag of the ! parent after performing full restriction. ! ! Arguments: ! ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine prepare_sibling_derefinement(status) ! import external procedures and variables ! use blocks , only : block_meta, list_meta #ifdef MPI use blocks , only : block_data #endif /* MPI */ use blocks , only : nchildren use blocks , only : set_neighbors_refine #ifdef MPI use blocks , only : append_datablock, remove_datablock, link_blocks #endif /* MPI */ use coordinates , only : toplev #ifdef MPI use coordinates , only : nn => bcells use equations , only : nv #endif /* MPI */ use iso_fortran_env, only : error_unit #ifdef MPI use mpitools , only : nprocs, nproc use mpitools , only : send_array, receive_array #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status ! local pointers ! type(block_meta), pointer :: pmeta, pparent, pchild #ifdef MPI type(block_data), pointer :: pdata #endif /* MPI */ ! local variables ! logical :: flag integer(kind=4) :: l, p #ifdef MPI ! tag for the MPI data exchange ! integer(kind=4) :: itag ! local buffer for data block exchange ! #if NDIMS == 3 real(kind=8), dimension(nv,nn,nn,nn) :: rbuf #else /* NDIMS == 3 */ real(kind=8), dimension(nv,nn,nn, 1) :: rbuf #endif /* NDIMS == 3 */ #endif /* MPI */ ! local parameters ! character(len=*), parameter :: loc = 'MESH::check_children_derefinement()' !------------------------------------------------------------------------------- ! ! reset the status flag ! status = 0 ! 1) check if the siblings of the block selected for derefinement, can be ! derefined as well, if not cancel the derefinement of all siblings ! ! iterate over levels and check sibling derefinement ! do l = 2, toplev ! iterate over all meta blocks ! pmeta => list_meta do while (associated(pmeta)) ! check only leafs at the current level ! if (pmeta%leaf .and. pmeta%level == l) then ! check if block is selected for derefinement ! if (pmeta%refine == -1) then ! assign pparent to the parent block of pmeta ! pparent => pmeta%parent #ifdef DEBUG ! check if pparent is associated ! if (associated(pparent)) then #endif /* DEBUG */ ! reset derefinement flag ! flag = .true. ! iterate over all children ! do p = 1, nchildren ! assign pchild to the current child ! pchild => pparent%child(p)%ptr #ifdef DEBUG ! check if pchild is associated ! if (associated(pchild)) then #endif /* DEBUG */ ! check if the current child is a leaf and selected for derefinement as well ! flag = flag .and. (pchild%leaf .and. pchild%refine == -1) #ifdef DEBUG else ! pchild is associated write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Child does not exist!" status = 1 go to 100 end if ! pparent is associated #endif /* DEBUG */ end do ! over all children ! if children can be derefined, set the refine flag of the parent to -1, ! otherwise, cancel the derefinement of all siblings ! if (flag) then pparent%refine = -1 else ! iterate over all children ! do p = 1, nchildren ! assign pchild to the current child ! pchild => pparent%child(p)%ptr ! reset its derefinement ! pchild%refine = max(0, pchild%refine) end do ! children end if ! ~flag #ifdef DEBUG else ! pparent is associated write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Current meta block has no parent!" status = 1 go to 100 end if ! pparent is associated #endif /* DEBUG */ end if ! %refine = -1 end if ! only leafs at level l pmeta => pmeta%next end do ! iterate over meta blocks end do ! levels #ifdef MPI ! 2) bring all siblings together to the same process ! ! assign pmeta to the first meta block on the list ! pmeta => list_meta ! iterate over all meta blocks ! do while (associated(pmeta)) ! process only parent blocks (not leafs) ! if (.not. pmeta%leaf) then ! check if the first child is selected for derefinement ! if (pmeta%refine == -1) then ! assign pchild with the first child ! pchild => pmeta%child(1)%ptr ! set the parent process to be the same as the first child ! pmeta%process = pchild%process ! iterate over remaining children and if they are not on the same process, ! bring them to the parent's one ! do p = 2, nchildren ! assign pchild to the current child ! pchild => pmeta%child(p)%ptr ! if pchild belongs to a different process move its data block to the process ! of its parent ! if (pchild%process /= pmeta%process) then ! generate the tag for communication ! itag = pchild%process * nprocs + pmeta%process + nprocs + p + 1 ! send data block from the current child to the parent process and deallocate it ! if (pchild%process == nproc) then ! assign pdata to the daba block of the current child ! pdata => pchild%data #ifdef DEBUG ! check if pdata is associated ! if (associated(pdata)) then #endif /* DEBUG */ ! copy data to the local buffer ! rbuf(:,:,:,:) = pdata%u(:,:,:,:) #ifdef DEBUG else ! pdata associated write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Child has no data block associated!" status = 1 go to 100 end if ! pdata associated #endif /* DEBUG */ ! send data ! call send_array(pmeta%process, itag, rbuf) ! deallocate the associated data block (it has to be pchild%data, and not pdata, ! otherwise, pchild%data won't be nullified) ! call remove_datablock(pchild%data, status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot remove data block!" go to 100 end if end if ! pchild%process == nproc ! allocate data block at the curent child, and receive its data from ! a different process ! if (pmeta%process == nproc) then ! receive the data ! call receive_array(pchild%process, itag, rbuf(:,:,:,:)) ! allocate data block for the current child ! call append_datablock(pdata, status) if (status == 0) then call link_blocks(pchild, pdata) ! copy buffer to data block ! pdata%u(:,:,:,:) = rbuf(:,:,:,:) else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot append data block!" go to 100 end if end if ! pmeta%process == nproc ! set the current processor of the block ! pchild%process = pmeta%process end if ! pchild belongs to a different process end do ! children end if ! pmeta children are selected for derefinement end if ! the block is parent pmeta => pmeta%next end do ! iterate over meta blocks #endif /* MPI */ ! error entry point ! 100 continue !------------------------------------------------------------------------------- ! end subroutine prepare_sibling_derefinement ! !=============================================================================== ! ! subroutine DEREFINE_SELECTED_BLOCKS: ! ----------------------------------- ! ! Subroutine scans over all blocks and restrict those selected. ! ! Note: This subroutine resets the flag %refine set in subroutine ! prepare_sibling_derefinement(). ! ! Arguments: ! ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine derefine_selected_blocks(status) ! import external procedures and variables ! use blocks , only : block_meta, block_data, list_meta use blocks , only : append_datablock, link_blocks, derefine_block use coordinates , only : toplev use iso_fortran_env, only : error_unit #ifdef MPI use mpitools , only : nproc #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! local variables ! integer(kind=4) :: l ! local parameters ! character(len=*), parameter :: loc = 'MESH::derefine_selected_blocks()' !------------------------------------------------------------------------------- ! ! reset the status flag ! status = 0 ! iterate over levels and restrict the blocks selected for restriction ! do l = toplev - 1, 1, -1 ! iterate over all meta blocks ! pmeta => list_meta do while (associated(pmeta)) ! process non-leafs at the current level selected for restriction ! if (.not. pmeta%leaf .and. pmeta%level == l & .and. pmeta%refine == -1) then #ifdef MPI ! check if pmeta belongs to the current process ! if (pmeta%process == nproc) then #endif /* MPI */ ! check if a data block is associated with pmeta, if not create one ! if (.not. associated(pmeta%data)) then ! append new data block and link it with the current pmeta ! call append_datablock(pdata, status) if (status == 0) then call link_blocks(pmeta, pdata) else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot append new data block!" go to 100 end if end if ! no data block associated ! perform the block restriction ! call restrict_block(pmeta) #ifdef MPI end if ! pmeta belongs to the current process #endif /* MPI */ ! perform the mesh derefinement, and reset the refinement flag of ! the current block ! call derefine_block(pmeta, status) if (status == 0) then pmeta%refine = 0 else write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot derefine current meta block!" go to 100 end if end if ! non-leaf at current level selected for derefinement pmeta => pmeta%next end do ! iterate over meta blocks end do ! levels ! error entry level ! 100 continue !------------------------------------------------------------------------------- ! end subroutine derefine_selected_blocks ! !=============================================================================== ! ! subroutine REFINE_SELECTED_BLOCKS: ! --------------------------------- ! ! Subroutine scans over all blocks and prolongates those selected. ! ! Arguments: ! ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine refine_selected_blocks(status) ! import external procedures and variables ! use blocks , only : block_meta, list_meta use blocks , only : refine_block, remove_datablock use coordinates , only : toplev use iso_fortran_env, only : error_unit #ifdef MPI use mpitools , only : nproc #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status ! local pointers ! type(block_meta), pointer :: pmeta, pparent ! local variables ! integer(kind=4) :: l ! local parameters ! character(len=*), parameter :: loc = 'MESH::refine_selected_blocks()' ! !------------------------------------------------------------------------------- ! ! reset the status flag ! status = 0 ! iterate over all levels and prolong those selected for prolongation ! do l = 1, toplev - 1 ! assign pmeta to the first meta block on the list ! pmeta => list_meta ! iterate over all meta blocks ! do while (associated(pmeta)) ! process leafs at the current level selected for prolongation ! if (pmeta%leaf .and. pmeta%level == l .and. pmeta%refine == 1) then ! assign pparent with the new parent block ! pparent => pmeta #ifdef MPI ! check if pmeta belongs to the current process ! if (pmeta%process == nproc) then #endif /* MPI */ ! prepare child blocks with allocating the data blocks ! call refine_block(pmeta, .true., status) ! perform the data prolongation ! call prolong_block(pparent, status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot prolong current data block!" go to 100 end if ! remove the data block associated with the new parent ! call remove_datablock(pparent%data, status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot remove current data block!" go to 100 end if #ifdef MPI else ! pmeta belongs to the current process ! prepare child blocks without actually allocating the data blocks ! call refine_block(pmeta, .false., status) if (status /= 0) then write(error_unit,"('[',a,']: ',a)") trim(loc) & , "Cannot refine current meta block!" go to 100 end if end if ! pmeta belongs to the current process #endif /* MPI */ end if ! leaf at current level selected for refinement ! assign pmeta to the next meta block ! pmeta => pmeta%next end do ! iterate over meta blocks end do ! levels ! error entry level ! 100 continue !------------------------------------------------------------------------------- ! end subroutine refine_selected_blocks !=============================================================================== ! end module