m_af_types.f90

#include "cpp_macros.h"
 
module m_af_types
  use iso_c_binding, only: c_ptr
 
  implicit none
  public
 
  ! dp stands for double precision (8 byte reals)
  integer, parameter :: dp = kind(0.0d0)
 
  integer, parameter :: af_max_lvl = 30
 
  integer, parameter :: af_min_lvl = 1
 
  integer, parameter :: af_max_num_vars = 1024
 
  integer, parameter :: af_rm_ref = -1
 
  integer, parameter :: af_keep_ref = 0
 
  integer, parameter :: af_do_ref = 1
 
  integer, parameter :: af_derefine = -2
 
  integer, parameter :: af_refine = 2
 
  integer, parameter :: af_no_box = 0
 
  integer, parameter :: af_phys_boundary = -1
 
  integer, parameter :: af_init_tag = -huge(1)
 
  integer, parameter :: af_xyz = 1
 
  integer, parameter :: af_cyl = 2
 
  character(len=*), parameter :: af_coord_names(2) = &
       ["Cartesian  ", "Cylindrical"]
 
  integer, parameter :: af_bc_dirichlet = -10
 
  integer, parameter :: af_bc_neumann = -11
 
  integer, parameter :: af_bc_continuous = -12
 
  integer, parameter :: af_bc_dirichlet_copy = -13
 
  integer, parameter :: af_nlen = 20
 
  type lvl_t
     integer, allocatable :: ids(:)
     integer, allocatable :: leaves(:)
     integer, allocatable :: parents(:)
  end type lvl_t
 
  type ref_lvl_t
     integer, allocatable :: add(:)
  end type ref_lvl_t
 
  type ref_info_t
     integer                      :: n_add = 0
     integer                      :: n_rm  = 0
     integer, allocatable         :: rm(:)
     type(ref_lvl_t), allocatable :: lvls(:)
  end type ref_info_t
 
#if NDIM == 1
 
  integer, parameter :: af_num_children = 2
 
  integer, parameter :: af_child_dix(1, 2) = reshape([0,1], [1,2])
  integer, parameter :: af_child_adj_nb(1, 2) = reshape([1,2], [1,2])
  logical, parameter :: af_child_low(1, 2) = reshape([.true., .false.], [1, 2])
 
  integer, parameter :: af_num_neighbors = 2
  integer, parameter :: af_neighb_lowx = 1
  integer, parameter :: af_neighb_highx = 2
 
  integer, parameter :: af_neighb_dix(1, 2) = reshape([-1,1], [1,2])
  logical, parameter :: af_neighb_low(2) = [.true., .false.]
  integer, parameter :: af_low_neighbs(1) = [1]
  integer, parameter :: af_high_neighbs(1) = [2]
  integer, parameter :: af_neighb_high_pm(2) = [-1, 1]
 
  integer, parameter :: af_neighb_rev(2) = [2, 1]
  integer, parameter :: af_neighb_dim(2) = [1, 1]
#elif NDIM == 2
 
  integer, parameter :: af_num_children = 4
 
  integer, parameter :: af_child_dix(2, 4) = reshape([0,0,1,0,0,1,1,1], [2,4])
  integer, parameter :: af_child_adj_nb(2, 4) = reshape([1,3,2,4,1,2,3,4], [2,4])
  logical, parameter :: af_child_low(2, 4) = reshape([.true., .true., &
       .false., .true., .true., .false., .false., .false.], [2, 4])
 
  integer, parameter :: af_num_neighbors = 4
  integer, parameter :: af_neighb_lowx = 1
  integer, parameter :: af_neighb_highx = 2
  integer, parameter :: af_neighb_lowy = 3
  integer, parameter :: af_neighb_highy = 4
 
  integer, parameter :: af_neighb_dix(2, 4) = reshape([-1,0,1,0,0,-1,0,1], [2,4])
  logical, parameter :: af_neighb_low(4) = [.true., .false., .true., .false.]
  integer, parameter :: af_low_neighbs(2) = [1, 3]
  integer, parameter :: af_high_neighbs(2) = [2, 4]
  integer, parameter :: af_neighb_high_pm(4) = [-1, 1, -1, 1]
 
  integer, parameter :: af_neighb_rev(4) = [2, 1, 4, 3]
  integer, parameter :: af_neighb_dim(4) = [1, 1, 2, 2]
