!!****************************************************************************** !! !! 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-2014 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 ! 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; ! verbose - flag determining if the subroutine should be verbose; ! iret - return flag of the procedure execution status; ! !=============================================================================== ! subroutine initialize_mesh(nrun, verbose, iret) ! import external procedures and variables ! use blocks , only : datablock_set_dims use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use coordinates , only : toplev, im, jm, km use equations , only : nv use mpitools , only : master, nprocs ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(in) :: nrun logical, intent(in) :: verbose integer, intent(inout) :: iret ! local variables ! character(len=64) :: fn integer(kind=4) :: i, j, k, l, n ! local arrays ! integer(kind=4), dimension(3) :: bm, rm, dm ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! set timer descriptions ! call set_timer('mesh initialization' , imi) call set_timer('mesh statistics' , ims) call set_timer('initial mesh generation', img) call set_timer('adaptive mesh update' , imu) call set_timer('block autobalancing' , ima) call set_timer('block restriction' , imr) call set_timer('block prolongation' , imp) ! start accounting time for module initialization/finalization ! call start_timer(imi) #endif /* PROFILE */ ! set data block dimensions ! call datablock_set_dims(nv, nv, im, jm, km) ! 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 open(unit = funit, file = fn, form = 'formatted', status = 'replace' & , buffered = 'yes') #else /* INTEL */ open(unit = funit, file = fn, form = 'formatted', status = 'replace') #endif /* INTEL */ ! 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 #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: ! ! iret - return flag of the procedure execution status; ! !=============================================================================== ! subroutine finalize_mesh(iret) ! import external procedures and variables ! use mpitools , only : master ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(inout) :: iret ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for module initialization/finalization ! call start_timer(imi) #endif /* PROFILE */ ! only master posses a handler to the mesh statistics file ! if (master) then ! close the statistics file ! close(funit) end if ! master #ifdef PROFILE ! stop accounting time for module initialization/finalization ! call stop_timer(imi) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine finalize_mesh ! !=============================================================================== ! ! subroutine STORE_MESH_STATS: ! --------------------------- ! ! Subroutine prepares and stores various mesh statistics. ! ! Arguments: ! ! step - the integration step; ! time - the physical time; ! !=============================================================================== ! subroutine store_mesh_stats(step, time) ! import external procedures and variables ! use blocks , only : ndims, block_meta, list_meta use blocks , only : get_mblocks, get_nleafs use coordinates , only : ng, im, jm, km, toplev, effres 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 real(kind=8) :: cv, ef ! local pointers ! type(block_meta), pointer :: pmeta ! local saved variables ! logical , save :: first = .true. integer(kind=4), save :: nm = 0, nl = 0, nt = 0, nc = 0 ! 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 */ ! store the mesh statistics only on master ! if (master) then ! store the mesh statistics only if they changed ! if (nm /= get_mblocks() .or. nl /= get_nleafs()) then ! get the mximum block number and maximum level ! if (first) then ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! determine the number of base blocks ! do while(associated(pmeta)) ! if the block is at the lowest level, count it ! if (pmeta%level == 1) nt = nt + 1 ! associate the pointer with the next block ! pmeta => pmeta%next end do ! pmeta ! calculate the maximum block number ! nt = nt * 2**(ndims*(toplev - 1)) ! calculate the number of cell in one block (including the ghost cells) ! nc = im * jm * km ! 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 ! block number) and the efficiency (the cells count for adaptive mesh ! divided by the cell count for corresponding uniform mesh) ! cv = (1.0d+00 * nl) / nt ef = (1.0d+00 * nl) * nc / product(effres(1:ndims) + 2 * ng) ! initialize the level and process block counter ! ldist(:) = 0 #ifdef MPI cdist(:) = 0 #endif /* MPI */ ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! scan all meta blocks and prepare get the block level and process ! distributions ! do while(associated(pmeta)) ! process only leafs ! if (pmeta%leaf) then ! increase the block level and process counts ! ldist(pmeta%level) = ldist(pmeta%level) + 1 #ifdef MPI cdist(pmeta%cpu+1) = cdist(pmeta%cpu+1) + 1 #endif /* MPI */ end if ! the leaf ! associate the pointer with the next block ! pmeta => pmeta%next end do ! pmeta ! write down the block statistics ! write(funit, "(i9,3e14.6,2(2x,i9))", advance="no") & step, time, cv, ef, nm, nl ! write down 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 ! write down 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 */ ! write the new line symbol ! 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. ! ! !=============================================================================== ! subroutine generate_mesh() ! import external procedures and variables ! use blocks , only : block_meta, block_data, list_meta, list_data use blocks , only : ndims, nchildren, nsides, nfaces use blocks , only : allocate_datablock, deallocate_datablock use blocks , only : append_datablock, remove_datablock use blocks , only : link_blocks, unlink_blocks, refine_block use blocks , only : get_mblocks, get_nleafs use coordinates , only : minlev, maxlev, res use domains , only : setup_domain use error , only : print_info, print_error use mpitools , only : master, nproc, nprocs use problems , only : setup_problem use refinement , only : check_refinement_criterion ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_meta), pointer :: pmeta, pneigh, pnext type(block_data), pointer :: pdata ! local variables ! integer(kind=4) :: level, lev, idir, iside, iface integer(kind=4) :: np, nl integer(kind=4), dimension(0:nprocs-1) :: lb !------------------------------------------------------------------------------- ! #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() ! 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 ! if (master) then ! print the levels while generating them ! write(*,*) write(*,"(1x,a)" ) "Generating the initial mesh:" write(*,"(4x,a)",advance="no") "generating level = " end if ! allocate temporary data block to determine the initial block structure ! call allocate_datablock(pdata) ! reset the currently processed level ! level = 1 ! iterate over all level up to the maximum one ! do while (level < maxlev) ! print the currently processed level ! if (master) write(*, "(1x,i2)", advance = "no") level !! DETERMINE THE REFINEMENT OF ALL BLOCKS AT CURRENT LEVEL !! ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks at the current level and initialize the problem ! for them ! do while (associated(pmeta)) ! set the initial conditions at the current block ! if (pmeta%level == level) then ! associate the temporary data block with the current meta block ! call link_blocks(pmeta, pdata) ! set up the initial conditions for the block ! call setup_problem(pdata) ! check the refinement criterion for the current block, and do not allow for ! the derefinement ! pmeta%refine = max(0, check_refinement_criterion(pdata)) ! if the level is lower then the minimum level set the block to be refined ! 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 ! unlink the data block from the current meta block ! call unlink_blocks(pmeta, pdata) end if ! pmeta%level == 1 ! assign pointer to the next block ! pmeta => pmeta%next end do ! pmeta !! 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 ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks at the level n and if the current block is ! selected for the refinement and its neighbors are at lower levels select them ! for refinement too; ! do while (associated(pmeta)) ! check if the current block is a leaf at the level n and selected for ! refinement ! if (pmeta%leaf .and. pmeta%level == lev .and. pmeta%refine == 1) then ! iterate over all neighbors ! do idir = 1, ndims do iside = 1, nsides do iface = 1, nfaces ! assign pointer to the neighbor ! pneigh => pmeta%neigh(idir,iside,iface)%ptr ! check if the neighbor is associated ! if (associated(pneigh)) then ! if the neighbor has lower level, select it to be refined too ! if (pneigh%level < pmeta%level) pneigh%refine = 1 end if ! neighbor's pointer is associated end do ! iface end do ! iside end do ! idir end if ! leaf at level n and marked for refinement ! assign pointer to the next block ! 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 ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks ! do while (associated(pmeta)) ! check if the current block is at the level lev and refine if if ! it is selected for refinement ! if (pmeta%level == lev .and. pmeta%refine == 1) then ! perform the refinement without creating new data blocks ! call refine_block(pmeta, res(pmeta%level + 1,:), .false.) end if ! selected for refinement and at the current level ! assign pointer to the next block ! pmeta => pmeta%next end do ! pmeta end do ! lev = 1, level ! increase the level number ! level = level + 1 end do ! level = 1, maxlev ! deallocate temporary data block ! call deallocate_datablock(pdata) ! print the currently processed level ! if (master) write(*, "(1x,i2)") level !! ASSOCIATE THE META BLOCKS WITH PROCESSES AND CREATE INITIAL DATA BLOCKS !! ! divide blocks between all processes, use the number of data blocks to do this, ! but keep blocks from the top level which have the same parent packed together ! 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 ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks ! do while (associated(pmeta)) ! assign the process number to the current block ! pmeta%cpu = 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%cpu == nproc) then ! append new data block ! call append_datablock(pdata) ! link the data block with the current meta block ! call link_blocks(pmeta, pdata) ! set up the initial conditions for the block ! call setup_problem(pdata) 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(nprocs - 1, np + 1) end if ! nl >= lb(np) end if ! leaf ! assign pointer to the next block ! pmeta => pmeta%next end do ! pmeta #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 data block, and ! refines or derefines it if necessary by restricting or prolongating ! the block data. In the MPI version the data blocks are redistributed ! among all processes after the mesh update. ! ! !=============================================================================== ! subroutine update_mesh() ! import external procedures and variables ! use blocks , only : block_meta, block_data, list_meta, list_data use blocks , only : nchildren, ndims, nsides, nfaces use blocks , only : get_nleafs use blocks , only : refine_block, derefine_block use blocks , only : append_datablock, remove_datablock, link_blocks use coordinates , only : minlev, maxlev, toplev, im, jm, km, res use equations , only : nv use error , only : print_info, print_error #ifdef MPI use mpitools , only : master, nprocs, nproc use mpitools , only : reduce_sum_integer_array use mpitools , only : send_real_array, receive_real_array #endif /* MPI */ use refinement , only : check_refinement_criterion ! local variables are not implicit by default ! implicit none ! local variables ! logical :: flag integer(kind=4) :: nl, i, j, k, l, n, p integer :: iret ! local pointers ! type(block_meta), pointer :: pmeta, pneigh, pchild, pparent type(block_data), pointer :: pdata #ifdef MPI ! tag for the MPI data exchange ! integer(kind=4) :: itag ! array for storing the refinement flags ! integer(kind=4), dimension(:), allocatable :: ibuf ! array for number of data block for autobalancing ! integer(kind=4), dimension(0:nprocs-1) :: lb ! local buffer for data block exchange ! real(kind=8) , dimension(nv,im,jm,km) :: rbuf #endif /* MPI */ !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the adaptive mesh refinement update ! call start_timer(imu) #endif /* PROFILE */ #ifdef DEBUG ! check the mesh when debugging ! call check_mesh('before update_mesh') #endif /* DEBUG */ !! DETERMINE THE REFINEMENT OF ALL DATA BLOCKS !! ! set the pointer to the first block on the data block list ! pdata => list_data ! iterate over all blocks in the data block list ! do while (associated(pdata)) ! assign a pointer to the meta block associated with the current data block ! pmeta => pdata%meta ! continue if the current data block has a meta block associated ! if (associated(pmeta)) then ! if the associated meta block is a leaf ! if (pmeta%leaf) then ! check the refinement criterion for the current data block ! pmeta%refine = check_refinement_criterion(pdata) ! correct the refinement of the block for the base and top 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 if ! pmeta is a leaf end if ! pmeta associated ! assign a pointer to the next data block ! pdata => pdata%next end do ! pdata #ifdef MPI !! EXCHANGE REFINEMENT FLAGS BETWEEN ALL PROCESSES !! ! get the number of leafs ! nl = get_nleafs() ! allocate a buffer for the refinement field values ! allocate(ibuf(nl)) ! reset the buffer ! ibuf(:) = 0 ! reset the leaf block counter ! l = 0 ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks ! do while (associated(pmeta)) ! process only leafs ! if (pmeta%leaf) then ! increase the leaf block counter ! l = l + 1 ! store the refinement flag for all blocks at the current process for ! exchange with other processors ! ibuf(l) = pmeta%refine end if ! pmeta is the leaf ! assign a pointer to the next meta block ! pmeta => pmeta%next end do ! pmeta ! update refinement flags across all processors ! call reduce_sum_integer_array(nl, ibuf(1:nl), iret) ! reset the leaf block counter ! l = 0 ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! iterate over all meta blocks ! do while (associated(pmeta)) ! process only leafs ! if (pmeta%leaf) then ! increase the leaf block counter ! l = l + 1 ! update non-local block refinement flags ! pmeta%refine = ibuf(l) end if ! pmeta is the leaf ! assign a pointer to the next meta block ! pmeta => pmeta%next end do ! pmeta ! deallocate the buffer ! if (allocated(ibuf)) deallocate(ibuf) #endif /* MPI */ !! SELECT NEIGHBORS OF REFINED BLOCKS TO BE REFINED IF NECESSARY !! ! iterate over all levels starting from top and correct the refinement ! of neighbor blocks ! do l = toplev, 1, -1 ! set the pointer to the first block on the meta block 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 ! iterte over all neighbors of the current leaf ! do i = 1, ndims do j = 1, nsides do k = 1, nfaces ! assign a pointer to the current neighbor ! pneigh => pmeta%neigh(i,j,k)%ptr ! check if the pointer is associated with any block ! if (associated(pneigh)) then !=== conditions for blocks selected to be refined ! if (pmeta%refine == 1) then ! if the neighbor is set to be derefined, reset its flag (this applies to ! blocks at the current or lower level) ! pneigh%refine = max(0, pneigh%refine) ! if the neighbor is at lower level, always set it to be refined ! if (pneigh%level < pmeta%level) pneigh%refine = 1 end if ! refine = 1 !=== conditions for blocks which stay at the same level ! if (pmeta%refine == 0) then ! if the neighbor lays at lower level and is set to be derefined, cancel its ! derefinement ! if (pneigh%level < pmeta%level) & pneigh%refine = max(0, pneigh%refine) end if ! refine = 0 !=== conditions for blocks which are selected to be derefined ! if (pmeta%refine == -1) then ! if the neighbor is at lower level and is set to be derefined, cancel its ! derefinement ! if (pneigh%level < pmeta%level) & pneigh%refine = max(0, pneigh%refine) ! if a neighbor is set to be refined, cancel the derefinement of current block ! if (pneigh%refine == 1) pmeta%refine = 0 end if ! refine = -1 end if ! associated(pneigh) end do ! nfaces end do ! nsides end do ! ndims end if ! the leaf at level l ! assign a pointer to the next block ! pmeta => pmeta%next end do ! meta blocks !! CHECK IF BLOCK CHILDREN CAN BE DEREFINED !! ! now check all derefined block if their siblings are set for derefinement too ! and are at the same level; check only levels > 1 ! if (l > 1) then ! set the pointer to the first block on the meta block list ! pmeta => list_meta ! iterate over all blocks ! do while (associated(pmeta)) ! check only leafs at the current level ! if (pmeta%leaf .and. pmeta%level == l) then ! check blocks which are selected to be derefined ! if (pmeta%refine == -1) then ! assign a pointer to the parent of the current block ! pparent => pmeta%parent ! check if parent is associated with any block ! if (associated(pparent)) then ! reset derefinement flag ! flag = .true. ! iterate over all children ! do p = 1, nchildren ! assign a pointer to the current child ! pchild => pparent%child(p)%ptr ! check if the current child is the leaf ! flag = flag .and. (pchild%leaf) ! check if the current child is set to be derefined ! flag = flag .and. (pchild%refine == -1) end do ! over all children ! if not all children can be derefined, cancel the derefinement ! if (.not. flag) then ! iterate over all children ! do p = 1, nchildren ! assign a pointer to the current child ! pchild => pparent%child(p)%ptr ! reset derefinement of the current child ! pchild%refine = max(0, pchild%refine) end do ! children end if ! ~flag end if ! pparent is associated end if ! refine = -1 end if ! leafs at level l ! assign a pointer to the next block ! pmeta => pmeta%next end do ! meta blocks end if ! l > 1 end do ! levels #ifdef MPI !! BRING BACK ALL CHILDREN SELECTED FOR DEREFINEMENT TO THE SAME PROCESS !! ! set the pointer to the first block on the meta block 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%child(1)%ptr%refine == -1) then ! check if the parent blocks is on the same processor as the next block, if not ! move it to the same processor ! if (pmeta%cpu /= pmeta%next%cpu) pmeta%cpu = pmeta%next%cpu ! find the case when child blocks are spread across at least 2 processors ! flag = .false. do p = 1, nchildren flag = flag .or. (pmeta%child(p)%ptr%cpu /= pmeta%cpu) end do if (flag) then ! iterate over all children ! do p = 1, nchildren ! generate the tag for communication ! itag = pmeta%child(p)%ptr%cpu * nprocs + pmeta%cpu & + nprocs + p + 1 ! if the current children is not on the same processor, then ... ! if (pmeta%child(p)%ptr%cpu /= pmeta%cpu) then ! if the meta block is on the same process ! if (pmeta%cpu == nproc) then ! allocate data blocks for children on the processor which will receive data ! call append_datablock(pdata) call link_blocks(pmeta%child(p)%ptr, pdata) ! receive the data ! call receive_real_array(size(rbuf), pmeta%child(p)%ptr%cpu & , itag, rbuf, iret) ! coppy buffer to data ! pmeta%child(p)%ptr%data%u(:,:,:,:) = rbuf(:,:,:,:) end if ! send data to the right processor and deallocate data block ! if (pmeta%child(p)%ptr%cpu == nproc) then ! copy data to buffer ! rbuf(:,:,:,:) = pmeta%child(p)%ptr%data%u(:,:,:,:) ! send data ! call send_real_array(size(rbuf), pmeta%cpu, itag, rbuf, iret) ! deallocate data block ! call remove_datablock(pmeta%child(p)%ptr%data) end if ! set the current processor of the block ! pmeta%child(p)%ptr%cpu = pmeta%cpu end if ! if child is are on different processes end do ! nchildren end if ! children spread over different processes end if ! children selected for derefinement end if ! the block is parent ! assign a pointer to the next block ! pmeta => pmeta%next end do ! over meta blocks #endif /* MPI */ !! DEREFINE SELECTED BLOCKS !! ! perform the actual derefinement ! do l = toplev, 2, -1 pmeta => list_meta do while (associated(pmeta)) if (pmeta%leaf) then if (pmeta%level .eq. l) then if (pmeta%refine .eq. -1) then pparent => pmeta%parent if (associated(pparent)) then #ifdef MPI if (pmeta%cpu .eq. nproc) then #endif /* MPI */ if (.not. associated(pparent%data)) then call append_datablock(pdata) call link_blocks(pparent, pdata) end if call restrict_block(pparent) #ifdef MPI end if #endif /* MPI */ call derefine_block(pparent) pmeta => pparent else call print_error("mesh::update_mesh" & , "Parent of derefined block is not associated!") end if end if end if end if pmeta => pmeta%next end do end do !! REFINE SELECTED BLOCKS !! ! perform the actual refinement starting from the lowest level ! do l = 1, toplev - 1 pmeta => list_meta do while (associated(pmeta)) if (pmeta%leaf) then if (pmeta%level .eq. l) then if (pmeta%refine .eq. 1) then pparent => pmeta #ifdef MPI if (pmeta%cpu .eq. nproc) then #endif /* MPI */ call refine_block(pmeta, res(pmeta%level + 1,:), .true.) call prolong_block(pparent) call remove_datablock(pparent%data) #ifdef MPI else call refine_block(pmeta, res(pmeta%level + 1,:), .false.) end if #endif /* MPI */ end if end if end if pmeta => pmeta%next end do end do ! l = 1, toplev - 1 #ifdef MPI ! redistribute blocks equally among all processors ! call redistribute_blocks() #endif /* MPI */ #ifdef DEBUG ! check mesh ! call check_mesh('after update_mesh') #endif /* DEBUG */ #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, list_data use blocks , only : get_nleafs use blocks , only : append_datablock, remove_datablock, link_blocks use coordinates , only : im, jm, km use equations , only : nv use mpitools , only : nprocs, nproc use mpitools , only : send_real_array, receive_real_array #endif /* MPI */ ! local variables are not implicit by default ! implicit none #ifdef MPI ! local variables ! integer :: iret 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:nprocs-1) :: lb ! local buffer for data block exchange ! real(kind=8) , dimension(nv,im,jm,km) :: rbuf #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)) ! check if the block belongs to another process ! if (pmeta%cpu /= np) then ! check if the block is the leaf ! if (pmeta%leaf) then ! generate a tag for communication ! itag = pmeta%cpu * nprocs + np + nprocs + 1 ! sends the block to the right process ! if (nproc == pmeta%cpu) then ! copy data to buffer ! rbuf(:,:,:,:) = pmeta%data%u(:,:,:,:) ! send data ! call send_real_array(size(rbuf), np, itag, rbuf, iret) ! remove data block from the current process ! call remove_datablock(pmeta%data) ! send data block ! end if ! nproc == pmeta%cpu ! 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) call link_blocks(pmeta, pdata) ! receive the data ! call receive_real_array(size(rbuf), pmeta%cpu, itag, rbuf, iret) ! coppy the buffer to data block ! pmeta%data%u(:,:,:,:) = rbuf(:,:,:,:) end if ! nproc == n end if ! leaf ! set new processor number ! pmeta%cpu = np end if ! pmeta%cpu /= 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(nprocs - 1, np + 1) end if ! l >= lb(n) end if ! leaf ! 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: ! ! pblock - the input meta block; ! !=============================================================================== ! subroutine prolong_block(pblock) ! import external procedures and variables ! use blocks , only : block_meta, block_data, nchildren use coordinates , only : ng, nh, in, jn, kn, im, jm, km use coordinates , only : ib, ie, jb, je, kb, ke use equations , only : nv use interpolations , only : limiter ! local variables are not implicit by default ! implicit none ! input arguments ! type(block_meta), pointer, intent(inout) :: pblock ! local variables ! integer :: i, j, k, q, p integer :: il, iu, jl, ju, kl, ku integer :: ic, jc, kc, ip, jp, kp real :: dul, dur, dux, duy, duz ! local pointers ! type(block_meta), pointer :: pchild type(block_data), pointer :: pdata ! local arrays ! integer, dimension(3) :: dm ! local allocatable arrays ! real, dimension(:,:,:,:), allocatable :: u !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the block prolongation ! call start_timer(imp) #endif /* PROFILE */ ! assign the pdata pointer ! pdata => pblock%data ! prepare dimensions ! dm(:) = 2 * ((/ in, jn, kn /) + ng) #if NDIMS == 2 dm(3) = 1 #endif /* NDIMS == 2 */ ! allocate array for the expanded arrays ! allocate(u(nv, dm(1), dm(2), dm(3))) ! prepare indices ! il = ib - nh iu = ie + nh jl = jb - nh ju = je + nh #if NDIMS == 3 kl = kb - nh ku = ke + nh #endif /* NDIMS == 3 */ ! expand the block variables ! #if NDIMS == 2 do k = 1, km 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) dux = limiter(0.25d+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) duy = limiter(0.25d+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 ) duz = limiter(0.25d+00, dul, dur) #endif /* NDIMS == 3 */ #if NDIMS == 2 u(p,ic,jc,kc) = pdata%u(p,i,j,k) - (dux + duy) u(p,ip,jc,kc) = pdata%u(p,i,j,k) + (dux - duy) u(p,ic,jp,kc) = pdata%u(p,i,j,k) + (duy - dux) u(p,ip,jp,kc) = pdata%u(p,i,j,k) + (dux + duy) #endif /* NDIMS == 2 */ #if NDIMS == 3 u(p,ic,jc,kc) = pdata%u(p,i,j,k) - dux - duy - duz u(p,ip,jc,kc) = pdata%u(p,i,j,k) + dux - duy - duz u(p,ic,jp,kc) = pdata%u(p,i,j,k) - dux + duy - duz u(p,ip,jp,kc) = pdata%u(p,i,j,k) + dux + duy - duz u(p,ic,jc,kp) = pdata%u(p,i,j,k) - dux - duy + duz u(p,ip,jc,kp) = pdata%u(p,i,j,k) + dux - duy + duz u(p,ic,jp,kp) = pdata%u(p,i,j,k) - dux + duy + duz u(p,ip,jp,kp) = pdata%u(p,i,j,k) + dux + duy + duz #endif /* NDIMS == 3 */ end do end do end do end do ! iterate over all children ! do p = 1, nchildren ! assign pointer to the current child ! pchild => pblock%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 * in jl = 1 + jc * jn #if NDIMS == 3 kl = 1 + kc * kn #endif /* NDIMS == 3 */ iu = il + im - 1 ju = jl + jm - 1 #if NDIMS == 3 ku = kl + km - 1 #endif /* NDIMS == 3 */ ! copy data to the current child ! #if NDIMS == 2 pchild%data%u(1:nv,1:im,1:jm,1:km) = u(1:nv,il:iu,jl:ju, 1:km) #endif /* NDIMS == 2 */ #if NDIMS == 3 pchild%data%u(1:nv,1:im,1:jm,1:km) = u(1:nv,il:iu,jl:ju,kl:ku) #endif /* NDIMS == 3 */ end do ! nchildren ! deallocate local arrays ! if (allocated(u)) deallocate(u) #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: ! ! pblock - the input meta block; ! !=============================================================================== ! subroutine restrict_block(pblock) ! import external procedures and variables ! use blocks , only : ndims use blocks , only : block_meta, block_data, nchildren use coordinates , only : ng, nh, in, jn, kn, im, jm, km use coordinates , only : ih, jh, kh, ib, jb, kb, ie, je, ke use equations , only : nv ! local variables are not implicit by default ! implicit none ! subroutine arguments ! type(block_meta), pointer, intent(inout) :: pblock ! local variables ! integer :: p integer :: il, jl, kl, iu, ju, ku, ip, jp, kp integer :: is, js, ks, it, jt, kt ! 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 => pblock%data ! iterate over all children ! do p = 1, nchildren ! assign a pointer to the current child ! pchild => pblock%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 = ie is = ib - nh it = ih else il = ib iu = im is = ih + 1 it = ie + nh end if ip = il + 1 if (pos(2) == 0) then jl = 1 ju = je js = jb - nh jt = jh else jl = jb ju = jm js = jh + 1 jt = je + nh end if jp = jl + 1 #if NDIMS == 3 if (pos(3) == 0) then kl = 1 ku = ke ks = kb - nh kt = kh else kl = kb ku = km ks = kh + 1 kt = ke + nh 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, 1 ) = & 2.50d-01 * ((pchild%u(1:nv,il:iu:2,jl:ju:2, 1 ) & + pchild%u(1:nv,ip:iu:2,jp:ju:2, 1 )) & + (pchild%u(1:nv,il:iu:2,jp:ju:2, 1 ) & + pchild%u(1:nv,ip:iu:2,jl:ju:2, 1 ))) #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 #ifdef PROFILE ! stop accounting time for the block restriction ! call stop_timer(imr) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine restrict_block #ifdef DEBUG ! !=============================================================================== ! ! check_mesh: subroutine checks if the block structure is correct ! subroutine CHECK_MESH: ! --------------------- ! ! Subroutine checks if the meta block structure is correct. ! ! Arguments: ! ! string - the identification string; ! !=============================================================================== ! subroutine check_mesh(string) ! import external procedures and variables ! use blocks , only : block_meta, list_meta use blocks , only : check_metablock ! local variables are not implicit by default ! implicit none ! input arguments ! character(len=*), intent(in) :: string ! local pointers ! type(block_meta), pointer :: pmeta !------------------------------------------------------------------------------- ! ! assign the pointer with the first block on the list ! pmeta => list_meta ! iterate over all meta blocks ! do while(associated(pmeta)) ! check the current block ! call check_metablock(pmeta, string) ! assign the pointer with the next block on the meta block list ! pmeta => pmeta%next end do ! over meta blocks !------------------------------------------------------------------------------- ! end subroutine check_mesh #endif /* DEBUG */ !=============================================================================== ! end module