1 #include "../src/cpp_macros.h"
19 public :: mg_fas_vcycle
22 public :: mg_lsf_dist_linear
23 public :: mg_lsf_dist_gss
26 public :: mg_set_box_tag
27 public :: mg_set_operators_lvl
28 public :: mg_set_operators_tree
29 public :: mg_update_operator_stencil
32 public :: mg_box_lpl_gradient
36 public :: mg_compute_phi_gradient
37 public :: mg_compute_field_norm
38 public :: mg_box_field_norm
43 subroutine mg_init(tree, mg)
45 type(
af_t),
intent(inout) :: tree
46 type(
mg_t),
intent(inout) :: mg
48 if (mg%i_phi < 0) stop
"mg_init: i_phi not set"
49 if (mg%i_tmp < 0) stop
"mg_init: i_tmp not set"
50 if (mg%i_rhs < 0) stop
"mg_init: i_rhs not set"
52 if (.not.
associated(mg%sides_bc)) stop
"mg_init: sides_bc not set"
53 if (.not. tree%ready) error stop
"mg_init: tree not initialized"
56 if (.not.
associated(mg%box_op)) mg%box_op => mg_auto_op
58 if (.not.
associated(mg%box_gsrb)) mg%box_gsrb => mg_auto_gsrb
59 if (.not.
associated(mg%box_corr)) mg%box_corr => mg_auto_corr
60 if (.not.
associated(mg%box_rstr)) mg%box_rstr => mg_auto_rstr
61 if (.not.
associated(mg%sides_rb)) mg%sides_rb => mg_auto_rb
64 if (mg%operator_key == af_stencil_none)
then
65 tree%n_stencil_keys_stored = tree%n_stencil_keys_stored + 1
66 mg%operator_key = tree%n_stencil_keys_stored
69 if (mg%prolongation_key == af_stencil_none)
then
70 tree%n_stencil_keys_stored = tree%n_stencil_keys_stored + 1
71 mg%prolongation_key = tree%n_stencil_keys_stored
74 if (mg%i_lsf /= -1)
then
75 error stop
"Set tree%mg_i_lsf instead of mg%i_lsf"
76 else if (iand(mg%operator_mask, mg_lsf_box) > 0)
then
77 mg%i_lsf = tree%mg_i_lsf
80 if (mg%i_eps /= -1)
then
81 error stop
"Set tree%mg_i_eps instead of mg%i_eps"
82 else if (iand(mg%operator_mask, mg_veps_box+mg_ceps_box) > 0)
then
83 mg%i_eps = tree%mg_i_eps
86 if (mg%i_lsf /= -1)
then
87 if (.not.
associated(mg%lsf_dist)) &
88 mg%lsf_dist => mg_lsf_dist_linear
89 if (
associated(mg%lsf_dist, mg_lsf_dist_gss) .and. .not. &
90 associated(mg%lsf))
then
91 error stop
"mg_lsf_dist_gss requires mg%lsf to be set"
95 call mg_set_operators_lvl(tree, mg, 1)
96 call coarse_solver_initialize(tree, mg)
98 if (mg%i_lsf /= -1)
then
99 call check_coarse_representation_lsf(tree, mg)
102 if (.not. tree%has_cc_method(mg%i_phi))
then
104 call af_set_cc_methods(tree, mg%i_phi, mg%sides_bc, mg%sides_rb)
107 mg%initialized = .true.
109 end subroutine mg_init
111 subroutine mg_destroy(mg)
112 type(
mg_t),
intent(inout) :: mg
113 if (.not. mg%initialized) error stop
"mg%initialized is false"
114 call coarse_solver_destroy(mg%csolver)
115 end subroutine mg_destroy
118 subroutine mg_use(tree, mg)
119 type(
af_t),
intent(inout) :: tree
120 type(
mg_t),
intent(in),
target :: mg
122 if (.not. mg%initialized) error stop
"mg%initialized is false"
123 tree%mg_current_operator_mask = mg%operator_mask
125 call mg_set_operators_tree(tree, mg)
126 end subroutine mg_use
129 subroutine done_with_mg(tree)
130 type(
af_t),
intent(inout) :: tree
131 tree%mg_current_operator_mask = -1
132 end subroutine done_with_mg
137 subroutine mg_fas_fmg(tree, mg, set_residual, have_guess)
139 type(
af_t),
intent(inout) :: tree
140 type(
mg_t),
intent(inout) :: mg
141 logical,
intent(in) :: set_residual
142 logical,
intent(in) :: have_guess
145 call mg_use(tree, mg)
148 do lvl = tree%highest_lvl, 2, -1
150 call set_coarse_phi_rhs(tree, lvl, mg)
153 call init_phi_rhs(tree, mg)
158 call af_boxes_copy_cc(tree%boxes, tree%lvls(lvl)%ids, &
160 call mg_fas_vcycle(tree, mg, set_residual .and. &
161 lvl == tree%highest_lvl, lvl, standalone=.false.)
163 do lvl = 2, tree%highest_lvl
165 call af_boxes_copy_cc(tree%boxes, tree%lvls(lvl)%ids, &
170 call correct_children(tree%boxes, tree%lvls(lvl-1)%parents, mg)
173 call af_gc_lvl(tree, lvl, [mg%i_phi])
175 call mg_fas_vcycle(tree, mg, set_residual .and. &
176 lvl == tree%highest_lvl, lvl, standalone=.false.)
179 call done_with_mg(tree)
180 end subroutine mg_fas_fmg
185 subroutine mg_fas_vcycle(tree, mg, set_residual, highest_lvl, standalone)
187 type(
af_t),
intent(inout) :: tree
188 type(
mg_t),
intent(in) :: mg
189 logical,
intent(in) :: set_residual
190 integer,
intent(in),
optional :: highest_lvl
191 logical,
intent(in),
optional :: standalone
192 integer :: lvl, i, id, max_lvl
194 real(dp) :: sum_phi, mean_phi
196 by_itself = .true.;
if (
present(standalone)) by_itself = standalone
199 call mg_use(tree, mg)
202 max_lvl = tree%highest_lvl
203 if (
present(highest_lvl)) max_lvl = highest_lvl
205 do lvl = max_lvl, 2, -1
207 call gsrb_boxes(tree, lvl, mg, mg_cycle_down)
211 call update_coarse(tree, lvl, mg)
214 call solve_coarse_grid(tree, mg)
220 call correct_children(tree%boxes, tree%lvls(lvl-1)%parents, mg)
223 call af_gc_lvl(tree, lvl, [mg%i_phi])
226 call gsrb_boxes(tree, lvl, mg, mg_cycle_up)
229 if (set_residual)
then
233 do i = 1,
size(tree%lvls(lvl)%ids)
234 id = tree%lvls(lvl)%ids(i)
235 call residual_box(tree%boxes(id), mg)
244 if (mg%subtract_mean)
then
245 call af_tree_sum_cc(tree, mg%i_phi, sum_phi)
246 mean_phi = sum_phi / af_total_volume(tree)
250 do i = 1,
size(tree%lvls(lvl)%ids)
251 id = tree%lvls(lvl)%ids(i)
252 tree%boxes(id)%cc(dtimes(:), mg%i_phi) = &
253 tree%boxes(id)%cc(dtimes(:), mg%i_phi) - mean_phi
261 call done_with_mg(tree)
264 end subroutine mg_fas_vcycle
266 subroutine solve_coarse_grid(tree, mg)
268 type(
af_t),
intent(inout) :: tree
269 type(
mg_t),
intent(in) :: mg
270 integer :: num_iterations, max_threads
271 integer :: n_unknowns
273 call coarse_solver_set_rhs_phi(tree, mg)
274 n_unknowns = product(tree%coarse_grid_size)
276 if (n_unknowns < mg%csolver%min_unknowns_for_openmp)
then
278 max_threads = omp_get_max_threads()
279 call omp_set_num_threads(1)
282 call coarse_solver(mg%csolver, num_iterations)
283 call coarse_solver_get_phi(tree, mg)
286 call af_gc_lvl(tree, 1, [mg%i_phi])
288 if (n_unknowns < mg%csolver%min_unknowns_for_openmp)
then
289 call omp_set_num_threads(max_threads)
291 end subroutine solve_coarse_grid
294 subroutine mg_sides_rb(boxes, id, nb, iv)
295 type(
box_t),
intent(inout) :: boxes(:)
296 integer,
intent(in) :: id
297 integer,
intent(in) :: nb
298 integer,
intent(in) :: iv
299 integer :: nc, ix, dix, ijk, di, co(ndim)
300 integer :: hnc, p_id, p_nb_id
306 real(dp) :: tmp(0:boxes(id)%n_cell/2+1)
307 real(dp) :: gc(boxes(id)%n_cell)
310 real(dp) :: tmp(0:boxes(id)%n_cell/2+1, 0:boxes(id)%n_cell/2+1)
311 real(dp) :: gc(boxes(id)%n_cell, boxes(id)%n_cell)
314 real(dp) :: grad(ndim-1)
317 nc = boxes(id)%n_cell
320 p_id = boxes(id)%parent
321 p_nb_id = boxes(p_id)%neighbors(nb)
322 co = af_get_child_offset(boxes(id))
324 associate(box => boxes(p_nb_id))
328 case (af_neighb_lowx)
330 case (af_neighb_highx)
333 case (af_neighb_lowx)
334 tmp = box%cc(nc, co(2):co(2)+hnc+1, iv)
335 case (af_neighb_highx)
336 tmp = box%cc(1, co(2):co(2)+hnc+1, iv)
337 case (af_neighb_lowy)
338 tmp = box%cc(co(1):co(1)+hnc+1, nc, iv)
339 case (af_neighb_highy)
340 tmp = box%cc(co(1):co(1)+hnc+1, 1, iv)
342 case (af_neighb_lowx)
343 tmp = box%cc(nc, co(2):co(2)+hnc+1, co(3):co(3)+hnc+1, iv)
344 case (af_neighb_highx)
345 tmp = box%cc(1, co(2):co(2)+hnc+1, co(3):co(3)+hnc+1, iv)
346 case (af_neighb_lowy)
347 tmp = box%cc(co(1):co(1)+hnc+1, nc, co(3):co(3)+hnc+1, iv)
348 case (af_neighb_highy)
349 tmp = box%cc(co(1):co(1)+hnc+1, 1, co(3):co(3)+hnc+1, iv)
350 case (af_neighb_lowz)
351 tmp = box%cc(co(1):co(1)+hnc+1, co(2):co(2)+hnc+1, nc, iv)
352 case (af_neighb_highz)
353 tmp = box%cc(co(1):co(1)+hnc+1, co(2):co(2)+hnc+1, 1, iv)
356 error stop
"mg_sides_rb: wrong argument for nb"
366 grad(1) = 0.125_dp * (tmp(i+1) - tmp(i-1))
367 gc(2*i-1) = tmp(i) - grad(1)
368 gc(2*i) = tmp(i) + grad(1)
373 grad(1) = 0.125_dp * (tmp(i+1, j) - tmp(i-1, j))
374 grad(2) = 0.125_dp * (tmp(i, j+1) - tmp(i, j-1))
375 gc(2*i-1, 2*j-1) = tmp(i, j) - grad(1) - grad(2)
376 gc(2*i, 2*j-1) = tmp(i, j) + grad(1) - grad(2)
377 gc(2*i-1, 2*j) = tmp(i, j) - grad(1) + grad(2)
378 gc(2*i, 2*j) = tmp(i, j) + grad(1) + grad(2)
383 if (af_neighb_low(nb))
then
391 select case (af_neighb_dim(nb))
396 boxes(id)%cc(i-di, iv) = 0.5_dp * gc &
397 + 0.75_dp * boxes(id)%cc(i, iv) &
398 - 0.25_dp * boxes(id)%cc(i+di, iv)
404 dj = -1 + 2 * iand(j, 1)
405 boxes(id)%cc(i-di, j, iv) = 0.5_dp * gc(j) &
406 + 0.75_dp * boxes(id)%cc(i, j, iv) &
407 - 0.25_dp * boxes(id)%cc(i+di, j, iv)
413 di = -1 + 2 * iand(i, 1)
414 boxes(id)%cc(i, j-dj, iv) = 0.5_dp * gc(i) &
415 + 0.75_dp * boxes(id)%cc(i, j, iv) &
416 - 0.25_dp * boxes(id)%cc(i, j+dj, iv)
423 dk = -1 + 2 * iand(k, 1)
425 dj = -1 + 2 * iand(j, 1)
426 boxes(id)%cc(i-di, j, k, iv) = &
427 0.5_dp * gc(j, k) + &
428 0.75_dp * boxes(id)%cc(i, j, k, iv) - &
429 0.25_dp * boxes(id)%cc(i+di, j, k, iv)
436 dk = -1 + 2 * iand(k, 1)
438 di = -1 + 2 * iand(i, 1)
439 boxes(id)%cc(i, j-dj, k, iv) = &
440 0.5_dp * gc(i, k) + &
441 0.75_dp * boxes(id)%cc(i, j, k, iv) - &
442 0.25_dp * boxes(id)%cc(i, j+dj, k, iv)
449 dj = -1 + 2 * iand(j, 1)
451 di = -1 + 2 * iand(i, 1)
452 boxes(id)%cc(i, j, k-dk, iv) = &
453 0.5_dp * gc(i, j) + &
454 0.75_dp * boxes(id)%cc(i, j, k, iv) - &
455 0.25_dp * boxes(id)%cc(i, j, k+dk, iv)
461 end subroutine mg_sides_rb
468 subroutine mg_sides_rb_extrap(boxes, id, nb, iv)
469 type(box_t),
intent(inout) :: boxes(:)
470 integer,
intent(in) :: id
471 integer,
intent(in) :: nb
472 integer,
intent(in) :: iv
473 integer :: nc, ix, dix, IJK, di
481 nc = boxes(id)%n_cell
483 if (af_neighb_low(nb))
then
491 call af_gc_prolong_copy(boxes, id, nb, iv, 0)
493 select case (af_neighb_dim(nb))
498 boxes(id)%cc(i-di, iv) = 0.5_dp * boxes(id)%cc(i-di, iv) &
499 + 0.75_dp * boxes(id)%cc(i, iv) &
500 - 0.25_dp * boxes(id)%cc(i+di, iv)
506 dj = -1 + 2 * iand(j, 1)
509 boxes(id)%cc(i-di, j, iv) = 0.5_dp * boxes(id)%cc(i-di, j, iv) + &
510 1.125_dp * boxes(id)%cc(i, j, iv) - 0.375_dp * &
511 (boxes(id)%cc(i+di, j, iv) + boxes(id)%cc(i, j+dj, iv)) &
512 + 0.125_dp * boxes(id)%cc(i+di, j+dj, iv)
528 di = -1 + 2 * iand(i, 1)
531 boxes(id)%cc(i, j-dj, iv) = 0.5_dp * boxes(id)%cc(i, j-dj, iv) + &
532 1.125_dp * boxes(id)%cc(i, j, iv) - 0.375_dp * &
533 (boxes(id)%cc(i+di, j, iv) + boxes(id)%cc(i, j+dj, iv)) &
534 + 0.125_dp * boxes(id)%cc(i+di, j+dj, iv)
546 dk = -1 + 2 * iand(k, 1)
548 dj = -1 + 2 * iand(j, 1)
562 boxes(id)%cc(i-di, j, k, iv) = &
563 0.5_dp * boxes(id)%cc(i-di, j, k, iv) + &
564 0.75_dp * boxes(id)%cc(i, j, k, iv) - &
565 0.25_dp * boxes(id)%cc(i+di, j+dj, k+dk, iv)
572 dk = -1 + 2 * iand(k, 1)
574 di = -1 + 2 * iand(i, 1)
587 boxes(id)%cc(i, j-dj, k, iv) = &
588 0.5_dp * boxes(id)%cc(i, j-dj, k, iv) + &
589 0.75_dp * boxes(id)%cc(i, j, k, iv) - &
590 0.25_dp * boxes(id)%cc(i+di, j+dj, k+dk, iv)
597 dj = -1 + 2 * iand(j, 1)
599 di = -1 + 2 * iand(i, 1)
612 boxes(id)%cc(i, j, k-dk, iv) = &
613 0.5_dp * boxes(id)%cc(i, j, k-dk, iv) + &
614 0.75_dp * boxes(id)%cc(i, j, k, iv) - &
615 0.25_dp * boxes(id)%cc(i+di, j+dj, k+dk, iv)
621 end subroutine mg_sides_rb_extrap
624 subroutine correct_children(boxes, ids, mg)
625 type(box_t),
intent(inout) :: boxes(:)
626 integer,
intent(in) :: ids(:)
627 type(mg_t),
intent(in) :: mg
628 integer :: i, id, nc, i_c, c_id
633 nc = boxes(id)%n_cell
636 boxes(id)%cc(dtimes(:), mg%i_tmp) = boxes(id)%cc(dtimes(:), mg%i_phi) - &
637 boxes(id)%cc(dtimes(:), mg%i_tmp)
639 do i_c = 1, af_num_children
640 c_id = boxes(id)%children(i_c)
641 if (c_id == af_no_box) cycle
642 call mg%box_corr(boxes(id), boxes(c_id), mg)
646 end subroutine correct_children
648 subroutine gsrb_boxes(tree, lvl, mg, type_cycle)
649 type(af_t),
intent(inout) :: tree
650 type(mg_t),
intent(in) :: mg
651 integer,
intent(in) :: lvl
652 integer,
intent(in) :: type_cycle
653 integer :: n, i, n_cycle
654 logical :: use_corners
656 select case (type_cycle)
658 n_cycle = mg%n_cycle_down
660 n_cycle = mg%n_cycle_up
662 error stop
"gsrb_boxes: invalid cycle type"
665 associate(ids => tree%lvls(lvl)%ids)
667 do n = 1, 2 * n_cycle
670 call mg%box_gsrb(tree%boxes(ids(i)), n, mg)
676 use_corners = mg%use_corners .or. &
677 (type_cycle /= mg_cycle_down .and. n == 2 * n_cycle)
681 call af_gc_box(tree, ids(i), [mg%i_phi], use_corners)
687 end subroutine gsrb_boxes
691 subroutine update_coarse(tree, lvl, mg)
692 use m_af_utils,
only: af_box_add_cc, af_box_copy_cc
693 type(af_t),
intent(inout) :: tree
694 integer,
intent(in) :: lvl
695 type(mg_t),
intent(in) :: mg
696 integer :: i, id, p_id, nc
697 real(dp),
allocatable :: tmp(DTIMES(:))
699 id = tree%lvls(lvl)%ids(1)
701 allocate(tmp(dtimes(1:nc)))
705 do i = 1,
size(tree%lvls(lvl)%ids)
706 id = tree%lvls(lvl)%ids(i)
707 p_id = tree%boxes(id)%parent
711 tmp = tree%boxes(id)%cc(dtimes(1:nc), mg%i_tmp)
712 call residual_box(tree%boxes(id), mg)
713 call mg%box_rstr(tree%boxes(id), tree%boxes(p_id), mg%i_tmp, mg)
714 call mg%box_rstr(tree%boxes(id), tree%boxes(p_id), mg%i_phi, mg)
715 tree%boxes(id)%cc(dtimes(1:nc), mg%i_tmp) = tmp
719 call af_gc_lvl(tree, lvl-1, [mg%i_phi])
725 do i = 1,
size(tree%lvls(lvl-1)%parents)
726 id = tree%lvls(lvl-1)%parents(i)
729 call mg%box_op(tree%boxes(id), mg%i_rhs, mg)
732 call af_box_add_cc(tree%boxes(id), mg%i_tmp, mg%i_rhs)
735 call af_box_copy_cc(tree%boxes(id), mg%i_phi, mg%i_tmp)
738 end subroutine update_coarse
742 subroutine set_coarse_phi_rhs(tree, lvl, mg)
744 type(af_t),
intent(inout) :: tree
745 integer,
intent(in) :: lvl
746 type(mg_t),
intent(in) :: mg
747 integer :: i, id, p_id
750 if (lvl == tree%highest_lvl)
then
751 call af_gc_lvl(tree, lvl, [mg%i_phi])
755 do i = 1,
size(tree%lvls(lvl)%ids)
756 id = tree%lvls(lvl)%ids(i)
757 p_id = tree%boxes(id)%parent
759 call residual_box(tree%boxes(id), mg)
760 call mg%box_rstr(tree%boxes(id), tree%boxes(p_id), mg%i_tmp, mg)
761 call mg%box_rstr(tree%boxes(id), tree%boxes(p_id), mg%i_phi, mg)
765 call af_gc_lvl(tree, lvl-1, [mg%i_phi])
770 do i = 1,
size(tree%lvls(lvl-1)%parents)
771 id = tree%lvls(lvl-1)%parents(i)
772 call mg%box_op(tree%boxes(id), mg%i_rhs, mg)
773 call af_box_add_cc(tree%boxes(id), mg%i_tmp, mg%i_rhs)
776 end subroutine set_coarse_phi_rhs
779 subroutine init_phi_rhs(tree, mg)
781 type(af_t),
intent(inout) :: tree
782 type(mg_t),
intent(in) :: mg
783 integer :: i, id, p_id, lvl
786 do lvl = tree%highest_lvl, 1+1, -1
788 do i = 1,
size(tree%lvls(lvl)%ids)
789 id = tree%lvls(lvl)%ids(i)
790 call af_box_clear_cc(tree%boxes(id), mg%i_phi)
792 p_id = tree%boxes(id)%parent
793 call mg%box_rstr(tree%boxes(id), tree%boxes(p_id), mg%i_rhs, mg)
799 end subroutine init_phi_rhs
801 subroutine residual_box(box, mg)
802 type(box_t),
intent(inout) :: box
803 type(mg_t),
intent(in) :: mg
806 call mg%box_op(box, mg%i_tmp, mg)
808 box%cc(dtimes(1:nc), mg%i_tmp) = box%cc(dtimes(1:nc), mg%i_rhs) &
809 - box%cc(dtimes(1:nc), mg%i_tmp)
810 end subroutine residual_box
813 subroutine mg_auto_gsrb(box, redblack_cntr, mg)
814 type(box_t),
intent(inout) :: box
815 integer,
intent(in) :: redblack_cntr
816 type(mg_t),
intent(in) :: mg
818 call af_stencil_gsrb_box(box, mg%operator_key, redblack_cntr, &
820 end subroutine mg_auto_gsrb
823 subroutine mg_store_operator_stencil(box, mg, ix)
824 type(box_t),
intent(inout) :: box
825 type(mg_t),
intent(in) :: mg
827 integer,
intent(inout) :: ix
829 if (ix == af_stencil_none)
then
830 call af_stencil_prepare_store(box, mg%operator_key, ix)
833 select case (mg%operator_type)
834 case (mg_normal_operator)
835 call mg_box_lpl_stencil(box, mg, ix)
836 case (mg_lsf_operator)
837 call mg_box_lsf_stencil(box, mg, ix)
838 case (mg_eps_operator)
839 call mg_box_lpld_stencil(box, mg, ix)
840 case (mg_auto_operator)
842 select case (iand(box%tag, mg%operator_mask))
844 call mg_box_lpl_stencil(box, mg, ix)
846 call mg_box_lsf_stencil(box, mg, ix)
847 case (mg_veps_box, mg_ceps_box)
848 call mg_box_lpld_stencil(box, mg, ix)
849 case (mg_veps_box+mg_lsf_box, mg_ceps_box+mg_lsf_box)
850 call mg_box_lpld_lsf_stencil(box, mg, ix)
852 error stop
"mg_store_operator_stencil: unknown box tag"
855 error stop
"mg_store_operator_stencil: unknown mg%operator_type"
858 call af_stencil_check_box(box, mg%operator_key, ix)
859 end subroutine mg_store_operator_stencil
862 subroutine mg_store_prolongation_stencil(tree, id, mg)
863 type(af_t),
intent(inout) :: tree
864 integer,
intent(in) :: id
865 type(mg_t),
intent(in) :: mg
868 call af_stencil_prepare_store(tree%boxes(id), mg%prolongation_key, ix)
870 p_id = tree%boxes(id)%parent
871 if (p_id <= af_no_box) error stop
"Box does not have parent"
873 select case (mg%prolongation_type)
874 case (mg_prolong_linear)
875 call mg_box_prolong_linear_stencil(tree%boxes(id), &
876 tree%boxes(p_id), mg, ix)
877 case (mg_prolong_sparse)
878 call mg_box_prolong_sparse_stencil(tree%boxes(id), &
879 tree%boxes(p_id), mg, ix)
880 case (mg_prolong_auto)
882 associate(box=>tree%boxes(id), box_p=>tree%boxes(p_id))
883 select case (iand(box%tag, mg%operator_mask))
884 case (mg_normal_box, mg_ceps_box)
885 call mg_box_prolong_linear_stencil(box, box_p, mg, ix)
886 case (mg_lsf_box, mg_ceps_box+mg_lsf_box)
887 if (mg%lsf_use_custom_prolongation)
then
888 call mg_box_prolong_lsf_stencil(box, box_p, mg, ix)
890 call mg_box_prolong_linear_stencil(box, box_p, mg, ix)
892 case (mg_veps_box, mg_veps_box + mg_lsf_box)
893 call mg_box_prolong_eps_stencil(box, box_p, mg, ix)
895 error stop
"mg_store_prolongation_stencil: unknown box tag"
899 error stop
"mg_store_prolongation_stencil: unknown mg%prolongation_type"
902 call af_stencil_check_box(tree%boxes(id), mg%prolongation_key, ix)
903 end subroutine mg_store_prolongation_stencil
906 subroutine mg_auto_op(box, i_out, mg)
907 type(box_t),
intent(inout) :: box
908 integer,
intent(in) :: i_out
909 type(mg_t),
intent(in) :: mg
911 call af_stencil_apply_box(box, mg%operator_key, mg%i_phi, i_out)
912 end subroutine mg_auto_op
915 subroutine mg_auto_rstr(box_c, box_p, iv, mg)
916 type(box_t),
intent(in) :: box_c
917 type(box_t),
intent(inout) :: box_p
918 integer,
intent(in) :: iv
919 type(mg_t),
intent(in) :: mg
922 call mg_box_rstr_lpl(box_c, box_p, iv, mg)
923 end subroutine mg_auto_rstr
926 subroutine mg_auto_rb(boxes, id, nb, iv, op_mask)
927 type(box_t),
intent(inout) :: boxes(:)
928 integer,
intent(in) :: id
929 integer,
intent(in) :: nb
930 integer,
intent(in) :: iv
931 integer,
intent(in) :: op_mask
933 select case (iand(boxes(id)%tag, op_mask))
936 call mg_sides_rb_extrap(boxes, id, nb, iv)
938 call mg_sides_rb(boxes, id, nb, iv)
940 end subroutine mg_auto_rb
943 subroutine mg_auto_corr(box_p, box_c, mg)
944 type(box_t),
intent(inout) :: box_c
945 type(box_t),
intent(in) :: box_p
946 type(mg_t),
intent(in) :: mg
948 call af_stencil_prolong_box(box_p, box_c, mg%prolongation_key, &
949 mg%i_tmp, mg%i_phi, .true.)
950 end subroutine mg_auto_corr
954 subroutine get_possible_lsf_root_mask(box, nc, dmax, mg, root_mask)
955 type(box_t),
intent(in) :: box
956 real(dp),
intent(in) :: dmax
957 integer,
intent(in) :: nc
958 type(mg_t),
intent(in) :: mg
960 logical,
intent(out) :: root_mask(DTIMES(nc))
962 real(dp) :: gradnorm, rr(NDIM)
964 if (.not.
associated(mg%lsf)) error stop
"mg%lsf not set"
968 rr = af_r_cc(box, [ijk])
969 gradnorm = norm2(numerical_gradient(mg%lsf, rr))
970 root_mask(ijk) = (abs(box%cc(ijk, mg%i_lsf)) < dmax * gradnorm * &
971 mg%lsf_gradient_safety_factor)
973 end subroutine get_possible_lsf_root_mask
977 subroutine store_lsf_distance_matrix(box, nc, mg, boundary)
978 type(box_t),
intent(inout) :: box
979 integer,
intent(in) :: nc
980 type(mg_t),
intent(in) :: mg
981 logical,
intent(out) :: boundary
982 logical :: root_mask(DTIMES(nc))
983 real(dp) :: dd(2*NDIM), a(NDIM)
984 integer :: ixs(NDIM, nc**NDIM), IJK, ix, n
986 real(dp) :: v(2*NDIM, nc**NDIM)
988 integer :: nb, dim, i_step, n_steps
989 real(dp) :: dist, gradient(NDIM), dvec(NDIM)
990 real(dp) :: x(NDIM), step_size(NDIM), min_dr
992 min_dr = minval(box%dr)
998 call get_possible_lsf_root_mask(box, nc, norm2(box%dr), mg, root_mask)
999 n_mask = count(root_mask)
1001 if (n_mask > 0)
then
1003 call af_stencil_prepare_store(box, mg_lsf_mask_key, ix)
1004 box%stencils(ix)%stype = stencil_sparse
1005 box%stencils(ix)%shape = af_stencil_mask
1006 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell, &
1011 if (root_mask(ijk))
then
1013 box%stencils(ix)%sparse_ix(:, m) = [ijk]
1022 if (root_mask(ijk))
then
1023 a = af_r_cc(box, [ijk])
1025 dd(1) = mg%lsf_dist(a, af_r_cc(box, [i-1]), mg)
1026 dd(2) = mg%lsf_dist(a, af_r_cc(box, [i+1]), mg)
1028 dd(1) = mg%lsf_dist(a, af_r_cc(box, [i-1, j]), mg)
1029 dd(2) = mg%lsf_dist(a, af_r_cc(box, [i+1, j]), mg)
1030 dd(3) = mg%lsf_dist(a, af_r_cc(box, [i, j-1]), mg)
1031 dd(4) = mg%lsf_dist(a, af_r_cc(box, [i, j+1]), mg)
1033 dd(1) = mg%lsf_dist(a, af_r_cc(box, [i-1, j, k]), mg)
1034 dd(2) = mg%lsf_dist(a, af_r_cc(box, [i+1, j, k]), mg)
1035 dd(3) = mg%lsf_dist(a, af_r_cc(box, [i, j-1, k]), mg)
1036 dd(4) = mg%lsf_dist(a, af_r_cc(box, [i, j+1, k]), mg)
1037 dd(5) = mg%lsf_dist(a, af_r_cc(box, [i, j, k-1]), mg)
1038 dd(6) = mg%lsf_dist(a, af_r_cc(box, [i, j, k+1]), mg)
1044 if (min_dr > mg%lsf_length_scale .and. all(dd >= 1))
then
1045 n_steps = ceiling(min_dr/mg%lsf_length_scale)
1046 step_size = sign(mg%lsf_length_scale, box%cc(ijk, mg%i_lsf))
1049 do i_step = 1, n_steps
1050 gradient = numerical_gradient(mg%lsf, x)
1051 gradient = gradient/max(norm2(gradient), 1e-50_dp)
1054 x = x - gradient * step_size
1055 if (mg%lsf(x) * box%cc(ijk, mg%i_lsf) <= 0)
exit
1058 dist = mg%lsf_dist(a, x, mg)
1062 dist = dist * norm2(x - a)/min_dr
1066 dim = maxloc(abs(dvec), dim=1)
1068 if (dvec(dim) > 0) nb = nb + 1
1079 if (any(dd < 1.0_dp))
then
1088 call af_stencil_prepare_store(box, mg_lsf_distance_key, ix)
1089 box%stencils(ix)%stype = stencil_sparse
1090 box%stencils(ix)%shape = af_stencil_246
1091 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell, &
1093 box%stencils(ix)%sparse_ix(:, :) = ixs(:, 1:n)
1094 box%stencils(ix)%sparse_v(:, :) = v(:, 1:n)
1097 end subroutine store_lsf_distance_matrix
1100 subroutine mg_set_box_tag(tree, id, mg)
1101 type(af_t),
intent(inout) :: tree
1102 integer,
intent(in) :: id
1103 type(mg_t),
intent(in) :: mg
1107 associate(box => tree%boxes(id))
1108 box%tag = mg_normal_box
1109 if (tree%mg_i_lsf /= -1)
then
1110 call store_lsf_distance_matrix(box, box%n_cell, mg, has_lsf)
1113 box%tag = box%tag + mg_lsf_box
1117 if (tree%mg_i_eps /= -1)
then
1118 a = minval(box%cc(dtimes(:), tree%mg_i_eps))
1119 b = maxval(box%cc(dtimes(:), tree%mg_i_eps))
1123 box%tag = box%tag + mg_veps_box
1124 else if (maxval(abs([a, b] - 1)) > 1e-8_dp)
then
1126 box%tag = box%tag + mg_ceps_box
1131 end subroutine mg_set_box_tag
1133 subroutine mg_set_operators_lvl(tree, mg, lvl)
1134 type(af_t),
intent(inout) :: tree
1135 type(mg_t),
intent(in) :: mg
1136 integer,
intent(in) :: lvl
1140 do i = 1,
size(tree%lvls(lvl)%ids)
1141 id = tree%lvls(lvl)%ids(i)
1142 associate(box => tree%boxes(id))
1144 if (box%tag == af_init_tag)
then
1147 call mg_set_box_tag(tree, id, mg)
1151 n = af_stencil_index(box, mg%operator_key)
1152 if (n == af_stencil_none)
then
1153 call mg_store_operator_stencil(box, mg, n)
1157 if (
allocated(box%stencils(n)%f))
then
1158 box%stencils(n)%bc_correction = box%stencils(n)%f * &
1159 mg_lsf_boundary_value(box, mg)
1163 n = af_stencil_index(box, mg%prolongation_key)
1164 if (n == af_stencil_none)
then
1165 call mg_store_prolongation_stencil(tree, id, mg)
1171 end subroutine mg_set_operators_lvl
1174 subroutine mg_update_operator_stencil(tree, mg)
1175 type(af_t),
intent(inout) :: tree
1176 type(mg_t),
intent(inout) :: mg
1177 integer :: lvl, i, id, n
1180 do lvl = 1, tree%highest_lvl
1182 do i = 1,
size(tree%lvls(lvl)%ids)
1183 id = tree%lvls(lvl)%ids(i)
1184 associate(box => tree%boxes(id))
1185 n = af_stencil_index(box, mg%operator_key)
1186 if (n == af_stencil_none) error stop
"Operator stencil not yet stored"
1187 call mg_store_operator_stencil(box, mg, n)
1194 call coarse_solver_update_matrix(tree, mg)
1195 end subroutine mg_update_operator_stencil
1197 subroutine mg_set_operators_tree(tree, mg)
1198 type(af_t),
intent(inout) :: tree
1199 type(mg_t),
intent(in) :: mg
1202 do lvl = 1, tree%highest_lvl
1203 call mg_set_operators_lvl(tree, mg, lvl)
1205 end subroutine mg_set_operators_tree
1208 subroutine mg_box_rstr_lpl(box_c, box_p, iv, mg)
1210 type(box_t),
intent(in) :: box_c
1211 type(box_t),
intent(inout) :: box_p
1212 integer,
intent(in) :: iv
1213 type(mg_t),
intent(in) :: mg
1215 if (iv == mg%i_phi)
then
1218 call af_restrict_box(box_c, box_p, [iv], use_geometry=.false.)
1221 call af_restrict_box(box_c, box_p, [iv], use_geometry=.true.)
1223 end subroutine mg_box_rstr_lpl
1227 subroutine mg_box_lpl_stencil(box, mg, ix)
1228 type(box_t),
intent(inout) :: box
1229 type(mg_t),
intent(in) :: mg
1230 integer,
intent(in) :: ix
1231 real(dp) :: inv_dr2(NDIM)
1234 box%stencils(ix)%shape = af_stencil_357
1235 box%stencils(ix)%stype = stencil_constant
1236 box%stencils(ix)%cylindrical_gradient = (box%coord_t == af_cyl)
1237 inv_dr2 = 1 / box%dr**2
1238 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell)
1241 box%stencils(ix)%c(2*idim:2*idim+1) = inv_dr2(idim)
1243 box%stencils(ix)%c(1) = -sum(box%stencils(ix)%c(2:)) - mg%helmholtz_lambda
1245 end subroutine mg_box_lpl_stencil
1248 subroutine mg_box_prolong_linear_stencil(box, box_p, mg, ix)
1249 type(box_t),
intent(inout) :: box
1250 type(box_t),
intent(in) :: box_p
1251 type(mg_t),
intent(in) :: mg
1252 integer,
intent(in) :: ix
1254 box%stencils(ix)%shape = af_stencil_p248
1255 box%stencils(ix)%stype = stencil_constant
1256 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell)
1259 box%stencils(ix)%c(:) = [0.75_dp, 0.25_dp]
1261 box%stencils(ix)%c(:) = [9, 3, 3, 1] / 16.0_dp
1263 box%stencils(ix)%c(:) = [27, 9, 9, 3, 9, 3, 3, 1] / 64.0_dp
1265 end subroutine mg_box_prolong_linear_stencil
1268 subroutine mg_box_prolong_sparse_stencil(box, box_p, mg, ix)
1269 type(box_t),
intent(inout) :: box
1270 type(box_t),
intent(in) :: box_p
1271 type(mg_t),
intent(in) :: mg
1272 integer,
intent(in) :: ix
1274 box%stencils(ix)%shape = af_stencil_p234
1275 box%stencils(ix)%stype = stencil_constant
1276 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell)
1279 box%stencils(ix)%c(:) = [0.75_dp, 0.25_dp]
1281 box%stencils(ix)%c(:) = [0.5_dp, 0.25_dp, 0.25_dp]
1283 box%stencils(ix)%c(:) = [0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp]
1285 end subroutine mg_box_prolong_sparse_stencil
1289 subroutine mg_box_prolong_eps_stencil(box, box_p, mg, ix)
1290 type(box_t),
intent(inout) :: box
1291 type(box_t),
intent(in) :: box_p
1292 type(mg_t),
intent(in) :: mg
1293 integer,
intent(in) :: ix
1294 real(dp) :: a0, a(NDIM)
1295 integer :: i_eps, nc
1297 integer :: IJK, IJK_(c1)
1298 integer :: IJK_(c2), ix_offset(NDIM)
1300 real(dp),
parameter :: third = 1/3.0_dp
1304 box%stencils(ix)%shape = af_stencil_p234
1305 box%stencils(ix)%stype = stencil_variable
1306 ix_offset = af_get_child_offset(box)
1307 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell)
1311 associate(v => box%stencils(ix)%v)
1316 i_c1 = ix_offset(1) + ishft(i+1, -1)
1317 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
1319 a0 = box_p%cc(i_c1, i_eps)
1320 a(1) = box_p%cc(i_c2, i_eps)
1323 v(:, ijk) = [a0, a(1)] / (a0 + a(1))
1327 j_c1 = ix_offset(2) + ishft(j+1, -1)
1328 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
1330 i_c1 = ix_offset(1) + ishft(i+1, -1)
1331 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
1333 a0 = box_p%cc(i_c1, j_c1, i_eps)
1334 a(1) = box_p%cc(i_c2, j_c1, i_eps)
1335 a(2) = box_p%cc(i_c1, j_c2, i_eps)
1338 v(1, ijk) = 0.5_dp * sum(a0 / (a0 + a(:)))
1339 v(2:, ijk) = 0.5_dp * a(:) / (a0 + a(:))
1344 k_c1 = ix_offset(3) + ishft(k+1, -1)
1345 k_c2 = k_c1 + 1 - 2 * iand(k, 1)
1347 j_c1 = ix_offset(2) + ishft(j+1, -1)
1348 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
1350 i_c1 = ix_offset(1) + ishft(i+1, -1)
1351 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
1353 a0 = box_p%cc(i_c1, j_c1, k_c1, i_eps)
1354 a(1) = box_p%cc(i_c2, j_c1, k_c1, i_eps)
1355 a(2) = box_p%cc(i_c1, j_c2, k_c1, i_eps)
1356 a(3) = box_p%cc(i_c1, j_c1, k_c2, i_eps)
1359 v(1, ijk) = third * sum((a0 - 0.5_dp * a(:))/(a0 + a(:)))
1360 v(2:, ijk) = 0.5_dp * a(:) / (a0 + a(:))
1367 call af_stencil_try_constant(box, ix, epsilon(1.0_dp), success)
1369 end subroutine mg_box_prolong_eps_stencil
1373 subroutine mg_box_prolong_lsf_stencil(box, box_p, mg, ix)
1374 type(box_t),
intent(inout) :: box
1375 type(box_t),
intent(in) :: box_p
1376 type(mg_t),
intent(in) :: mg
1377 integer,
intent(in) :: ix
1378 real(dp) :: dd(NDIM+1), a(NDIM)
1379 integer :: i_lsf, nc, ix_mask
1381 integer :: IJK, IJK_(c1), n, n_mask
1382 integer :: IJK_(c2), ix_offset(NDIM)
1385 box%stencils(ix)%shape = af_stencil_p234
1386 box%stencils(ix)%stype = stencil_variable
1387 ix_offset = af_get_child_offset(box)
1389 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell)
1391 ix_mask = af_stencil_index(box, mg_lsf_mask_key)
1392 if (ix_mask == af_stencil_none) error stop
"No LSF root mask stored"
1393 n_mask =
size(box%stencils(ix_mask)%sparse_ix, 2)
1395 associate(v => box%stencils(ix)%v)
1403 i = box%stencils(ix_mask)%sparse_ix(1, n)
1404 i_c1 = ix_offset(1) + ishft(i+1, -1)
1405 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
1407 a = af_r_cc(box, [ijk])
1408 dd(1) = mg%lsf_dist(a, af_r_cc(box_p, [i_c1]), mg)
1409 dd(2) = mg%lsf_dist(a, af_r_cc(box_p, [i_c2]), mg)
1411 v(:, ijk) = [3 * dd(2), dd(1)]
1412 v(:, ijk) = v(:, ijk) / sum(v(:, ijk))
1413 where (dd < 1) v(:, ijk) = 0
1417 v(2, :, :) = 0.25_dp
1418 v(3, :, :) = 0.25_dp
1421 i = box%stencils(ix_mask)%sparse_ix(1, n)
1422 j = box%stencils(ix_mask)%sparse_ix(2, n)
1424 j_c1 = ix_offset(2) + ishft(j+1, -1)
1425 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
1426 i_c1 = ix_offset(1) + ishft(i+1, -1)
1427 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
1429 a = af_r_cc(box, [ijk])
1430 dd(1) = mg%lsf_dist(a, af_r_cc(box_p, [i_c1, j_c1]), mg)
1431 dd(2) = mg%lsf_dist(a, af_r_cc(box_p, [i_c2, j_c1]), mg)
1432 dd(3) = mg%lsf_dist(a, af_r_cc(box_p, [i_c1, j_c2]), mg)
1434 v(:, ijk) = [2 * dd(2) * dd(3), dd(1) * dd(3), dd(1) * dd(2)]
1435 v(:, ijk) = v(:, ijk) / sum(v(:, ijk))
1436 where (dd < 1) v(:, ijk) = 0
1439 v(:, :, :, :) = 0.25_dp
1442 i = box%stencils(ix_mask)%sparse_ix(1, n)
1443 j = box%stencils(ix_mask)%sparse_ix(2, n)
1444 k = box%stencils(ix_mask)%sparse_ix(3, n)
1446 k_c1 = ix_offset(3) + ishft(k+1, -1)
1447 k_c2 = k_c1 + 1 - 2 * iand(k, 1)
1448 j_c1 = ix_offset(2) + ishft(j+1, -1)
1449 j_c2 = j_c1 + 1 - 2 * iand(j, 1)
1450 i_c1 = ix_offset(1) + ishft(i+1, -1)
1451 i_c2 = i_c1 + 1 - 2 * iand(i, 1)
1453 a = af_r_cc(box, [ijk])
1454 dd(1) = mg%lsf_dist(a, af_r_cc(box_p, [i_c1, j_c1, k_c1]), mg)
1455 dd(2) = mg%lsf_dist(a, af_r_cc(box_p, [i_c2, j_c1, k_c1]), mg)
1456 dd(3) = mg%lsf_dist(a, af_r_cc(box_p, [i_c1, j_c2, k_c1]), mg)
1457 dd(4) = mg%lsf_dist(a, af_r_cc(box_p, [i_c1, j_c1, k_c2]), mg)
1459 v(:, ijk) = [dd(2) * dd(3) * dd(4), &
1460 dd(1) * dd(3) * dd(4), dd(1) * dd(2) * dd(4), &
1461 dd(1) * dd(2) * dd(3)]
1462 v(:, ijk) = v(:, ijk) / sum(v(:, ijk))
1463 where (dd < 1) v(:, ijk) = 0
1468 call af_stencil_try_constant(box, ix, epsilon(1.0_dp), success)
1470 end subroutine mg_box_prolong_lsf_stencil
1474 subroutine mg_box_lpld_stencil(box, mg, ix)
1475 type(box_t),
intent(inout) :: box
1476 type(mg_t),
intent(in) :: mg
1477 integer,
intent(in) :: ix
1479 real(dp) :: idr2(2*NDIM), a0, a(2*NDIM)
1483 idr2(1:2*ndim:2) = 1/box%dr**2
1484 idr2(2:2*ndim:2) = idr2(1:2*ndim:2)
1486 box%stencils(ix)%shape = af_stencil_357
1487 box%stencils(ix)%stype = stencil_variable
1488 box%stencils(ix)%cylindrical_gradient = (box%coord_t == af_cyl)
1489 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell)
1491 associate(cc => box%cc, n => mg%i_phi, i_eps => mg%i_eps)
1494 a0 = box%cc(i, i_eps)
1495 a(1:2) = box%cc(i-1:i+1:2, i_eps)
1497 a0 = box%cc(i, j, i_eps)
1498 a(1:2) = box%cc(i-1:i+1:2, j, i_eps)
1499 a(3:4) = box%cc(i, j-1:j+1:2, i_eps)
1501 a0 = box%cc(i, j, k, i_eps)
1502 a(1:2) = box%cc(i-1:i+1:2, j, k, i_eps)
1503 a(3:4) = box%cc(i, j-1:j+1:2, k, i_eps)
1504 a(5:6) = box%cc(i, j, k-1:k+1:2, i_eps)
1506 box%stencils(ix)%v(2:, ijk) = idr2 * 2 * a0*a(:)/(a0 + a(:))
1507 box%stencils(ix)%v(1, ijk) = -sum(box%stencils(ix)%v(2:, ijk))
1511 call af_stencil_try_constant(box, ix, epsilon(1.0_dp), success)
1513 end subroutine mg_box_lpld_stencil
1516 subroutine mg_box_lpld_lsf_stencil(box, mg, ix)
1517 type(box_t),
intent(inout) :: box
1518 type(mg_t),
intent(in) :: mg
1519 integer,
intent(in) :: ix
1521 real(dp) :: idr2(2*NDIM), a0, a(2*NDIM), dr2(NDIM)
1522 integer :: n, m, idim, s_ix(NDIM), ix_dist
1523 real(dp) :: dd(2*NDIM)
1524 real(dp),
allocatable :: all_distances(:, DTIMES(:))
1528 idr2(1:2*ndim:2) = 1/box%dr**2
1529 idr2(2:2*ndim:2) = idr2(1:2*ndim:2)
1532 box%stencils(ix)%shape = af_stencil_357
1533 box%stencils(ix)%stype = stencil_variable
1537 box%stencils(ix)%cylindrical_gradient = (box%coord_t == af_cyl)
1538 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell, use_f=.true.)
1539 allocate(all_distances(2*ndim, dtimes(nc)))
1541 box%stencils(ix)%f = 0.0_dp
1543 all_distances = 1.0_dp
1545 ix_dist = af_stencil_index(box, mg_lsf_distance_key)
1546 if (ix_dist == af_stencil_none) error stop
"No distances stored"
1549 do n = 1,
size(box%stencils(ix_dist)%sparse_ix, 2)
1550 s_ix = box%stencils(ix_dist)%sparse_ix(:, n)
1551 all_distances(:, dindex(s_ix)) = &
1552 box%stencils(ix_dist)%sparse_v(:, n)
1555 associate(cc => box%cc, n => mg%i_phi, i_eps => mg%i_eps)
1557 dd = all_distances(:, ijk)
1560 a0 = box%cc(i, i_eps)
1561 a(1:2) = box%cc(i-1:i+1:2, i_eps)
1563 a0 = box%cc(i, j, i_eps)
1564 a(1:2) = box%cc(i-1:i+1:2, j, i_eps)
1565 a(3:4) = box%cc(i, j-1:j+1:2, i_eps)
1567 a0 = box%cc(i, j, k, i_eps)
1568 a(1:2) = box%cc(i-1:i+1:2, j, k, i_eps)
1569 a(3:4) = box%cc(i, j-1:j+1:2, k, i_eps)
1570 a(5:6) = box%cc(i, j, k-1:k+1:2, i_eps)
1574 box%stencils(ix)%v(1+2*idim-1, ijk) = 1 / &
1575 (0.5_dp * dr2(idim) * (dd(2*idim-1) + dd(2*idim)) * &
1577 box%stencils(ix)%v(1+2*idim, ijk) = 1 / &
1578 (0.5_dp * dr2(idim) * (dd(2*idim-1) + dd(2*idim)) * &
1585 box%stencils(ix)%v(2:, ijk) = box%stencils(ix)%v(2:, ijk) * &
1586 2 * a0*a(:)/(a0 + a(:))
1588 box%stencils(ix)%v(1, ijk) = -sum(box%stencils(ix)%v(2:, ijk))
1592 if (dd(m) < 1.0_dp)
then
1593 box%stencils(ix)%f(ijk) = box%stencils(ix)%f(ijk) - &
1594 box%stencils(ix)%v(m+1, ijk)
1595 box%stencils(ix)%v(m+1, ijk) = 0.0_dp
1601 call af_stencil_try_constant(box, ix, epsilon(1.0_dp), success)
1603 end subroutine mg_box_lpld_lsf_stencil
1607 function mg_lsf_dist_linear(a, b, mg)
result(dist)
1608 real(dp),
intent(in) :: a(ndim)
1609 real(dp),
intent(in) :: b(ndim)
1610 type(mg_t),
intent(in) :: mg
1611 real(dp) :: dist, lsf_a, lsf_b
1616 if (lsf_a * lsf_b < 0)
then
1618 dist = lsf_a / (lsf_a - lsf_b)
1619 dist = max(dist, mg%lsf_min_rel_distance)
1623 end function mg_lsf_dist_linear
1629 function mg_lsf_dist_gss(a, b, mg)
result(dist)
1630 real(dp),
intent(in) :: a(ndim)
1631 real(dp),
intent(in) :: b(ndim)
1632 type(mg_t),
intent(in) :: mg
1633 real(dp) :: bracket(ndim, 2), b_new(ndim)
1634 real(dp) :: dist, r_root(ndim), lsf_a, lsf_b
1635 integer,
parameter :: max_iter = 100
1641 if (lsf_a * lsf_b <= 0)
then
1642 r_root = bisection(mg%lsf, a, b, mg%lsf_tol, max_iter)
1645 bracket = gss(mg%lsf, a, b, minimization=(lsf_a >= 0), &
1646 tol=mg%lsf_tol, find_bracket=.true.)
1649 if (mg%lsf(bracket(:, 1)) * lsf_a <= 0)
then
1650 b_new = bracket(:, 1)
1652 b_new = bracket(:, 2)
1655 if (mg%lsf(b_new) * lsf_a > 0)
then
1658 r_root = bisection(mg%lsf, a, b_new, mg%lsf_tol, max_iter)
1662 dist = norm2(r_root - a)/norm2(b-a)
1663 dist = max(dist, mg%lsf_min_rel_distance)
1664 end function mg_lsf_dist_gss
1667 function bisection(f, in_a, in_b, tol, max_iter)
result(c)
1668 procedure(mg_func_lsf) :: f
1669 real(dp),
intent(in) :: in_a(ndim), in_b(ndim), tol
1670 integer,
intent(in) :: max_iter
1671 real(dp) :: a(ndim), b(ndim), c(ndim), fc
1680 if (0.5 * norm2(b-a) < tol .or. abs(fc) <= 0)
exit
1682 if (fc * f(a) >= 0)
then
1688 end function bisection
1694 function gss(f, in_a, in_b, minimization, tol, find_bracket)
result(bracket)
1695 procedure(mg_func_lsf) :: f
1696 real(dp),
intent(in) :: in_a(ndim)
1697 real(dp),
intent(in) :: in_b(ndim)
1699 logical,
intent(in) :: minimization
1700 real(dp),
intent(in) :: tol
1701 logical,
intent(in) :: find_bracket
1702 real(dp) :: bracket(ndim, 2)
1703 real(dp) :: a(ndim), b(ndim), c(ndim), d(ndim)
1704 real(dp) :: h(ndim), yc, yd, ya
1705 real(dp) :: invphi, invphi2
1708 invphi = (sqrt(5.0_dp) - 1) / 2
1709 invphi2 = (3 - sqrt(5.0_dp)) / 2
1715 if (norm2(h) <= tol)
then
1722 n = int(ceiling(log(tol / norm2(h)) / log(invphi)))
1731 if ((yc < yd) .eqv. minimization)
then
1748 if (find_bracket .and. ya * yc <= 0 .and. ya * yd <= 0)
exit
1751 if ((yc < yd) .eqv. minimization)
then
1762 subroutine mg_box_lsf_stencil(box, mg, ix)
1763 type(box_t),
intent(inout) :: box
1764 type(mg_t),
intent(in) :: mg
1765 integer,
intent(in) :: ix
1766 integer :: IJK, n, nc, idim
1767 integer :: s_ix(NDIM), ix_dist
1768 real(dp) :: dd(2*NDIM), dr2(NDIM)
1769 real(dp),
allocatable :: all_distances(:, DTIMES(:))
1777 ix_dist = af_stencil_index(box, mg_lsf_distance_key)
1778 if (ix_dist == af_stencil_none) error stop
"No distances stored"
1781 box%stencils(ix)%shape = af_stencil_357
1782 box%stencils(ix)%stype = stencil_variable
1784 box%stencils(ix)%cylindrical_gradient = .false.
1785 call af_stencil_allocate_coeff(box%stencils(ix), box%n_cell, use_f=.true.)
1786 box%stencils(ix)%f = 0.0_dp
1788 allocate(all_distances(2*ndim, dtimes(nc)))
1791 all_distances = 1.0_dp
1794 do n = 1,
size(box%stencils(ix_dist)%sparse_ix, 2)
1795 s_ix = box%stencils(ix_dist)%sparse_ix(:, n)
1796 all_distances(:, dindex(s_ix)) = &
1797 box%stencils(ix_dist)%sparse_v(:, n)
1801 dd = all_distances(:, ijk)
1805 box%stencils(ix)%v(1+2*idim-1, ijk) = 1 / &
1806 (0.5_dp * dr2(idim) * (dd(2*idim-1) + dd(2*idim)) * &
1808 box%stencils(ix)%v(1+2*idim, ijk) = 1 / &
1809 (0.5_dp * dr2(idim) * (dd(2*idim-1) + dd(2*idim)) * &
1814 if (box%coord_t == af_cyl)
then
1817 tmp = 1/(box%dr(1) * (dd(1) + dd(2)) * af_cyl_radius_cc(box, i))
1818 box%stencils(ix)%v(2, ijk) = box%stencils(ix)%v(2, ijk) - tmp
1819 box%stencils(ix)%v(3, ijk) = box%stencils(ix)%v(3, ijk) + tmp
1822 box%stencils(ix)%v(1, ijk) = -sum(box%stencils(ix)%v(2:, ijk))
1826 if (dd(n) < 1.0_dp)
then
1827 box%stencils(ix)%f(ijk) = box%stencils(ix)%f(ijk) - &
1828 box%stencils(ix)%v(n+1, ijk)
1829 box%stencils(ix)%v(n+1, ijk) = 0.0_dp
1834 end subroutine mg_box_lsf_stencil
1837 subroutine mg_compute_phi_gradient(tree, mg, i_fc, fac, i_norm)
1838 type(af_t),
intent(inout) :: tree
1839 type(mg_t),
intent(in) :: mg
1840 integer,
intent(in) :: i_fc
1841 real(dp),
intent(in) :: fac
1843 integer,
intent(in),
optional :: i_norm
1844 integer :: lvl, i, id
1847 do lvl = 1, tree%highest_lvl
1849 do i = 1,
size(tree%lvls(lvl)%ids)
1850 id = tree%lvls(lvl)%ids(i)
1851 select case (iand(tree%boxes(id)%tag, mg%operator_mask))
1852 case (mg_normal_box, mg_ceps_box)
1853 call mg_box_lpl_gradient(tree, id, mg, i_fc, fac)
1854 case (mg_lsf_box, mg_ceps_box+mg_lsf_box, mg_veps_box+mg_lsf_box)
1856 if (af_has_children(tree%boxes(id)))
then
1859 call mg_box_lpl_gradient(tree, id, mg, i_fc, fac)
1861 call mg_box_lpllsf_gradient(tree, id, mg, i_fc, fac)
1865 call mg_box_lpl_gradient(tree, id, mg, i_fc, fac)
1867 error stop
"box tag not set"
1869 error stop
"unknown box tag"
1872 if (
present(i_norm))
then
1873 call mg_box_field_norm(tree, id, i_fc, i_norm)
1879 end subroutine mg_compute_phi_gradient
1882 subroutine mg_box_lpl_gradient(tree, id, mg, i_fc, fac)
1883 type(af_t),
intent(inout) :: tree
1884 integer,
intent(in) :: id
1885 type(mg_t),
intent(in) :: mg
1886 integer,
intent(in) :: i_fc
1887 real(dp),
intent(in) :: fac
1888 integer :: nc, i_phi, i_eps
1889 real(dp) :: inv_dr(ndim)
1891 associate(box => tree%boxes(id), cc => tree%boxes(id)%cc)
1894 inv_dr = fac / box%dr
1897 box%fc(1:nc+1, 1, i_fc) = inv_dr(1) * &
1898 (cc(1:nc+1, i_phi) - cc(0:nc, i_phi))
1900 box%fc(1:nc+1, 1:nc, 1, i_fc) = inv_dr(1) * &
1901 (cc(1:nc+1, 1:nc, i_phi) - cc(0:nc, 1:nc, i_phi))
1902 box%fc(1:nc, 1:nc+1, 2, i_fc) = inv_dr(2) * &
1903 (cc(1:nc, 1:nc+1, i_phi) - cc(1:nc, 0:nc, i_phi))
1905 box%fc(1:nc+1, 1:nc, 1:nc, 1, i_fc) = inv_dr(1) * &
1906 (cc(1:nc+1, 1:nc, 1:nc, i_phi) - &
1907 cc(0:nc, 1:nc, 1:nc, i_phi))
1908 box%fc(1:nc, 1:nc+1, 1:nc, 2, i_fc) = inv_dr(2) * &
1909 (cc(1:nc, 1:nc+1, 1:nc, i_phi) - &
1910 cc(1:nc, 0:nc, 1:nc, i_phi))
1911 box%fc(1:nc, 1:nc, 1:nc+1, 3, i_fc) = inv_dr(3) * &
1912 (cc(1:nc, 1:nc, 1:nc+1, i_phi) - &
1913 cc(1:nc, 1:nc, 0:nc, i_phi))
1916 if (iand(box%tag, mg%operator_mask) == mg_veps_box)
then
1921 box%fc(1, 1, i_fc) = 2 * inv_dr(1) * &
1922 (cc(1, i_phi) - cc(0, i_phi)) * &
1924 (cc(1, i_eps) + cc(0, i_eps))
1925 box%fc(nc+1, 1, i_fc) = 2 * inv_dr(1) * &
1926 (cc(nc+1, i_phi) - cc(nc, i_phi)) * &
1928 (cc(nc+1, i_eps) + cc(nc, i_eps))
1930 box%fc(1, 1:nc, 1, i_fc) = 2 * inv_dr(1) * &
1931 (cc(1, 1:nc, i_phi) - cc(0, 1:nc, i_phi)) * &
1932 cc(0, 1:nc, i_eps) / &
1933 (cc(1, 1:nc, i_eps) + cc(0, 1:nc, i_eps))
1934 box%fc(nc+1, 1:nc, 1, i_fc) = 2 * inv_dr(1) * &
1935 (cc(nc+1, 1:nc, i_phi) - cc(nc, 1:nc, i_phi)) * &
1936 cc(nc+1, 1:nc, i_eps) / &
1937 (cc(nc+1, 1:nc, i_eps) + cc(nc, 1:nc, i_eps))
1938 box%fc(1:nc, 1, 2, i_fc) = 2 * inv_dr(2) * &
1939 (cc(1:nc, 1, i_phi) - cc(1:nc, 0, i_phi)) * &
1940 cc(1:nc, 0, i_eps) / &
1941 (cc(1:nc, 1, i_eps) + cc(1:nc, 0, i_eps))
1942 box%fc(1:nc, nc+1, 2, i_fc) = 2 * inv_dr(2) * &
1943 (cc(1:nc, nc+1, i_phi) - cc(1:nc, nc, i_phi)) * &
1944 cc(1:nc, nc+1, i_eps) / &
1945 (cc(1:nc, nc+1, i_eps) + cc(1:nc, nc, i_eps))
1947 box%fc(1, 1:nc, 1:nc, 1, i_fc) = 2 * inv_dr(1) * &
1948 (cc(1, 1:nc, 1:nc, i_phi) - cc(0, 1:nc, 1:nc, i_phi)) * &
1949 cc(0, 1:nc, 1:nc, i_eps) / &
1950 (cc(1, 1:nc, 1:nc, i_eps) + cc(0, 1:nc, 1:nc, i_eps))
1951 box%fc(nc+1, 1:nc, 1:nc, 1, i_fc) = 2 * inv_dr(1) * &
1952 (cc(nc+1, 1:nc, 1:nc, i_phi) - cc(nc, 1:nc, 1:nc, i_phi)) * &
1953 cc(nc+1, 1:nc, 1:nc, i_eps) / &
1954 (cc(nc+1, 1:nc, 1:nc, i_eps) + cc(nc, 1:nc, 1:nc, i_eps))
1955 box%fc(1:nc, 1, 1:nc, 2, i_fc) = 2 * inv_dr(2) * &
1956 (cc(1:nc, 1, 1:nc, i_phi) - cc(1:nc, 0, 1:nc, i_phi)) * &
1957 cc(1:nc, 0, 1:nc, i_eps) / &
1958 (cc(1:nc, 1, 1:nc, i_eps) + cc(1:nc, 0, 1:nc, i_eps))
1959 box%fc(1:nc, nc+1, 1:nc, 2, i_fc) = 2 * inv_dr(2) * &
1960 (cc(1:nc, nc+1, 1:nc, i_phi) - cc(1:nc, nc, 1:nc, i_phi)) * &
1961 cc(1:nc, nc+1, 1:nc, i_eps) / &
1962 (cc(1:nc, nc+1, 1:nc, i_eps) + cc(1:nc, nc, 1:nc, i_eps))
1963 box%fc(1:nc, 1:nc, 1, 3, i_fc) = 2 * inv_dr(3) * &
1964 (cc(1:nc, 1:nc, 1, i_phi) - cc(1:nc, 1:nc, 0, i_phi)) * &
1965 cc(1:nc, 1:nc, 0, i_eps) / &
1966 (cc(1:nc, 1:nc, 1, i_eps) + cc(1:nc, 1:nc, 0, i_eps))
1967 box%fc(1:nc, 1:nc, nc+1, 3, i_fc) = 2 * inv_dr(3) * &
1968 (cc(1:nc, 1:nc, nc+1, i_phi) - cc(1:nc, 1:nc, nc, i_phi)) * &
1969 cc(1:nc, 1:nc, nc+1, i_eps) / &
1970 (cc(1:nc, 1:nc, nc+1, i_eps) + cc(1:nc, 1:nc, nc, i_eps))
1974 end subroutine mg_box_lpl_gradient
1977 subroutine mg_compute_field_norm(tree, i_fc, i_norm)
1978 type(af_t),
intent(inout) :: tree
1979 integer,
intent(in) :: i_fc
1980 integer,
intent(in) :: i_norm
1981 integer :: lvl, i, id
1984 do lvl = 1, tree%highest_lvl
1986 do i = 1,
size(tree%lvls(lvl)%ids)
1987 id = tree%lvls(lvl)%ids(i)
1988 call mg_box_field_norm(tree, id, i_fc, i_norm)
1993 end subroutine mg_compute_field_norm
1995 subroutine mg_box_field_norm(tree, id, i_fc, i_norm)
1996 type(af_t),
intent(inout) :: tree
1997 integer,
intent(in) :: id
1998 integer,
intent(in) :: i_fc
2000 integer,
intent(in) :: i_norm
2003 associate(box => tree%boxes(id))
2006 box%cc(1:nc, i_norm) = 0.5_dp * sqrt(&
2007 (box%fc(1:nc, 1, i_fc) + &
2008 box%fc(2:nc+1, 1, i_fc))**2)
2010 box%cc(1:nc, 1:nc, i_norm) = 0.5_dp * sqrt(&
2011 (box%fc(1:nc, 1:nc, 1, i_fc) + &
2012 box%fc(2:nc+1, 1:nc, 1, i_fc))**2 + &
2013 (box%fc(1:nc, 1:nc, 2, i_fc) + &
2014 box%fc(1:nc, 2:nc+1, 2, i_fc))**2)
2016 box%cc(1:nc, 1:nc, 1:nc, i_norm) = 0.5_dp * sqrt(&
2017 (box%fc(1:nc, 1:nc, 1:nc, 1, i_fc) + &
2018 box%fc(2:nc+1, 1:nc, 1:nc, 1, i_fc))**2 + &
2019 (box%fc(1:nc, 1:nc, 1:nc, 2, i_fc) + &
2020 box%fc(1:nc, 2:nc+1, 1:nc, 2, i_fc))**2 + &
2021 (box%fc(1:nc, 1:nc, 1:nc, 3, i_fc) + &
2022 box%fc(1:nc, 1:nc, 2:nc+1, 3, i_fc))**2)
2025 end subroutine mg_box_field_norm
2030 subroutine mg_box_lpllsf_gradient(tree, id, mg, i_fc, fac)
2031 type(af_t),
intent(inout) :: tree
2032 integer,
intent(in) :: id
2033 type(mg_t),
intent(in) :: mg
2034 integer,
intent(in) :: i_fc
2035 real(dp),
intent(in) :: fac
2037 integer :: n, nc, i_phi, ix_dist
2038 real(dp) :: inv_dr(NDIM), dd(2*NDIM)
2039 real(dp) :: bc(DTIMES(tree%n_cell))
2042 call mg_box_lpl_gradient(tree, id, mg, i_fc, fac)
2044 ix_dist = af_stencil_index(tree%boxes(id), mg_lsf_distance_key)
2045 if (ix_dist == af_stencil_none) error stop
"No distances stored"
2047 associate(box => tree%boxes(id), cc => tree%boxes(id)%cc)
2050 inv_dr = fac / box%dr
2051 bc = mg_lsf_boundary_value(box, mg)
2054 do n = 1,
size(box%stencils(ix_dist)%sparse_ix, 2)
2055 dd = box%stencils(ix_dist)%sparse_v(:, n)
2058 i = box%stencils(ix_dist)%sparse_ix(1, n)
2060 if (dd(1) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2061 box%fc(i, 1, i_fc) = inv_dr(1) * &
2062 (cc(i, i_phi) - bc(ijk)) / dd(1)
2064 if (dd(2) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2065 box%fc(i+1, 1, i_fc) = inv_dr(1) * &
2066 (bc(ijk) - cc(i, i_phi)) / dd(2)
2069 i = box%stencils(ix_dist)%sparse_ix(1, n)
2070 j = box%stencils(ix_dist)%sparse_ix(2, n)
2072 if (dd(1) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2073 box%fc(i, j, 1, i_fc) = inv_dr(1) * &
2074 (cc(i, j, i_phi) - bc(ijk)) / dd(1)
2076 if (dd(2) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2077 box%fc(i+1, j, 1, i_fc) = inv_dr(1) * &
2078 (bc(ijk) - cc(i, j, i_phi)) / dd(2)
2080 if (dd(3) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2081 box%fc(i, j, 2, i_fc) = inv_dr(2) * &
2082 (cc(i, j, i_phi) - bc(ijk)) / dd(3)
2084 if (dd(4) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2085 box%fc(i, j+1, 2, i_fc) = inv_dr(2) * &
2086 (bc(ijk) - cc(i, j, i_phi)) / dd(4)
2089 i = box%stencils(ix_dist)%sparse_ix(1, n)
2090 j = box%stencils(ix_dist)%sparse_ix(2, n)
2091 k = box%stencils(ix_dist)%sparse_ix(3, n)
2093 if (dd(1) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2094 box%fc(i, j, k, 1, i_fc) = inv_dr(1) * &
2095 (cc(ijk, i_phi) - bc(ijk)) / dd(1)
2097 if (dd(2) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2098 box%fc(i+1, j, k, 1, i_fc) = inv_dr(1) * &
2099 (bc(ijk) - cc(ijk, i_phi)) / dd(2)
2101 if (dd(3) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2102 box%fc(i, j, k, 2, i_fc) = inv_dr(2) * &
2103 (cc(ijk, i_phi) - bc(ijk)) / dd(3)
2105 if (dd(4) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2106 box%fc(i, j+1, k, 2, i_fc) = inv_dr(2) * &
2107 (bc(ijk) - cc(ijk, i_phi)) / dd(4)
2109 if (dd(5) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2110 box%fc(i, j, k, 3, i_fc) = inv_dr(3) * &
2111 (cc(ijk, i_phi) - bc(ijk)) / dd(5)
2113 if (dd(6) < 1 .and. cc(ijk, mg%i_lsf) >= 0)
then
2114 box%fc(i, j, k+1, 3, i_fc) = inv_dr(3) * &
2115 (bc(ijk) - cc(ijk, i_phi)) / dd(6)
2120 end subroutine mg_box_lpllsf_gradient
2124 subroutine check_coarse_representation_lsf(tree, mg)
2125 type(af_t),
intent(in) :: tree
2126 type(mg_t),
intent(in) :: mg
2127 integer :: i, id, ix
2129 do i = 1,
size(tree%lvls(1)%ids)
2130 id = tree%lvls(1)%ids(i)
2131 ix = af_stencil_index(tree%boxes(id), mg_lsf_distance_key)
2132 if (ix > af_stencil_none)
exit
2135 if (i ==
size(tree%lvls(1)%ids)+1)
then
2136 print *,
"No roots found for level set function on coarse grid"
2137 print *,
"you should probably use a finer coarse grid"
2138 error stop
"level set function not resolved on coarse grid"
2141 end subroutine check_coarse_representation_lsf
2144 function numerical_gradient(f, r)
result(gradient)
2145 procedure(mg_func_lsf) :: f
2146 real(dp),
intent(in) :: r(ndim)
2147 real(dp),
parameter :: sqrteps = sqrt(epsilon(1.0_dp))
2148 real(dp),
parameter :: min_stepsize = epsilon(1.0_dp)
2149 real(dp) :: r_eval(ndim), gradient(ndim)
2150 real(dp) :: stepsize(ndim), flo, fhi
2153 stepsize = max(min_stepsize, sqrteps * abs(r))
2158 r_eval(idim) = r(idim) - stepsize(idim)
2161 r_eval(idim) = r(idim) + stepsize(idim)
2165 gradient(idim) = (fhi - flo)/(2 * stepsize(idim))
2168 r_eval(idim) = r(idim)
2170 end function numerical_gradient
This module contains the core routines of Afivo, namely those that deal with initializing and changin...
This module contains routines related to the filling of ghost cells. Note that corner ghost cells are...
This module contains the geometric multigrid routines that come with Afivo.
This module contains routines for restriction: going from fine to coarse variables.
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...
This module contains all kinds of different 'helper' routines for Afivo. If the number of routines fo...
Module to solve elliptic PDEs on the coarse grid. This module contains an interface to Hypre,...
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 to store multigrid options in.