#elif NDIM == 3
 
  integer, parameter :: af_num_children = 8
 
  integer, parameter :: af_child_dix(3, 8) = reshape( &
       [0,0,0, 1,0,0, 0,1,0, 1,1,0, &
       0,0,1, 1,0,1, 0,1,1, 1,1,1], [3,8])
  integer, parameter :: af_child_adj_nb(4, 6) = reshape( &
       [1,3,5,7, 2,4,6,8, 1,2,5,6, 3,4,7,8, 1,2,3,4, 5,6,7,8], [4,6])
  logical, parameter :: af_child_low(3, 8) = reshape([ &
       .true., .true., .true., .false., .true., .true., &
       .true., .false., .true., .false., .false., .true., &
       .true., .true., .false., .false., .true., .false., &
       .true., .false., .false., .false., .false., .false.], [3, 8])
 
  integer, parameter :: af_num_neighbors = 6
  integer, parameter :: af_neighb_lowx = 1
  integer, parameter :: af_neighb_highx = 2
  integer, parameter :: af_neighb_lowy = 3
  integer, parameter :: af_neighb_highy = 4
  integer, parameter :: af_neighb_lowz = 5
  integer, parameter :: af_neighb_highz = 6
  integer, parameter :: af_neighb_dix(3, 6) = reshape( &
       [-1,0,0, 1,0,0, 0,-1,0, 0,1,0, 0,0,-1, 0,0,1], [3,6])
  logical, parameter :: af_neighb_low(6) = &
       [.true., .false., .true., .false., .true., .false.]
  integer, parameter :: af_low_neighbs(3) = [1, 3, 5]
  integer, parameter :: af_high_neighbs(3) = [2, 4, 6]
  integer, parameter :: af_neighb_high_pm(6) = [-1, 1, -1, 1, -1, 1]
  integer, parameter :: af_neighb_rev(6) = [2, 1, 4, 3, 6, 5]
  integer, parameter :: af_neighb_dim(6) = [1, 1, 2, 2, 3, 3]
 
  integer, parameter :: af_num_edges = 12
  integer, parameter :: af_edge_dim(12) = &
       [1,1,1,1, 2,2,2,2, 3,3,3,3]
  integer, parameter :: af_edge_dir(3, 12) = reshape( &
       [0,-1,-1, 0,1,-1, 0,-1,1, 0,1,1, &
       -1,0,-1, 1,0,-1, -1,0,1, 1,0,1, &
       -1,-1,0, 1,-1,0, -1,1,0, 1,1,0], [3, 12])
  integer, parameter :: af_nb_adj_edge(2, 12) = reshape( &
       [3,5, 4,5, 3,6, 4,6, &
       1,5, 2,5, 1,6, 2,6, &
       1,3, 2,3, 1,4, 2,4], [2,12])
  integer, parameter :: af_edge_min_ix(3, 12) = reshape( &
       [0,0,0, 0,1,0, 0,0,1, 0,1,1, &
       0,0,0, 1,0,0, 0,0,1, 1,0,1, &
       0,0,0, 1,0,0, 0,1,0, 1,1,0], [3,12])
