!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2021 Grzegorz Kowal !! !! 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: BOUNDARIES !! !! This module handles the boundary synchronization. !! !!****************************************************************************** ! module boundaries #ifdef MPI use blocks, only : pointer_info #endif /* MPI */ implicit none ! parameters corresponding to the boundary type ! integer, parameter :: bnd_periodic = 0 integer, parameter :: bnd_open = 1 integer, parameter :: bnd_outflow = 2 integer, parameter :: bnd_reflective = 3 integer, parameter :: bnd_gravity = 4 integer, parameter :: bnd_user = 5 ! variable to store boundary type flags ! integer, dimension(3,2), save :: bnd_type = bnd_periodic #ifdef MPI ! arrays to store information about blocks which need to be exchange between ! processes ! type(pointer_info), dimension(:,:), allocatable, save :: barray integer , dimension(:,:), allocatable, save :: bcount #endif /* MPI */ ! by default everything is private ! private ! declare public subroutines ! public :: initialize_boundaries, finalize_boundaries, print_boundaries public :: boundary_variables, boundary_fluxes public :: bnd_type, bnd_periodic !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_BOUNDARIES: ! -------------------------------- ! ! Subroutine initializes the module BOUNDARIES by setting its parameters. ! ! Arguments: ! ! status - return flag of the procedure execution status; ! !=============================================================================== ! subroutine initialize_boundaries(status) ! import external procedures and variables ! use coordinates, only : periodic #ifdef MPI use mpitools , only : npmax #endif /* MPI */ use parameters , only : get_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status ! module parameters for the boundary update order and boundary type ! character(len = 32) :: xlbndry = "periodic" character(len = 32) :: xubndry = "periodic" character(len = 32) :: ylbndry = "periodic" character(len = 32) :: yubndry = "periodic" character(len = 32) :: zlbndry = "periodic" character(len = 32) :: zubndry = "periodic" ! local variables ! integer :: n ! !------------------------------------------------------------------------------- ! status = 0 ! get runtime values for the boundary types ! call get_parameter("xlbndry", xlbndry) call get_parameter("xubndry", xubndry) call get_parameter("ylbndry", ylbndry) call get_parameter("yubndry", yubndry) call get_parameter("zlbndry", zlbndry) call get_parameter("zubndry", zubndry) ! fill the boundary type flags ! select case(xlbndry) case("open") bnd_type(1,1) = bnd_open case("outflow", "out") bnd_type(1,1) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(1,1) = bnd_reflective case("hydrostatic", "gravity") bnd_type(1,1) = bnd_gravity case("user", "custom") bnd_type(1,1) = bnd_user case default bnd_type(1,1) = bnd_periodic end select select case(xubndry) case("open") bnd_type(1,2) = bnd_open case("outflow", "out") bnd_type(1,2) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(1,2) = bnd_reflective case("hydrostatic", "gravity") bnd_type(1,2) = bnd_gravity case("user", "custom") bnd_type(1,2) = bnd_user case default bnd_type(1,2) = bnd_periodic end select select case(ylbndry) case("open") bnd_type(2,1) = bnd_open case("outflow", "out") bnd_type(2,1) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(2,1) = bnd_reflective case("hydrostatic", "gravity") bnd_type(2,1) = bnd_gravity case("user", "custom") bnd_type(2,1) = bnd_user case default bnd_type(2,1) = bnd_periodic end select select case(yubndry) case("open") bnd_type(2,2) = bnd_open case("outflow", "out") bnd_type(2,2) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(2,2) = bnd_reflective case("hydrostatic", "gravity") bnd_type(2,2) = bnd_gravity case("user", "custom") bnd_type(2,2) = bnd_user case default bnd_type(2,2) = bnd_periodic end select select case(zlbndry) case("open") bnd_type(3,1) = bnd_open case("outflow", "out") bnd_type(3,1) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(3,1) = bnd_reflective case("hydrostatic", "gravity") bnd_type(3,1) = bnd_gravity case("user", "custom") bnd_type(3,1) = bnd_user case default bnd_type(3,1) = bnd_periodic end select select case(zubndry) case("open") bnd_type(3,2) = bnd_open case("outflow", "out") bnd_type(3,2) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(3,2) = bnd_reflective case("hydrostatic", "gravity") bnd_type(3,2) = bnd_gravity case("user", "custom") bnd_type(3,2) = bnd_user case default bnd_type(3,2) = bnd_periodic end select ! set domain periodicity ! do n = 1, NDIMS periodic(n) = (bnd_type(n,1) == bnd_periodic) .and. & (bnd_type(n,2) == bnd_periodic) end do #ifdef MPI ! allocate the exchange arrays ! allocate(barray(0:npmax,0:npmax), bcount(0:npmax,0:npmax), stat = status) ! prepare the exchange arrays ! if (status == 0) call prepare_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine initialize_boundaries ! !=============================================================================== ! ! subroutine FINALIZE_BOUNDARIES: ! ------------------------------ ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine finalize_boundaries(status) ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(out) :: status ! !------------------------------------------------------------------------------- ! status = 0 #ifdef MPI ! release the exchange arrays ! call release_exchange_array() ! deallocate the exchange arrays ! deallocate(barray, bcount, stat = status) #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine finalize_boundaries ! !=============================================================================== ! ! subroutine PRINT_BOUNDARIES: ! --------------------------- ! ! Subroutine prints module parameters and setup. ! ! Arguments: ! ! verbose - flag determining if the subroutine should be verbose; ! !=============================================================================== ! subroutine print_boundaries(verbose) use helpers , only : print_section, print_parameter use parameters, only : get_parameter implicit none logical, intent(in) :: verbose character(len=32) :: xlbndry = "periodic" character(len=32) :: xubndry = "periodic" character(len=32) :: ylbndry = "periodic" character(len=32) :: yubndry = "periodic" character(len=32) :: zlbndry = "periodic" character(len=32) :: zubndry = "periodic" integer :: w !------------------------------------------------------------------------------- ! if (verbose) then call get_parameter("xlbndry", xlbndry) call get_parameter("xubndry", xubndry) call get_parameter("ylbndry", ylbndry) call get_parameter("yubndry", yubndry) call get_parameter("zlbndry", zlbndry) call get_parameter("zubndry", zubndry) call print_section(verbose, "Boundaries") w = max(0, max(len(trim(xlbndry)), len(trim(xubndry)))) w = max(w, max(len(trim(ylbndry)), len(trim(yubndry)))) #if NDIMS == 3 w = max(w, max(len(trim(zlbndry)), len(trim(zubndry)))) #endif /* NDIMS == 3 */ call print_parameter(verbose, "X-boundary", xlbndry, xubndry, w) call print_parameter(verbose, "Y-boundary", ylbndry, yubndry, w) #if NDIMS == 3 call print_parameter(verbose, "Z-boundary", zlbndry, zubndry, w) #endif /* NDIMS == 3 */ end if !------------------------------------------------------------------------------- ! end subroutine print_boundaries ! !=============================================================================== ! ! subroutine BOUNDARY_VARIABLES: ! ----------------------------- ! ! Subroutine updates the ghost zones of the data blocks from their neighbors ! or applies the specific boundary conditions. ! ! Arguments: ! ! t, dt - time and time increment; ! !=============================================================================== ! subroutine boundary_variables(t, dt) ! import external procedures and variables ! use blocks , only : ndims use coordinates, only : minlev, maxlev ! local variables are not implicit by default ! implicit none ! subroutine arguments ! real(kind=8), intent(in) :: t, dt ! local variables ! integer :: idir ! !------------------------------------------------------------------------------- ! #if NDIMS == 3 ! update face boundaries between blocks at the same levels ! do idir = 1, ndims call boundaries_face_copy(idir) end do ! idir #endif /* NDIMS == 3 */ ! update edge boundaries between blocks at the same levels ! do idir = 1, ndims call boundaries_edge_copy(idir) end do ! idir ! update corner boundaries between blocks at the same levels ! call boundaries_corner_copy() ! do prolongation and restriction only if blocks are at different levels ! if (minlev /= maxlev) then #if NDIMS == 3 ! restrict face boundaries from higher level blocks ! do idir = 1, ndims call boundaries_face_restrict(idir) end do ! idir #endif /* NDIMS == 3 */ ! restricts edge boundaries from block at higher level ! do idir = 1, ndims call boundaries_edge_restrict(idir) end do ! idir ! restricts corner boundaries from blocks at higher levels ! call boundaries_corner_restrict() ! update specific boundaries ! call boundaries_specific(t, dt) #if NDIMS == 3 ! prolong face boundaries from lower level blocks ! do idir = 1, ndims call boundaries_face_prolong(idir) end do ! idir #endif /* NDIMS == 3 */ ! prolongs edge boundaries from block at lower level ! do idir = 1, ndims call boundaries_edge_prolong(idir) end do ! idir ! prolong corner boundaries from blocks at lower levels ! call boundaries_corner_prolong() end if ! minlev /= maxlev ! update specific boundaries ! call boundaries_specific(t, dt) ! convert updated primitive variables to conservative ones in all ghost cells ! call update_ghost_cells() !------------------------------------------------------------------------------- ! end subroutine boundary_variables ! !=============================================================================== ! ! subroutine BOUNDARY_FLUXES: ! -------------------------- ! ! Subroutine updates the numerical fluxes from neighbors which lay on ! higher level. At higher levels the numerical fluxes are calculated with ! smaller error, since the resolution is higher, therefore we take those ! fluxes and restrict them down to the level of the updated block. ! ! !=============================================================================== ! subroutine boundary_fluxes() ! import external procedures and variables ! use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf #ifdef MPI use blocks , only : block_info, pointer_info #endif /* MPI */ use blocks , only : ndims, nsides use coordinates, only : minlev, maxlev use coordinates, only : nh => ncells_half use coordinates, only : nb, ne, nbm, nbp, nep use coordinates, only : adxi, adyi #if NDIMS == 3 use coordinates, only : adzi #endif /* NDIMS == 3 */ use equations , only : nf, ns use equations , only : idn, isl, isu #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_data), pointer :: pdata type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: n #if NDIMS == 2 integer :: m #endif /* NDIMS == 2 */ integer :: i , il , iu integer :: j , jl , ju integer :: k = 1, kl = 1, ku = 1 integer :: s #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: scount, rcount integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ #if NDIMS == 3 real(kind=8), dimension(nf,nh,nh) :: fh, df real(kind=8), dimension(nh,nh) :: sl, sr, ps #else /* NDIMS == 3 */ real(kind=8), dimension(nf,nh, 1) :: fh, df real(kind=8), dimension(nh, 1) :: sl, sr, ps #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! ! quit if all blocks are at the same level ! if (minlev == maxlev) return #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update the fluxes between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta and data block ! pmeta => pleaf%meta pdata => pmeta%data ! iterate over all dimensions ! do n = 1, ndims #if NDIMS == 2 m = 3 - n #endif /* NDIMS == 2 */ ! iterate over all corners ! #if NDIMS == 3 do k = 1, nsides #endif /* NDIMS == 3 */ do j = 1, nsides do i = 1, nsides ! associate pneigh with the neighbor ! #if NDIMS == 2 pneigh => pmeta%edges(i,j,m)%ptr #endif /* NDIMS == 2 */ #if NDIMS == 3 pneigh => pmeta%faces(i,j,k,n)%ptr #endif /* NDIMS == 3 */ ! process only if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor lays at higher level ! if (pneigh%level > pmeta%level) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! update the flux depending on the direction ! select case(n) case(1) ! prepare the boundary layer indices for X-direction flux ! if (j == 1) then jl = nb ju = nb + nh - 1 else jl = ne - nh + 1 ju = ne end if #if NDIMS == 3 if (k == 1) then kl = nb ku = nb + nh - 1 else kl = ne - nh + 1 ku = ne end if #endif /* NDIMS == 3 */ ! update the flux at the X-face of the block ! if (i == 1) then #if NDIMS == 3 fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fx(:,2,nb :ne:2,nb :ne:2) & + pneigh%data%fx(:,2,nbp:ne:2,nbp:ne:2)) & + (pneigh%data%fx(:,2,nbp:ne:2,nb :ne:2) & + pneigh%data%fx(:,2,nb :ne:2,nbp:ne:2))) #else /* NDIMS == 3 */ fh(:,:,:) = & 5.0d-01 * (pneigh%data%fx(:,2,nb :ne:2,:) & + pneigh%data%fx(:,2,nbp:ne:2,:)) #endif /* NDIMS == 3 */ df(:,:,:) = (fh(:,:,:) - pdata%fx(:,1,jl:ju,kl:ku)) & * adxi(pmeta%level) pdata%du(:,nbm,jl:ju,kl:ku) = & pdata%du(:,nbm,jl:ju,kl:ku) - df(:,:,:) pdata%du(:,nb ,jl:ju,kl:ku) = & pdata%du(:,nb ,jl:ju,kl:ku) + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,nbm,jl:ju,kl:ku) & + sr(:,:) * pdata%q(s,nb ,jl:ju,kl:ku) pdata%du(s,nbm,jl:ju,kl:ku) = & pdata%du(s,nbm,jl:ju,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,nb ,jl:ju,kl:ku) = & pdata%du(s,nb ,jl:ju,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if else #if NDIMS == 3 fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fx(:,1,nb :ne:2,nb :ne:2) & + pneigh%data%fx(:,1,nbp:ne:2,nbp:ne:2)) & + (pneigh%data%fx(:,1,nbp:ne:2,nb :ne:2) & + pneigh%data%fx(:,1,nb :ne:2,nbp:ne:2))) #else /* NDIMS == 3 */ fh(:,:,:) = & 5.0d-01 * (pneigh%data%fx(:,1,nb :ne:2,:) & + pneigh%data%fx(:,1,nbp:ne:2,:)) #endif /* NDIMS == 3 */ df(:,:,:) = (fh(:,:,:) - pdata%fx(:,2,jl:ju,kl:ku)) & * adxi(pmeta%level) pdata%du(:,ne ,jl:ju,kl:ku) = & pdata%du(:,ne ,jl:ju,kl:ku) - df(:,:,:) pdata%du(:,nep,jl:ju,kl:ku) = & pdata%du(:,nep,jl:ju,kl:ku) + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,ne ,jl:ju,kl:ku) & + sr(:,:) * pdata%q(s,nep,jl:ju,kl:ku) pdata%du(s,ne ,jl:ju,kl:ku) = & pdata%du(s,ne ,jl:ju,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,nep,jl:ju,kl:ku) = & pdata%du(s,nep,jl:ju,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if end if case(2) ! prepare the boundary layer indices for Y-direction flux ! if (i == 1) then il = nb iu = nb + nh - 1 else il = ne - nh + 1 iu = ne end if #if NDIMS == 3 if (k == 1) then kl = nb ku = nb + nh - 1 else kl = ne - nh + 1 ku = ne end if #endif /* NDIMS == 3 */ ! update the flux at the Y-face of the block ! if (j == 1) then #if NDIMS == 3 fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fy(:,nb :ne:2,2,nb :ne:2) & + pneigh%data%fy(:,nbp:ne:2,2,nbp:ne:2)) & + (pneigh%data%fy(:,nbp:ne:2,2,nb :ne:2) & + pneigh%data%fy(:,nb :ne:2,2,nbp:ne:2))) #else /* NDIMS == 3 */ fh(:,:,:) = & 5.0d-01 * (pneigh%data%fy(:,nb :ne:2,2,:) & + pneigh%data%fy(:,nbp:ne:2,2,:)) #endif /* NDIMS == 3 */ df(:,:,:) = (fh(:,:,:) - pdata%fy(:,il:iu,1,kl:ku)) & * adyi(pmeta%level) pdata%du(:,il:iu,nbm,kl:ku) = & pdata%du(:,il:iu,nbm,kl:ku) - df(:,:,:) pdata%du(:,il:iu,nb ,kl:ku) = & pdata%du(:,il:iu,nb ,kl:ku) + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,il:iu,nbm,kl:ku) & + sr(:,:) * pdata%q(s,il:iu,nb ,kl:ku) pdata%du(s,il:iu,nbm,kl:ku) = & pdata%du(s,il:iu,nbm,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,nb ,kl:ku) = & pdata%du(s,il:iu,nb ,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if else #if NDIMS == 3 fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fy(:,nb :ne:2,1,nb :ne:2) & + pneigh%data%fy(:,nbp:ne:2,1,nbp:ne:2)) & + (pneigh%data%fy(:,nbp:ne:2,1,nb :ne:2) & + pneigh%data%fy(:,nb :ne:2,1,nbp:ne:2))) #else /* NDIMS == 3 */ fh(:,:,:) = & 5.0d-01 * (pneigh%data%fy(:,nb :ne:2,1,:) & + pneigh%data%fy(:,nbp:ne:2,1,:)) #endif /* NDIMS == 3 */ df(:,:,:) = (fh(:,:,:) - pdata%fy(:,il:iu,2,kl:ku)) & * adyi(pmeta%level) pdata%du(:,il:iu,ne ,kl:ku) = & pdata%du(:,il:iu,ne ,kl:ku) - df(:,:,:) pdata%du(:,il:iu,nep,kl:ku) = & pdata%du(:,il:iu,nep,kl:ku) + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,il:iu,ne ,kl:ku) & + sr(:,:) * pdata%q(s,il:iu,nep,kl:ku) pdata%du(s,il:iu,ne ,kl:ku) = & pdata%du(s,il:iu,ne ,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,nep,kl:ku) = & pdata%du(s,il:iu,nep,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if end if #if NDIMS == 3 case(3) ! prepare the boundary layer indices for Z-direction flux ! if (i == 1) then il = nb iu = nb + nh - 1 else il = ne - nh + 1 iu = ne end if if (j == 1) then jl = nb ju = nb + nh - 1 else jl = ne - nh + 1 ju = ne end if ! update the flux at the Z-face of the block ! if (k == 1) then fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fz(:,nb :ne:2,nb :ne:2,2) & + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,2)) & + (pneigh%data%fz(:,nbp:ne:2,nb :ne:2,2) & + pneigh%data%fz(:,nb :ne:2,nbp:ne:2,2))) df(:,:,:) = (fh(:,:,:) - pdata%fz(:,il:iu,jl:ju,1)) & * adzi(pmeta%level) pdata%du(:,il:iu,jl:ju,nbm) = & pdata%du(:,il:iu,jl:ju,nbm) - df(:,:,:) pdata%du(:,il:iu,jl:ju,nb ) = & pdata%du(:,il:iu,jl:ju,nb ) + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,nbm) & + sr(:,:) * pdata%q(s,il:iu,jl:ju,nb ) pdata%du(s,il:iu,jl:ju,nbm) = & pdata%du(s,il:iu,jl:ju,nbm) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,jl:ju,nb ) = & pdata%du(s,il:iu,jl:ju,nb ) & + df(idn,:,:) * ps(:,:) end do end if else fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fz(:,nb :ne:2,nb :ne:2,1) & + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,1)) & + (pneigh%data%fz(:,nbp:ne:2,nb :ne:2,1) & + pneigh%data%fz(:,nb :ne:2,nbp:ne:2,1))) df(:,:,:) = (fh(:,:,:) - pdata%fz(:,il:iu,jl:ju,2)) & * adzi(pmeta%level) pdata%du(:,il:iu,jl:ju,ne ) = & pdata%du(:,il:iu,jl:ju,ne ) - df(:,:,:) pdata%du(:,il:iu,jl:ju,nep) = & pdata%du(:,il:iu,jl:ju,nep) + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, fh(idn,:,:)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,ne ) & + sr(:,:) * pdata%q(s,il:iu,jl:ju,nep) pdata%du(s,il:iu,jl:ju,ne ) = & pdata%du(s,il:iu,jl:ju,ne ) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,jl:ju,nep) = & pdata%du(s,il:iu,jl:ju,nep) & + df(idn,:,:) * ps(:,:) end do end if end if #endif /* NDIMS == 3 */ end select #ifdef MPI end if ! pneigh on the current process else ! pneigh%proc /= pmeta%proc ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, n, (/ i, j, k /)) end if ! pneigh%proc /= pmeta%proc #endif /* MPI */ end if ! pmeta level < pneigh level end if ! pneigh associated end do ! i = 1, nsides end do ! j = 1, nsides #if NDIMS == 3 end do ! k = 1, nsides #endif /* NDIMS == 3 */ end do ! n = 1, ndims ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI ! update flux boundaries between neighbors laying on different processes ! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member) ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) end if ! get the number of blocks to exchange ! scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) ! process only pairs which have anything to exchange ! if ((scount + rcount) > 0) then ! allocate buffers for variable exchange ! #if NDIMS == 3 allocate(sbuf(nf,nh,nh,scount)) allocate(rbuf(nf,nh,nh,rcount)) #else /* NDIMS == 3 */ allocate(sbuf(nf,nh, 1,scount)) allocate(rbuf(nf,nh, 1,rcount)) #endif /* NDIMS == 3 */ !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate pinfo with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan all blocks on the list ! do while(associated(pinfo)) ! increase the block count ! l = l + 1 ! associate pneigh pointer ! pneigh => pinfo%neigh ! get neighbor direction and corner coordinates ! n = pinfo%direction i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! update directional flux from the neighbor ! select case(n) case(1) ! update the flux edge from the neighbor at higher level ! if (i == 1) then #if NDIMS == 3 sbuf(:,:,:,l) = & 2.5d-01 * ((pneigh%data%fx(:,2,nb :ne:2,nb :ne:2) & + pneigh%data%fx(:,2,nbp:ne:2,nbp:ne:2)) & + (pneigh%data%fx(:,2,nbp:ne:2,nb :ne:2) & + pneigh%data%fx(:,2,nb :ne:2,nbp:ne:2))) #else /* NDIMS == 3 */ sbuf(:,:,:,l) = & 5.0d-01 * (pneigh%data%fx(:,2,nb :ne:2,:) & + pneigh%data%fx(:,2,nbp:ne:2,:)) #endif /* NDIMS == 3 */ else #if NDIMS == 3 sbuf(:,:,:,l) = & 2.5d-01 * ((pneigh%data%fx(:,1,nb :ne:2,nb :ne:2) & + pneigh%data%fx(:,1,nbp:ne:2,nbp:ne:2)) & + (pneigh%data%fx(:,1,nbp:ne:2,nb :ne:2) & + pneigh%data%fx(:,1,nb :ne:2,nbp:ne:2))) #else /* NDIMS == 3 */ sbuf(:,:,:,l) = & 5.0d-01 * (pneigh%data%fx(:,1,nb :ne:2,:) & + pneigh%data%fx(:,1,nbp:ne:2,:)) #endif /* NDIMS == 3 */ end if case(2) ! update the flux edge from the neighbor at higher level ! if (j == 1) then #if NDIMS == 3 sbuf(:,:,:,l) = & 2.5d-01 * ((pneigh%data%fy(:,nb :ne:2,2,nb :ne:2) & + pneigh%data%fy(:,nbp:ne:2,2,nbp:ne:2)) & + (pneigh%data%fy(:,nbp:ne:2,2,nb :ne:2) & + pneigh%data%fy(:,nb :ne:2,2,nbp:ne:2))) #else /* NDIMS == 3 */ sbuf(:,:,:,l) = & 5.0d-01 * (pneigh%data%fy(:,nb :ne:2,2,:) & + pneigh%data%fy(:,nbp:ne:2,2,:)) #endif /* NDIMS == 3 */ else #if NDIMS == 3 sbuf(:,:,:,l) = & 2.5d-01 * ((pneigh%data%fy(:,nb :ne:2,1,nb :ne:2) & + pneigh%data%fy(:,nbp:ne:2,1,nbp:ne:2)) & + (pneigh%data%fy(:,nbp:ne:2,1,nb :ne:2) & + pneigh%data%fy(:,nb :ne:2,1,nbp:ne:2))) #else /* NDIMS == 3 */ sbuf(:,:,:,l) = & 5.0d-01 * (pneigh%data%fy(:,nb :ne:2,1,:) & + pneigh%data%fy(:,nbp:ne:2,1,:)) #endif /* NDIMS == 3 */ end if #if NDIMS == 3 case(3) ! update the flux edge from the neighbor at higher level ! if (k == 1) then sbuf(:,:,:,l) = & 2.5d-01 * ((pneigh%data%fz(:,nb :ne:2,nb :ne:2,2) & + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,2)) & + (pneigh%data%fz(:,nbp:ne:2,nb :ne:2,2) & + pneigh%data%fz(:,nb :ne:2,nbp:ne:2,2))) else sbuf(:,:,:,l) = & 2.5d-01 * ((pneigh%data%fz(:,nb :ne:2,nb :ne:2,1) & + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,1)) & + (pneigh%data%fz(:,nbp:ne:2,nb :ne:2,1) & + pneigh%data%fz(:,nb :ne:2,nbp:ne:2,1))) end if #endif /* NDIMS == 3 */ end select ! associate pinfo with the next block ! pinfo => pinfo%prev end do ! %ptr blocks !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, sbuf, rbuf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate pinfo with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! scan all blocks on the list ! do while(associated(pinfo)) ! increase the block count ! l = l + 1 ! associate meta and data block pointers ! pmeta => pinfo%meta pdata => pmeta%data ! get neighbor direction and corner indices ! n = pinfo%direction i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! update directional flux from the neighbor ! select case(n) case(1) ! prepare the boundary layer indices depending on the corner position ! if (j == 1) then jl = nb ju = nb + nh - 1 else jl = ne - nh + 1 ju = ne end if #if NDIMS == 3 if (k == 1) then kl = nb ku = nb + nh - 1 else kl = ne - nh + 1 ku = ne end if #endif /* NDIMS == 3 */ ! update the flux edge from the neighbor at higher level ! if (i == 1) then df(:,:,:) = (rbuf(:,:,:,l) - pdata%fx(:,1,jl:ju,kl:ku)) & * adxi(pmeta%level) pdata%du(:,nbm,jl:ju,kl:ku) = & pdata%du(:,nbm,jl:ju,kl:ku) - df(:,:,:) pdata%du(:,nb ,jl:ju,kl:ku) = & pdata%du(:,nb ,jl:ju,kl:ku) + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,nbm,jl:ju,kl:ku) & + sr(:,:) * pdata%q(s,nb ,jl:ju,kl:ku) pdata%du(s,nbm,jl:ju,kl:ku) = pdata%du(s,nbm,jl:ju,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,nb ,jl:ju,kl:ku) = pdata%du(s,nb ,jl:ju,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if else df(:,:,:) = (rbuf(:,:,:,l) - pdata%fx(:,2,jl:ju,kl:ku)) & * adxi(pmeta%level) pdata%du(:,ne ,jl:ju,kl:ku) = & pdata%du(:,ne ,jl:ju,kl:ku) - df(:,:,:) pdata%du(:,nep,jl:ju,kl:ku) = & pdata%du(:,nep,jl:ju,kl:ku) + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,ne ,jl:ju,kl:ku) & + sr(:,:) * pdata%q(s,nep,jl:ju,kl:ku) pdata%du(s,ne ,jl:ju,kl:ku) = pdata%du(s,ne ,jl:ju,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,nep,jl:ju,kl:ku) = pdata%du(s,nep,jl:ju,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if end if case(2) ! prepare the boundary layer indices depending on the corner position ! if (i == 1) then il = nb iu = nb + nh - 1 else il = ne - nh + 1 iu = ne end if #if NDIMS == 3 if (k == 1) then kl = nb ku = nb + nh - 1 else kl = ne - nh + 1 ku = ne end if #endif /* NDIMS == 3 */ ! update the flux edge from the neighbor at higher level ! if (j == 1) then df(:,:,:) = (rbuf(:,:,:,l) - pdata%fy(:,il:iu,1,kl:ku)) & * adyi(pmeta%level) pdata%du(:,il:iu,nbm,kl:ku) = pdata%du(:,il:iu,nbm,kl:ku) & - df(:,:,:) pdata%du(:,il:iu,nb ,kl:ku) = pdata%du(:,il:iu,nb ,kl:ku) & + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,il:iu,nbm,kl:ku) & + sr(:,:) * pdata%q(s,il:iu,nb ,kl:ku) pdata%du(s,il:iu,nbm,kl:ku) = pdata%du(s,il:iu,nbm,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,nb ,kl:ku) = pdata%du(s,il:iu,nb ,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if else df(:,:,:) = (rbuf(:,:,:,l) - pdata%fy(:,il:iu,2,kl:ku)) & * adyi(pmeta%level) pdata%du(:,il:iu,ne ,kl:ku) = pdata%du(:,il:iu,ne ,kl:ku) & - df(:,:,:) pdata%du(:,il:iu,nep,kl:ku) = pdata%du(:,il:iu,nep,kl:ku) & + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,il:iu,ne ,kl:ku) & + sr(:,:) * pdata%q(s,il:iu,nep,kl:ku) pdata%du(s,il:iu,ne ,kl:ku) = pdata%du(s,il:iu,ne ,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,nep,kl:ku) = pdata%du(s,il:iu,nep,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if end if #if NDIMS == 3 case(3) ! prepare the boundary layer indices depending on the corner position ! if (i == 1) then il = nb iu = nb + nh - 1 else il = ne - nh + 1 iu = ne end if if (j == 1) then jl = nb ju = nb + nh - 1 else jl = ne - nh + 1 ju = ne end if ! update the flux edge from the neighbor at higher level ! if (k == 1) then df(:,:,:) = (rbuf(:,:,:,l) - pdata%fz(:,il:iu,jl:ju,1)) & * adzi(pmeta%level) pdata%du(:,il:iu,jl:ju,nbm) = pdata%du(:,il:iu,jl:ju,nbm) & - df(:,:,:) pdata%du(:,il:iu,jl:ju,nb ) = pdata%du(:,il:iu,jl:ju,nb ) & + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,nbm) & + sr(:,:) * pdata%q(s,il:iu,jl:ju,nb ) pdata%du(s,il:iu,jl:ju,nbm) = pdata%du(s,il:iu,jl:ju,nbm) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,jl:ju,nb ) = pdata%du(s,il:iu,jl:ju,nb ) & + df(idn,:,:) * ps(:,:) end do end if else df(:,:,:) = (rbuf(:,:,:,l) - pdata%fz(:,il:iu,jl:ju,2)) & * adzi(pmeta%level) pdata%du(:,il:iu,jl:ju,ne ) = pdata%du(:,il:iu,jl:ju,ne ) & - df(:,:,:) pdata%du(:,il:iu,jl:ju,nep) = pdata%du(:,il:iu,jl:ju,nep) & + df(:,:,:) if (ns > 0) then sl(:,:) = sign(5.0d-01, rbuf(idn,:,:,l)) + 5.0d-01 sr(:,:) = 1.0d+00 - sl(:,:) do s = isl, isu ps(:,:) = sl(:,:) * pdata%q(s,il:iu,jl:ju,ne ) & + sr(:,:) * pdata%q(s,il:iu,jl:ju,nep) pdata%du(s,il:iu,jl:ju,ne ) = pdata%du(s,il:iu,jl:ju,ne ) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,jl:ju,nep) = pdata%du(s,il:iu,jl:ju,nep) & + df(idn,:,:) * ps(:,:) end do end if end if #endif /* NDIMS == 3 */ end select ! associate pinfo with the next block ! pinfo => pinfo%prev end do ! %ptr blocks ! deallocate data buffer ! deallocate(sbuf, rbuf) end if ! (scount + rcount) > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundary_fluxes ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! DOMAIN SPECIFIC BOUNDARY SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BOUNDARIES_SPECIFIC: ! ------------------------------ ! ! Subroutine scans over all leaf blocks in order to find blocks without ! neighbors and update the corresponding boundaries for the selected ! boundary type. ! ! Arguments: ! ! t, dt - time and time increment; ! !=============================================================================== ! subroutine boundaries_specific(t, dt) ! import external procedures and variables ! use blocks , only : block_meta, block_leaf use blocks , only : list_leaf use blocks , only : ndims, nsides use coordinates, only : nn => bcells use coordinates, only : ax, ay #if NDIMS == 3 use coordinates, only : az #endif /* NDIMS == 3 */ use coordinates, only : periodic #ifdef MPI use mpitools , only : nproc #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! real(kind=8), intent(in) :: t, dt ! local pointers ! type(block_meta), pointer :: pmeta type(block_leaf), pointer :: pleaf ! local variables ! integer :: i, j, k = 1, n #if NDIMS == 2 integer :: m #endif /* NDIMS == 2 */ ! local arrays ! real(kind=8), dimension(nn) :: x real(kind=8), dimension(nn) :: y #if NDIMS == 3 real(kind=8), dimension(nn) :: z #else /* NDIMS == 3 */ real(kind=8), dimension( 1) :: z #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! process only if this block is marked for update ! if (pmeta%update) then #ifdef MPI ! check if the block belongs to the local process ! if (pmeta%process == nproc) then #endif /* MPI */ ! prepare block coordinates ! x(:) = pmeta%xmin + ax(pmeta%level,:) y(:) = pmeta%ymin + ay(pmeta%level,:) #if NDIMS == 3 z(:) = pmeta%zmin + az(pmeta%level,:) #else /* NDIMS == 3 */ z( : ) = 0.0d+00 #endif /* NDIMS == 3 */ #if NDIMS == 2 ! iterate over all directions ! do n = 1, ndims ! process boundaries only if they are not periodic along the given direction ! if (.not. periodic(n)) then ! calculate the edge direction (in 2D we don't have face neighbors, so we have ! to use edge neighbors) ! m = 3 - n ! iterate over all corners ! do j = 1, nsides do i = 1, nsides ! if the face neighbor is not associated, apply specific boundaries ! if (.not. associated(pmeta%edges(i,j,m)%ptr)) & call block_boundary_specific(n, (/ i, j, k /) & , t, dt, x(:), y(:), z(:) & , pmeta%data%q(:,:,:,:)) end do ! i = 1, sides end do ! j = 1, sides end if ! not periodic end do ! n = 1, ndims #endif /* NDIMS == 2 */ #if NDIMS == 3 ! iterate over all directions ! do n = 1, ndims ! process boundaries only if they are not periodic along the given direction ! if (.not. periodic(n)) then ! iterate over all corners ! do k = 1, nsides do j = 1, nsides do i = 1, nsides ! if the face neighbor is not associated, apply specific boundaries ! if (.not. associated(pmeta%faces(i,j,k,n)%ptr)) & call block_boundary_specific(n, (/ i, j, k /) & , t, dt, x(:), y(:), z(:) & , pmeta%data%q(:,:,:,:)) end do ! i = 1, sides end do ! j = 1, sides end do ! k = 1, sides end if ! not periodic end do ! n = 1, ndims #endif /* NDIMS == 3 */ #ifdef MPI end if ! block belongs to the local process #endif /* MPI */ end if ! if pmeta marked for update ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks !------------------------------------------------------------------------------- ! end subroutine boundaries_specific #if NDIMS == 3 ! !=============================================================================== ! ! DOMAIN FACE BOUNDARY UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BOUNDARIES_FACE_COPY: ! ------------------------------- ! ! Subroutine updates the face boundaries between blocks on the same level. ! ! Arguments: ! ! idir - the direction to be processed; ! !=============================================================================== ! subroutine boundaries_face_copy(idir) ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : nh => ncells_half, ng => nghosts use coordinates, only : faces_gc, faces_dc use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(in) :: idir ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i , j , k = 1 integer :: il, jl, kl integer :: iu, ju, ku integer :: is, js, ks integer :: it, jt, kt #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: ecount = 0 integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: buf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! do k = 1, nsides do j = 1, nsides do i = 1, nsides ! associate pneigh with the current neighbor ! pneigh => pmeta%faces(i,j,k,idir)%ptr ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor is at the same level ! if (pneigh%level == pmeta%level) then ! process only blocks and neighbors which are marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare region indices of the block and its neighbor for the face boundary ! update ! il = faces_gc(i,j,k,idir)%l(1) jl = faces_gc(i,j,k,idir)%l(2) kl = faces_gc(i,j,k,idir)%l(3) iu = faces_gc(i,j,k,idir)%u(1) ju = faces_gc(i,j,k,idir)%u(2) ku = faces_gc(i,j,k,idir)%u(3) is = faces_dc(i,j,k,idir)%l(1) js = faces_dc(i,j,k,idir)%l(2) ks = faces_dc(i,j,k,idir)%l(3) it = faces_dc(i,j,k,idir)%u(1) jt = faces_dc(i,j,k,idir)%u(2) kt = faces_dc(i,j,k,idir)%u(3) ! copy the corresponding face region from the neighbor to the current data ! block ! pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & pneigh%data%q(1:nv,is:it,js:jt,ks:kt) #ifdef MPI end if ! pneigh on the current process else ! block and neighbor belong to different processes ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, idir, (/ i, j, k /)) end if ! block and neighbor belong to different processes #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at the same level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides end do ! k = 1, nsides ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member), ! and the number of blocks to exchange ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) ecount = bcount(sproc,rproc) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) ecount = bcount(rproc,sproc) end if ! process only pairs which have anything to exchange ! if (ecount > 0) then ! allocate data buffer for variables to exchange ! select case(idir) case(1) allocate(buf(ecount,nv,ng,nh,nh)) case(2) allocate(buf(ecount,nv,nh,ng,nh)) case(3) allocate(buf(ecount,nv,nh,nh,ng)) end select !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pneigh to the associated neighbor block ! pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) k = pinfo%corner(3) ! prepare region indices for the face boundary update ! is = faces_dc(i,j,k,idir)%l(1) js = faces_dc(i,j,k,idir)%l(2) ks = faces_dc(i,j,k,idir)%l(3) it = faces_dc(i,j,k,idir)%u(1) jt = faces_dc(i,j,k,idir)%u(2) kt = faces_dc(i,j,k,idir)%u(3) ! copy the corresponding face region from the neighbor and insert it to ! the buffer ! select case(idir) case(1) buf(l,1:nv,1:ng,1:nh,1:nh) = pneigh%data%q(1:nv,is:it,js:jt,ks:kt) case(2) buf(l,1:nv,1:nh,1:ng,1:nh) = pneigh%data%q(1:nv,is:it,js:jt,ks:kt) case(3) buf(l,1:nv,1:nh,1:nh,1:ng) = pneigh%data%q(1:nv,is:it,js:jt,ks:kt) end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, buf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) k = pinfo%corner(3) ! prepare region indices for the face boundary update ! il = faces_gc(i,j,k,idir)%l(1) jl = faces_gc(i,j,k,idir)%l(2) kl = faces_gc(i,j,k,idir)%l(3) iu = faces_gc(i,j,k,idir)%u(1) ju = faces_gc(i,j,k,idir)%u(2) ku = faces_gc(i,j,k,idir)%u(3) ! update the corresponding face region of the current block ! select case(idir) case(1) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = buf(l,1:nv,1:ng,1:nh,1:nh) case(2) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = buf(l,1:nv,1:nh,1:ng,1:nh) case(3) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = buf(l,1:nv,1:nh,1:nh,1:ng) end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(buf) end if ! ecount > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_face_copy ! !=============================================================================== ! ! subroutine BOUNDARIES_FACE_RESTRICT: ! ----------------------------------- ! ! Subroutine updates the face boundaries from blocks on higher level. ! ! Arguments: ! ! idir - the direction to be processed; ! !=============================================================================== ! subroutine boundaries_face_restrict(idir) ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : nh => ncells_half, ng => nghosts use coordinates, only : faces_gr use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(in) :: idir ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i , j , k = 1 integer :: il, jl, kl integer :: iu, ju, ku #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: scount, rcount integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! do k = 1, nsides do j = 1, nsides do i = 1, nsides ! associate pneigh with the current neighbor ! pneigh => pmeta%faces(i,j,k,idir)%ptr ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor is at higher level ! if (pneigh%level > pmeta%level) then ! process only blocks and neighbors which are marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare region indices of the block and its neighbor for the face boundary ! update ! il = faces_gr(i,j,k,idir)%l(1) jl = faces_gr(i,j,k,idir)%l(2) kl = faces_gr(i,j,k,idir)%l(3) iu = faces_gr(i,j,k,idir)%u(1) ju = faces_gr(i,j,k,idir)%u(2) ku = faces_gr(i,j,k,idir)%u(3) ! extract the corresponding face region from the neighbor and insert it in ! the current data block ! call block_face_restrict(idir, i, j, k & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku)) #ifdef MPI end if ! pneigh on the current process else ! block and neighbor belong to different processes ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, idir, (/ i, j, k /)) end if ! block and neighbor belong to different processes #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at the same level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides end do ! k = 1, nsides ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member) ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) end if ! get the number of blocks to exchange ! scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) ! process only pairs which have anything to exchange ! if ((scount + rcount) > 0) then ! allocate data buffer for variables to exchange ! select case(idir) case(1) allocate(sbuf(scount,nv,ng,nh,nh)) allocate(rbuf(rcount,nv,ng,nh,nh)) case(2) allocate(sbuf(scount,nv,nh,ng,nh)) allocate(rbuf(rcount,nv,nh,ng,nh)) case(3) allocate(sbuf(scount,nv,nh,nh,ng)) allocate(rbuf(rcount,nv,nh,nh,ng)) end select !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pneigh to the associated neighbor block ! pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) k = pinfo%corner(3) ! restrict the corresponding face region from the neighbor and insert it ! to the buffer ! select case(idir) case(1) call block_face_restrict(idir, i, j, k & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:nh,1:nh)) case(2) call block_face_restrict(idir, i, j, k & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:nh,1:ng,1:nh)) case(3) call block_face_restrict(idir, i, j, k & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:nh,1:nh,1:ng)) end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, sbuf, rbuf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) k = pinfo%corner(3) ! prepare region indices for the face boundary update ! il = faces_gr(i,j,k,idir)%l(1) jl = faces_gr(i,j,k,idir)%l(2) kl = faces_gr(i,j,k,idir)%l(3) iu = faces_gr(i,j,k,idir)%u(1) ju = faces_gr(i,j,k,idir)%u(2) ku = faces_gr(i,j,k,idir)%u(3) ! update the corresponding face region of the current block ! select case(idir) case(1) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:ng,1:nh,1:nh) case(2) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:nh,1:ng,1:nh) case(3) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:nh,1:nh,1:ng) end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(sbuf, rbuf) end if ! (scount + rcount) > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_face_restrict ! !=============================================================================== ! ! subroutine BOUNDARIES_FACE_PROLONG: ! ---------------------------------- ! ! Subroutine updates the face boundaries from blocks on lower level. ! ! Arguments: ! ! idir - the direction to be processed; ! !=============================================================================== ! subroutine boundaries_face_prolong(idir) ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : ni => ncells, ng => nghosts use coordinates, only : faces_gp use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(in) :: idir ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i , j , k = 1 integer :: ic, jc, kc integer :: ih, jh, kh integer :: il = 1, jl = 1, kl = 1 integer :: iu = 1, ju = 1, ku = 1 #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: scount, rcount integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! ! calculate the sizes ! ih = ni + ng jh = ni + ng kh = ni + ng #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! do k = 1, nsides kc = k do j = 1, nsides jc = j do i = 1, nsides ic = i ! associate pneigh with the current neighbor ! pneigh => pmeta%faces(i,j,k,idir)%ptr ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor lays at lower level ! if (pneigh%level < pmeta%level) then ! process only blocks and neighbors which are marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare indices of the region in which the boundaries should be updated ! select case(idir) case(1) jc = pmeta%pos(2) + 1 kc = pmeta%pos(3) + 1 il = faces_gp(i ,jc,kc,idir)%l(1) jl = faces_gp(i ,jc,kc,idir)%l(2) kl = faces_gp(i ,jc,kc,idir)%l(3) iu = faces_gp(i ,jc,kc,idir)%u(1) ju = faces_gp(i ,jc,kc,idir)%u(2) ku = faces_gp(i ,jc,kc,idir)%u(3) case(2) ic = pmeta%pos(1) + 1 kc = pmeta%pos(3) + 1 il = faces_gp(ic,j ,kc,idir)%l(1) jl = faces_gp(ic,j ,kc,idir)%l(2) kl = faces_gp(ic,j ,kc,idir)%l(3) iu = faces_gp(ic,j ,kc,idir)%u(1) ju = faces_gp(ic,j ,kc,idir)%u(2) ku = faces_gp(ic,j ,kc,idir)%u(3) case(3) ic = pmeta%pos(1) + 1 jc = pmeta%pos(2) + 1 il = faces_gp(ic,jc,k ,idir)%l(1) jl = faces_gp(ic,jc,k ,idir)%l(2) kl = faces_gp(ic,jc,k ,idir)%l(3) iu = faces_gp(ic,jc,k ,idir)%u(1) ju = faces_gp(ic,jc,k ,idir)%u(2) ku = faces_gp(ic,jc,k ,idir)%u(3) end select ! take the neighbor volume, extract the corresponding face region and insert ! it in the current data block ! call block_face_prolong(idir, ic, jc, kc & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku)) #ifdef MPI end if ! pneigh on the current process else ! block and neighbor belong to different processes ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, idir, (/ i, j, k /)) end if ! block and neighbor belong to different processes #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at lower level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides end do ! k = 1, nsides ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member) ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) end if ! get the number of blocks to exchange ! scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) ! process only pairs which have anything to exchange ! if ((scount + rcount) > 0) then ! allocate data buffer for variables to exchange ! select case(idir) case(1) allocate(sbuf(scount,nv,ng,jh,kh)) allocate(rbuf(rcount,nv,ng,jh,kh)) case(2) allocate(sbuf(scount,nv,ih,ng,kh)) allocate(rbuf(rcount,nv,ih,ng,kh)) case(3) allocate(sbuf(scount,nv,ih,jh,ng)) allocate(rbuf(rcount,nv,ih,jh,ng)) end select !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pmeta and pneigh to the right blocks ! pmeta => pinfo%meta pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) k = pinfo%corner(3) ! extract the corresponding face region from the neighbor and insert it ! to the buffer ! select case(idir) case(1) j = pmeta%pos(2) + 1 k = pmeta%pos(3) + 1 call block_face_prolong(idir, i, j, k & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:jh,1:kh)) case(2) i = pmeta%pos(1) + 1 k = pmeta%pos(3) + 1 call block_face_prolong(idir, i, j, k & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ih,1:ng,1:kh)) case(3) i = pmeta%pos(1) + 1 j = pmeta%pos(2) + 1 call block_face_prolong(idir, i, j, k & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ih,1:jh,1:ng)) end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, sbuf, rbuf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) k = pinfo%corner(3) ! update the corresponding face region of the current block ! select case(idir) case(1) jc = pmeta%pos(2) + 1 kc = pmeta%pos(3) + 1 il = faces_gp(i ,jc,kc,idir)%l(1) jl = faces_gp(i ,jc,kc,idir)%l(2) kl = faces_gp(i ,jc,kc,idir)%l(3) iu = faces_gp(i ,jc,kc,idir)%u(1) ju = faces_gp(i ,jc,kc,idir)%u(2) ku = faces_gp(i ,jc,kc,idir)%u(3) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = rbuf(l,1:nv,1:ng,1:jh,1:kh) case(2) ic = pmeta%pos(1) + 1 kc = pmeta%pos(3) + 1 il = faces_gp(ic,j ,kc,idir)%l(1) jl = faces_gp(ic,j ,kc,idir)%l(2) kl = faces_gp(ic,j ,kc,idir)%l(3) iu = faces_gp(ic,j ,kc,idir)%u(1) ju = faces_gp(ic,j ,kc,idir)%u(2) ku = faces_gp(ic,j ,kc,idir)%u(3) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = rbuf(l,1:nv,1:ih,1:ng,1:kh) case(3) ic = pmeta%pos(1) + 1 jc = pmeta%pos(2) + 1 il = faces_gp(ic,jc,k ,idir)%l(1) jl = faces_gp(ic,jc,k ,idir)%l(2) kl = faces_gp(ic,jc,k ,idir)%l(3) iu = faces_gp(ic,jc,k ,idir)%u(1) ju = faces_gp(ic,jc,k ,idir)%u(2) ku = faces_gp(ic,jc,k ,idir)%u(3) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = rbuf(l,1:nv,1:ih,1:jh,1:ng) end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(sbuf, rbuf) end if ! (scount + rcount) > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_face_prolong #endif /* NDIMS == 3 */ ! !=============================================================================== ! ! DOMAIN EDGE BOUNDARY UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BOUNDARIES_EDGE_COPY: ! ------------------------------- ! ! Subroutine updates the edge boundaries from blocks on the same level. ! ! Arguments: ! ! idir - the direction to be processed; ! !=============================================================================== ! subroutine boundaries_edge_copy(idir) ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : nh => ncells_half, ng => nghosts use coordinates, only : edges_gc, edges_dc use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(in) :: idir ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i, il, iu, is, it integer :: j, jl, ju, js, jt integer :: k = 1 #if NDIMS == 3 integer :: kl, ku, ks, kt #endif /* NDIMS == 3 */ #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: ecount = 0 integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: buf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! #if NDIMS == 3 do k = 1, nsides #endif /* NDIMS == 3 */ do j = 1, nsides do i = 1, nsides ! associate pneigh with the current neighbor ! #if NDIMS == 2 pneigh => pmeta%edges(i,j,idir)%ptr #endif /* NDIMS == 2 */ #if NDIMS == 3 pneigh => pmeta%edges(i,j,k,idir)%ptr #endif /* NDIMS == 3 */ ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor is at the same level ! if (pneigh%level == pmeta%level) then ! process only blocks and neighbors which are marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare region indices of the block and its neighbor for the edge boundary ! update #if NDIMS == 2 il = edges_gc(i,j ,idir)%l(1) jl = edges_gc(i,j ,idir)%l(2) iu = edges_gc(i,j ,idir)%u(1) ju = edges_gc(i,j ,idir)%u(2) is = edges_dc(i,j ,idir)%l(1) js = edges_dc(i,j ,idir)%l(2) it = edges_dc(i,j ,idir)%u(1) jt = edges_dc(i,j ,idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_gc(i,j,k,idir)%l(1) jl = edges_gc(i,j,k,idir)%l(2) kl = edges_gc(i,j,k,idir)%l(3) iu = edges_gc(i,j,k,idir)%u(1) ju = edges_gc(i,j,k,idir)%u(2) ku = edges_gc(i,j,k,idir)%u(3) is = edges_dc(i,j,k,idir)%l(1) js = edges_dc(i,j,k,idir)%l(2) ks = edges_dc(i,j,k,idir)%l(3) it = edges_dc(i,j,k,idir)%u(1) jt = edges_dc(i,j,k,idir)%u(2) kt = edges_dc(i,j,k,idir)%u(3) #endif /* NDIMS == 3 */ ! copy the corresponding edge region from the neighbor and insert it in ! the current data block ! #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = & pneigh%data%q(1:nv,is:it,js:jt, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & pneigh%data%q(1:nv,is:it,js:jt,ks:kt) #endif /* NDIMS == 3 */ #ifdef MPI end if ! pneigh on the current process else ! block and neighbor belong to different processes ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, idir, (/ i, j, k /)) end if ! block and neighbor belong to different processes #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at the same level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides #if NDIMS == 3 end do ! k = 1, nsides #endif /* NDIMS == 3 */ ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member), ! and the number of blocks to exchange ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) ecount = bcount(sproc,rproc) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) ecount = bcount(rproc,sproc) end if ! process only pairs which have anything to exchange ! if (ecount > 0) then ! allocate buffers for variable exchange ! select case(idir) #if NDIMS == 2 case(1) allocate(buf(ecount,nv,nh,ng, 1)) case(2) allocate(buf(ecount,nv,ng,nh, 1)) #endif /* NDIMS == 2 */ #if NDIMS == 3 case(1) allocate(buf(ecount,nv,nh,ng,ng)) case(2) allocate(buf(ecount,nv,ng,nh,ng)) case(3) allocate(buf(ecount,nv,ng,ng,nh)) #endif /* NDIMS == 3 */ end select !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pneigh to the associated neighbor block ! pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! prepare indices of the region for edge boundary update ! #if NDIMS == 2 is = edges_dc(i,j ,idir)%l(1) js = edges_dc(i,j ,idir)%l(2) it = edges_dc(i,j ,idir)%u(1) jt = edges_dc(i,j ,idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 is = edges_dc(i,j,k,idir)%l(1) js = edges_dc(i,j,k,idir)%l(2) ks = edges_dc(i,j,k,idir)%l(3) it = edges_dc(i,j,k,idir)%u(1) jt = edges_dc(i,j,k,idir)%u(2) kt = edges_dc(i,j,k,idir)%u(3) #endif /* NDIMS == 3 */ ! copy the corresponding edge region from the neighbor and insert it in ! the buffer ! select case(idir) case(1) #if NDIMS == 2 buf(l,1:nv,1:nh,1:ng, : ) = pneigh%data%q(1:nv,is:it,js:jt, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 buf(l,1:nv,1:nh,1:ng,1:ng) = pneigh%data%q(1:nv,is:it,js:jt,ks:kt) #endif /* NDIMS == 3 */ case(2) #if NDIMS == 2 buf(l,1:nv,1:ng,1:nh, : ) = pneigh%data%q(1:nv,is:it,js:jt, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 buf(l,1:nv,1:ng,1:nh,1:ng) = pneigh%data%q(1:nv,is:it,js:jt,ks:kt) #endif /* NDIMS == 3 */ #if NDIMS == 3 case(3) buf(l,1:nv,1:ng,1:ng,1:nh) = pneigh%data%q(1:nv,is:it,js:jt,ks:kt) #endif /* NDIMS == 3 */ end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, buf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! prepare indices of the region for the edge update ! #if NDIMS == 2 il = edges_gc(i,j ,idir)%l(1) jl = edges_gc(i,j ,idir)%l(2) iu = edges_gc(i,j ,idir)%u(1) ju = edges_gc(i,j ,idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_gc(i,j,k,idir)%l(1) jl = edges_gc(i,j,k,idir)%l(2) kl = edges_gc(i,j,k,idir)%l(3) iu = edges_gc(i,j,k,idir)%u(1) ju = edges_gc(i,j,k,idir)%u(2) ku = edges_gc(i,j,k,idir)%u(3) #endif /* NDIMS == 3 */ ! update the corresponding edge region of the current block ! select case(idir) case(1) #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = buf(l,1:nv,1:nh,1:ng, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = buf(l,1:nv,1:nh,1:ng,1:ng) #endif /* NDIMS == 3 */ case(2) #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = buf(l,1:nv,1:ng,1:nh, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = buf(l,1:nv,1:ng,1:nh,1:ng) #endif /* NDIMS == 3 */ #if NDIMS == 3 case(3) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = buf(l,1:nv,1:ng,1:ng,1:nh) #endif /* NDIMS == 3 */ end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(buf) end if ! ecount > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_edge_copy ! !=============================================================================== ! ! subroutine BOUNDARIES_EDGE_RESTRICT: ! ----------------------------------- ! ! Subroutine updates the edge boundaries from blocks on higher level. ! ! Arguments: ! ! idir - the direction to be processed; ! !=============================================================================== ! subroutine boundaries_edge_restrict(idir) ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : nh => ncells_half, ng => nghosts use coordinates, only : edges_gr use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(in) :: idir ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i, il, iu integer :: j, jl, ju integer :: k = 1 #if NDIMS == 3 integer :: kl, ku #endif /* NDIMS == 3 */ #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: scount, rcount integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! #if NDIMS == 3 do k = 1, nsides #endif /* NDIMS == 3 */ do j = 1, nsides do i = 1, nsides ! assign pneigh to the current neighbor ! #if NDIMS == 2 pneigh => pmeta%edges(i,j,idir)%ptr #endif /* NDIMS == 2 */ #if NDIMS == 3 pneigh => pmeta%edges(i,j,k,idir)%ptr #endif /* NDIMS == 3 */ ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor is at higher level ! if (pneigh%level > pmeta%level) then ! process only blocks and neighbors which are marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare the region indices for edge boundary update ! #if NDIMS == 2 il = edges_gr(i,j, idir)%l(1) jl = edges_gr(i,j, idir)%l(2) iu = edges_gr(i,j, idir)%u(1) ju = edges_gr(i,j, idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_gr(i,j,k,idir)%l(1) jl = edges_gr(i,j,k,idir)%l(2) kl = edges_gr(i,j,k,idir)%l(3) iu = edges_gr(i,j,k,idir)%u(1) ju = edges_gr(i,j,k,idir)%u(2) ku = edges_gr(i,j,k,idir)%u(3) #endif /* NDIMS == 3 */ ! extract the corresponding edge region from the neighbor to the current ! data block ! #if NDIMS == 2 call block_edge_restrict(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_edge_restrict(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku)) #endif /* NDIMS == 3 */ #ifdef MPI end if ! pneigh on the current process else ! block and neighbor belong to different processes ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, idir, (/ i, j, k /)) end if ! block and neighbor belong to different processes #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at the same level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides #if NDIMS == 3 end do ! k = 1, nsides #endif /* NDIMS == 3 */ ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member) ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) end if ! get the number of blocks to exchange ! scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) ! process only pairs which have anything to exchange ! if ((scount + rcount) > 0) then ! allocate buffers for variable exchange ! select case(idir) #if NDIMS == 2 case(1) allocate(sbuf(scount,nv,nh,ng, 1)) allocate(rbuf(rcount,nv,nh,ng, 1)) case(2) allocate(sbuf(scount,nv,ng,nh, 1)) allocate(rbuf(rcount,nv,ng,nh, 1)) #endif /* NDIMS == 2 */ #if NDIMS == 3 case(1) allocate(sbuf(scount,nv,nh,ng,ng)) allocate(rbuf(rcount,nv,nh,ng,ng)) case(2) allocate(sbuf(scount,nv,ng,nh,ng)) allocate(rbuf(rcount,nv,ng,nh,ng)) case(3) allocate(sbuf(scount,nv,ng,ng,nh)) allocate(rbuf(rcount,nv,ng,ng,nh)) #endif /* NDIMS == 3 */ end select !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pneigh to the associated neighbor block ! pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! extract the corresponding edge region from the neighbor and insert it ! to the buffer ! select case(idir) case(1) #if NDIMS == 2 call block_edge_restrict(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:nh,1:ng, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_edge_restrict(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:nh,1:ng,1:ng)) #endif /* NDIMS == 3 */ case(2) #if NDIMS == 2 call block_edge_restrict(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:nh, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_edge_restrict(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:nh,1:ng)) #endif /* NDIMS == 3 */ #if NDIMS == 3 case(3) call block_edge_restrict(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:ng,1:nh)) #endif /* NDIMS == 3 */ end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, sbuf, rbuf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! prepare the region indices for edge boundary update ! #if NDIMS == 2 il = edges_gr(i,j, idir)%l(1) jl = edges_gr(i,j, idir)%l(2) iu = edges_gr(i,j, idir)%u(1) ju = edges_gr(i,j, idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_gr(i,j,k,idir)%l(1) jl = edges_gr(i,j,k,idir)%l(2) kl = edges_gr(i,j,k,idir)%l(3) iu = edges_gr(i,j,k,idir)%u(1) ju = edges_gr(i,j,k,idir)%u(2) ku = edges_gr(i,j,k,idir)%u(3) #endif /* NDIMS == 3 */ ! update the corresponding corner region of the current block ! select case(idir) case(1) #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = & rbuf(l,1:nv,1:nh,1:ng, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:nh,1:ng,1:ng) #endif /* NDIMS == 3 */ case(2) #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = & rbuf(l,1:nv,1:ng,1:nh, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:ng,1:nh,1:ng) #endif /* NDIMS == 3 */ #if NDIMS == 3 case(3) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:ng,1:ng,1:nh) #endif /* NDIMS == 3 */ end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(sbuf, rbuf) end if ! (scount + rcount) > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_edge_restrict ! !=============================================================================== ! ! subroutine BOUNDARIES_EDGE_PROLONG: ! ---------------------------------- ! ! Subroutine updates the edge boundaries from blocks on lower level. ! ! Arguments: ! ! idir - the direction to be processed; ! !=============================================================================== ! subroutine boundaries_edge_prolong(idir) ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : ni => ncells, ng => nghosts use coordinates, only : edges_gp use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(in) :: idir ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i, ic, ih, il = 1, iu = 1 integer :: j, jc, jh, jl = 1, ju = 1 integer :: k = 1, kc = 1 #if NDIMS == 3 integer :: kh, kl = 1, ku = 1 #endif /* NDIMS == 3 */ #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: scount, rcount integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! ! calculate the sizes ! ih = ni + ng jh = ni + ng #if NDIMS == 3 kh = ni + ng #endif /* NDIMS == 3 */ #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! #if NDIMS == 3 do k = 1, nsides kc = k #endif /* NDIMS == 3 */ do j = 1, nsides jc = j do i = 1, nsides ic = i ! assign pneigh to the current neighbor ! #if NDIMS == 2 pneigh => pmeta%edges(i,j,idir)%ptr #endif /* NDIMS == 2 */ #if NDIMS == 3 pneigh => pmeta%edges(i,j,k,idir)%ptr #endif /* NDIMS == 3 */ ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor lays at lower level ! if (pneigh%level < pmeta%level) then ! process only blocks and neighbors which are marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare the region indices for edge boundary update ! select case(idir) case(1) ic = pmeta%pos(1) + 1 #if NDIMS == 2 il = edges_gp(ic,j ,idir)%l(1) iu = edges_gp(ic,j ,idir)%u(1) jl = edges_gp(i ,j ,idir)%l(2) ju = edges_gp(i ,j ,idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_gp(ic,j,k,idir)%l(1) iu = edges_gp(ic,j,k,idir)%u(1) jl = edges_gp(i ,j,k,idir)%l(2) ju = edges_gp(i ,j,k,idir)%u(2) kl = edges_gp(i ,j,k,idir)%l(3) ku = edges_gp(i ,j,k,idir)%u(3) #endif /* NDIMS == 3 */ case(2) jc = pmeta%pos(2) + 1 #if NDIMS == 2 il = edges_gp(i,j ,idir)%l(1) iu = edges_gp(i,j ,idir)%u(1) jl = edges_gp(i,jc ,idir)%l(2) ju = edges_gp(i,jc ,idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_gp(i,j ,k,idir)%l(1) iu = edges_gp(i,j ,k,idir)%u(1) jl = edges_gp(i,jc,k,idir)%l(2) ju = edges_gp(i,jc,k,idir)%u(2) kl = edges_gp(i,j ,k,idir)%l(3) ku = edges_gp(i,j ,k,idir)%u(3) case(3) kc = pmeta%pos(3) + 1 il = edges_gp(i,j,k ,idir)%l(1) iu = edges_gp(i,j,k ,idir)%u(1) jl = edges_gp(i,j,k ,idir)%l(2) ju = edges_gp(i,j,k ,idir)%u(2) kl = edges_gp(i,j,kc,idir)%l(3) ku = edges_gp(i,j,kc,idir)%u(3) #endif /* NDIMS == 3 */ end select ! extract the corresponding edge region from the neighbor and insert it in ! the current data block ! #if NDIMS == 2 call block_edge_prolong(idir, (/ ic, jc, kc /) & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_edge_prolong(idir, (/ ic, jc, kc /) & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku)) #endif /* NDIMS == 3 */ #ifdef MPI end if ! pneigh on the current process else ! block and neighbor belong to different processes ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, idir, (/ i, j, k /)) end if ! block and neighbor belong to different processes #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at lower level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides #if NDIMS == 3 end do ! k = 1, nsides #endif /* NDIMS == 3 */ ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member) ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) end if ! get the number of blocks to exchange ! scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) ! process only pairs which have anything to exchange ! if ((scount + rcount) > 0) then ! allocate buffers for variable exchange ! select case(idir) #if NDIMS == 2 case(1) allocate(sbuf(scount,nv,ih,ng, 1)) allocate(rbuf(rcount,nv,ih,ng, 1)) case(2) allocate(sbuf(scount,nv,ng,jh, 1)) allocate(rbuf(rcount,nv,ng,jh, 1)) #endif /* NDIMS == 2 */ #if NDIMS == 3 case(1) allocate(sbuf(scount,nv,ih,ng,ng)) allocate(rbuf(rcount,nv,ih,ng,ng)) case(2) allocate(sbuf(scount,nv,ng,jh,ng)) allocate(rbuf(rcount,nv,ng,jh,ng)) case(3) allocate(sbuf(scount,nv,ng,ng,kh)) allocate(rbuf(rcount,nv,ng,ng,kh)) #endif /* NDIMS == 3 */ end select !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pmeta and pneigh to the associated blocks ! pmeta => pinfo%meta pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! extract the corresponding edge region from the neighbor and insert it ! to the buffer ! select case(idir) case(1) i = pmeta%pos(1) + 1 #if NDIMS == 2 call block_edge_prolong(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ih,1:ng, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_edge_prolong(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ih,1:ng,1:ng)) #endif /* NDIMS == 3 */ case(2) j = pmeta%pos(2) + 1 #if NDIMS == 2 call block_edge_prolong(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:jh, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_edge_prolong(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:jh,1:ng)) #endif /* NDIMS == 3 */ #if NDIMS == 3 case(3) k = pmeta%pos(3) + 1 call block_edge_prolong(idir, (/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:ng,1:kh)) #endif /* NDIMS == 3 */ end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, sbuf, rbuf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! update the corresponding corner region of the current block ! select case(idir) case(1) ic = pmeta%pos(1) + 1 #if NDIMS == 2 il = edges_gp(ic,j ,idir)%l(1) iu = edges_gp(ic,j ,idir)%u(1) jl = edges_gp(i ,j ,idir)%l(2) ju = edges_gp(i ,j ,idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_gp(ic,j,k,idir)%l(1) iu = edges_gp(ic,j,k,idir)%u(1) jl = edges_gp(i ,j,k,idir)%l(2) ju = edges_gp(i ,j,k,idir)%u(2) kl = edges_gp(i ,j,k,idir)%l(3) ku = edges_gp(i ,j,k,idir)%u(3) #endif /* NDIMS == 3 */ #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = & rbuf(l,1:nv,1:ih,1:ng, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:ih,1:ng,1:ng) #endif /* NDIMS == 3 */ case(2) jc = pmeta%pos(2) + 1 #if NDIMS == 2 il = edges_gp(i,j ,idir)%l(1) iu = edges_gp(i,j ,idir)%u(1) jl = edges_gp(i,jc ,idir)%l(2) ju = edges_gp(i,jc ,idir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_gp(i,j ,k,idir)%l(1) iu = edges_gp(i,j ,k,idir)%u(1) jl = edges_gp(i,jc,k,idir)%l(2) ju = edges_gp(i,jc,k,idir)%u(2) kl = edges_gp(i,j ,k,idir)%l(3) ku = edges_gp(i,j ,k,idir)%u(3) #endif /* NDIMS == 3 */ #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = & rbuf(l,1:nv,1:ng,1:jh, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:ng,1:jh,1:ng) case(3) kc = pmeta%pos(3) + 1 il = edges_gp(i,j,k ,idir)%l(1) iu = edges_gp(i,j,k ,idir)%u(1) jl = edges_gp(i,j,k ,idir)%l(2) ju = edges_gp(i,j,k ,idir)%u(2) kl = edges_gp(i,j,kc,idir)%l(3) ku = edges_gp(i,j,kc,idir)%u(3) pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = & rbuf(l,1:nv,1:ng,1:ng,1:kh) #endif /* NDIMS == 3 */ end select ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(sbuf, rbuf) end if ! (scount + rcount) > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_edge_prolong ! !=============================================================================== ! ! DOMAIN CORNER BOUNDARY UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BOUNDARIES_CORNER_COPY: ! --------------------------------- ! ! Subroutine updates the corner boundaries from blocks on the same level. ! ! !=============================================================================== ! subroutine boundaries_corner_copy() ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : ng => nghosts use coordinates, only : corners_gc, corners_dc use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i, il, iu, is, it integer :: j, jl, ju, js, jt integer :: k = 1 #if NDIMS == 3 integer :: kl, ku, ks, kt #endif /* NDIMS == 3 */ #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: ecount = 0 integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: buf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! #if NDIMS == 3 do k = 1, nsides #endif /* NDIMS == 3 */ do j = 1, nsides do i = 1, nsides ! assign pneigh to the current neighbor ! #if NDIMS == 2 pneigh => pmeta%corners(i,j)%ptr #endif /* NDIMS == 2 */ #if NDIMS == 3 pneigh => pmeta%corners(i,j,k)%ptr #endif /* NDIMS == 3 */ ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor is at the same level ! if (pneigh%level == pmeta%level) then ! skip if the block and its neighbor are not marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare region indices of the block and its neighbor for the corner boundary ! update #if NDIMS == 2 il = corners_gc(i,j )%l(1) jl = corners_gc(i,j )%l(2) iu = corners_gc(i,j )%u(1) ju = corners_gc(i,j )%u(2) is = corners_dc(i,j )%l(1) js = corners_dc(i,j )%l(2) it = corners_dc(i,j )%u(1) jt = corners_dc(i,j )%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = corners_gc(i,j,k)%l(1) jl = corners_gc(i,j,k)%l(2) kl = corners_gc(i,j,k)%l(3) iu = corners_gc(i,j,k)%u(1) ju = corners_gc(i,j,k)%u(2) ku = corners_gc(i,j,k)%u(3) is = corners_dc(i,j,k)%l(1) js = corners_dc(i,j,k)%l(2) ks = corners_dc(i,j,k)%l(3) it = corners_dc(i,j,k)%u(1) jt = corners_dc(i,j,k)%u(2) kt = corners_dc(i,j,k)%u(3) #endif /* NDIMS == 3 */ ! copy the corresponding corner region from the neighbor to the current ! data block ! #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) & = pneigh%data%q(1:nv,is:it,js:jt, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) & = pneigh%data%q(1:nv,is:it,js:jt,ks:kt) #endif /* NDIMS == 3 */ #ifdef MPI end if ! pneigh on the current process else ! block and neighbor belong to different processes ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, -1, (/ i, j, k /)) end if ! block and neighbor belong to different processes #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at the same level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides #if NDIMS == 3 end do ! k = 1, nsides #endif /* NDIMS == 3 */ ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member), ! and the number of blocks to exchange ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) ecount = bcount(sproc,rproc) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) ecount = bcount(rproc,sproc) end if ! process only pairs which have anything to exchange ! if (ecount > 0) then ! allocate buffers for variable exchange ! #if NDIMS == 2 allocate(buf(ecount,nv,ng,ng, 1)) #endif /* NDIMS == 2 */ #if NDIMS == 3 allocate(buf(ecount,nv,ng,ng,ng)) #endif /* NDIMS == 3 */ !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pneigh to the associated neighbor block ! pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! prepare the corner region indices for the neighbor ! #if NDIMS == 2 is = corners_dc(i,j )%l(1) js = corners_dc(i,j )%l(2) it = corners_dc(i,j )%u(1) jt = corners_dc(i,j )%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 is = corners_dc(i,j,k)%l(1) js = corners_dc(i,j,k)%l(2) ks = corners_dc(i,j,k)%l(3) it = corners_dc(i,j,k)%u(1) jt = corners_dc(i,j,k)%u(2) kt = corners_dc(i,j,k)%u(3) #endif /* NDIMS == 3 */ ! copy the corresponding corner region from the neighbor to the buffer ! #if NDIMS == 2 buf(l,1:nv,1:ng,1:ng, : ) = pneigh%data%q(1:nv,is:it,js:jt, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 buf(l,1:nv,1:ng,1:ng,1:ng) = pneigh%data%q(1:nv,is:it,js:jt,ks:kt) #endif /* NDIMS == 3 */ ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, buf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! prepare the corner region indices for the block ! #if NDIMS == 2 il = corners_gc(i,j )%l(1) jl = corners_gc(i,j )%l(2) iu = corners_gc(i,j )%u(1) ju = corners_gc(i,j )%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = corners_gc(i,j,k)%l(1) jl = corners_gc(i,j,k)%l(2) kl = corners_gc(i,j,k)%l(3) iu = corners_gc(i,j,k)%u(1) ju = corners_gc(i,j,k)%u(2) ku = corners_gc(i,j,k)%u(3) #endif /* NDIMS == 3 */ ! update the corresponding corner region of the current block ! #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = buf(l,1:nv,1:ng,1:ng, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = buf(l,1:nv,1:ng,1:ng,1:ng) #endif /* NDIMS == 3 */ ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(buf) end if ! ecount > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_corner_copy ! !=============================================================================== ! ! subroutine BOUNDARIES_CORNER_RESTRICT: ! ------------------------------------- ! ! Subroutine updates the corner boundaries from blocks on higher level. ! ! !=============================================================================== ! subroutine boundaries_corner_restrict() ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : ng => nghosts use coordinates, only : corners_gr use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i, il, iu integer :: j, jl, ju integer :: k = 1 #if NDIMS == 3 integer :: kl, ku #endif /* NDIMS == 3 */ #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: scount, rcount integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! #if NDIMS == 3 do k = 1, nsides #endif /* NDIMS == 3 */ do j = 1, nsides do i = 1, nsides ! assign pneigh to the current neighbor ! #if NDIMS == 2 pneigh => pmeta%corners(i,j)%ptr #endif /* NDIMS == 2 */ #if NDIMS == 3 pneigh => pmeta%corners(i,j,k)%ptr #endif /* NDIMS == 3 */ ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor is at higher level ! if (pneigh%level > pmeta%level) then ! skip if the block and its neighbor are not marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare the region indices for corner boundary update ! #if NDIMS == 2 il = corners_gr(i,j )%l(1) jl = corners_gr(i,j )%l(2) iu = corners_gr(i,j )%u(1) ju = corners_gr(i,j )%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = corners_gr(i,j,k)%l(1) jl = corners_gr(i,j,k)%l(2) kl = corners_gr(i,j,k)%l(3) iu = corners_gr(i,j,k)%u(1) ju = corners_gr(i,j,k)%u(2) ku = corners_gr(i,j,k)%u(3) #endif /* NDIMS == 3 */ ! extract and restrict the corresponding corner region from the neighbor and ! insert it in the current data block ! #if NDIMS == 2 call block_corner_restrict((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_corner_restrict((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku)) #endif /* NDIMS == 3 */ #ifdef MPI end if ! block on the current processor else ! block and neighbor on different processors ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, -1, (/ i, j, k /)) end if ! block and neighbor on different processors #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at higher level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides #if NDIMS == 3 end do ! k = 1, nsides #endif /* NDIMS == 3 */ ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member) ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) end if ! get the number of blocks to exchange ! scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) ! process only pairs which have anything to exchange ! if ((scount + rcount) > 0) then ! allocate buffers for variable exchange ! #if NDIMS == 2 allocate(sbuf(scount,nv,ng,ng, 1)) allocate(rbuf(rcount,nv,ng,ng, 1)) #endif /* NDIMS == 2 */ #if NDIMS == 3 allocate(sbuf(scount,nv,ng,ng,ng)) allocate(rbuf(rcount,nv,ng,ng,ng)) #endif /* NDIMS == 3 */ !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pneigh to the associated neighbor block ! pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! restrict and extract the corresponding corner region from the neighbor and ! insert it to the buffer ! #if NDIMS == 2 call block_corner_restrict((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:ng,1: )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_corner_restrict((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:ng,1:ng)) #endif /* NDIMS == 3 */ ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, sbuf, rbuf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! prepare the region indices for corner boundary update ! #if NDIMS == 2 il = corners_gr(i,j )%l(1) jl = corners_gr(i,j )%l(2) iu = corners_gr(i,j )%u(1) ju = corners_gr(i,j )%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = corners_gr(i,j,k)%l(1) jl = corners_gr(i,j,k)%l(2) kl = corners_gr(i,j,k)%l(3) iu = corners_gr(i,j,k)%u(1) ju = corners_gr(i,j,k)%u(2) ku = corners_gr(i,j,k)%u(3) #endif /* NDIMS == 3 */ ! update the corresponding corner region of the current block ! #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = rbuf(l,1:nv,1:ng,1:ng, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = rbuf(l,1:nv,1:ng,1:ng,1:ng) #endif /* NDIMS == 3 */ ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(sbuf, rbuf) end if ! (scount + rcount) > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_corner_restrict ! !=============================================================================== ! ! subroutine BOUNDARIES_CORNER_PROLONG: ! ------------------------------------ ! ! Subroutine updates the corner boundaries from blocks on lower level. ! ! !=============================================================================== ! subroutine boundaries_corner_prolong() ! import external procedures and variables ! use blocks , only : nsides use blocks , only : block_meta, block_data, block_leaf use blocks , only : list_leaf use blocks , only : block_info, pointer_info use coordinates, only : ng => nghosts use coordinates, only : corners_gp use equations , only : nv #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ ! local variables are not implicit by default ! implicit none ! local pointers ! type(block_meta), pointer :: pmeta, pneigh type(block_leaf), pointer :: pleaf #ifdef MPI type(block_info), pointer :: pinfo #endif /* MPI */ ! local variables ! integer :: i, il, iu integer :: j, jl, ju integer :: k = 1 #if NDIMS == 3 integer :: kl, ku #endif /* NDIMS == 3 */ #ifdef MPI integer :: sproc = 0, rproc = 0 integer :: scount, rcount integer :: l, p ! local arrays ! real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ ! !------------------------------------------------------------------------------- ! #ifdef MPI ! prepare the block exchange structures ! call prepare_exchange_array() #endif /* MPI */ ! update boundaries between blocks on the same process ! ! associate pleaf with the first block on the leaf list ! pleaf => list_leaf ! scan all leaf meta blocks in the list ! do while(associated(pleaf)) ! get the associated meta block ! pmeta => pleaf%meta ! scan over all block corners ! #if NDIMS == 3 do k = 1, nsides #endif /* NDIMS == 3 */ do j = 1, nsides do i = 1, nsides ! assign pneigh to the current neighbor ! #if NDIMS == 2 pneigh => pmeta%corners(i,j)%ptr #endif /* NDIMS == 2 */ #if NDIMS == 3 pneigh => pmeta%corners(i,j,k)%ptr #endif /* NDIMS == 3 */ ! check if the neighbor is associated ! if (associated(pneigh)) then ! check if the neighbor lays at lower level ! if (pneigh%level < pmeta%level) then ! skip if the block and its neighbor are not marked for update ! if (pmeta%update .or. pneigh%update) then #ifdef MPI ! check if the block and its neighbor belong to the same process ! if (pmeta%process == pneigh%process) then ! check if the neighbor belongs to the current process ! if (pneigh%process == nproc) then #endif /* MPI */ ! prepare the region indices for corner boundary update ! #if NDIMS == 2 il = corners_gp(i,j )%l(1) jl = corners_gp(i,j )%l(2) iu = corners_gp(i,j )%u(1) ju = corners_gp(i,j )%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = corners_gp(i,j,k)%l(1) jl = corners_gp(i,j,k)%l(2) kl = corners_gp(i,j,k)%l(3) iu = corners_gp(i,j,k)%u(1) ju = corners_gp(i,j,k)%u(2) ku = corners_gp(i,j,k)%u(3) #endif /* NDIMS == 3 */ ! restrict and extract the corresponding corner region from the neighbor and ! insert it in the current data block ! #if NDIMS == 2 call block_corner_prolong((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_corner_prolong((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku)) #endif /* NDIMS == 3 */ #ifdef MPI end if ! block on the current processor else ! block and neighbor on different processors ! append the block to the exchange list ! call append_exchange_block(pmeta, pneigh, -1, (/ i, j, k /)) end if ! block and neighbor on different processors #endif /* MPI */ end if ! pmeta and pneigh marked for update end if ! neighbor at lower level end if ! neighbor associated end do ! i = 1, nsides end do ! j = 1, nsides #if NDIMS == 3 end do ! k = 1, nsides #endif /* NDIMS == 3 */ ! associate pleaf with the next leaf on the list ! pleaf => pleaf%next end do ! over leaf blocks #ifdef MPI !! 3. UPDATE VARIABLE BOUNDARIES BETWEEN BLOCKS BELONGING TO DIFFERENT PROCESSES !! ! iterate over all process pairs ! do p = 1, npairs ! process only pairs related to this process ! if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then ! get sending and receiving process identifiers (depending on pair member) ! if (pairs(p,1) == nproc) then sproc = pairs(p,1) rproc = pairs(p,2) end if if (pairs(p,2) == nproc) then sproc = pairs(p,2) rproc = pairs(p,1) end if ! get the number of blocks to exchange ! scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) ! process only pairs which have anything to exchange ! if ((scount + rcount) > 0) then ! allocate buffers for variable exchange ! #if NDIMS == 2 allocate(sbuf(scount,nv,ng,ng, 1)) allocate(rbuf(rcount,nv,ng,ng, 1)) #endif /* NDIMS == 2 */ #if NDIMS == 3 allocate(sbuf(scount,nv,ng,ng,ng)) allocate(rbuf(rcount,nv,ng,ng,ng)) #endif /* NDIMS == 3 */ !! PREPARE BLOCKS FOR SENDING !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(sproc,rproc)%ptr ! scan over all blocks on the block exchange list ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign pneigh to the associated neighbor block ! pneigh => pinfo%neigh ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! prolong the corresponding corner region from the neighbor and insert it in ! the buffer ! #if NDIMS == 2 call block_corner_prolong((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:ng, : )) #endif /* NDIMS == 2 */ #if NDIMS == 3 call block_corner_prolong((/ i, j, k /) & , pneigh%data%q(1:nv, : , : , : ) & , sbuf(l,1:nv,1:ng,1:ng,1:ng)) #endif /* NDIMS == 3 */ ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list !! SEND PREPARED BLOCKS AND RECEIVE NEW ONES !! ! exchange data ! call exchange_arrays(rproc, p, sbuf, rbuf) !! PROCESS RECEIVED BLOCKS !! ! reset the block counter ! l = 0 ! associate the pointer with the first block in the exchange list ! pinfo => barray(rproc,sproc)%ptr ! iterate over all received blocks and update boundaries of the corresponding ! data blocks ! do while(associated(pinfo)) ! increase the block counter ! l = l + 1 ! assign a pointer to the associated data block ! pmeta => pinfo%meta ! get the corner coordinates ! i = pinfo%corner(1) j = pinfo%corner(2) #if NDIMS == 3 k = pinfo%corner(3) #endif /* NDIMS == 3 */ ! prepare the region indices for corner boundary update ! #if NDIMS == 2 il = corners_gp(i,j )%l(1) jl = corners_gp(i,j )%l(2) iu = corners_gp(i,j )%u(1) ju = corners_gp(i,j )%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = corners_gp(i,j,k)%l(1) jl = corners_gp(i,j,k)%l(2) kl = corners_gp(i,j,k)%l(3) iu = corners_gp(i,j,k)%u(1) ju = corners_gp(i,j,k)%u(2) ku = corners_gp(i,j,k)%u(3) #endif /* NDIMS == 3 */ ! update the corresponding corner region of the current block ! #if NDIMS == 2 pmeta%data%q(1:nv,il:iu,jl:ju, : ) = rbuf(l,1:nv,1:ng,1:ng, : ) #endif /* NDIMS == 2 */ #if NDIMS == 3 pmeta%data%q(1:nv,il:iu,jl:ju,kl:ku) = rbuf(l,1:nv,1:ng,1:ng,1:ng) #endif /* NDIMS == 3 */ ! associate the pointer with the next block ! pinfo => pinfo%prev end do ! %ptr block list ! deallocate data buffer ! deallocate(sbuf, rbuf) end if ! (scount + rcount) > 0 end if ! pairs(p,1) == nproc || pairs(p,2) == nproc end do ! p = 1, npairs ! release the memory used by the array of exchange block lists ! call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_corner_prolong ! !=============================================================================== ! ! BLOCK SPECIFIC BOUNDARY SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_BOUNDARY_SPECIFIC: ! ---------------------------------- ! ! Subroutine applies specific boundary conditions to the pointed data block. ! ! Arguments: ! ! nc - the direction; ! side - the side of the boundary; ! t, dt - time and time increment; ! x, y, z - the block coordinates; ! qn - the variable array; ! !=============================================================================== ! subroutine block_boundary_specific(nc, side, t, dt, x, y, z, qn) use coordinates , only : nn => bcells, nh => bcells_half, ng => nghosts use coordinates , only : nb, ne, nbl, neu use equations , only : nv use equations , only : idn, ipr, ivx, ivy, ibx, iby #if NDIMS == 3 use equations , only : ivz, ibz #endif /* NDIMS == 3 */ use equations , only : csnd2 use gravity , only : gravitational_acceleration use helpers , only : print_message use user_problem, only : user_boundary_x, user_boundary_y, & user_boundary_z implicit none integer , intent(in) :: nc integer , dimension(3) , intent(in) :: side real(kind=8) , intent(in) :: t, dt real(kind=8), dimension(:) , intent(inout) :: x real(kind=8), dimension(:) , intent(inout) :: y real(kind=8), dimension(:) , intent(inout) :: z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn integer :: i, il, iu, is, it, im1, ip1 integer :: j, jl, ju, js, jt, jm1, jp1 integer :: k, kl, ku #if NDIMS == 3 integer :: ks, kt, km1, kp1 real(kind=8) :: dz, dzh, zi #endif /* NDIMS == 3 */ real(kind=8) :: dx, dxh, xi, dy, dyh, yi real(kind=8), dimension(3) :: ga character(len=*), parameter :: loc = 'BOUNDARIES::block_boundary_specific()' !------------------------------------------------------------------------------- ! ! apply specific boundaries depending on the direction ! select case(nc) case(1) ! prepare indices for the boundaries ! if (side(2) == 1) then jl = 1 ju = nh - 1 else jl = nh ju = nn end if #if NDIMS == 3 if (side(3) == 1) then kl = 1 ku = nh - 1 else kl = nh ku = nn end if #else /* NDIMS == 3 */ kl = 1 ku = 1 #endif /* NDIMS == 3 */ ! apply selected boundary condition ! select case(bnd_type(nc,side(1))) ! "open" boundary conditions ! case(bnd_open) if (side(1) == 1) then do i = nbl, 1, -1 qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,nb,jl:ju,kl:ku) end do else do i = neu, nn qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,ne,jl:ju,kl:ku) end do end if ! "outflow" boundary conditions ! case(bnd_outflow) if (side(1) == 1) then do i = nbl, 1, -1 qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,nb,jl:ju,kl:ku) qn(ivx ,i,jl:ju,kl:ku) = min(0.0d+00, qn(ivx,nb,jl:ju,kl:ku)) end do ! i = nbl, 1, -1 else do i = neu, nn qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,ne,jl:ju,kl:ku) qn(ivx ,i,jl:ju,kl:ku) = max(0.0d+00, qn(ivx,ne,jl:ju,kl:ku)) end do ! i = neu, nn end if ! "reflective" boundary conditions ! case(bnd_reflective) if (side(1) == 1) then do i = 1, ng it = nb - i is = nbl + i qn(1:nv,it,jl:ju,kl:ku) = qn(1:nv,is,jl:ju,kl:ku) qn(ivx ,it,jl:ju,kl:ku) = - qn(ivx ,is,jl:ju,kl:ku) if (ibx > 0) then qn(ibx ,it,jl:ju,kl:ku) = - qn(ibx ,is,jl:ju,kl:ku) end if end do else do i = 1, ng it = ne + i is = neu - i qn(1:nv,it,jl:ju,kl:ku) = qn(1:nv,is,jl:ju,kl:ku) qn(ivx ,it,jl:ju,kl:ku) = - qn(ivx ,is,jl:ju,kl:ku) if (ibx > 0) then qn(ibx ,it,jl:ju,kl:ku) = - qn(ibx ,is,jl:ju,kl:ku) end if end do end if ! "gravity" or "hydrostatic" boundary conditions ! case(bnd_gravity) dx = x(nb) - x(nbl) dxh = 0.5d+00 * dx if (ipr > 0) then if (side(1) == 1) then do i = nbl, 1, -1 ip1 = i + 1 xi = x(i) + dxh do k = kl, ku do j = jl, ju qn(1:nv,i,j,k) = qn(1:nv,nb,j,k) call gravitational_acceleration(t, dt, xi, y(j), z(k), ga(:)) qn(ipr,i,j,k) = qn(ipr,ip1,j,k) & - (qn(idn,ip1,j,k) + qn(idn,i,j,k)) * ga(1) * dxh end do end do end do else do i = neu, nn im1 = i - 1 xi = x(i) - dxh do k = kl, ku do j = jl, ju qn(1:nv,i,j,k) = qn(1:nv,ne,j,k) call gravitational_acceleration(t, dt, xi, y(j), z(k), ga(:)) qn(ipr,i,j,k) = qn(ipr,im1,j,k) & + (qn(idn,im1,j,k) + qn(idn,i,j,k)) * ga(1) * dxh end do end do end do end if else if (side(1) == 1) then do i = nbl, 1, -1 ip1 = i + 1 xi = x(i) + dxh do k = kl, ku do j = jl, ju qn(1:nv,i,j,k) = qn(1:nv,nb,j,k) call gravitational_acceleration(t, dt, xi, y(j), z(k), ga(:)) qn(idn,i,j,k) = qn(idn,ip1,j,k) * exp(- ga(1) * dx / csnd2) end do end do end do else do i = neu, nn im1 = i - 1 xi = x(i) - dxh do k = kl, ku do j = jl, ju qn(1:nv,i,j,k) = qn(1:nv,ne,j,k) call gravitational_acceleration(t, dt, xi, y(j), z(k), ga(:)) qn(idn,i,j,k) = qn(idn,im1,j,k) * exp( ga(1) * dx / csnd2) end do end do end do end if end if ! user specific boundary conditions ! case(bnd_user) call user_boundary_x(side(1), jl, ju, kl, ku & , t, dt, x(:), y(:), z(:), qn(:,:,:,:)) ! wrong boundary conditions ! case default if (side(1) == 1) then call print_message(loc, "Wrong left X boundary type!") else call print_message(loc, "Wrong right X boundary type!") end if end select case(2) ! prepare indices for the boundaries ! if (side(1) == 1) then il = 1 iu = nh - 1 else il = nh iu = nn end if #if NDIMS == 3 if (side(3) == 1) then kl = 1 ku = nh - 1 else kl = nh ku = nn end if #else /* NDIMS == 3 */ kl = 1 ku = 1 #endif /* NDIMS == 3 */ ! apply selected boundary condition ! select case(bnd_type(nc,side(2))) ! "open" boundary conditions ! case(bnd_open) if (side(2) == 1) then do j = nbl, 1, -1 qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,nb,kl:ku) end do else do j = neu, nn qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,ne,kl:ku) end do end if ! "outflow" boundary conditions ! case(bnd_outflow) if (side(2) == 1) then do j = nbl, 1, -1 qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,nb,kl:ku) qn(ivy ,il:iu,j,kl:ku) = min(0.0d+00, qn(ivy,il:iu,nb,kl:ku)) end do ! j = nbl, 1, -1 else do j = neu, nn qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,ne,kl:ku) qn(ivy ,il:iu,j,kl:ku) = max(0.0d+00, qn(ivy,il:iu,ne,kl:ku)) end do ! j = neu, ng end if ! "reflective" boundary conditions ! case(bnd_reflective) if (side(2) == 1) then do j = 1, ng jt = nb - j js = nbl + j qn(1:nv,il:iu,jt,kl:ku) = qn(1:nv,il:iu,js,kl:ku) qn(ivy ,il:iu,jt,kl:ku) = - qn(ivy ,il:iu,js,kl:ku) if (iby > 0) then qn(iby ,il:iu,jt,kl:ku) = - qn(iby ,il:iu,js,kl:ku) end if end do else do j = 1, ng jt = ne + j js = neu - j qn(1:nv,il:iu,jt,kl:ku) = qn(1:nv,il:iu,js,kl:ku) qn(ivy ,il:iu,jt,kl:ku) = - qn(ivy ,il:iu,js,kl:ku) if (iby > 0) then qn(iby ,il:iu,jt,kl:ku) = - qn(iby ,il:iu,js,kl:ku) end if end do end if ! "gravity" or "hydrostatic" boundary conditions ! case(bnd_gravity) dy = y(nb) - y(nbl) dyh = 0.5d+00 * dy if (ipr > 0) then if (side(2) == 1) then do j = nbl, 1, -1 jp1 = j + 1 yi = y(j) + dyh do k = kl, ku do i = il, iu qn(1:nv,i,j,k) = qn(1:nv,i,nb,k) call gravitational_acceleration(t, dt, x(i), yi, z(k), ga(:)) qn(ipr,i,j,k) = qn(ipr,i,jp1,k) & - (qn(idn,i,jp1,k) + qn(idn,i,j,k)) * ga(2) * dyh end do end do end do else do j = neu, nn jm1 = j - 1 yi = y(j) - dyh do k = kl, ku do i = il, iu qn(1:nv,i,j,k) = qn(1:nv,i,ne,k) call gravitational_acceleration(t, dt, x(i), yi, z(k), ga(:)) qn(ipr,i,j,k) = qn(ipr,i,jm1,k) & + (qn(idn,i,jm1,k) + qn(idn,i,j,k)) * ga(2) * dyh end do end do end do end if else if (side(2) == 1) then do j = nbl, 1, -1 jp1 = j + 1 yi = y(j) + dyh do k = kl, ku do i = il, iu qn(1:nv,i,j,k) = qn(1:nv,i,nb,k) call gravitational_acceleration(t, dt, x(i), yi, z(k), ga(:)) qn(idn,i,j,k) = qn(idn,i,jp1,k) * exp(- ga(2) * dy / csnd2) end do end do end do else do j = neu, nn jm1 = j - 1 yi = y(j) - dyh do k = kl, ku do i = il, iu qn(1:nv,i,j,k) = qn(1:nv,i,ne,k) call gravitational_acceleration(t, dt, x(i), yi, z(k), ga(:)) qn(idn,i,j,k) = qn(idn,i,jm1,k) * exp( ga(2) * dy / csnd2) end do end do end do end if end if ! user specific boundary conditions ! case(bnd_user) call user_boundary_y(side(2), il, iu, kl, ku & , t, dt, x(:), y(:), z(:), qn(:,:,:,:)) ! wrong boundary conditions ! case default if (side(2) == 1) then call print_message(loc, "Wrong left Y boundary type!") else call print_message(loc, "Wrong right Y boundary type!") end if end select #if NDIMS == 3 case(3) ! prepare indices for the boundaries ! if (side(1) == 1) then il = 1 iu = nh - 1 else il = nh iu = nn end if if (side(2) == 1) then jl = 1 ju = nh - 1 else jl = nh ju = nn end if ! apply selected boundary condition ! select case(bnd_type(nc,side(3))) ! "open" boundary conditions ! case(bnd_open) if (side(3) == 1) then do k = nbl, 1, -1 qn(1:nv,il:iu,jl:ju,k) = qn(1:nv,il:iu,jl:ju,nb) end do else do k = neu, nn qn(1:nv,il:iu,jl:ju,k) = qn(1:nv,il:iu,jl:ju,ne) end do end if ! "outflow" boundary conditions ! case(bnd_outflow) if (side(3) == 1) then do k = nbl, 1, -1 qn(1:nv,il:iu,jl:ju,k) = qn(1:nv,il:iu,jl:ju,nb) qn(ivz ,il:iu,jl:ju,k) = min(0.0d+00, qn(ivz,il:iu,jl:ju,nb)) end do ! k = nbl, 1, -1 else do k = neu, nn qn(1:nv,il:iu,jl:ju,k) = qn(1:nv,il:iu,jl:ju,ne) qn(ivz ,il:iu,jl:ju,k) = max(0.0d+00, qn(ivz,il:iu,jl:ju,ne)) end do ! k = neu, nn end if ! "reflective" boundary conditions ! case(bnd_reflective) if (side(3) == 1) then do k = 1, ng kt = nb - k ks = nbl + k qn(1:nv,il:iu,jl:ju,kt) = qn(1:nv,il:iu,jl:ju,ks) qn(ivz ,il:iu,jl:ju,kt) = - qn(ivz ,il:iu,jl:ju,ks) if (ibz > 0) then qn(ibz ,il:iu,jl:ju,kt) = - qn(ibz ,il:iu,jl:ju,ks) end if end do else do k = 1, ng kt = ne + k ks = neu - k qn(1:nv,il:iu,jl:ju,kt) = qn(1:nv,il:iu,jl:ju,ks) qn(ivz ,il:iu,jl:ju,kt) = - qn(ivz ,il:iu,jl:ju,ks) if (ibz > 0) then qn(ibz ,il:iu,jl:ju,kt) = - qn(ibz ,il:iu,jl:ju,ks) end if end do end if ! "gravity" or "hydrostatic" boundary conditions ! case(bnd_gravity) dz = z(nb) - z(nbl) dzh = 0.5d+00 * dz if (ipr > 0) then if (side(3) == 1) then do k = nbl, 1, -1 kp1 = k + 1 zi = z(k) + dzh do j = jl, ju do i = il, iu qn(1:nv,i,j,k) = qn(1:nv,i,j,nb) call gravitational_acceleration(t, dt, x(i), y(j), zi, ga(:)) qn(ipr,i,j,k) = qn(ipr,i,j,kp1) & - (qn(idn,i,j,kp1) + qn(idn,i,j,k)) * ga(3) * dzh end do end do end do else do k = neu, nn km1 = k - 1 zi = z(k) - dzh do j = jl, ju do i = il, iu qn(1:nv,i,j,k) = qn(1:nv,i,j,ne) call gravitational_acceleration(t, dt, x(i), y(j), zi, ga(:)) qn(ipr,i,j,k) = qn(ipr,i,j,km1) & + (qn(idn,i,j,km1) + qn(idn,i,j,k)) * ga(3) * dzh end do end do end do end if else if (side(3) == 1) then do k = nbl, 1, -1 kp1 = k + 1 zi = z(k) + dzh do j = jl, ju do i = il, iu qn(1:nv,i,j,k) = qn(1:nv,i,j,nb) call gravitational_acceleration(t, dt, x(i), y(j), zi, ga(:)) qn(idn,i,j,k) = qn(idn,i,j,kp1) * exp(- ga(3) * dz / csnd2) end do end do end do else do k = neu, nn km1 = k - 1 zi = z(k) - dzh do j = jl, ju do i = il, iu qn(1:nv,i,j,k) = qn(1:nv,i,j,ne) call gravitational_acceleration(t, dt, x(i), y(j), zi, ga(:)) qn(idn,i,j,k) = qn(idn,i,j,km1) * exp( ga(3) * dz / csnd2) end do end do end do end if end if ! user specific boundary conditions ! case(bnd_user) call user_boundary_z(side(3), il, iu, jl, ju & , t, dt, x(:), y(:), z(:), qn(:,:,:,:)) ! wrong boundary conditions ! case default if (side(3) == 1) then call print_message(loc, "Wrong left Z boundary type!") else call print_message(loc, "Wrong right Z boundary type!") end if end select #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine block_boundary_specific #if NDIMS == 3 ! !=============================================================================== ! ! BLOCK FACE UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_FACE_RESTRICT: ! ------------------------------ ! ! Subroutine returns the face boundary region restricted from the provided ! input variable array. ! ! Arguments: ! ! nc - the face direction; ! ic, jc, kc - the corner position; ! qn - the input neighbor variable array; ! qb - the output face boundary array; ! !=============================================================================== ! subroutine block_face_restrict(nc, ic, jc, kc, qn, qb) ! import external procedures and variables ! use coordinates, only : nh => ncells_half, ng => nghosts use coordinates, only : faces_dr use equations , only : nv ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer , intent(in) :: nc, ic, jc, kc real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(out) :: qb ! local variables ! integer :: il, jl, kl integer :: ip, jp, kp integer :: iu, ju, ku ! !------------------------------------------------------------------------------- ! ! prepare indices for the face region ! il = faces_dr(ic,jc,kc,nc)%l(1) jl = faces_dr(ic,jc,kc,nc)%l(2) kl = faces_dr(ic,jc,kc,nc)%l(3) ip = il + 1 jp = jl + 1 kp = kl + 1 iu = faces_dr(ic,jc,kc,nc)%u(1) ju = faces_dr(ic,jc,kc,nc)%u(2) ku = faces_dr(ic,jc,kc,nc)%u(3) ! process depending on the direction ! select case(nc) case(1) ! restrict the face region to the output array ! qb(1:nv,1:ng,1:nh,1:nh) = & 1.25d-01 * (((qn(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + ((qn(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + (qn(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) case(2) ! restrict the face region to the output array ! qb(1:nv,1:nh,1:ng,1:nh) = & 1.25d-01 * (((qn(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + ((qn(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + (qn(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) case(3) ! restrict the face region to the output array ! qb(1:nv,1:nh,1:nh,1:ng) = & 1.25d-01 * (((qn(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + ((qn(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + (qn(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) end select !------------------------------------------------------------------------------- ! end subroutine block_face_restrict ! !=============================================================================== ! ! subroutine BLOCK_FACE_PROLONG: ! ----------------------------- ! ! Subroutine returns the face boundary region prolongated from the provided ! input variable array. ! ! Arguments: ! ! nc - the face direction; ! ic, jc, kc - the corner position; ! qn - the input neighbor variable array; ! qb - the output face boundary array; ! !=============================================================================== ! subroutine block_face_prolong(nc, ic, jc, kc, qn, qb) use coordinates , only : faces_dp use equations , only : nv, positive use helpers , only : print_message use interpolations, only : limiter_prol implicit none integer , intent(in) :: nc, ic, jc, kc real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(out) :: qb integer :: i, j, k, p integer :: il, jl, kl integer :: iu, ju, ku integer :: is, js, ks integer :: it, jt, kt integer :: im1, jm1, km1 integer :: ip1, jp1, kp1 real(kind=8) :: dql, dqr real(kind=8) :: dq1, dq2, dq3, dq4 real(kind=8), dimension(3) :: dq character(len=80) :: msg character(len=*), parameter :: loc = 'BOUNDARIES::block_face_prolong()' !------------------------------------------------------------------------------- ! il = faces_dp(ic,jc,kc,nc)%l(1) jl = faces_dp(ic,jc,kc,nc)%l(2) kl = faces_dp(ic,jc,kc,nc)%l(3) iu = faces_dp(ic,jc,kc,nc)%u(1) ju = faces_dp(ic,jc,kc,nc)%u(2) ku = faces_dp(ic,jc,kc,nc)%u(3) do k = kl, ku km1 = k - 1 kp1 = k + 1 ks = 2 * (k - kl) + 1 kt = ks + 1 do j = jl, ju jm1 = j - 1 jp1 = j + 1 js = 2 * (j - jl) + 1 jt = js + 1 do i = il, iu im1 = i - 1 ip1 = i + 1 is = 2 * (i - il) + 1 it = is + 1 do p = 1, nv dql = qn(p,i ,j,k) - qn(p,im1,j,k) dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) dq(1) = limiter_prol(2.5d-01, dql, dqr) dql = qn(p,i,j ,k) - qn(p,i,jm1,k) dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) dq(2) = limiter_prol(2.5d-01, dql, dqr) dql = qn(p,i,j,k ) - qn(p,i,j,km1) dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) dq(3) = limiter_prol(2.5d-01, dql, dqr) if (positive(p) .and. qn(p,i,j,k) < sum(abs(dq(1:NDIMS)))) then if (qn(p,i,j,k) > 0.0d+00) then do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) dq(:) = 0.5d+00 * dq(:) end do else write(msg,"(a,3i4,a)") & "Positive variable is not positive at (", i, j, k, " )" call print_message(loc, msg) dq(:) = 0.0d+00 end if end if dq1 = dq(1) + dq(2) + dq(3) dq2 = dq(1) - dq(2) - dq(3) dq3 = dq(1) - dq(2) + dq(3) dq4 = dq(1) + dq(2) - dq(3) qb(p,is,js,ks) = qn(p,i,j,k) - dq1 qb(p,it,js,ks) = qn(p,i,j,k) + dq2 qb(p,is,jt,ks) = qn(p,i,j,k) - dq3 qb(p,it,jt,ks) = qn(p,i,j,k) + dq4 qb(p,is,js,kt) = qn(p,i,j,k) - dq4 qb(p,it,js,kt) = qn(p,i,j,k) + dq3 qb(p,is,jt,kt) = qn(p,i,j,k) - dq2 qb(p,it,jt,kt) = qn(p,i,j,k) + dq1 end do ! p end do ! i end do ! j end do ! k !------------------------------------------------------------------------------- ! end subroutine block_face_prolong #endif /* NDIMS == 3 */ ! !=============================================================================== ! ! BLOCK EDGE UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_EDGE_RESTRICT: ! ------------------------------ ! ! Subroutine returns the edge boundary region by restricting the corresponding ! region from the provided input variable array. ! ! Arguments: ! ! dir - the edge direction; ! pos - the edge position; ! qn - the input neighbor variable array; ! qb - the output edge boundary array; ! !=============================================================================== ! subroutine block_edge_restrict(dir, pos, qn, qb) ! import external procedures and variables ! use coordinates, only : nh => ncells_half, ng => nghosts use coordinates, only : edges_dr use equations , only : nv ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer , intent(in) :: dir integer , dimension(3) , intent(in) :: pos real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(out) :: qb ! local variables ! integer :: il, ip, iu integer :: jl, jp, ju #if NDIMS == 3 integer :: kl, kp, ku #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! ! prepare indices for the edge region ! #if NDIMS == 2 il = edges_dr(pos(1),pos(2),dir)%l(1) jl = edges_dr(pos(1),pos(2),dir)%l(2) ip = il + 1 jp = jl + 1 iu = edges_dr(pos(1),pos(2),dir)%u(1) ju = edges_dr(pos(1),pos(2),dir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_dr(pos(1),pos(2),pos(3),dir)%l(1) jl = edges_dr(pos(1),pos(2),pos(3),dir)%l(2) kl = edges_dr(pos(1),pos(2),pos(3),dir)%l(3) ip = il + 1 jp = jl + 1 kp = kl + 1 iu = edges_dr(pos(1),pos(2),pos(3),dir)%u(1) ju = edges_dr(pos(1),pos(2),pos(3),dir)%u(2) ku = edges_dr(pos(1),pos(2),pos(3),dir)%u(3) #endif /* NDIMS == 3 */ ! process depending on the direction ! select case(dir) case(1) ! restrict the edge region to the output array ! #if NDIMS == 2 qb(1:nv,1:nh,1:ng, : ) = & 2.50d-01 * ((qn(1:nv,il:iu:2,jl:ju:2, : ) & + qn(1:nv,ip:iu:2,jp:ju:2, : )) & + (qn(1:nv,il:iu:2,jp:ju:2, : ) & + qn(1:nv,ip:iu:2,jl:ju:2, : ))) #endif /* NDIMS == 2 */ #if NDIMS == 3 qb(1:nv,1:nh,1:ng,1:ng) = & 1.25d-01 * (((qn(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + ((qn(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + (qn(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) #endif /* NDIMS == 3 */ case(2) ! restrict the edge region to the output array ! #if NDIMS == 2 qb(1:nv,1:ng,1:nh, : ) = & 2.50d-01 * ((qn(1:nv,il:iu:2,jl:ju:2, : ) & + qn(1:nv,ip:iu:2,jp:ju:2, : )) & + (qn(1:nv,il:iu:2,jp:ju:2, : ) & + qn(1:nv,ip:iu:2,jl:ju:2, : ))) #endif /* NDIMS == 2 */ #if NDIMS == 3 qb(1:nv,1:ng,1:nh,1:ng) = & 1.25d-01 * (((qn(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + ((qn(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + (qn(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) #endif /* NDIMS == 3 */ #if NDIMS == 3 case(3) ! restrict the edge region to the output array ! qb(1:nv,1:ng,1:ng,1:nh) = & 1.25d-01 * (((qn(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + ((qn(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + (qn(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine block_edge_restrict ! !=============================================================================== ! ! subroutine BLOCK_EDGE_PROLONG: ! ----------------------------- ! ! Subroutine returns the edge boundary region by prolongating ! the corresponding region from the provided input variable array. ! ! Arguments: ! ! dir - the edge direction; ! pos - the edge position; ! qn - the input neighbor variable array; ! qb - the output edge boundary array; ! !=============================================================================== ! subroutine block_edge_prolong(dir, pos, qn, qb) use coordinates , only : edges_dp use equations , only : nv, positive use helpers , only : print_message use interpolations, only : limiter_prol implicit none integer , intent(in) :: dir integer , dimension(3) , intent(in) :: pos real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(out) :: qb integer :: p integer :: i, il, iu, is, it, im1, ip1 integer :: j, jl, ju, js, jt, jm1, jp1 integer :: k, kt #if NDIMS == 3 integer :: kl, ku, ks, km1, kp1 #endif /* NDIMS == 3 */ real(kind=8) :: dql, dqr real(kind=8) :: dq1, dq2 #if NDIMS == 3 real(kind=8) :: dq3, dq4 #endif /* NDIMS == 3 */ real(kind=8), dimension(NDIMS) :: dq character(len=80) :: msg character(len=*), parameter :: loc = 'BOUNDARIES::block_edge_prolong()' ! !------------------------------------------------------------------------------- ! #if NDIMS == 2 il = edges_dp(pos(1),pos(2),dir)%l(1) jl = edges_dp(pos(1),pos(2),dir)%l(2) iu = edges_dp(pos(1),pos(2),dir)%u(1) ju = edges_dp(pos(1),pos(2),dir)%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = edges_dp(pos(1),pos(2),pos(3),dir)%l(1) jl = edges_dp(pos(1),pos(2),pos(3),dir)%l(2) kl = edges_dp(pos(1),pos(2),pos(3),dir)%l(3) iu = edges_dp(pos(1),pos(2),pos(3),dir)%u(1) ju = edges_dp(pos(1),pos(2),pos(3),dir)%u(2) ku = edges_dp(pos(1),pos(2),pos(3),dir)%u(3) #endif /* NDIMS == 3 */ #if NDIMS == 2 k = 1 kt = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = kl, ku km1 = k - 1 kp1 = k + 1 ks = 2 * (k - kl) + 1 kt = ks + 1 #endif /* NDIMS == 3 */ do j = jl, ju jm1 = j - 1 jp1 = j + 1 js = 2 * (j - jl) + 1 jt = js + 1 do i = il, iu im1 = i - 1 ip1 = i + 1 is = 2 * (i - il) + 1 it = is + 1 do p = 1, nv dql = qn(p,i ,j,k) - qn(p,im1,j,k) dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) dq(1) = limiter_prol(2.5d-01, dql, dqr) dql = qn(p,i,j ,k) - qn(p,i,jm1,k) dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) dq(2) = limiter_prol(2.5d-01, dql, dqr) #if NDIMS == 3 dql = qn(p,i,j,k ) - qn(p,i,j,km1) dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) dq(3) = limiter_prol(2.5d-01, dql, dqr) #endif /* NDIMS == 3 */ if (positive(p) .and. qn(p,i,j,k) < sum(abs(dq(1:NDIMS)))) then if (qn(p,i,j,k) > 0.0d+00) then do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) dq(:) = 0.5d+00 * dq(:) end do else write(msg,"(a,3i4,a)") & "Positive variable is not positive at (", i, j, k, " )" call print_message(loc, msg) dq(:) = 0.0d+00 end if end if #if NDIMS == 2 dq1 = dq(1) + dq(2) dq2 = dq(1) - dq(2) qb(p,is,js,k ) = qn(p,i,j,k) - dq1 qb(p,it,js,k ) = qn(p,i,j,k) + dq2 qb(p,is,jt,k ) = qn(p,i,j,k) - dq2 qb(p,it,jt,k ) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 2 */ #if NDIMS == 3 dq1 = dq(1) + dq(2) + dq(3) dq2 = dq(1) - dq(2) - dq(3) dq3 = dq(1) - dq(2) + dq(3) dq4 = dq(1) + dq(2) - dq(3) qb(p,is,js,ks) = qn(p,i,j,k) - dq1 qb(p,it,js,ks) = qn(p,i,j,k) + dq2 qb(p,is,jt,ks) = qn(p,i,j,k) - dq3 qb(p,it,jt,ks) = qn(p,i,j,k) + dq4 qb(p,is,js,kt) = qn(p,i,j,k) - dq4 qb(p,it,js,kt) = qn(p,i,j,k) + dq3 qb(p,is,jt,kt) = qn(p,i,j,k) - dq2 qb(p,it,jt,kt) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 3 */ end do ! p end do ! i end do ! j #if NDIMS == 3 end do ! k #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine block_edge_prolong ! !=============================================================================== ! ! BLOCK CORNER UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_CORNER_RESTRICT: ! -------------------------------- ! ! Subroutine returns the corner boundary region by restricting ! the corresponding region from the provided input variable array. ! ! Arguments: ! ! pos - the corner position; ! qn - the input neighbor variable array; ! qb - the output corner boundary array; ! !=============================================================================== ! subroutine block_corner_restrict(pos, qn, qb) ! import external procedures and variables ! use coordinates, only : ng => nghosts use coordinates, only : corners_dr use equations , only : nv ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer , dimension(3) , intent(in) :: pos real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(out) :: qb ! local variables ! integer :: il, ip, iu integer :: jl, jp, ju #if NDIMS == 3 integer :: kl, kp, ku #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! ! prepare indices for the corner region ! #if NDIMS == 2 il = corners_dr(pos(1),pos(2))%l(1) jl = corners_dr(pos(1),pos(2))%l(2) ip = il + 1 jp = jl + 1 iu = corners_dr(pos(1),pos(2))%u(1) ju = corners_dr(pos(1),pos(2))%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = corners_dr(pos(1),pos(2),pos(3))%l(1) jl = corners_dr(pos(1),pos(2),pos(3))%l(2) kl = corners_dr(pos(1),pos(2),pos(3))%l(3) ip = il + 1 jp = jl + 1 kp = kl + 1 iu = corners_dr(pos(1),pos(2),pos(3))%u(1) ju = corners_dr(pos(1),pos(2),pos(3))%u(2) ku = corners_dr(pos(1),pos(2),pos(3))%u(3) #endif /* NDIMS == 3 */ ! restrict the corner region to the output array ! #if NDIMS == 2 qb(1:nv,1:ng,1:ng, : ) = & 2.50d-01 * ((qn(1:nv,il:iu:2,jl:ju:2, : ) & + qn(1:nv,ip:iu:2,jp:ju:2, : )) & + (qn(1:nv,il:iu:2,jp:ju:2, : ) & + qn(1:nv,ip:iu:2,jl:ju:2, : ))) #endif /* NDIMS == 2 */ #if NDIMS == 3 qb(1:nv,1:ng,1:ng,1:ng) = & 1.25d-01 * (((qn(1:nv,il:iu:2,jl:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(1:nv,il:iu:2,jl:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jp:ju:2,kl:ku:2))) & + ((qn(1:nv,il:iu:2,jp:ju:2,kp:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)) & + (qn(1:nv,il:iu:2,jp:ju:2,kl:ku:2) & + qn(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)))) #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine block_corner_restrict ! !=============================================================================== ! ! subroutine BLOCK_CORNER_PROLONG: ! ------------------------------- ! ! Subroutine returns the corner boundary region by prolongating ! the corresponding region from the provided input variable array. ! ! Arguments: ! ! pos - the corner position; ! qn - the input neighbor variable array; ! qb - the output corner boundary array; ! !=============================================================================== ! subroutine block_corner_prolong(pos, qn, qb) use coordinates , only : corners_dp use equations , only : nv, positive use helpers , only : print_message use interpolations, only : limiter_prol implicit none integer , dimension(3) , intent(in) :: pos real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(out) :: qb integer :: p integer :: i, il, iu, is, it, im1, ip1 integer :: j, jl, ju, js, jt, jm1, jp1 integer :: k, kt #if NDIMS == 3 integer :: kl, ku, ks, km1, kp1 #endif /* NDIMS == 3 */ real(kind=8) :: dql, dqr real(kind=8) :: dq1, dq2 #if NDIMS == 3 real(kind=8) :: dq3, dq4 #endif /* NDIMS == 3 */ real(kind=8), dimension(NDIMS) :: dq character(len=80) :: msg character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong()' ! !------------------------------------------------------------------------------- ! #if NDIMS == 2 il = corners_dp(pos(1),pos(2))%l(1) jl = corners_dp(pos(1),pos(2))%l(2) iu = corners_dp(pos(1),pos(2))%u(1) ju = corners_dp(pos(1),pos(2))%u(2) #endif /* NDIMS == 2 */ #if NDIMS == 3 il = corners_dp(pos(1),pos(2),pos(3))%l(1) jl = corners_dp(pos(1),pos(2),pos(3))%l(2) kl = corners_dp(pos(1),pos(2),pos(3))%l(3) iu = corners_dp(pos(1),pos(2),pos(3))%u(1) ju = corners_dp(pos(1),pos(2),pos(3))%u(2) ku = corners_dp(pos(1),pos(2),pos(3))%u(3) #endif /* NDIMS == 3 */ #if NDIMS == 2 k = 1 kt = 1 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = kl, ku km1 = k - 1 kp1 = k + 1 ks = 2 * (k - kl) + 1 kt = ks + 1 #endif /* NDIMS == 3 */ do j = jl, ju jm1 = j - 1 jp1 = j + 1 js = 2 * (j - jl) + 1 jt = js + 1 do i = il, iu im1 = i - 1 ip1 = i + 1 is = 2 * (i - il) + 1 it = is + 1 do p = 1, nv dql = qn(p,i ,j,k) - qn(p,im1,j,k) dqr = qn(p,ip1,j,k) - qn(p,i ,j,k) dq(1) = limiter_prol(2.5d-01, dql, dqr) dql = qn(p,i,j ,k) - qn(p,i,jm1,k) dqr = qn(p,i,jp1,k) - qn(p,i,j ,k) dq(2) = limiter_prol(2.5d-01, dql, dqr) #if NDIMS == 3 dql = qn(p,i,j,k ) - qn(p,i,j,km1) dqr = qn(p,i,j,kp1) - qn(p,i,j,k ) dq(3) = limiter_prol(2.5d-01, dql, dqr) #endif /* NDIMS == 3 */ if (positive(p) .and. qn(p,i,j,k) < sum(abs(dq(1:NDIMS)))) then if (qn(p,i,j,k) > 0.0d+00) then do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) dq(:) = 0.5d+00 * dq(:) end do else write(msg,"(a,3i4,a)") & "Positive variable is not positive at (", i, j, k, " )" call print_message(loc, msg) dq(:) = 0.0d+00 end if end if #if NDIMS == 2 dq1 = dq(1) + dq(2) dq2 = dq(1) - dq(2) qb(p,is,js,k ) = qn(p,i,j,k) - dq1 qb(p,it,js,k ) = qn(p,i,j,k) + dq2 qb(p,is,jt,k ) = qn(p,i,j,k) - dq2 qb(p,it,jt,k ) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 2 */ #if NDIMS == 3 dq1 = dq(1) + dq(2) + dq(3) dq2 = dq(1) - dq(2) - dq(3) dq3 = dq(1) - dq(2) + dq(3) dq4 = dq(1) + dq(2) - dq(3) qb(p,is,js,ks) = qn(p,i,j,k) - dq1 qb(p,it,js,ks) = qn(p,i,j,k) + dq2 qb(p,is,jt,ks) = qn(p,i,j,k) - dq3 qb(p,it,jt,ks) = qn(p,i,j,k) + dq4 qb(p,is,js,kt) = qn(p,i,j,k) - dq4 qb(p,it,js,kt) = qn(p,i,j,k) + dq3 qb(p,is,jt,kt) = qn(p,i,j,k) - dq2 qb(p,it,jt,kt) = qn(p,i,j,k) + dq1 #endif /* NDIMS == 3 */ end do ! p end do ! i end do ! j #if NDIMS == 3 end do ! k #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine block_corner_prolong ! !=============================================================================== ! ! OTHER BOUNDARY SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine UPDATE_GHOST_CELLS: ! ----------------------------- ! ! Subroutine updates conservative variables in all ghost cells from ! already updated primitive variables. ! ! !=============================================================================== ! subroutine update_ghost_cells() ! include external variables ! use blocks , only : block_data, list_data use coordinates, only : nn => bcells, nb, ne, nbl, neu use equations , only : prim2cons ! local variables are not implicit by default ! implicit none ! local variables ! integer :: i, j, k = 1 ! local pointers ! type(block_data), pointer :: pdata ! !------------------------------------------------------------------------------- ! ! assign the pointer to the first block on the list ! pdata => list_data ! scan all data blocks until the last is reached ! do while(associated(pdata)) ! update the X and Y boundary ghost cells ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ ! update lower layers of the Y boundary ! do j = 1, nbl call prim2cons(pdata%q(:,:,j,k), pdata%u(:,:,j,k), .true.) end do ! j = 1, nbl ! update upper layers of the Y boundary ! do j = neu, nn call prim2cons(pdata%q(:,:,j,k), pdata%u(:,:,j,k), .true.) end do ! j = neu, nn ! update remaining left layers of the X boundary ! do i = 1, nbl call prim2cons(pdata%q(:,i,nb:ne,k), pdata%u(:,i,nb:ne,k), .true.) end do ! i = 1, nbl ! update remaining right layers of the X boundary ! do i = neu, nn call prim2cons(pdata%q(:,i,nb:ne,k), pdata%u(:,i,nb:ne,k), .true.) end do ! i = neu, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ #if NDIMS == 3 ! update the Z boundary ghost cells ! do j = nb, ne ! update the remaining front layers of the Z boundary ! do k = 1, nbl call prim2cons(pdata%q(:,nb:ne,j,k), pdata%u(:,nb:ne,j,k), .true.) end do ! k = 1, nbl ! update the remaining back layers of the Z boundary ! do k = neu, nn call prim2cons(pdata%q(:,nb:ne,j,k), pdata%u(:,nb:ne,j,k), .true.) end do ! k = neu, nn end do ! j = nb, ne #endif /* NDIMS == 3 */ ! assign the pointer to the next block on the list ! pdata => pdata%next end do ! data blocks !------------------------------------------------------------------------------- ! end subroutine update_ghost_cells #ifdef MPI ! !=============================================================================== ! ! subroutine PREPARE_EXCHANGE_ARRAY: ! --------------------------------- ! ! Subroutine prepares the arrays for block exchange lists and their counters. ! ! !=============================================================================== ! subroutine prepare_exchange_array() ! include external variables ! use mpitools, only : npmax ! local variables are not implicit by default ! implicit none ! local variables ! integer :: icol, irow ! !------------------------------------------------------------------------------- ! ! iterate over all elements of the block exchange array ! do irow = 0, npmax do icol = 0, npmax ! nullify the array element pointer ! nullify(barray(irow,icol)%ptr) ! reset the corresponding counter ! bcount(irow,icol) = 0 end do ! icol = 0, npmax end do ! irow = 0, npmax !------------------------------------------------------------------------------- ! end subroutine prepare_exchange_array ! !=============================================================================== ! ! subroutine RELEASE_EXCHANGE_ARRAY: ! --------------------------------- ! ! Subroutine releases objects on the array of block exchange lists. ! ! !=============================================================================== ! subroutine release_exchange_array() ! include external variables ! use blocks , only : block_info, pointer_info use mpitools, only : npmax ! local variables are not implicit by default ! implicit none ! local variables ! integer :: icol, irow ! local pointers ! type(block_info), pointer :: pinfo ! !------------------------------------------------------------------------------- ! ! iterate over all elements of the block exchange array ! do irow = 0, npmax do icol = 0, npmax ! associate pinfo with the first block in the exchange list ! pinfo => barray(irow,icol)%ptr ! scan all elements on the exchange list ! do while(associated(pinfo)) ! associate the exchange list pointer ! barray(irow,icol)%ptr => pinfo%prev ! nullify pointer fields ! nullify(pinfo%prev) nullify(pinfo%next) nullify(pinfo%meta) nullify(pinfo%neigh) ! deallocate info block ! deallocate(pinfo) ! associate pinfo with the next block ! pinfo => barray(irow,icol)%ptr end do ! %ptr blocks end do ! icol = 0, npmax end do ! irow = 0, npmax !------------------------------------------------------------------------------- ! end subroutine release_exchange_array ! !=============================================================================== ! ! subroutine APPEND_EXCHANGE_BLOCK: ! --------------------------------- ! ! Subroutine appends an info block to the element of array of block ! exchange lists. The element is determined by the processes of the meta ! and neighbor blocks. ! ! Arguments: ! ! pmeta - the pointer to meta block; ! pneigh - the pointer to the neighbor of pmeta; ! dir - the direction of the neighbor; ! pos - the position of the neighbor; ! !=============================================================================== ! subroutine append_exchange_block(pmeta, pneigh, dir, pos) ! include external variables ! use blocks, only : block_info, block_meta ! local variables are not implicit by default ! implicit none ! subroutine arguments ! type(block_meta), pointer, intent(inout) :: pmeta, pneigh integer , intent(in) :: dir integer, dimension(3) , intent(in) :: pos ! local variables ! integer :: icol, irow ! local pointers ! type(block_info), pointer :: pinfo ! !------------------------------------------------------------------------------- ! ! get the column and row indices ! irow = pneigh%process icol = pmeta%process ! increase the counter for the number of blocks to exchange ! bcount(irow,icol) = bcount(irow,icol) + 1 ! allocate a new info object ! allocate(pinfo) ! fill out its fields ! pinfo%meta => pmeta pinfo%neigh => pneigh pinfo%direction = dir pinfo%corner(1:NDIMS) = pos(1:NDIMS) pinfo%level_difference = pmeta%level - pneigh%level ! nullify pointer fields ! nullify(pinfo%prev) nullify(pinfo%next) ! check if the list is empty ! if (associated(barray(irow,icol)%ptr)) then ! if it is, associate the newly created block with it ! pinfo%prev => barray(irow,icol)%ptr end if ! %ptr associated ! point the list to the newly created block ! barray(irow,icol)%ptr => pinfo !------------------------------------------------------------------------------- ! end subroutine append_exchange_block #endif /* MPI */ !=============================================================================== ! end module