Afivo  0.3
m_af_prolong.f90
1 !> This module contains the routines related to prolongation: going from
2 !> coarse to fine variables.
3 module 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 
22 contains
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 
819 end 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 Jan 11 2025 05:02:04 for Afivo by doxygen 1.9.1