#endif
 
  type af_cc_methods
     procedure(af_subr_prolong), pointer, nopass :: prolong => null()
     procedure(af_subr_restrict), pointer, nopass :: restrict => null()
     procedure(af_subr_bc), pointer, nopass :: bc => null()
     procedure(af_subr_rb), pointer, nopass :: rb => null()
     procedure(af_subr_bc_custom), pointer, nopass :: bc_custom => null()
     procedure(af_subr_funcval), pointer, nopass :: funcval => null()
     integer :: prolong_limiter = -1
  end type af_cc_methods
 
  integer, parameter :: af_stencil_none = 0
 
  type stencil_t
     integer               :: key = af_stencil_none
     integer               :: shape = af_stencil_none
     integer               :: stype = -1
     logical               :: cylindrical_gradient = .false.
     real(dp), allocatable :: c(:)
     real(dp), allocatable :: v(:, dtimes(:))
     real(dp), allocatable :: f(dtimes(:))
     real(dp), allocatable :: bc_correction(dtimes(:))
     integer, allocatable  :: sparse_ix(:, :)
     real(dp), allocatable :: sparse_v(:, :)
  end type stencil_t
 
  type box_t
     integer               :: lvl
     logical               :: in_use=.false.  
     integer               :: tag=af_init_tag 
     integer               :: ix(ndim)
     integer               :: parent
     integer               :: children(af_num_children)
     integer               :: neighbors(af_num_neighbors)
     integer               :: neighbor_mat(dtimes(-1:1))
     integer               :: n_cell
     real(dp)              :: dr(ndim)
     real(dp)              :: r_min(ndim)
     integer               :: coord_t
     real(dp), allocatable :: cc(dtimes(:), :)
     real(dp), allocatable :: fc(dtimes(:), :, :)
 
     integer               :: n_bc = 0
     integer, allocatable  :: bc_index_to_nb(:)
     integer               :: nb_to_bc_index(af_num_neighbors)
     integer, allocatable :: bc_type(:, :)
     real(dp), allocatable :: bc_val(:, :, :)
     real(dp), allocatable :: bc_coords(:, :, :)
 
     integer         :: n_stencils = 0
     type(stencil_t), allocatable :: stencils(:)
  end type box_t
 
  type af_t
     logical  :: ready       = .false. 
     integer  :: box_limit
     integer  :: highest_lvl = 0
     integer  :: highest_id  = 0
     integer  :: n_cell
     integer  :: n_var_cell  = 0
     integer  :: n_var_face  = 0
     integer  :: coord_t
     integer  :: coarse_grid_size(ndim) = -1
     logical  :: periodic(ndim) = .false.  
     real(dp) :: r_base(ndim)
     real(dp) :: dr_base(ndim)
 
     character(len=af_nlen) :: cc_names(af_max_num_vars)
 
     character(len=af_nlen) :: fc_names(af_max_num_vars)
 
     integer :: cc_num_copies(af_max_num_vars) = 1
 
     logical :: cc_write_output(af_max_num_vars) = .true.
 
     logical :: cc_write_binary(af_max_num_vars) = .true.
 
     logical :: fc_write_binary(af_max_num_vars) = .true.
 
     type(af_cc_methods) :: cc_methods(af_max_num_vars)
 
     integer :: n_stencil_keys_stored = 0
 
     logical :: has_cc_method(af_max_num_vars) = .false.
 
     integer, allocatable :: cc_auto_vars(:)
 
     integer, allocatable :: cc_func_vars(:)
 
     type(lvl_t) :: lvls(af_min_lvl:af_max_lvl)
 
     type(box_t), allocatable :: boxes(:)
 
     integer, allocatable :: removed_ids(:)
 
     integer :: n_removed_ids = 0
 
     integer :: mg_i_eps = -1
 
     integer :: mg_i_lsf = -1
 
     integer :: mg_current_operator_mask = -1
  end type af_t
 
  type af_loc_t
     integer :: id = -1
     integer :: ix(ndim) = -1
  end type af_loc_t
 
    abstract interface
 
     subroutine af_subr_ref(box, cell_flags)
       import
       type(box_t), intent(in) :: box
       integer, intent(out)      :: cell_flags(DTIMES(box%n_cell))
     end subroutine af_subr_ref
 
     subroutine af_subr(box)
       import
       type(box_t), intent(inout) :: box
     end subroutine af_subr
 
     subroutine af_subr_arg(box, rarg)
       import
       type(box_t), intent(inout) :: box
       real(dp), intent(in)        :: rarg(:)
     end subroutine af_subr_arg
 
     subroutine af_subr_boxes(boxes, id)
       import
       type(box_t), intent(inout) :: boxes(:)
       integer, intent(in)         :: id
     end subroutine af_subr_boxes
 
     subroutine af_subr_boxes_arg(boxes, id, rarg)
       import
       type(box_t), intent(inout) :: boxes(:)
       integer, intent(in)         :: id
       real(dp), intent(in)        :: rarg(:)
     end subroutine af_subr_boxes_arg
 
     subroutine af_subr_tree(tree, id)
       import
       type(af_t), intent(inout) :: tree
       integer, intent(in)       :: id
     end subroutine af_subr_tree
 
     subroutine af_subr_tree_arg(tree, id, rarg)
       import
       type(af_t), intent(inout) :: tree
       integer, intent(in)       :: id
       real(dp), intent(in)      :: rarg(:)
     end subroutine af_subr_tree_arg
 
     subroutine af_subr_rb(boxes, id, nb, iv, op_mask)
       import
       type(box_t), intent(inout) :: boxes(:)
       integer, intent(in)        :: id
       integer, intent(in)        :: nb
       integer, intent(in)        :: iv
       integer, intent(in)        :: op_mask
     end subroutine af_subr_rb
 
     subroutine af_subr_bc(box, nb, iv, coords, bc_val, bc_type)
       import
       type(box_t), intent(in) :: box
       integer, intent(in)     :: nb
       integer, intent(in)     :: iv
       real(dp), intent(in)    :: coords(NDIM, box%n_cell**(NDIM-1))
       real(dp), intent(out)   :: bc_val(box%n_cell**(NDIM-1))
       integer, intent(out)    :: bc_type
     end subroutine af_subr_bc
 
     subroutine af_subr_bc_custom(box, nb, iv, n_gc, cc)
       import
       type(box_t), intent(inout) :: box
       integer, intent(in)     :: nb
       integer, intent(in)     :: iv
       integer, intent(in)     :: n_gc
       real(dp), intent(inout), optional :: &
            cc(DTIMES(1-n_gc:box%n_cell+n_gc))
     end subroutine af_subr_bc_custom
 
     subroutine af_subr_funcval(box, iv)
       import
       type(box_t), intent(inout) :: box
       integer, intent(in)        :: iv
     end subroutine af_subr_funcval
 
     subroutine af_subr_prolong(box_p, box_c, iv, iv_to, add, limiter)
       import
       type(box_t), intent(in)     :: box_p
       type(box_t), intent(inout)  :: box_c
       integer, intent(in)           :: iv
       integer, intent(in), optional :: iv_to
       logical, intent(in), optional :: add
       integer, intent(in), optional :: limiter
     end subroutine af_subr_prolong
 
     subroutine af_subr_restrict(box_c, box_p, ivs, use_geometry)
       import
       type(box_t), intent(in)       :: box_c
       type(box_t), intent(inout)    :: box_p
       integer, intent(in)           :: ivs(:)
       logical, intent(in), optional :: use_geometry
     end subroutine af_subr_restrict
  end interface
 
  ! *** Types related to multigrid ***
 
  ! The mg module supports different multigrid operators, and uses these tags to
  ! identify boxes / operators
  integer, parameter :: mg_normal_operator = 1
  integer, parameter :: mg_lsf_operator    = 2
  integer, parameter :: mg_eps_operator    = 3
  integer, parameter :: mg_auto_operator   = 4
 
  integer, parameter :: mg_normal_box = 0
  integer, parameter :: mg_lsf_box    = 1
  integer, parameter :: mg_veps_box   = 2
  integer, parameter :: mg_ceps_box   = 4
 
  integer, parameter :: mg_prolong_linear = 17
  integer, parameter :: mg_prolong_sparse = 18
  integer, parameter :: mg_prolong_auto   = 19
 
  integer, parameter :: mg_lsf_distance_key = 31
  integer, parameter :: mg_lsf_mask_key = 32
 
  ! Labels for the different steps of a multigrid cycle
  integer, parameter :: mg_cycle_down = 1
  integer, parameter :: mg_cycle_up   = 3
 
  type coarse_solve_t
     type(c_ptr) :: matrix
     type(c_ptr) :: rhs
     type(c_ptr) :: phi
     type(c_ptr) :: solver
     type(c_ptr) :: grid
 
     real(dp), allocatable :: bc_to_rhs(:, :, :)
 
     real(dp), allocatable :: lsf_fac(dtimes(:), :)
 
     integer  :: symmetric      = 1
     integer  :: solver_type    = -1
     integer  :: max_iterations = 50
     integer  :: n_cycle_down   = 1
     integer  :: n_cycle_up     = 1
     real(dp) :: tolerance      = 1e-6_dp
 
     integer  :: min_unknowns_for_openmp = 10*1000
  end type coarse_solve_t
 
  type :: mg_t
     integer :: i_phi        = -1
     integer :: i_rhs        = -1
     integer :: i_tmp        = -1
 
     integer :: operator_mask = -1
     integer :: i_eps = -1
     integer :: i_lsf = -1
 
     integer :: n_cycle_down = 2
     integer :: n_cycle_up   = 2
 
     logical :: initialized = .false.
     logical :: use_corners = .false.
     logical :: subtract_mean = .false.
 
     real(dp) :: helmholtz_lambda = 0.0_dp
 
     real(dp) :: lsf_boundary_value = 0.0_dp
 
     real(dp) :: lsf_gradient_safety_factor = 1.5_dp
 
     real(dp) :: lsf_length_scale = 1e100_dp
 
     real(dp) :: lsf_tol = 1e-8_dp
 
     real(dp) :: lsf_min_rel_distance = 1e-4_dp
 
     logical :: lsf_use_custom_prolongation = .false.
 
     procedure(mg_func_lsf), pointer, nopass :: lsf => null()
 
     procedure(mg_lsf_distf), pointer, nopass :: lsf_dist => null()
 
     procedure(mg_func_lsf), pointer, nopass :: lsf_boundary_function => null()
 
     procedure(af_subr_bc), pointer, nopass   :: sides_bc => null()
 
     procedure(af_subr_rb), pointer, nopass   :: sides_rb => null()
 
     procedure(mg_box_op), pointer, nopass   :: box_op => null()
 
     integer :: operator_type = mg_auto_operator
 
     integer :: operator_key = af_stencil_none
 
     integer :: prolongation_type = mg_prolong_auto
 
     integer :: prolongation_key = af_stencil_none
 
     procedure(mg_box_gsrb), pointer, nopass :: box_gsrb => null()
 
     procedure(mg_box_corr), pointer, nopass :: box_corr => null()
 
     procedure(mg_box_rstr), pointer, nopass :: box_rstr => null()
 
     procedure(mg_box_stencil), pointer, nopass :: box_stencil => null()
 
     type(coarse_solve_t) :: csolver
  end type mg_t
 
  abstract interface
 
     subroutine mg_box_op(box, i_out, mg)
       import
       type(box_t), intent(inout) :: box
       type(mg_t), intent(in)     :: mg
       integer, intent(in)        :: i_out
     end subroutine mg_box_op
 
     subroutine mg_box_gsrb(box, redblack_cntr, mg)
       import
       type(box_t), intent(inout) :: box
       type(mg_t), intent(in)     :: mg
       integer, intent(in)        :: redblack_cntr
     end subroutine mg_box_gsrb
 
     subroutine mg_box_corr(box_p, box_c, mg)
       import
       type(box_t), intent(inout) :: box_c
       type(box_t), intent(in)    :: box_p
       type(mg_t), intent(in)     :: mg
     end subroutine mg_box_corr
 
     subroutine mg_box_rstr(box_c, box_p, iv, mg)
       import
       type(box_t), intent(in)    :: box_c
       type(box_t), intent(inout) :: box_p
       integer, intent(in)        :: iv
       type(mg_t), intent(in)     :: mg
     end subroutine mg_box_rstr
 
     subroutine mg_box_stencil(box, mg, stencil, bc_to_rhs, lsf_fac)
       import
       type(box_t), intent(in) :: box
       type(mg_t), intent(in)  :: mg
       real(dp), intent(inout) :: stencil(2*NDIM+1, DTIMES(box%n_cell))
       real(dp), intent(inout) :: bc_to_rhs(box%n_cell**(NDIM-1), af_num_neighbors)
       real(dp), intent(inout) :: lsf_fac(DTIMES(box%n_cell))
     end subroutine mg_box_stencil
 
     real(dp) function mg_func_lsf(rr)
       import
       real(dp), intent(in) :: rr(ndim)
     end function mg_func_lsf
 
     real(dp) function mg_lsf_distf(a, b, mg)
       import
       real(dp), intent(in)   :: a(ndim)
       real(dp), intent(in)   :: b(ndim)
       type(mg_t), intent(in) :: mg
     end function mg_lsf_distf
  end interface
 
