Commit 9f6a9dbd authored by Jannis Teunissen's avatar Jannis Teunissen

Allow to set ghost cells manually, for more flexible b.c.

This is done by passing a 'bc_custom' routine to af_set_cc_methods instead of
the standard boundary condition method.
parent 62a80f1e
......@@ -9,7 +9,7 @@ program boundary_conditions
integer, parameter :: box_size = 8
integer, parameter :: n_iterations = 20
integer :: i_phi
integer :: i_phi, i_rho
type(af_t) :: tree
integer :: iter
......@@ -19,7 +19,9 @@ program boundary_conditions
print *, "Number of threads", af_get_max_threads()
call af_add_cc_variable(tree, "phi", ix=i_phi)
call af_add_cc_variable(tree, "rho", ix=i_rho)
call af_set_cc_methods(tree, i_phi, boundary_method)
call af_set_cc_methods(tree, i_rho, bc_custom=custom_boundary_method)
! Initialize tree
call af_init(tree, & ! Tree to initialize
......@@ -32,13 +34,13 @@ program boundary_conditions
do iter = 1, n_iterations
if (iter == 1) then
! Set initial conditions on first iteration
call af_loop_box(tree, set_phi_zero)
call af_loop_box(tree, set_initial_cond)
else
! Replace phi by average over neighbors
call af_loop_box(tree, average_phi)
call af_loop_box(tree, average_vars)
end if
call af_gc_tree(tree, [i_phi])
call af_gc_tree(tree, [i_phi, i_rho])
write(fname, "(A,I0)") "boundary_conditions_" // DIMNAME // "_", iter
call af_write_silo(tree, trim(fname), dir="output")
......@@ -47,45 +49,47 @@ program boundary_conditions
contains
! This routine sets the initial conditions for each box
subroutine set_phi_zero(box)
subroutine set_initial_cond(box)
type(box_t), intent(inout) :: box
box%cc(DTIMES(:), i_phi) = 0.0_dp
end subroutine set_phi_zero
box%cc(DTIMES(:), i_rho) = 1.0_dp
end subroutine set_initial_cond
subroutine average_phi(box)
subroutine average_vars(box)
type(box_t), intent(inout) :: box
integer :: IJK, nc
real(dp) :: tmp(DTIMES(box%n_cell))
integer :: IJK, nc, ivs(2)
real(dp) :: tmp(DTIMES(box%n_cell), 2)
nc = box%n_cell
ivs = [i_phi, i_rho]
do KJI_DO(1,nc)
#if NDIM == 1
tmp(i) = 0.5_dp * ( &
box%cc(i+1, i_phi) + &
box%cc(i-1, i_phi))
tmp(i, :) = 0.5_dp * ( &
box%cc(i+1, ivs) + &
box%cc(i-1, ivs))
#elif NDIM == 2
tmp(i, j) = 0.25_dp * ( &
box%cc(i+1, j, i_phi) + &
box%cc(i-1, j, i_phi) + &
box%cc(i, j+1, i_phi) + &
box%cc(i, j-1, i_phi))
tmp(i, j, :) = 0.25_dp * ( &
box%cc(i+1, j, ivs) + &
box%cc(i-1, j, ivs) + &
box%cc(i, j+1, ivs) + &
box%cc(i, j-1, ivs))
#elif NDIM == 3
tmp(i, j, k) = 1/6.0_dp * ( &
box%cc(i+1, j, k, i_phi) + &
box%cc(i-1, j, k, i_phi) + &
box%cc(i, j+1, k, i_phi) + &
box%cc(i, j-1, k, i_phi) + &
box%cc(i, j, k+1, i_phi) + &
box%cc(i, j, k-1, i_phi))
tmp(i, j, k, :) = 1/6.0_dp * ( &
box%cc(i+1, j, k, ivs) + &
box%cc(i-1, j, k, ivs) + &
box%cc(i, j+1, k, ivs) + &
box%cc(i, j-1, k, ivs) + &
box%cc(i, j, k+1, ivs) + &
box%cc(i, j, k-1, ivs))
#endif
end do; CLOSE_DO
! Average new and old value
box%cc(DTIMES(1:nc), i_phi) = 0.5_dp * (&
tmp(DTIMES(:)) + box%cc(DTIMES(1:nc), i_phi))
end subroutine average_phi
box%cc(DTIMES(1:nc), ivs) = 0.5_dp * (&
tmp(DTIMES(:), :) + box%cc(DTIMES(1:nc), ivs))
end subroutine average_vars
!> [boundary_method]
subroutine boundary_method(box, nb, iv, coords, bc_val, bc_type)
......@@ -112,4 +116,25 @@ contains
end subroutine boundary_method
!> [boundary_method]
end program boundary_conditions
!> With this method we can set ghost cells manually
subroutine custom_boundary_method(box, nb, iv, n_gc, cc)
type(box_t), intent(inout) :: box !< Box that needs b.c.
integer, intent(in) :: nb !< Direction
integer, intent(in) :: iv !< Index of variable
integer, intent(in) :: n_gc !< Number of ghost cells to fill
!> If n_gc > 1, fill ghost cell values in this array instead of box%cc
real(dp), intent(inout), optional :: &
cc(DTIMES(1-n_gc:box%n_cell+n_gc))
integer :: lo(NDIM), hi(NDIM)
if (n_gc /= 1) error stop "not implemented"
! Get ghost cell index range
call af_get_index_bc_inside(nb, box%n_cell, 1, lo, hi)
! Set all ghost cells to a scalar value
box%cc(DSLICE(lo,hi), iv) = 0.0_dp
end subroutine custom_boundary_method
end program
......@@ -311,16 +311,18 @@ contains
end subroutine af_set_coarse_grid
!> Set the methods for a cell-centered variable
subroutine af_set_cc_methods(tree, iv, bc, rb, prolong, restrict)
subroutine af_set_cc_methods(tree, iv, bc, rb, prolong, restrict, &
bc_custom)
use m_af_ghostcell, only: af_gc_interp
use m_af_prolong, only: af_prolong_linear
use m_af_restrict, only: af_restrict_box
type(af_t), intent(inout) :: tree !< Tree to operate on
integer, intent(in) :: iv !< Index of variable
procedure(af_subr_bc) :: bc !< Boundary condition method
procedure(af_subr_bc), optional :: bc !< Boundary condition method
procedure(af_subr_rb), optional :: rb !< Refinement boundary method
procedure(af_subr_prolong), optional :: prolong !< Prolongation method
procedure(af_subr_restrict), optional :: restrict !< Restriction method
procedure(af_subr_bc_custom), optional :: bc_custom !< Custom b.c. method
integer :: i
if (tree%has_cc_method(iv)) then
......@@ -331,7 +333,13 @@ contains
! Set methods for the variable and its copies
do i = iv, iv + tree%cc_num_copies(iv) - 1
tree%cc_methods(i)%bc => bc
if (present(bc)) then
tree%cc_methods(i)%bc => bc
else if (present(bc_custom)) then
tree%cc_methods(i)%bc_custom => bc_custom
else
error stop "af_set_cc_methods: either bc or bc_custom required"
end if
if (present(rb)) then
tree%cc_methods(i)%rb => rb
......
......@@ -65,63 +65,62 @@ contains
!> Fill ghost cells for variables ivs
subroutine af_gc_box(tree, id, ivs, corners)
type(af_t), intent(inout) :: tree !< Tree to fill ghost cells on
integer, intent(in) :: id !< Id of box for which we set ghost cells
integer, intent(in) :: ivs(:) !< Variables for which ghost cells are set
logical, intent(in), optional :: corners !< Fill corner ghost cells (default: yes)
logical :: do_corners
integer :: i, iv
type(af_t), intent(inout) :: tree !< Tree to fill ghost cells on
integer, intent(in) :: id !< Id of box for which we set ghost cells
integer, intent(in) :: ivs(:) !< Variables for which ghost cells are set
logical, intent(in), optional :: corners !< Fill corner ghost cells (default: yes)
logical :: do_corners
integer :: i, iv
integer :: nb, nb_id, bc_type
integer :: lo(NDIM), hi(NDIM), dnb(NDIM)
real(dp) :: coords(NDIM, tree%n_cell**(NDIM-1))
real(dp) :: bc_val(tree%n_cell**(NDIM-1))
do_corners = .true.
if (present(corners)) do_corners = corners
do i = 1, size(ivs)
iv = ivs(i)
if (.not. tree%has_cc_method(iv)) then
print *, "For variable ", trim(tree%cc_names(iv))
error stop "af_gc_box: no methods stored"
end if
call af_gc_box_sides(tree%boxes, id, iv, tree%cc_methods(iv)%rb, &
tree%cc_methods(iv)%bc)
if (do_corners) call af_gc_box_corner(tree%boxes, id, iv)
end do
end subroutine af_gc_box
!> Fill ghost cells for variable iv on the sides of a box, using subr_rb on
!> refinement boundaries and subr_bc on physical boundaries.
subroutine af_gc_box_sides(boxes, id, iv, subr_rb, subr_bc)
type(box_t), intent(inout) :: boxes(:) !< List of all the boxes
integer, intent(in) :: id !< Id of box for which we set ghost cells
integer, intent(in) :: iv !< Variable for which ghost cells are set
procedure(af_subr_rb) :: subr_rb !< Procedure called at refinement boundaries
procedure(af_subr_bc) :: subr_bc !< Procedure called at physical boundaries
integer :: nb, nb_id, bc_type
integer :: lo(NDIM), hi(NDIM), dnb(NDIM)
real(dp) :: coords(NDIM, boxes(id)%n_cell**(NDIM-1))
real(dp) :: bc_val(boxes(id)%n_cell**(NDIM-1))
do nb = 1, af_num_neighbors
nb_id = boxes(id)%neighbors(nb)
nb_id = tree%boxes(id)%neighbors(nb)
if (nb_id > af_no_box) then
! There is a neighbor
call af_get_index_bc_outside(nb, boxes(id)%n_cell, 1, lo, hi)
call af_get_index_bc_outside(nb, tree%boxes(id)%n_cell, 1, lo, hi)
dnb = af_neighb_offset([nb])
call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, hi, iv)
call copy_from_nb(tree%boxes(id), tree%boxes(nb_id), dnb, lo, hi, ivs)
else if (nb_id == af_no_box) then
! Refinement boundary
call subr_rb(boxes, id, nb, iv)
do i = 1, size(ivs)
iv = ivs(i)
call tree%cc_methods(iv)%rb(tree%boxes, id, nb, iv)
end do
else
! Physical boundary
call af_gc_get_boundary_coords(boxes(id), nb, coords)
call subr_bc(boxes(id), nb, iv, coords, bc_val, bc_type)
call bc_to_gc(boxes(id), nb, iv, bc_val, bc_type)
call af_gc_get_boundary_coords(tree%boxes(id), nb, coords)
do i = 1, size(ivs)
iv = ivs(i)
if (associated(tree%cc_methods(iv)%bc_custom)) then
call tree%cc_methods(iv)%bc_custom(tree%boxes(id), &
nb, iv, 1)
else
call tree%cc_methods(iv)%bc(tree%boxes(id), nb, iv, &
coords, bc_val, bc_type)
call bc_to_gc(tree%boxes(id), nb, iv, bc_val, bc_type)
end if
end do
end if
end do
end subroutine af_gc_box_sides
if (do_corners) call af_gc_box_corner(tree%boxes, id, ivs)
end subroutine af_gc_box
!> Get coordinates at the faces of a box boundary
subroutine af_gc_get_boundary_coords(box, nb, coords)
......@@ -167,10 +166,10 @@ contains
!> Fill corner ghost cells for variable iv on corners/edges of a box. If there
!> is no box to copy the data from, use linear extrapolation. This routine
!> assumes ghost cells on the sides of the box are available.
subroutine af_gc_box_corner(boxes, id, iv)
subroutine af_gc_box_corner(boxes, id, ivs)
type(box_t), intent(inout) :: boxes(:) !< List of all the boxes
integer, intent(in) :: id !< Id of box for which we set ghost cells
integer, intent(in) :: iv !< Variable for which ghost cells are set
integer, intent(in) :: ivs(:) !< Variable for which ghost cells are set
integer :: n, nb_id, dnb(NDIM), lo(NDIM)
#if NDIM == 3
integer :: hi(NDIM), dim
......@@ -192,9 +191,9 @@ contains
hi = lo
hi(dim) = boxes(id)%n_cell
dnb = af_neighb_offset(af_nb_adj_edge(:, n))
call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, hi, iv)
call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, hi, ivs)
else
call af_edge_gc_extrap(boxes(id), lo, dim, iv)
call af_edge_gc_extrap(boxes(id), lo, dim, ivs)
end if
end do
#endif
......@@ -207,9 +206,9 @@ contains
lo = af_child_dix(:, n) * (boxes(id)%n_cell + 1)
if (nb_id > af_no_box) then
call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, lo, iv)
call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, lo, ivs)
else
call af_corner_gc_extrap(boxes(id), lo, iv)
call af_corner_gc_extrap(boxes(id), lo, ivs)
end if
end do
end subroutine af_gc_box_corner
......@@ -693,21 +692,21 @@ contains
bc_val = 0.0_dp
end subroutine af_bc_set_continuous
subroutine copy_from_nb(box, box_nb, dnb, lo, hi, iv)
subroutine copy_from_nb(box, box_nb, dnb, lo, hi, ivs)
type(box_t), intent(inout) :: box !< Box on which to fill ghost cells
type(box_t), intent(in) :: box_nb !< Neighbouring box
integer, intent(in) :: dnb(NDIM) !< Neighbor spatial index offset
integer, intent(in) :: lo(NDIM) !< Ghost cell low index
integer, intent(in) :: hi(NDIM) !< Ghost cell high index
integer, intent(in) :: iv !< Ghost cell variable
integer, intent(in) :: ivs(:) !< Ghost cell variable
integer :: nlo(NDIM), nhi(NDIM)
! Get indices on neighbor
nlo = lo - dnb * box%n_cell
nhi = hi - dnb * box%n_cell
box%cc(DSLICE(lo, hi), iv) = &
box_nb%cc(DSLICE(nlo, nhi), iv)
box%cc(DSLICE(lo, hi), ivs) = &
box_nb%cc(DSLICE(nlo, nhi), ivs)
end subroutine copy_from_nb
!> Get array of cell-centered variables with multiple ghost cells, excluding corners
......@@ -762,10 +761,15 @@ contains
call af_gc_get_boundary_coords(tree%boxes(id), nb, coords)
do i = 1, size(ivs)
iv = ivs(i)
call tree%cc_methods(iv)%bc(tree%boxes(id), &
nb, iv, coords, bc_val, bc_type)
call bc_to_gc2(nc, cc(DTIMES(:), i), nb, bc_val, &
bc_type, tree%boxes(id)%dr)
if (associated(tree%cc_methods(iv)%bc_custom)) then
call tree%cc_methods(iv)%bc_custom(tree%boxes(id), &
nb, iv, 2, cc(DTIMES(:), i))
else
call tree%cc_methods(iv)%bc(tree%boxes(id), &
nb, iv, coords, bc_val, bc_type)
call bc_to_gc2(nc, cc(DTIMES(:), i), nb, bc_val, &
bc_type, tree%boxes(id)%dr)
end if
end do
end if
end do
......@@ -900,24 +904,24 @@ contains
!> This fills corner ghost cells using linear extrapolation. The ghost cells
!> on the sides already need to be filled.
subroutine af_corner_gc_extrap(box, ix, iv)
subroutine af_corner_gc_extrap(box, ix, ivs)
type(box_t), intent(inout) :: box !< Box to fill ghost cells for
integer, intent(in) :: ix(NDIM) !< Cell-index of corner
integer, intent(in) :: iv !< Variable to fill
integer, intent(in) :: ivs(:) !< Variable to fill
integer :: di(NDIM)
di = 1 - 2 * iand(ix, 1) ! 0 -> di = 1, nc+1 -> di = -1
#if NDIM == 2
box%cc(ix(1), ix(2), iv) = box%cc(ix(1)+di(1), ix(2), iv) &
+ box%cc(ix(1), ix(2)+di(2), iv) &
- box%cc(ix(1)+di(1), ix(2)+di(2), iv)
box%cc(ix(1), ix(2), ivs) = box%cc(ix(1)+di(1), ix(2), ivs) &
+ box%cc(ix(1), ix(2)+di(2), ivs) &
- box%cc(ix(1)+di(1), ix(2)+di(2), ivs)
#elif NDIM == 3
box%cc(ix(1), ix(2), ix(3), iv) = &
box%cc(ix(1), ix(2)+di(2), ix(3)+di(3), iv) + &
box%cc(ix(1)+di(1), ix(2), ix(3)+di(3), iv) + &
box%cc(ix(1)+di(1), ix(2)+di(2), ix(3), iv) - 2 * &
box%cc(ix(1)+di(1), ix(2)+di(2), ix(3)+di(3), iv)
box%cc(ix(1), ix(2), ix(3), ivs) = &
box%cc(ix(1), ix(2)+di(2), ix(3)+di(3), ivs) + &
box%cc(ix(1)+di(1), ix(2), ix(3)+di(3), ivs) + &
box%cc(ix(1)+di(1), ix(2)+di(2), ix(3), ivs) - 2 * &
box%cc(ix(1)+di(1), ix(2)+di(2), ix(3)+di(3), ivs)
#endif
end subroutine af_corner_gc_extrap
......@@ -925,11 +929,11 @@ contains
!> This fills edge ghost cells using linear extrapolation. The ghost cells on
!> the sides already need to be filled. This routine basically performs the
!> same operation as af_corner_gc_extrap does in 2D.
subroutine af_edge_gc_extrap(box, lo, dim, iv)
subroutine af_edge_gc_extrap(box, lo, dim, ivs)
type(box_t), intent(inout) :: box !< Box to operate on
integer, intent(in) :: lo(NDIM) !< Lowest index of edge ghost cells
integer, intent(in) :: dim !< Dimension parallel to edge
integer, intent(in) :: iv !< Variable to fill
integer, intent(in) :: ivs(:) !< Variable to fill
integer :: di(NDIM), ix(NDIM), ia(NDIM), ib(NDIM), ic(NDIM)
integer :: n, o_dims(NDIM-1)
......@@ -958,10 +962,10 @@ contains
ic(dim) = n
ix(dim) = n
box%cc(ix(1), ix(2), ix(3), iv) = &
box%cc(ia(1), ia(2), ia(3), iv) + &
box%cc(ib(1), ib(2), ib(3), iv) - &
box%cc(ic(1), ic(2), ic(3), iv)
box%cc(ix(1), ix(2), ix(3), ivs) = &
box%cc(ia(1), ia(2), ia(3), ivs) + &
box%cc(ib(1), ib(2), ib(3), ivs) - &
box%cc(ic(1), ic(2), ic(3), ivs)
end do
end subroutine af_edge_gc_extrap
......
......@@ -245,6 +245,8 @@ module m_af_types
procedure(af_subr_bc), pointer, nopass :: bc => null()
!> Refinement boundary routine
procedure(af_subr_rb), pointer, nopass :: rb => null()
!> Custom boundary condition routine
procedure(af_subr_bc_custom), pointer, nopass :: bc_custom => null()
end type af_cc_methods
!> The basic building block of afivo: a box with cell-centered and face
......@@ -404,6 +406,20 @@ module m_af_types
integer, intent(out) :: bc_type !< Type of b.c.
end subroutine af_subr_bc
!> To fill ghost cells near physical boundaries in a custom way. If the
!> number of ghost cells to fill is greater than one (n_gc > 1), fill ghost
!> cells in the optional argument cc.
subroutine af_subr_bc_custom(box, nb, iv, n_gc, cc)
import
type(box_t), intent(inout) :: box !< Box that needs b.c.
integer, intent(in) :: nb !< Direction
integer, intent(in) :: iv !< Index of variable
integer, intent(in) :: n_gc !< Number of ghost cells to fill
!> If n_gc > 1, fill ghost cell values in this array instead of box%cc
real(dp), intent(inout), optional :: &
cc(DTIMES(1-n_gc:box%n_cell+n_gc))
end subroutine af_subr_bc_custom
!> Subroutine for prolongation
subroutine af_subr_prolong(box_p, box_c, iv, iv_to, add)
import
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment