!!****************************************************************************** !! !! 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-2023 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 ! interfaces for custom boundary conditions ! abstract interface subroutine custom_boundary_iface(idir, il, iu, jl, ju, t, dt, x, y, z, qn) implicit none integer , intent(in) :: idir, il, iu, jl, ju real(kind=8) , intent(in) :: t, dt real(kind=8), dimension(:) , intent(in) :: x, y, z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn end subroutine end interface ! pointers to custom boundary conditions, one per direction handling both sides ! procedure(custom_boundary_iface), pointer, save :: & custom_boundary_x => null(), & custom_boundary_y => null(), & custom_boundary_z => null() ! supported boundary types ! enum, bind(c) enumerator bnd_custom enumerator bnd_periodic enumerator bnd_open enumerator bnd_outflow enumerator bnd_reflective enumerator bnd_gravity end enum ! the boundary type array for both sides along all directions ! integer, dimension(3,2), save :: bnd_type = bnd_periodic #ifdef MPI ! arrays to store information about blocks which need to exchange ghost zones ! between processes ! type(pointer_info), dimension(:,:), allocatable, save :: barray integer , dimension(:,:), allocatable, save :: bcount #endif /* MPI */ private public :: initialize_boundaries, finalize_boundaries, print_boundaries public :: boundary_fluxes, boundary_variables public :: custom_boundary_x, custom_boundary_y, custom_boundary_z !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_BOUNDARIES: ! -------------------------------- ! ! Subroutine initializes the module BOUNDARIES. ! ! Arguments: ! ! status - the subroutine call status; ! !=============================================================================== ! subroutine initialize_boundaries(status) use coordinates, only : periodic use helpers , only : print_message #ifdef MPI use mpitools , only : npmax #endif /* MPI */ use parameters , only : get_parameter implicit none integer, intent(out) :: status character(len=64) :: str, bnd integer :: n, s character(len=*), parameter :: loc = 'BOUNDARIES::initialize_boundaries()' !------------------------------------------------------------------------------- ! status = 0 do n = 1, NDIMS do s = 1, 2 bnd = "periodic" if (s == 1) then str = char(119+n) // "lbndry" else str = char(119+n) // "ubndry" end if call get_parameter(str, bnd) if (s == 1) then str = char(119+n) // "_lower_boundary" else str = char(119+n) // "_upper_boundary" end if call get_parameter(str, bnd) select case(bnd) case("open") bnd_type(n,s) = bnd_open case("outflow", "out") bnd_type(n,s) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(n,s) = bnd_reflective case("hydrostatic", "gravity") bnd_type(n,s) = bnd_gravity case("user", "custom") bnd_type(n,s) = bnd_custom case default bnd_type(n,s) = bnd_periodic end select end do if ((bnd_type(n,1) == bnd_periodic .and. & bnd_type(n,2) /= bnd_periodic) .or. & (bnd_type(n,2) == bnd_periodic .and. & bnd_type(n,1) /= bnd_periodic)) then call print_message(loc, char(87+n) // & "-boundary cannot be periodic on one " // & "side and non-periodic on another!") status = 1 return end if periodic(n) = (bnd_type(n,1) == bnd_periodic) .and. & (bnd_type(n,2) == bnd_periodic) end do #ifdef MPI ! allocate the arrays for the list of data blocks at different processes which ! need to exchange boundaries, and for the count of the blocks in each list ! allocate(barray(0:npmax,0:npmax), bcount(0:npmax,0:npmax), stat=status) ! generate the list of blocks on different processes ! if (status == 0) call prepare_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine initialize_boundaries ! !=============================================================================== ! ! subroutine FINALIZE_BOUNDARIES: ! ------------------------------ ! ! Subroutine finalizes module. ! ! Arguments: ! ! status - the subroutine call status; ! !=============================================================================== ! subroutine finalize_boundaries(status) use helpers, only : print_message implicit none integer, intent(out) :: status character(len=*), parameter :: loc = 'BOUNDARIES::finalize_boundaries()' !------------------------------------------------------------------------------- ! status = 0 #ifdef MPI call release_exchange_array() deallocate(barray, bcount, stat=status) if (status /= 0) then call print_message(loc, "Could not deallocate block exchange arrays!") end if #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine finalize_boundaries ! !=============================================================================== ! ! subroutine PRINT_BOUNDARIES: ! --------------------------- ! ! Subroutine prints module parameters and setup. ! ! Arguments: ! ! verbose - the verbose flag; ! !=============================================================================== ! subroutine print_boundaries(verbose) use helpers, only : print_section, print_parameter implicit none logical, intent(in) :: verbose integer :: n, w character(len=64), dimension(0:8) :: bnd_names !------------------------------------------------------------------------------- ! if (verbose) then bnd_names(:) = "unknown" bnd_names(bnd_custom) = "custom" bnd_names(bnd_periodic) = "periodic" bnd_names(bnd_open) = "open" bnd_names(bnd_outflow) = "outflow" bnd_names(bnd_reflective) = "reflective" bnd_names(bnd_gravity) = "hydrostatic" w = 4 do n = 1, NDIMS w = max(w, max(len(trim(bnd_names(bnd_type(n,1)))), & len(trim(bnd_names(bnd_type(n,2)))))) end do call print_section(verbose, "Boundaries") do n = 1, NDIMS call print_parameter(verbose, char(87+n) // "-boundary", & bnd_names(bnd_type(n,1)), & bnd_names(bnd_type(n,2)), w) end do end if !------------------------------------------------------------------------------- ! end subroutine print_boundaries ! !=============================================================================== ! ! subroutine BOUNDARY_FLUXES: ! -------------------------- ! ! Subroutine corrects the block interface fluxes from neighbors at ! higher levels. ! ! Arguments: ! ! status - the subroutine call status; ! !=============================================================================== ! subroutine boundary_fluxes(status) 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 : nsides use coordinates, only : minlev, maxlev use coordinates, only : nh => ncells_half use coordinates, only : nb, ne, nbm, nbp, nem, 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 use helpers , only : print_message #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(out) :: status 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 */ integer :: n #if NDIMS == 2 integer :: m #endif /* NDIMS == 2 */ integer :: i, il, iu integer :: j, jl, ju integer :: k, kl, ku integer :: s #ifdef MPI integer :: sproc, rproc integer :: scount, rcount integer :: l, p #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 */ #ifdef MPI real(kind=8), dimension(:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims character(len=*), parameter :: loc = 'BOUNDARIES::boundary_fluxes()' !------------------------------------------------------------------------------- ! status = 0 if (minlev == maxlev) return if (first) then nlims(1,1) = nb nlims(2,1) = nb + nh nlims(:,2) = nlims(:,1) + nh - 1 first = .false. end if #if NDIMS == 2 k = 1 kl = 1 ku = 1 #endif /* NDIMS == 2 */ do n = 1, NDIMS #if NDIMS == 2 m = 3 - n #endif /* NDIMS == 2 */ #ifdef MPI sproc = 0 rproc = 0 scount = 0 rcount = 0 call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta pdata => pmeta%data #if NDIMS == 3 do k = 1, nsides #endif /* NDIMS == 3 */ do j = 1, nsides do i = 1, nsides #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 */ if (associated(pneigh)) then if (pneigh%level > pmeta%level) then #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ select case(n) case(1) s = 3 - i if (i == 1) then il = nbm iu = nb else il = ne iu = nep end if jl = nlims(j,1) ju = nlims(j,2) #if NDIMS == 3 kl = nlims(k,1) ku = nlims(k,2) fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fx(:,nb:nem:2,nb:nem:2,s) & + pneigh%data%fx(:,nbp:ne:2,nbp:ne:2,s)) & + (pneigh%data%fx(:,nbp:ne:2,nb:nem:2,s) & + pneigh%data%fx(:,nb:nem:2,nbp:ne:2,s))) #else /* NDIMS == 3 */ fh(:,:,:) = & 5.0d-01 * (pneigh%data%fx(:,nb:nem:2,:,s) & + pneigh%data%fx(:,nbp:ne:2,:,s)) #endif /* NDIMS == 3 */ df(:,:,:) = (fh(:,:,:) - pdata%fx(:,jl:ju,kl:ku,i)) & * adxi(pmeta%level) pdata%du(:,il,jl:ju,kl:ku) = & pdata%du(:,il,jl:ju,kl:ku) - df(:,:,:) pdata%du(:,iu,jl:ju,kl:ku) = & pdata%du(:,iu,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%u(s,il,jl:ju,kl:ku) & + sr(:,:) * pdata%u(s,iu,jl:ju,kl:ku) pdata%du(s,il,jl:ju,kl:ku) = & pdata%du(s,il,jl:ju,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,iu,jl:ju,kl:ku) = & pdata%du(s,iu,jl:ju,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if case(2) s = 3 - j if (j == 1) then jl = nbm ju = nb else jl = ne ju = nep end if il = nlims(i,1) iu = nlims(i,2) #if NDIMS == 3 kl = nlims(k,1) ku = nlims(k,2) fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fy(:,nb:nem:2,nb:nem:2,s) & + pneigh%data%fy(:,nbp:ne:2,nbp:ne:2,s)) & + (pneigh%data%fy(:,nbp:ne:2,nb:nem:2,s) & + pneigh%data%fy(:,nb:nem:2,nbp:ne:2,s))) #else /* NDIMS == 3 */ fh(:,:,:) = & 5.0d-01 * (pneigh%data%fy(:,nb:nem:2,:,s) & + pneigh%data%fy(:,nbp:ne:2,:,s)) #endif /* NDIMS == 3 */ df(:,:,:) = (fh(:,:,:) - pdata%fy(:,il:iu,kl:ku,j)) & * adyi(pmeta%level) pdata%du(:,il:iu,jl,kl:ku) = & pdata%du(:,il:iu,jl,kl:ku) - df(:,:,:) pdata%du(:,il:iu,ju,kl:ku) = & pdata%du(:,il:iu,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%u(s,il:iu,jl,kl:ku) & + sr(:,:) * pdata%u(s,il:iu,ju,kl:ku) pdata%du(s,il:iu,jl,kl:ku) = & pdata%du(s,il:iu,jl,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,ju,kl:ku) = & pdata%du(s,il:iu,ju,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if #if NDIMS == 3 case(3) s = 3 - k if (k == 1) then kl = nbm ku = nb else kl = ne ku = nep end if il = nlims(i,1) iu = nlims(i,2) jl = nlims(j,1) ju = nlims(j,2) fh(:,:,:) = & 2.5d-01 * ((pneigh%data%fz(:,nb:nem:2,nb:nem:2,s) & + pneigh%data%fz(:,nbp:ne:2,nbp:ne:2,s)) & + (pneigh%data%fz(:,nbp:ne:2,nb:nem:2,s) & + pneigh%data%fz(:,nb:nem:2,nbp:ne:2,s))) df(:,:,:) = (fh(:,:,:) - pdata%fz(:,il:iu,jl:ju,k)) & * adzi(pmeta%level) pdata%du(:,il:iu,jl:ju,kl) = & pdata%du(:,il:iu,jl:ju,kl) - df(:,:,:) pdata%du(:,il:iu,jl:ju,ku) = & pdata%du(:,il:iu,jl:ju,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%u(s,il:iu,jl:ju,kl) & + sr(:,:) * pdata%u(s,il:iu,jl:ju,ku) pdata%du(s,il:iu,jl:ju,kl) = & pdata%du(s,il:iu,jl:ju,kl) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,jl:ju,ku) = & pdata%du(s,il:iu,jl:ju,ku) & + df(idn,:,:) * ps(:,:) end do end if #endif /* NDIMS == 3 */ end select #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, n, [ i, j, k ]) end if #endif /* MPI */ end if end if end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) if (scount > 0 .or. rcount > 0) then #if NDIMS == 3 allocate(sbuf(nf,nh,nh,scount), rbuf(nf,nh,nh,rcount), stat=status) #else /* NDIMS == 3 */ allocate(sbuf(nf,nh, 1,scount), rbuf(nf,nh, 1,rcount), stat=status) #endif /* NDIMS == 3 */ if (status /= 0) & call print_message(loc, "Could not allocate exchange arrays!") if (scount > 0) then l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%neigh%data select case(pinfo%direction) case(1) s = 3 - pinfo%location(1) #if NDIMS == 3 sbuf(:,:,:,l) = 2.5d-01 * ((pdata%fx(:,nb:nem:2,nb:nem:2,s) & + pdata%fx(:,nbp:ne:2,nbp:ne:2,s)) & + (pdata%fx(:,nbp:ne:2,nb:nem:2,s) & + pdata%fx(:,nb:nem:2,nbp:ne:2,s))) #else /* NDIMS == 3 */ sbuf(:,:,:,l) = 5.0d-01 * (pdata%fx(:,nb:nem:2,:,s) & + pdata%fx(:,nbp:ne:2,:,s)) #endif /* NDIMS == 3 */ case(2) s = 3 - pinfo%location(2) #if NDIMS == 3 sbuf(:,:,:,l) = 2.5d-01 * ((pdata%fy(:,nb:nem:2,nb:nem:2,s) & + pdata%fy(:,nbp:ne:2,nbp:ne:2,s)) & + (pdata%fy(:,nbp:ne:2,nb:nem:2,s) & + pdata%fy(:,nb:nem:2,nbp:ne:2,s))) #else /* NDIMS == 3 */ sbuf(:,:,:,l) = 5.0d-01 * (pdata%fy(:,nb:nem:2,:,s) & + pdata%fy(:,nbp:ne:2,:,s)) #endif /* NDIMS == 3 */ #if NDIMS == 3 case(3) s = 3 - pinfo%location(3) sbuf(:,:,:,l) = 2.5d-01 * ((pdata%fz(:,nb:nem:2,nb:nem:2,s) & + pdata%fz(:,nbp:ne:2,nbp:ne:2,s)) & + (pdata%fz(:,nbp:ne:2,nb:nem:2,s) & + pdata%fz(:,nb:nem:2,nbp:ne:2,s))) #endif /* NDIMS == 3 */ end select pinfo => pinfo%prev end do ! %ptr blocks end if call exchange_arrays(rproc, p, sbuf, rbuf) if (rcount > 0) then l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pmeta => pinfo%meta pdata => pmeta%data i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ select case(n) case(1) if (i == 1) then il = nbm iu = nb else il = ne iu = nep end if jl = nlims(j,1) ju = nlims(j,2) #if NDIMS == 3 kl = nlims(k,1) ku = nlims(k,2) #endif /* NDIMS == 3 */ df(:,:,:) = (rbuf(:,:,:,l) - pdata%fx(:,jl:ju,kl:ku,i)) & * adxi(pmeta%level) pdata%du(:,il,jl:ju,kl:ku) = pdata%du(:,il,jl:ju,kl:ku) & - df(:,:,:) pdata%du(:,iu,jl:ju,kl:ku) = pdata%du(:,iu,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%u(s,il,jl:ju,kl:ku) & + sr(:,:) * pdata%u(s,iu,jl:ju,kl:ku) pdata%du(s,il,jl:ju,kl:ku) = pdata%du(s,il,jl:ju,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,iu,jl:ju,kl:ku) = pdata%du(s,iu,jl:ju,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if case(2) if (j == 1) then jl = nbm ju = nb else jl = ne ju = nep end if il = nlims(i,1) iu = nlims(i,2) #if NDIMS == 3 kl = nlims(k,1) ku = nlims(k,2) #endif /* NDIMS == 3 */ df(:,:,:) = (rbuf(:,:,:,l) - pdata%fy(:,il:iu,kl:ku,j)) & * adyi(pmeta%level) pdata%du(:,il:iu,jl,kl:ku) = pdata%du(:,il:iu,jl,kl:ku) & - df(:,:,:) pdata%du(:,il:iu,ju,kl:ku) = pdata%du(:,il:iu,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%u(s,il:iu,jl,kl:ku) & + sr(:,:) * pdata%u(s,il:iu,ju,kl:ku) pdata%du(s,il:iu,jl,kl:ku) = pdata%du(s,il:iu,jl,kl:ku) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,ju,kl:ku) = pdata%du(s,il:iu,ju,kl:ku) & + df(idn,:,:) * ps(:,:) end do end if #if NDIMS == 3 case(3) if (k == 1) then kl = nbm ku = nb else kl = ne ku = nep end if il = nlims(i,1) iu = nlims(i,2) jl = nlims(j,1) ju = nlims(j,2) df(:,:,:) = (rbuf(:,:,:,l) - pdata%fz(:,il:iu,jl:ju,k)) & * adzi(pmeta%level) pdata%du(:,il:iu,jl:ju,kl) = pdata%du(:,il:iu,jl:ju,kl) & - df(:,:,:) pdata%du(:,il:iu,jl:ju,ku) = pdata%du(:,il:iu,jl:ju,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%u(s,il:iu,jl:ju,kl) & + sr(:,:) * pdata%u(s,il:iu,jl:ju,ku) pdata%du(s,il:iu,jl:ju,kl) = pdata%du(s,il:iu,jl:ju,kl) & - df(idn,:,:) * ps(:,:) pdata%du(s,il:iu,jl:ju,ku) = pdata%du(s,il:iu,jl:ju,ku) & + df(idn,:,:) * ps(:,:) end do end if #endif /* NDIMS == 3 */ end select pinfo => pinfo%prev end do end if deallocate(sbuf, rbuf, stat=status) if (status /= 0) & call print_message(loc, "Could not deallocate exchange arrays!") end if end if end do call release_exchange_array() #endif /* MPI */ end do !------------------------------------------------------------------------------- ! end subroutine boundary_fluxes ! !=============================================================================== ! ! 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; ! status - the call status; ! !=============================================================================== ! subroutine boundary_variables(t, dt, status) use coordinates, only : minlev, maxlev implicit none real(kind=8), intent(in) :: t, dt integer , intent(out) :: status integer :: level !------------------------------------------------------------------------------- ! #if NDIMS == 3 call boundaries_face_copy(status) #endif /* NDIMS == 3 */ call boundaries_edge_copy(status) call boundaries_corner_copy(status) if (minlev /= maxlev) then #if NDIMS == 3 call boundaries_face_restrict(status) #endif /* NDIMS == 3 */ call boundaries_edge_restrict(status) call boundaries_corner_restrict(status) call boundaries_specific(t, dt, status) do level = minlev, maxlev #if NDIMS == 3 call boundaries_face_prolong(level, status) #endif /* NDIMS == 3 */ call boundaries_edge_prolong(level, status) call boundaries_corner_prolong(level, status) end do end if ! minlev /= maxlev call boundaries_specific(t, dt, status) !------------------------------------------------------------------------------- ! end subroutine boundary_variables ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== #if NDIMS == 3 ! !=============================================================================== ! ! DOMAIN FACE BOUNDARY UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BOUNDARIES_FACE_COPY: ! ------------------------------- ! ! Subroutine updates the blocks' face ghost zones from blocks at ! the same level. ! ! Arguments: ! ! status - the call status; ! !=============================================================================== ! subroutine boundaries_face_copy(status) 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 : bcells, ncells_half, nghosts, nb, ne #ifdef MPI use equations , only : nv #endif /* MPI */ use helpers , only : print_message #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(out) :: status 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 */ integer :: n integer :: i, il, iu, is, it integer :: j, jl, ju, js, jt integer :: k, kl, ku, ks, kt #ifdef MPI integer :: sproc, rproc integer :: ecount integer :: l, p integer, dimension(1) :: dm integer, dimension(4) :: pm real(kind=8), dimension(:,:), allocatable :: buf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims, dlims, tlims character(len=*), parameter :: loc = 'BOUNDARIES::boundaries_face_copy()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = nb - 1 nlims(2,1) = ne + 1 nlims(2,2) = bcells dlims(1,1) = ne - nghosts + 1 dlims(1,2) = ne dlims(2,1) = nb dlims(2,2) = nb + nghosts - 1 tlims(1,1) = nb tlims(1,2) = nb + ncells_half - 1 tlims(2,1) = ne - ncells_half + 1 tlims(2,2) = ne first = .false. end if status = 0 #ifdef MPI sproc = 0 rproc = 0 ecount = 0 dm(1) = nv * ncells_half**2 * nghosts pm(1) = nv call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta pdata => pmeta%data do n = 1, 3 do k = 1, nsides if (n == 3) then kl = nlims(k,1) ku = nlims(k,2) ks = dlims(k,1) kt = dlims(k,2) else kl = tlims(k,1) ku = tlims(k,2) ks = tlims(k,1) kt = tlims(k,2) end if do j = 1, nsides if (n == 2) then jl = nlims(j,1) ju = nlims(j,2) js = dlims(j,1) jt = dlims(j,2) else jl = tlims(j,1) ju = tlims(j,2) js = tlims(j,1) jt = tlims(j,2) end if do i = 1, nsides if (n == 1) then il = nlims(i,1) iu = nlims(i,2) is = dlims(i,1) it = dlims(i,2) else il = tlims(i,1) iu = tlims(i,2) is = tlims(i,1) it = tlims(i,2) end if pneigh => pmeta%faces(i,j,k,n)%ptr if (associated(pneigh)) then if (pneigh%level == pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ pdata%u(:,il:iu,jl:ju,kl:ku) = & pneigh%data%u(:,is:it,js:jt,ks:kt) #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, n, [ i, j, k ]) end if #endif /* MPI */ end if end if end if end do end do end do end do pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 if (ecount > 0) then allocate(buf(dm(1),ecount)) if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%neigh%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) k = pinfo%location(3) if (n == 1) then is = dlims(i,1) it = dlims(i,2) else is = tlims(i,1) it = tlims(i,2) end if if (n == 2) then js = dlims(j,1) jt = dlims(j,2) else js = tlims(j,1) jt = tlims(j,2) end if if (n == 3) then ks = dlims(k,1) kt = dlims(k,2) else ks = tlims(k,1) kt = tlims(k,2) end if buf(:,l) = reshape(pdata%u(:,is:it,js:jt,ks:kt), dm) pinfo => pinfo%prev end do call exchange_arrays(rproc, p, buf) l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%meta%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) k = pinfo%location(3) if (n == 1) then pm(2) = nghosts il = nlims(i,1) iu = nlims(i,2) else pm(2) = ncells_half il = tlims(i,1) iu = tlims(i,2) end if if (n == 2) then pm(3) = nghosts jl = nlims(j,1) ju = nlims(j,2) else pm(3) = ncells_half jl = tlims(j,1) ju = tlims(j,2) end if if (n == 3) then pm(4) = nghosts kl = nlims(k,1) ku = nlims(k,2) else pm(4) = ncells_half kl = tlims(k,1) ku = tlims(k,2) end if pdata%u(:,il:iu,jl:ju,kl:ku) = reshape(buf(:,l), pm) pinfo => pinfo%prev end do deallocate(buf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_face_copy ! !=============================================================================== ! ! subroutine BOUNDARIES_FACE_RESTRICT: ! ----------------------------------- ! ! Subroutine updates the blocks' face ghost zones from blocks at ! higher levels. ! ! Arguments: ! ! status - the call status; ! !=============================================================================== ! subroutine boundaries_face_restrict(status) 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 : bcells, ncells_half, nb, ne #ifdef MPI use coordinates, only : nghosts use equations , only : nv #endif /* MPI */ use helpers , only : print_message #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(out) :: status 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 */ integer :: n integer :: i, il, iu integer :: j, jl, ju integer :: k, kl, ku #ifdef MPI integer :: sproc, rproc integer :: scount, rcount integer :: l, p integer, dimension(1) :: dm integer, dimension(3) :: bm integer, dimension(4) :: pm real(kind=8), dimension(:,:,:,:), allocatable :: tmp real(kind=8), dimension(:,:) , allocatable :: sbuf, rbuf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims, tlims character(len=*), parameter :: loc = & 'BOUNDARIES::boundaries_face_restrict()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = nb - 1 nlims(2,1) = ne + 1 nlims(2,2) = bcells tlims(1,1) = nb tlims(1,2) = nb + ncells_half - 1 tlims(2,1) = ne - ncells_half + 1 tlims(2,2) = ne first = .false. end if status = 0 #ifdef MPI sproc = 0 rproc = 0 scount = 0 rcount = 0 dm(1) = nv * ncells_half**2 * nghosts pm(1) = nv call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta pdata => pmeta%data do n = 1, 3 do k = 1, nsides if (n == 3) then kl = nlims(k,1) ku = nlims(k,2) else kl = tlims(k,1) ku = tlims(k,2) end if do j = 1, nsides if (n == 2) then jl = nlims(j,1) ju = nlims(j,2) else jl = tlims(j,1) ju = tlims(j,2) end if do i = 1, nsides if (n == 1) then il = nlims(i,1) iu = nlims(i,2) else il = tlims(i,1) iu = tlims(i,2) end if pneigh => pmeta%faces(i,j,k,n)%ptr if (associated(pneigh)) then if (pneigh%level > pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ call block_face_restrict(n, [ i, j, k ], pneigh%data%u,& pdata%u(:,il:iu,jl:ju,kl:ku)) #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, n, [ i, j, k ]) end if #endif /* MPI */ end if end if end if end do end do end do end do pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) if (scount > 0 .or. rcount > 0) then allocate(sbuf(dm(1),scount), rbuf(dm(1),rcount), stat=status) if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") if (scount > 0) then l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%neigh%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ bm(:) = ncells_half bm(n) = nghosts allocate(tmp(nv,bm(1),bm(2),bm(3)), stat=status) if (status /= 0) & call print_message(loc, & "Could not allocate the temporary buffer!") call block_face_restrict(n, [ i, j, k ], pdata%u, tmp) sbuf(:,l) = reshape(tmp, dm) deallocate(tmp, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the temporary buffer!") pinfo => pinfo%prev end do end if call exchange_arrays(rproc, p, sbuf, rbuf) if (rcount > 0) then l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%meta%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) k = pinfo%location(3) if (n == 1) then pm(2) = nghosts il = nlims(i,1) iu = nlims(i,2) else pm(2) = ncells_half il = tlims(i,1) iu = tlims(i,2) end if if (n == 2) then pm(3) = nghosts jl = nlims(j,1) ju = nlims(j,2) else pm(3) = ncells_half jl = tlims(j,1) ju = tlims(j,2) end if if (n == 3) then pm(4) = nghosts kl = nlims(k,1) ku = nlims(k,2) else pm(4) = ncells_half kl = tlims(k,1) ku = tlims(k,2) end if pdata%u(:,il:iu,jl:ju,kl:ku) = reshape(rbuf(:,l), pm) pinfo => pinfo%prev end do end if deallocate(sbuf, rbuf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_face_restrict ! !=============================================================================== ! ! subroutine BOUNDARIES_FACE_PROLONG: ! ---------------------------------- ! ! Subroutine updates the blocks' face ghost zones from blocks at lower levels. ! ! Arguments: ! ! level - the level to process; ! status - the call status; ! !=============================================================================== ! subroutine boundaries_face_prolong(level, status) 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 : bcells, nb, ne #ifdef MPI use coordinates, only : ncells, nghosts use equations , only : nv #endif /* MPI */ use helpers , only : print_message #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(in) :: level integer, intent(out) :: status 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 */ integer :: n integer :: i, il, iu integer :: j, jl, ju integer :: k, kl, ku integer, dimension(3) :: bm #ifdef MPI integer :: sproc, rproc integer :: scount, rcount integer :: l, p integer, dimension(1) :: dm integer, dimension(4) :: pm real(kind=8), dimension(:,:,:,:), allocatable :: tmp real(kind=8), dimension(:,:) , allocatable :: sbuf, rbuf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims, tlims character(len=*), parameter :: loc = 'BOUNDARIES::boundaries_face_prolong()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = nb - 1 nlims(2,1) = ne + 1 nlims(2,2) = bcells tlims(1,1) = nb tlims(1,2) = bcells tlims(2,1) = 1 tlims(2,2) = ne first = .false. end if status = 0 #ifdef MPI sproc = 0 rproc = 0 scount = 0 rcount = 0 dm(1) = nv * (ncells + nghosts)**2 * nghosts pm(1) = nv call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta if (pmeta%level == level) then pdata => pmeta%data do n = 1, 3 do k = 1, nsides if (n == 3) then bm(3) = k kl = nlims(k,1) ku = nlims(k,2) else bm(3) = pmeta%pos(3) + 1 kl = tlims(bm(3),1) ku = tlims(bm(3),2) end if do j = 1, nsides if (n == 2) then bm(2) = j jl = nlims(j,1) ju = nlims(j,2) else bm(2) = pmeta%pos(2) + 1 jl = tlims(bm(2),1) ju = tlims(bm(2),2) end if do i = 1, nsides if (n == 1) then bm(1) = i il = nlims(i,1) iu = nlims(i,2) else bm(1) = pmeta%pos(1) + 1 il = tlims(bm(1),1) iu = tlims(bm(1),2) end if pneigh => pmeta%faces(i,j,k,n)%ptr if (associated(pneigh)) then if (pneigh%level < pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ call block_face_prolong(n, bm, pneigh%data%u, & pdata%u(:,il:iu,jl:ju,kl:ku), status) #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, n, bm) end if #endif /* MPI */ end if end if end if end do end do end do end do end if pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) if (scount > 0 .or. rcount > 0) then allocate(sbuf(dm(1),scount), rbuf(dm(1),rcount), stat=status) if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") if (scount > 0) then l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%neigh%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ bm(:) = ncells + nghosts bm(n) = nghosts allocate(tmp(nv,bm(1),bm(2),bm(3)), stat=status) if (status /= 0) & call print_message(loc, & "Could not allocate the temporary buffer!") call block_face_prolong(n, [ i, j, k ], pdata%u, tmp, status) sbuf(:,l) = reshape(tmp, dm) deallocate(tmp, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the temporary buffer!") pinfo => pinfo%prev end do end if call exchange_arrays(rproc, p, sbuf, rbuf) if (rcount > 0) then l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%meta%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) k = pinfo%location(3) if (n == 1) then pm(2) = nghosts il = nlims(i,1) iu = nlims(i,2) else pm(2) = ncells + nghosts il = tlims(i,1) iu = tlims(i,2) end if if (n == 2) then pm(3) = nghosts jl = nlims(j,1) ju = nlims(j,2) else pm(3) = ncells + nghosts jl = tlims(j,1) ju = tlims(j,2) end if if (n == 3) then pm(4) = nghosts kl = nlims(k,1) ku = nlims(k,2) else pm(4) = ncells + nghosts kl = tlims(k,1) ku = tlims(k,2) end if pdata%u(:,il:iu,jl:ju,kl:ku) = reshape(rbuf(:,l), pm) pinfo => pinfo%prev end do end if deallocate(sbuf, rbuf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do 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 blocks' edge ghost zones from blocks at ! the same level. ! ! Arguments: ! ! status - the call status; ! !=============================================================================== ! subroutine boundaries_edge_copy(status) 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 : bcells, ncells_half, nghosts, nb, ne #ifdef MPI use equations , only : nv #endif /* MPI */ use helpers , only : print_message #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(out) :: status 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 */ integer :: n integer :: i, il, iu, is, it integer :: j, jl, ju, js, jt integer :: k, kl, ku, ks, kt #ifdef MPI integer :: sproc, rproc integer :: ecount integer :: l, p integer, dimension(1) :: dm integer, dimension(4) :: pm real(kind=8), dimension(:,:), allocatable :: buf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims, dlims, tlims character(len=*), parameter :: loc = 'BOUNDARIES::boundaries_edge_copy()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = nb - 1 nlims(2,1) = ne + 1 nlims(2,2) = bcells dlims(1,1) = ne - nghosts + 1 dlims(1,2) = ne dlims(2,1) = nb dlims(2,2) = nb + nghosts - 1 tlims(1,1) = nb tlims(1,2) = nb + ncells_half - 1 tlims(2,1) = ne - ncells_half + 1 tlims(2,2) = ne first = .false. end if status = 0 #if NDIMS == 2 k = 1 kl = 1 ku = 1 ks = 1 kt = 1 #endif /* NDIMS == 2 */ #ifdef MPI sproc = 0 rproc = 0 ecount = 0 dm(1) = nv * nghosts**(NDIMS - 1) * ncells_half pm(1) = nv #if NDIMS == 2 pm(4) = 1 #endif /* NDIMS == 2 */ call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta pdata => pmeta%data do n = 1, NDIMS #if NDIMS == 3 do k = 1, nsides if (n == 3) then kl = tlims(k,1) ku = tlims(k,2) ks = tlims(k,1) kt = tlims(k,2) else kl = nlims(k,1) ku = nlims(k,2) ks = dlims(k,1) kt = dlims(k,2) end if #endif /* NDIMS == 3 */ do j = 1, nsides if (n == 2) then jl = tlims(j,1) ju = tlims(j,2) js = tlims(j,1) jt = tlims(j,2) else jl = nlims(j,1) ju = nlims(j,2) js = dlims(j,1) jt = dlims(j,2) end if do i = 1, nsides if (n == 1) then il = tlims(i,1) iu = tlims(i,2) is = tlims(i,1) it = tlims(i,2) else il = nlims(i,1) iu = nlims(i,2) is = dlims(i,1) it = dlims(i,2) end if #if NDIMS == 3 pneigh => pmeta%edges(i,j,k,n)%ptr #else /* NDIMS == 3 */ pneigh => pmeta%edges(i,j,n)%ptr #endif /* NDIMS == 3 */ if (associated(pneigh)) then if (pneigh%level == pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ pdata%u(:,il:iu,jl:ju,kl:ku) = & pneigh%data%u(:,is:it,js:jt,ks:kt) #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, n, [ i, j, k ]) end if #endif /* MPI */ end if end if end if end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ end do pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 if (ecount > 0) then allocate(buf(dm(1),ecount)) if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%neigh%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ if (n == 1) then is = tlims(i,1) it = tlims(i,2) else is = dlims(i,1) it = dlims(i,2) end if if (n == 2) then js = tlims(j,1) jt = tlims(j,2) else js = dlims(j,1) jt = dlims(j,2) end if #if NDIMS == 3 if (n == 3) then ks = tlims(k,1) kt = tlims(k,2) else ks = dlims(k,1) kt = dlims(k,2) end if #endif /* NDIMS == 3 */ buf(:,l) = reshape(pdata%u(:,is:it,js:jt,ks:kt), dm) pinfo => pinfo%prev end do call exchange_arrays(rproc, p, buf) l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%meta%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ if (n == 1) then pm(2) = ncells_half il = tlims(i,1) iu = tlims(i,2) else pm(2) = nghosts il = nlims(i,1) iu = nlims(i,2) end if if (n == 2) then pm(3) = ncells_half jl = tlims(j,1) ju = tlims(j,2) else pm(3) = nghosts jl = nlims(j,1) ju = nlims(j,2) end if #if NDIMS == 3 if (n == 3) then pm(4) = ncells_half kl = tlims(k,1) ku = tlims(k,2) else pm(4) = nghosts kl = nlims(k,1) ku = nlims(k,2) end if #endif /* NDIMS == 3 */ pdata%u(:,il:iu,jl:ju,kl:ku) = reshape(buf(:,l), pm) pinfo => pinfo%prev end do deallocate(buf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_edge_copy ! !=============================================================================== ! ! subroutine BOUNDARIES_EDGE_RESTRICT: ! ----------------------------------- ! ! Subroutine updates the blocks' edge ghost zones from blocks at ! higher levels. ! ! Arguments: ! ! status - the call status; ! !=============================================================================== ! subroutine boundaries_edge_restrict(status) 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 : bcells, ncells_half, nb, ne #ifdef MPI use coordinates, only : nghosts use equations , only : nv #endif /* MPI */ use helpers , only : print_message #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(out) :: status 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 */ integer :: n integer :: i, il, iu integer :: j, jl, ju integer :: k, kl, ku #ifdef MPI integer :: sproc, rproc integer :: scount, rcount integer :: l, p integer, dimension(1) :: dm integer, dimension(3) :: bm integer, dimension(4) :: pm real(kind=8), dimension(:,:,:,:), allocatable :: tmp real(kind=8), dimension(:,:) , allocatable :: sbuf, rbuf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims, tlims character(len=*), parameter :: loc = & 'BOUNDARIES::boundaries_edge_restrict()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = nb - 1 nlims(2,1) = ne + 1 nlims(2,2) = bcells tlims(1,1) = nb tlims(1,2) = nb + ncells_half - 1 tlims(2,1) = ne - ncells_half + 1 tlims(2,2) = ne first = .false. end if status = 0 #if NDIMS == 2 k = 1 kl = 1 ku = 1 #endif /* NDIMS == 2 */ #ifdef MPI sproc = 0 rproc = 0 dm(1) = nv * nghosts**(NDIMS - 1) * ncells_half pm(1) = nv #if NDIMS == 2 pm(4) = 1 bm(3) = 1 #endif /* NDIMS == 2 */ call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta pdata => pmeta%data do n = 1, NDIMS #if NDIMS == 3 do k = 1, nsides if (n == 3) then kl = tlims(k,1) ku = tlims(k,2) else kl = nlims(k,1) ku = nlims(k,2) end if #endif /* NDIMS == 3 */ do j = 1, nsides if (n == 2) then jl = tlims(j,1) ju = tlims(j,2) else jl = nlims(j,1) ju = nlims(j,2) end if do i = 1, nsides if (n == 1) then il = tlims(i,1) iu = tlims(i,2) else il = nlims(i,1) iu = nlims(i,2) end if #if NDIMS == 3 pneigh => pmeta%edges(i,j,k,n)%ptr #else /* NDIMS == 3 */ pneigh => pmeta%edges(i,j,n)%ptr #endif /* NDIMS == 3 */ if (associated(pneigh)) then if (pneigh%level > pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ call block_edge_restrict(n, [ i, j, k ], & pneigh%data%u, pdata%u(:,il:iu,jl:ju,kl:ku)) #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, n, [ i, j, k ]) end if #endif /* MPI */ end if end if end if end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ end do pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) if (scount > 0 .or. rcount > 0) then allocate(sbuf(dm(1),scount), rbuf(dm(1),rcount), stat=status) if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") if (scount > 0) then l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%neigh%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ bm(1:NDIMS) = nghosts bm(n) = ncells_half allocate(tmp(nv,bm(1),bm(2),bm(3)), stat=status) if (status /= 0) & call print_message(loc, & "Could not allocate the temporary buffer!") call block_edge_restrict(n, [ i, j, k ], pdata%u, tmp) sbuf(:,l) = reshape(tmp, dm) deallocate(tmp, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the temporary buffer!") pinfo => pinfo%prev end do end if call exchange_arrays(rproc, p, sbuf, rbuf) if (rcount > 0) then l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%meta%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ if (n == 1) then pm(2) = ncells_half il = tlims(i,1) iu = tlims(i,2) else pm(2) = nghosts il = nlims(i,1) iu = nlims(i,2) end if if (n == 2) then pm(3) = ncells_half jl = tlims(j,1) ju = tlims(j,2) else pm(3) = nghosts jl = nlims(j,1) ju = nlims(j,2) end if #if NDIMS == 3 if (n == 3) then pm(4) = ncells_half kl = tlims(k,1) ku = tlims(k,2) else pm(4) = nghosts kl = nlims(k,1) ku = nlims(k,2) end if #endif /* NDIMS == 3 */ pdata%u(:,il:iu,jl:ju,kl:ku) = reshape(rbuf(:,l), pm) pinfo => pinfo%prev end do end if deallocate(sbuf, rbuf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_edge_restrict ! !=============================================================================== ! ! subroutine BOUNDARIES_EDGE_PROLONG: ! ---------------------------------- ! ! Subroutine updates the blocks' edge ghost zones from blocks at lower levels. ! ! Arguments: ! ! level - the level to process; ! status - the call status; ! !=============================================================================== ! subroutine boundaries_edge_prolong(level, status) 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 : bcells, nb, ne #ifdef MPI use coordinates, only : ncells, nghosts use equations , only : nv #endif /* MPI */ use helpers , only : print_message #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(in) :: level integer, intent(out) :: status 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 */ integer :: n integer :: i, il, iu integer :: j, jl, ju integer :: k, kl, ku integer, dimension(3) :: bm #ifdef MPI integer :: sproc, rproc integer :: scount, rcount integer :: l, p integer, dimension(1) :: dm integer, dimension(4) :: pm real(kind=8), dimension(:,:,:,:), allocatable :: tmp real(kind=8), dimension(:,:) , allocatable :: sbuf, rbuf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims, tlims character(len=*), parameter :: loc = 'BOUNDARIES::boundaries_edge_prolong()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = nb - 1 nlims(2,1) = ne + 1 nlims(2,2) = bcells tlims(1,1) = nb tlims(1,2) = bcells tlims(2,1) = 1 tlims(2,2) = ne first = .false. end if status = 0 #if NDIMS == 2 k = 1 kl = 1 ku = 1 #endif /* NDIMS == 2 */ #ifdef MPI sproc = 0 rproc = 0 scount = 0 rcount = 0 dm(1) = nv * nghosts**(NDIMS - 1) * (ncells + nghosts) pm(1) = nv #if NDIMS == 2 pm(4) = 1 bm(3) = 1 #endif /* NDIMS == 2 */ call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta if (pmeta%level == level) then pdata => pmeta%data do n = 1, NDIMS #if NDIMS == 3 do k = 1, nsides if (n == 3) then bm(3) = pmeta%pos(3) + 1 kl = tlims(bm(3),1) ku = tlims(bm(3),2) else bm(3) = k kl = nlims(k,1) ku = nlims(k,2) end if #endif /* NDIMS == 3 */ do j = 1, nsides if (n == 2) then bm(2) = pmeta%pos(2) + 1 jl = tlims(bm(2),1) ju = tlims(bm(2),2) else bm(2) = j jl = nlims(j,1) ju = nlims(j,2) end if do i = 1, nsides if (n == 1) then bm(1) = pmeta%pos(1) + 1 il = tlims(bm(1),1) iu = tlims(bm(1),2) else bm(1) = i il = nlims(i,1) iu = nlims(i,2) end if #if NDIMS == 3 pneigh => pmeta%edges(i,j,k,n)%ptr #else /* NDIMS == 3 */ pneigh => pmeta%edges(i,j,n)%ptr #endif /* NDIMS == 3 */ if (associated(pneigh)) then if (pneigh%level < pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ call block_edge_prolong(n, bm, pneigh%data%u, & pdata%u(:,il:iu,jl:ju,kl:ku), status) #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, n, bm) end if #endif /* MPI */ end if end if end if end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ end do end if pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) if (scount > 0 .or. rcount > 0) then allocate(sbuf(dm(1),scount), rbuf(dm(1),rcount), stat=status) if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") if (scount > 0) then l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%neigh%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ bm(1:NDIMS) = nghosts bm(n) = ncells + nghosts allocate(tmp(nv,bm(1),bm(2),bm(3)), stat=status) if (status /= 0) & call print_message(loc, & "Could not allocate the temporary buffer!") call block_edge_prolong(n, [ i, j, k ], pdata%u, tmp, status) sbuf(:,l) = reshape(tmp, dm) deallocate(tmp, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the temporary buffer!") pinfo => pinfo%prev end do end if call exchange_arrays(rproc, p, sbuf, rbuf) if (rcount > 0) then l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%meta%data n = pinfo%direction i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ if (n == 1) then pm(2) = ncells + nghosts il = tlims(i,1) iu = tlims(i,2) else pm(2) = nghosts il = nlims(i,1) iu = nlims(i,2) end if if (n == 2) then pm(3) = ncells + nghosts jl = tlims(j,1) ju = tlims(j,2) else pm(3) = nghosts jl = nlims(j,1) ju = nlims(j,2) end if #if NDIMS == 3 if (n == 3) then pm(4) = ncells + nghosts kl = tlims(k,1) ku = tlims(k,2) else pm(4) = nghosts kl = nlims(k,1) ku = nlims(k,2) end if #endif /* NDIMS == 3 */ pdata%u(:,il:iu,jl:ju,kl:ku) = reshape(rbuf(:,l), pm) pinfo => pinfo%prev end do end if deallocate(sbuf, rbuf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_edge_prolong ! !=============================================================================== ! ! DOMAIN CORNER BOUNDARY UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BOUNDARIES_CORNER_COPY: ! --------------------------------- ! ! Subroutine updates the blocks' corner ghost zones from blocks ! at the same level. ! ! Arguments: ! ! status - the call status; ! !=============================================================================== ! subroutine boundaries_corner_copy(status) 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 : nn => bcells, ng => nghosts, nb, ne #ifdef MPI use equations , only : nv use helpers , only : print_message use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(out) :: status 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 */ integer :: i, il, iu, is, it integer :: j, jl, ju, js, jt integer :: k, kl, ku, ks, kt #ifdef MPI integer :: sproc, rproc integer :: ecount integer :: l, p real(kind=8), dimension(:,:,:,:,:), allocatable :: buf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: dlims, blims #ifdef MPI character(len=*), parameter :: loc = 'BOUNDARIES::boundaries_corner_copy()' #endif /* MPI */ !------------------------------------------------------------------------------- ! if (first) then dlims(1,1) = ne - ng + 1 dlims(1,2) = ne dlims(2,1) = nb dlims(2,2) = nb + ng - 1 blims(1,1) = 1 blims(1,2) = ng blims(2,1) = nn - ng + 1 blims(2,2) = nn first = .false. end if status = 0 #if NDIMS == 2 k = 1 kl = 1 ku = 1 ks = 1 kt = 1 #endif /* NDIMS == 2 */ #ifdef MPI sproc = 0 rproc = 0 ecount = 0 call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta pdata => pmeta%data #if NDIMS == 3 do k = 1, nsides kl = blims(k,1) ku = blims(k,2) ks = dlims(k,1) kt = dlims(k,2) #endif /* NDIMS == 3 */ do j = 1, nsides jl = blims(j,1) ju = blims(j,2) js = dlims(j,1) jt = dlims(j,2) do i = 1, nsides il = blims(i,1) iu = blims(i,2) is = dlims(i,1) it = dlims(i,2) #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 */ if (associated(pneigh)) then if (pneigh%level == pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ pdata%u(:,il:iu,jl:ju,kl:ku) = & pneigh%data%u(:,is:it,js:jt,ks:kt) #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, -1, [ i, j, k ]) end if #endif /* MPI */ end if end if end if end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 if (ecount > 0) then #if NDIMS == 3 allocate(buf(nv,ng,ng,ng,ecount), stat=status) #else /* NDIMS == 3 */ allocate(buf(nv,ng,ng, 1,ecount), stat=status) #endif /* NDIMS == 3 */ if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pneigh => pinfo%neigh is = dlims(pinfo%location(1),1) it = dlims(pinfo%location(1),2) js = dlims(pinfo%location(2),1) jt = dlims(pinfo%location(2),2) #if NDIMS == 3 ks = dlims(pinfo%location(3),1) kt = dlims(pinfo%location(3),2) #endif /* NDIMS == 3 */ buf(:,:,:,:,l) = pneigh%data%u(:,is:it,js:jt,ks:kt) pinfo => pinfo%prev end do call exchange_arrays(rproc, p, buf) l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pmeta => pinfo%meta il = blims(pinfo%location(1),1) iu = blims(pinfo%location(1),2) jl = blims(pinfo%location(2),1) ju = blims(pinfo%location(2),2) #if NDIMS == 3 kl = blims(pinfo%location(3),1) ku = blims(pinfo%location(3),2) #endif /* NDIMS == 3 */ pmeta%data%u(:,il:iu,jl:ju,kl:ku) = buf(:,:,:,:,l) pinfo => pinfo%prev end do deallocate(buf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do ! p = 1, npairs call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_corner_copy ! !=============================================================================== ! ! subroutine BOUNDARIES_CORNER_RESTRICT: ! ------------------------------------- ! ! Subroutine updates the blocks' corner ghost zones from blocks ! at higher levels. ! ! Arguments: ! ! status - the call status; ! !=============================================================================== ! subroutine boundaries_corner_restrict(status) 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 : nn => bcells, ng => nghosts #ifdef MPI use equations , only : nv use helpers , only : print_message use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(out) :: status 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 */ integer :: i, il, iu integer :: j, jl, ju integer :: k, kl, ku #ifdef MPI integer :: sproc, rproc integer :: scount, rcount integer :: l, p real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims #ifdef MPI character(len=*), parameter :: loc = & 'BOUNDARIES::boundaries_corner_restrict()' #endif /* MPI */ !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = ng nlims(2,1) = nn - ng + 1 nlims(2,2) = nn first = .false. end if status = 0 #if NDIMS == 2 k = 1 kl = 1 ku = 1 #endif /* NDIMS == 2 */ #ifdef MPI sproc = 0 rproc = 0 scount = 0 rcount = 0 call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta pdata => pmeta%data #if NDIMS == 3 do k = 1, nsides kl = nlims(k,1) ku = nlims(k,2) #endif /* NDIMS == 3 */ do j = 1, nsides jl = nlims(j,1) ju = nlims(j,2) do i = 1, nsides il = nlims(i,1) iu = nlims(i,2) #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 */ if (associated(pneigh)) then if (pneigh%level > pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ call block_corner_restrict([ i, j, k ], pneigh%data%u, & pdata%u(:,il:iu,jl:ju,kl:ku)) #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, -1, [ i, j, k ]) end if #endif /* MPI */ end if end if end if end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) if (scount > 0 .or. rcount > 0) then #if NDIMS == 3 allocate(sbuf(nv,ng,ng,ng,scount), & rbuf(nv,ng,ng,ng,rcount), stat=status) #else /* NDIMS == 3 */ allocate(sbuf(nv,ng,ng, 1,scount), & rbuf(nv,ng,ng, 1,rcount), stat=status) #endif /* NDIMS == 3 */ if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") if (scount > 0) then l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pneigh => pinfo%neigh i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ call block_corner_restrict([ i, j, k ], pneigh%data%u, & sbuf(:,:,:,:,l)) pinfo => pinfo%prev end do end if call exchange_arrays(rproc, p, sbuf, rbuf) if (rcount > 0) then l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pmeta => pinfo%meta il = nlims(pinfo%location(1),1) iu = nlims(pinfo%location(1),2) jl = nlims(pinfo%location(2),1) ju = nlims(pinfo%location(2),2) #if NDIMS == 3 kl = nlims(pinfo%location(3),1) ku = nlims(pinfo%location(3),2) #endif /* NDIMS == 3 */ pmeta%data%u(:,il:iu,jl:ju,kl:ku) = rbuf(:,:,:,:,l) pinfo => pinfo%prev end do end if deallocate(sbuf, rbuf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do ! p = 1, npairs call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_corner_restrict ! !=============================================================================== ! ! subroutine BOUNDARIES_CORNER_PROLONG: ! ------------------------------------ ! ! Subroutine updates the blocks' corner ghost zones from blocks ! at lower levels. ! ! Arguments: ! ! level - the level to process; ! status - the call status; ! !=============================================================================== ! subroutine boundaries_corner_prolong(level, status) 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 : nn => bcells, ng => nghosts #ifdef MPI use equations , only : nv #endif /* MPI */ use helpers , only : print_message #ifdef MPI use mpitools , only : nproc, npairs, pairs use mpitools , only : exchange_arrays #endif /* MPI */ implicit none integer, intent(in) :: level integer, intent(out) :: status 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 */ integer :: i, il, iu integer :: j, jl, ju integer :: k, kl, ku #ifdef MPI integer :: sproc, rproc integer :: scount, rcount integer :: l, p real(kind=8), dimension(:,:,:,:,:), allocatable :: sbuf, rbuf #endif /* MPI */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims character(len=*), parameter :: loc = & 'BOUNDARIES::boundaries_corner_prolong()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = ng nlims(2,1) = nn - ng + 1 nlims(2,2) = nn first = .false. end if status = 0 #if NDIMS == 2 k = 1 kl = 1 ku = 1 #endif /* NDIMS == 2 */ #ifdef MPI sproc = 0 rproc = 0 scount = 0 rcount = 0 call prepare_exchange_array() #endif /* MPI */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta if (pmeta%level == level) then pdata => pmeta%data #if NDIMS == 3 do k = 1, nsides kl = nlims(k,1) ku = nlims(k,2) #endif /* NDIMS == 3 */ do j = 1, nsides jl = nlims(j,1) ju = nlims(j,2) do i = 1, nsides il = nlims(i,1) iu = nlims(i,2) #if NDIMS == 3 pneigh => pmeta%corners(i,j,k)%ptr #else /* NDIMS == 3 */ pneigh => pmeta%corners(i,j)%ptr #endif /* NDIMS == 3 */ if (associated(pneigh)) then if (pneigh%level < pmeta%level) then if (pmeta%update .or. pneigh%update) then pmeta%boundary = .true. #ifdef MPI if (pmeta%process == pneigh%process) then if (pneigh%process == nproc) then #endif /* MPI */ call block_corner_prolong([ i, j, k ], pneigh%data%u, & pdata%u(:,il:iu,jl:ju,kl:ku), status) if (status /= 0) & call print_message(loc, & "Block's corner prolongation failed!") #ifdef MPI end if else call append_exchange_block(pmeta, pneigh, -1, [ i, j, k ]) end if #endif /* MPI */ end if end if end if end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ end if pleaf => pleaf%next end do #ifdef MPI do p = 1, npairs if (pairs(p,1) == nproc .or. pairs(p,2) == nproc) then 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 scount = bcount(sproc,rproc) rcount = bcount(rproc,sproc) if (scount > 0 .or. rcount > 0) then #if NDIMS == 3 allocate(sbuf(nv,ng,ng,ng,scount), & rbuf(nv,ng,ng,ng,rcount), stat=status) #else /* NDIMS == 3 */ allocate(sbuf(nv,ng,ng, 1,scount), & rbuf(nv,ng,ng, 1,rcount), stat=status) #endif /* NDIMS == 3 */ if (status /= 0) & call print_message(loc, "Could not allocate the exchange buffers!") if (scount > 0) then l = 0 pinfo => barray(sproc,rproc)%ptr do while(associated(pinfo)) l = l + 1 pneigh => pinfo%neigh i = pinfo%location(1) j = pinfo%location(2) #if NDIMS == 3 k = pinfo%location(3) #endif /* NDIMS == 3 */ call block_corner_prolong([ i, j, k ], pneigh%data%u, & sbuf(:,:,:,:,l), status) pinfo => pinfo%prev end do end if call exchange_arrays(rproc, p, sbuf, rbuf) if (rcount > 0) then l = 0 pinfo => barray(rproc,sproc)%ptr do while(associated(pinfo)) l = l + 1 pdata => pinfo%meta%data il = nlims(pinfo%location(1),1) iu = nlims(pinfo%location(1),2) jl = nlims(pinfo%location(2),1) ju = nlims(pinfo%location(2),2) #if NDIMS == 3 kl = nlims(pinfo%location(3),1) ku = nlims(pinfo%location(3),2) #endif /* NDIMS == 3 */ pdata%u(:,il:iu,jl:ju,kl:ku) = rbuf(:,:,:,:,l) pinfo => pinfo%prev end do end if deallocate(sbuf, rbuf, stat=status) if (status /= 0) & call print_message(loc, & "Could not deallocate the exchange buffers!") end if end if end do ! p = 1, npairs call release_exchange_array() #endif /* MPI */ !------------------------------------------------------------------------------- ! end subroutine boundaries_corner_prolong ! !=============================================================================== ! ! 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; ! status - the call status; ! !=============================================================================== ! subroutine boundaries_specific(t, dt, status) 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 */ implicit none real(kind=8), intent(in) :: t, dt integer , intent(out) :: status type(block_meta), pointer :: pmeta type(block_leaf), pointer :: pleaf integer :: i, j, k, n #if NDIMS == 2 integer :: m #endif /* NDIMS == 2 */ integer, dimension(NDIMS) :: s real(kind=8), dimension(nn) :: x, y #if NDIMS == 3 real(kind=8), dimension(nn) :: z #else /* NDIMS == 3 */ real(kind=8), dimension( 1) :: z #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! #if NDIMS == 2 k = 1 z = 0.0d+00 #endif /* NDIMS == 2 */ pleaf => list_leaf do while(associated(pleaf)) pmeta => pleaf%meta if (pmeta%update) then #ifdef MPI if (pmeta%process == nproc) then #endif /* MPI */ x(:) = pmeta%xmin + ax(pmeta%level,:) y(:) = pmeta%ymin + ay(pmeta%level,:) #if NDIMS == 3 z(:) = pmeta%zmin + az(pmeta%level,:) #endif /* NDIMS == 3 */ #if NDIMS == 2 do n = 1, ndims 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 do j = 1, nsides s(2) = j do i = 1, nsides s(1) = i if (.not. associated(pmeta%edges(i,j,m)%ptr)) then pmeta%boundary = .true. call block_primitive_variables(n, pmeta%data, status) call block_boundary_specific(n, [ i, j, k ], t, dt, & x, y, z, pmeta%data%q, status) call block_conservative_variables(n, s(n), pmeta%data) end if end do end do end if end do ! n = 1, ndims #endif /* NDIMS == 2 */ #if NDIMS == 3 do n = 1, ndims if (.not. periodic(n)) then do k = 1, nsides s(3) = k do j = 1, nsides s(2) = j do i = 1, nsides s(1) = i if (.not. associated(pmeta%faces(i,j,k,n)%ptr)) then pmeta%boundary = .true. call block_primitive_variables(n, pmeta%data, status) call block_boundary_specific(n, [ i, j, k ], t, dt, & x, y, z, pmeta%data%q, status) call block_conservative_variables(n, s(n), pmeta%data) end if end do end do end do end if end do #endif /* NDIMS == 3 */ #ifdef MPI end if #endif /* MPI */ end if pleaf => pleaf%next end do !------------------------------------------------------------------------------- ! end subroutine boundaries_specific #if NDIMS == 3 ! !=============================================================================== ! ! BLOCK FACE UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_FACE_RESTRICT: ! ------------------------------ ! ! Subroutine restricts the region of qn specified by the face position ! and inserts it in the corresponding face boundary region of qb. ! ! Arguments: ! ! dir - the face direction; ! pos - the face position; ! qn - the array to restrict from; ! qb - the array of the face boundary region; ! !=============================================================================== ! subroutine block_face_restrict(dir, pos, qn, qb) use coordinates, only : nghosts_double, nb, ne implicit none integer , intent(in) :: dir integer , dimension(3) , intent(in) :: pos real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(inout) :: qb integer :: il, iu, im, ip integer :: jl, ju, jm, jp integer :: kl, ku, km, kp logical , save :: first = .true. integer, dimension(2,2), save :: nlims !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = ne - nghosts_double + 1 nlims(1,2) = ne nlims(2,1) = nb nlims(2,2) = nb + nghosts_double - 1 first = .false. end if select case(dir) case(1) il = nlims(pos(1),1) iu = nlims(pos(1),2) jl = nb ju = ne kl = nb ku = ne case(2) il = nb iu = ne jl = nlims(pos(2),1) ju = nlims(pos(2),2) kl = nb ku = ne case default il = nb iu = ne jl = nb ju = ne kl = nlims(pos(3),1) ku = nlims(pos(3),2) end select ip = il + 1 im = iu - 1 jp = jl + 1 jm = ju - 1 kp = kl + 1 km = ku - 1 qb(:,:,:,:) = 1.25d-01 * (((qn(:,il:im:2,jl:jm:2,kl:km:2) & + qn(:,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(:,il:im:2,jl:jm:2,kp:ku:2) & + qn(:,ip:iu:2,jp:ju:2,kl:km:2))) & + ((qn(:,il:im:2,jp:ju:2,kp:ku:2) & + qn(:,ip:iu:2,jl:jm:2,kl:km:2)) & + (qn(:,il:im:2,jp:ju:2,kl:km:2) & + qn(:,ip:iu:2,jl:jm:2,kp:ku:2)))) !------------------------------------------------------------------------------- ! end subroutine block_face_restrict ! !=============================================================================== ! ! subroutine BLOCK_FACE_PROLONG: ! ----------------------------- ! ! Subroutine prolongs the region of qn specified by the face position ! and inserts it in the corresponding face boundary region of qb. ! ! Arguments: ! ! dir - the face direction; ! pos - the face position; ! qn - the array to prolong from; ! qb - the array of the face boundary region; ! status - the call status; ! !=============================================================================== ! subroutine block_face_prolong(dir, pos, qn, qb, status) use coordinates , only : nb, ne, ncells_half, nghosts_half 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(inout) :: qb integer , intent(out) :: status integer :: p integer :: i, il, iu, is, it, im, ip integer :: j, jl, ju, js, jt, jm, jp integer :: k, kl, ku, ks, kt, km, kp real(kind=8) :: dq1, dq2, dq3, dq4 real(kind=8), dimension(3) :: dq logical , save :: first = .true. integer, dimension(2,2), save :: nlims, tlims character(len=*), parameter :: loc = 'BOUNDARIES::block_face_prolong()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = ne - nghosts_half + 1 nlims(1,2) = ne nlims(2,1) = nb nlims(2,2) = nb + nghosts_half - 1 tlims(1,1) = nb tlims(1,2) = nb + ncells_half + nghosts_half - 1 tlims(2,1) = ne - ncells_half - nghosts_half + 1 tlims(2,2) = ne first = .false. end if status = 0 select case(dir) case(1) il = nlims(pos(1),1) iu = nlims(pos(1),2) jl = tlims(pos(2),1) ju = tlims(pos(2),2) kl = tlims(pos(3),1) ku = tlims(pos(3),2) case(2) il = tlims(pos(1),1) iu = tlims(pos(1),2) jl = nlims(pos(2),1) ju = nlims(pos(2),2) kl = tlims(pos(3),1) ku = tlims(pos(3),2) case default il = tlims(pos(1),1) iu = tlims(pos(1),2) jl = tlims(pos(2),1) ju = tlims(pos(2),2) kl = nlims(pos(3),1) ku = nlims(pos(3),2) end select do k = kl, ku km = k - 1 kp = k + 1 ks = 2 * (k - kl) + 1 kt = ks + 1 do j = jl, ju jm = j - 1 jp = j + 1 js = 2 * (j - jl) + 1 jt = js + 1 do i = il, iu im = i - 1 ip = i + 1 is = 2 * (i - il) + 1 it = is + 1 do p = 1, nv dq1 = qn(p,i ,j,k) - qn(p,im,j,k) dq2 = qn(p,ip,j,k) - qn(p,i ,j,k) dq(1) = limiter_prol(2.5d-01, dq1, dq2) dq1 = qn(p,i,j ,k) - qn(p,i,jm,k) dq2 = qn(p,i,jp,k) - qn(p,i,j ,k) dq(2) = limiter_prol(2.5d-01, dq1, dq2) dq1 = qn(p,i,j,k ) - qn(p,i,j,km) dq2 = qn(p,i,j,kp) - qn(p,i,j,k ) dq(3) = limiter_prol(2.5d-01, dq1, dq2) if (positive(p)) then if (qn(p,i,j,k) > 0.0d+00) then dq1 = sum(abs(dq(1:NDIMS))) if (qn(p,i,j,k) <= dq1) & dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:) else call print_message(loc, "Positive variable is not positive!") status = 1 return 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 end do end do end do !------------------------------------------------------------------------------- ! end subroutine block_face_prolong #endif /* NDIMS == 3 */ ! !=============================================================================== ! ! BLOCK EDGE UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_EDGE_RESTRICT: ! ------------------------------ ! ! Subroutine restricts the region of qn specified by the edge position ! and inserts it in the corresponding edge boundary region of qb. ! ! Arguments: ! ! dir - the edge direction; ! pos - the edge position; ! qn - the array to restrict from; ! qb - the array of the edge boundary region; ! !=============================================================================== ! subroutine block_edge_restrict(dir, pos, qn, qb) use coordinates, only : nb, ne, nghosts_double implicit none integer , intent(in) :: dir integer , dimension(3) , intent(in) :: pos real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(inout) :: qb integer :: il, iu, im, ip integer :: jl, ju, jm, jp #if NDIMS == 3 integer :: kl, ku, km, kp #endif /* NDIMS == 3 */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = ne - nghosts_double + 1 nlims(1,2) = ne nlims(2,1) = nb nlims(2,2) = nb + nghosts_double - 1 first = .false. end if il = nlims(pos(1),1) iu = nlims(pos(1),2) jl = nlims(pos(2),1) ju = nlims(pos(2),2) #if NDIMS == 3 kl = nlims(pos(3),1) ku = nlims(pos(3),2) #endif /* NDIMS == 3 */ select case(dir) case(1) il = nb iu = ne case(2) jl = nb ju = ne #if NDIMS == 3 case(3) kl = nb ku = ne #endif /* NDIMS == 3 */ end select ip = il + 1 im = iu - 1 jp = jl + 1 jm = ju - 1 #if NDIMS == 3 kp = kl + 1 km = ku - 1 qb(:,:,:,:) = 1.25d-01 * (((qn(:,il:im:2,jl:jm:2,kl:km:2) & + qn(:,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(:,il:im:2,jl:jm:2,kp:ku:2) & + qn(:,ip:iu:2,jp:ju:2,kl:km:2))) & + ((qn(:,il:im:2,jp:ju:2,kp:ku:2) & + qn(:,ip:iu:2,jl:jm:2,kl:km:2)) & + (qn(:,il:im:2,jp:ju:2,kl:km:2) & + qn(:,ip:iu:2,jl:jm:2,kp:ku:2)))) #else /* NDIMS == 3 */ qb(:,:,:,:) = 2.50d-01 * ((qn(:,il:im:2,jl:jm:2,:) & + qn(:,ip:iu:2,jp:ju:2,:)) & + (qn(:,il:im:2,jp:ju:2,:) & + qn(:,ip:iu:2,jl:jm:2,:))) #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine block_edge_restrict ! !=============================================================================== ! ! subroutine BLOCK_EDGE_PROLONG: ! ----------------------------- ! ! Subroutine prolongss the region of qn specified by the edge position ! and inserts it in the corresponding edge boundary region of qb. ! ! Arguments: ! ! dir - the edge direction; ! pos - the edge position; ! qn - the array to prolong from; ! qb - the array of the edge boundary region; ! status - the call status; ! !=============================================================================== ! subroutine block_edge_prolong(dir, pos, qn, qb, status) use coordinates , only : nb, ne, ncells_half, nghosts_half 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(inout) :: qb integer , intent(out) :: status integer :: p integer :: i, il, iu, is, it, im, ip integer :: j, jl, ju, js, jt, jm, jp #if NDIMS == 3 integer :: k, kl, ku, ks, kt, km, kp #else /* NDIMS == 3 */ integer :: k #endif /* NDIMS == 3 */ real(kind=8) :: dq1, dq2 #if NDIMS == 3 real(kind=8) :: dq3, dq4 #endif /* NDIMS == 3 */ real(kind=8), dimension(NDIMS) :: dq logical , save :: first = .true. integer, dimension(2,2), save :: nlims, tlims character(len=*), parameter :: loc = 'BOUNDARIES::block_edge_prolong()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = ne - nghosts_half + 1 nlims(1,2) = ne nlims(2,1) = nb nlims(2,2) = nb + nghosts_half - 1 tlims(1,1) = nb tlims(1,2) = nb + ncells_half + nghosts_half - 1 tlims(2,1) = ne - ncells_half - nghosts_half + 1 tlims(2,2) = ne first = .false. end if status = 0 il = nlims(pos(1),1) iu = nlims(pos(1),2) jl = nlims(pos(2),1) ju = nlims(pos(2),2) #if NDIMS == 3 kl = nlims(pos(3),1) ku = nlims(pos(3),2) #endif /* NDIMS == 3 */ select case(dir) case(1) il = tlims(pos(1),1) iu = tlims(pos(1),2) case(2) jl = tlims(pos(2),1) ju = tlims(pos(2),2) #if NDIMS == 3 case(3) kl = tlims(pos(3),1) ku = tlims(pos(3),2) #endif /* NDIMS == 3 */ end select #if NDIMS == 3 do k = kl, ku km = k - 1 kp = k + 1 ks = 2 * (k - kl) + 1 kt = ks + 1 #else /* NDIMS == 3 */ k = 1 #endif /* NDIMS == 3 */ do j = jl, ju jm = j - 1 jp = j + 1 js = 2 * (j - jl) + 1 jt = js + 1 do i = il, iu im = i - 1 ip = i + 1 is = 2 * (i - il) + 1 it = is + 1 do p = 1, nv dq1 = qn(p,i ,j,k) - qn(p,im,j,k) dq2 = qn(p,ip,j,k) - qn(p,i ,j,k) dq(1) = limiter_prol(2.5d-01, dq1, dq2) dq1 = qn(p,i,j ,k) - qn(p,i,jm,k) dq2 = qn(p,i,jp,k) - qn(p,i,j ,k) dq(2) = limiter_prol(2.5d-01, dq1, dq2) #if NDIMS == 3 dq1 = qn(p,i,j,k ) - qn(p,i,j,km) dq2 = qn(p,i,j,kp) - qn(p,i,j,k ) dq(3) = limiter_prol(2.5d-01, dq1, dq2) #endif /* NDIMS == 3 */ if (positive(p)) then if (qn(p,i,j,k) > 0.0d+00) then dq1 = sum(abs(dq(1:NDIMS))) if (qn(p,i,j,k) <= dq1) & dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:) else call print_message(loc, "Positive variable is not positive!") status = 1 return 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 end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine block_edge_prolong ! !=============================================================================== ! ! BLOCK CORNER UPDATE SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_CORNER_RESTRICT: ! -------------------------------- ! ! Subroutine restricts the region of qn specified by the corner position ! and inserts it in the corresponding corner boundary region of qb. ! ! Arguments: ! ! pos - the corner position; ! qn - the array to restrict from; ! qb - the array of the corner boundary region; ! !=============================================================================== ! subroutine block_corner_restrict(pos, qn, qb) use coordinates, only : nb, ne, nghosts_double implicit none integer , dimension(3) , intent(in) :: pos real(kind=8), dimension(:,:,:,:), intent(in) :: qn real(kind=8), dimension(:,:,:,:), intent(inout) :: qb integer :: il, iu, ip, im integer :: jl, ju, jp, jm #if NDIMS == 3 integer :: kl, ku, kp, km #endif /* NDIMS == 3 */ logical , save :: first = .true. integer, dimension(2,2), save :: nlims !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = ne - nghosts_double + 1 nlims(1,2) = ne nlims(2,1) = nb nlims(2,2) = nb + nghosts_double - 1 first = .false. end if il = nlims(pos(1),1) iu = nlims(pos(1),2) jl = nlims(pos(2),1) ju = nlims(pos(2),2) #if NDIMS == 3 kl = nlims(pos(3),1) ku = nlims(pos(3),2) #endif /* NDIMS == 3 */ ip = il + 1 jp = jl + 1 #if NDIMS == 3 kp = kl + 1 #endif /* NDIMS == 3 */ im = iu - 1 jm = ju - 1 #if NDIMS == 3 km = ku - 1 #endif /* NDIMS == 3 */ #if NDIMS == 3 qb(:,:,:,:) = 1.25d-01 * (((qn(:,il:im:2,jl:jm:2,kl:km:2) & + qn(:,ip:iu:2,jp:ju:2,kp:ku:2)) & + (qn(:,il:im:2,jl:jm:2,kp:ku:2) & + qn(:,ip:iu:2,jp:ju:2,kl:km:2))) & + ((qn(:,il:im:2,jp:ju:2,kp:ku:2) & + qn(:,ip:iu:2,jl:jm:2,kl:km:2)) & + (qn(:,il:im:2,jp:ju:2,kl:km:2) & + qn(:,ip:iu:2,jl:jm:2,kp:ku:2)))) #else /* NDIMS == 3 */ qb(:,:,:,:) = 2.50d-01 * ((qn(:,il:im:2,jl:jm:2,:) & + qn(:,ip:iu:2,jp:ju:2,:)) & + (qn(:,il:im:2,jp:ju:2,:) & + qn(:,ip:iu:2,jl:jm:2,:))) #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine block_corner_restrict ! !=============================================================================== ! ! subroutine BLOCK_CORNER_PROLONG: ! ------------------------------- ! ! Subroutine prolongs the region of qn specified by the corner position ! and inserts it in the corresponding corner boundary region of qb. ! ! Remarks: ! ! This subroutine requires the ghost zone regions of the qn already ! updated in order to perform a consistent interpolation. ! ! Arguments: ! ! pos - the corner position; ! qn - the array to prolong from; ! qb - the array of the corner boundary region; ! status - the call status; ! !=============================================================================== ! subroutine block_corner_prolong(pos, qn, qb, status) use coordinates , only : nb, ne, nghosts_half 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(inout) :: qb integer , intent(out) :: status integer :: p integer :: i, il, iu, im, ip, is, it integer :: j, jl, ju, jm, jp, js, jt #if NDIMS == 3 integer :: k, kl, ku, km, kp, ks, kt #else /* NDIMS == 3 */ integer :: k #endif /* NDIMS == 3 */ real(kind=8) :: dq1, dq2 #if NDIMS == 3 real(kind=8) :: dq3, dq4 #endif /* NDIMS == 3 */ real(kind=8), dimension(NDIMS) :: dq logical , save :: first = .true. integer, dimension(2,2), save :: nlims character(len=*), parameter :: loc = 'BOUNDARIES::block_corner_prolong()' !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = ne - nghosts_half + 1 nlims(1,2) = ne nlims(2,1) = nb nlims(2,2) = nb + nghosts_half - 1 first = .false. end if status = 0 il = nlims(pos(1),1) iu = nlims(pos(1),2) jl = nlims(pos(2),1) ju = nlims(pos(2),2) #if NDIMS == 3 kl = nlims(pos(3),1) ku = nlims(pos(3),2) do k = kl, ku km = k - 1 kp = k + 1 ks = 2 * (k - kl) + 1 kt = ks + 1 #else /* NDIMS == 3 */ k = 1 #endif /* NDIMS == 3 */ do j = jl, ju jm = j - 1 jp = j + 1 js = 2 * (j - jl) + 1 jt = js + 1 do i = il, iu im = i - 1 ip = i + 1 is = 2 * (i - il) + 1 it = is + 1 do p = 1, nv dq1 = qn(p,i ,j,k) - qn(p,im,j,k) dq2 = qn(p,ip,j,k) - qn(p,i ,j,k) dq(1) = limiter_prol(2.5d-01, dq1, dq2) dq1 = qn(p,i,j ,k) - qn(p,i,jm,k) dq2 = qn(p,i,jp,k) - qn(p,i,j ,k) dq(2) = limiter_prol(2.5d-01, dq1, dq2) #if NDIMS == 3 dq1 = qn(p,i,j,k ) - qn(p,i,j,km) dq2 = qn(p,i,j,kp) - qn(p,i,j,k ) dq(3) = limiter_prol(2.5d-01, dq1, dq2) #endif /* NDIMS == 3 */ if (positive(p)) then if (qn(p,i,j,k) > 0.0d+00) then dq1 = sum(abs(dq(1:NDIMS))) if (qn(p,i,j,k) <= dq1) & dq(:) = (5.0d-01 * qn(p,i,j,k) / dq1) * dq(:) else call print_message(loc, "Positive variable is not positive!") status = 1 return 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 end do end do #if NDIMS == 3 end do #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine block_corner_prolong ! !=============================================================================== ! ! BLOCK SPECIFIC BOUNDARY SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_BOUNDARY_SPECIFIC: ! ---------------------------------- ! ! Subroutine applies specific boundary conditions to the pointed data block. ! ! Arguments: ! ! dir - the direction; ! side - the side of the boundary; ! t, dt - time and time increment; ! x, y, z - the block coordinates; ! qn - the variable array; ! status - the call status; ! !=============================================================================== ! subroutine block_boundary_specific(dir, side, t, dt, x, y, z, qn, status) use coordinates , only : nn => bcells, nh => bcells_half, ng => nghosts use coordinates , only : nb, ne, nbl, neu use equations , only : magnetized 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 implicit none integer , intent(in) :: dir integer , dimension(3) , intent(in) :: side real(kind=8) , intent(in) :: t, dt real(kind=8), dimension(:) , intent(inout) :: x, y, z real(kind=8), dimension(:,:,:,:), intent(inout) :: qn integer , intent(out) :: status integer :: i, il, iu, is, it, im, ip integer :: j, jl, ju, js, jt, jm, jp integer :: k, kl, ku #if NDIMS == 3 integer :: ks, kt, km, kp 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()' !------------------------------------------------------------------------------- ! status = 0 select case(dir) case(1) 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 */ select case(bnd_type(dir,side(1))) case(bnd_outflow) if (side(1) == 1) then do i = nbl, 1, -1 qn( : ,i,jl:ju,kl:ku) = qn( : ,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 else do i = neu, nn qn( : ,i,jl:ju,kl:ku) = qn( : ,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 end if case(bnd_reflective) if (side(1) == 1) then do i = 1, ng it = nb - i is = nbl + i qn( : ,it,jl:ju,kl:ku) = qn( : ,is,jl:ju,kl:ku) qn(ivx,it,jl:ju,kl:ku) = - qn(ivx,is,jl:ju,kl:ku) if (magnetized) 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( : ,it,jl:ju,kl:ku) = qn( : ,is,jl:ju,kl:ku) qn(ivx,it,jl:ju,kl:ku) = - qn(ivx,is,jl:ju,kl:ku) if (magnetized) then qn(ibx,it,jl:ju,kl:ku) = - qn(ibx,is,jl:ju,kl:ku) end if end do end if 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 ip = i + 1 xi = x(i) + dxh do k = kl, ku do j = jl, ju qn(:,i,j,k) = qn(:,nb,j,k) call gravitational_acceleration(t, dt, xi, y(j), z(k), ga(:)) qn(ipr,i,j,k) = qn(ipr,ip,j,k) & - (qn(idn,ip,j,k) + qn(idn,i,j,k)) * ga(1) * dxh end do end do end do else do i = neu, nn im = i - 1 xi = x(i) - dxh do k = kl, ku do j = jl, ju qn(:,i,j,k) = qn(:,ne,j,k) call gravitational_acceleration(t, dt, xi, y(j), z(k), ga(:)) qn(ipr,i,j,k) = qn(ipr,im,j,k) & + (qn(idn,im,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 ip = i + 1 xi = x(i) + dxh do k = kl, ku do j = jl, ju qn(:,i,j,k) = qn(:,nb,j,k) call gravitational_acceleration(t, dt, xi, y(j), z(k), ga(:)) qn(idn,i,j,k) = qn(idn,ip,j,k) * exp(- ga(1) * dx / csnd2) end do end do end do else do i = neu, nn im = i - 1 xi = x(i) - dxh do k = kl, ku do j = jl, ju qn(:,i,j,k) = qn(:,ne,j,k) call gravitational_acceleration(t, dt, xi, y(j), z(k), ga(:)) qn(idn,i,j,k) = qn(idn,im,j,k) * exp( ga(1) * dx / csnd2) end do end do end do end if end if case(bnd_custom) if (associated(custom_boundary_x)) then call custom_boundary_x(side(1), jl, ju, kl, ku, & t, dt, x(:), y(:), z(:), qn(:,:,:,:)) else call print_message(loc, & "No subroutine associated with custom_boundary_x()!") status = 1 return end if case default if (side(1) == 1) then do i = nbl, 1, -1 qn(:,i,jl:ju,kl:ku) = qn(:,nb,jl:ju,kl:ku) end do else do i = neu, nn qn(:,i,jl:ju,kl:ku) = qn(:,ne,jl:ju,kl:ku) end do end if end select case(2) 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 */ select case(bnd_type(dir,side(2))) case(bnd_outflow) if (side(2) == 1) then do j = nbl, 1, -1 qn( : ,il:iu,j,kl:ku) = qn( : ,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 else do j = neu, nn qn( : ,il:iu,j,kl:ku) = qn( : ,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 end if case(bnd_reflective) if (side(2) == 1) then do j = 1, ng jt = nb - j js = nbl + j qn( : ,il:iu,jt,kl:ku) = qn( : ,il:iu,js,kl:ku) qn(ivy,il:iu,jt,kl:ku) = - qn(ivy,il:iu,js,kl:ku) if (magnetized) 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( : ,il:iu,jt,kl:ku) = qn( : ,il:iu,js,kl:ku) qn(ivy,il:iu,jt,kl:ku) = - qn(ivy,il:iu,js,kl:ku) if (magnetized) then qn(iby,il:iu,jt,kl:ku) = - qn(iby,il:iu,js,kl:ku) end if end do end if 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 jp = j + 1 yi = y(j) + dyh do k = kl, ku do i = il, iu qn(:,i,j,k) = qn(:,i,nb,k) call gravitational_acceleration(t, dt, x(i), yi, z(k), ga(:)) qn(ipr,i,j,k) = qn(ipr,i,jp,k) & - (qn(idn,i,jp,k) + qn(idn,i,j,k)) * ga(2) * dyh end do end do end do else do j = neu, nn jm = j - 1 yi = y(j) - dyh do k = kl, ku do i = il, iu qn(:,i,j,k) = qn(:,i,ne,k) call gravitational_acceleration(t, dt, x(i), yi, z(k), ga(:)) qn(ipr,i,j,k) = qn(ipr,i,jm,k) & + (qn(idn,i,jm,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 jp = j + 1 yi = y(j) + dyh do k = kl, ku do i = il, iu qn(:,i,j,k) = qn(:,i,nb,k) call gravitational_acceleration(t, dt, x(i), yi, z(k), ga(:)) qn(idn,i,j,k) = qn(idn,i,jp,k) * exp(- ga(2) * dy / csnd2) end do end do end do else do j = neu, nn jm = j - 1 yi = y(j) - dyh do k = kl, ku do i = il, iu qn(:,i,j,k) = qn(:,i,ne,k) call gravitational_acceleration(t, dt, x(i), yi, z(k), ga(:)) qn(idn,i,j,k) = qn(idn,i,jm,k) * exp( ga(2) * dy / csnd2) end do end do end do end if end if case(bnd_custom) if (associated(custom_boundary_y)) then call custom_boundary_y(side(2), il, iu, kl, ku, & t, dt, x(:), y(:), z(:), qn(:,:,:,:)) else call print_message(loc, & "No subroutine associated with custom_boundary_y()!") status = 1 return end if case default if (side(2) == 1) then do j = nbl, 1, -1 qn(:,il:iu,j,kl:ku) = qn(:,il:iu,nb,kl:ku) end do else do j = neu, nn qn(:,il:iu,j,kl:ku) = qn(:,il:iu,ne,kl:ku) end do end if end select #if NDIMS == 3 case(3) 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 select case(bnd_type(dir,side(3))) case(bnd_outflow) if (side(3) == 1) then do k = nbl, 1, -1 qn( : ,il:iu,jl:ju,k) = qn( : ,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 else do k = neu, nn qn( : ,il:iu,jl:ju,k) = qn( : ,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 end if case(bnd_reflective) if (side(3) == 1) then do k = 1, ng kt = nb - k ks = nbl + k qn( : ,il:iu,jl:ju,kt) = qn( : ,il:iu,jl:ju,ks) qn(ivz,il:iu,jl:ju,kt) = - qn(ivz,il:iu,jl:ju,ks) if (magnetized) 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( : ,il:iu,jl:ju,kt) = qn( : ,il:iu,jl:ju,ks) qn(ivz,il:iu,jl:ju,kt) = - qn(ivz,il:iu,jl:ju,ks) if (magnetized) then qn(ibz,il:iu,jl:ju,kt) = - qn(ibz,il:iu,jl:ju,ks) end if end do end if 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 kp = k + 1 zi = z(k) + dzh do j = jl, ju do i = il, iu qn(:,i,j,k) = qn(:,i,j,nb) call gravitational_acceleration(t, dt, x(i), y(j), zi, ga(:)) qn(ipr,i,j,k) = qn(ipr,i,j,kp) & - (qn(idn,i,j,kp) + qn(idn,i,j,k)) * ga(3) * dzh end do end do end do else do k = neu, nn km = k - 1 zi = z(k) - dzh do j = jl, ju do i = il, iu qn(:,i,j,k) = qn(:,i,j,ne) call gravitational_acceleration(t, dt, x(i), y(j), zi, ga(:)) qn(ipr,i,j,k) = qn(ipr,i,j,km) & + (qn(idn,i,j,km) + 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 kp = k + 1 zi = z(k) + dzh do j = jl, ju do i = il, iu qn(:,i,j,k) = qn(:,i,j,nb) call gravitational_acceleration(t, dt, x(i), y(j), zi, ga(:)) qn(idn,i,j,k) = qn(idn,i,j,kp) * exp(- ga(3) * dz / csnd2) end do end do end do else do k = neu, nn km = k - 1 zi = z(k) - dzh do j = jl, ju do i = il, iu qn(:,i,j,k) = qn(:,i,j,ne) call gravitational_acceleration(t, dt, x(i), y(j), zi, ga(:)) qn(idn,i,j,k) = qn(idn,i,j,km) * exp( ga(3) * dz / csnd2) end do end do end do end if end if case(bnd_custom) if (associated(custom_boundary_z)) then call custom_boundary_z(side(3), il, iu, jl, ju, & t, dt, x(:), y(:), z(:), qn(:,:,:,:)) else call print_message(loc, & "No subroutine associated with custom_boundary_z()!") status = 1 return end if case default if (side(3) == 1) then do k = nbl, 1, -1 qn(:,il:iu,jl:ju,k) = qn(:,il:iu,jl:ju,nb) end do else do k = neu, nn qn(:,il:iu,jl:ju,k) = qn(:,il:iu,jl:ju,ne) end do end if end select #endif /* NDIMS == 3 */ end select !------------------------------------------------------------------------------- ! end subroutine block_boundary_specific ! !=============================================================================== ! ! OTHER BOUNDARY SUBROUTINES ! !=============================================================================== ! !=============================================================================== ! ! subroutine BLOCK_PRIMITIVE_VARIABLES: ! ------------------------------------ ! ! Subroutine updates primitive variables in the block interior excluding the ! ghost zones along the specified direction. ! ! Arguments: ! ! dir - the direction; ! side - the side of the boundary; ! pdata - the data block pointer; ! status - the call status; ! !=============================================================================== ! subroutine block_primitive_variables(dir, pdata, status) use blocks , only : block_data use coordinates, only : nn => bcells, nb, ne use equations , only : cons2prim implicit none integer , intent(in) :: dir integer , intent(out) :: status type(block_data), pointer, intent(inout) :: pdata integer :: il, iu integer :: jl, ju, j integer :: kl, ku, k !------------------------------------------------------------------------------- ! if (dir == 1) then il = nb iu = ne else il = 1 iu = nn end if if (dir == 2) then jl = nb ju = ne else jl = 1 ju = nn end if #if NDIMS == 3 if (dir == 3) then kl = nb ku = ne else kl = 1 ku = nn end if #else /* NDIMS == 3 */ kl = 1 ku = 1 #endif /* NDIMS == 3 */ do k = kl, ku do j = jl, ju call cons2prim(pdata%u(:,il:iu,j,k), pdata%q(:,il:iu,j,k), status) end do end do !------------------------------------------------------------------------------- ! end subroutine block_primitive_variables ! !=============================================================================== ! ! subroutine BLOCK_CONSERVATIVE_VARIABLES: ! --------------------------------------- ! ! Subroutine updates conservative variables in the ghost region specified by ! the direction and side. ! ! Arguments: ! ! dir - the direction; ! side - the side of the boundary; ! pdata - the data block pointer; ! !=============================================================================== ! subroutine block_conservative_variables(dir, side, pdata) use blocks , only : block_data use coordinates, only : nn => bcells, nb, ne use equations , only : prim2cons implicit none integer , intent(in) :: dir, side type(block_data), pointer, intent(inout) :: pdata integer :: il, iu integer :: jl, ju, j integer :: kl, ku, k logical , save :: first = .true. integer, dimension(2,2), save :: nlims !------------------------------------------------------------------------------- ! if (first) then nlims(1,1) = 1 nlims(1,2) = nb - 1 nlims(2,1) = ne + 1 nlims(2,2) = nn first = .false. end if if (dir == 1) then il = nlims(side,1) iu = nlims(side,2) else il = 1 iu = nn end if if (dir == 2) then jl = nlims(side,1) ju = nlims(side,2) else jl = 1 ju = nn end if #if NDIMS == 3 if (dir == 3) then kl = nlims(side,1) ku = nlims(side,2) else kl = 1 ku = nn end if #else /* NDIMS == 3 */ kl = 1 ku = 1 #endif /* NDIMS == 3 */ do k = kl, ku do j = jl, ju call prim2cons(pdata%q(:,il:iu,j,k), pdata%u(:,il:iu,j,k), .true.) end do end do !------------------------------------------------------------------------------- ! end subroutine block_conservative_variables #ifdef MPI ! !=============================================================================== ! ! subroutine PREPARE_EXCHANGE_ARRAY: ! --------------------------------- ! ! Subroutine prepares the arrays for block exchange lists and their counters. ! ! !=============================================================================== ! subroutine prepare_exchange_array() use mpitools, only : npmax implicit none integer :: icol, irow !------------------------------------------------------------------------------- ! do irow = 0, npmax do icol = 0, npmax nullify(barray(irow,icol)%ptr) bcount(irow,icol) = 0 end do end do !------------------------------------------------------------------------------- ! end subroutine prepare_exchange_array ! !=============================================================================== ! ! subroutine RELEASE_EXCHANGE_ARRAY: ! --------------------------------- ! ! Subroutine releases objects on the array of block exchange lists. ! ! !=============================================================================== ! subroutine release_exchange_array() use blocks , only : block_info, pointer_info use mpitools, only : npmax implicit none integer :: icol, irow type(block_info), pointer :: pinfo !------------------------------------------------------------------------------- ! do irow = 0, npmax do icol = 0, npmax pinfo => barray(irow,icol)%ptr do while(associated(pinfo)) barray(irow,icol)%ptr => pinfo%prev nullify(pinfo%prev) nullify(pinfo%next) nullify(pinfo%meta) nullify(pinfo%neigh) deallocate(pinfo) pinfo => barray(irow,icol)%ptr end do end do end do !------------------------------------------------------------------------------- ! end subroutine release_exchange_array ! !=============================================================================== ! ! subroutine APPEND_EXCHANGE_BLOCK: ! --------------------------------- ! ! Subroutine appends an info block to the array of the lists of the block ! exchange. The element is determined by the process numbers 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) use blocks, only : block_info, block_meta implicit none type(block_meta), pointer, intent(inout) :: pmeta, pneigh integer , intent(in) :: dir integer, dimension(3) , intent(in) :: pos integer :: icol, irow type(block_info), pointer :: pinfo !------------------------------------------------------------------------------- ! irow = pneigh%process icol = pmeta%process bcount(irow,icol) = bcount(irow,icol) + 1 allocate(pinfo) pinfo%meta => pmeta pinfo%neigh => pneigh pinfo%direction = dir pinfo%location(1:NDIMS) = pos(1:NDIMS) nullify(pinfo%prev) nullify(pinfo%next) if (associated(barray(irow,icol)%ptr)) & pinfo%prev => barray(irow,icol)%ptr barray(irow,icol)%ptr => pinfo !------------------------------------------------------------------------------- ! end subroutine append_exchange_block #endif /* MPI */ !=============================================================================== ! end module