contains
 
  integer function af_get_max_threads()
    use omp_lib, only: omp_get_max_threads
    af_get_max_threads = omp_get_max_threads()
  end function af_get_max_threads
 
  subroutine af_print_info(tree)
    type(af_t), intent(in) :: tree
    character(len=15) :: fmt_string
 
    if (.not. allocated(tree%boxes)) then
       print *, "af_init has not been called for this tree"
    else if (.not. tree%ready) then
       print *, "af_set_base has not been called for this tree"
    else
       write(*, "(A)") ""
       write(*, "(A)") " *** af_print_info(tree) ***"
       write(*, "(A,I10)") " Current maximum level:  ", tree%highest_lvl
       write(*, "(A,I10)") " Number of boxes used:   ", af_num_boxes_used(tree)
       write(*, "(A,I10)") " Number of leaves used:  ", af_num_leaves_used(tree)
       write(*, "(A,I10)") " Memory limit(boxes):    ", tree%box_limit
       write(*, "(A,E10.2)") " Memory limit(GByte):    ", &
            tree%box_limit * 0.5_dp**30 * &
            af_box_bytes(tree%n_cell, tree%n_var_cell, tree%n_var_face)
       write(*, "(A,I10)") " Highest id in box list: ", tree%highest_id
       write(*, "(A,I10)") " Size of box list:       ", size(tree%boxes)
       write(*, "(A,I10)") " Box size (cells):       ", tree%n_cell
       write(*, "(A,I10)") " Number of cc variables: ", tree%n_var_cell
       write(*, "(A,I10)") " Number of fc variables: ", tree%n_var_face
       write(*, "(A,A15)") " Type of coordinates:    ", &
            af_coord_names(tree%coord_t)
       write(fmt_string, "(A,I0,A)") "(A,", ndim, "E10.2)"
       write(*, fmt_string) " min. coords:            ", tree%r_base
       write(*, fmt_string) " dr at level one:        ", tree%dr_base
       write(*, "(A,E10.2)") " minval(dr):             ", af_min_dr(tree)
       write(*, "(A)") ""
    end if
  end subroutine af_print_info
 
  pure function af_box_bytes(n_cell, n_var_cell, n_var_face) result(box_bytes)
    integer, intent(in) :: n_cell
    integer, intent(in) :: n_var_cell
    integer, intent(in) :: n_var_face
    integer             :: box_bytes
    type(box_t)       :: dummy_box
 
    box_bytes = 8 * n_var_cell * (n_cell + 2)**ndim + &
         8 * n_var_face * (n_cell + 1) * n_cell**(ndim-1) + &
         int(storage_size(dummy_box) / 8)
  end function af_box_bytes
 
  pure function af_num_boxes_used(tree) result(n_boxes)
    type(af_t), intent(in) :: tree
    integer                 :: n_boxes, lvl
 
    n_boxes = 0
    do lvl = 1, tree%highest_lvl
       n_boxes = n_boxes + size(tree%lvls(lvl)%ids)
    end do
  end function af_num_boxes_used
 
  pure function af_num_leaves_used(tree) result(n_boxes)
    type(af_t), intent(in) :: tree
    integer                 :: n_boxes, lvl
 
    n_boxes = 0
    do lvl = 1, tree%highest_lvl
       n_boxes = n_boxes + size(tree%lvls(lvl)%leaves)
    end do
  end function af_num_leaves_used
 
  pure function af_num_cells_used(tree) result(n_cells)
    type(af_t), intent(in) :: tree
    integer                 :: n_cells
 
    n_cells = af_num_leaves_used(tree) * tree%n_cell**ndim
  end function af_num_cells_used
 
  pure function af_total_volume(tree) result(vol)
    type(af_t), intent(in) :: tree
    real(dp)                :: vol
    integer                 :: i, id
    real(dp)                :: r0, r1, box_len(ndim)
    real(dp), parameter     :: pi = acos(-1.0_dp)
 
    box_len = tree%n_cell * tree%dr_base
 
    if (ndim == 2 .and. tree%coord_t == af_cyl) then
       vol     = 0.0_dp
       do i = 1, size(tree%lvls(1)%ids)
          id  = tree%lvls(1)%ids(i)
          r0  = tree%boxes(id)%r_min(1)
          r1  = r0 + box_len(1)
          vol = vol + pi * (r1**2 - r0**2) * box_len(ndim)
       end do
    else
       vol = size(tree%lvls(1)%ids) * product(box_len)
    end if
  end function af_total_volume
 
  elemental logical function af_has_children(box)
    type(box_t), intent(in) :: box
 
    ! Boxes are either fully refined or not, so we only need to check one of the
    ! children
    af_has_children = (box%children(1) /= af_no_box)
  end function af_has_children
 
  pure function af_is_phys_boundary(boxes, id, nbs) result(boundary)
    type(box_t), intent(in) :: boxes(:)
    integer, intent(in)       :: id
    integer, intent(in)       :: nbs(:)
    logical                   :: boundary(size(nbs))
    integer                   :: n, nb, p_id, nb_id, dim, dix
 
    do n = 1, size(nbs)
       nb = nbs(n)
       nb_id = boxes(id)%neighbors(nb)
 
       if (nb_id < af_no_box) then
          boundary(n) = .true.  ! Physical boundary
       else                     ! There could be a periodic boundary
          dim = af_neighb_dim(nb)
 
          if (nb_id == af_no_box) then
             ! Refinement boundary, compute index difference on parent (which is
             ! guaranteed to be there)
             p_id = boxes(id)%parent
             nb_id = boxes(p_id)%neighbors(nb)
             dix = boxes(nb_id)%ix(dim) - boxes(p_id)%ix(dim)
          else
             ! Normal neighbor, compute index difference
             dix = boxes(nb_id)%ix(dim) - boxes(id)%ix(dim)
          end if
 
          if (dix /= af_neighb_high_pm(nb)) then
             ! The difference in index is not the expected +-1, so a periodic
             ! boundary
             boundary(n) = .true.
          else
             boundary(n) = .false.
          end if
       end if
    end do
 
  end function af_is_phys_boundary
 
  pure logical function af_is_ref_boundary(boxes, id, nb)
    type(box_t), intent(in) :: boxes(:)
    integer, intent(in)     :: id
    integer, intent(in)     :: nb
    integer                 :: nb_id
 
    af_is_ref_boundary = .false.
    nb_id = boxes(id)%neighbors(nb)
 
    if (nb_id == af_no_box) then
       af_is_ref_boundary = .true.
    else if (nb_id > af_no_box) then
       if (.not. af_has_children(boxes(id)) .and. &
            af_has_children(boxes(nb_id))) then
          af_is_ref_boundary = .true.
       end if
    end if
  end function af_is_ref_boundary
 
  pure function af_get_child_offset(box, nb) result(ix_offset)
    type(box_t), intent(in)           :: box
    integer, intent(in), optional      :: nb
    integer                            :: ix_offset(ndim)
 
    ix_offset = iand(box%ix-1, 1) * ishft(box%n_cell, -1) ! * n_cell / 2
    if (present(nb)) ix_offset = ix_offset - af_neighb_dix(:, nb) * box%n_cell
  end function af_get_child_offset
 
  pure function af_get_ix_on_parent(box, ix) result(p_ix)
    type(box_t), intent(in) :: box
    integer, intent(in)       :: ix(ndim)
    integer                   :: p_ix(ndim)
    p_ix = af_get_child_offset(box) + ishft(ix+1, -1)
  end function af_get_ix_on_parent
 
  pure function af_get_ix_on_neighb(box, ix, nb) result(nb_ix)
    type(box_t), intent(in) :: box
    integer, intent(in)       :: ix(ndim)
    integer, intent(in)       :: nb
    integer                   :: nb_ix(ndim), nb_dim
 
    nb_dim        = af_neighb_dim(nb)
    nb_ix         = ix
    nb_ix(nb_dim) = nb_ix(nb_dim) - af_neighb_high_pm(nb) * box%n_cell
  end function af_get_ix_on_neighb
 
  pure function af_neighb_offset(nbs) result(dix)
    integer, intent(in) :: nbs(:)
    integer             :: n, dim, dix(ndim)
 
    dix = 0
    do n = 1, size(nbs)
       dim = af_neighb_dim(nbs(n))
       dix(dim) = dix(dim) + af_neighb_high_pm(nbs(n))
    end do
  end function af_neighb_offset
 
  subroutine af_get_index_bc_inside(nb, nc, n_gc, lo, hi)
    integer, intent(in)  :: nb
    integer, intent(in)  :: nc
    integer, intent(in)  :: n_gc
    integer, intent(out) :: lo(NDIM)
    integer, intent(out) :: hi(NDIM)
    integer              :: nb_dim
 
    ! Determine index range next to boundary
    nb_dim     = af_neighb_dim(nb)
    lo(:)      = 1
    hi(:)      = nc
    if (af_neighb_low(nb)) then
       lo(nb_dim) = 1
       hi(nb_dim) = n_gc
    else
       lo(nb_dim) = nc - n_gc + 1
       hi(nb_dim) = nc
    end if
  end subroutine af_get_index_bc_inside
 
  subroutine af_get_index_bc_outside(nb, nc, n_gc, lo, hi)
    integer, intent(in)  :: nb
    integer, intent(in)  :: nc
    integer, intent(in)  :: n_gc
    integer, intent(out) :: lo(NDIM)
    integer, intent(out) :: hi(NDIM)
    integer              :: nb_dim
 
    ! Determine index range next to boundary
    nb_dim     = af_neighb_dim(nb)
    lo(:)      = 1
    hi(:)      = nc
    if (af_neighb_low(nb)) then
       lo(nb_dim) = 1 - n_gc
       hi(nb_dim) = 0
    else
       lo(nb_dim) = nc + 1
       hi(nb_dim) = nc + n_gc
    end if
  end subroutine af_get_index_bc_outside
 
  subroutine af_get_index_bface_inside(nb, nc, n_gc, lo, hi)
    integer, intent(in)  :: nb
    integer, intent(in)  :: nc
    integer, intent(in)  :: n_gc
    integer, intent(out) :: lo(NDIM)
    integer, intent(out) :: hi(NDIM)
    integer              :: nb_dim
 
    ! Determine index range next to boundary
    nb_dim     = af_neighb_dim(nb)
    lo(:)      = 1
    hi(:)      = nc
    if (af_neighb_low(nb)) then
       lo(nb_dim) = 1
       hi(nb_dim) = n_gc
    else
       lo(nb_dim) = nc - n_gc + 2
       hi(nb_dim) = nc + 1
    end if
  end subroutine af_get_index_bface_inside
 
  pure integer function af_ix_to_ichild(ix)
    integer, intent(in) :: ix(ndim)
    ! The index can range from 1 (all ix odd) and 2**NDIM (all ix even)
