subroutine mg_build_rectangle(mg, domain_size, box_size, dx, r_min, periodic,&
n_finer)
type(mg_t), intent(inout) :: mg
integer, intent(in) :: domain_size(3)
integer, intent(in) :: box_size
real(dp), intent(in) :: dx(3)
real(dp), intent(in) :: r_min(3)
logical, intent(in) :: periodic(3)
integer, intent(in) :: n_finer
integer :: i, j, k, lvl, n, id, nx(3), IJK_vec(3), idim
integer :: boxes_per_dim(3, mg_lvl_lo:1)
integer :: periodic_offset(3)
if (modulo(box_size, 2) /= 0) error stop "box_size should be even"
if (any(modulo(domain_size, box_size) /= 0)) error stop "box_size does not divide domain_size"
if (all(periodic)) then
mg%subtract_mean = .true.
end if
nx = domain_size
mg%box_size = box_size
mg%box_size_lvl(1) = box_size
mg%domain_size_lvl(:, 1) = domain_size
mg%first_normal_lvl = 1
mg%dr(:, 1) = dx
mg%r_min(:) = r_min
mg%periodic = periodic
boxes_per_dim(:, :) = 0
boxes_per_dim(:, 1) = domain_size / box_size
do lvl = 1, mg_lvl_lo+1, -1
! For a Gauss-Seidel (non red-black) smoother, we should avoid boxes
! containing a single cell
if (any(modulo(nx, 2) == 1 .or. nx == mg%coarsest_grid .or. &
(mg%box_size_lvl(lvl) == mg%coarsest_grid .and. mg%smoother_type == &
mg_smoother_gs))) exit
if (all(modulo(nx/mg%box_size_lvl(lvl), 2) == 0)) then
mg%box_size_lvl(lvl-1) = mg%box_size_lvl(lvl)
boxes_per_dim(:, lvl-1) = boxes_per_dim(:, lvl)/2
mg%first_normal_lvl = lvl-1
else
mg%box_size_lvl(lvl-1) = mg%box_size_lvl(lvl)/2
boxes_per_dim(:, lvl-1) = boxes_per_dim(:, lvl)
end if
mg%dr(:, lvl-1) = mg%dr(:, lvl) * 2
nx = nx / 2
mg%domain_size_lvl(:, lvl-1) = nx
end do
mg%lowest_lvl = lvl
mg%highest_lvl = 1
do lvl = 2, mg_lvl_hi
mg%dr(:, lvl) = mg%dr(:, lvl-1) * 0.5_dp
mg%box_size_lvl(lvl) = box_size
mg%domain_size_lvl(:, lvl) = 2 * mg%domain_size_lvl(:, lvl-1)
end do
n = sum(product(boxes_per_dim, dim=1)) + n_finer
allocate(mg%boxes(n))
! Create lowest level
nx = boxes_per_dim(:, mg%lowest_lvl)
periodic_offset = [nx(1)-1, (nx(2)-1)*nx(1), (nx(3)-1) * nx(2) * nx(1)]
do k=1,nx(3); do j=1,nx(2); do i=1,nx(1)
mg%n_boxes = mg%n_boxes + 1
n = mg%n_boxes
mg%boxes(n)%rank = 0
mg%boxes(n)%id = n
mg%boxes(n)%lvl = mg%lowest_lvl
mg%boxes(n)%ix(:) = [i, j, k]
mg%boxes(n)%r_min(:) = r_min + (mg%boxes(n)%ix(:) - 1) * &
mg%box_size_lvl(mg%lowest_lvl) * mg%dr(:, mg%lowest_lvl)
mg%boxes(n)%dr(:) = mg%dr(:, mg%lowest_lvl)
mg%boxes(n)%parent = mg_no_box
mg%boxes(n)%children(:) = mg_no_box
! Set default neighbors
mg%boxes(n)%neighbors(:) = [n-1, n+1, n-nx(1), n+nx(1), n-nx(1)*nx(2),&
n+nx(1)*nx(2)]
! Handle boundaries
IJK_vec = [i, j, k]
do idim = 1, 3
if (IJK_vec(idim) == 1) then
if (periodic(idim)) then
mg%boxes(n)%neighbors(2*idim-1) = n + periodic_offset(idim)
else
mg%boxes(n)%neighbors(2*idim-1) = mg_physical_boundary
end if
end if
if (IJK_vec(idim) == nx(idim)) then
if (periodic(idim)) then
mg%boxes(n)%neighbors(2*idim) = n - periodic_offset(idim)
else
mg%boxes(n)%neighbors(2*idim) = mg_physical_boundary
end if
end if
end do
end do; end do; end do
mg%lvls(mg%lowest_lvl)%ids = [(n, n=1, mg%n_boxes)]
! Add higher levels
do lvl = mg%lowest_lvl, 0
if (mg%box_size_lvl(lvl+1) == mg%box_size_lvl(lvl)) then
do i = 1, size(mg%lvls(lvl)%ids)
id = mg%lvls(lvl)%ids(i)
call mg_add_children(mg, id)
end do
call mg_set_leaves_parents(mg%boxes, mg%lvls(lvl))
call mg_set_next_level_ids(mg, lvl)
call mg_set_neighbors_lvl(mg, lvl+1)
else
do i = 1, size(mg%lvls(lvl)%ids)
id = mg%lvls(lvl)%ids(i)
call add_single_child(mg, id, size(mg%lvls(lvl)%ids))
end do
call mg_set_leaves_parents(mg%boxes, mg%lvls(lvl))
call mg_set_next_level_ids(mg, lvl)
end if
end do
call mg_set_leaves_parents(mg%boxes, mg%lvls(1))
! No refinement boundaries
do lvl = mg%lowest_lvl, 1
if (allocated(mg%lvls(lvl)%ref_bnds)) deallocate(mg%lvls(lvl)%ref_bnds)
allocate(mg%lvls(lvl)%ref_bnds(0))
end do
mg%tree_created = .true.
end subroutine mg_build_rectangle