1 #include "cpp_macros.h"
10 integer,
parameter :: stencil_list_initial_size = 10
12 integer,
parameter,
public :: stencil_constant = 1
13 integer,
parameter,
public :: stencil_variable = 2
14 integer,
parameter,
public :: stencil_sparse = 3
17 integer,
parameter,
private :: num_shapes = 5
20 integer,
parameter,
public :: af_stencil_357 = 1
23 integer,
parameter,
public :: af_stencil_p234 = 2
26 integer,
parameter,
public :: af_stencil_p248 = 3
29 integer,
parameter,
public :: af_stencil_246 = 4
32 integer,
parameter,
public :: af_stencil_mask = 5
35 integer,
parameter,
public :: af_stencil_sizes(num_shapes) = &
36 [2*ndim+1, ndim+1, 2**ndim, 2*ndim, 1]
42 type(
box_t),
intent(in) :: box
47 public :: af_stencil_print_info
48 public :: af_stencil_index
49 public :: af_stencil_prepare_store
50 public :: af_stencil_remove
51 public :: af_stencil_store_box
52 public :: af_stencil_check_box
53 public :: af_stencil_allocate_coeff
54 public :: af_stencil_store
55 public :: af_stencil_apply
56 public :: af_stencil_apply_box
57 public :: af_stencil_gsrb_box
58 public :: af_stencil_prolong_box
59 public :: af_stencil_get_box
60 public :: af_stencil_try_constant
65 subroutine af_stencil_print_info(tree)
66 type(
af_t),
intent(in) :: tree
67 integer :: id, n_stencils_stored
68 integer :: n_boxes_with_stencils
69 integer :: n_constant_stored, n_variable_stored
70 integer :: n, n_stencils, n_sparse_stored
76 n_boxes_with_stencils = 0
78 do id = 1, tree%highest_id
79 if (.not. tree%boxes(id)%in_use) cycle
81 n_stencils = tree%boxes(id)%n_stencils
82 if (n_stencils > 0)
then
83 n_boxes_with_stencils = n_boxes_with_stencils + 1
84 n_stencils_stored = n_stencils_stored + n_stencils
87 select case (tree%boxes(id)%stencils(n)%stype)
88 case (stencil_constant)
89 n_constant_stored = n_constant_stored + 1
90 case (stencil_variable)
91 n_variable_stored = n_variable_stored + 1
93 n_sparse_stored = n_sparse_stored + 1
100 write(*,
'(A)')
' ** Information about stencils **'
101 write(*,
'(A, I0)')
' #boxes with stencils: ', n_boxes_with_stencils
102 write(*,
'(A, I0)')
' #stencils stored: ', n_stencils_stored
103 write(*,
'(A, I0)')
' #constant stencils: ', n_constant_stored
104 write(*,
'(A, I0)')
' #variable stencils: ', n_variable_stored
105 write(*,
'(A, I0)')
' #sparse stencils: ', n_sparse_stored
106 end subroutine af_stencil_print_info
109 pure integer function af_stencil_index(box, key)
110 type(
box_t),
intent(in) :: box
111 integer,
intent(in) :: key
114 af_stencil_index = af_stencil_none
115 do i = 1, box%n_stencils
116 if (box%stencils(i)%key == key)
then
121 end function af_stencil_index
124 subroutine af_stencil_store(tree, key, set_stencil)
125 type(
af_t),
intent(inout) :: tree
126 integer,
intent(in) :: key
129 integer :: lvl, i, id
132 do lvl = 1, tree%highest_lvl
134 do i = 1,
size(tree%lvls(lvl)%ids)
135 id = tree%lvls(lvl)%ids(i)
136 call af_stencil_store_box(tree%boxes(id), key, set_stencil)
141 end subroutine af_stencil_store
144 subroutine af_stencil_prepare_store(box, key, ix)
145 type(
box_t),
intent(inout) :: box
146 integer,
intent(in) :: key
147 integer,
intent(out) :: ix
150 ix = af_stencil_index(box, key)
152 if (ix == af_stencil_none)
then
154 if (.not.
allocated(box%stencils))
then
155 allocate(box%stencils(stencil_list_initial_size))
156 else if (box%n_stencils ==
size(box%stencils))
then
157 allocate(tmp(2 * box%n_stencils))
158 tmp(1:box%n_stencils) = box%stencils
159 call move_alloc(tmp, box%stencils)
162 box%n_stencils = box%n_stencils + 1
164 box%stencils(ix)%key = key
166 error stop
"Stencil with this key has already been stored"
168 end subroutine af_stencil_prepare_store
171 subroutine af_stencil_allocate_coeff(stencil, nc, use_f, n_sparse)
172 type(
stencil_t),
intent(inout) :: stencil
174 integer,
intent(in) :: nc
177 logical,
intent(in),
optional :: use_f
179 integer,
intent(in),
optional :: n_sparse
180 logical :: allocate_f
183 allocate_f = .false.;
if (
present(use_f)) allocate_f = use_f
184 if (stencil%shape < 1 .or. stencil%shape > num_shapes) &
185 error stop
"Unknown stencil shape"
187 n_coeff = af_stencil_sizes(stencil%shape)
189 select case (stencil%stype)
190 case (stencil_constant)
191 if (.not.
allocated(stencil%c))
then
192 allocate(stencil%c(n_coeff))
193 else if (
size(stencil%c) /= n_coeff)
then
194 deallocate(stencil%c)
195 allocate(stencil%c(n_coeff))
197 case (stencil_variable)
198 if (.not.
allocated(stencil%v))
then
199 allocate(stencil%v(n_coeff, dtimes(nc)))
200 else if (
size(stencil%v, 1) /= n_coeff)
then
201 deallocate(stencil%v)
202 allocate(stencil%v(n_coeff, dtimes(nc)))
205 if (allocate_f .and. .not.
allocated(stencil%f))
then
206 allocate(stencil%f(dtimes(nc)))
207 allocate(stencil%bc_correction(dtimes(nc)))
209 case (stencil_sparse)
210 if (.not.
present(n_sparse)) error stop
"n_sparse required"
212 if (.not.
allocated(stencil%sparse_v))
then
213 allocate(stencil%sparse_ix(ndim, n_sparse))
214 allocate(stencil%sparse_v(n_coeff, n_sparse))
215 else if (any(
size(stencil%sparse_v) /= [n_coeff, n_sparse]))
then
216 deallocate(stencil%sparse_v)
217 deallocate(stencil%sparse_ix)
218 allocate(stencil%sparse_ix(ndim, n_sparse))
219 allocate(stencil%sparse_v(n_coeff, n_sparse))
222 error stop
"Unknow stencil%stype"
224 end subroutine af_stencil_allocate_coeff
227 subroutine af_stencil_remove(tree, key)
228 type(
af_t),
intent(inout) :: tree
229 integer,
intent(in) :: key
230 integer :: lvl, i, id
233 do lvl = 1, tree%highest_lvl
235 do i = 1,
size(tree%lvls(lvl)%ids)
236 id = tree%lvls(lvl)%ids(i)
237 call af_stencil_remove_box(tree%boxes(id), key)
242 end subroutine af_stencil_remove
245 subroutine af_stencil_store_box(box, key, set_stencil)
246 type(
box_t),
intent(inout) :: box
247 integer,
intent(in) :: key
252 call af_stencil_prepare_store(box, key, ix)
253 call set_stencil(box, box%stencils(ix))
254 end subroutine af_stencil_store_box
257 subroutine af_stencil_remove_box(box, key)
258 type(
box_t),
intent(inout) :: box
259 integer,
intent(in) :: key
263 ix = af_stencil_index(box, key)
264 if (ix == af_stencil_none) error stop
"Stencil not present"
267 if (ix == box%n_stencils)
then
268 box%stencils(ix) = empty_stencil
270 box%stencils(ix) = box%stencils(box%n_stencils)
271 box%stencils(box%n_stencils) = empty_stencil
274 box%n_stencils = box%n_stencils - 1
275 end subroutine af_stencil_remove_box
278 subroutine af_stencil_check_box(box, key, ix)
279 type(
box_t),
intent(in) :: box
280 integer,
intent(in) :: key
281 integer,
intent(in) :: ix
283 if (key == af_stencil_none) &
284 error stop
"key == af_stencil_none"
286 if (af_stencil_index(box, key) /= ix) &
287 error stop
"Stencil with key not found at index"
289 associate(stencil => box%stencils(ix))
290 if (stencil%shape < 1 .or. stencil%shape > num_shapes) &
291 error stop
"Unknown stencil shape"
292 select case (stencil%stype)
293 case (stencil_constant)
294 if (.not.
allocated(stencil%c)) &
295 error stop
"stencil%c not allocated"
296 case (stencil_variable)
297 if (.not.
allocated(stencil%v)) &
298 error stop
"stencil%v not allocated"
299 case (stencil_sparse)
300 if (.not.
allocated(stencil%sparse_ix)) &
301 error stop
"stencil%sparse_ix not allocated"
302 if (.not.
allocated(stencil%sparse_v)) &
303 error stop
"stencil%sparse_v not allocated"
305 error stop
"Unknow stencil%stype"
308 end subroutine af_stencil_check_box
310 subroutine af_stencil_apply(tree, key, iv, i_out)
311 type(af_t),
intent(inout) :: tree
312 integer,
intent(in) :: key
313 integer,
intent(in) :: iv
314 integer,
intent(in) :: i_out
316 integer :: lvl, i, id
319 do lvl = 1, tree%highest_lvl
321 do i = 1,
size(tree%lvls(lvl)%ids)
322 id = tree%lvls(lvl)%ids(i)
323 call af_stencil_apply_box(tree%boxes(id), key, iv, i_out)
328 end subroutine af_stencil_apply
330 subroutine af_stencil_apply_box(box, key, iv, i_out)
331 type(box_t),
intent(inout) :: box
332 integer,
intent(in) :: key
333 integer,
intent(in) :: iv
334 integer,
intent(in) :: i_out
337 i = af_stencil_index(box, key)
338 if (i == af_stencil_none) error stop
"Stencil not stored"
340 select case (box%stencils(i)%shape)
341 case (af_stencil_357)
342 call stencil_apply_357(box, box%stencils(i), iv, i_out)
344 error stop
"Unknown stencil"
346 end subroutine af_stencil_apply_box
349 subroutine stencil_apply_357(box, stencil, iv, i_out)
350 type(box_t),
intent(inout) :: box
351 type(stencil_t),
intent(in) :: stencil
352 integer,
intent(in) :: iv
353 integer,
intent(in) :: i_out
354 real(dp) :: c(2*NDIM+1)
357 real(dp) :: rfac(2, box%n_cell), c_cyl(2*NDIM+1)
358 real(dp) :: cc_cyl(2*NDIM+1, box%n_cell)
361 if (iv == i_out) error stop
"Cannot have iv == i_out"
362 if (stencil%stype == stencil_sparse) error stop
"sparse not implemented"
364 associate(cc => box%cc, nc => box%n_cell)
367 if (stencil%stype == stencil_constant)
then
371 c(1) * cc(i, iv) + c(2) * cc(i-1, iv) + c(3) * cc(i+1, iv)
375 c = stencil%v(:, ijk)
377 c(1) * cc(i, iv) + c(2) * cc(i-1, iv) + c(3) * cc(i+1, iv)
381 if (stencil%cylindrical_gradient)
then
384 call af_cyl_flux_factors(box, rfac)
386 if (stencil%stype == stencil_constant)
then
391 cc_cyl(2:3, i) = rfac(1:2, i) * c(2:3)
392 cc_cyl(1, i) = c(1) - (cc_cyl(2, i) - c(2)) &
393 - (cc_cyl(3, i) - c(3))
394 cc_cyl(4:, i) = c(4:)
399 cc_cyl(1, i) * cc(i, j, iv) + &
400 cc_cyl(2, i) * cc(i-1, j, iv) + &
401 cc_cyl(3, i) * cc(i+1, j, iv) + &
402 cc_cyl(4, i) * cc(i, j-1, iv) + &
403 cc_cyl(5, i) * cc(i, j+1, iv)
408 c = stencil%v(:, ijk)
409 c_cyl(2:3) = rfac(1:2, i) * c(2:3)
410 c_cyl(1) = c(1) - (c_cyl(2) - c(2)) - (c_cyl(3) - c(3))
413 c_cyl(1) * cc(i, j, iv) + &
414 c_cyl(2) * cc(i-1, j, iv) + &
415 c_cyl(3) * cc(i+1, j, iv) + &
416 c_cyl(4) * cc(i, j-1, iv) + &
417 c_cyl(5) * cc(i, j+1, iv)
422 if (stencil%stype == stencil_constant)
then
426 c(1) * cc(i, j, iv) + &
427 c(2) * cc(i-1, j, iv) + &
428 c(3) * cc(i+1, j, iv) + &
429 c(4) * cc(i, j-1, iv) + &
430 c(5) * cc(i, j+1, iv)
434 c = stencil%v(:, ijk)
436 c(1) * cc(i, j, iv) + &
437 c(2) * cc(i-1, j, iv) + &
438 c(3) * cc(i+1, j, iv) + &
439 c(4) * cc(i, j-1, iv) + &
440 c(5) * cc(i, j+1, iv)
445 if (stencil%stype == stencil_constant)
then
448 cc(i, j, k, i_out) = &
449 c(1) * cc(i, j, k, iv) + &
450 c(2) * cc(i-1, j, k, iv) + &
451 c(3) * cc(i+1, j, k, iv) + &
452 c(4) * cc(i, j-1, k, iv) + &
453 c(5) * cc(i, j+1, k, iv) + &
454 c(6) * cc(i, j, k-1, iv) + &
455 c(7) * cc(i, j, k+1, iv)
459 c = stencil%v(:, ijk)
460 cc(i, j, k, i_out) = &
461 c(1) * cc(i, j, k, iv) + &
462 c(2) * cc(i-1, j, k, iv) + &
463 c(3) * cc(i+1, j, k, iv) + &
464 c(4) * cc(i, j-1, k, iv) + &
465 c(5) * cc(i, j+1, k, iv) + &
466 c(6) * cc(i, j, k-1, iv) + &
467 c(7) * cc(i, j, k+1, iv)
472 if (
allocated(stencil%bc_correction))
then
473 cc(dtimes(1:nc), i_out) = cc(dtimes(1:nc), i_out) - &
474 stencil%bc_correction
478 end subroutine stencil_apply_357
534 subroutine af_stencil_prolong_box(box_p, box_c, key, iv, iv_to, add)
535 type(box_t),
intent(in) :: box_p
536 type(box_t),
intent(inout) :: box_c
537 integer,
intent(in) :: key
538 integer,
intent(in) :: iv
539 integer,
intent(in) :: iv_to
540 logical,
intent(in),
optional :: add
545 if (
present(add)) add_to = add
547 i = af_stencil_index(box_c, key)
548 if (i == af_stencil_none) error stop
"Stencil not stored"
550 select case (box_c%stencils(i)%shape)
551 case (af_stencil_p234)
552 call stencil_prolong_234(box_p, box_c, box_c%stencils(i), &
554 case (af_stencil_p248)
555 call stencil_prolong_248(box_p, box_c, box_c%stencils(i), &
558 print *,
"box_c%stencils(i)%shape: ", box_c%stencils(i)%shape
559 error stop
"Unknown stencil"
561 end subroutine af_stencil_prolong_box
564 subroutine stencil_prolong_234(box_p, box_c, stencil, iv, iv_to, add)
565 type(box_t),
intent(in) :: box_p
566 type(box_t),
intent(inout) :: box_c
567 type(stencil_t),
intent(in) :: stencil
568 integer,
intent(in) :: iv
569 integer,
intent(in) :: iv_to
570 logical,
intent(in) :: add
572 integer :: nc, ix_offset(NDIM), IJK
573 integer :: IJK_(c1), IJK_(c2)
574 real(dp) :: c(NDIM+1)
577 ix_offset = af_get_child_offset(box_c)
578 if (stencil%stype == stencil_sparse) error stop
"sparse not implemented"
579 if (.not. add) box_c%cc(dtimes(1:nc), iv_to) = 0
584 if (stencil%stype == stencil_constant)
then
587 i_c1 = ix_offset(1) + ishft(i+1, -1)
588 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
589 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
590 c(1) * box_p%cc(i_c1, iv) + &
591 c(2) * box_p%cc(i_c2, iv)
595 c = stencil%v(:, ijk)
596 i_c1 = ix_offset(1) + ishft(i+1, -1)
597 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
598 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
599 c(1) * box_p%cc(i_c1, iv) + &
600 c(2) * box_p%cc(i_c2, iv)
604 if (stencil%stype == stencil_constant)
then
607 j_c1 = ix_offset(2) + ishft(j+1, -1)
608 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
610 i_c1 = ix_offset(1) + ishft(i+1, -1)
611 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
612 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
613 c(1) * box_p%cc(i_c1, j_c1, iv) + &
614 c(2) * box_p%cc(i_c2, j_c1, iv) + &
615 c(3) * box_p%cc(i_c1, j_c2, iv)
620 j_c1 = ix_offset(2) + ishft(j+1, -1)
621 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
623 i_c1 = ix_offset(1) + ishft(i+1, -1)
624 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
625 c = stencil%v(:, ijk)
626 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
627 c(1) * box_p%cc(i_c1, j_c1, iv) + &
628 c(2) * box_p%cc(i_c2, j_c1, iv) + &
629 c(3) * box_p%cc(i_c1, j_c2, iv)
634 if (stencil%stype == stencil_constant)
then
637 k_c1 = ix_offset(3) + ishft(k+1, -1)
638 k_c2 = k_c1 + 1 - 2 * iand(k, 1)
640 j_c1 = ix_offset(2) + ishft(j+1, -1)
641 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
643 i_c1 = ix_offset(1) + ishft(i+1, -1)
644 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
645 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
646 c(1) * box_p%cc(i_c1, j_c1, k_c1, iv) + &
647 c(2) * box_p%cc(i_c2, j_c1, k_c1, iv) + &
648 c(3) * box_p%cc(i_c1, j_c2, k_c1, iv) + &
649 c(4) * box_p%cc(i_c1, j_c1, k_c2, iv)
655 k_c1 = ix_offset(3) + ishft(k+1, -1)
656 k_c2 = k_c1 + 1 - 2 * iand(k, 1)
658 j_c1 = ix_offset(2) + ishft(j+1, -1)
659 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
661 i_c1 = ix_offset(1) + ishft(i+1, -1)
662 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
663 c = stencil%v(:, ijk)
664 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
665 c(1) * box_p%cc(i_c1, j_c1, k_c1, iv) + &
666 c(2) * box_p%cc(i_c2, j_c1, k_c1, iv) + &
667 c(3) * box_p%cc(i_c1, j_c2, k_c1, iv) + &
668 c(4) * box_p%cc(i_c1, j_c1, k_c2, iv)
674 end subroutine stencil_prolong_234
677 subroutine stencil_prolong_248(box_p, box_c, stencil, iv, iv_to, add)
678 type(box_t),
intent(in) :: box_p
679 type(box_t),
intent(inout) :: box_c
680 type(stencil_t),
intent(in) :: stencil
681 integer,
intent(in) :: iv
682 integer,
intent(in) :: iv_to
683 logical,
intent(in) :: add
685 integer :: nc, ix_offset(NDIM), IJK
686 integer :: IJK_(c1), IJK_(c2)
687 real(dp) :: c(2**NDIM)
690 ix_offset = af_get_child_offset(box_c)
691 if (stencil%stype == stencil_sparse) error stop
"sparse not implemented"
692 if (.not. add) box_c%cc(dtimes(1:nc), iv_to) = 0
697 if (stencil%stype == stencil_constant)
then
700 i_c1 = ix_offset(1) + ishft(i+1, -1)
701 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
702 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
703 c(1) * box_p%cc(i_c1, iv) + &
704 c(2) * box_p%cc(i_c2, iv)
708 c = stencil%v(:, ijk)
709 i_c1 = ix_offset(1) + ishft(i+1, -1)
710 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
711 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
712 c(1) * box_p%cc(i_c1, iv) + &
713 c(2) * box_p%cc(i_c2, iv)
717 if (stencil%stype == stencil_constant)
then
720 j_c1 = ix_offset(2) + ishft(j+1, -1)
721 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
723 i_c1 = ix_offset(1) + ishft(i+1, -1)
724 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
725 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
726 c(1) * box_p%cc(i_c1, j_c1, iv) + &
727 c(2) * box_p%cc(i_c2, j_c1, iv) + &
728 c(3) * box_p%cc(i_c1, j_c2, iv) + &
729 c(4) * box_p%cc(i_c2, j_c2, iv)
734 j_c1 = ix_offset(2) + ishft(j+1, -1)
735 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
737 i_c1 = ix_offset(1) + ishft(i+1, -1)
738 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
739 c = stencil%v(:, ijk)
740 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
741 c(1) * box_p%cc(i_c1, j_c1, iv) + &
742 c(2) * box_p%cc(i_c2, j_c1, iv) + &
743 c(3) * box_p%cc(i_c1, j_c2, iv) + &
744 c(4) * box_p%cc(i_c2, j_c2, iv)
749 if (stencil%stype == stencil_constant)
then
752 k_c1 = ix_offset(3) + ishft(k+1, -1)
753 k_c2 = k_c1 + 1 - 2 * iand(k, 1)
755 j_c1 = ix_offset(2) + ishft(j+1, -1)
756 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
758 i_c1 = ix_offset(1) + ishft(i+1, -1)
759 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
760 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
761 c(1) * box_p%cc(i_c1, j_c1, k_c1, iv) + &
762 c(2) * box_p%cc(i_c2, j_c1, k_c1, iv) + &
763 c(3) * box_p%cc(i_c1, j_c2, k_c1, iv) + &
764 c(4) * box_p%cc(i_c2, j_c2, k_c1, iv) + &
765 c(5) * box_p%cc(i_c1, j_c1, k_c2, iv) + &
766 c(6) * box_p%cc(i_c2, j_c1, k_c2, iv) + &
767 c(7) * box_p%cc(i_c1, j_c2, k_c2, iv) + &
768 c(8) * box_p%cc(i_c2, j_c2, k_c2, iv)
774 k_c1 = ix_offset(3) + ishft(k+1, -1)
775 k_c2 = k_c1 + 1 - 2 * iand(k, 1)
777 j_c1 = ix_offset(2) + ishft(j+1, -1)
778 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
780 i_c1 = ix_offset(1) + ishft(i+1, -1)
781 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
782 c = stencil%v(:, ijk)
783 box_c%cc(ijk, iv_to) = box_c%cc(ijk, iv_to) + &
784 c(1) * box_p%cc(i_c1, j_c1, k_c1, iv) + &
785 c(2) * box_p%cc(i_c2, j_c1, k_c1, iv) + &
786 c(3) * box_p%cc(i_c1, j_c2, k_c1, iv) + &
787 c(4) * box_p%cc(i_c2, j_c2, k_c1, iv) + &
788 c(5) * box_p%cc(i_c1, j_c1, k_c2, iv) + &
789 c(6) * box_p%cc(i_c2, j_c1, k_c2, iv) + &
790 c(7) * box_p%cc(i_c1, j_c2, k_c2, iv) + &
791 c(8) * box_p%cc(i_c2, j_c2, k_c2, iv)
797 end subroutine stencil_prolong_248
800 subroutine af_stencil_gsrb_box(box, key, redblack, iv, i_rhs)
801 type(box_t),
intent(inout) :: box
802 integer,
intent(in) :: key
803 integer,
intent(in) :: redblack
804 integer,
intent(in) :: iv
805 integer,
intent(in) :: i_rhs
808 i = af_stencil_index(box, key)
809 if (i == af_stencil_none) error stop
"Stencil not stored"
811 select case (box%stencils(i)%shape)
812 case (af_stencil_357)
813 call stencil_gsrb_357(box, box%stencils(i), redblack, iv, i_rhs)
815 error stop
"Unknown stencil"
817 end subroutine af_stencil_gsrb_box
820 subroutine stencil_gsrb_357(box, stencil, redblack, iv, i_rhs)
821 type(box_t),
intent(inout) :: box
822 type(stencil_t),
intent(in) :: stencil
823 integer,
intent(in) :: redblack
824 integer,
intent(in) :: iv
825 integer,
intent(in) :: i_rhs
827 real(dp) :: c(2*NDIM+1), inv_c1
828 integer :: IJK, i0, nc
830 real(dp) :: rfac(2, box%n_cell), c_cyl(2*NDIM+1)
831 real(dp) :: cc_cyl(2*NDIM+1, box%n_cell), inv_cc1(box%n_cell)
834 if (stencil%stype == stencil_sparse) error stop
"sparse not implemented"
837 associate(cc => box%cc, nc => box%n_cell)
838 if (
allocated(stencil%bc_correction))
then
839 cc(dtimes(1:nc), i_rhs) = cc(dtimes(1:nc), i_rhs) + &
840 stencil%bc_correction
844 i0 = 2 - iand(redblack, 1)
845 if (stencil%stype == stencil_constant)
then
850 cc(ijk, iv) = (cc(ijk, i_rhs) &
851 - c(2) * cc(i-1, iv) &
852 - c(3) * cc(i+1, iv)) * inv_c1
856 c = stencil%v(:, ijk)
857 cc(ijk, iv) = (cc(ijk, i_rhs) &
858 - c(2) * cc(i-1, iv) &
859 - c(3) * cc(i+1, iv)) / c(1)
863 if (stencil%cylindrical_gradient)
then
866 call af_cyl_flux_factors(box, rfac)
868 if (stencil%stype == stencil_constant)
then
873 cc_cyl(2:3, i) = rfac(1:2, i) * c(2:3)
874 cc_cyl(1, i) = c(1) - (cc_cyl(2, i) - c(2)) &
875 - (cc_cyl(3, i) - c(3))
876 cc_cyl(4:, i) = c(4:)
877 inv_cc1(i) = 1 / cc_cyl(1, i)
881 i0 = 2 - iand(ieor(redblack, j), 1)
883 cc(ijk, iv) = (cc(ijk, i_rhs) &
884 - cc_cyl(2, i) * cc(i-1, j, iv) &
885 - cc_cyl(3, i) * cc(i+1, j, iv) &
886 - cc_cyl(4, i) * cc(i, j-1, iv) &
887 - cc_cyl(5, i) * cc(i, j+1, iv)) * inv_cc1(i)
893 i0 = 2 - iand(ieor(redblack, j), 1)
895 c = stencil%v(:, ijk)
896 c_cyl(2:3) = rfac(1:2, i) * c(2:3)
897 c_cyl(1) = c(1) - (c_cyl(2) - c(2)) - (c_cyl(3) - c(3))
900 cc(ijk, iv) = (cc(ijk, i_rhs) &
901 - c_cyl(2) * cc(i-1, j, iv) &
902 - c_cyl(3) * cc(i+1, j, iv) &
903 - c_cyl(4) * cc(i, j-1, iv) &
904 - c_cyl(5) * cc(i, j+1, iv)) / c_cyl(1)
909 if (stencil%stype == stencil_constant)
then
914 i0 = 2 - iand(ieor(redblack, j), 1)
916 cc(ijk, iv) = (cc(ijk, i_rhs) &
917 - c(2) * cc(i-1, j, iv) &
918 - c(3) * cc(i+1, j, iv) &
919 - c(4) * cc(i, j-1, iv) &
920 - c(5) * cc(i, j+1, iv)) * inv_c1
925 i0 = 2 - iand(ieor(redblack, j), 1)
927 c = stencil%v(:, ijk)
928 cc(ijk, iv) = (cc(ijk, i_rhs) &
929 - c(2) * cc(i-1, j, iv) &
930 - c(3) * cc(i+1, j, iv) &
931 - c(4) * cc(i, j-1, iv) &
932 - c(5) * cc(i, j+1, iv)) / c(1)
938 if (stencil%stype == stencil_constant)
then
944 i0 = 2 - iand(ieor(redblack, k+j), 1)
946 cc(ijk, iv) = (cc(ijk, i_rhs) &
947 - c(2) * cc(i-1, j, k, iv) &
948 - c(3) * cc(i+1, j, k, iv) &
949 - c(4) * cc(i, j-1, k, iv) &
950 - c(5) * cc(i, j+1, k, iv) &
951 - c(6) * cc(i, j, k-1, iv) &
952 - c(7) * cc(i, j, k+1, iv)) * inv_c1
959 i0 = 2 - iand(ieor(redblack, k+j), 1)
961 c = stencil%v(:, ijk)
962 cc(ijk, iv) = (cc(ijk, i_rhs) &
963 - c(2) * cc(i-1, j, k, iv) &
964 - c(3) * cc(i+1, j, k, iv) &
965 - c(4) * cc(i, j-1, k, iv) &
966 - c(5) * cc(i, j+1, k, iv) &
967 - c(6) * cc(i, j, k-1, iv) &
968 - c(7) * cc(i, j, k+1, iv)) / c(1)
975 if (
allocated(stencil%bc_correction))
then
976 cc(dtimes(1:nc), i_rhs) = cc(dtimes(1:nc), i_rhs) - &
977 stencil%bc_correction
980 end subroutine stencil_gsrb_357
983 subroutine af_stencil_try_constant(box, ix, abs_tol, success)
984 type(box_t),
intent(inout) :: box
985 integer,
intent(in) :: ix
986 real(dp),
intent(in) :: abs_tol
988 integer :: ijk, nc, n_coeff
990 if (box%stencils(ix)%stype == stencil_sparse) &
991 error stop
"sparse not implemented"
992 if (.not.
allocated(box%stencils(ix)%v)) &
993 error stop
"No variable stencil present"
1000 if (any(abs(box%stencils(ix)%v(:, ijk) - &
1001 box%stencils(ix)%v(:, dtimes(1))) > abs_tol))
return
1004 box%stencils(ix)%stype = stencil_constant
1005 n_coeff =
size(box%stencils(ix)%v, 1)
1006 allocate(box%stencils(ix)%c(n_coeff))
1007 box%stencils(ix)%c = box%stencils(ix)%v(:, dtimes(1))
1008 deallocate(box%stencils(ix)%v)
1010 end subroutine af_stencil_try_constant
1013 subroutine af_stencil_get_box(box, key, v)
1014 type(box_t),
intent(inout) :: box
1015 integer,
intent(in) :: key
1017 real(dp),
intent(inout) :: v(:, dtimes(:))
1020 i = af_stencil_index(box, key)
1021 if (i == af_stencil_none) error stop
"Stencil not stored"
1023 select case (box%stencils(i)%shape)
1024 case (af_stencil_357)
1025 call stencil_get_357(box, box%stencils(i), v)
1027 error stop
"Unknown stencil"
1029 end subroutine af_stencil_get_box
1032 subroutine stencil_get_357(box, stencil, v)
1033 type(box_t),
intent(inout) :: box
1034 type(stencil_t),
intent(in) :: stencil
1036 real(dp),
intent(inout) :: v(:, DTIMES(:))
1038 real(dp) :: c_cyl(size(v, 1)), rfac(2, box%n_cell)
1044 if (
size(v, 1) /= af_stencil_sizes(stencil%shape)) &
1045 error stop
"Argument v has wrong size"
1047 if (stencil%stype == stencil_constant)
then
1049 v(:, ijk) = stencil%c
1056 if (stencil%cylindrical_gradient)
then
1059 call af_cyl_flux_factors(box, rfac)
1061 c_cyl(2:3) = rfac(1:2, i) * v(2:3, ijk)
1062 c_cyl(1) = v(1, ijk) - (c_cyl(2) - v(2, ijk)) &
1063 - (c_cyl(3) - v(3, ijk))
1064 c_cyl(4:) = v(4:, ijk)
1069 end subroutine stencil_get_357
Subroutine for setting a stencil on a box.
This module contains functionality for dealing with numerical stencils.
This module contains the basic types and constants that are used in the NDIM-dimensional version of A...
Type which stores all the boxes and levels, as well as some information about the number of boxes,...
The basic building block of afivo: a box with cell-centered and face centered data,...
Type for storing a numerical stencil for a box.