#if NDIM == 1
    af_ix_to_ichild = 2 - iand(ix(1), 1)
#elif NDIM == 2
    af_ix_to_ichild = 4 - 2 * iand(ix(2), 1) - iand(ix(1), 1)
#elif NDIM == 3
    af_ix_to_ichild = &
         8 - 4 * iand(ix(3), 1) - 2 * iand(ix(2), 1) - iand(ix(1), 1)
#endif
  end function af_ix_to_ichild
 
  pure function af_cc_ix(box, r) result(cc_ix)
    type(box_t), intent(in) :: box
    real(dp), intent(in)     :: r(ndim)
    integer                  :: cc_ix(ndim)
    cc_ix = ceiling((r - box%r_min) / box%dr)
  end function af_cc_ix
 
  pure function af_r_cc(box, cc_ix) result(r)
    type(box_t), intent(in) :: box
    integer, intent(in)      :: cc_ix(ndim)
    real(dp)                 :: r(ndim)
    r = box%r_min + (cc_ix-0.5_dp) * box%dr
  end function af_r_cc
 
  pure function af_r_fc(box, dim, fc_ix) result(r)
    type(box_t), intent(in) :: box
    integer, intent(in)     :: dim
    integer, intent(in)     :: fc_ix(ndim)
    real(dp)                :: r(ndim)
    r      = box%r_min + (fc_ix-0.5_dp) * box%dr
    r(dim) = r(dim) - 0.5_dp * box%dr(dim)
  end function af_r_fc
 
  pure function af_rr_cc(box, cc_ix) result(r)
    type(box_t), intent(in) :: box
    real(dp), intent(in)     :: cc_ix(ndim)
    real(dp)                 :: r(ndim)
    r = box%r_min + (cc_ix-0.5_dp) * box%dr
  end function af_rr_cc
 
  pure function af_r_center(box) result(r_center)
    type(box_t), intent(in) :: box
    real(dp)                 :: r_center(ndim)
    r_center = box%r_min + 0.5_dp * box%n_cell * box%dr
  end function af_r_center
 
  pure function af_min_dr(tree) result(dr)
    type(af_t), intent(in) :: tree
    real(dp)               :: dr
    dr = minval(af_lvl_dr(tree, tree%highest_lvl))
  end function af_min_dr
 
  pure function af_lvl_dr(tree, lvl) result(dr)
    type(af_t), intent(in) :: tree
    integer, intent(in)    :: lvl
    real(dp)               :: dr(ndim)
    dr = tree%dr_base * 0.5_dp**(lvl-1)
  end function af_lvl_dr
 
  subroutine af_set_box_gc(box, nb, iv, gc_scalar, gc_array)
    type(box_t), intent(inout) :: box
    integer, intent(in)         :: nb
    integer, intent(in)         :: iv
    real(dp), optional, intent(in) :: gc_scalar
