!!****************************************************************************** !! !! 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-2012 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 implicit none ! maximum number of covering blocks ! integer , save :: tblocks = 1 ! log file for the mesh statistics ! character(len=32), save :: fname = "mesh.log" integer , save :: funit = 11 contains ! !=============================================================================== ! ! init_mesh: subroutine initializes mesh module and coordinate variables ! !=============================================================================== ! subroutine init_mesh(flag) use blocks , only : datablock_set_dims use config , only : toplev, in, jn, kn, im, jm, km, ncells, rdims, ng & , xmin, xmax, ymin, ymax, zmin, zmax use mpitools , only : master, nprocs use variables, only : nqt, nvr implicit none ! input arguments ! logical, intent(in) :: flag ! local variables ! integer(kind=4) :: i, j, k, l, n logical :: info ! local arrays ! integer(kind=4), dimension(3) :: bm, rm, dm ! !------------------------------------------------------------------------------- ! ! set data block dimensions ! call datablock_set_dims(nqt, nvr, im, jm, km) ! print general information about resolutions ! if (master) then ! prepare the file for logging mesh statistics; only the master process handles ! this part ! ! check if the mesh statistics file exists ! inquire(file = fname, exist = info) ! if flag is set to .true. or the original mesh statistics file does not exist, ! create a new one and write the header in it, otherwise open the original file ! and move the writing position to the end to allow for appending ! if (flag .or. .not. info) then ! create a new mesh statistics file ! open(unit=funit, file=fname, form='formatted', status='replace') ! write the mesh statistics file header ! write(funit,"('#')") write(funit,"('#',4x,'step',5x,'time',11x,'leaf',4x,'meta'," // & "6x,'coverage',8x,'AMR',11x,'block distribution')") write(funit,"('#',27x,'blocks',2x,'blocks',4x,'efficiency'," // & "6x,'efficiency')",advance='no') write(funit,"(1x,a12)",advance="no") 'level = 1' do l = 2, toplev write(funit,"(2x,i6)",advance="no") l end do #ifdef MPI write(funit,"(1x,a10)",advance="no") 'cpu = 1' do n = 2, nprocs write(funit,"(2x,i6)",advance="no") n end do #endif /* MPI */ write(funit,"('' )") write(funit,"('#')") else ! open the mesh statistics file and set the position at the end of file ! open(unit=funit, file=fname, form='formatted', position='append') ! write a marker that the job has been restarted from here ! write(funit,"('#',1x,a)") "job restarted from this point" end if end if ! master !------------------------------------------------------------------------------- ! end subroutine init_mesh ! !=============================================================================== ! ! generate_mesh: subroutine generate the initial block structure by creating ! blocks according to the geometry, initial problem and ! refinement criterium ! !=============================================================================== ! subroutine generate_mesh() use blocks , only : block_meta, block_data, list_meta, list_data use blocks , only : refine_block, deallocate_datablock use blocks , only : nchild, nsides, nfaces use blocks , only : get_mblocks, get_nleafs use config , only : minlev, maxlev, rdims use coordinates, only : res use error , only : print_info, print_error use mpitools, only : master, nproc, nprocs use problem , only : init_domain, init_problem, check_ref implicit none ! local pointers ! type(block_meta), pointer :: pmeta, pneigh, pnext type(block_data), pointer :: pdata ! local variables ! integer(kind=4) :: i, j, k, l, n #ifdef MPI integer(kind=4), dimension(0:nprocs-1) :: lb #endif /* MPI */ !------------------------------------------------------------------------------- ! ! allocate the initial structure of blocks according to the problem ! call init_domain(res(1,:)) ! 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 write(*,*) write(*,"(1x,a)" ) "Generating the initial mesh:" write(*,"(4x,a)",advance="no") "generating level = " end if l = 1 do while (l .le. maxlev) ! print the level currently processed ! if (master) & write(*,"(1x,i2)",advance="no") l ! iterate over all data blocks at the current level and initialize the problem ! pdata => list_data do while (associated(pdata)) ! set the initial conditions at the current block ! if (pdata%meta%level .le. l) & call init_problem(pdata) ! assign pointer to the next block ! pdata => pdata%next end do ! at the maximum level we only initialize the problem (without checking the ! refinement criterion) ! if (l .lt. maxlev) then ! iterate over all data blocks at the current level and check the refinement ! criterion; do not allow for derefinement ! pdata => list_data do while (associated(pdata)) if (pdata%meta%level .eq. l) then pdata%meta%refine = max(0, check_ref(pdata)) ! if there is only one block, and it is set not to be refined, refine it anyway ! because the resolution for the problem initiation may be too small ! if (get_mblocks() .eq. 1 .and. l .eq. 1) & pdata%meta%refine = 1 ! if the level is lower then the minimum level set the block to be refined ! anyway ! if (l .lt. minlev) pdata%meta%refine = 1 end if ! assign pointer to the next block ! pdata => pdata%next end do ! walk through all levels down from the current level and check if select all ! neighbors for the refinement if they are at lower level; there is no need for ! checking the blocks at the lowest level; ! do n = l, 2, -1 ! 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; ! pmeta => list_meta do while (associated(pmeta)) ! check if the current block is at the level n, is a leaf, and selected for ! refinement ! if (pmeta%level .eq. n) then if (pmeta%leaf) then if (pmeta%refine .eq. 1) then ! iterate over all neighbors ! do i = 1, NDIMS do j = 1, nsides do k = 1, nfaces ! assign pointer to the neighbor ! pneigh => pmeta%neigh(i,j,k)%ptr ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor is a leaf, if not something wrong is going on ! if (pneigh%leaf) then ! if the neighbor has lower level, select it to be refined too ! if (pneigh%level .lt. pmeta%level) & pneigh%refine = 1 else call print_error("mesh::init_mesh", "Neighbor is not a leaf!") end if end if end do end do end do end if end if end if ! assign pointer to the next block ! pmeta => pmeta%next end do end do !! refine all selected blocks starting from the lowest level !! ! walk through the levels starting from the lowest to the current level ! do n = 1, l ! iterate over all meta blocks ! pmeta => list_meta do while (associated(pmeta)) ! check if the current block is at the level n and, if it is selected for ! refinement and if so, perform the refinement on this block ! if (pmeta%level .eq. n .and. pmeta%refine .eq. 1) then ! perform the refinement ! call refine_block(pmeta, res(pmeta%level + 1,:), .true.) end if ! assign pointer to the next block ! pmeta => pmeta%next end do end do end if l = l + 1 end do ! deallocate data blocks of non leafs ! pmeta => list_meta do while (associated(pmeta)) if (.not. pmeta%leaf) & call deallocate_datablock(pmeta%data) ! assign pointer to the next block ! pmeta => pmeta%next end do #ifdef MPI ! 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 ! l = mod(get_nleafs(), nprocs) - 1 lb( : ) = get_nleafs() / nprocs lb(0:l) = lb(0:l) + 1 ! reset the processor and block numbers n = 0 l = 0 pmeta => list_meta do while (associated(pmeta)) ! assign the cpu to the current block ! pmeta%cpu = n ! 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 l = l + 1 if (l .ge. lb(n)) then n = min(nprocs - 1, n + 1) l = 0 end if end if ! assign pointer to the next block ! pmeta => pmeta%next end do ! remove all data blocks which do not belong to the current process ! pmeta => list_meta do while (associated(pmeta)) pnext => pmeta%next ! if the current block belongs to another process and its data field is ! associated, deallocate its data field ! if (pmeta%cpu .ne. nproc .and. associated(pmeta%data)) & call deallocate_datablock(pmeta%data) ! assign pointer to the next block ! pmeta => pnext end do #endif /* MPI */ ! go to a new line after generating levels ! if (master) then write(*,*) end if !------------------------------------------------------------------------------- ! end subroutine generate_mesh ! !=============================================================================== ! ! update_mesh: subroutine check the refinement criterion for each block, ! refines or derefines it if necessary, and restricts or ! prolongates all data to the newly created blocks ! !=============================================================================== ! subroutine update_mesh() use blocks , only : block_meta, block_data, list_meta, list_data & , nchild, ndims, nsides, nfaces & , refine_block, derefine_block, append_datablock & , associate_blocks, deallocate_datablock use blocks , only : get_nleafs use config , only : minlev, maxlev, toplev, im, jm, km use coordinates, only : res use error , only : print_info, print_error #ifdef MPI use mpitools , only : reduce_sum_integer_array use mpitools , only : send_real_array, receive_real_array use mpitools , only : master, nprocs, nproc #endif /* MPI */ use problem , only : check_ref use variables, only : nqt implicit none ! local variables ! logical :: flag integer(kind=4) :: i, j, k, l, n, p integer :: iret #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(nqt,im,jm,km) :: rbuf #endif /* MPI */ ! local pointers ! type(block_meta), pointer :: pmeta, pneigh, pchild, pparent type(block_data), pointer :: pdata !------------------------------------------------------------------------------- ! #ifdef DEBUG ! check mesh ! call check_mesh('before update_mesh') #endif /* DEBUG */ ! iterate over elements of the data block list ! pdata => list_data do while (associated(pdata)) ! assign a pointer to the meta block associated with the current data block ! pmeta => pdata%meta ! 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_ref(pdata) ! correct the refinement of the block for the base and top levels ! if (pmeta%level .lt. minlev) pmeta%refine = 1 if (pmeta%level .eq. minlev) pmeta%refine = max(0, pmeta%refine) if (pmeta%level .eq. maxlev) pmeta%refine = min(0, pmeta%refine) if (pmeta%level .gt. 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 #ifdef MPI ! allocate buffer for the refinement field values ! allocate(ibuf(get_nleafs())) ! reset the buffer ! ibuf(:) = 0 ! store refinement flags for all blocks for exchange between processors ! l = 1 pmeta => list_meta do while (associated(pmeta)) if (pmeta%leaf) then ibuf(l) = pmeta%refine l = l + 1 end if pmeta => pmeta%next end do ! update refinement flags across all processors ! call reduce_sum_integer_array(get_nleafs(), ibuf(1:get_nleafs()), iret) ! update non-local block refinement flags ! l = 1 pmeta => list_meta do while (associated(pmeta)) if (pmeta%leaf) then pmeta%refine = ibuf(l) l = l + 1 end if pmeta => pmeta%next end do ! deallocate the buffer ! if (allocated(ibuf)) deallocate(ibuf) #endif /* MPI */ ! iterate over all levels starting from top and correct the refinement of blocks ! do l = toplev, 1, -1 ! 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 .eq. 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 which are selected to be refined ! if (pmeta%refine .eq. 1) then ! if the neighbor is set to be derefined, reset its flags (this applies to ! blocks at the current and lower levels) ! pneigh%refine = max(0, pneigh%refine) ! if the neighbor is at lower level, always set it to be refined ! if (pneigh%level .lt. pmeta%level) pneigh%refine = 1 end if ! refine = 1 != conditions for blocks which stay at the same level ! if (pmeta%refine .eq. 0) then ! if the neighbor lays at lower level and is set to be derefined, cancel its ! derefinement ! if (pneigh%level .lt. pmeta%level) & pneigh%refine = max(0, pneigh%refine) end if ! refine = 0 != conditions for blocks which are selected to be derefined ! if (pmeta%refine .eq. -1) then ! if the neighbor is at lower level and is set to be derefined, cancel its ! derefinement ! if (pneigh%level .lt. pmeta%level) & pneigh%refine = max(0, pneigh%refine) ! if the neighbor is set to be refined, cancel derefinement of the current block ! if (pneigh%refine .eq. 1) pmeta%refine = 0 end if ! refine = -1 end if ! associated(pneigh) end do end do end do end if ! leafs at level l ! assign a pointer to the next block ! pmeta => pmeta%next end do ! meta blocks != now check all derefined block if their siblings are set for derefinement too ! and are at the same level; check ony levels >= 2 ! if (l .ge. 2) then ! iterate over all blocks ! pmeta => list_meta do while (associated(pmeta)) ! check only leafs at the current level ! if (pmeta%leaf .and. pmeta%level .eq. l) then ! check blocks which are selected to be derefined ! if (pmeta%refine .eq. -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, nchild ! assign a pointer to the current child ! pchild => pparent%child(p)%ptr ! check if the current child is a leaf ! flag = flag .and. (pchild%leaf) ! check if the current child is set to be derefined ! flag = flag .and. (pchild%refine .eq. -1) end do ! over all children ! if not all children are proper for derefinement, cancel derefinement of all ! children ! if (.not. flag) then ! iterate over all children ! do p = 1, nchild ! 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 >= 2 end do ! levels #ifdef MPI ! find all sibling blocks which are spread over different processors ! pmeta => list_meta do while (associated(pmeta)) if (.not. pmeta%leaf) then if (pmeta%child(1)%ptr%refine .eq. -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 .ne. 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, nchild flag = flag .or. (pmeta%child(p)%ptr%cpu .ne. pmeta%cpu) end do if (flag) then ! iterate over all children ! do p = 1, nchild ! 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 .ne. pmeta%cpu) then ! allocate data blocks for children on the processor which will receive data ! if (pmeta%cpu .eq. nproc) then call append_datablock(pdata) call associate_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 .eq. 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 deallocate_datablock(pmeta%child(p)%ptr%data) end if ! set the current processor of the block ! pmeta%child(p)%ptr%cpu = pmeta%cpu end if end do end if end if end if pmeta => pmeta%next end do #endif /* MPI */ ! 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 associate_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 ! perform the actual refinement ! 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 deallocate_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 #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 */ !------------------------------------------------------------------------------- ! end subroutine update_mesh #ifdef MPI ! !=============================================================================== ! ! redistribute_blocks: subroutine redistributes blocks equally between ! processors ! !=============================================================================== ! subroutine redistribute_blocks() use blocks , only : block_meta, block_data, list_meta, list_data use blocks , only : get_nleafs, append_datablock, deallocate_datablock & , associate_blocks use config , only : im, jm, km use mpitools , only : send_real_array, receive_real_array use mpitools , only : nprocs, nproc use variables, only : nqt implicit none ! local variables ! integer :: iret integer(kind=4) :: l, n ! ! 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(nqt,im,jm,km) :: rbuf ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata !------------------------------------------------------------------------------- ! !! AUTO BALANCING !! ! calculate the new division ! l = mod(get_nleafs(), nprocs) - 1 lb( : ) = get_nleafs() / nprocs lb(0:l) = lb(0:l) + 1 ! iterate over all metablocks and reassign the processor numbers ! n = 0 l = 0 pmeta => list_meta do while (associated(pmeta)) ! assign the cpu to the current block ! if (pmeta%cpu .ne. n) then if (pmeta%leaf) then ! generate the tag for communication ! itag = pmeta%cpu * nprocs + n + nprocs + 1 if (nproc .eq. pmeta%cpu) then ! copy data to buffer ! rbuf(:,:,:,:) = pmeta%data%u(:,:,:,:) ! send data ! call send_real_array(size(rbuf), n, itag, rbuf, iret) ! deallocate data block ! call deallocate_datablock(pmeta%data) ! send data block ! end if if (nproc .eq. n) then ! allocate data block for the current block ! call append_datablock(pdata) call associate_blocks(pmeta, pdata) ! receive the data ! call receive_real_array(size(rbuf), pmeta%cpu, itag, rbuf, iret) ! coppy buffer to data ! pmeta%data%u(:,:,:,:) = rbuf(:,:,:,:) ! receive data block ! end if end if ! set new processor number ! pmeta%cpu = n end if ! 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 l = l + 1 if (l .ge. lb(n)) then n = min(nprocs - 1, n + 1) l = 0 end if end if ! assign pointer to the next block ! pmeta => pmeta%next end do !------------------------------------------------------------------------------- ! end subroutine redistribute_blocks #endif /* MPI */ ! !=============================================================================== ! ! prolong_block: subroutine expands the block data and copy them to children ! !=============================================================================== ! subroutine prolong_block(pblock) use blocks , only : block_meta, block_data, nchild use config , only : ng, nh, in, jn, kn, im, jm, km use config , only : ib, ie, jb, je, kb, ke use interpolation, only : minmod use variables , only : nqt 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 arrays ! integer, dimension(3) :: dm ! local allocatable arrays ! real, dimension(:,:,:,:), allocatable :: u ! local pointers ! type(block_meta), pointer :: pchild type(block_data), pointer :: pdata !------------------------------------------------------------------------------- ! ! 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(nqt, 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, nqt dur = pdata%u(p,i+1,j,k) - pdata%u(p,i ,j,k) dul = pdata%u(p,i ,j,k) - pdata%u(p,i-1,j,k) dux = 0.125d0 * minmod(dur, dul) dur = pdata%u(p,i,j+1,k) - pdata%u(p,i,j ,k) dul = pdata%u(p,i,j ,k) - pdata%u(p,i,j-1,k) duy = 0.125d0 * minmod(dur, dul) #if NDIMS == 3 dur = pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k ) dul = pdata%u(p,i,j,k ) - pdata%u(p,i,j,k-1) duz = 0.125d0 * minmod(dur, dul) #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) - dux + duy 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, nchild ! 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:nqt,1:im,1:jm,1:km) = u(1:nqt,il:iu,jl:ju, 1:km) #endif /* NDIMS == 2 */ #if NDIMS == 3 pchild%data%u(1:nqt,1:im,1:jm,1:km) = u(1:nqt,il:iu,jl:ju,kl:ku) #endif /* NDIMS == 3 */ end do ! deallocate local arrays ! if (allocated(u)) deallocate(u) !------------------------------------------------------------------------------- ! end subroutine prolong_block ! !=============================================================================== ! ! restrict_block: subroutine shrinks the block data and copy them from children ! !=============================================================================== ! subroutine restrict_block(pblock) use blocks , only : block_meta, block_data, nchild use config , only : ng, in, ih, im, ib, ie, nh, jn, jh, jm, jb, je & , kn, kh, km, kb, ke use variables , only : nfl #ifdef MHD use variables , only : ibx, iby, ibz #ifdef GLM use variables , only : iph #endif /* GLM */ #endif /* MHD */ implicit none ! input arguments ! type(block_meta), pointer, intent(inout) :: pblock ! local variables ! integer :: p integer :: if, jf, kf integer :: il, jl, kl, iu, ju, ku integer :: ip, jp, kp integer :: is, js, ks, it, jt, kt ! local pointers ! type(block_data), pointer :: pparent, pchild !------------------------------------------------------------------------------- ! ! assign pointers ! pparent => pblock%data ! iterate over all children ! do p = 1, nchild ! assign pointer to the current child ! pchild => pblock%child(p)%ptr%data ! obtain the position of the current child in the parent block ! if = pchild%meta%pos(1) jf = pchild%meta%pos(2) #if NDIMS == 3 kf = pchild%meta%pos(3) #endif /* NDIMS == 3 */ ! calculate the bound indices of the source nad destination arrays ! if (if .eq. 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 (jf .eq. 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 (kf .eq. 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 */ ! copy the variables from the current child to the proper location of ! the parent block ! #if NDIMS == 2 pparent%u(1:nfl,is:it,js:jt,1) = & 0.25 * (pchild%u(1:nfl,il:iu:2,jl:ju:2,1) & + pchild%u(1:nfl,ip:iu:2,jl:ju:2,1) & + pchild%u(1:nfl,il:iu:2,jp:ju:2,1) & + pchild%u(1:nfl,ip:iu:2,jp:ju:2,1)) #ifdef MHD pparent%u(ibx:ibz,is:it,js:jt,1) = & 0.25 * (pchild%u(ibx:ibz,il:iu:2,jl:ju:2,1) & + pchild%u(ibx:ibz,ip:iu:2,jl:ju:2,1) & + pchild%u(ibx:ibz,il:iu:2,jp:ju:2,1) & + pchild%u(ibx:ibz,ip:iu:2,jp:ju:2,1)) #ifdef GLM pparent%u(iph ,is:it,js:jt,1) = & 0.25 * (pchild%u(iph ,il:iu:2,jl:ju:2,1) & + pchild%u(iph ,ip:iu:2,jl:ju:2,1) & + pchild%u(iph ,il:iu:2,jp:ju:2,1) & + pchild%u(iph ,ip:iu:2,jp:ju:2,1)) #endif /* GLM */ #endif /* MHD */ #endif /* NDIMS == 2 */ #if NDIMS == 3 pparent%u(1:nfl,is:it,js:jt,ks:kt) = & 0.125 * (pchild%u(1:nfl,il:iu:2,jl:ju:2,kl:ku:2) & + pchild%u(1:nfl,ip:iu:2,jl:ju:2,kl:ku:2) & + pchild%u(1:nfl,il:iu:2,jp:ju:2,kl:ku:2) & + pchild%u(1:nfl,ip:iu:2,jp:ju:2,kl:ku:2) & + pchild%u(1:nfl,il:iu:2,jl:ju:2,kp:ku:2) & + pchild%u(1:nfl,ip:iu:2,jl:ju:2,kp:ku:2) & + pchild%u(1:nfl,il:iu:2,jp:ju:2,kp:ku:2) & + pchild%u(1:nfl,ip:iu:2,jp:ju:2,kp:ku:2)) #ifdef MHD pparent%u(ibx:ibz,is:it,js:jt,ks:kt) = & 0.125 * (pchild%u(ibx:ibz,il:iu:2,jl:ju:2,kl:ku:2) & + pchild%u(ibx:ibz,ip:iu:2,jl:ju:2,kl:ku:2) & + pchild%u(ibx:ibz,il:iu:2,jp:ju:2,kl:ku:2) & + pchild%u(ibx:ibz,ip:iu:2,jp:ju:2,kl:ku:2) & + pchild%u(ibx:ibz,il:iu:2,jl:ju:2,kp:ku:2) & + pchild%u(ibx:ibz,ip:iu:2,jl:ju:2,kp:ku:2) & + pchild%u(ibx:ibz,il:iu:2,jp:ju:2,kp:ku:2) & + pchild%u(ibx:ibz,ip:iu:2,jp:ju:2,kp:ku:2)) #ifdef GLM pparent%u(iph ,is:it,js:jt,ks:kt) = & 0.125 * (pchild%u(iph ,il:iu:2,jl:ju:2,kl:ku:2) & + pchild%u(iph ,ip:iu:2,jl:ju:2,kl:ku:2) & + pchild%u(iph ,il:iu:2,jp:ju:2,kl:ku:2) & + pchild%u(iph ,ip:iu:2,jp:ju:2,kl:ku:2) & + pchild%u(iph ,il:iu:2,jl:ju:2,kp:ku:2) & + pchild%u(iph ,ip:iu:2,jl:ju:2,kp:ku:2) & + pchild%u(iph ,il:iu:2,jp:ju:2,kp:ku:2) & + pchild%u(iph ,ip:iu:2,jp:ju:2,kp:ku:2)) #endif /* GLM */ #endif /* MHD */ #endif /* NDIMS == 3 */ end do ! !------------------------------------------------------------------------------- ! end subroutine restrict_block ! !=============================================================================== ! ! clears_mesh: subroutine deallocates mesh, removing blocks ! !=============================================================================== ! subroutine clear_mesh() use error , only : print_info use mpitools, only : master implicit none !------------------------------------------------------------------------------- ! ! close the handler of the mesh statistics file ! if (master) close(funit) !------------------------------------------------------------------------------- ! end subroutine clear_mesh ! !=============================================================================== ! ! store_mesh_stats: subroutine stores mesh statistics ! !=============================================================================== ! subroutine store_mesh_stats(n, t) use blocks , only : block_meta, list_meta use blocks , only : get_mblocks, get_nleafs use config , only : ncells, nghost, toplev use coordinates, only : effres use mpitools, only : master, nprocs implicit none ! arguments ! integer , intent(in) :: n real(kind=8), intent(in) :: t ! local variables ! integer(kind=4) :: l real(kind=4) :: cov, eff ! local arrays ! integer(kind=4), dimension(toplev) :: ldist #ifdef MPI integer(kind=4), dimension(nprocs) :: cdist #endif /* MPI */ ! local pointers ! type(block_meta), pointer :: pmeta ! local variables ! integer(kind=4), save :: nm = 0, nl = 0 !------------------------------------------------------------------------------- ! ! store the statistics about mesh ! if (master) then if (nm .ne. get_mblocks() .or. nl .ne. get_nleafs()) then ! set new numbers of meta blocks and leafs ! nm = get_mblocks() nl = get_nleafs() ! calculate the coverage ! cov = (1.0 * nl) / tblocks eff = (1.0 * nl * (ncells + 2 * nghost)**NDIMS) & / product(effres(1:NDIMS) + 2 * nghost) ! get the block level distribution ! ldist(:) = 0 #ifdef MPI cdist(:) = 0 #endif /* MPI */ pmeta => list_meta do while(associated(pmeta)) if (pmeta%leaf) then ldist(pmeta%level) = ldist(pmeta%level) + 1 #ifdef MPI cdist(pmeta%cpu+1) = cdist(pmeta%cpu+1) + 1 #endif /* MPI */ end if pmeta => pmeta%next end do ! write down the statistics ! write(funit,"(2x,i8,2x,1pe14.8,2(2x,i6),2(2x,1pe14.8))",advance="no") & n, t, nl, nm, cov, eff write(funit,"(' ')",advance="no") do l = 1, toplev write(funit,"(2x,i6)",advance="no") ldist(l) end do #ifdef MPI write(funit,"(' ')",advance="no") do l = 1, nprocs write(funit,"(2x,i6)",advance="no") cdist(l) end do #endif /* MPI */ write(funit,"('')") end if ! number of blocks or leafs changed end if ! is master ! !------------------------------------------------------------------------------- ! end subroutine store_mesh_stats #ifdef DEBUG ! !=============================================================================== ! ! check_mesh: subroutine checks if the block structure is correct ! !=============================================================================== ! subroutine check_mesh(string) use blocks, only : block_meta, list_meta use blocks, only : check_metablock implicit none ! input arguments ! character(len=*), intent(in) :: string ! local pointers ! type(block_meta), pointer :: pmeta !------------------------------------------------------------------------------- ! ! check meta blocks ! pmeta => list_meta do while(associated(pmeta)) ! check the current block ! call check_metablock(pmeta, string) pmeta => pmeta%next end do !------------------------------------------------------------------------------- ! end subroutine check_mesh #endif /* DEBUG */ !=============================================================================== ! end module