Afivo 0.3
•All Classes Namespaces Functions Variables Pages
m_af_prolong.f90
1!> This module contains the routines related to prolongation: going from
2!> coarse to fine variables.
3module m_af_prolong
4#include "cpp_macros.h"
5 use m_af_types
6
7 implicit none
8 private
9
10 public :: af_prolong_copy_from
11 public :: af_prolong_copy
12 public :: af_prolong_zeroth
13 public :: af_prolong_linear_from
14 public :: af_prolong_sparse
15 public :: af_prolong_linear
16 ! public :: af_prolong_quadratic_from
17 ! public :: af_prolong_quadratic
18
19 public :: af_prolong_limit
20 public :: af_prolong_linear_cons
21
22contains
23
24 !> Zeroth-order prolongation to children.
25 subroutine af_prolong_copy_from(boxes, id, iv, iv_to, add)
26 type(box_t), intent(inout) :: boxes(:) !< List of all boxes
27 integer, intent(in) :: id !< Box whose children we will fill
28 integer, intent(in) :: iv !< Variable that is prolonged
29 integer, intent(in), optional :: iv_to !< Destination variable
30 logical, intent(in), optional :: add !< Add to old values
31 integer :: i_c, c_id
32
33 do i_c = 1, af_num_children
34 c_id = boxes(id)%children(i_c)
35 if (c_id == af_no_box) cycle
36 call af_prolong_copy(boxes(id), boxes(c_id), iv, iv_to=iv_to, add=add)
37 end do
38 end subroutine af_prolong_copy_from
39
40 !> Partial prolongation to a child (from parent) using injection (simply copy value)
41 subroutine af_prolong_copy(box_p, box_c, iv, iv_to, low, high, add)
42 type(box_t), intent(in) :: box_p !< Parent box
43 type(box_t), intent(inout) :: box_c !< Child box
44 integer, intent(in) :: iv !< Variable to fill
45 integer, intent(in), optional :: iv_to !< Destination variable
46 integer, intent(in), optional :: low(ndim) !< Min cell index at child
47 integer, intent(in), optional :: high(ndim) !< Max cell index at child
48 logical, intent(in), optional :: add !< Add to old values
49 logical :: add_to
50 integer :: nc, ix_offset(ndim), ivc
51 integer :: ijk, ijk_(c1), lo(ndim), hi(ndim)
52
53 nc = box_c%n_cell
54 add_to = .false.; if (present(add)) add_to = add
55 ivc = iv; if (present(iv_to)) ivc = iv_to
56 lo = 1; if (present(low)) lo = low
57 hi = nc; if (present(high)) hi = high
58
59 ! Offset of child w.r.t. parent
60 ix_offset = af_get_child_offset(box_c)
61
62 if (add_to) then
63#if NDIM == 1
64 do i = lo(1), hi(1)
65 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
66 box_c%cc(i, ivc) = box_c%cc(i, ivc) + &
67 box_p%cc(i_c1, iv)
68 end do
69#elif NDIM == 2
70 do j = lo(2), hi(2)
71 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
72 do i = lo(1), hi(1)
73 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
74 box_c%cc(i, j, ivc) = box_c%cc(i, j, ivc) + &
75 box_p%cc(i_c1, j_c1, iv)
76 end do
77 end do
78#elif NDIM == 3
79 do k = lo(3), hi(3)
80 k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
81 do j = lo(2), hi(2)
82 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
83 do i = lo(1), hi(1)
84 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
85 box_c%cc(i, j, k, ivc) = box_c%cc(i, j, k, ivc) + &
86 box_p%cc(i_c1, j_c1, k_c1, iv)
87 end do
88 end do
89 end do
90#endif
91 else
92#if NDIM == 1
93 do i = lo(1), hi(1)
94 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
95 box_c%cc(i, ivc) = box_p%cc(i_c1, iv)
96 end do
97#elif NDIM == 2
98 do j = lo(2), hi(2)
99 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
100 do i = lo(1), hi(1)
101 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
102 box_c%cc(i, j, ivc) = box_p%cc(i_c1, j_c1, iv)
103 end do
104 end do
105#elif NDIM == 3
106 do k = lo(3), hi(3)
107 k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
108 do j = lo(2), hi(2)
109 j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
110 do i = lo(1), hi(1)
111 i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
112 box_c%cc(i, j, k, ivc) = box_p%cc(i_c1, j_c1, k_c1, iv)
113 end do
114 end do
115 end do
116#endif
117 end if
118 end subroutine af_prolong_copy
119
120 !> Zeroth order prolongation
121 subroutine af_prolong_zeroth(box_p, box_c, iv, iv_to, add, limiter)
122 type(box_t), intent(in) :: box_p !< Parent box
123 type(box_t), intent(inout) :: box_c !< Child box
124 integer, intent(in) :: iv !< Variable to fill
125 integer, intent(in), optional :: iv_to !< Destination variable
126 logical, intent(in), optional :: add !< Add to old values
127 integer, intent(in), optional :: limiter !< What kind of limiter to use
128 integer :: hnc, nc, ix_offset(ndim), ivc
129 integer :: ijk, ijk_(c), ijk_(f)
130 real(dp) :: f0
131 logical :: add_to
132
133 nc = box_c%n_cell
134 hnc = ishft(box_c%n_cell, -1)
135 ix_offset = af_get_child_offset(box_c)
136 add_to = .false.; if (present(add)) add_to = add
137 ivc = iv; if (present(iv_to)) ivc = iv_to
138
139 if (.not. add_to) then
140 box_c%cc(dtimes(1:nc), ivc) = 0
141 end if
142
143#if NDIM == 1
144 do i = 1, hnc
145 i_c = i + ix_offset(1)
146 i_f = 2 * i - 1
147
148 f0 = box_p%cc(i_c, iv)
149 box_c%cc(i_f:i_f+1, ivc) = &
150 box_c%cc(i_f:i_f+1, ivc) + f0
151 end do
152#elif NDIM == 2
153 do j = 1, hnc
154 j_c = j + ix_offset(2)
155 j_f = 2 * j - 1
156 do i = 1, hnc
157 i_c = i + ix_offset(1)
158 i_f = 2 * i - 1
159
160 f0 = box_p%cc(i_c, j_c, iv)
161 box_c%cc(i_f:i_f+1, j_f:j_f+1, ivc) = &
162 box_c%cc(i_f:i_f+1, j_f:j_f+1, ivc) + f0
163 end do
164 end do
165#elif NDIM == 3
166 do k = 1, hnc
167 k_c = k + ix_offset(3)
168 k_f = 2 * k - 1
169 do j = 1, hnc
170 j_c = j + ix_offset(2)
171 j_f = 2 * j - 1
172 do i = 1, hnc
173 i_c = i + ix_offset(1)
174 i_f = 2 * i - 1
175
176 f0 = box_p%cc(i_c, j_c, k_c, iv)
177 box_c%cc(i_f:i_f+1, j_f:j_f+1, k_f:k_f+1, ivc) = &
178 box_c%cc(i_f:i_f+1, j_f:j_f+1, k_f:k_f+1, ivc) + f0
179 end do
180 end do
181 end do
182#endif
183 end subroutine af_prolong_zeroth
184
185 !> Linear prolongation to children. We use 2-1-1 interpolation (2d) and
186 !> 1-1-1-1 interpolation (3D), which do not require corner ghost cells.
187 subroutine af_prolong_linear_from(boxes, id, iv, iv_to, add)
188 type(box_t), intent(inout) :: boxes(:) !< List of all boxes
189 integer, intent(in) :: id !< Box whose children we will fill
190 integer, intent(in) :: iv !< Variable that is filled
191 integer, intent(in), optional :: iv_to !< Destination variable
192 logical, intent(in), optional :: add !< Add to old values
193 integer :: i_c, c_id
194
195 do i_c = 1, af_num_children
196 c_id = boxes(id)%children(i_c)
197 if (c_id == af_no_box) cycle
198 call af_prolong_linear(boxes(id), boxes(c_id), iv, iv_to, add)
199 end do
200 end subroutine af_prolong_linear_from
201
202 !> Prolongation to a child (from parent) using linear interpolation. We use
203 !> 2-1-1 interpolation (2D) and 1-1-1-1 interpolation (3D) which do not need
204 !> corner ghost cells.
205 subroutine af_prolong_sparse(box_p, box_c, iv, iv_to, add, limiter)
206 type(box_t), intent(in) :: box_p !< Parent box
207 type(box_t), intent(inout) :: box_c !< Child box
208 integer, intent(in) :: iv !< Variable to fill
209 integer, intent(in), optional :: iv_to !< Destination variable
210 logical, intent(in), optional :: add !< Add to old values
211 integer, intent(in), optional :: limiter !< What kind of limiter to use
212 integer :: hnc, nc, ix_offset(ndim), ivc
213 integer :: ijk, ijk_(c), ijk_(f)
214 real(dp) :: f0, flx, fhx
215 logical :: add_to
216#if NDIM > 1
217 real(dp) :: fly, fhy
218#endif
219#if NDIM > 2
220 real(dp) :: flz, fhz
221#endif
222
223 nc = box_c%n_cell
224 hnc = ishft(box_c%n_cell, -1)
225 ix_offset = af_get_child_offset(box_c)
226 add_to = .false.; if (present(add)) add_to = add
227 ivc = iv; if (present(iv_to)) ivc = iv_to
228
229 if (.not. add_to) then
230 box_c%cc(dtimes(1:nc), ivc) = 0
231 end if
232
233#if NDIM == 1
234 do i = 1, hnc
235 i_c = i + ix_offset(1)
236 i_f = 2 * i - 1
237
238 f0 = 0.75_dp * box_p%cc(i_c, iv)
239 flx = 0.25_dp * box_p%cc(i_c-1, iv)
240 fhx = 0.25_dp * box_p%cc(i_c+1, iv)
241
242 box_c%cc(i_f, ivc) = f0 + flx + box_c%cc(i_f, ivc)
243 box_c%cc(i_f+1, ivc) = f0 + fhx + box_c%cc(i_f+1, ivc)
244 end do
245#elif NDIM == 2
246 do j = 1, hnc
247 j_c = j + ix_offset(2)
248 j_f = 2 * j - 1
249 do i = 1, hnc
250 i_c = i + ix_offset(1)
251 i_f = 2 * i - 1
252
253 f0 = 0.5_dp * box_p%cc(i_c, j_c, iv)
254 flx = 0.25_dp * box_p%cc(i_c-1, j_c, iv)
255 fhx = 0.25_dp * box_p%cc(i_c+1, j_c, iv)
256 fly = 0.25_dp * box_p%cc(i_c, j_c-1, iv)
257 fhy = 0.25_dp * box_p%cc(i_c, j_c+1, iv)
258
259 box_c%cc(i_f, j_f, ivc) = f0 + flx + fly &
260 + box_c%cc(i_f, j_f, ivc)
261 box_c%cc(i_f+1, j_f, ivc) = f0 + fhx + fly &
262 + box_c%cc(i_f+1, j_f, ivc)
263 box_c%cc(i_f, j_f+1, ivc) = f0 + flx + fhy &
264 + box_c%cc(i_f, j_f+1, ivc)
265 box_c%cc(i_f+1, j_f+1, ivc) = f0 + fhx + fhy &
266 + box_c%cc(i_f+1, j_f+1, ivc)
267 end do
268 end do
269#elif NDIM == 3
270 do k = 1, hnc
271 k_c = k + ix_offset(3)
272 k_f = 2 * k - 1
273 do j = 1, hnc
274 j_c = j + ix_offset(2)
275 j_f = 2 * j - 1
276 do i = 1, hnc
277 i_c = i + ix_offset(1)
278 i_f = 2 * i - 1
279
280 f0 = 0.25_dp * box_p%cc(i_c, j_c, k_c, iv)
281 flx = 0.25_dp * box_p%cc(i_c-1, j_c, k_c, iv)
282 fhx = 0.25_dp * box_p%cc(i_c+1, j_c, k_c, iv)
283 fly = 0.25_dp * box_p%cc(i_c, j_c-1, k_c, iv)
284 fhy = 0.25_dp * box_p%cc(i_c, j_c+1, k_c, iv)
285 flz = 0.25_dp * box_p%cc(i_c, j_c, k_c-1, iv)
286 fhz = 0.25_dp * box_p%cc(i_c, j_c, k_c+1, iv)
287
288 box_c%cc(i_f, j_f, k_f, ivc) = f0 + flx + &
289 fly + flz + box_c%cc(i_f, j_f, k_f, ivc)
290 box_c%cc(i_f+1, j_f, k_f, ivc) = f0 + fhx + &
291 fly + flz + box_c%cc(i_f+1, j_f, k_f, ivc)
292 box_c%cc(i_f, j_f+1, k_f, ivc) = f0 + flx + &
293 fhy + flz + box_c%cc(i_f, j_f+1, k_f, ivc)
294 box_c%cc(i_f+1, j_f+1, k_f, ivc) = f0 + fhx + &
295 fhy + flz + box_c%cc(i_f+1, j_f+1, k_f, ivc)
296 box_c%cc(i_f, j_f, k_f+1, ivc) = f0 + flx + &
297 fly + fhz + box_c%cc(i_f, j_f, k_f+1, ivc)
298 box_c%cc(i_f+1, j_f, k_f+1, ivc) = f0 + fhx + &
299 fly + fhz + box_c%cc(i_f+1, j_f, k_f+1, ivc)
300 box_c%cc(i_f, j_f+1, k_f+1, ivc) = f0 + flx + &
301 fhy + fhz + box_c%cc(i_f, j_f+1, k_f+1, ivc)
302 box_c%cc(i_f+1, j_f+1, k_f+1, ivc) = f0 + fhx + &
303 fhy + fhz + box_c%cc(i_f+1, j_f+1, k_f+1, ivc)
304 end do
305 end do
306 end do
307#endif
308 end subroutine af_prolong_sparse
309
310 !> Conservative prolongation using the limited gradient from the coarse cells
311 subroutine af_prolong_limit(box_p, box_c, iv, iv_to, add, limiter)
312 use m_af_limiters
313 type(box_t), intent(in) :: box_p !< Parent box
314 type(box_t), intent(inout) :: box_c !< Child box
315 integer, intent(in) :: iv !< Variable to fill
316 integer, intent(in), optional :: iv_to !< Destination variable
317 logical, intent(in), optional :: add !< Add to old values
318 integer, intent(in), optional :: limiter !< What kind of limiter to use
319 integer :: hnc, nc, ix_offset(ndim), ivc
320 integer :: ijk, ijk_(c), ijk_(f)
321 real(dp) :: f(0:ndim), slopes_a(ndim), slopes_b(ndim)
322 logical :: add_to
323
324 nc = box_c%n_cell
325 hnc = ishft(box_c%n_cell, -1)
326 ix_offset = af_get_child_offset(box_c)
327 add_to = .false.; if (present(add)) add_to = add
328 ivc = iv; if (present(iv_to)) ivc = iv_to
329
330 if (.not. present(limiter)) error stop "limiter argument required"
331
332 if (.not. add_to) then
333 box_c%cc(dtimes(1:nc), ivc) = 0
334 end if
335
336 associate(cc_p => box_p%cc, cc_c => box_c%cc)
337#if NDIM == 1
338 do i = 1, hnc
339 i_c = i + ix_offset(1)
340 i_f = 2 * i - 1
341
342 f(0) = cc_p(i_c, iv)
343 slopes_a = [cc_p(i_c, iv) - cc_p(i_c-1, iv)]
344 slopes_b = [cc_p(i_c+1, iv) - cc_p(i_c, iv)]
345 f(1:1) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
346
347 cc_c(i_f, ivc) = f(0) - f(1) + cc_c(i_f, ivc)
348 cc_c(i_f+1, ivc) = f(0) + f(1) + cc_c(i_f+1, ivc)
349 end do
350#elif NDIM == 2
351 do j = 1, hnc
352 j_c = j + ix_offset(2)
353 j_f = 2 * j - 1
354 do i = 1, hnc
355 i_c = i + ix_offset(1)
356 i_f = 2 * i - 1
357
358 f(0) = cc_p(i_c, j_c, iv)
359 slopes_a = [&
360 cc_p(i_c, j_c, iv) - cc_p(i_c-1, j_c, iv), &
361 cc_p(i_c, j_c, iv) - cc_p(i_c, j_c-1, iv)]
362 slopes_b = [&
363 cc_p(i_c+1, j_c, iv) - cc_p(i_c, j_c, iv), &
364 cc_p(i_c, j_c+1, iv) - cc_p(i_c, j_c, iv)]
365 f(1:2) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
366
367 cc_c(i_f, j_f, ivc) = f(0) - f(1) - f(2) &
368 + cc_c(i_f, j_f, ivc)
369 cc_c(i_f+1, j_f, ivc) = f(0) + f(1) - f(2) &
370 + cc_c(i_f+1, j_f, ivc)
371 cc_c(i_f, j_f+1, ivc) = f(0) - f(1) + f(2) &
372 + cc_c(i_f, j_f+1, ivc)
373 cc_c(i_f+1, j_f+1, ivc) = f(0) + f(1) + f(2) &
374 + cc_c(i_f+1, j_f+1, ivc)
375 end do
376 end do
377#elif NDIM == 3
378 do k = 1, hnc
379 k_c = k + ix_offset(3)
380 k_f = 2 * k - 1
381 do j = 1, hnc
382 j_c = j + ix_offset(2)
383 j_f = 2 * j - 1
384 do i = 1, hnc
385 i_c = i + ix_offset(1)
386 i_f = 2 * i - 1
387
388 f(0) = cc_p(i_c, j_c, k_c, iv)
389 slopes_a = [&
390 cc_p(i_c, j_c, k_c, iv) - cc_p(i_c-1, j_c, k_c, iv), &
391 cc_p(i_c, j_c, k_c, iv) - cc_p(i_c, j_c-1, k_c, iv), &
392 cc_p(i_c, j_c, k_c, iv) - cc_p(i_c, j_c, k_c-1, iv)]
393 slopes_b = [&
394 cc_p(i_c+1, j_c, k_c, iv) - cc_p(i_c, j_c, k_c, iv), &
395 cc_p(i_c, j_c+1, k_c, iv) - cc_p(i_c, j_c, k_c, iv), &
396 cc_p(i_c, j_c, k_c+1, iv) - cc_p(i_c, j_c, k_c, iv)]
397 f(1:3) = 0.25_dp * af_limiter_apply(slopes_a, slopes_b, limiter)
398
399 cc_c(i_f, j_f, k_f, ivc) = f(0) - f(1) - &
400 f(2) - f(3) + cc_c(i_f, j_f, k_f, ivc)
401 cc_c(i_f+1, j_f, k_f, ivc) = f(0) + f(1) - &
402 f(2) - f(3) + cc_c(i_f+1, j_f, k_f, ivc)
403 cc_c(i_f, j_f+1, k_f, ivc) = f(0) - f(1) + &
404 f(2) - f(3) + cc_c(i_f, j_f+1, k_f, ivc)
405 cc_c(i_f+1, j_f+1, k_f, ivc) = f(0) + f(1) + &
406 f(2) - f(3) + cc_c(i_f+1, j_f+1, k_f, ivc)
407 cc_c(i_f, j_f, k_f+1, ivc) = f(0) - f(1) - &
408 f(2) + f(3) + cc_c(i_f, j_f, k_f+1, ivc)
409 cc_c(i_f+1, j_f, k_f+1, ivc) = f(0) + f(1) - &
410 f(2) + f(3) + cc_c(i_f+1, j_f, k_f+1, ivc)
411 cc_c(i_f, j_f+1, k_f+1, ivc) = f(0) - f(1) + &
412 f(2) + f(3) + cc_c(i_f, j_f+1, k_f+1, ivc)
413 cc_c(i_f+1, j_f+1, k_f+1, ivc) = f(0) + f(1) + &
414 f(2) + f(3) + cc_c(i_f+1, j_f+1, k_f+1, ivc)
415 end do
416 end do
417 end do
418#endif
419 end associate
420 end subroutine af_prolong_limit
421
422 ! Prolong with a linear (unlimited) slope in the coarse cells, which can
423 ! result in negative densities. This procedure is conservative.
424 subroutine af_prolong_linear_cons(box_p, box_c, iv, iv_to, add, limiter)
425 type(box_t), intent(in) :: box_p !< Parent box
426 type(box_t), intent(inout) :: box_c !< Child box
427 integer, intent(in) :: iv !< Variable to fill
428 integer, intent(in), optional :: iv_to !< Destination variable
429 logical, intent(in), optional :: add !< Add to old values
430 integer, intent(in), optional :: limiter !< What kind of limiter to use
431 integer :: hnc, nc, ix_offset(ndim), ivc
432 integer :: ijk, ijk_(c), ijk_(f)
433 real(dp) :: f(0:ndim)
434 logical :: add_to
435
436 nc = box_c%n_cell
437 hnc = ishft(box_c%n_cell, -1)
438 ix_offset = af_get_child_offset(box_c)
439 add_to = .false.; if (present(add)) add_to = add
440 ivc = iv; if (present(iv_to)) ivc = iv_to
441
442 if (.not. add_to) then
443 box_c%cc(dtimes(1:nc), ivc) = 0
444 end if
445
446#if NDIM == 1
447 do i = 1, hnc
448 i_c = i + ix_offset(1)
449 i_f = 2 * i - 1
450
451 f(0) = box_p%cc(i_c, iv)
452 f(1) = 0.125_dp * (box_p%cc(i_c+1, iv) - &
453 box_p%cc(i_c-1, iv))
454
455 box_c%cc(i_f, ivc) = f(0) - f(1) + box_c%cc(i_f, ivc)
456 box_c%cc(i_f+1, ivc) = f(0) + f(1) + box_c%cc(i_f+1, ivc)
457 end do
458#elif NDIM == 2
459 do j = 1, hnc
460 j_c = j + ix_offset(2)
461 j_f = 2 * j - 1
462 do i = 1, hnc
463 i_c = i + ix_offset(1)
464 i_f = 2 * i - 1
465
466 f(0) = box_p%cc(i_c, j_c, iv)
467 f(1) = 0.125_dp * (box_p%cc(i_c+1, j_c, iv) - &
468 box_p%cc(i_c-1, j_c, iv))
469 f(2) = 0.125_dp * (box_p%cc(i_c, j_c+1, iv) - &
470 box_p%cc(i_c, j_c-1, iv))
471
472 if (box_p%coord_t == af_cyl) then
473 ! Conservative prolongation for cylindrical coords
474 f(0) = f(0) - 0.25_dp * box_p%dr(1) * f(1) / &
475 af_cyl_radius_cc(box_p, i_c)
476 end if
477
478 box_c%cc(i_f, j_f, ivc) = f(0) - f(1) - f(2) &
479 + box_c%cc(i_f, j_f, ivc)
480 box_c%cc(i_f+1, j_f, ivc) = f(0) + f(1) - f(2) &
481 + box_c%cc(i_f+1, j_f, ivc)
482 box_c%cc(i_f, j_f+1, ivc) = f(0) - f(1) + f(2) &
483 + box_c%cc(i_f, j_f+1, ivc)
484 box_c%cc(i_f+1, j_f+1, ivc) = f(0) + f(1) + f(2) &
485 + box_c%cc(i_f+1, j_f+1, ivc)
486 end do
487 end do
488#elif NDIM == 3
489 do k = 1, hnc
490 k_c = k + ix_offset(3)
491 k_f = 2 * k - 1
492 do j = 1, hnc
493 j_c = j + ix_offset(2)
494 j_f = 2 * j - 1
495 do i = 1, hnc
496 i_c = i + ix_offset(1)
497 i_f = 2 * i - 1
498
499 f(0) = box_p%cc(i_c, j_c, k_c, iv)
500 f(1) = 0.125_dp * (box_p%cc(i_c+1, j_c, k_c, iv) - &
501 box_p%cc(i_c-1, j_c, k_c, iv))
502 f(2) = 0.125_dp * (box_p%cc(i_c, j_c+1, k_c, iv) - &
503 box_p%cc(i_c, j_c-1, k_c, iv))
504 f(3) = 0.125_dp * (box_p%cc(i_c, j_c, k_c+1, iv) - &
505 box_p%cc(i_c, j_c, k_c-1, iv))
506
507 box_c%cc(i_f, j_f, k_f, ivc) = f(0) - f(1) - &
508 f(2) - f(3) + box_c%cc(i_f, j_f, k_f, ivc)
509 box_c%cc(i_f+1, j_f, k_f, ivc) = f(0) + f(1) - &
510 f(2) - f(3) + box_c%cc(i_f+1, j_f, k_f, ivc)
511 box_c%cc(i_f, j_f+1, k_f, ivc) = f(0) - f(1) + &
512 f(2) - f(3) + box_c%cc(i_f, j_f+1, k_f, ivc)
513 box_c%cc(i_f+1, j_f+1, k_f, ivc) = f(0) + f(1) + &
514 f(2) - f(3) + box_c%cc(i_f+1, j_f+1, k_f, ivc)
515 box_c%cc(i_f, j_f, k_f+1, ivc) = f(0) - f(1) - &
516 f(2) + f(3) + box_c%cc(i_f, j_f, k_f+1, ivc)
517 box_c%cc(i_f+1, j_f, k_f+1, ivc) = f(0) + f(1) - &
518 f(2) + f(3) + box_c%cc(i_f+1, j_f, k_f+1, ivc)
519 box_c%cc(i_f, j_f+1, k_f+1, ivc) = f(0) - f(1) + &
520 f(2) + f(3) + box_c%cc(i_f, j_f+1, k_f+1, ivc)
521 box_c%cc(i_f+1, j_f+1, k_f+1, ivc) = f(0) + f(1) + &
522 f(2) + f(3) + box_c%cc(i_f+1, j_f+1, k_f+1, ivc)
523 end do
524 end do
525 end do
526#endif
527
528 end subroutine af_prolong_linear_cons
529
530 !> Bi/trilinear prolongation to a child (from parent)
531 subroutine af_prolong_linear(box_p, box_c, iv, iv_to, add, limiter)
532 type(box_t), intent(in) :: box_p !< Parent box
533 type(box_t), intent(inout) :: box_c !< Child box
534 integer, intent(in) :: iv !< Variable to fill
535 integer, intent(in), optional :: iv_to !< Destination variable
536 logical, intent(in), optional :: add !< Add to old values
537 integer, intent(in), optional :: limiter !< What kind of limiter to use
538 integer :: hnc, nc, ix_offset(ndim), ivc
539 integer :: ijk, ijk_(c), ijk_(f)
540 logical :: add_to
541#if NDIM == 1
542 real(dp) :: f0, flx, fhx
543 real(dp), parameter :: f1=1/4.0_dp, f3=3/4.0_dp
544#elif NDIM == 2
545 real(dp) :: f0, flx, fhx, fly, fhy
546 real(dp) :: fll, fhl, flh, fhh
547 real(dp), parameter :: f1 = 1/16.0_dp, f3=3/16.0_dp, f9=9/16.0_dp
548#elif NDIM == 3
549 real(dp) :: f000, f00l, f0l0, f0ll, fl00, fl0l, fll0
550 real(dp) :: flll, f00h, f0h0, f0hh, fh00, fh0h, fhh0
551 real(dp) :: fhhh, f0lh, f0hl, fl0h, fh0l, flh0, fhl0
552 real(dp) :: fllh, flhl, fhll, fhhl, fhlh, flhh
553 real(dp), parameter :: f1 = 1/64.0_dp, f3=3/64.0_dp, f9=9/64.0_dp
554 real(dp), parameter :: f27 = 27/64.0_dp
555#endif
556
557 nc = box_c%n_cell
558 hnc = ishft(box_c%n_cell, -1)
559 ix_offset = af_get_child_offset(box_c)
560 add_to = .false.; if (present(add)) add_to = add
561 ivc = iv; if (present(iv_to)) ivc = iv_to
562
563 if (.not. add_to) then
564 box_c%cc(dtimes(1:nc), ivc) = 0
565 end if
566
567#if NDIM == 1
568 do i = 1, hnc
569 i_c = i + ix_offset(1)
570 i_f = 2 * i - 1
571
572 f0 = f3 * box_p%cc(i_c, iv)
573 flx = f1 * box_p%cc(i_c-1, iv)
574 fhx = f1 * box_p%cc(i_c+1, iv)
575
576 box_c%cc(i_f, ivc) = f0 + flx + box_c%cc(i_f, ivc)
577 box_c%cc(i_f+1, ivc) = f0 + fhx + box_c%cc(i_f+1, ivc)
578 end do
579#elif NDIM == 2
580 do j = 1, hnc
581 j_c = j + ix_offset(2)
582 j_f = 2 * j - 1
583 do i = 1, hnc
584 i_c = i + ix_offset(1)
585 i_f = 2 * i - 1
586
587 f0 = f9 * box_p%cc(i_c, j_c, iv)
588 flx = f3 * box_p%cc(i_c-1, j_c, iv)
589 fhx = f3 * box_p%cc(i_c+1, j_c, iv)
590 fly = f3 * box_p%cc(i_c, j_c-1, iv)
591 fhy = f3 * box_p%cc(i_c, j_c+1, iv)
592 fll = f1 * box_p%cc(i_c-1, j_c-1, iv)
593 fhl = f1 * box_p%cc(i_c+1, j_c-1, iv)
594 flh = f1 * box_p%cc(i_c-1, j_c+1, iv)
595 fhh = f1 * box_p%cc(i_c+1, j_c+1, iv)
596
597 box_c%cc(i_f, j_f, ivc) = f0 + flx + fly + fll &
598 + box_c%cc(i_f, j_f, ivc)
599 box_c%cc(i_f+1, j_f, ivc) = f0 + fhx + fly + fhl &
600 + box_c%cc(i_f+1, j_f, ivc)
601 box_c%cc(i_f, j_f+1, ivc) = f0 + flx + fhy + flh &
602 + box_c%cc(i_f, j_f+1, ivc)
603 box_c%cc(i_f+1, j_f+1, ivc) = f0 + fhx + fhy + fhh &
604 + box_c%cc(i_f+1, j_f+1, ivc)
605 end do
606 end do
607#elif NDIM == 3
608 do k = 1, hnc
609 k_c = k + ix_offset(3)
610 k_f = 2 * k - 1
611 do j = 1, hnc
612 j_c = j + ix_offset(2)
613 j_f = 2 * j - 1
614 do i = 1, hnc
615 i_c = i + ix_offset(1)
616 i_f = 2 * i - 1
617
618 f000 = f27 * box_p%cc(i_c, j_c, k_c, iv)
619
620 f00l = f9 * box_p%cc(i_c, j_c, k_c-1, iv)
621 f0l0 = f9 * box_p%cc(i_c, j_c-1, k_c, iv)
622 f0ll = f3 * box_p%cc(i_c, j_c-1, k_c-1, iv)
623 fl00 = f9 * box_p%cc(i_c-1, j_c, k_c, iv)
624 fl0l = f3 * box_p%cc(i_c-1, j_c, k_c-1, iv)
625 fll0 = f3 * box_p%cc(i_c-1, j_c-1, k_c, iv)
626 flll = f1 * box_p%cc(i_c-1, j_c-1, k_c-1, iv)
627
628 f00h = f9 * box_p%cc(i_c, j_c, k_c+1, iv)
629 f0h0 = f9 * box_p%cc(i_c, j_c+1, k_c, iv)
630 f0hh = f3 * box_p%cc(i_c, j_c+1, k_c+1, iv)
631 fh00 = f9 * box_p%cc(i_c+1, j_c, k_c, iv)
632 fh0h = f3 * box_p%cc(i_c+1, j_c, k_c+1, iv)
633 fhh0 = f3 * box_p%cc(i_c+1, j_c+1, k_c, iv)
634 fhhh = f1 * box_p%cc(i_c+1, j_c+1, k_c+1, iv)
635
636 fl0h = f3 * box_p%cc(i_c-1, j_c, k_c+1, iv)
637 fh0l = f3 * box_p%cc(i_c+1, j_c, k_c-1, iv)
638 flh0 = f3 * box_p%cc(i_c-1, j_c+1, k_c, iv)
639 fhl0 = f3 * box_p%cc(i_c+1, j_c-1, k_c, iv)
640 f0lh = f3 * box_p%cc(i_c, j_c-1, k_c+1, iv)
641 f0hl = f3 * box_p%cc(i_c, j_c+1, k_c-1, iv)
642
643 fllh = f1 * box_p%cc(i_c-1, j_c-1, k_c+1, iv)
644 flhl = f1 * box_p%cc(i_c-1, j_c+1, k_c-1, iv)
645 fhll = f1 * box_p%cc(i_c+1, j_c-1, k_c-1, iv)
646
647 fhhl = f1 * box_p%cc(i_c+1, j_c+1, k_c-1, iv)
648 fhlh = f1 * box_p%cc(i_c+1, j_c-1, k_c+1, iv)
649 flhh = f1 * box_p%cc(i_c-1, j_c+1, k_c+1, iv)
650
651 box_c%cc(i_f, j_f, k_f, ivc) = f000 + fl00 + &
652 f0l0 + f00l + fll0 + fl0l + f0ll + flll + &
653 box_c%cc(i_f, j_f, k_f, ivc)
654 box_c%cc(i_f+1, j_f, k_f, ivc) = f000 + fh00 + &
655 f0l0 + f00l + fhl0 + fh0l + f0ll + fhll + &
656 box_c%cc(i_f+1, j_f, k_f, ivc)
657 box_c%cc(i_f, j_f+1, k_f, ivc) = f000 + fl00 + &
658 f0h0 + f00l + flh0 + fl0l + f0hl + flhl + &
659 box_c%cc(i_f, j_f+1, k_f, ivc)
660 box_c%cc(i_f+1, j_f+1, k_f, ivc) = f000 + fh00 + &
661 f0h0 + f00l + fhh0 + fh0l + f0hl + fhhl + &
662 box_c%cc(i_f+1, j_f+1, k_f, ivc)
663 box_c%cc(i_f, j_f, k_f+1, ivc) = f000 + fl00 + &
664 f0l0 + f00h + fll0 + fl0h + f0lh + fllh + &
665 box_c%cc(i_f, j_f, k_f+1, ivc)
666 box_c%cc(i_f+1, j_f, k_f+1, ivc) = f000 + fh00 + &
667 f0l0 + f00h + fhl0 + fh0h + f0lh + fhlh + &
668 box_c%cc(i_f+1, j_f, k_f+1, ivc)
669 box_c%cc(i_f, j_f+1, k_f+1, ivc) = f000 + fl00 + &
670 f0h0 + f00h + flh0 + fl0h + f0hh + flhh + &
671 box_c%cc(i_f, j_f+1, k_f+1, ivc)
672 box_c%cc(i_f+1, j_f+1, k_f+1, ivc) = f000 + fh00 + &
673 f0h0 + f00h + fhh0 + fh0h + f0hh + fhhh + &
674 box_c%cc(i_f+1, j_f+1, k_f+1, ivc)
675 end do
676 end do
677 end do
678#endif
679 end subroutine af_prolong_linear
680
681 ! !> Quadratic prolongation to children. We use stencils that do not require
682 ! !> corner ghost cells.
683 ! subroutine af_prolong_quadratic_from(boxes, id, iv, iv_to, add)
684 ! type(box_t), intent(inout) :: boxes(:) !< List of all boxes
685 ! integer, intent(in) :: id !< Box whose children we will fill
686 ! integer, intent(in) :: iv !< Variable that is filled
687 ! integer, intent(in), optional :: iv_to !< Destination variable
688 ! logical, intent(in), optional :: add !< Add to old values
689 ! integer :: i_c, c_id
690
691 ! do i_c = 1, af_num_children
692 ! c_id = boxes(id)%children(i_c)
693 ! if (c_id == af_no_box) cycle
694 ! call af_prolong_quadratic(boxes(id), boxes(c_id), iv, iv_to, add)
695 ! end do
696 ! end subroutine af_prolong_quadratic_from
697
698! !> Prolongation to a child (from parent) using quadratic interpolation (using
699! !> a local Taylor approximation)
700! !> \todo Mixed derivatives in 3D version
701! subroutine af_prolong_quadratic(box_p, box_c, iv, iv_to, add)
702! type(box_t), intent(in) :: box_p !< Parent box
703! type(box_t), intent(inout) :: box_c !< Child box
704! integer, intent(in) :: iv !< Variable to fill
705! integer, intent(in), optional :: iv_to !< Destination variable
706! logical, intent(in), optional :: add !< Add to old values
707! logical :: add_to
708! integer :: hnc, ix_offset(NDIM)
709! integer :: i, j, nc, ivc
710! integer :: i_c, i_f, j_c, j_f
711! real(dp) :: f0, fx, fy, fxx, fyy, f2
712! #if NDIM == 2
713! real(dp) :: fxy(2**NDIM)
714! #elif NDIM == 3
715! real(dp) :: fz, fzz
716! integer :: k, k_c, k_f
717! #endif
718
719! nc = box_c%n_cell
720! hnc = ishft(box_c%n_cell, -1)
721! ix_offset = af_get_child_offset(box_c)
722! add_to = .false.; if (present(add)) add_to = add
723! ivc = iv; if (present(iv_to)) ivc = iv_to
724
725! if (.not. add_to) then
726! #if NDIM == 2
727! box_c%cc(1:nc, 1:nc, ivc) = 0
728! #elif NDIM == 3
729! box_c%cc(1:nc, 1:nc, 1:nc, ivc) = 0
730! #endif
731! end if
732
733! #if NDIM == 2
734! do j = 1, hnc
735! j_c = j + ix_offset(2)
736! j_f = 2 * j - 1
737! do i = 1, hnc
738! i_c = i + ix_offset(1)
739! i_f = 2 * i - 1
740
741! f0 = box_p%cc(i_c, j_c, iv)
742! fx = 0.125_dp * (box_p%cc(i_c+1, j_c, iv) - &
743! box_p%cc(i_c-1, j_c, iv))
744! fy = 0.125_dp * (box_p%cc(i_c, j_c+1, iv) - &
745! box_p%cc(i_c, j_c-1, iv))
746! fxx = 0.03125_dp * (box_p%cc(i_c-1, j_c, iv) - &
747! 2 * f0 + box_p%cc(i_c+1, j_c, iv))
748! fyy = 0.03125_dp * (box_p%cc(i_c, j_c-1, iv) - &
749! 2 * f0 + box_p%cc(i_c, j_c+1, iv))
750! f2 = fxx + fyy
751
752! fxy(1) = 0.0625_dp * (box_p%cc(i_c-1, j_c-1, iv) + f0 - &
753! box_p%cc(i_c-1, j_c, iv) - box_p%cc(i_c, j_c-1, iv))
754! fxy(2) = 0.0625_dp * (box_p%cc(i_c+1, j_c-1, iv) + f0 - &
755! box_p%cc(i_c+1, j_c, iv) - box_p%cc(i_c, j_c-1, iv))
756! fxy(3) = 0.0625_dp * (box_p%cc(i_c-1, j_c+1, iv) + f0 - &
757! box_p%cc(i_c-1, j_c, iv) - box_p%cc(i_c, j_c+1, iv))
758! fxy(4) = 0.0625_dp * (box_p%cc(i_c+1, j_c+1, iv) + f0 - &
759! box_p%cc(i_c+1, j_c, iv) - box_p%cc(i_c, j_c+1, iv))
760
761! box_c%cc(i_f, j_f, ivc) = f0 - fx - fy + f2 + fxy(1) + &
762! box_c%cc(i_f, j_f, ivc)
763! box_c%cc(i_f+1, j_f, ivc) = f0 + fx - fy + f2 + fxy(2) + &
764! box_c%cc(i_f+1, j_f, ivc)
765! box_c%cc(i_f, j_f+1, ivc) = f0 - fx + fy + f2 + fxy(3) + &
766! box_c%cc(i_f, j_f+1, ivc)
767! box_c%cc(i_f+1, j_f+1, ivc) = f0 + fx + fy + f2 + fxy(4) + &
768! box_c%cc(i_f+1, j_f+1, ivc)
769! end do
770! end do
771! #elif NDIM == 3
772! do k = 1, hnc
773! k_c = k + ix_offset(3)
774! k_f = 2 * k - 1
775! do j = 1, hnc
776! j_c = j + ix_offset(2)
777! j_f = 2 * j - 1
778! do i = 1, hnc
779! i_c = i + ix_offset(1)
780! i_f = 2 * i - 1
781
782! f0 = box_p%cc(i_c, j_c, k_c, iv)
783! fx = 0.125_dp * (box_p%cc(i_c+1, j_c, k_c, iv) - &
784! box_p%cc(i_c-1, j_c, k_c, iv))
785! fy = 0.125_dp * (box_p%cc(i_c, j_c+1, k_c, iv) - &
786! box_p%cc(i_c, j_c-1, k_c, iv))
787! fz = 0.125_dp * (box_p%cc(i_c, j_c, k_c+1, iv) - &
788! box_p%cc(i_c, j_c, k_c-1, iv))
789! fxx = 0.03125_dp * (box_p%cc(i_c-1, j_c, k_c, iv) - &
790! 2 * f0 + box_p%cc(i_c+1, j_c, k_c, iv))
791! fyy = 0.03125_dp * (box_p%cc(i_c, j_c-1, k_c, iv) - &
792! 2 * f0 + box_p%cc(i_c, j_c+1, k_c, iv))
793! fzz = 0.03125_dp * (box_p%cc(i_c, j_c, k_c-1, iv) - &
794! 2 * f0 + box_p%cc(i_c, j_c, k_c+1, iv))
795! f2 = fxx + fyy + fzz
796
797! box_c%cc(i_f, j_f, k_f, ivc) = f0 - fx - fy - fz + f2 + &
798! box_c%cc(i_f, j_f, k_f, ivc)
799! box_c%cc(i_f+1, j_f, k_f, ivc) = f0 + fx - fy - fz + f2 + &
800! box_c%cc(i_f+1, j_f, k_f, ivc)
801! box_c%cc(i_f, j_f+1, k_f, ivc) = f0 - fx + fy - fz + f2 + &
802! box_c%cc(i_f, j_f+1, k_f, ivc)
803! box_c%cc(i_f+1, j_f+1, k_f, ivc) = f0 + fx + fy - fz + f2 + &
804! box_c%cc(i_f+1, j_f+1, k_f, ivc)
805! box_c%cc(i_f, j_f, k_f+1, ivc) = f0 - fx - fy + fz + f2 + &
806! box_c%cc(i_f, j_f, k_f+1, ivc)
807! box_c%cc(i_f+1, j_f, k_f+1, ivc) = f0 + fx - fy + fz + f2 + &
808! box_c%cc(i_f+1, j_f, k_f+1, ivc)
809! box_c%cc(i_f, j_f+1, k_f+1, ivc) = f0 - fx + fy + fz + f2 + &
810! box_c%cc(i_f, j_f+1, k_f+1, ivc)
811! box_c%cc(i_f+1, j_f+1, k_f+1, ivc) = f0 + fx + fy + fz + f2 + &
812! box_c%cc(i_f+1, j_f+1, k_f+1, ivc)
813! end do
814! end do
815! end do
816! #endif
817! end subroutine af_prolong_quadratic
818
819end module m_af_prolong
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::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 Mar 1 2025 05:02:04 for Afivo by doxygen 1.9.4