Afivo  0.3
m_af_ghostcell.f90
1 !> This module contains routines related to the filling of ghost cells. Note that
2 !> corner ghost cells are not used in Afivo.
3 module m_af_ghostcell
4 #include "cpp_macros.h"
5  use m_af_types
6 
7  implicit none
8  private
9 
10  public :: af_gc_tree
11  public :: af_gc_lvl
12  public :: af_gc_box
13  public :: af_bc_dirichlet_zero
14  public :: af_bc_neumann_zero
15  public :: af_bc_set_continuous
16  public :: af_gc_interp
17  public :: af_gc_prolong_copy
18  public :: af_gc_interp_lim
19  public :: af_gc2_box
20 
21 contains
22 
23  !> Fill ghost cells for variables ivs on the sides of all boxes, using
24  !> subr_rb on refinement boundaries and subr_bc on physical boundaries
25  subroutine af_gc_tree(tree, ivs, corners, leaves_only)
26  type(af_t), intent(inout) :: tree !< Tree to fill ghost cells on
27  integer, intent(in) :: ivs(:) !< Variables for which ghost cells are set
28  logical, intent(in), optional :: corners !< Fill corner ghost cells (default: yes)
29  logical, intent(in), optional :: leaves_only !< Fill only leaves' ghost cells (default: false)
30  integer :: lvl
31  logical :: all_ids
32 
33  if (.not. tree%ready) error stop "af_gc_tree: tree not ready"
34  all_ids = .true.
35  if (present(leaves_only)) all_ids = .not. leaves_only
36 
37  if (all_ids) then
38  do lvl = 1, tree%highest_lvl
39  call af_gc_lvl(tree, lvl, ivs, corners)
40  end do
41  else
42  do lvl = 1, tree%highest_lvl
43  call af_gc_lvl(tree, lvl, ivs, corners)
44  end do
45  end if
46  end subroutine af_gc_tree
47 
48  !> Fill ghost cells for variables ivs on the sides of all boxes
49  subroutine af_gc_lvl(tree, lvl, ivs, corners)
50  type(af_t), intent(inout) :: tree !< Tree to fill ghost cells on
51  integer, intent(in) :: lvl !< Fill on this refinement level
52  integer, intent(in) :: ivs(:) !< Variables for which ghost cells are set
53  logical, intent(in), optional :: corners !< Fill corner ghost cells (default: yes)
54  integer :: i
55 
56  !$omp parallel do
57  do i = 1, size(tree%lvls(lvl)%ids)
58  call af_gc_box(tree, tree%lvls(lvl)%ids(i), ivs, corners)
59  end do
60  !$omp end parallel do
61  end subroutine af_gc_lvl
62 
63  !> Fill ghost cells for variables ivs
64  subroutine af_gc_box(tree, id, ivs, corners)
65  type(af_t), intent(inout) :: tree !< Tree to fill ghost cells on
66  integer, intent(in) :: id !< Id of box for which we set ghost cells
67  integer, intent(in) :: ivs(:) !< Variables for which ghost cells are set
68  logical, intent(in), optional :: corners !< Fill corner ghost cells (default: yes)
69  logical :: do_corners
70  integer :: i, iv
71  integer :: nb, nb_id, bc_type, ix
72  integer :: lo(ndim), hi(ndim), dnb(ndim)
73  real(dp) :: bc_val(tree%n_cell**(ndim-1))
74 
75  do_corners = .true.
76  if (present(corners)) do_corners = corners
77 
78  do i = 1, size(ivs)
79  iv = ivs(i)
80  if (.not. tree%has_cc_method(iv)) then
81  print *, "For variable ", trim(tree%cc_names(iv))
82  error stop "af_gc_box: no methods stored"
83  end if
84  end do
85 
86  do nb = 1, af_num_neighbors
87  nb_id = tree%boxes(id)%neighbors(nb)
88 
89  if (nb_id > af_no_box) then
90  ! There is a neighbor
91  call af_get_index_bc_outside(nb, tree%boxes(id)%n_cell, 1, lo, hi)
92  dnb = af_neighb_offset([nb])
93  call copy_from_nb(tree%boxes(id), tree%boxes(nb_id), dnb, lo, hi, ivs)
94  else if (nb_id == af_no_box) then
95  ! Refinement boundary
96  do i = 1, size(ivs)
97  iv = ivs(i)
98  call tree%cc_methods(iv)%rb(tree%boxes, id, nb, iv, &
99  tree%mg_current_operator_mask)
100  end do
101  else
102  do i = 1, size(ivs)
103  iv = ivs(i)
104  if (associated(tree%cc_methods(iv)%bc_custom)) then
105  call tree%cc_methods(iv)%bc_custom(tree%boxes(id), &
106  nb, iv, 1)
107  else
108  ix = tree%boxes(id)%nb_to_bc_index(nb)
109  call tree%cc_methods(iv)%bc(tree%boxes(id), nb, iv, &
110  tree%boxes(id)%bc_coords(:, :, ix), bc_val, bc_type)
111  call bc_to_gc(tree%boxes(id), nb, iv, bc_val, bc_type)
112  tree%boxes(id)%bc_val(:, iv, ix) = bc_val
113  tree%boxes(id)%bc_type(iv, ix) = bc_type
114  end if
115  end do
116  end if
117  end do
118 
119  if (do_corners) call af_gc_box_corner(tree%boxes, id, ivs)
120  end subroutine af_gc_box
121 
122  !> Fill corner ghost cells for variable iv on corners/edges of a box. If there
123  !> is no box to copy the data from, use linear extrapolation. This routine
124  !> assumes ghost cells on the sides of the box are available.
125  subroutine af_gc_box_corner(boxes, id, ivs)
126  type(box_t), intent(inout) :: boxes(:) !< List of all the boxes
127  integer, intent(in) :: id !< Id of box for which we set ghost cells
128  integer, intent(in) :: ivs(:) !< Variable for which ghost cells are set
129  integer :: n, nb_id, dnb(NDIM), lo(NDIM)
130 #if NDIM == 3
131  integer :: hi(NDIM), dim
132 #endif
133 
134 #if NDIM == 3
135  ! Have to do edges before corners (since extrapolation can depend on edge values)
136  do n = 1, af_num_edges
137  dim = af_edge_dim(n)
138 
139  ! Check whether there is a neighbor, and find its index
140  nb_id = boxes(id)%neighbor_mat(af_edge_dir(1, n), &
141  af_edge_dir(2, n), af_edge_dir(3, n))
142 
143  lo = af_edge_min_ix(:, n) * (boxes(id)%n_cell + 1)
144  lo(dim) = 1
145 
146  if (nb_id > af_no_box) then
147  hi = lo
148  hi(dim) = boxes(id)%n_cell
149  dnb = af_neighb_offset(af_nb_adj_edge(:, n))
150  call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, hi, ivs)
151  else
152  call af_edge_gc_extrap(boxes(id), lo, dim, ivs)
153  end if
154  end do
155 #endif
156 
157  do n = 1, af_num_children
158  ! Check whether there is a neighbor, and find its index
159  dnb = 2 * af_child_dix(:, n) - 1
160 
161  nb_id = boxes(id)%neighbor_mat(dindex(dnb))
162  lo = af_child_dix(:, n) * (boxes(id)%n_cell + 1)
163 
164  if (nb_id > af_no_box) then
165  call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, lo, ivs)
166  else
167  call af_corner_gc_extrap(boxes(id), lo, ivs)
168  end if
169  end do
170  end subroutine af_gc_box_corner
171 
172  !> Convert a boundary condition to ghost cell data
173  subroutine bc_to_gc(box, nb, iv, bc_val, bc_type)
174  type(box_t), intent(inout) :: box
175  integer, intent(in) :: iv !< Variable to fill
176  integer, intent(in) :: nb !< Neighbor direction
177  real(dp), intent(in) :: bc_val(box%n_cell**(NDIM-1))
178  integer, intent(in) :: bc_type !< Type of b.c.
179  real(dp) :: c0, c1, c2
180  integer :: nc
181 
182  nc = box%n_cell
183 
184  ! If we call the interior point x1, x2 and the ghost point x0, then a
185  ! Dirichlet boundary value b can be imposed as:
186  ! x0 = -x1 + 2*b
187  ! A Neumann b.c. can be imposed as:
188  ! x0 = x1 +/- dx * b
189  ! A continuous boundary (same slope) as:
190  ! x0 = 2 * x1 - x2
191  ! Below, we set coefficients to handle these cases
192  select case (bc_type)
193  case (af_bc_dirichlet)
194  c0 = 2
195  c1 = -1
196  c2 = 0
197  case (af_bc_neumann)
198  c0 = box%dr(af_neighb_dim(nb)) * af_neighb_high_pm(nb) ! This gives a + or - sign
199  c1 = 1
200  c2 = 0
201  case (af_bc_continuous)
202  c0 = 0
203  c1 = 2
204  c2 = -1
205  case (af_bc_dirichlet_copy)
206  c0 = 1
207  c1 = 0
208  c2 = 0
209  case default
210  stop "fill_bc: unknown boundary condition"
211  end select
212 
213  select case (nb)
214 #if NDIM == 1
215  case (af_neighb_lowx)
216  box%cc(0, iv) = &
217  c0 * bc_val(1) + &
218  c1 * box%cc(1, iv) + &
219  c2 * box%cc(2, iv)
220  case (af_neighb_highx)
221  box%cc(nc+1, iv) = &
222  c0 * bc_val(1) + &
223  c1 * box%cc(nc, iv) + &
224  c2 * box%cc(nc-1, iv)
225 #elif NDIM == 2
226  case (af_neighb_lowx)
227  box%cc(0, 1:nc, iv) = &
228  c0 * bc_val + &
229  c1 * box%cc(1, 1:nc, iv) + &
230  c2 * box%cc(2, 1:nc, iv)
231  case (af_neighb_highx)
232  box%cc(nc+1, 1:nc, iv) = &
233  c0 * bc_val + &
234  c1 * box%cc(nc, 1:nc, iv) + &
235  c2 * box%cc(nc-1, 1:nc, iv)
236  case (af_neighb_lowy)
237  box%cc(1:nc, 0, iv) = &
238  c0 * bc_val + &
239  c1 * box%cc(1:nc, 1, iv) + &
240  c2 * box%cc(1:nc, 2, iv)
241  case (af_neighb_highy)
242  box%cc(1:nc, nc+1, iv) = &
243  c0 * bc_val + &
244  c1 * box%cc(1:nc, nc, iv) + &
245  c2 * box%cc(1:nc, nc-1, iv)
246 #elif NDIM == 3
247  case (af_neighb_lowx)
248  box%cc(0, 1:nc, 1:nc, iv) = &
249  c0 * reshape(bc_val, [nc, nc]) + &
250  c1 * box%cc(1, 1:nc, 1:nc, iv) + &
251  c2 * box%cc(2, 1:nc, 1:nc, iv)
252  case (af_neighb_highx)
253  box%cc(nc+1, 1:nc, 1:nc, iv) = &
254  c0 * reshape(bc_val, [nc, nc]) + &
255  c1 * box%cc(nc, 1:nc, 1:nc, iv) + &
256  c2 * box%cc(nc-1, 1:nc, 1:nc, iv)
257  case (af_neighb_lowy)
258  box%cc(1:nc, 0, 1:nc, iv) = &
259  c0 * reshape(bc_val, [nc, nc]) + &
260  c1 * box%cc(1:nc, 1, 1:nc, iv) + &
261  c2 * box%cc(1:nc, 2, 1:nc, iv)
262  case (af_neighb_highy)
263  box%cc(1:nc, nc+1, 1:nc, iv) = &
264  c0 * reshape(bc_val, [nc, nc]) + &
265  c1 * box%cc(1:nc, nc, 1:nc, iv) + &
266  c2 * box%cc(1:nc, nc-1, 1:nc, iv)
267  case (af_neighb_lowz)
268  box%cc(1:nc, 1:nc, 0, iv) = &
269  c0 * reshape(bc_val, [nc, nc]) + &
270  c1 * box%cc(1:nc, 1:nc, 1, iv) + &
271  c2 * box%cc(1:nc, 1:nc, 2, iv)
272  case (af_neighb_highz)
273  box%cc(1:nc, 1:nc, nc+1, iv) = &
274  c0 * reshape(bc_val, [nc, nc]) + &
275  c1 * box%cc(1:nc, 1:nc, nc, iv) + &
276  c2 * box%cc(1:nc, 1:nc, nc-1, iv)
277 #endif
278  end select
279  end subroutine bc_to_gc
280 
281  !> Convert a boundary condition to two layers of ghost cell data
282  subroutine bc_to_gc2(nc, cc, nb, bc_val, bc_type, dr)
283  integer, intent(in) :: nc !< Number of cells
284  real(dp), intent(inout) :: cc(DTIMES(-1:nc+2)) !< Cell-centered data
285  integer, intent(in) :: nb !< Neighbor direction
286  real(dp), intent(in) :: bc_val(nc**(NDIM-1)) !< Boundary condition
287  integer, intent(in) :: bc_type !< Type of b.c.
288  real(dp), intent(in) :: dr(NDIM) !< Grid spacing
289  real(dp) :: c0, c1, c2
290 
291  ! If we call the interior point x1, x2 and the ghost point x0, then a
292  ! Dirichlet boundary value b can be imposed as:
293  ! x0 = -x1 + 2*b = c0 * b + c1 * x1
294  ! x-1 = -x2 + 2*b = c2 * b + c1 * x2
295  ! A Neumann b.c. can be imposed as:
296  ! x0 = x1 +/- dx * b = c0 * b + c1 * x1
297  ! x-1 = x2 +/- 3 * dx * b = c2 * b + c1 * x2
298  !
299  ! The second ghost cell here is a copy of the first one, this might not
300  ! always be ideal, but it ensures the af_bc_dirichlet_copy variant does not
301  ! introduce negative values
302  !
303  ! Below, we set coefficients to handle these cases
304  select case (bc_type)
305  case (af_bc_dirichlet)
306  c0 = 2
307  c1 = -1
308  c2 = c0
309  case (af_bc_neumann)
310  c0 = dr(af_neighb_dim(nb)) * af_neighb_high_pm(nb) ! This gives a + or - sign
311  c1 = 1
312  c2 = 3 * c0
313  case (af_bc_dirichlet_copy)
314  c0 = 1
315  c1 = 0
316  c2 = c0
317  case default
318  stop "fill_bc: unknown boundary condition"
319  end select
320 
321  select case (nb)
322 #if NDIM == 1
323  case (af_neighb_lowx)
324  cc(0) = c0 * bc_val(1) + c1 * cc(1)
325  cc(-1) = c2 * bc_val(1) + c1 * cc(2)
326  case (af_neighb_highx)
327  cc(nc+1) = c0 * bc_val(1) + c1 * cc(nc)
328  cc(nc+2) = c2 * bc_val(1) + c1 * cc(nc-1)
329 #elif NDIM == 2
330  case (af_neighb_lowx)
331  cc(0, 1:nc) = c0 * bc_val + c1 * cc(1, 1:nc)
332  cc(-1, 1:nc) = c2 * bc_val + c1 * cc(2, 1:nc)
333  case (af_neighb_highx)
334  cc(nc+1, 1:nc) = c0 * bc_val + c1 * cc(nc, 1:nc)
335  cc(nc+2, 1:nc) = c2 * bc_val + c1 * cc(nc-1, 1:nc)
336  case (af_neighb_lowy)
337  cc(1:nc, 0) = c0 * bc_val + c1 * cc(1:nc, 1)
338  cc(1:nc, -1) = c2 * bc_val + c1 * cc(1:nc, 2)
339  case (af_neighb_highy)
340  cc(1:nc, nc+1) = c0 * bc_val + c1 * cc(1:nc, nc)
341  cc(1:nc, nc+2) = c2 * bc_val + c1 * cc(1:nc, nc-1)
342 #elif NDIM == 3
343  case (af_neighb_lowx)
344  cc(0, 1:nc, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
345  c1 * cc(1, 1:nc, 1:nc)
346  cc(-1, 1:nc, 1:nc) = c2 * reshape(bc_val, [nc, nc]) + &
347  c1 * cc(2, 1:nc, 1:nc)
348  case (af_neighb_highx)
349  cc(nc+1, 1:nc, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
350  c1 * cc(nc, 1:nc, 1:nc)
351  cc(nc+2, 1:nc, 1:nc) = c2 * reshape(bc_val, [nc, nc]) + &
352  c1 * cc(nc-1, 1:nc, 1:nc)
353  case (af_neighb_lowy)
354  cc(1:nc, 0, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
355  c1 * cc(1:nc, 1, 1:nc)
356  cc(1:nc, -1, 1:nc) = c2 * reshape(bc_val, [nc, nc]) + &
357  c1 * cc(1:nc, 2, 1:nc)
358  case (af_neighb_highy)
359  cc(1:nc, nc+1, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
360  c1 * cc(1:nc, nc, 1:nc)
361  cc(1:nc, nc+2, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
362  c1 * cc(1:nc, nc-1, 1:nc)
363  case (af_neighb_lowz)
364  cc(1:nc, 1:nc, 0) = c0 * reshape(bc_val, [nc, nc]) + &
365  c1 * cc(1:nc, 1:nc, 1)
366  cc(1:nc, 1:nc, -1) = c2 * reshape(bc_val, [nc, nc]) + &
367  c1 * cc(1:nc, 1:nc, 2)
368  case (af_neighb_highz)
369  cc(1:nc, 1:nc, nc+1) = c0 * reshape(bc_val, [nc, nc]) + &
370  c1 * cc(1:nc, 1:nc, nc)
371  cc(1:nc, 1:nc, nc+2) = c2 * reshape(bc_val, [nc, nc]) + &
372  c1 * cc(1:nc, 1:nc, nc-1)
373 #endif
374  end select
375  end subroutine bc_to_gc2
376 
377  !> Partial prolongation to the ghost cells of box id from parent
378  subroutine af_gc_prolong_copy(boxes, id, nb, iv, op_mask)
379  use m_af_prolong, only: af_prolong_copy
380  type(box_t), intent(inout) :: boxes(:) !< List of all boxes
381  integer, intent(in) :: id !< Id of child
382  integer, intent(in) :: iv !< Variable to fill
383  integer, intent(in) :: nb !< Neighbor to get data from
384  integer, intent(in) :: op_mask !< Multigrid operator mask
385  integer :: p_id, lo(ndim), hi(ndim)
386 
387  p_id = boxes(id)%parent
388  call af_get_index_bc_outside(nb, boxes(id)%n_cell, 1, lo, hi)
389  call af_prolong_copy(boxes(p_id), boxes(id), iv, low=lo, high=hi)
390  end subroutine af_gc_prolong_copy
391 
392  !> Interpolate between fine points and coarse neighbors to fill ghost cells
393  !> near refinement boundaries
394  subroutine af_gc_interp(boxes, id, nb, iv, op_mask)
395  type(box_t), intent(inout) :: boxes(:) !< List of all boxes
396  integer, intent(in) :: id !< Id of box
397  integer, intent(in) :: nb !< Ghost cell direction
398  integer, intent(in) :: iv !< Ghost cell variable
399  integer, intent(in) :: op_mask !< Multigrid operator mask
400  integer :: nc, ix, ix_c, ix_f
401  integer :: p_nb_id
402  integer :: p_id, ix_offset(ndim)
403  real(dp), parameter :: third = 1/3.0_dp
404 #if NDIM > 1
405  real(dp), parameter :: sixth = 1/6.0_dp
406  integer :: i_c1, i_c2, j_c1, j_c2, i, j
407 #endif
408 #if NDIM > 2
409  integer :: k_c1, k_c2, k
410 #endif
411 
412  nc = boxes(id)%n_cell
413  p_id = boxes(id)%parent
414  p_nb_id = boxes(p_id)%neighbors(nb)
415  ix_offset = af_get_child_offset(boxes(id), nb)
416 
417  if (af_neighb_low(nb)) then
418  ix = 0
419  ix_f = 1
420  ix_c = nc
421  else
422  ix = nc+1
423  ix_f = nc
424  ix_c = 1
425  end if
426 
427  select case (af_neighb_dim(nb))
428 #if NDIM == 1
429  case (1)
430  boxes(id)%cc(ix, iv) = &
431  2 * third * boxes(p_nb_id)%cc(ix_c, iv) + &
432  third * boxes(id)%cc(ix_f, iv)
433 #elif NDIM == 2
434  case (1)
435  do j = 1, nc
436  j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
437  j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
438  boxes(id)%cc(ix, j, iv) = &
439  0.5_dp * boxes(p_nb_id)%cc(ix_c, j_c1, iv) + &
440  sixth * boxes(p_nb_id)%cc(ix_c, j_c2, iv) + &
441  third * boxes(id)%cc(ix_f, j, iv)
442  end do
443  case (2)
444  do i = 1, nc
445  i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
446  i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
447  boxes(id)%cc(i, ix, iv) = &
448  0.5_dp * boxes(p_nb_id)%cc(i_c1, ix_c, iv) + &
449  sixth * boxes(p_nb_id)%cc(i_c2, ix_c, iv) + &
450  third * boxes(id)%cc(i, ix_f, iv)
451  end do
452 #elif NDIM==3
453  case (1)
454  do k = 1, nc
455  k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
456  k_c2 = k_c1 + 1 - 2 * iand(k, 1) ! even: +1, odd: -1
457  do j = 1, nc
458  j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
459  j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
460  boxes(id)%cc(ix, j, k, iv) = &
461  third * boxes(p_nb_id)%cc(ix_c, j_c1, k_c1, iv) + &
462  sixth * boxes(p_nb_id)%cc(ix_c, j_c2, k_c1, iv) + &
463  sixth * boxes(p_nb_id)%cc(ix_c, j_c1, k_c2, iv) + &
464  third * boxes(id)%cc(ix_f, j, k, iv)
465  end do
466  end do
467  case (2)
468  do k = 1, nc
469  k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
470  k_c2 = k_c1 + 1 - 2 * iand(k, 1) ! even: +1, odd: -1
471  do i = 1, nc
472  i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
473  i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
474  boxes(id)%cc(i, ix, k, iv) = &
475  third * boxes(p_nb_id)%cc(i_c1, ix_c, k_c1, iv) + &
476  sixth * boxes(p_nb_id)%cc(i_c2, ix_c, k_c1, iv) + &
477  sixth * boxes(p_nb_id)%cc(i_c1, ix_c, k_c2, iv) + &
478  third * boxes(id)%cc(i, ix_f, k, iv)
479  end do
480  end do
481  case (3)
482  do j = 1, nc
483  j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
484  j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
485  do i = 1, nc
486  i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
487  i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
488  boxes(id)%cc(i, j, ix, iv) = &
489  third * boxes(p_nb_id)%cc(i_c1, j_c1, ix_c, iv) + &
490  sixth * boxes(p_nb_id)%cc(i_c1, j_c2, ix_c, iv) + &
491  sixth * boxes(p_nb_id)%cc(i_c2, j_c1, ix_c, iv) + &
492  third * boxes(id)%cc(i, j, ix_f, iv)
493  end do
494  end do
495 #endif
496  end select
497 
498  end subroutine af_gc_interp
499 
500  !> Interpolate between fine points and coarse neighbors to fill ghost cells
501  !> near refinement boundaries. The ghost values are less than twice the coarse
502  !> values.
503  subroutine af_gc_interp_lim(boxes, id, nb, iv, op_mask)
504  type(box_t), intent(inout) :: boxes(:) !< List of all boxes
505  integer, intent(in) :: id !< Id of box
506  integer, intent(in) :: nb !< Ghost cell direction
507  integer, intent(in) :: iv !< Ghost cell variable
508  integer, intent(in) :: op_mask !< Multigrid operator mask
509  integer :: nc, ix, ix_c, ix_f
510  integer :: p_id, ix_offset(ndim), p_nb_id
511  real(dp) :: c(ndim)
512  real(dp), parameter :: third = 1/3.0_dp
513 #if NDIM > 1
514  integer :: ijk, ijk_(c1), ijk_(c2)
515  real(dp), parameter :: sixth = 1/6.0_dp
516 #endif
517 
518  nc = boxes(id)%n_cell
519  p_id = boxes(id)%parent
520  p_nb_id = boxes(p_id)%neighbors(nb)
521  ix_offset = af_get_child_offset(boxes(id), nb)
522 
523  if (af_neighb_low(nb)) then
524  ix = 0
525  ix_f = 1
526  ix_c = nc
527  else
528  ix = nc+1
529  ix_f = nc
530  ix_c = 1
531  end if
532 
533  select case (af_neighb_dim(nb))
534 #if NDIM == 1
535  case (1)
536  c(1) = boxes(p_nb_id)%cc(ix_c, iv)
537  boxes(id)%cc(ix, iv) = (2 * c(1) + boxes(id)%cc(ix_f, iv)) * third
538  if (boxes(id)%cc(ix, iv) > 2 * c(1)) boxes(id)%cc(ix, iv) = 2 * c(1)
539 #elif NDIM == 2
540  case (1)
541  do j = 1, nc
542  j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
543  j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
544  c(1) = boxes(p_nb_id)%cc(ix_c, j_c1, iv)
545  c(2) = boxes(p_nb_id)%cc(ix_c, j_c2, iv)
546  boxes(id)%cc(ix, j, iv) = 0.5_dp * c(1) + sixth * c(2) + &
547  third * boxes(id)%cc(ix_f, j, iv)
548  if (boxes(id)%cc(ix, j, iv) > 2 * c(1)) boxes(id)%cc(ix, j, iv) = 2 * c(1)
549  end do
550  case (2)
551  do i = 1, nc
552  i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
553  i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
554  c(1) = boxes(p_nb_id)%cc(i_c1, ix_c, iv)
555  c(2) = boxes(p_nb_id)%cc(i_c2, ix_c, iv)
556  boxes(id)%cc(i, ix, iv) = 0.5_dp * c(1) + sixth * c(2) + &
557  third * boxes(id)%cc(i, ix_f, iv)
558  if (boxes(id)%cc(i, ix, iv) > 2 * c(1)) boxes(id)%cc(i, ix, iv) = 2 * c(1)
559  end do
560 #elif NDIM==3
561  case (1)
562  do k = 1, nc
563  k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
564  k_c2 = k_c1 + 1 - 2 * iand(k, 1) ! even: +1, odd: -1
565  do j = 1, nc
566  j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
567  j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
568  c(1) = boxes(p_nb_id)%cc(ix_c, j_c1, k_c1, iv)
569  c(2) = boxes(p_nb_id)%cc(ix_c, j_c2, k_c1, iv)
570  c(3) = boxes(p_nb_id)%cc(ix_c, j_c1, k_c2, iv)
571  boxes(id)%cc(ix, j, k, iv) = third * c(1) + sixth * c(2) + &
572  sixth * c(3) + third * boxes(id)%cc(ix_f, j, k, iv)
573  if (boxes(id)%cc(ix, j, k, iv) > 2 * c(1)) &
574  boxes(id)%cc(ix, j, k, iv) = 2 * c(1)
575  end do
576  end do
577  case (2)
578  do k = 1, nc
579  k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
580  k_c2 = k_c1 + 1 - 2 * iand(k, 1) ! even: +1, odd: -1
581  do i = 1, nc
582  i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
583  i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
584  c(1) = boxes(p_nb_id)%cc(i_c1, ix_c, k_c1, iv)
585  c(2) = boxes(p_nb_id)%cc(i_c2, ix_c, k_c1, iv)
586  c(3) = boxes(p_nb_id)%cc(i_c1, ix_c, k_c2, iv)
587  boxes(id)%cc(i, ix, k, iv) = third * c(1) + sixth * c(2) + &
588  sixth * c(3) + third * boxes(id)%cc(i, ix_f, k, iv)
589  if (boxes(id)%cc(i, ix, k, iv) > 2 * c(1)) &
590  boxes(id)%cc(i, ix, k, iv) = 2 * c(1)
591  end do
592  end do
593  case (3)
594  do j = 1, nc
595  j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
596  j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
597  do i = 1, nc
598  i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
599  i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
600  c(1) = boxes(p_nb_id)%cc(i_c1, j_c1, ix_c, iv)
601  c(2) = boxes(p_nb_id)%cc(i_c1, j_c2, ix_c, iv)
602  c(3) = boxes(p_nb_id)%cc(i_c2, j_c1, ix_c, iv)
603  boxes(id)%cc(i, j, ix, iv) = third * c(1) + sixth * c(2) + &
604  sixth * c(3) + third * boxes(id)%cc(i, j, ix_f, iv)
605  if (boxes(id)%cc(i, j, ix, iv) > 2 * c(1)) &
606  boxes(id)%cc(i, j, ix, iv) = 2 * c(1)
607  end do
608  end do
609 #endif
610  end select
611 
612  end subroutine af_gc_interp_lim
613 
614  !> This fills ghost cells near physical boundaries using Neumann zero
615  subroutine af_bc_neumann_zero(box, nb, iv, coords, bc_val, bc_type)
616  type(box_t), intent(in) :: box
617  integer, intent(in) :: nb
618  integer, intent(in) :: iv
619  real(dp), intent(in) :: coords(ndim, box%n_cell**(ndim-1))
620  real(dp), intent(out) :: bc_val(box%n_cell**(ndim-1))
621  integer, intent(out) :: bc_type
622 
623  bc_type = af_bc_neumann
624  bc_val = 0.0_dp
625  end subroutine af_bc_neumann_zero
626 
627  !> This fills ghost cells near physical boundaries using Dirichlet zero
628  subroutine af_bc_dirichlet_zero(box, nb, iv, coords, bc_val, bc_type)
629  type(box_t), intent(in) :: box
630  integer, intent(in) :: nb
631  integer, intent(in) :: iv
632  real(dp), intent(in) :: coords(ndim, box%n_cell**(ndim-1))
633  real(dp), intent(out) :: bc_val(box%n_cell**(ndim-1))
634  integer, intent(out) :: bc_type
635 
636  bc_type = af_bc_dirichlet
637  bc_val = 0.0_dp
638  end subroutine af_bc_dirichlet_zero
639 
640  !> This fills ghost cells near physical boundaries using the same slope
641  subroutine af_bc_set_continuous(box, nb, iv, coords, bc_val, bc_type)
642  type(box_t), intent(in) :: box
643  integer, intent(in) :: nb
644  integer, intent(in) :: iv
645  real(dp), intent(in) :: coords(ndim, box%n_cell**(ndim-1))
646  real(dp), intent(out) :: bc_val(box%n_cell**(ndim-1))
647  integer, intent(out) :: bc_type
648 
649  bc_type = af_bc_continuous
650  ! Set values to zero (to prevent problems with NaN)
651  bc_val = 0.0_dp
652  end subroutine af_bc_set_continuous
653 
654  subroutine copy_from_nb(box, box_nb, dnb, lo, hi, ivs)
655  type(box_t), intent(inout) :: box !< Box on which to fill ghost cells
656  type(box_t), intent(in) :: box_nb !< Neighbouring box
657  integer, intent(in) :: dnb(NDIM) !< Neighbor spatial index offset
658  integer, intent(in) :: lo(NDIM) !< Ghost cell low index
659  integer, intent(in) :: hi(NDIM) !< Ghost cell high index
660  integer, intent(in) :: ivs(:) !< Ghost cell variable
661  integer :: nlo(NDIM), nhi(NDIM)
662 
663  ! Get indices on neighbor
664  nlo = lo - dnb * box%n_cell
665  nhi = hi - dnb * box%n_cell
666 
667  box%cc(dslice(lo, hi), ivs) = &
668  box_nb%cc(dslice(nlo, nhi), ivs)
669  end subroutine copy_from_nb
670 
671  !> Get array of cell-centered variables with multiple ghost cells, excluding corners
672  subroutine af_gc2_box(tree, id, ivs, cc)
673  type(af_t), intent(inout) :: tree !< Tree to fill ghost cells on
674  integer, intent(in) :: id !< Id of box for which we set ghost cells
675  integer, intent(in) :: ivs(:) !< Variables for which ghost cells are set
676  real(dp), intent(inout) :: cc(dtimes(-1:tree%n_cell+2), size(ivs))
677  integer :: i, iv, nb, nc, nb_id, bc_type, ix
678  integer :: lo(ndim), hi(ndim), dnb(ndim)
679  integer :: nlo(ndim), nhi(ndim)
680  real(dp) :: bc_val(tree%n_cell**(ndim-1))
681 
682  nc = tree%n_cell
683 
684  do i = 1, size(ivs)
685  iv = ivs(i)
686 
687  ! Copy interior region
688  cc(dtimes(0:nc+1), i) = tree%boxes(id)%cc(dtimes(:), iv)
689 
690  if (.not. tree%has_cc_method(iv)) then
691  print *, "For variable ", trim(tree%cc_names(iv))
692  error stop "No methods stored"
693  end if
694  end do
695 
696  do nb = 1, af_num_neighbors
697  nb_id = tree%boxes(id)%neighbors(nb)
698 
699  if (nb_id > af_no_box) then
700  ! There is a neighbor
701  call af_get_index_bc_outside(nb, tree%n_cell, 2, lo, hi)
702 
703  dnb = af_neighb_offset([nb])
704  nlo = lo - dnb * tree%n_cell
705  nhi = hi - dnb * tree%n_cell
706 
707  cc(dslice(lo, hi), :) = &
708  tree%boxes(nb_id)%cc(dslice(nlo, nhi), ivs)
709 
710  else if (nb_id == af_no_box) then
711  ! Refinement boundary
712  do i = 1, size(ivs)
713  iv = ivs(i)
714  call gc2_prolong_rb(tree, id, nb, iv, tree%n_cell, &
715  cc(dtimes(:), i))
716  end do
717  else
718  ! Physical boundary
719  do i = 1, size(ivs)
720  iv = ivs(i)
721  if (associated(tree%cc_methods(iv)%bc_custom)) then
722  call tree%cc_methods(iv)%bc_custom(tree%boxes(id), &
723  nb, iv, 2, cc(dtimes(:), i))
724  else
725  ix = tree%boxes(id)%nb_to_bc_index(nb)
726  call tree%cc_methods(iv)%bc(tree%boxes(id), nb, iv, &
727  tree%boxes(id)%bc_coords(:, :, ix), bc_val, bc_type)
728  call bc_to_gc2(nc, cc(dtimes(:), i), nb, bc_val, bc_type, &
729  tree%boxes(id)%dr)
730  tree%boxes(id)%bc_val(:, iv, ix) = bc_val
731  tree%boxes(id)%bc_type(iv, ix) = bc_type
732  end if
733  end do
734  end if
735  end do
736 
737  do i = 1, size(ivs)
738  iv = ivs(i)
739 
740  ! Copy back updated ghost cells
741  tree%boxes(id)%cc(dtimes(:), iv) = cc(dtimes(0:nc+1), i)
742  end do
743 
744  end subroutine af_gc2_box
745 
746  !> Conservative prolongation with a limited slope for ghost cells near
747  !> refinement boundaries
748  !>
749  !> This method assumes that the ghost cells for leaves on a coarser level have
750  !> already been set, which is typically the case.
751  !>
752  !> @todo make compatible with arbitrary number of ghost cells
753  subroutine gc2_prolong_rb(tree, id, nb, iv, nc, cc)
754  use m_af_limiters
755  type(af_t), intent(in) :: tree !< Tree
756  integer, intent(in) :: id !< Index of box to fill ghost cells for
757  integer, intent(in) :: nb !< Neighbor direction
758  integer, intent(in) :: iv !< Index of variable
759  integer, intent(in) :: nc !< Number of cells
760  real(dp), intent(inout) :: cc(DTIMES(-1:nc+2)) !< Enlarged array
761  integer :: IJK, IJK_(f), p_id, p_nb_id
762  integer :: lo_c(NDIM), hi_c(NDIM), ix_offset(NDIM)
763  integer :: lo(NDIM), hi(NDIM)
764  real(dp) :: f(0:NDIM), slopes_a(NDIM), slopes_b(NDIM)
765 
766  p_id = tree%boxes(id)%parent
767  p_nb_id = tree%boxes(p_id)%neighbors(nb)
768 
769  ! Get index in current box
770  call af_get_index_bc_outside(nb, nc, 2, lo, hi)
771 
772  ! Convert to index on parent box
773  ix_offset = af_get_child_offset(tree%boxes(id))
774  lo_c = ix_offset + (lo+1)/2
775  hi_c = ix_offset + (hi+1)/2
776 
777  ! Convert to index on neighbor of parent box. These indices should then lie
778  ! in the range [1, nc]
779  lo_c = lo_c - af_neighb_dix(:, nb) * nc
780  hi_c = hi_c - af_neighb_dix(:, nb) * nc
781 
782  associate(cc_p => tree%boxes(p_nb_id)%cc, &
783  limiter => tree%cc_methods(iv)%prolong_limiter)
784 #if NDIM == 1
785  do i = lo_c(1), hi_c(1)
786  i_f = lo(1) + 2 * (i - lo_c(1))
787 
788  ! Compute slopes on parent
789  f(0) = cc_p(i, iv)
790  slopes_a = [cc_p(i, iv) - cc_p(i-1, iv)]
791  slopes_b = [cc_p(i+1, iv) - cc_p(i, iv)]
792  f(1:1) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
793 
794  ! Prolong to fine cells
795  cc(i_f) = f(0) - f(1)
796  cc(i_f+1) = f(0) + f(1)
797  end do
798 #elif NDIM == 2
799  do j = lo_c(2), hi_c(2)
800  j_f = lo(2) + 2 * (j - lo_c(2))
801  do i = lo_c(1), hi_c(1)
802  i_f = lo(1) + 2 * (i - lo_c(1))
803 
804  ! Compute slopes on parent
805  f(0) = cc_p(i, j, iv)
806  slopes_a = [cc_p(i, j, iv) - cc_p(i-1, j, iv), &
807  cc_p(i, j, iv) - cc_p(i, j-1, iv)]
808  slopes_b = [cc_p(i+1, j, iv) - cc_p(i, j, iv), &
809  cc_p(i, j+1, iv) - cc_p(i, j, iv)]
810  f(1:2) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
811 
812  ! Prolong to fine cells
813  cc(i_f, j_f) = f(0) - f(1) - f(2)
814  cc(i_f, j_f+1) = f(0) - f(1) + f(2)
815  cc(i_f+1, j_f) = f(0) + f(1) - f(2)
816  cc(i_f+1, j_f+1) = f(0) + f(1) + f(2)
817  end do
818  end do
819 #elif NDIM == 3
820  do k = lo_c(3), hi_c(3)
821  k_f = lo(3) + 2 * (k - lo_c(3))
822  do j = lo_c(2), hi_c(2)
823  j_f = lo(2) + 2 * (j - lo_c(2))
824  do i = lo_c(1), hi_c(1)
825  i_f = lo(1) + 2 * (i - lo_c(1))
826 
827  ! Compute slopes on parent
828  f(0) = cc_p(i, j, k, iv)
829  slopes_a = [cc_p(i, j, k, iv) - cc_p(i-1, j, k, iv), &
830  cc_p(i, j, k, iv) - cc_p(i, j-1, k, iv), &
831  cc_p(i, j, k, iv) - cc_p(i, j, k-1, iv)]
832  slopes_b = [cc_p(i+1, j, k, iv) - cc_p(i, j, k, iv), &
833  cc_p(i, j+1, k, iv) - cc_p(i, j, k, iv), &
834  cc_p(i, j, k+1, iv) - cc_p(i, j, k, iv)]
835  f(1:3) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
836 
837  ! Prolong to fine cells
838  cc(i_f, j_f, k_f) = f(0) - f(1) - f(2) - f(3)
839  cc(i_f, j_f, k_f+1) = f(0) - f(1) - f(2) + f(3)
840  cc(i_f, j_f+1, k_f) = f(0) - f(1) + f(2) - f(3)
841  cc(i_f, j_f+1, k_f+1) = f(0) - f(1) + f(2) + f(3)
842  cc(i_f+1, j_f, k_f) = f(0) + f(1) - f(2) - f(3)
843  cc(i_f+1, j_f, k_f+1) = f(0) + f(1) - f(2) + f(3)
844  cc(i_f+1, j_f+1, k_f) = f(0) + f(1) + f(2) - f(3)
845  cc(i_f+1, j_f+1, k_f+1) = f(0) + f(1) + f(2) + f(3)
846  end do
847  end do
848  end do
849 #endif
850  end associate
851 
852  contains
853 
854 
855 
856  end subroutine gc2_prolong_rb
857 
858  !> This fills corner ghost cells using linear extrapolation. The ghost cells
859  !> on the sides already need to be filled.
860  subroutine af_corner_gc_extrap(box, ix, ivs)
861  type(box_t), intent(inout) :: box !< Box to fill ghost cells for
862  integer, intent(in) :: ix(NDIM) !< Cell-index of corner
863  integer, intent(in) :: ivs(:) !< Variable to fill
864  integer :: di(NDIM)
865 
866  di = 1 - 2 * iand(ix, 1) ! 0 -> di = 1, nc+1 -> di = -1
867 
868 #if NDIM == 2
869  box%cc(ix(1), ix(2), ivs) = box%cc(ix(1)+di(1), ix(2), ivs) &
870  + box%cc(ix(1), ix(2)+di(2), ivs) &
871  - box%cc(ix(1)+di(1), ix(2)+di(2), ivs)
872 #elif NDIM == 3
873  box%cc(ix(1), ix(2), ix(3), ivs) = &
874  box%cc(ix(1), ix(2)+di(2), ix(3)+di(3), ivs) + &
875  box%cc(ix(1)+di(1), ix(2), ix(3)+di(3), ivs) + &
876  box%cc(ix(1)+di(1), ix(2)+di(2), ix(3), ivs) - 2 * &
877  box%cc(ix(1)+di(1), ix(2)+di(2), ix(3)+di(3), ivs)
878 #endif
879  end subroutine af_corner_gc_extrap
880 
881 #if NDIM == 3
882  !> This fills edge ghost cells using linear extrapolation. The ghost cells on
883  !> the sides already need to be filled. This routine basically performs the
884  !> same operation as af_corner_gc_extrap does in 2D.
885  subroutine af_edge_gc_extrap(box, lo, dim, ivs)
886  type(box_t), intent(inout) :: box !< Box to operate on
887  integer, intent(in) :: lo(NDIM) !< Lowest index of edge ghost cells
888  integer, intent(in) :: dim !< Dimension parallel to edge
889  integer, intent(in) :: ivs(:) !< Variable to fill
890  integer :: di(NDIM), ix(NDIM), ia(NDIM), ib(NDIM), ic(NDIM)
891  integer :: n, o_dims(NDIM-1)
892 
893  ! Dimensions other than/perpendicular to dim
894  o_dims = [1 + mod(dim, ndim), 1 + mod(dim + 1, ndim)]
895 
896  ! Index offsets
897  di = 1 - 2 * iand(lo, 1) ! 0 -> di = 1, nc+1 -> di = -1
898  di(dim) = 0
899 
900  ! Neighbor index in direction one
901  ia = lo
902  ia(o_dims(1)) = ia(o_dims(1)) + di(o_dims(1))
903 
904  ! Neighbor index in direction two
905  ib = lo
906  ib(o_dims(2)) = ib(o_dims(2)) + di(o_dims(2))
907 
908  ! Diagional neighbor index
909  ic = lo + di
910  ix = lo
911 
912  do n = 1, box%n_cell
913  ia(dim) = n
914  ib(dim) = n
915  ic(dim) = n
916  ix(dim) = n
917 
918  box%cc(ix(1), ix(2), ix(3), ivs) = &
919  box%cc(ia(1), ia(2), ia(3), ivs) + &
920  box%cc(ib(1), ib(2), ib(3), ivs) - &
921  box%cc(ic(1), ic(2), ic(3), ivs)
922  end do
923 
924  end subroutine af_edge_gc_extrap
925 #endif
926 
927 end module m_af_ghostcell
m_af_ghostcell
This module contains routines related to the filling of ghost cells. Note that corner ghost cells are...
Definition: m_af_ghostcell.f90:3
m_af_limiters
Module containing slope limiters.
Definition: m_af_limiters.f90:2
m_af_prolong
This module contains the routines related to prolongation: going from coarse to fine variables.
Definition: m_af_prolong.f90:3
m_af_types
This module contains the basic types and constants that are used in the NDIM-dimensional version of A...
Definition: m_af_types.f90:3
m_af_types::af_t
Type which stores all the boxes and levels, as well as some information about the number of boxes,...
Definition: m_af_types.f90:326
m_af_types::box_t
The basic building block of afivo: a box with cell-centered and face centered data,...
Definition: m_af_types.f90:286
Generated on Sat Jan 11 2025 05:02:04 for Afivo by doxygen 1.9.1