Afivo 0.3
•All Classes Namespaces Functions Variables Pages
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.
3module 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
21contains
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 ! Below, we set coefficients to handle these cases
300 select case (bc_type)
301 case (af_bc_dirichlet)
302 c0 = 2
303 c1 = -1
304 c2 = c0
305 case (af_bc_neumann)
306 c0 = dr(af_neighb_dim(nb)) * af_neighb_high_pm(nb) ! This gives a + or - sign
307 c1 = 1
308 c2 = 3 * c0
309 case (af_bc_dirichlet_copy)
310 c0 = 1
311 c1 = 0
312 c2 = c0
313 case default
314 stop "fill_bc: unknown boundary condition"
315 end select
316
317 select case (nb)
318#if NDIM == 1
319 case (af_neighb_lowx)
320 cc(0) = c0 * bc_val(1) + c1 * cc(1)
321 cc(-1) = c2 * bc_val(1) + c1 * cc(2)
322 case (af_neighb_highx)
323 cc(nc+1) = c0 * bc_val(1) + c1 * cc(nc)
324 cc(nc+2) = c2 * bc_val(1) + c1 * cc(nc-1)
325#elif NDIM == 2
326 case (af_neighb_lowx)
327 cc(0, 1:nc) = c0 * bc_val + c1 * cc(1, 1:nc)
328 cc(-1, 1:nc) = c2 * bc_val + c1 * cc(2, 1:nc)
329 case (af_neighb_highx)
330 cc(nc+1, 1:nc) = c0 * bc_val + c1 * cc(nc, 1:nc)
331 cc(nc+2, 1:nc) = c2 * bc_val + c1 * cc(nc-1, 1:nc)
332 case (af_neighb_lowy)
333 cc(1:nc, 0) = c0 * bc_val + c1 * cc(1:nc, 1)
334 cc(1:nc, -1) = c2 * bc_val + c1 * cc(1:nc, 2)
335 case (af_neighb_highy)
336 cc(1:nc, nc+1) = c0 * bc_val + c1 * cc(1:nc, nc)
337 cc(1:nc, nc+2) = c2 * bc_val + c1 * cc(1:nc, nc-1)
338#elif NDIM == 3
339 case (af_neighb_lowx)
340 cc(0, 1:nc, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
341 c1 * cc(1, 1:nc, 1:nc)
342 cc(-1, 1:nc, 1:nc) = c2 * reshape(bc_val, [nc, nc]) + &
343 c1 * cc(2, 1:nc, 1:nc)
344 case (af_neighb_highx)
345 cc(nc+1, 1:nc, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
346 c1 * cc(nc, 1:nc, 1:nc)
347 cc(nc+2, 1:nc, 1:nc) = c2 * reshape(bc_val, [nc, nc]) + &
348 c1 * cc(nc-1, 1:nc, 1:nc)
349 case (af_neighb_lowy)
350 cc(1:nc, 0, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
351 c1 * cc(1:nc, 1, 1:nc)
352 cc(1:nc, -1, 1:nc) = c2 * reshape(bc_val, [nc, nc]) + &
353 c1 * cc(1:nc, 2, 1:nc)
354 case (af_neighb_highy)
355 cc(1:nc, nc+1, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
356 c1 * cc(1:nc, nc, 1:nc)
357 cc(1:nc, nc+2, 1:nc) = c0 * reshape(bc_val, [nc, nc]) + &
358 c1 * cc(1:nc, nc-1, 1:nc)
359 case (af_neighb_lowz)
360 cc(1:nc, 1:nc, 0) = c0 * reshape(bc_val, [nc, nc]) + &
361 c1 * cc(1:nc, 1:nc, 1)
362 cc(1:nc, 1:nc, -1) = c2 * reshape(bc_val, [nc, nc]) + &
363 c1 * cc(1:nc, 1:nc, 2)
364 case (af_neighb_highz)
365 cc(1:nc, 1:nc, nc+1) = c0 * reshape(bc_val, [nc, nc]) + &
366 c1 * cc(1:nc, 1:nc, nc)
367 cc(1:nc, 1:nc, nc+2) = c2 * reshape(bc_val, [nc, nc]) + &
368 c1 * cc(1:nc, 1:nc, nc-1)
369#endif
370 end select
371 end subroutine bc_to_gc2
372
373 !> Partial prolongation to the ghost cells of box id from parent
374 subroutine af_gc_prolong_copy(boxes, id, nb, iv, op_mask)
375 use m_af_prolong, only: af_prolong_copy
376 type(box_t), intent(inout) :: boxes(:) !< List of all boxes
377 integer, intent(in) :: id !< Id of child
378 integer, intent(in) :: iv !< Variable to fill
379 integer, intent(in) :: nb !< Neighbor to get data from
380 integer, intent(in) :: op_mask !< Multigrid operator mask
381 integer :: p_id, lo(ndim), hi(ndim)
382
383 p_id = boxes(id)%parent
384 call af_get_index_bc_outside(nb, boxes(id)%n_cell, 1, lo, hi)
385 call af_prolong_copy(boxes(p_id), boxes(id), iv, low=lo, high=hi)
386 end subroutine af_gc_prolong_copy
387
388 !> Interpolate between fine points and coarse neighbors to fill ghost cells
389 !> near refinement boundaries
390 subroutine af_gc_interp(boxes, id, nb, iv, op_mask)
391 type(box_t), intent(inout) :: boxes(:) !< List of all boxes
392 integer, intent(in) :: id !< Id of box
393 integer, intent(in) :: nb !< Ghost cell direction
394 integer, intent(in) :: iv !< Ghost cell variable
395 integer, intent(in) :: op_mask !< Multigrid operator mask
396 integer :: nc, ix, ix_c, ix_f
397 integer :: p_nb_id
398 integer :: p_id, ix_offset(ndim)
399 real(dp), parameter :: third = 1/3.0_dp
400#if NDIM > 1
401 real(dp), parameter :: sixth = 1/6.0_dp
402 integer :: i_c1, i_c2, j_c1, j_c2, i, j
403#endif
404#if NDIM > 2
405 integer :: k_c1, k_c2, k
406#endif
407
408 nc = boxes(id)%n_cell
409 p_id = boxes(id)%parent
410 p_nb_id = boxes(p_id)%neighbors(nb)
411 ix_offset = af_get_child_offset(boxes(id), nb)
412
413 if (af_neighb_low(nb)) then
414 ix = 0
415 ix_f = 1
416 ix_c = nc
417 else
418 ix = nc+1
419 ix_f = nc
420 ix_c = 1
421 end if
422
423 select case (af_neighb_dim(nb))
424#if NDIM == 1
425 case (1)
426 boxes(id)%cc(ix, iv) = &
427 2 * third * boxes(p_nb_id)%cc(ix_c, iv) + &
428 third * boxes(id)%cc(ix_f, iv)
429#elif NDIM == 2
430 case (1)
431 do j = 1, nc
432 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
433 j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
434 boxes(id)%cc(ix, j, iv) = &
435 0.5_dp * boxes(p_nb_id)%cc(ix_c, j_c1, iv) + &
436 sixth * boxes(p_nb_id)%cc(ix_c, j_c2, iv) + &
437 third * boxes(id)%cc(ix_f, j, iv)
438 end do
439 case (2)
440 do i = 1, nc
441 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
442 i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
443 boxes(id)%cc(i, ix, iv) = &
444 0.5_dp * boxes(p_nb_id)%cc(i_c1, ix_c, iv) + &
445 sixth * boxes(p_nb_id)%cc(i_c2, ix_c, iv) + &
446 third * boxes(id)%cc(i, ix_f, iv)
447 end do
448#elif NDIM==3
449 case (1)
450 do k = 1, nc
451 k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
452 k_c2 = k_c1 + 1 - 2 * iand(k, 1) ! even: +1, odd: -1
453 do j = 1, nc
454 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
455 j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
456 boxes(id)%cc(ix, j, k, iv) = &
457 third * boxes(p_nb_id)%cc(ix_c, j_c1, k_c1, iv) + &
458 sixth * boxes(p_nb_id)%cc(ix_c, j_c2, k_c1, iv) + &
459 sixth * boxes(p_nb_id)%cc(ix_c, j_c1, k_c2, iv) + &
460 third * boxes(id)%cc(ix_f, j, k, iv)
461 end do
462 end do
463 case (2)
464 do k = 1, nc
465 k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
466 k_c2 = k_c1 + 1 - 2 * iand(k, 1) ! even: +1, odd: -1
467 do i = 1, nc
468 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
469 i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
470 boxes(id)%cc(i, ix, k, iv) = &
471 third * boxes(p_nb_id)%cc(i_c1, ix_c, k_c1, iv) + &
472 sixth * boxes(p_nb_id)%cc(i_c2, ix_c, k_c1, iv) + &
473 sixth * boxes(p_nb_id)%cc(i_c1, ix_c, k_c2, iv) + &
474 third * boxes(id)%cc(i, ix_f, k, iv)
475 end do
476 end do
477 case (3)
478 do j = 1, nc
479 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
480 j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
481 do i = 1, nc
482 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
483 i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
484 boxes(id)%cc(i, j, ix, iv) = &
485 third * boxes(p_nb_id)%cc(i_c1, j_c1, ix_c, iv) + &
486 sixth * boxes(p_nb_id)%cc(i_c1, j_c2, ix_c, iv) + &
487 sixth * boxes(p_nb_id)%cc(i_c2, j_c1, ix_c, iv) + &
488 third * boxes(id)%cc(i, j, ix_f, iv)
489 end do
490 end do
491#endif
492 end select
493
494 end subroutine af_gc_interp
495
496 !> Interpolate between fine points and coarse neighbors to fill ghost cells
497 !> near refinement boundaries. The ghost values are less than twice the coarse
498 !> values.
499 subroutine af_gc_interp_lim(boxes, id, nb, iv, op_mask)
500 type(box_t), intent(inout) :: boxes(:) !< List of all boxes
501 integer, intent(in) :: id !< Id of box
502 integer, intent(in) :: nb !< Ghost cell direction
503 integer, intent(in) :: iv !< Ghost cell variable
504 integer, intent(in) :: op_mask !< Multigrid operator mask
505 integer :: nc, ix, ix_c, ix_f
506 integer :: p_id, ix_offset(ndim), p_nb_id
507 real(dp) :: c(ndim)
508 real(dp), parameter :: third = 1/3.0_dp
509#if NDIM > 1
510 integer :: ijk, ijk_(c1), ijk_(c2)
511 real(dp), parameter :: sixth = 1/6.0_dp
512#endif
513
514 nc = boxes(id)%n_cell
515 p_id = boxes(id)%parent
516 p_nb_id = boxes(p_id)%neighbors(nb)
517 ix_offset = af_get_child_offset(boxes(id), nb)
518
519 if (af_neighb_low(nb)) then
520 ix = 0
521 ix_f = 1
522 ix_c = nc
523 else
524 ix = nc+1
525 ix_f = nc
526 ix_c = 1
527 end if
528
529 select case (af_neighb_dim(nb))
530#if NDIM == 1
531 case (1)
532 c(1) = boxes(p_nb_id)%cc(ix_c, iv)
533 boxes(id)%cc(ix, iv) = (2 * c(1) + boxes(id)%cc(ix_f, iv)) * third
534 if (boxes(id)%cc(ix, iv) > 2 * c(1)) boxes(id)%cc(ix, iv) = 2 * c(1)
535#elif NDIM == 2
536 case (1)
537 do j = 1, nc
538 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
539 j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
540 c(1) = boxes(p_nb_id)%cc(ix_c, j_c1, iv)
541 c(2) = boxes(p_nb_id)%cc(ix_c, j_c2, iv)
542 boxes(id)%cc(ix, j, iv) = 0.5_dp * c(1) + sixth * c(2) + &
543 third * boxes(id)%cc(ix_f, j, iv)
544 if (boxes(id)%cc(ix, j, iv) > 2 * c(1)) boxes(id)%cc(ix, j, iv) = 2 * c(1)
545 end do
546 case (2)
547 do i = 1, nc
548 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
549 i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
550 c(1) = boxes(p_nb_id)%cc(i_c1, ix_c, iv)
551 c(2) = boxes(p_nb_id)%cc(i_c2, ix_c, iv)
552 boxes(id)%cc(i, ix, iv) = 0.5_dp * c(1) + sixth * c(2) + &
553 third * boxes(id)%cc(i, ix_f, iv)
554 if (boxes(id)%cc(i, ix, iv) > 2 * c(1)) boxes(id)%cc(i, ix, iv) = 2 * c(1)
555 end do
556#elif NDIM==3
557 case (1)
558 do k = 1, nc
559 k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
560 k_c2 = k_c1 + 1 - 2 * iand(k, 1) ! even: +1, odd: -1
561 do j = 1, nc
562 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
563 j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
564 c(1) = boxes(p_nb_id)%cc(ix_c, j_c1, k_c1, iv)
565 c(2) = boxes(p_nb_id)%cc(ix_c, j_c2, k_c1, iv)
566 c(3) = boxes(p_nb_id)%cc(ix_c, j_c1, k_c2, iv)
567 boxes(id)%cc(ix, j, k, iv) = third * c(1) + sixth * c(2) + &
568 sixth * c(3) + third * boxes(id)%cc(ix_f, j, k, iv)
569 if (boxes(id)%cc(ix, j, k, iv) > 2 * c(1)) &
570 boxes(id)%cc(ix, j, k, iv) = 2 * c(1)
571 end do
572 end do
573 case (2)
574 do k = 1, nc
575 k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
576 k_c2 = k_c1 + 1 - 2 * iand(k, 1) ! even: +1, odd: -1
577 do i = 1, nc
578 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
579 i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
580 c(1) = boxes(p_nb_id)%cc(i_c1, ix_c, k_c1, iv)
581 c(2) = boxes(p_nb_id)%cc(i_c2, ix_c, k_c1, iv)
582 c(3) = boxes(p_nb_id)%cc(i_c1, ix_c, k_c2, iv)
583 boxes(id)%cc(i, ix, k, iv) = third * c(1) + sixth * c(2) + &
584 sixth * c(3) + third * boxes(id)%cc(i, ix_f, k, iv)
585 if (boxes(id)%cc(i, ix, k, iv) > 2 * c(1)) &
586 boxes(id)%cc(i, ix, k, iv) = 2 * c(1)
587 end do
588 end do
589 case (3)
590 do j = 1, nc
591 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
592 j_c2 = j_c1 + 1 - 2 * iand(j, 1) ! even: +1, odd: -1
593 do i = 1, nc
594 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
595 i_c2 = i_c1 + 1 - 2 * iand(i, 1) ! even: +1, odd: -1
596 c(1) = boxes(p_nb_id)%cc(i_c1, j_c1, ix_c, iv)
597 c(2) = boxes(p_nb_id)%cc(i_c1, j_c2, ix_c, iv)
598 c(3) = boxes(p_nb_id)%cc(i_c2, j_c1, ix_c, iv)
599 boxes(id)%cc(i, j, ix, iv) = third * c(1) + sixth * c(2) + &
600 sixth * c(3) + third * boxes(id)%cc(i, j, ix_f, iv)
601 if (boxes(id)%cc(i, j, ix, iv) > 2 * c(1)) &
602 boxes(id)%cc(i, j, ix, iv) = 2 * c(1)
603 end do
604 end do
605#endif
606 end select
607
608 end subroutine af_gc_interp_lim
609
610 !> This fills ghost cells near physical boundaries using Neumann zero
611 subroutine af_bc_neumann_zero(box, nb, iv, coords, bc_val, bc_type)
612 type(box_t), intent(in) :: box
613 integer, intent(in) :: nb
614 integer, intent(in) :: iv
615 real(dp), intent(in) :: coords(ndim, box%n_cell**(ndim-1))
616 real(dp), intent(out) :: bc_val(box%n_cell**(ndim-1))
617 integer, intent(out) :: bc_type
618
619 bc_type = af_bc_neumann
620 bc_val = 0.0_dp
621 end subroutine af_bc_neumann_zero
622
623 !> This fills ghost cells near physical boundaries using Dirichlet zero
624 subroutine af_bc_dirichlet_zero(box, nb, iv, coords, bc_val, bc_type)
625 type(box_t), intent(in) :: box
626 integer, intent(in) :: nb
627 integer, intent(in) :: iv
628 real(dp), intent(in) :: coords(ndim, box%n_cell**(ndim-1))
629 real(dp), intent(out) :: bc_val(box%n_cell**(ndim-1))
630 integer, intent(out) :: bc_type
631
632 bc_type = af_bc_dirichlet
633 bc_val = 0.0_dp
634 end subroutine af_bc_dirichlet_zero
635
636 !> This fills ghost cells near physical boundaries using the same slope
637 subroutine af_bc_set_continuous(box, nb, iv, coords, bc_val, bc_type)
638 type(box_t), intent(in) :: box
639 integer, intent(in) :: nb
640 integer, intent(in) :: iv
641 real(dp), intent(in) :: coords(ndim, box%n_cell**(ndim-1))
642 real(dp), intent(out) :: bc_val(box%n_cell**(ndim-1))
643 integer, intent(out) :: bc_type
644
645 bc_type = af_bc_continuous
646 ! Set values to zero (to prevent problems with NaN)
647 bc_val = 0.0_dp
648 end subroutine af_bc_set_continuous
649
650 subroutine copy_from_nb(box, box_nb, dnb, lo, hi, ivs)
651 type(box_t), intent(inout) :: box !< Box on which to fill ghost cells
652 type(box_t), intent(in) :: box_nb !< Neighbouring box
653 integer, intent(in) :: dnb(NDIM) !< Neighbor spatial index offset
654 integer, intent(in) :: lo(NDIM) !< Ghost cell low index
655 integer, intent(in) :: hi(NDIM) !< Ghost cell high index
656 integer, intent(in) :: ivs(:) !< Ghost cell variable
657 integer :: nlo(NDIM), nhi(NDIM)
658
659 ! Get indices on neighbor
660 nlo = lo - dnb * box%n_cell
661 nhi = hi - dnb * box%n_cell
662
663 box%cc(dslice(lo, hi), ivs) = &
664 box_nb%cc(dslice(nlo, nhi), ivs)
665 end subroutine copy_from_nb
666
667 !> Get array of cell-centered variables with multiple ghost cells, excluding corners
668 subroutine af_gc2_box(tree, id, ivs, cc)
669 type(af_t), intent(inout) :: tree !< Tree to fill ghost cells on
670 integer, intent(in) :: id !< Id of box for which we set ghost cells
671 integer, intent(in) :: ivs(:) !< Variables for which ghost cells are set
672 real(dp), intent(inout) :: cc(dtimes(-1:tree%n_cell+2), size(ivs))
673 integer :: i, iv, nb, nc, nb_id, bc_type, ix
674 integer :: lo(ndim), hi(ndim), dnb(ndim)
675 integer :: nlo(ndim), nhi(ndim)
676 real(dp) :: bc_val(tree%n_cell**(ndim-1))
677
678 nc = tree%n_cell
679
680 do i = 1, size(ivs)
681 iv = ivs(i)
682
683 ! Copy interior region
684 cc(dtimes(0:nc+1), i) = tree%boxes(id)%cc(dtimes(:), iv)
685
686 if (.not. tree%has_cc_method(iv)) then
687 print *, "For variable ", trim(tree%cc_names(iv))
688 error stop "No methods stored"
689 end if
690 end do
691
692 do nb = 1, af_num_neighbors
693 nb_id = tree%boxes(id)%neighbors(nb)
694
695 if (nb_id > af_no_box) then
696 ! There is a neighbor
697 call af_get_index_bc_outside(nb, tree%n_cell, 2, lo, hi)
698
699 dnb = af_neighb_offset([nb])
700 nlo = lo - dnb * tree%n_cell
701 nhi = hi - dnb * tree%n_cell
702
703 cc(dslice(lo, hi), :) = &
704 tree%boxes(nb_id)%cc(dslice(nlo, nhi), ivs)
705
706 else if (nb_id == af_no_box) then
707 ! Refinement boundary
708 do i = 1, size(ivs)
709 iv = ivs(i)
710 call gc2_prolong_rb(tree, id, nb, iv, tree%n_cell, &
711 cc(dtimes(:), i))
712 end do
713 else
714 ! Physical boundary
715 do i = 1, size(ivs)
716 iv = ivs(i)
717 if (associated(tree%cc_methods(iv)%bc_custom)) then
718 call tree%cc_methods(iv)%bc_custom(tree%boxes(id), &
719 nb, iv, 2, cc(dtimes(:), i))
720 else
721 ix = tree%boxes(id)%nb_to_bc_index(nb)
722 call tree%cc_methods(iv)%bc(tree%boxes(id), nb, iv, &
723 tree%boxes(id)%bc_coords(:, :, ix), bc_val, bc_type)
724 call bc_to_gc2(nc, cc(dtimes(:), i), nb, bc_val, bc_type, &
725 tree%boxes(id)%dr)
726 tree%boxes(id)%bc_val(:, iv, ix) = bc_val
727 tree%boxes(id)%bc_type(iv, ix) = bc_type
728 end if
729 end do
730 end if
731 end do
732
733 do i = 1, size(ivs)
734 iv = ivs(i)
735
736 ! Copy back updated ghost cells
737 tree%boxes(id)%cc(dtimes(:), iv) = cc(dtimes(0:nc+1), i)
738 end do
739
740 end subroutine af_gc2_box
741
742 !> Conservative prolongation with a limited slope for ghost cells near
743 !> refinement boundaries
744 !>
745 !> This method assumes that the ghost cells for leaves on a coarser level have
746 !> already been set, which is typically the case.
747 !>
748 !> @todo make compatible with arbitrary number of ghost cells
749 subroutine gc2_prolong_rb(tree, id, nb, iv, nc, cc)
750 use m_af_limiters
751 type(af_t), intent(in) :: tree !< Tree
752 integer, intent(in) :: id !< Index of box to fill ghost cells for
753 integer, intent(in) :: nb !< Neighbor direction
754 integer, intent(in) :: iv !< Index of variable
755 integer, intent(in) :: nc !< Number of cells
756 real(dp), intent(inout) :: cc(DTIMES(-1:nc+2)) !< Enlarged array
757 integer :: IJK, IJK_(f), p_id, p_nb_id
758 integer :: lo_c(NDIM), hi_c(NDIM), ix_offset(NDIM)
759 integer :: lo(NDIM), hi(NDIM)
760 real(dp) :: f(0:NDIM), slopes_a(NDIM), slopes_b(NDIM)
761
762 p_id = tree%boxes(id)%parent
763 p_nb_id = tree%boxes(p_id)%neighbors(nb)
764
765 ! Get index in current box
766 call af_get_index_bc_outside(nb, nc, 2, lo, hi)
767
768 ! Convert to index on parent box
769 ix_offset = af_get_child_offset(tree%boxes(id))
770 lo_c = ix_offset + (lo+1)/2
771 hi_c = ix_offset + (hi+1)/2
772
773 ! Convert to index on neighbor of parent box. These indices should then lie
774 ! in the range [1, nc]
775 lo_c = lo_c - af_neighb_dix(:, nb) * nc
776 hi_c = hi_c - af_neighb_dix(:, nb) * nc
777
778 associate(cc_p => tree%boxes(p_nb_id)%cc, &
779 limiter => tree%cc_methods(iv)%prolong_limiter)
780#if NDIM == 1
781 do i = lo_c(1), hi_c(1)
782 i_f = lo(1) + 2 * (i - lo_c(1))
783
784 ! Compute slopes on parent
785 f(0) = cc_p(i, iv)
786 slopes_a = [cc_p(i, iv) - cc_p(i-1, iv)]
787 slopes_b = [cc_p(i+1, iv) - cc_p(i, iv)]
788 f(1:1) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
789
790 ! Prolong to fine cells
791 cc(i_f) = f(0) - f(1)
792 cc(i_f+1) = f(0) + f(1)
793 end do
794#elif NDIM == 2
795 do j = lo_c(2), hi_c(2)
796 j_f = lo(2) + 2 * (j - lo_c(2))
797 do i = lo_c(1), hi_c(1)
798 i_f = lo(1) + 2 * (i - lo_c(1))
799
800 ! Compute slopes on parent
801 f(0) = cc_p(i, j, iv)
802 slopes_a = [cc_p(i, j, iv) - cc_p(i-1, j, iv), &
803 cc_p(i, j, iv) - cc_p(i, j-1, iv)]
804 slopes_b = [cc_p(i+1, j, iv) - cc_p(i, j, iv), &
805 cc_p(i, j+1, iv) - cc_p(i, j, iv)]
806 f(1:2) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
807
808 ! Prolong to fine cells
809 cc(i_f, j_f) = f(0) - f(1) - f(2)
810 cc(i_f, j_f+1) = f(0) - f(1) + f(2)
811 cc(i_f+1, j_f) = f(0) + f(1) - f(2)
812 cc(i_f+1, j_f+1) = f(0) + f(1) + f(2)
813 end do
814 end do
815#elif NDIM == 3
816 do k = lo_c(3), hi_c(3)
817 k_f = lo(3) + 2 * (k - lo_c(3))
818 do j = lo_c(2), hi_c(2)
819 j_f = lo(2) + 2 * (j - lo_c(2))
820 do i = lo_c(1), hi_c(1)
821 i_f = lo(1) + 2 * (i - lo_c(1))
822
823 ! Compute slopes on parent
824 f(0) = cc_p(i, j, k, iv)
825 slopes_a = [cc_p(i, j, k, iv) - cc_p(i-1, j, k, iv), &
826 cc_p(i, j, k, iv) - cc_p(i, j-1, k, iv), &
827 cc_p(i, j, k, iv) - cc_p(i, j, k-1, iv)]
828 slopes_b = [cc_p(i+1, j, k, iv) - cc_p(i, j, k, iv), &
829 cc_p(i, j+1, k, iv) - cc_p(i, j, k, iv), &
830 cc_p(i, j, k+1, iv) - cc_p(i, j, k, iv)]
831 f(1:3) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
832
833 ! Prolong to fine cells
834 cc(i_f, j_f, k_f) = f(0) - f(1) - f(2) - f(3)
835 cc(i_f, j_f, k_f+1) = f(0) - f(1) - f(2) + f(3)
836 cc(i_f, j_f+1, k_f) = f(0) - f(1) + f(2) - f(3)
837 cc(i_f, j_f+1, k_f+1) = f(0) - f(1) + f(2) + f(3)
838 cc(i_f+1, j_f, k_f) = f(0) + f(1) - f(2) - f(3)
839 cc(i_f+1, j_f, k_f+1) = f(0) + f(1) - f(2) + f(3)
840 cc(i_f+1, j_f+1, k_f) = f(0) + f(1) + f(2) - f(3)
841 cc(i_f+1, j_f+1, k_f+1) = f(0) + f(1) + f(2) + f(3)
842 end do
843 end do
844 end do
845#endif
846 end associate
847
848 contains
849
850
851
852 end subroutine gc2_prolong_rb
853
854 !> This fills corner ghost cells using linear extrapolation. The ghost cells
855 !> on the sides already need to be filled.
856 subroutine af_corner_gc_extrap(box, ix, ivs)
857 type(box_t), intent(inout) :: box !< Box to fill ghost cells for
858 integer, intent(in) :: ix(NDIM) !< Cell-index of corner
859 integer, intent(in) :: ivs(:) !< Variable to fill
860 integer :: di(NDIM)
861
862 di = 1 - 2 * iand(ix, 1) ! 0 -> di = 1, nc+1 -> di = -1
863
864#if NDIM == 2
865 box%cc(ix(1), ix(2), ivs) = box%cc(ix(1)+di(1), ix(2), ivs) &
866 + box%cc(ix(1), ix(2)+di(2), ivs) &
867 - box%cc(ix(1)+di(1), ix(2)+di(2), ivs)
868#elif NDIM == 3
869 box%cc(ix(1), ix(2), ix(3), ivs) = &
870 box%cc(ix(1), ix(2)+di(2), ix(3)+di(3), ivs) + &
871 box%cc(ix(1)+di(1), ix(2), ix(3)+di(3), ivs) + &
872 box%cc(ix(1)+di(1), ix(2)+di(2), ix(3), ivs) - 2 * &
873 box%cc(ix(1)+di(1), ix(2)+di(2), ix(3)+di(3), ivs)
874#endif
875 end subroutine af_corner_gc_extrap
876
877#if NDIM == 3
878 !> This fills edge ghost cells using linear extrapolation. The ghost cells on
879 !> the sides already need to be filled. This routine basically performs the
880 !> same operation as af_corner_gc_extrap does in 2D.
881 subroutine af_edge_gc_extrap(box, lo, dim, ivs)
882 type(box_t), intent(inout) :: box !< Box to operate on
883 integer, intent(in) :: lo(NDIM) !< Lowest index of edge ghost cells
884 integer, intent(in) :: dim !< Dimension parallel to edge
885 integer, intent(in) :: ivs(:) !< Variable to fill
886 integer :: di(NDIM), ix(NDIM), ia(NDIM), ib(NDIM), ic(NDIM)
887 integer :: n, o_dims(NDIM-1)
888
889 ! Dimensions other than/perpendicular to dim
890 o_dims = [1 + mod(dim, ndim), 1 + mod(dim + 1, ndim)]
891
892 ! Index offsets
893 di = 1 - 2 * iand(lo, 1) ! 0 -> di = 1, nc+1 -> di = -1
894 di(dim) = 0
895
896 ! Neighbor index in direction one
897 ia = lo
898 ia(o_dims(1)) = ia(o_dims(1)) + di(o_dims(1))
899
900 ! Neighbor index in direction two
901 ib = lo
902 ib(o_dims(2)) = ib(o_dims(2)) + di(o_dims(2))
903
904 ! Diagional neighbor index
905 ic = lo + di
906 ix = lo
907
908 do n = 1, box%n_cell
909 ia(dim) = n
910 ib(dim) = n
911 ic(dim) = n
912 ix(dim) = n
913
914 box%cc(ix(1), ix(2), ix(3), ivs) = &
915 box%cc(ia(1), ia(2), ia(3), ivs) + &
916 box%cc(ib(1), ib(2), ib(3), ivs) - &
917 box%cc(ic(1), ic(2), ic(3), ivs)
918 end do
919
920 end subroutine af_edge_gc_extrap
921#endif
922
923end 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 Apr 19 2025 05:02:04 for Afivo by doxygen 1.9.4