#if NDIM == 1
    real(dp), optional, intent(in) :: gc_array(1)
#elif NDIM == 2
    real(dp), optional, intent(in) :: gc_array(box%n_cell)
#elif NDIM == 3
 
    real(dp), optional, intent(in) :: gc_array(box%n_cell, box%n_cell)
#endif
    integer                     :: nc
 
    nc = box%n_cell
 
    if (present(gc_array)) then
       select case (nb)
#if NDIM == 1
       case (af_neighb_lowx)
          box%cc(0, iv)    = gc_array(1)
       case (af_neighb_highx)
          box%cc(nc+1, iv) = gc_array(1)
#elif NDIM == 2
       case (af_neighb_lowx)
          box%cc(0, 1:nc, iv)    = gc_array
       case (af_neighb_highx)
          box%cc(nc+1, 1:nc, iv) = gc_array
       case (af_neighb_lowy)
          box%cc(1:nc, 0, iv)    = gc_array
       case (af_neighb_highy)
          box%cc(1:nc, nc+1, iv) = gc_array
#elif NDIM == 3
       case (af_neighb_lowx)
          box%cc(0, 1:nc, 1:nc, iv)    = gc_array
       case (af_neighb_highx)
          box%cc(nc+1, 1:nc, 1:nc, iv) = gc_array
       case (af_neighb_lowy)
          box%cc(1:nc, 0, 1:nc, iv)    = gc_array
       case (af_neighb_highy)
          box%cc(1:nc, nc+1, 1:nc, iv) = gc_array
       case (af_neighb_lowz)
          box%cc(1:nc, 1:nc, 0, iv)    = gc_array
       case (af_neighb_highz)
          box%cc(1:nc, 1:nc, nc+1, iv) = gc_array
