Afivo  0.3
m_a2_prolong.f90
1 
4 module m_a2_prolong
5  use m_a2_types
6 
7  implicit none
8  private
9 
10  public :: a2_prolong_copy_from
11  public :: a2_prolong_copy
12  public :: a2_prolong_linear_from
13  public :: a2_prolong_sparse
14  public :: a2_prolong_linear
15  ! public :: a2_prolong_quadratic_from
16  ! public :: a2_prolong_quadratic
17 
18  public :: a2_prolong_limit
19  public :: a2_prolong_linear_cons
20 
21 contains
22 
24  subroutine a2_prolong_copy_from(boxes, id, iv, iv_to, add)
25  type(box2_t), intent(inout) :: boxes(:)
26  integer, intent(in) :: id
27  integer, intent(in) :: iv
28  integer, intent(in), optional :: iv_to
29  logical, intent(in), optional :: add
30  integer :: i_c, c_id
31 
32  do i_c = 1, a2_num_children
33  c_id = boxes(id)%children(i_c)
34  if (c_id == af_no_box) cycle
35  call a2_prolong_copy(boxes(id), boxes(c_id), iv, iv_to=iv_to, add=add)
36  end do
37  end subroutine a2_prolong_copy_from
38 
40  subroutine a2_prolong_copy(box_p, box_c, iv, iv_to, low, high, add)
41  type(box2_t), intent(in) :: box_p
42  type(box2_t), intent(inout) :: box_c
43  integer, intent(in) :: iv
44  integer, intent(in), optional :: iv_to
45  integer, intent(in), optional :: low(2)
46  integer, intent(in), optional :: high(2)
47  logical, intent(in), optional :: add
48  logical :: add_to
49  integer :: nc, ix_offset(2), ivc
50  integer :: i, j, i_c1, j_c1, lo(2), hi(2)
51 
52  nc = box_c%n_cell
53  add_to = .false.; if (present(add)) add_to = add
54  ivc = iv; if (present(iv_to)) ivc = iv_to
55  lo = 1; if (present(low)) lo = low
56  hi = nc; if (present(high)) hi = high
57 
58  ! Offset of child w.r.t. parent
59  ix_offset = a2_get_child_offset(box_c)
60 
61  if (add_to) then
62  do j = lo(2), hi(2)
63  j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
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, j, ivc) = box_c%cc(i, j, ivc) + &
67  box_p%cc(i_c1, j_c1, iv)
68  end do
69  end do
70  else
71  do j = lo(2), hi(2)
72  j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
73  do i = lo(1), hi(1)
74  i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
75  box_c%cc(i, j, ivc) = box_p%cc(i_c1, j_c1, iv)
76  end do
77  end do
78  end if
79  end subroutine a2_prolong_copy
80 
83  subroutine a2_prolong_linear_from(boxes, id, iv, iv_to, add)
84  type(box2_t), intent(inout) :: boxes(:)
85  integer, intent(in) :: id
86  integer, intent(in) :: iv
87  integer, intent(in), optional :: iv_to
88  logical, intent(in), optional :: add
89  integer :: i_c, c_id
90 
91  do i_c = 1, a2_num_children
92  c_id = boxes(id)%children(i_c)
93  if (c_id == af_no_box) cycle
94  call a2_prolong_linear(boxes(id), boxes(c_id), iv, iv_to, add)
95  end do
96  end subroutine a2_prolong_linear_from
97 
101  subroutine a2_prolong_sparse(box_p, box_c, iv, iv_to, add)
102  type(box2_t), intent(in) :: box_p
103  type(box2_t), intent(inout) :: box_c
104  integer, intent(in) :: iv
105  integer, intent(in), optional :: iv_to
106  logical, intent(in), optional :: add
107  integer :: hnc, nc, ix_offset(2), ivc
108  integer :: i, j, i_c, i_f, j_c, j_f
109  real(dp) :: f0, flx, fhx, fly, fhy
110  logical :: add_to
111 
112  nc = box_c%n_cell
113  hnc = ishft(box_c%n_cell, -1)
114  ix_offset = a2_get_child_offset(box_c)
115  add_to = .false.; if (present(add)) add_to = add
116  ivc = iv; if (present(iv_to)) ivc = iv_to
117 
118  if (.not. add_to) then
119  box_c%cc(1:nc, 1:nc, ivc) = 0
120  end if
121 
122  do j = 1, hnc
123  j_c = j + ix_offset(2)
124  j_f = 2 * j - 1
125  do i = 1, hnc
126  i_c = i + ix_offset(1)
127  i_f = 2 * i - 1
128 
129  f0 = 0.5_dp * box_p%cc(i_c, j_c, iv)
130  flx = 0.25_dp * box_p%cc(i_c-1, j_c, iv)
131  fhx = 0.25_dp * box_p%cc(i_c+1, j_c, iv)
132  fly = 0.25_dp * box_p%cc(i_c, j_c-1, iv)
133  fhy = 0.25_dp * box_p%cc(i_c, j_c+1, iv)
134 
135  box_c%cc(i_f, j_f, ivc) = f0 + flx + fly &
136  + box_c%cc(i_f, j_f, ivc)
137  box_c%cc(i_f+1, j_f, ivc) = f0 + fhx + fly &
138  + box_c%cc(i_f+1, j_f, ivc)
139  box_c%cc(i_f, j_f+1, ivc) = f0 + flx + fhy &
140  + box_c%cc(i_f, j_f+1, ivc)
141  box_c%cc(i_f+1, j_f+1, ivc) = f0 + fhx + fhy &
142  + box_c%cc(i_f+1, j_f+1, ivc)
143  end do
144  end do
145  end subroutine a2_prolong_sparse
146 
147  ! Prolong with a limited slope in the coarse cells, which takes the minimum of
148  ! the left/right slopes if they have the same sign (else it is zero). This
149  ! procedure is conservative.
150  subroutine a2_prolong_limit(box_p, box_c, iv, iv_to, add)
151  type(box2_t), intent(in) :: box_p
152  type(box2_t), intent(inout) :: box_c
153  integer, intent(in) :: iv
154  integer, intent(in), optional :: iv_to
155  logical, intent(in), optional :: add
156  integer :: hnc, nc, ix_offset(2), ivc
157  integer :: i, j, i_c, i_f, j_c, j_f
158  real(dp) :: f0, fxL, fxR, fyL, fyR, fx, fy
159  logical :: add_to
160 
161  nc = box_c%n_cell
162  hnc = ishft(box_c%n_cell, -1)
163  ix_offset = a2_get_child_offset(box_c)
164  add_to = .false.; if (present(add)) add_to = add
165  ivc = iv; if (present(iv_to)) ivc = iv_to
166 
167  if (.not. add_to) then
168  box_c%cc(1:nc, 1:nc, ivc) = 0
169  end if
170 
171  do j = 1, hnc
172  j_c = j + ix_offset(2)
173  j_f = 2 * j - 1
174  do i = 1, hnc
175  i_c = i + ix_offset(1)
176  i_f = 2 * i - 1
177 
178  f0 = box_p%cc(i_c, j_c, iv)
179  fxl = f0 - box_p%cc(i_c-1, j_c, iv)
180  fxr = box_p%cc(i_c+1, j_c, iv) - f0
181  fyl = f0 - box_p%cc(i_c, j_c-1, iv)
182  fyr = box_p%cc(i_c, j_c+1, iv) - f0
183 
184  fx = 0.25_dp * limit_slope(fxl, fxr)
185  fy = 0.25_dp * limit_slope(fyl, fyr)
186 
187  box_c%cc(i_f, j_f, ivc) = f0 - fx - fy &
188  + box_c%cc(i_f, j_f, ivc)
189  box_c%cc(i_f+1, j_f, ivc) = f0 + fx - fy &
190  + box_c%cc(i_f+1, j_f, ivc)
191  box_c%cc(i_f, j_f+1, ivc) = f0 - fx + fy &
192  + box_c%cc(i_f, j_f+1, ivc)
193  box_c%cc(i_f+1, j_f+1, ivc) = f0 + fx + fy &
194  + box_c%cc(i_f+1, j_f+1, ivc)
195  end do
196  end do
197  contains
198 
199  ! Take minimum of two slopes if they have the same sign, else take zero
200  elemental function limit_slope(ll, rr) result(slope)
201  real(dp), intent(in) :: ll, rr
202  real(dp) :: slope
203 
204  if (ll * rr < 0) then
205  slope = 0.0_dp
206  else if (ll * ll < rr * rr) then
207  slope = ll
208  else
209  slope = rr
210  end if
211  end function limit_slope
212 
213  end subroutine a2_prolong_limit
214 
215  ! Prolong with a linear (unlimited) slope in the coarse cells, which can
216  ! result in negative densities. This procedure is conservative.
217  subroutine a2_prolong_linear_cons(box_p, box_c, iv, iv_to, add)
218  type(box2_t), intent(in) :: box_p
219  type(box2_t), intent(inout) :: box_c
220  integer, intent(in) :: iv
221  integer, intent(in), optional :: iv_to
222  logical, intent(in), optional :: add
223  integer :: hnc, nc, ix_offset(2), ivc
224  integer :: i, j, i_c, i_f, j_c, j_f
225  real(dp) :: f0, fx, fy
226  logical :: add_to
227 
228  nc = box_c%n_cell
229  hnc = ishft(box_c%n_cell, -1)
230  ix_offset = a2_get_child_offset(box_c)
231  add_to = .false.; if (present(add)) add_to = add
232  ivc = iv; if (present(iv_to)) ivc = iv_to
233 
234  if (.not. add_to) then
235  box_c%cc(1:nc, 1:nc, ivc) = 0
236  end if
237 
238  do j = 1, hnc
239  j_c = j + ix_offset(2)
240  j_f = 2 * j - 1
241  do i = 1, hnc
242  i_c = i + ix_offset(1)
243  i_f = 2 * i - 1
244 
245  f0 = box_p%cc(i_c, j_c, iv)
246  fx = 0.125_dp * (box_p%cc(i_c+1, j_c, iv) - &
247  box_p%cc(i_c-1, j_c, iv))
248  fy = 0.125_dp * (box_p%cc(i_c, j_c+1, iv) - &
249  box_p%cc(i_c, j_c-1, iv))
250 
251  box_c%cc(i_f, j_f, ivc) = f0 - fx - fy &
252  + box_c%cc(i_f, j_f, ivc)
253  box_c%cc(i_f+1, j_f, ivc) = f0 + fx - fy &
254  + box_c%cc(i_f+1, j_f, ivc)
255  box_c%cc(i_f, j_f+1, ivc) = f0 - fx + fy &
256  + box_c%cc(i_f, j_f+1, ivc)
257  box_c%cc(i_f+1, j_f+1, ivc) = f0 + fx + fy &
258  + box_c%cc(i_f+1, j_f+1, ivc)
259  end do
260  end do
261 
262  end subroutine a2_prolong_linear_cons
263 
265  subroutine a2_prolong_linear(box_p, box_c, iv, iv_to, add)
266  type(box2_t), intent(in) :: box_p
267  type(box2_t), intent(inout) :: box_c
268  integer, intent(in) :: iv
269  integer, intent(in), optional :: iv_to
270  logical, intent(in), optional :: add
271  integer :: hnc, nc, ix_offset(2), ivc
272  integer :: i, j, i_c, i_f, j_c, j_f
273  logical :: add_to
274  real(dp) :: f0, flx, fhx, fly, fhy
275  real(dp) :: fll, fhl, flh, fhh
276  real(dp), parameter :: f1 = 1/16.0_dp, f3=3/16.0_dp, f9=9/16.0_dp
277 
278  nc = box_c%n_cell
279  hnc = ishft(box_c%n_cell, -1)
280  ix_offset = a2_get_child_offset(box_c)
281  add_to = .false.; if (present(add)) add_to = add
282  ivc = iv; if (present(iv_to)) ivc = iv_to
283 
284  if (.not. add_to) then
285  box_c%cc(1:nc, 1:nc, ivc) = 0
286  end if
287 
288  do j = 1, hnc
289  j_c = j + ix_offset(2)
290  j_f = 2 * j - 1
291  do i = 1, hnc
292  i_c = i + ix_offset(1)
293  i_f = 2 * i - 1
294 
295  f0 = f9 * box_p%cc(i_c, j_c, iv)
296  flx = f3 * box_p%cc(i_c-1, j_c, iv)
297  fhx = f3 * box_p%cc(i_c+1, j_c, iv)
298  fly = f3 * box_p%cc(i_c, j_c-1, iv)
299  fhy = f3 * box_p%cc(i_c, j_c+1, iv)
300  fll = f1 * box_p%cc(i_c-1, j_c-1, iv)
301  fhl = f1 * box_p%cc(i_c+1, j_c-1, iv)
302  flh = f1 * box_p%cc(i_c-1, j_c+1, iv)
303  fhh = f1 * box_p%cc(i_c+1, j_c+1, iv)
304 
305  box_c%cc(i_f, j_f, ivc) = f0 + flx + fly + fll &
306  + box_c%cc(i_f, j_f, ivc)
307  box_c%cc(i_f+1, j_f, ivc) = f0 + fhx + fly + fhl &
308  + box_c%cc(i_f+1, j_f, ivc)
309  box_c%cc(i_f, j_f+1, ivc) = f0 + flx + fhy + flh &
310  + box_c%cc(i_f, j_f+1, ivc)
311  box_c%cc(i_f+1, j_f+1, ivc) = f0 + fhx + fhy + fhh &
312  + box_c%cc(i_f+1, j_f+1, ivc)
313  end do
314  end do
315  end subroutine a2_prolong_linear
316 
317  ! !> Quadratic prolongation to children. We use stencils that do not require
318  ! !> corner ghost cells.
319  ! subroutine a2_prolong_quadratic_from(boxes, id, iv, iv_to, add)
320  ! type(box2_t), intent(inout) :: boxes(:) !< List of all boxes
321  ! integer, intent(in) :: id !< Box whose children we will fill
322  ! integer, intent(in) :: iv !< Variable that is filled
323  ! integer, intent(in), optional :: iv_to !< Destination variable
324  ! logical, intent(in), optional :: add !< Add to old values
325  ! integer :: i_c, c_id
326 
327  ! do i_c = 1, a2_num_children
328  ! c_id = boxes(id)%children(i_c)
329  ! if (c_id == af_no_box) cycle
330  ! call a2_prolong_quadratic(boxes(id), boxes(c_id), iv, iv_to, add)
331  ! end do
332  ! end subroutine a2_prolong_quadratic_from
333 
334 ! !> Prolongation to a child (from parent) using quadratic interpolation (using
335 ! !> a local Taylor approximation)
336 ! !> \todo Mixed derivatives in 3D version
337 ! subroutine a2_prolong_quadratic(box_p, box_c, iv, iv_to, add)
338 ! type(box2_t), intent(in) :: box_p !< Parent box
339 ! type(box2_t), intent(inout) :: box_c !< Child box
340 ! integer, intent(in) :: iv !< Variable to fill
341 ! integer, intent(in), optional :: iv_to !< Destination variable
342 ! logical, intent(in), optional :: add !< Add to old values
343 ! logical :: add_to
344 ! integer :: hnc, ix_offset(2)
345 ! integer :: i, j, nc, ivc
346 ! integer :: i_c, i_f, j_c, j_f
347 ! real(dp) :: f0, fx, fy, fxx, fyy, f2
348 ! #if 2 == 2
349 ! real(dp) :: fxy(2**2)
350 ! #elif 2 == 3
351 ! real(dp) :: fz, fzz
352 ! integer :: k, k_c, k_f
353 ! #endif
354 
355 ! nc = box_c%n_cell
356 ! hnc = ishft(box_c%n_cell, -1)
357 ! ix_offset = a2_get_child_offset(box_c)
358 ! add_to = .false.; if (present(add)) add_to = add
359 ! ivc = iv; if (present(iv_to)) ivc = iv_to
360 
361 ! if (.not. add_to) then
362 ! #if 2 == 2
363 ! box_c%cc(1:nc, 1:nc, ivc) = 0
364 ! #elif 2 == 3
365 ! box_c%cc(1:nc, 1:nc, 1:nc, ivc) = 0
366 ! #endif
367 ! end if
368 
369 ! #if 2 == 2
370 ! do j = 1, hnc
371 ! j_c = j + ix_offset(2)
372 ! j_f = 2 * j - 1
373 ! do i = 1, hnc
374 ! i_c = i + ix_offset(1)
375 ! i_f = 2 * i - 1
376 
377 ! f0 = box_p%cc(i_c, j_c, iv)
378 ! fx = 0.125_dp * (box_p%cc(i_c+1, j_c, iv) - &
379 ! box_p%cc(i_c-1, j_c, iv))
380 ! fy = 0.125_dp * (box_p%cc(i_c, j_c+1, iv) - &
381 ! box_p%cc(i_c, j_c-1, iv))
382 ! fxx = 0.03125_dp * (box_p%cc(i_c-1, j_c, iv) - &
383 ! 2 * f0 + box_p%cc(i_c+1, j_c, iv))
384 ! fyy = 0.03125_dp * (box_p%cc(i_c, j_c-1, iv) - &
385 ! 2 * f0 + box_p%cc(i_c, j_c+1, iv))
386 ! f2 = fxx + fyy
387 
388 ! fxy(1) = 0.0625_dp * (box_p%cc(i_c-1, j_c-1, iv) + f0 - &
389 ! box_p%cc(i_c-1, j_c, iv) - box_p%cc(i_c, j_c-1, iv))
390 ! fxy(2) = 0.0625_dp * (box_p%cc(i_c+1, j_c-1, iv) + f0 - &
391 ! box_p%cc(i_c+1, j_c, iv) - box_p%cc(i_c, j_c-1, iv))
392 ! fxy(3) = 0.0625_dp * (box_p%cc(i_c-1, j_c+1, iv) + f0 - &
393 ! box_p%cc(i_c-1, j_c, iv) - box_p%cc(i_c, j_c+1, iv))
394 ! fxy(4) = 0.0625_dp * (box_p%cc(i_c+1, j_c+1, iv) + f0 - &
395 ! box_p%cc(i_c+1, j_c, iv) - box_p%cc(i_c, j_c+1, iv))
396 
397 ! box_c%cc(i_f, j_f, ivc) = f0 - fx - fy + f2 + fxy(1) + &
398 ! box_c%cc(i_f, j_f, ivc)
399 ! box_c%cc(i_f+1, j_f, ivc) = f0 + fx - fy + f2 + fxy(2) + &
400 ! box_c%cc(i_f+1, j_f, ivc)
401 ! box_c%cc(i_f, j_f+1, ivc) = f0 - fx + fy + f2 + fxy(3) + &
402 ! box_c%cc(i_f, j_f+1, ivc)
403 ! box_c%cc(i_f+1, j_f+1, ivc) = f0 + fx + fy + f2 + fxy(4) + &
404 ! box_c%cc(i_f+1, j_f+1, ivc)
405 ! end do
406 ! end do
407 ! #elif 2 == 3
408 ! do k = 1, hnc
409 ! k_c = k + ix_offset(3)
410 ! k_f = 2 * k - 1
411 ! do j = 1, hnc
412 ! j_c = j + ix_offset(2)
413 ! j_f = 2 * j - 1
414 ! do i = 1, hnc
415 ! i_c = i + ix_offset(1)
416 ! i_f = 2 * i - 1
417 
418 ! f0 = box_p%cc(i_c, j_c, k_c, iv)
419 ! fx = 0.125_dp * (box_p%cc(i_c+1, j_c, k_c, iv) - &
420 ! box_p%cc(i_c-1, j_c, k_c, iv))
421 ! fy = 0.125_dp * (box_p%cc(i_c, j_c+1, k_c, iv) - &
422 ! box_p%cc(i_c, j_c-1, k_c, iv))
423 ! fz = 0.125_dp * (box_p%cc(i_c, j_c, k_c+1, iv) - &
424 ! box_p%cc(i_c, j_c, k_c-1, iv))
425 ! fxx = 0.03125_dp * (box_p%cc(i_c-1, j_c, k_c, iv) - &
426 ! 2 * f0 + box_p%cc(i_c+1, j_c, k_c, iv))
427 ! fyy = 0.03125_dp * (box_p%cc(i_c, j_c-1, k_c, iv) - &
428 ! 2 * f0 + box_p%cc(i_c, j_c+1, k_c, iv))
429 ! fzz = 0.03125_dp * (box_p%cc(i_c, j_c, k_c-1, iv) - &
430 ! 2 * f0 + box_p%cc(i_c, j_c, k_c+1, iv))
431 ! f2 = fxx + fyy + fzz
432 
433 ! box_c%cc(i_f, j_f, k_f, ivc) = f0 - fx - fy - fz + f2 + &
434 ! box_c%cc(i_f, j_f, k_f, ivc)
435 ! box_c%cc(i_f+1, j_f, k_f, ivc) = f0 + fx - fy - fz + f2 + &
436 ! box_c%cc(i_f+1, j_f, k_f, ivc)
437 ! box_c%cc(i_f, j_f+1, k_f, ivc) = f0 - fx + fy - fz + f2 + &
438 ! box_c%cc(i_f, j_f+1, k_f, ivc)
439 ! box_c%cc(i_f+1, j_f+1, k_f, ivc) = f0 + fx + fy - fz + f2 + &
440 ! box_c%cc(i_f+1, j_f+1, k_f, ivc)
441 ! box_c%cc(i_f, j_f, k_f+1, ivc) = f0 - fx - fy + fz + f2 + &
442 ! box_c%cc(i_f, j_f, k_f+1, ivc)
443 ! box_c%cc(i_f+1, j_f, k_f+1, ivc) = f0 + fx - fy + fz + f2 + &
444 ! box_c%cc(i_f+1, j_f, k_f+1, ivc)
445 ! box_c%cc(i_f, j_f+1, k_f+1, ivc) = f0 - fx + fy + fz + f2 + &
446 ! box_c%cc(i_f, j_f+1, k_f+1, ivc)
447 ! box_c%cc(i_f+1, j_f+1, k_f+1, ivc) = f0 + fx + fy + fz + f2 + &
448 ! box_c%cc(i_f+1, j_f+1, k_f+1, ivc)
449 ! end do
450 ! end do
451 ! end do
452 ! #endif
453 ! end subroutine a2_prolong_quadratic
454 
455 end module m_a2_prolong
m_a2_types
This module contains the basic types and constants that are used in the 2-dimensional version of Afiv...
Definition: m_a2_types.f90:5
m_a2_prolong
This module contains the routines related to prolongation: going from coarse to fine variables...
Definition: m_a2_prolong.f90:4
m_a2_types::box2_t
The basic building block of afivo: a box with cell-centered and face centered data, and information about its position, neighbors, children etc.
Definition: m_a2_types.f90:71
Generated on Wed Aug 29 2018 17:01:16 for Afivo by   doxygen 1.8.12