#endif
       end select
    else if (present(gc_scalar)) then
       select case (nb)
#if NDIM == 1
       case (af_neighb_lowx)
          box%cc(0, iv)    = gc_scalar
       case (af_neighb_highx)
          box%cc(nc+1, iv) = gc_scalar
#elif NDIM == 2
       case (af_neighb_lowx)
          box%cc(0, 1:nc, iv)    = gc_scalar
       case (af_neighb_highx)
          box%cc(nc+1, 1:nc, iv) = gc_scalar
       case (af_neighb_lowy)
          box%cc(1:nc, 0, iv)    = gc_scalar
       case (af_neighb_highy)
          box%cc(1:nc, nc+1, iv) = gc_scalar
#elif NDIM == 3
       case (af_neighb_lowx)
          box%cc(0, 1:nc, 1:nc, iv)    = gc_scalar
       case (af_neighb_highx)
          box%cc(nc+1, 1:nc, 1:nc, iv) = gc_scalar
       case (af_neighb_lowy)
          box%cc(1:nc, 0, 1:nc, iv)    = gc_scalar
       case (af_neighb_highy)
          box%cc(1:nc, nc+1, 1:nc, iv) = gc_scalar
       case (af_neighb_lowz)
          box%cc(1:nc, 1:nc, 0, iv)    = gc_scalar
       case (af_neighb_highz)
          box%cc(1:nc, 1:nc, nc+1, iv) = gc_scalar
#endif
       end select
    else
       stop "af_set_box_gc: requires gc_scalar or gc_array"
    end if
  end subroutine af_set_box_gc
 
#if NDIM == 2
 
  pure function af_cyl_radius_cc(box, i) result(r)
    type(box_t), intent(in) :: box
    integer, intent(in)      :: i
    real(dp)                 :: r
    r = box%r_min(1) + (i-0.5_dp) * box%dr(1)
  end function af_cyl_radius_cc
 
  pure function af_cyl_volume_cc(box, i) result(dvol)
    type(box_t), intent(in) :: box
    integer, intent(in)     :: i
    real(dp), parameter     :: two_pi = 2 * acos(-1.0_dp)
    real(dp)                :: dvol
    dvol = two_pi * (box%r_min(1) + (i-0.5_dp) * box%dr(1)) * product(box%dr)
  end function af_cyl_volume_cc
 
  subroutine af_cyl_child_weights(box, i, inner, outer)
    type(box_t), intent(in) :: box
    integer, intent(in)      :: i
    real(dp), intent(out)    :: inner, outer
    real(dp)                 :: tmp
 
    tmp = 0.25_dp * box%dr(1) / af_cyl_radius_cc(box, i)
    inner = 1 - tmp
    outer = 1 + tmp
  end subroutine af_cyl_child_weights
 
  subroutine af_cyl_flux_factors(box, flux_factors)
    type(box_t), intent(in) :: box
    real(dp), intent(out)   :: flux_factors(2, box%n_cell)
    integer                 :: i
    real(dp)                :: r, inv_r
 
    do i = 1, box%n_cell
       r                  = af_cyl_radius_cc(box, i)
       inv_r              = 1/r
       flux_factors(1, i) = (r - 0.5_dp * box%dr(1)) * inv_r
       flux_factors(2, i) = (r + 0.5_dp * box%dr(1)) * inv_r
    end do
  end subroutine af_cyl_flux_factors
#endif
 
  subroutine af_get_face_coords(box, nb, coords)
    type(box_t), intent(in) :: box
    integer, intent(in)     :: nb
    real(dp), intent(out)   :: coords(NDIM, box%n_cell**(NDIM-1))
    integer                 :: i, nb_dim, bc_dim(NDIM-1)
    integer, parameter      :: all_dims(NDIM) = [(i, i = 1, ndim)]
    real(dp)                :: rmin(NDIM)
#if NDIM == 3
    integer                 :: j, ix
#endif
 
    nb_dim       = af_neighb_dim(nb)
    bc_dim       = pack(all_dims, all_dims /= nb_dim)
    rmin(bc_dim) = box%r_min(bc_dim) + 0.5_dp * box%dr(bc_dim)
 
    if (af_neighb_low(nb)) then
       rmin(nb_dim) = box%r_min(nb_dim)
    else
       rmin(nb_dim) = box%r_min(nb_dim) + box%n_cell * box%dr(nb_dim)
    end if
 
#if NDIM == 1
    coords(nb_dim, 1) = rmin(nb_dim)
#elif NDIM == 2
    do i = 1, box%n_cell
       coords(bc_dim, i) = rmin(bc_dim) + (i-1) * box%dr(bc_dim)
       coords(nb_dim, i) = rmin(nb_dim)
    end do
#elif NDIM == 3
    ix = 0
    do j = 1, box%n_cell
       do i = 1, box%n_cell
          ix = ix + 1
          coords(bc_dim, ix) = rmin(bc_dim) + [i-1, j-1] * box%dr(bc_dim)
          coords(nb_dim, ix) = rmin(nb_dim)
       end do
    end do
#endif
  end subroutine af_get_face_coords
 
end module m_af_types