30 integer,
parameter,
public ::
dp = kind(0.0d0)
33 integer,
parameter,
public ::
i8 = selected_int_kind(18)
111 integer,
parameter,
public ::
mg_child_dix(2, 4) = reshape([0,0,1,0,0,1,1,1], [2,4])
113 integer,
parameter,
public ::
mg_child_rev(4, 2) = reshape([2,1,4,3,3,4,1,2], [4,2])
115 integer,
parameter,
public ::
mg_child_adj_nb(2, 4) = reshape([1,3,2,4,1,2,3,4], [2,4])
117 logical,
parameter,
public ::
mg_child_low(2, 4) = reshape([.true., .true., &
118 .false., .true., .true., .false., .false., .false.], [2, 4])
128 integer,
parameter,
public ::
mg_neighb_dix(2, 4) = reshape([-1,0,1,0,0,-1,0,1], [2,4])
130 logical,
parameter,
public ::
mg_neighb_low(4) = [.true., .false., .true., .false.]
141 integer,
allocatable :: leaves(:)
142 integer,
allocatable :: parents(:)
143 integer,
allocatable :: ref_bnds(:)
144 integer,
allocatable :: ids(:)
145 integer,
allocatable :: my_leaves(:)
146 integer,
allocatable :: my_parents(:)
147 integer,
allocatable :: my_ref_bnds(:)
148 integer,
allocatable :: my_ids(:)
158 integer :: children(2**2)
159 integer :: neighbors(2*2)
163 real(
dp),
allocatable :: cc(:, :, :)
171 integer,
allocatable :: ix(:)
172 real(
dp),
allocatable :: send(:)
173 real(
dp),
allocatable :: recv(:)
177 integer,
allocatable :: n_send(:, :)
178 integer,
allocatable :: n_recv(:, :)
183 real(
dp) :: bc_value = 0.0_dp
185 procedure(
mg_subr_bc),
pointer,
nopass :: boundary_cond => null()
187 procedure(mg_subr_rb),
pointer,
nopass :: refinement_bnd => null()
191 character(len=20) :: name
192 real(
dp) :: t = 0.0_dp
198 logical :: tree_created = .false.
200 logical :: is_allocated = .false.
202 integer :: n_extra_vars = 0
206 integer :: n_cpu = -1
208 integer :: my_rank = -1
210 integer :: box_size = -1
212 integer :: highest_lvl = -1
214 integer :: lowest_lvl = -1
217 integer :: first_normal_lvl = -1
219 integer :: n_boxes = 0
244 logical :: phi_bc_data_stored = .false.
247 logical :: periodic(2) = .false.
262 integer :: n_smoother_substeps = 1
264 integer :: n_cycle_down = 2
266 integer :: n_cycle_up = 2
268 integer :: max_coarse_cycles = 1000
269 integer :: coarsest_grid(2) = 2
271 real(
dp) :: residual_coarse_abs = 1e-8_dp
273 real(
dp) :: residual_coarse_rel = 1e-8_dp
276 procedure(mg_box_op),
pointer,
nopass :: box_op => null()
279 procedure(mg_box_gsrb),
pointer,
nopass :: box_smoother => null()
282 procedure(mg_box_prolong),
pointer,
nopass :: box_prolong => null()
285 integer :: n_timers = 0
294 type(mg_box_t),
intent(in) :: box
295 integer,
intent(in) :: nc
296 integer,
intent(in) :: iv
297 integer,
intent(in) :: nb
298 integer,
intent(out) :: bc_type
300 real(dp),
intent(out) :: bc(nc)
304 subroutine mg_subr_rb(box, nc, iv, nb, cgc)
306 type(mg_box_t),
intent(inout) :: box
307 integer,
intent(in) :: nc
308 integer,
intent(in) :: iv
309 integer,
intent(in) :: nb
311 real(dp),
intent(in) :: cgc(nc)
312 end subroutine mg_subr_rb
315 subroutine mg_box_op(mg, id, nc, i_out)
317 type(mg_t),
intent(inout) :: mg
318 integer,
intent(in) :: id
319 integer,
intent(in) :: nc
320 integer,
intent(in) :: i_out
321 end subroutine mg_box_op
324 subroutine mg_box_gsrb(mg, id, nc, redblack_cntr)
326 type(mg_t),
intent(inout) :: mg
327 integer,
intent(in) :: id
328 integer,
intent(in) :: nc
329 integer,
intent(in) :: redblack_cntr
330 end subroutine mg_box_gsrb
333 subroutine mg_box_prolong(mg, p_id, dix, nc, iv, fine)
335 type(mg_t),
intent(inout) :: mg
336 integer,
intent(in) :: p_id
337 integer,
intent(in) :: dix(2)
338 integer,
intent(in) :: nc
339 integer,
intent(in) :: iv
340 real(dp),
intent(out) :: fine(nc, nc)
341 end subroutine mg_box_prolong
348 public :: mg_box_gsrb
349 public :: mg_box_prolong
453 integer :: timer_total_vcycle = -1
454 integer :: timer_total_fmg = -1
455 integer :: timer_smoother = -1
456 integer :: timer_smoother_gc = -1
457 integer :: timer_coarse = -1
458 integer :: timer_correct = -1
459 integer :: timer_update_coarse = -1
476 Integer,
Parameter :: kdp = selected_real_kind(15)
481 module procedure i_mrgrnk
511 integer,
intent(in) :: ix(2)
520 type(
mg_t),
intent(in) :: mg
521 integer,
intent(in) :: id
522 integer :: ix_offset(2)
524 if (mg%boxes(id)%lvl <= mg%first_normal_lvl)
then
527 ix_offset = iand(mg%boxes(id)%ix-1, 1) * &
528 ishft(mg%box_size, -1)
533 type(
mg_t),
intent(in) :: mg
536 do lvl = mg%first_normal_lvl, mg%highest_lvl-1
538 if (
size(mg%lvls(lvl)%leaves) /= 0 .and. &
539 size(mg%lvls(lvl)%parents) /= 0)
exit
546 type(
mg_t),
intent(in) :: mg
548 integer(i8) :: n_unknowns
551 do lvl = mg%first_normal_lvl, mg%highest_lvl
552 n_unknowns = n_unknowns +
size(mg%lvls(lvl)%leaves)
554 n_unknowns = n_unknowns * int(mg%box_size**3,
i8)
560 integer,
intent(in) :: nb
561 integer,
intent(in) :: nc
562 real(
dp),
intent(out) :: x(nc, 2)
563 integer :: i, ixs(2-1)
570 ixs = [(i, i = 1, 2-1)]
571 ixs(nb_dim:) = ixs(nb_dim:) + 1
575 rmin(nb_dim) = rmin(nb_dim) + box%dr(nb_dim) * nc
580 x(i, ixs(1)) = x(i, ixs(1)) + (i-0.5d0) * box%dr(ixs(1))
585 type(
mg_t),
intent(inout) :: mg
586 character(len=*),
intent(in) :: name
588 mg%n_timers = mg%n_timers + 1
596 timer%t0 = mpi_wtime()
602 timer%t = timer%t + mpi_wtime() - timer%t0
607 type(
mg_t),
intent(in) :: mg
609 real(
dp) :: tmin(mg%n_timers)
610 real(
dp) :: tmax(mg%n_timers)
612 call mpi_reduce(mg%timers(1:mg%n_timers)%t, tmin, mg%n_timers, &
613 mpi_double, mpi_min, 0, mg%comm, ierr)
614 call mpi_reduce(mg%timers(1:mg%n_timers)%t, tmax, mg%n_timers, &
615 mpi_double, mpi_max, 0, mg%comm, ierr)
617 if (mg%my_rank == 0)
then
618 write(*,
"(A20,2A16)")
"name ",
"min(s)",
"max(s)"
619 do n = 1, mg%n_timers
620 write(*,
"(A20,2F16.6)") mg%timers(n)%name, &
630 type(
mg_t),
intent(inout) :: mg
633 if (.not. mg%is_allocated) &
634 error stop
"deallocate_storage: tree is not allocated"
639 deallocate(mg%comm_restrict%n_send)
640 deallocate(mg%comm_restrict%n_recv)
642 deallocate(mg%comm_prolong%n_send)
643 deallocate(mg%comm_prolong%n_recv)
645 deallocate(mg%comm_ghostcell%n_send)
646 deallocate(mg%comm_ghostcell%n_recv)
648 do lvl = mg%lowest_lvl, mg%highest_lvl
649 deallocate(mg%lvls(lvl)%ids)
650 deallocate(mg%lvls(lvl)%leaves)
651 deallocate(mg%lvls(lvl)%parents)
652 deallocate(mg%lvls(lvl)%ref_bnds)
653 deallocate(mg%lvls(lvl)%my_ids)
654 deallocate(mg%lvls(lvl)%my_leaves)
655 deallocate(mg%lvls(lvl)%my_parents)
656 deallocate(mg%lvls(lvl)%my_ref_bnds)
659 mg%is_allocated = .false.
661 mg%phi_bc_data_stored = .false.
668 type(
mg_t),
intent(inout) :: mg
669 integer :: i, id, lvl, nc
670 integer :: n_send(0:mg%n_cpu-1, 3)
671 integer :: n_recv(0:mg%n_cpu-1, 3)
673 integer :: n_in, n_out, n_id
675 if (.not. mg%tree_created) &
676 error stop
"allocate_storage: tree is not yet created"
678 if (mg%is_allocated) &
679 error stop
"allocate_storage: tree is already allocated"
681 do lvl = mg%lowest_lvl, mg%highest_lvl
682 nc = mg%box_size_lvl(lvl)
683 do i = 1,
size(mg%lvls(lvl)%my_ids)
684 id = mg%lvls(lvl)%my_ids(i)
685 allocate(mg%boxes(id)%cc(0:nc+1, 0:nc+1, &
689 mg%boxes(id)%cc(:, :, :) = 0.0_dp
693 allocate(mg%buf(0:mg%n_cpu-1))
696 n_recv(:, 1), dsize(1))
698 n_recv(:, 2), dsize(2))
700 n_recv(:, 3), dsize(3))
703 n_out = maxval(n_send(i, :) * dsize(:))
704 n_in = maxval(n_recv(i, :) * dsize(:))
705 n_id = maxval(n_send(i, :))
706 allocate(mg%buf(i)%send(n_out))
707 allocate(mg%buf(i)%recv(n_in))
708 allocate(mg%buf(i)%ix(n_id))
711 mg%is_allocated = .true.
721 type(
mg_t),
intent(inout) :: mg
722 real(
dp),
intent(in) :: dt
723 real(
dp),
intent(in) :: diffusion_coeff
724 integer,
intent(in) :: order
725 real(
dp),
intent(in) :: max_res
726 integer,
parameter :: max_its = 10
736 call set_rhs(mg, -1/(dt * diffusion_coeff), 0.0_dp)
741 call set_rhs(mg, -2/(dt * diffusion_coeff), -1.0_dp)
743 error stop
"diffusion_solve: order should be 1 or 2"
751 if (res <= max_res)
exit
755 if (n == max_its + 1)
then
756 if (mg%my_rank == 0) &
757 print *,
"Did you specify boundary conditions correctly?"
758 error stop
"diffusion_solve: no convergence"
768 type(
mg_t),
intent(inout) :: mg
769 real(
dp),
intent(in) :: dt
770 integer,
intent(in) :: order
771 real(
dp),
intent(in) :: max_res
772 integer,
parameter :: max_its = 10
782 call set_rhs(mg, -1/dt, 0.0_dp)
787 call set_rhs(mg, -2/dt, -1.0_dp)
789 error stop
"diffusion_solve: order should be 1 or 2"
797 if (res <= max_res)
exit
801 if (n == max_its + 1)
then
802 if (mg%my_rank == 0)
then
803 print *,
"Did you specify boundary conditions correctly?"
804 print *,
"Or is the variation in diffusion too large?"
806 error stop
"diffusion_solve: no convergence"
817 type(
mg_t),
intent(inout) :: mg
818 real(
dp),
intent(in) :: dt
819 integer,
intent(in) :: order
820 real(
dp),
intent(in) :: max_res
821 integer,
parameter :: max_its = 100
831 call set_rhs(mg, -1/dt, 0.0_dp)
836 call set_rhs(mg, -2/dt, -1.0_dp)
838 error stop
"diffusion_solve: order should be 1 or 2"
847 if (res <= max_res)
exit
851 if (n == max_its + 1)
then
852 if (mg%my_rank == 0)
then
853 print *,
"Did you specify boundary conditions correctly?"
854 print *,
"Or is the variation in diffusion too large?"
856 error stop
"diffusion_solve: no convergence"
860 subroutine set_rhs(mg, f1, f2)
861 type(
mg_t),
intent(inout) :: mg
862 real(
dp),
intent(in) :: f1, f2
863 integer :: lvl, i, id, nc
865 do lvl = 1, mg%highest_lvl
866 nc = mg%box_size_lvl(lvl)
867 do i = 1,
size(mg%lvls(lvl)%my_leaves)
868 id = mg%lvls(lvl)%my_leaves(i)
869 mg%boxes(id)%cc(1:nc, 1:nc,
mg_irhs) = &
870 f1 * mg%boxes(id)%cc(1:nc, 1:nc,
mg_iphi) + &
871 f2 * mg%boxes(id)%cc(1:nc, 1:nc,
mg_irhs)
874 end subroutine set_rhs
879 type(
mg_t),
intent(inout) :: mg
881 if (all(mg%periodic))
then
884 mg%subtract_mean = .true.
887 select case (mg%geometry_type)
891 select case (mg%smoother_type)
895 mg%box_smoother => box_gs_lpl
897 error stop
"laplacian_set_methods: unsupported smoother type"
900 mg%box_op => box_clpl
902 select case (mg%smoother_type)
904 mg%box_smoother => box_gs_clpl
906 error stop
"laplacian_set_methods: unsupported smoother type"
909 error stop
"laplacian_set_methods: unsupported geometry"
915 subroutine box_gs_lpl(mg, id, nc, redblack_cntr)
916 type(
mg_t),
intent(inout) :: mg
917 integer,
intent(in) :: id
918 integer,
intent(in) :: nc
919 integer,
intent(in) :: redblack_cntr
920 integer :: i, j, i0, di
921 real(
dp) :: idr2(2), fac
924 idr2 = 1/mg%dr(:, mg%boxes(id)%lvl)**2
925 fac = 0.5_dp / sum(idr2)
937 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
940 i0 = 2 - iand(ieor(redblack_cntr, j), 1)
943 cc(i, j, n) = fac * ( &
944 idr2(1) * (cc(i+1, j, n) + cc(i-1, j, n)) + &
945 idr2(2) * (cc(i, j+1, n) + cc(i, j-1, n)) - &
950 end subroutine box_gs_lpl
991 subroutine box_lpl(mg, id, nc, i_out)
992 type(
mg_t),
intent(inout) :: mg
993 integer,
intent(in) :: id
994 integer,
intent(in) :: nc
995 integer,
intent(in) :: i_out
999 idr2 = 1 / mg%dr(:, mg%boxes(id)%lvl)**2
1001 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
1005 idr2(1) * (cc(i-1, j, n) + cc(i+1, j, n) - 2 * cc(i, j, n)) + &
1006 idr2(2) * (cc(i, j-1, n) + cc(i, j+1, n) - 2 * cc(i, j, n))
1010 end subroutine box_lpl
1014 subroutine box_clpl(mg, id, nc, i_out)
1015 type(
mg_t),
intent(inout) :: mg
1016 integer,
intent(in) :: id
1017 integer,
intent(in) :: nc
1018 integer,
intent(in) :: i_out
1020 real(
dp) :: dr(2), idr2(2)
1021 real(
dp) :: r_face(nc+1), r_inv(nc)
1023 dr = mg%dr(:, mg%boxes(id)%lvl)
1025 r_face = mg%boxes(id)%r_min(1) + dr(1) * [(i, i=0,nc)]
1026 r_inv = 1/(mg%boxes(id)%r_min(1) + dr(1) * [(i-0.5_dp, i=1,nc)])
1028 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
1031 cc(i, j, i_out) = idr2(1) * (&
1032 r_face(i) * r_inv(i) * cc(i-1, j, n) + &
1033 r_face(i+1) * r_inv(i) * cc(i+1, j, n) &
1034 - 2 * cc(i, j, n)) &
1035 + idr2(2) * (cc(i, j-1, n) + cc(i, j+1, n) - 2 * cc(i, j, n))
1039 end subroutine box_clpl
1044 subroutine box_gs_clpl(mg, id, nc, redblack_cntr)
1045 type(
mg_t),
intent(inout) :: mg
1046 integer,
intent(in) :: id
1047 integer,
intent(in) :: nc
1048 integer,
intent(in) :: redblack_cntr
1049 integer :: i, j, i0, di
1051 real(
dp) :: idr2(2), dr(2), dr2(2), fac
1052 real(
dp) :: r_face(nc+1), r_inv(nc)
1054 dr = mg%dr(:, mg%boxes(id)%lvl)
1057 fac = 0.5_dp / sum(idr2)
1058 r_face = mg%boxes(id)%r_min(1) + dr(1) * [(i, i=0,nc)]
1059 r_inv = 1/(mg%boxes(id)%r_min(1) + dr(1) * [(i-0.5_dp, i=1,nc)])
1071 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
1074 i0 = 2 - iand(ieor(redblack_cntr, j), 1)
1077 cc(i, j, n) = fac * (idr2(1) * &
1078 (r_face(i+1) * r_inv(i) * cc(i+1, j, n) + &
1079 r_face(i) * r_inv(i) * cc(i-1, j, n)) + &
1080 idr2(2) * (cc(i, j+1, n) + cc(i, j-1, n)) &
1085 end subroutine box_gs_clpl
1090 type(
mg_t),
intent(inout) :: mg
1092 if (mg%n_extra_vars == 0 .and. mg%is_allocated)
then
1093 error stop
"vhelmholtz_set_methods: mg%n_extra_vars == 0"
1095 mg%n_extra_vars = max(1, mg%n_extra_vars)
1098 mg%subtract_mean = .false.
1103 mg%bc(:,
mg_iveps)%bc_value = 0.0_dp
1105 select case (mg%geometry_type)
1107 mg%box_op => box_vhelmh
1109 select case (mg%smoother_type)
1111 mg%box_smoother => box_gs_vhelmh
1113 error stop
"vhelmholtz_set_methods: unsupported smoother type"
1116 error stop
"vhelmholtz_set_methods: unsupported geometry"
1122 real(
dp),
intent(in) :: lambda
1125 error stop
"vhelmholtz_set_lambda: lambda < 0 not allowed"
1131 subroutine box_gs_vhelmh(mg, id, nc, redblack_cntr)
1132 type(
mg_t),
intent(inout) :: mg
1133 integer,
intent(in) :: id
1134 integer,
intent(in) :: nc
1135 integer,
intent(in) :: redblack_cntr
1136 integer :: i, j, i0, di
1138 real(
dp) :: idr2(2*2), u(2*2)
1139 real(
dp) :: a0, a(2*2), c(2*2)
1142 idr2(1:2*2:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1143 idr2(2:2*2:2) = idr2(1:2*2:2)
1155 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
1159 i0 = 2 - iand(ieor(redblack_cntr, j), 1)
1162 a0 = cc(i, j, i_eps)
1163 u(1:2) = cc(i-1:i+1:2, j, n)
1164 a(1:2) = cc(i-1:i+1:2, j, i_eps)
1165 u(3:4) = cc(i, j-1:j+1:2, n)
1166 a(3:4) = cc(i, j-1:j+1:2, i_eps)
1167 c(:) = 2 * a0 * a(:) / (a0 + a(:)) * idr2
1170 (sum(c(:) * u(:)) - cc(i, j,
mg_irhs)) / &
1175 end subroutine box_gs_vhelmh
1178 subroutine box_vhelmh(mg, id, nc, i_out)
1179 type(
mg_t),
intent(inout) :: mg
1180 integer,
intent(in) :: id
1181 integer,
intent(in) :: nc
1182 integer,
intent(in) :: i_out
1184 real(
dp) :: idr2(2*2), a0, u0, u(2*2), a(2*2)
1187 idr2(1:2*2:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1188 idr2(2:2*2:2) = idr2(1:2*2:2)
1190 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
1194 a0 = cc(i, j, i_eps)
1195 a(1:2) = cc(i-1:i+1:2, j, i_eps)
1196 a(3:4) = cc(i, j-1:j+1:2, i_eps)
1198 u(1:2) = cc(i-1:i+1:2, j, n)
1199 u(3:4) = cc(i, j-1:j+1:2, n)
1201 cc(i, j, i_out) = sum(2 * idr2 * &
1202 a0*a(:)/(a0 + a(:)) * (u(:) - u0)) - &
1207 end subroutine box_vhelmh
1213 type(
mg_t),
intent(inout) :: mg
1214 integer,
intent(in) :: domain_size(2)
1215 integer,
intent(in) :: box_size
1216 real(
dp),
intent(in) :: dx(2)
1217 real(
dp),
intent(in) :: r_min(2)
1218 logical,
intent(in) :: periodic(2)
1219 integer,
intent(in) :: n_finer
1220 integer :: i, j, lvl, n, id, nx(2)
1221 integer :: boxes_per_dim(2,
mg_lvl_lo:1)
1222 integer :: periodic_offset(2)
1224 if (modulo(box_size, 2) /= 0) &
1225 error stop
"box_size should be even"
1226 if (any(modulo(domain_size, box_size) /= 0)) &
1227 error stop
"box_size does not divide domain_size"
1229 if (all(periodic))
then
1230 mg%subtract_mean = .true.
1234 mg%box_size = box_size
1235 mg%box_size_lvl(1) = box_size
1236 mg%domain_size_lvl(:, 1) = domain_size
1237 mg%first_normal_lvl = 1
1240 mg%periodic = periodic
1241 boxes_per_dim(:, :) = 0
1242 boxes_per_dim(:, 1) = domain_size / box_size
1247 if (any(modulo(nx, 2) == 1 .or. nx == mg%coarsest_grid .or. &
1248 (mg%box_size_lvl(lvl) == mg%coarsest_grid .and. &
1251 if (all(modulo(nx/mg%box_size_lvl(lvl), 2) == 0))
then
1252 mg%box_size_lvl(lvl-1) = mg%box_size_lvl(lvl)
1253 boxes_per_dim(:, lvl-1) = boxes_per_dim(:, lvl)/2
1254 mg%first_normal_lvl = lvl-1
1256 mg%box_size_lvl(lvl-1) = mg%box_size_lvl(lvl)/2
1257 boxes_per_dim(:, lvl-1) = boxes_per_dim(:, lvl)
1260 mg%dr(:, lvl-1) = mg%dr(:, lvl) * 2
1262 mg%domain_size_lvl(:, lvl-1) = nx
1269 mg%dr(:, lvl) = mg%dr(:, lvl-1) * 0.5_dp
1270 mg%box_size_lvl(lvl) = box_size
1271 mg%domain_size_lvl(:, lvl) = 2 * mg%domain_size_lvl(:, lvl-1)
1274 n = sum(product(boxes_per_dim, dim=1)) + n_finer
1275 allocate(mg%boxes(n))
1278 nx = boxes_per_dim(:, mg%lowest_lvl)
1279 periodic_offset = [nx(1)-1, (nx(2)-1)*nx(1)]
1281 do j=1,nx(2);
do i=1,nx(1)
1282 mg%n_boxes = mg%n_boxes + 1
1285 mg%boxes(n)%rank = 0
1287 mg%boxes(n)%lvl = mg%lowest_lvl
1288 mg%boxes(n)%ix(:) = [i, j]
1289 mg%boxes(n)%r_min(:) = r_min + (mg%boxes(n)%ix(:) - 1) * &
1290 mg%box_size_lvl(mg%lowest_lvl) * mg%dr(:, mg%lowest_lvl)
1291 mg%boxes(n)%dr(:) = mg%dr(:, mg%lowest_lvl)
1296 mg%boxes(n)%neighbors(:) = [n-1, n+1, n-nx(1), n+nx(1)]
1299 where ([i, j] == 1 .and. .not. periodic)
1303 where ([i, j] == 1 .and. periodic)
1308 where ([i, j] == nx .and. .not. periodic)
1312 where ([i, j] == nx .and. periodic)
1318 mg%lvls(mg%lowest_lvl)%ids = [(n, n=1, mg%n_boxes)]
1321 do lvl = mg%lowest_lvl, 0
1322 if (mg%box_size_lvl(lvl+1) == mg%box_size_lvl(lvl))
then
1323 do i = 1,
size(mg%lvls(lvl)%ids)
1324 id = mg%lvls(lvl)%ids(i)
1332 do i = 1,
size(mg%lvls(lvl)%ids)
1333 id = mg%lvls(lvl)%ids(i)
1334 call add_single_child(mg, id,
size(mg%lvls(lvl)%ids))
1345 do lvl = mg%lowest_lvl, 1
1346 if (
allocated(mg%lvls(lvl)%ref_bnds)) &
1347 deallocate(mg%lvls(lvl)%ref_bnds)
1348 allocate(mg%lvls(lvl)%ref_bnds(0))
1351 mg%tree_created = .true.
1355 type(
mg_t),
intent(inout) :: mg
1356 integer,
intent(in) :: lvl
1359 do i = 1,
size(mg%lvls(lvl)%ids)
1360 id = mg%lvls(lvl)%ids(i)
1361 call set_neighbs(mg%boxes, id)
1366 type(
mg_t),
intent(inout) :: mg
1367 integer,
intent(in) :: lvl
1370 if (
allocated(mg%lvls(lvl+1)%ids)) &
1371 deallocate(mg%lvls(lvl+1)%ids)
1374 if (mg%box_size_lvl(lvl+1) == mg%box_size_lvl(lvl))
then
1376 allocate(mg%lvls(lvl+1)%ids(n))
1379 do i = 1,
size(mg%lvls(lvl)%parents)
1380 id = mg%lvls(lvl)%parents(i)
1381 mg%lvls(lvl+1)%ids(n*(i-1)+1:n*i) = mg%boxes(id)%children
1384 n =
size(mg%lvls(lvl)%parents)
1385 allocate(mg%lvls(lvl+1)%ids(n))
1388 do i = 1,
size(mg%lvls(lvl)%parents)
1389 id = mg%lvls(lvl)%parents(i)
1390 mg%lvls(lvl+1)%ids(i) = mg%boxes(id)%children(1)
1397 subroutine set_neighbs(boxes, id)
1398 type(
mg_box_t),
intent(inout) :: boxes(:)
1399 integer,
intent(in) :: id
1400 integer :: nb, nb_id
1403 if (boxes(id)%neighbors(nb) ==
mg_no_box)
then
1404 nb_id = find_neighb(boxes, id, nb)
1406 boxes(id)%neighbors(nb) = nb_id
1411 end subroutine set_neighbs
1414 function find_neighb(boxes, id, nb)
result(nb_id)
1415 type(
mg_box_t),
intent(in) :: boxes(:)
1416 integer,
intent(in) :: id
1417 integer,
intent(in) :: nb
1418 integer :: nb_id, p_id, c_ix, d, old_pid
1420 p_id = boxes(id)%parent
1428 p_id = boxes(p_id)%neighbors(nb)
1433 end function find_neighb
1437 type(
mg_box_t),
intent(in) :: boxes(:)
1438 type(
mg_lvl_t),
intent(inout) :: level
1439 integer :: i, id, i_leaf, i_parent
1440 integer :: n_parents, n_leaves
1443 n_leaves =
size(level%ids) - n_parents
1445 if (.not.
allocated(level%parents))
then
1446 allocate(level%parents(n_parents))
1447 else if (n_parents /=
size(level%parents))
then
1448 deallocate(level%parents)
1449 allocate(level%parents(n_parents))
1452 if (.not.
allocated(level%leaves))
then
1453 allocate(level%leaves(n_leaves))
1454 else if (n_leaves /=
size(level%leaves))
then
1455 deallocate(level%leaves)
1456 allocate(level%leaves(n_leaves))
1461 do i = 1,
size(level%ids)
1464 i_parent = i_parent + 1
1465 level%parents(i_parent) = id
1468 level%leaves(i_leaf) = id
1475 type(
mg_box_t),
intent(in) :: boxes(:)
1476 type(
mg_lvl_t),
intent(inout) :: level
1477 integer,
allocatable :: tmp(:)
1478 integer :: i, id, nb, nb_id, ix
1480 if (
allocated(level%ref_bnds))
deallocate(level%ref_bnds)
1482 if (
size(level%parents) == 0)
then
1484 allocate(level%ref_bnds(0))
1486 allocate(tmp(
size(level%leaves)))
1488 do i = 1,
size(level%leaves)
1489 id = level%leaves(i)
1492 nb_id = boxes(id)%neighbors(nb)
1503 allocate(level%ref_bnds(ix))
1504 level%ref_bnds(:) = tmp(1:ix)
1509 type(
mg_t),
intent(inout) :: mg
1510 integer,
intent(in) :: id
1511 integer :: lvl, i, nb, child_nb(2**(2-1))
1513 integer :: c_id, c_ix_base(2)
1516 error stop
"mg_add_children: not enough space"
1521 mg%boxes(id)%children = c_ids
1522 c_ix_base = 2 * mg%boxes(id)%ix - 1
1523 lvl = mg%boxes(id)%lvl+1
1527 mg%boxes(c_id)%rank = mg%boxes(id)%rank
1529 mg%boxes(c_id)%lvl = lvl
1530 mg%boxes(c_id)%parent = id
1533 mg%boxes(c_id)%r_min = mg%boxes(id)%r_min + &
1535 mg%boxes(c_id)%dr(:) = mg%dr(:, lvl)
1540 if (mg%boxes(id)%neighbors(nb) <
mg_no_box)
then
1542 mg%boxes(child_nb)%neighbors(nb) = mg%boxes(id)%neighbors(nb)
1547 subroutine add_single_child(mg, id, n_boxes_lvl)
1548 type(
mg_t),
intent(inout) :: mg
1549 integer,
intent(in) :: id
1550 integer,
intent(in) :: n_boxes_lvl
1551 integer :: lvl, c_id
1553 c_id = mg%n_boxes + 1
1554 mg%n_boxes = mg%n_boxes + 1
1555 mg%boxes(id)%children(1) = c_id
1556 lvl = mg%boxes(id)%lvl+1
1558 mg%boxes(c_id)%rank = mg%boxes(id)%rank
1559 mg%boxes(c_id)%ix = mg%boxes(id)%ix
1560 mg%boxes(c_id)%lvl = lvl
1561 mg%boxes(c_id)%parent = id
1564 mg%boxes(c_id)%neighbors = mg%boxes(id)%neighbors
1566 mg%boxes(c_id)%neighbors = mg%boxes(id)%neighbors + n_boxes_lvl
1568 mg%boxes(c_id)%r_min = mg%boxes(id)%r_min
1569 mg%boxes(c_id)%dr(:) = mg%dr(:, lvl)
1571 end subroutine add_single_child
1581 type(
mg_t),
intent(inout) :: mg
1582 integer :: i, id, lvl, single_cpu_lvl
1583 integer :: work_left, my_work, i_cpu
1587 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1589 do lvl = mg%lowest_lvl, single_cpu_lvl
1590 do i = 1,
size(mg%lvls(lvl)%ids)
1591 id = mg%lvls(lvl)%ids(i)
1592 mg%boxes(id)%rank = 0
1598 do lvl = single_cpu_lvl+1, mg%highest_lvl
1599 work_left =
size(mg%lvls(lvl)%ids)
1603 do i = 1,
size(mg%lvls(lvl)%ids)
1604 if ((mg%n_cpu - i_cpu - 1) * my_work >= work_left)
then
1609 my_work = my_work + 1
1610 work_left = work_left - 1
1612 id = mg%lvls(lvl)%ids(i)
1613 mg%boxes(id)%rank = i_cpu
1617 do lvl = mg%lowest_lvl, mg%highest_lvl
1618 call update_lvl_info(mg, mg%lvls(lvl))
1630 type(
mg_t),
intent(inout) :: mg
1631 integer :: i, id, lvl, single_cpu_lvl
1632 integer :: work_left, my_work(0:mg%n_cpu), i_cpu
1635 integer :: coarse_rank
1639 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1643 do lvl = mg%highest_lvl, single_cpu_lvl+1, -1
1647 do i = 1,
size(mg%lvls(lvl)%parents)
1648 id = mg%lvls(lvl)%parents(i)
1650 c_ids = mg%boxes(id)%children
1651 c_ranks = mg%boxes(c_ids)%rank
1652 i_cpu = most_popular(c_ranks, my_work, mg%n_cpu)
1653 mg%boxes(id)%rank = i_cpu
1654 my_work(i_cpu) = my_work(i_cpu) + 1
1657 work_left =
size(mg%lvls(lvl)%leaves)
1660 do i = 1,
size(mg%lvls(lvl)%leaves)
1662 if ((mg%n_cpu - i_cpu - 1) * my_work(i_cpu) >= &
1663 work_left + sum(my_work(i_cpu+1:)))
then
1667 my_work(i_cpu) = my_work(i_cpu) + 1
1668 work_left = work_left - 1
1670 id = mg%lvls(lvl)%leaves(i)
1671 mg%boxes(id)%rank = i_cpu
1676 if (single_cpu_lvl < mg%highest_lvl)
then
1677 coarse_rank = most_popular(mg%boxes(&
1678 mg%lvls(single_cpu_lvl+1)%ids)%rank, my_work, mg%n_cpu)
1683 do lvl = mg%lowest_lvl, single_cpu_lvl
1684 do i = 1,
size(mg%lvls(lvl)%ids)
1685 id = mg%lvls(lvl)%ids(i)
1686 mg%boxes(id)%rank = coarse_rank
1690 do lvl = mg%lowest_lvl, mg%highest_lvl
1691 call update_lvl_info(mg, mg%lvls(lvl))
1699 type(
mg_t),
intent(inout) :: mg
1700 integer :: i, id, lvl
1703 integer :: single_cpu_lvl, coarse_rank
1704 integer :: my_work(0:mg%n_cpu), i_cpu
1708 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1710 do lvl = mg%highest_lvl-1, single_cpu_lvl+1, -1
1714 do i = 1,
size(mg%lvls(lvl)%leaves)
1715 id = mg%lvls(lvl)%leaves(i)
1716 i_cpu = mg%boxes(id)%rank
1717 my_work(i_cpu) = my_work(i_cpu) + 1
1720 do i = 1,
size(mg%lvls(lvl)%parents)
1721 id = mg%lvls(lvl)%parents(i)
1723 c_ids = mg%boxes(id)%children
1724 c_ranks = mg%boxes(c_ids)%rank
1725 i_cpu = most_popular(c_ranks, my_work, mg%n_cpu)
1726 mg%boxes(id)%rank = i_cpu
1727 my_work(i_cpu) = my_work(i_cpu) + 1
1733 if (single_cpu_lvl < mg%highest_lvl)
then
1734 coarse_rank = most_popular(mg%boxes(&
1735 mg%lvls(single_cpu_lvl+1)%ids)%rank, my_work, mg%n_cpu)
1740 do lvl = mg%lowest_lvl, single_cpu_lvl
1741 do i = 1,
size(mg%lvls(lvl)%ids)
1742 id = mg%lvls(lvl)%ids(i)
1743 mg%boxes(id)%rank = coarse_rank
1747 do lvl = mg%lowest_lvl, mg%highest_lvl
1748 call update_lvl_info(mg, mg%lvls(lvl))
1755 pure integer function most_popular(list, work, n_cpu)
1756 integer,
intent(in) :: list(:)
1757 integer,
intent(in) :: n_cpu
1758 integer,
intent(in) :: work(0:n_cpu-1)
1759 integer :: i, best_count, current_count
1760 integer :: my_work, best_work
1766 do i = 1,
size(list)
1767 current_count = count(list == list(i))
1768 my_work = work(list(i))
1771 if (current_count > best_count .or. &
1772 current_count == best_count .and. my_work < best_work)
then
1773 best_count = current_count
1775 most_popular = list(i)
1779 end function most_popular
1781 subroutine update_lvl_info(mg, lvl)
1782 type(
mg_t),
intent(inout) :: mg
1783 type(
mg_lvl_t),
intent(inout) :: lvl
1785 lvl%my_ids = pack(lvl%ids, &
1786 mg%boxes(lvl%ids)%rank == mg%my_rank)
1787 lvl%my_leaves = pack(lvl%leaves, &
1788 mg%boxes(lvl%leaves)%rank == mg%my_rank)
1789 lvl%my_parents = pack(lvl%parents, &
1790 mg%boxes(lvl%parents)%rank == mg%my_rank)
1791 lvl%my_ref_bnds = pack(lvl%ref_bnds, &
1792 mg%boxes(lvl%ref_bnds)%rank == mg%my_rank)
1793 end subroutine update_lvl_info
1798 type(
mg_t),
intent(inout) :: mg
1800 if (mg%n_extra_vars == 0 .and. mg%is_allocated)
then
1801 error stop
"vlaplacian_set_methods: mg%n_extra_vars == 0"
1803 mg%n_extra_vars = max(1, mg%n_extra_vars)
1806 mg%subtract_mean = .false.
1811 mg%bc(:,
mg_iveps)%bc_value = 0.0_dp
1813 select case (mg%geometry_type)
1815 mg%box_op => box_vlpl
1820 select case (mg%smoother_type)
1822 mg%box_smoother => box_gs_vlpl
1824 error stop
"vlaplacian_set_methods: unsupported smoother type"
1828 error stop
"vlaplacian_set_methods: unsupported geometry"
1834 subroutine box_gs_vlpl(mg, id, nc, redblack_cntr)
1835 type(
mg_t),
intent(inout) :: mg
1836 integer,
intent(in) :: id
1837 integer,
intent(in) :: nc
1838 integer,
intent(in) :: redblack_cntr
1839 integer :: i, j, i0, di
1841 real(
dp) :: idr2(2*2), u(2*2)
1842 real(
dp) :: a0, a(2*2), c(2*2)
1845 idr2(1:2*2:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1846 idr2(2:2*2:2) = idr2(1:2*2:2)
1858 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
1862 i0 = 2 - iand(ieor(redblack_cntr, j), 1)
1865 a0 = cc(i, j, i_eps)
1866 u(1:2) = cc(i-1:i+1:2, j, n)
1867 a(1:2) = cc(i-1:i+1:2, j, i_eps)
1868 u(3:4) = cc(i, j-1:j+1:2, n)
1869 a(3:4) = cc(i, j-1:j+1:2, i_eps)
1870 c(:) = 2 * a0 * a(:) / (a0 + a(:)) * idr2
1873 (sum(c(:) * u(:)) - cc(i, j,
mg_irhs)) / sum(c(:))
1877 end subroutine box_gs_vlpl
1880 subroutine box_vlpl(mg, id, nc, i_out)
1881 type(
mg_t),
intent(inout) :: mg
1882 integer,
intent(in) :: id
1883 integer,
intent(in) :: nc
1884 integer,
intent(in) :: i_out
1886 real(
dp) :: idr2(2*2), a0, u0, u(2*2), a(2*2)
1889 idr2(1:2*2:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1890 idr2(2:2*2:2) = idr2(1:2*2:2)
1892 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
1896 a0 = cc(i, j, i_eps)
1897 a(1:2) = cc(i-1:i+1:2, j, i_eps)
1898 a(3:4) = cc(i, j-1:j+1:2, i_eps)
1900 u(1:2) = cc(i-1:i+1:2, j, n)
1901 u(3:4) = cc(i, j-1:j+1:2, n)
1903 cc(i, j, i_out) = sum(2 * idr2 * &
1904 a0*a(:)/(a0 + a(:)) * (u(:) - u0))
1908 end subroutine box_vlpl
2016 type(
mg_t),
intent(inout) :: mg
2018 integer,
intent(in),
optional :: comm
2020 logical :: initialized
2022 call mpi_initialized(initialized, ierr)
2023 if (.not. initialized)
then
2027 if (
present(comm))
then
2030 mg%comm = mpi_comm_world
2033 call mpi_comm_rank(mg%comm, mg%my_rank, ierr)
2034 call mpi_comm_size(mg%comm, mg%n_cpu, ierr)
2039 type(
mg_t),
intent(inout) :: mg
2040 integer,
intent(in) :: dsize
2041 integer :: i, n_send, n_recv
2042 integer :: send_req(mg%n_cpu)
2043 integer :: recv_req(mg%n_cpu)
2049 do i = 0, mg%n_cpu - 1
2050 if (mg%buf(i)%i_send > 0)
then
2053 call mpi_isend(mg%buf(i)%send, mg%buf(i)%i_send, mpi_double, &
2054 i, 0, mg%comm, send_req(n_send), ierr)
2056 if (mg%buf(i)%i_recv > 0)
then
2058 call mpi_irecv(mg%buf(i)%recv, mg%buf(i)%i_recv, mpi_double, &
2059 i, 0, mg%comm, recv_req(n_recv), ierr)
2063 call mpi_waitall(n_recv, recv_req(1:n_recv), mpi_statuses_ignore, ierr)
2064 call mpi_waitall(n_send, send_req(1:n_send), mpi_statuses_ignore, ierr)
2071 type(
mg_buf_t),
intent(inout) :: gc
2072 integer,
intent(in) :: dsize
2073 integer :: ix_sort(gc%i_ix)
2074 real(dp) :: buf_cpy(gc%i_send)
2077 call mrgrnk(gc%ix(1:gc%i_ix), ix_sort)
2079 buf_cpy = gc%send(1:gc%i_send)
2083 j = (ix_sort(n)-1) * dsize
2084 gc%send(i+1:i+dsize) = buf_cpy(j+1:j+dsize)
2086 gc%ix(1:gc%i_ix) = gc%ix(ix_sort)
2094 type(
mg_t),
intent(inout) :: mg
2095 integer,
intent(out) :: n_send(0:mg%n_cpu-1)
2096 integer,
intent(out) :: n_recv(0:mg%n_cpu-1)
2097 integer,
intent(out) :: dsize
2098 integer :: i, id, lvl, nc
2100 allocate(mg%comm_ghostcell%n_send(0:mg%n_cpu-1, &
2101 mg%first_normal_lvl:mg%highest_lvl))
2102 allocate(mg%comm_ghostcell%n_recv(0:mg%n_cpu-1, &
2103 mg%first_normal_lvl:mg%highest_lvl))
2105 dsize = mg%box_size**(2-1)
2107 do lvl = mg%first_normal_lvl, mg%highest_lvl
2108 nc = mg%box_size_lvl(lvl)
2109 mg%buf(:)%i_send = 0
2110 mg%buf(:)%i_recv = 0
2113 do i = 1,
size(mg%lvls(lvl)%my_ids)
2114 id = mg%lvls(lvl)%my_ids(i)
2115 call buffer_ghost_cells(mg, id, nc, 1, dry_run=.true.)
2119 do i = 1,
size(mg%lvls(lvl-1)%my_ref_bnds)
2120 id = mg%lvls(lvl-1)%my_ref_bnds(i)
2121 call buffer_refinement_boundaries(mg, id, nc, 1, dry_run=.true.)
2126 mg%buf(:)%i_recv = 0
2127 do i = 1,
size(mg%lvls(lvl)%my_ids)
2128 id = mg%lvls(lvl)%my_ids(i)
2129 call set_ghost_cells(mg, id, nc, 1, dry_run=.true.)
2132 mg%comm_ghostcell%n_send(:, lvl) = mg%buf(:)%i_send/dsize
2133 mg%comm_ghostcell%n_recv(:, lvl) = mg%buf(:)%i_recv/dsize
2136 n_send = maxval(mg%comm_ghostcell%n_send, dim=2)
2137 n_recv = maxval(mg%comm_ghostcell%n_recv, dim=2)
2143 type(
mg_t),
intent(inout) :: mg
2148 do lvl = mg%lowest_lvl, mg%highest_lvl
2149 nc = mg%box_size_lvl(lvl)
2150 call mg_phi_bc_store_lvl(mg, lvl, nc)
2153 mg%phi_bc_data_stored = .true.
2156 subroutine mg_phi_bc_store_lvl(mg, lvl, nc)
2157 type(
mg_t),
intent(inout) :: mg
2158 integer,
intent(in) :: lvl
2159 integer,
intent(in) :: nc
2161 integer :: i, id, nb, nb_id, bc_type
2163 do i = 1,
size(mg%lvls(lvl)%my_ids)
2164 id = mg%lvls(lvl)%my_ids(i)
2166 nb_id = mg%boxes(id)%neighbors(nb)
2169 if (
associated(mg%bc(nb,
mg_iphi)%boundary_cond))
then
2170 call mg%bc(nb,
mg_iphi)%boundary_cond(mg%boxes(id), nc, &
2173 bc_type = mg%bc(nb,
mg_iphi)%bc_type
2174 bc = mg%bc(nb,
mg_iphi)%bc_value
2180 mg%boxes(id)%neighbors(nb) = bc_type
2183 call box_set_gc(mg%boxes(id), nb, nc,
mg_irhs, bc)
2187 end subroutine mg_phi_bc_store_lvl
2192 integer,
intent(in) :: iv
2195 do lvl = mg%lowest_lvl, mg%highest_lvl
2204 integer,
intent(in) :: lvl
2205 integer,
intent(in) :: iv
2206 integer :: i, id, dsize, nc
2208 if (lvl < mg%lowest_lvl) &
2209 error stop
"fill_ghost_cells_lvl: lvl < lowest_lvl"
2210 if (lvl > mg%highest_lvl) &
2211 error stop
"fill_ghost_cells_lvl: lvl > highest_lvl"
2213 nc = mg%box_size_lvl(lvl)
2215 if (lvl >= mg%first_normal_lvl)
then
2217 mg%buf(:)%i_send = 0
2218 mg%buf(:)%i_recv = 0
2221 do i = 1,
size(mg%lvls(lvl)%my_ids)
2222 id = mg%lvls(lvl)%my_ids(i)
2223 call buffer_ghost_cells(mg, id, nc, iv, .false.)
2227 do i = 1,
size(mg%lvls(lvl-1)%my_ref_bnds)
2228 id = mg%lvls(lvl-1)%my_ref_bnds(i)
2229 call buffer_refinement_boundaries(mg, id, nc, iv, .false.)
2234 mg%buf(:)%i_recv = mg%comm_ghostcell%n_recv(:, lvl) * dsize
2238 mg%buf(:)%i_recv = 0
2241 do i = 1,
size(mg%lvls(lvl)%my_ids)
2242 id = mg%lvls(lvl)%my_ids(i)
2243 call set_ghost_cells(mg, id, nc, iv, .false.)
2247 subroutine buffer_ghost_cells(mg, id, nc, iv, dry_run)
2248 type(
mg_t),
intent(inout) :: mg
2249 integer,
intent(in) :: id
2250 integer,
intent(in) :: nc
2251 integer,
intent(in) :: iv
2252 logical,
intent(in) :: dry_run
2253 integer :: nb, nb_id, nb_rank
2256 nb_id = mg%boxes(id)%neighbors(nb)
2260 nb_rank = mg%boxes(nb_id)%rank
2262 if (nb_rank /= mg%my_rank)
then
2263 call buffer_for_nb(mg, mg%boxes(id), nc, iv, nb_id, nb_rank, &
2268 end subroutine buffer_ghost_cells
2270 subroutine buffer_refinement_boundaries(mg, id, nc, iv, dry_run)
2271 type(
mg_t),
intent(inout) :: mg
2272 integer,
intent(in) :: id
2273 integer,
intent(in) :: nc
2274 integer,
intent(in) :: iv
2275 logical,
intent(in) :: dry_run
2276 integer :: nb, nb_id, c_ids(2**(2-1))
2277 integer :: n, c_id, c_rank
2280 nb_id = mg%boxes(id)%neighbors(nb)
2283 c_ids = mg%boxes(nb_id)%children(&
2288 c_rank = mg%boxes(c_id)%rank
2290 if (c_rank /= mg%my_rank)
then
2292 call buffer_for_fine_nb(mg, mg%boxes(id), nc, iv, c_id, &
2293 c_rank, nb, dry_run)
2299 end subroutine buffer_refinement_boundaries
2302 subroutine set_ghost_cells(mg, id, nc, iv, dry_run)
2303 type(
mg_t),
intent(inout) :: mg
2304 integer,
intent(in) :: id
2305 integer,
intent(in) :: nc
2306 integer,
intent(in) :: iv
2307 logical,
intent(in) :: dry_run
2309 integer :: nb, nb_id, nb_rank, bc_type
2312 nb_id = mg%boxes(id)%neighbors(nb)
2316 nb_rank = mg%boxes(nb_id)%rank
2318 if (nb_rank /= mg%my_rank)
then
2319 call fill_buffered_nb(mg, mg%boxes(id), nb_rank, &
2320 nb, nc, iv, dry_run)
2321 else if (.not. dry_run)
then
2322 call copy_from_nb(mg%boxes(id), mg%boxes(nb_id), &
2327 call fill_refinement_bnd(mg, id, nb, nc, iv, dry_run)
2328 else if (.not. dry_run)
then
2330 if (mg%phi_bc_data_stored .and. iv ==
mg_iphi)
then
2333 call box_get_gc(mg%boxes(id), nb, nc,
mg_irhs, bc)
2336 if (
associated(mg%bc(nb, iv)%boundary_cond))
then
2337 call mg%bc(nb, iv)%boundary_cond(mg%boxes(id), nc, iv, &
2340 bc_type = mg%bc(nb, iv)%bc_type
2341 bc = mg%bc(nb, iv)%bc_value
2345 call box_set_gc(mg%boxes(id), nb, nc, iv, bc)
2346 call bc_to_gc(mg, id, nc, iv, nb, bc_type)
2349 end subroutine set_ghost_cells
2351 subroutine fill_refinement_bnd(mg, id, nb, nc, iv, dry_run)
2352 type(
mg_t),
intent(inout) :: mg
2353 integer,
intent(in) :: id
2354 integer,
intent(in) :: nc
2355 integer,
intent(in) :: iv
2356 integer,
intent(in) :: nb
2357 logical,
intent(in) :: dry_run
2359 integer :: p_id, p_nb_id, ix_offset(2)
2360 integer :: i, dsize, p_nb_rank
2363 p_id = mg%boxes(id)%parent
2364 p_nb_id = mg%boxes(p_id)%neighbors(nb)
2365 p_nb_rank = mg%boxes(p_nb_id)%rank
2367 if (p_nb_rank /= mg%my_rank)
then
2368 i = mg%buf(p_nb_rank)%i_recv
2369 if (.not. dry_run)
then
2370 gc = reshape(mg%buf(p_nb_rank)%recv(i+1:i+dsize), shape(gc))
2372 mg%buf(p_nb_rank)%i_recv = mg%buf(p_nb_rank)%i_recv + dsize
2373 else if (.not. dry_run)
then
2375 call box_gc_for_fine_neighbor(mg%boxes(p_nb_id),
mg_neighb_rev(nb), &
2376 ix_offset, nc, iv, gc)
2379 if (.not. dry_run)
then
2380 if (
associated(mg%bc(nb, iv)%refinement_bnd))
then
2381 call mg%bc(nb, iv)%refinement_bnd(mg%boxes(id), nc, iv, nb, gc)
2383 call sides_rb(mg%boxes(id), nc, iv, nb, gc)
2386 end subroutine fill_refinement_bnd
2388 subroutine copy_from_nb(box, box_nb, nb, nc, iv)
2389 type(
mg_box_t),
intent(inout) :: box
2390 type(
mg_box_t),
intent(in) :: box_nb
2391 integer,
intent(in) :: nb
2392 integer,
intent(in) :: nc
2393 integer,
intent(in) :: iv
2396 call box_gc_for_neighbor(box_nb,
mg_neighb_rev(nb), nc, iv, gc)
2397 call box_set_gc(box, nb, nc, iv, gc)
2398 end subroutine copy_from_nb
2400 subroutine buffer_for_nb(mg, box, nc, iv, nb_id, nb_rank, nb, dry_run)
2401 type(
mg_t),
intent(inout) :: mg
2402 type(
mg_box_t),
intent(inout) :: box
2403 integer,
intent(in) :: nc
2404 integer,
intent(in) :: iv
2405 integer,
intent(in) :: nb_id
2406 integer,
intent(in) :: nb_rank
2407 integer,
intent(in) :: nb
2408 logical,
intent(in) :: dry_run
2412 i = mg%buf(nb_rank)%i_send
2415 if (.not. dry_run)
then
2416 call box_gc_for_neighbor(box, nb, nc, iv, gc)
2417 mg%buf(nb_rank)%send(i+1:i+dsize) = pack(gc, .true.)
2422 i = mg%buf(nb_rank)%i_ix
2423 if (.not. dry_run)
then
2427 mg%buf(nb_rank)%i_send = mg%buf(nb_rank)%i_send + dsize
2428 mg%buf(nb_rank)%i_ix = mg%buf(nb_rank)%i_ix + 1
2429 end subroutine buffer_for_nb
2431 subroutine buffer_for_fine_nb(mg, box, nc, iv, fine_id, fine_rank, nb, dry_run)
2432 type(
mg_t),
intent(inout) :: mg
2433 type(
mg_box_t),
intent(inout) :: box
2434 integer,
intent(in) :: nc
2435 integer,
intent(in) :: iv
2436 integer,
intent(in) :: fine_id
2437 integer,
intent(in) :: fine_rank
2438 integer,
intent(in) :: nb
2439 logical,
intent(in) :: dry_run
2440 integer :: i, dsize, ix_offset(2)
2443 i = mg%buf(fine_rank)%i_send
2446 if (.not. dry_run)
then
2448 call box_gc_for_fine_neighbor(box, nb, ix_offset, nc, iv, gc)
2449 mg%buf(fine_rank)%send(i+1:i+dsize) = pack(gc, .true.)
2454 i = mg%buf(fine_rank)%i_ix
2455 if (.not. dry_run)
then
2460 mg%buf(fine_rank)%i_send = mg%buf(fine_rank)%i_send + dsize
2461 mg%buf(fine_rank)%i_ix = mg%buf(fine_rank)%i_ix + 1
2462 end subroutine buffer_for_fine_nb
2464 subroutine fill_buffered_nb(mg, box, nb_rank, nb, nc, iv, dry_run)
2465 type(
mg_t),
intent(inout) :: mg
2466 type(
mg_box_t),
intent(inout) :: box
2467 integer,
intent(in) :: nb_rank
2468 integer,
intent(in) :: nb
2469 integer,
intent(in) :: nc
2470 integer,
intent(in) :: iv
2471 logical,
intent(in) :: dry_run
2475 i = mg%buf(nb_rank)%i_recv
2478 if (.not. dry_run)
then
2479 gc = reshape(mg%buf(nb_rank)%recv(i+1:i+dsize), shape(gc))
2480 call box_set_gc(box, nb, nc, iv, gc)
2482 mg%buf(nb_rank)%i_recv = mg%buf(nb_rank)%i_recv + dsize
2484 end subroutine fill_buffered_nb
2486 subroutine box_gc_for_neighbor(box, nb, nc, iv, gc)
2488 integer,
intent(in) :: nb, nc, iv
2489 real(
dp),
intent(out) :: gc(nc)
2493 gc = box%cc(1, 1:nc, iv)
2495 gc = box%cc(nc, 1:nc, iv)
2497 gc = box%cc(1:nc, 1, iv)
2499 gc = box%cc(1:nc, nc, iv)
2501 end subroutine box_gc_for_neighbor
2504 subroutine box_gc_for_fine_neighbor(box, nb, di, nc, iv, gc)
2506 integer,
intent(in) :: nb
2507 integer,
intent(in) :: di(2)
2508 integer,
intent(in) :: nc, iv
2509 real(
dp),
intent(out) :: gc(nc)
2510 real(
dp) :: tmp(0:nc/2+1), grad(2-1)
2518 tmp = box%cc(1, di(2):di(2)+hnc+1, iv)
2520 tmp = box%cc(nc, di(2):di(2)+hnc+1, iv)
2522 tmp = box%cc(di(1):di(1)+hnc+1, 1, iv)
2524 tmp = box%cc(di(1):di(1)+hnc+1, nc, iv)
2532 grad(1) = 0.125_dp * (tmp(i+1) - tmp(i-1))
2533 gc(2*i-1) = tmp(i) - grad(1)
2534 gc(2*i) = tmp(i) + grad(1)
2536 end subroutine box_gc_for_fine_neighbor
2538 subroutine box_get_gc(box, nb, nc, iv, gc)
2540 integer,
intent(in) :: nb, nc, iv
2541 real(
dp),
intent(out) :: gc(nc)
2545 gc = box%cc(0, 1:nc, iv)
2547 gc = box%cc(nc+1, 1:nc, iv)
2549 gc = box%cc(1:nc, 0, iv)
2551 gc = box%cc(1:nc, nc+1, iv)
2553 end subroutine box_get_gc
2555 subroutine box_set_gc(box, nb, nc, iv, gc)
2556 type(
mg_box_t),
intent(inout) :: box
2557 integer,
intent(in) :: nb, nc, iv
2558 real(
dp),
intent(in) :: gc(nc)
2562 box%cc(0, 1:nc, iv) = gc
2564 box%cc(nc+1, 1:nc, iv) = gc
2566 box%cc(1:nc, 0, iv) = gc
2568 box%cc(1:nc, nc+1, iv) = gc
2570 end subroutine box_set_gc
2572 subroutine bc_to_gc(mg, id, nc, iv, nb, bc_type)
2573 type(
mg_t),
intent(inout) :: mg
2574 integer,
intent(in) :: id
2575 integer,
intent(in) :: nc
2576 integer,
intent(in) :: iv
2577 integer,
intent(in) :: nb
2578 integer,
intent(in) :: bc_type
2579 real(
dp) :: c0, c1, c2, dr
2589 select case (bc_type)
2604 error stop
"bc_to_gc: unknown boundary condition"
2609 mg%boxes(id)%cc(0, 1:nc, iv) = &
2610 c0 * mg%boxes(id)%cc(0, 1:nc, iv) + &
2611 c1 * mg%boxes(id)%cc(1, 1:nc, iv) + &
2612 c2 * mg%boxes(id)%cc(2, 1:nc, iv)
2614 mg%boxes(id)%cc(nc+1, 1:nc, iv) = &
2615 c0 * mg%boxes(id)%cc(nc+1, 1:nc, iv) + &
2616 c1 * mg%boxes(id)%cc(nc, 1:nc, iv) + &
2617 c2 * mg%boxes(id)%cc(nc-1, 1:nc, iv)
2619 mg%boxes(id)%cc(1:nc, 0, iv) = &
2620 c0 * mg%boxes(id)%cc(1:nc, 0, iv) + &
2621 c1 * mg%boxes(id)%cc(1:nc, 1, iv) + &
2622 c2 * mg%boxes(id)%cc(1:nc, 2, iv)
2624 mg%boxes(id)%cc(1:nc, nc+1, iv) = &
2625 c0 * mg%boxes(id)%cc(1:nc, nc+1, iv) + &
2626 c1 * mg%boxes(id)%cc(1:nc, nc, iv) + &
2627 c2 * mg%boxes(id)%cc(1:nc, nc-1, iv)
2629 end subroutine bc_to_gc
2632 subroutine sides_rb(box, nc, iv, nb, gc)
2633 type(
mg_box_t),
intent(inout) :: box
2634 integer,
intent(in) :: nc
2635 integer,
intent(in) :: iv
2636 integer,
intent(in) :: nb
2638 real(
dp),
intent(in) :: gc(nc)
2640 integer :: i, j, ix, dix
2656 dj = -1 + 2 * iand(j, 1)
2657 box%cc(i-di, j, iv) = 0.5_dp * gc(j) &
2658 + 0.75_dp * box%cc(i, j, iv) &
2659 - 0.25_dp * box%cc(i+di, j, iv)
2665 di = -1 + 2 * iand(i, 1)
2666 box%cc(i, j-dj, iv) = 0.5_dp * gc(i) &
2667 + 0.75_dp * box%cc(i, j, iv) &
2668 - 0.25_dp * box%cc(i, j+dj, iv)
2672 end subroutine sides_rb
2678 type(
mg_t),
intent(inout) :: mg
2679 integer,
intent(out) :: n_send(0:mg%n_cpu-1)
2680 integer,
intent(out) :: n_recv(0:mg%n_cpu-1)
2681 integer,
intent(out) :: dsize
2682 integer :: lvl, min_lvl
2684 if (.not.
allocated(mg%comm_restrict%n_send))
then
2685 error stop
"Call restrict_buffer_size before prolong_buffer_size"
2688 min_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
2689 allocate(mg%comm_prolong%n_send(0:mg%n_cpu-1, &
2690 min_lvl:mg%highest_lvl))
2691 allocate(mg%comm_prolong%n_recv(0:mg%n_cpu-1, &
2692 min_lvl:mg%highest_lvl))
2694 mg%comm_prolong%n_recv(:, mg%highest_lvl) = 0
2695 mg%comm_prolong%n_send(:, mg%highest_lvl) = 0
2697 do lvl = min_lvl, mg%highest_lvl-1
2698 mg%comm_prolong%n_recv(:, lvl) = &
2699 mg%comm_restrict%n_send(:, lvl+1)
2700 mg%comm_prolong%n_send(:, lvl) = &
2701 mg%comm_restrict%n_recv(:, lvl+1)
2706 dsize = (mg%box_size)**2
2707 n_send = maxval(mg%comm_prolong%n_send, dim=2)
2708 n_recv = maxval(mg%comm_prolong%n_recv, dim=2)
2714 type(
mg_t),
intent(inout) :: mg
2715 integer,
intent(in) :: lvl
2716 integer,
intent(in) :: iv
2717 integer,
intent(in) :: iv_to
2718 procedure(mg_box_prolong) :: method
2719 logical,
intent(in) :: add
2720 integer :: i, id, dsize, nc
2722 if (lvl == mg%highest_lvl) error stop
"cannot prolong highest level"
2723 if (lvl < mg%lowest_lvl) error stop
"cannot prolong below lowest level"
2726 if (lvl >= mg%first_normal_lvl-1)
then
2727 dsize = mg%box_size**2
2728 mg%buf(:)%i_send = 0
2731 do i = 1,
size(mg%lvls(lvl)%my_ids)
2732 id = mg%lvls(lvl)%my_ids(i)
2733 call prolong_set_buffer(mg, id, mg%box_size, iv, method)
2736 mg%buf(:)%i_recv = mg%comm_prolong%n_recv(:, lvl) * dsize
2738 mg%buf(:)%i_recv = 0
2741 nc = mg%box_size_lvl(lvl+1)
2742 do i = 1,
size(mg%lvls(lvl+1)%my_ids)
2743 id = mg%lvls(lvl+1)%my_ids(i)
2744 call prolong_onto(mg, id, nc, iv, iv_to, add, method)
2751 subroutine prolong_set_buffer(mg, id, nc, iv, method)
2752 type(
mg_t),
intent(inout) :: mg
2753 integer,
intent(in) :: id
2754 integer,
intent(in) :: nc
2755 integer,
intent(in) :: iv
2756 procedure(mg_box_prolong) :: method
2757 integer :: i, dix(2)
2758 integer :: i_c, c_id, c_rank, dsize
2759 real(
dp) :: tmp(nc, nc)
2764 c_id = mg%boxes(id)%children(i_c)
2766 c_rank = mg%boxes(c_id)%rank
2767 if (c_rank /= mg%my_rank)
then
2769 call method(mg, id, dix, nc, iv, tmp)
2771 i = mg%buf(c_rank)%i_send
2772 mg%buf(c_rank)%send(i+1:i+dsize) = pack(tmp, .true.)
2773 mg%buf(c_rank)%i_send = mg%buf(c_rank)%i_send + dsize
2775 i = mg%buf(c_rank)%i_ix
2776 mg%buf(c_rank)%ix(i+1) = c_id
2777 mg%buf(c_rank)%i_ix = mg%buf(c_rank)%i_ix + 1
2781 end subroutine prolong_set_buffer
2784 subroutine prolong_onto(mg, id, nc, iv, iv_to, add, method)
2785 type(
mg_t),
intent(inout) :: mg
2786 integer,
intent(in) :: id
2787 integer,
intent(in) :: nc
2788 integer,
intent(in) :: iv
2789 integer,
intent(in) :: iv_to
2790 logical,
intent(in) :: add
2791 procedure(mg_box_prolong) :: method
2792 integer :: hnc, p_id, p_rank, i, dix(2), dsize
2793 real(
dp) :: tmp(nc, nc)
2796 p_id = mg%boxes(id)%parent
2797 p_rank = mg%boxes(p_id)%rank
2799 if (p_rank == mg%my_rank)
then
2801 call method(mg, p_id, dix, nc, iv, tmp)
2804 i = mg%buf(p_rank)%i_recv
2805 tmp = reshape(mg%buf(p_rank)%recv(i+1:i+dsize), [nc, nc])
2806 mg%buf(p_rank)%i_recv = mg%buf(p_rank)%i_recv + dsize
2810 mg%boxes(id)%cc(1:nc, 1:nc, iv_to) = &
2811 mg%boxes(id)%cc(1:nc, 1:nc, iv_to) + tmp
2813 mg%boxes(id)%cc(1:nc, 1:nc, iv_to) = tmp
2816 end subroutine prolong_onto
2820 type(
mg_t),
intent(inout) :: mg
2821 integer,
intent(in) :: p_id
2822 integer,
intent(in) :: dix(2)
2823 integer,
intent(in) :: nc
2824 integer,
intent(in) :: iv
2825 real(
dp),
intent(out) :: fine(nc, nc)
2827 integer :: i, j, hnc
2829 real(
dp) :: f0, flx, fhx, fly, fhy
2833 associate(crs => mg%boxes(p_id)%cc)
2839 f0 = 0.5_dp * crs(ic, jc, iv)
2840 flx = 0.25_dp * crs(ic-1, jc, iv)
2841 fhx = 0.25_dp * crs(ic+1, jc, iv)
2842 fly = 0.25_dp * crs(ic, jc-1, iv)
2843 fhy = 0.25_dp * crs(ic, jc+1, iv)
2845 fine(2*i-1, 2*j-1) = f0 + flx + fly
2846 fine(2*i , 2*j-1) = f0 + fhx + fly
2847 fine(2*i-1, 2*j) = f0 + flx + fhy
2848 fine(2*i , 2*j) = f0 + fhx + fhy
2857 type(
mg_t),
intent(inout) :: mg
2859 mg%subtract_mean = .false.
2861 select case (mg%geometry_type)
2863 mg%box_op => box_helmh
2865 select case (mg%smoother_type)
2867 mg%box_smoother => box_gs_helmh
2869 error stop
"helmholtz_set_methods: unsupported smoother type"
2872 error stop
"helmholtz_set_methods: unsupported geometry"
2878 real(
dp),
intent(in) :: lambda
2881 error stop
"helmholtz_set_lambda: lambda < 0 not allowed"
2887 subroutine box_gs_helmh(mg, id, nc, redblack_cntr)
2888 type(
mg_t),
intent(inout) :: mg
2889 integer,
intent(in) :: id
2890 integer,
intent(in) :: nc
2891 integer,
intent(in) :: redblack_cntr
2892 integer :: i, j, i0, di
2893 real(
dp) :: idr2(2), fac
2896 idr2 = 1/mg%dr(:, mg%boxes(id)%lvl)**2
2909 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
2912 i0 = 2 - iand(ieor(redblack_cntr, j), 1)
2915 cc(i, j, n) = fac * ( &
2916 idr2(1) * (cc(i+1, j, n) + cc(i-1, j, n)) + &
2917 idr2(2) * (cc(i, j+1, n) + cc(i, j-1, n)) - &
2922 end subroutine box_gs_helmh
2925 subroutine box_helmh(mg, id, nc, i_out)
2926 type(
mg_t),
intent(inout) :: mg
2927 integer,
intent(in) :: id
2928 integer,
intent(in) :: nc
2929 integer,
intent(in) :: i_out
2933 idr2 = 1 / mg%dr(:, mg%boxes(id)%lvl)**2
2935 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
2939 idr2(1) * (cc(i-1, j, n) + cc(i+1, j, n) - 2 * cc(i, j, n)) + &
2940 idr2(2) * (cc(i, j-1, n) + cc(i, j+1, n) - 2 * cc(i, j, n)) - &
2945 end subroutine box_helmh
2951 type(
mg_t),
intent(inout) :: mg
2956 select case (mg%operator_type)
2968 error stop
"mg_set_methods: unknown operator"
2974 mg%n_smoother_substeps = 2
2976 mg%n_smoother_substeps = 1
2980 subroutine check_methods(mg)
2981 type(
mg_t),
intent(inout) :: mg
2983 if (.not.
associated(mg%box_op) .or. &
2984 .not.
associated(mg%box_smoother))
then
2988 end subroutine check_methods
2990 subroutine mg_add_timers(mg)
2991 type(
mg_t),
intent(inout) :: mg
2992 timer_total_vcycle =
mg_add_timer(mg,
"mg total V-cycle")
2993 timer_total_fmg =
mg_add_timer(mg,
"mg total FMG cycle")
2995 timer_smoother_gc =
mg_add_timer(mg,
"mg smoother g.c.")
2998 timer_update_coarse =
mg_add_timer(mg,
"mg update coarse")
2999 end subroutine mg_add_timers
3003 type(
mg_t),
intent(inout) :: mg
3004 logical,
intent(in) :: have_guess
3005 real(
dp),
intent(out),
optional :: max_res
3006 integer :: lvl, i, id
3008 call check_methods(mg)
3009 if (timer_smoother == -1)
call mg_add_timers(mg)
3013 if (.not. have_guess)
then
3014 do lvl = mg%highest_lvl, mg%lowest_lvl, -1
3015 do i = 1,
size(mg%lvls(lvl)%my_ids)
3016 id = mg%lvls(lvl)%my_ids(i)
3017 mg%boxes(id)%cc(:, :,
mg_iphi) = 0.0_dp
3025 do lvl = mg%highest_lvl, mg%lowest_lvl+1, -1
3028 call update_coarse(mg, lvl)
3032 if (mg%subtract_mean)
then
3037 do lvl = mg%lowest_lvl, mg%highest_lvl
3039 do i = 1,
size(mg%lvls(lvl)%my_ids)
3040 id = mg%lvls(lvl)%my_ids(i)
3041 mg%boxes(id)%cc(:, :,
mg_iold) = &
3042 mg%boxes(id)%cc(:, :,
mg_iphi)
3045 if (lvl > mg%lowest_lvl)
then
3049 call correct_children(mg, lvl-1)
3057 if (lvl == mg%highest_lvl)
then
3070 type(
mg_t),
intent(inout) :: mg
3071 integer,
intent(in),
optional :: highest_lvl
3072 real(
dp),
intent(out),
optional :: max_res
3074 logical,
intent(in),
optional :: standalone
3075 integer :: lvl, min_lvl, i, max_lvl, ierr
3076 real(
dp) :: init_res, res
3077 logical :: is_standalone
3079 is_standalone = .true.
3080 if (
present(standalone)) is_standalone = standalone
3082 call check_methods(mg)
3083 if (timer_smoother == -1)
call mg_add_timers(mg)
3085 if (is_standalone) &
3088 if (mg%subtract_mean .and. .not.
present(highest_lvl))
then
3094 min_lvl = mg%lowest_lvl
3095 max_lvl = mg%highest_lvl
3096 if (
present(highest_lvl)) max_lvl = highest_lvl
3099 if (is_standalone)
then
3103 do lvl = max_lvl, min_lvl+1, -1
3105 call smooth_boxes(mg, lvl, mg%n_cycle_down)
3110 call update_coarse(mg, lvl)
3115 if (.not. all(mg%boxes(mg%lvls(min_lvl)%ids)%rank == &
3116 mg%boxes(mg%lvls(min_lvl)%ids(1))%rank))
then
3117 error stop
"Multiple CPUs for coarse grid (not implemented yet)"
3120 init_res = max_residual_lvl(mg, min_lvl)
3121 do i = 1, mg%max_coarse_cycles
3122 call smooth_boxes(mg, min_lvl, mg%n_cycle_up+mg%n_cycle_down)
3123 res = max_residual_lvl(mg, min_lvl)
3124 if (res < mg%residual_coarse_rel * init_res .or. &
3125 res < mg%residual_coarse_abs)
exit
3130 do lvl = min_lvl+1, max_lvl
3134 call correct_children(mg, lvl-1)
3141 call smooth_boxes(mg, lvl, mg%n_cycle_up)
3144 if (
present(max_res))
then
3146 do lvl = min_lvl, max_lvl
3147 res = max_residual_lvl(mg, lvl)
3148 init_res = max(res, init_res)
3150 call mpi_allreduce(init_res, max_res, 1, &
3151 mpi_double, mpi_max, mg%comm, ierr)
3155 if (mg%subtract_mean)
then
3159 if (is_standalone) &
3165 type(
mg_t),
intent(inout) :: mg
3166 integer,
intent(in) :: iv
3167 logical,
intent(in) :: include_ghostcells
3168 integer :: i, id, lvl, nc, ierr
3169 real(dp) :: sum_iv, mean_iv, volume
3173 call mpi_allreduce(sum_iv, mean_iv, 1, &
3174 mpi_double, mpi_sum, mg%comm, ierr)
3177 volume = nc**2 * product(mg%dr(:, 1)) *
size(mg%lvls(1)%ids)
3178 mean_iv = mean_iv / volume
3180 do lvl = mg%lowest_lvl, mg%highest_lvl
3181 nc = mg%box_size_lvl(lvl)
3183 do i = 1,
size(mg%lvls(lvl)%my_ids)
3184 id = mg%lvls(lvl)%my_ids(i)
3185 if (include_ghostcells)
then
3186 mg%boxes(id)%cc(:, :, iv) = &
3187 mg%boxes(id)%cc(:, :, iv) - mean_iv
3189 mg%boxes(id)%cc(1:nc, 1:nc, iv) = &
3190 mg%boxes(id)%cc(1:nc, 1:nc, iv) - mean_iv
3197 type(
mg_t),
intent(in) :: mg
3198 integer,
intent(in) :: iv
3199 integer :: lvl, i, id, nc
3203 do lvl = 1, mg%highest_lvl
3204 nc = mg%box_size_lvl(lvl)
3205 w = product(mg%dr(:, lvl))
3206 do i = 1,
size(mg%lvls(lvl)%my_leaves)
3207 id = mg%lvls(lvl)%my_leaves(i)
3209 sum(mg%boxes(id)%cc(1:nc, 1:nc, iv))
3214 real(
dp) function max_residual_lvl(mg, lvl)
3215 type(
mg_t),
intent(inout) :: mg
3216 integer,
intent(in) :: lvl
3217 integer :: i, id, nc
3220 nc = mg%box_size_lvl(lvl)
3221 max_residual_lvl = 0.0_dp
3223 do i = 1,
size(mg%lvls(lvl)%my_ids)
3224 id = mg%lvls(lvl)%my_ids(i)
3225 call residual_box(mg, id, nc)
3226 res = maxval(abs(mg%boxes(id)%cc(1:nc, 1:nc,
mg_ires)))
3227 max_residual_lvl = max(max_residual_lvl, res)
3229 end function max_residual_lvl
3265 subroutine update_coarse(mg, lvl)
3266 type(
mg_t),
intent(inout) :: mg
3267 integer,
intent(in) :: lvl
3268 integer :: i, id, nc, nc_c
3270 nc = mg%box_size_lvl(lvl)
3271 nc_c = mg%box_size_lvl(lvl-1)
3274 do i = 1,
size(mg%lvls(lvl)%my_ids)
3275 id = mg%lvls(lvl)%my_ids(i)
3276 call residual_box(mg, id, nc)
3287 do i = 1,
size(mg%lvls(lvl-1)%my_parents)
3288 id = mg%lvls(lvl-1)%my_parents(i)
3291 call mg%box_op(mg, id, nc_c,
mg_irhs)
3294 mg%boxes(id)%cc(1:nc_c, 1:nc_c,
mg_irhs) = &
3295 mg%boxes(id)%cc(1:nc_c, 1:nc_c,
mg_irhs) + &
3296 mg%boxes(id)%cc(1:nc_c, 1:nc_c,
mg_ires)
3299 mg%boxes(id)%cc(:, :,
mg_iold) = &
3300 mg%boxes(id)%cc(:, :,
mg_iphi)
3302 end subroutine update_coarse
3305 subroutine correct_children(mg, lvl)
3306 type(
mg_t),
intent(inout) :: mg
3307 integer,
intent(in) :: lvl
3310 do i = 1,
size(mg%lvls(lvl)%my_parents)
3311 id = mg%lvls(lvl)%my_parents(i)
3314 mg%boxes(id)%cc(:, :,
mg_ires) = &
3315 mg%boxes(id)%cc(:, :,
mg_iphi) - &
3316 mg%boxes(id)%cc(:, :,
mg_iold)
3320 end subroutine correct_children
3322 subroutine smooth_boxes(mg, lvl, n_cycle)
3323 type(
mg_t),
intent(inout) :: mg
3324 integer,
intent(in) :: lvl
3325 integer,
intent(in) :: n_cycle
3326 integer :: n, i, id, nc
3328 nc = mg%box_size_lvl(lvl)
3330 do n = 1, n_cycle * mg%n_smoother_substeps
3332 do i = 1,
size(mg%lvls(lvl)%my_ids)
3333 id = mg%lvls(lvl)%my_ids(i)
3334 call mg%box_smoother(mg, id, nc, n)
3342 end subroutine smooth_boxes
3344 subroutine residual_box(mg, id, nc)
3345 type(
mg_t),
intent(inout) :: mg
3346 integer,
intent(in) :: id
3347 integer,
intent(in) :: nc
3349 call mg%box_op(mg, id, nc,
mg_ires)
3351 mg%boxes(id)%cc(1:nc, 1:nc,
mg_ires) = &
3352 mg%boxes(id)%cc(1:nc, 1:nc,
mg_irhs) &
3353 - mg%boxes(id)%cc(1:nc, 1:nc,
mg_ires)
3354 end subroutine residual_box
3358 type(
mg_t),
intent(inout) :: mg
3359 integer,
intent(in) :: i_out
3360 procedure(mg_box_op),
optional :: op
3361 integer :: lvl, i, id, nc
3363 do lvl = mg%lowest_lvl, mg%highest_lvl
3364 nc = mg%box_size_lvl(lvl)
3365 do i = 1,
size(mg%lvls(lvl)%my_ids)
3366 id = mg%lvls(lvl)%my_ids(i)
3367 if (
present(op))
then
3368 call op(mg, id, nc, i_out)
3370 call mg%box_op(mg, id, nc, i_out)
3380 type(
mg_t),
intent(inout) :: mg
3381 integer,
intent(out) :: n_send(0:mg%n_cpu-1)
3382 integer,
intent(out) :: n_recv(0:mg%n_cpu-1)
3383 integer,
intent(out) :: dsize
3384 integer :: n_out(0:mg%n_cpu-1, mg%first_normal_lvl:mg%highest_lvl)
3385 integer :: n_in(0:mg%n_cpu-1, mg%first_normal_lvl:mg%highest_lvl)
3386 integer :: lvl, i, id, p_id, p_rank
3387 integer :: i_c, c_id, c_rank, min_lvl
3391 min_lvl = max(mg%lowest_lvl+1, mg%first_normal_lvl)
3393 do lvl = min_lvl, mg%highest_lvl
3395 do i = 1,
size(mg%lvls(lvl-1)%my_parents)
3396 id = mg%lvls(lvl-1)%my_parents(i)
3398 c_id = mg%boxes(id)%children(i_c)
3401 c_rank = mg%boxes(c_id)%rank
3402 if (c_rank /= mg%my_rank)
then
3403 n_in(c_rank, lvl) = n_in(c_rank, lvl) + 1
3410 do i = 1,
size(mg%lvls(lvl)%my_ids)
3411 id = mg%lvls(lvl)%my_ids(i)
3413 p_id = mg%boxes(id)%parent
3414 p_rank = mg%boxes(p_id)%rank
3415 if (p_rank /= mg%my_rank)
then
3416 n_out(p_rank, lvl) = n_out(p_rank, lvl) + 1
3422 allocate(mg%comm_restrict%n_send(0:mg%n_cpu-1, &
3423 mg%first_normal_lvl:mg%highest_lvl))
3424 allocate(mg%comm_restrict%n_recv(0:mg%n_cpu-1, &
3425 mg%first_normal_lvl:mg%highest_lvl))
3426 mg%comm_restrict%n_send = n_out
3427 mg%comm_restrict%n_recv = n_in
3429 dsize = (mg%box_size/2)**2
3430 n_send = maxval(n_out, dim=2)
3431 n_recv = maxval(n_in, dim=2)
3436 type(
mg_t),
intent(inout) :: mg
3437 integer,
intent(in) :: iv
3440 do lvl = mg%highest_lvl, mg%lowest_lvl+1, -1
3448 type(
mg_t),
intent(inout) :: mg
3449 integer,
intent(in) :: iv
3450 integer,
intent(in) :: lvl
3451 integer :: i, id, dsize, nc
3453 if (lvl <= mg%lowest_lvl) error stop
"cannot restrict lvl <= lowest_lvl"
3455 nc = mg%box_size_lvl(lvl)
3457 if (lvl >= mg%first_normal_lvl)
then
3460 mg%buf(:)%i_send = 0
3463 do i = 1,
size(mg%lvls(lvl)%my_ids)
3464 id = mg%lvls(lvl)%my_ids(i)
3465 call restrict_set_buffer(mg, id, iv)
3468 mg%buf(:)%i_recv = mg%comm_restrict%n_recv(:, lvl) * dsize
3470 mg%buf(:)%i_recv = 0
3473 do i = 1,
size(mg%lvls(lvl-1)%my_parents)
3474 id = mg%lvls(lvl-1)%my_parents(i)
3475 call restrict_onto(mg, id, nc, iv)
3479 subroutine restrict_set_buffer(mg, id, iv)
3480 type(
mg_t),
intent(inout) :: mg
3481 integer,
intent(in) :: id
3482 integer,
intent(in) :: iv
3483 integer :: i, j, n, hnc, p_id, p_rank
3484 real(
dp) :: tmp(mg%box_size/2, mg%box_size/2)
3487 p_id = mg%boxes(id)%parent
3488 p_rank = mg%boxes(p_id)%rank
3490 if (p_rank /= mg%my_rank)
then
3493 tmp(i, j) = 0.25_dp * &
3494 sum(mg%boxes(id)%cc(2*i-1:2*i, 2*j-1:2*j, iv))
3500 i = mg%buf(p_rank)%i_send
3501 mg%buf(p_rank)%send(i+1:i+n) = pack(tmp, .true.)
3502 mg%buf(p_rank)%i_send = mg%buf(p_rank)%i_send + n
3505 i = mg%buf(p_rank)%i_ix
3508 mg%buf(p_rank)%i_ix = mg%buf(p_rank)%i_ix + 1
3510 end subroutine restrict_set_buffer
3512 subroutine restrict_onto(mg, id, nc, iv)
3513 type(
mg_t),
intent(inout) :: mg
3514 integer,
intent(in) :: id
3515 integer,
intent(in) :: nc
3516 integer,
intent(in) :: iv
3517 integer :: i, j, hnc, dsize, i_c, c_id
3518 integer :: c_rank, dix(2)
3524 c_id = mg%boxes(id)%children(i_c)
3526 c_rank = mg%boxes(c_id)%rank
3529 if (c_rank == mg%my_rank)
then
3530 do j=1, hnc;
do i=1, hnc
3531 mg%boxes(id)%cc(dix(1)+i, dix(2)+j, iv) = 0.25_dp * &
3532 sum(mg%boxes(c_id)%cc(2*i-1:2*i, 2*j-1:2*j, iv))
3535 i = mg%buf(c_rank)%i_recv
3536 mg%boxes(id)%cc(dix(1)+1:dix(1)+hnc, &
3537 dix(2)+1:dix(2)+hnc, iv) = &
3538 reshape(mg%buf(c_rank)%recv(i+1:i+dsize), [hnc, hnc])
3539 mg%buf(c_rank)%i_recv = mg%buf(c_rank)%i_recv + dsize
3543 end subroutine restrict_onto
3547 Subroutine i_mrgrnk (XDONT, IRNGT)
3554 Integer,
Dimension (:),
Intent (In) :: xdont
3555 Integer,
Dimension (:),
Intent (Out) :: irngt
3557 Integer :: xvala, xvalb
3559 Integer,
Dimension (SIZE(IRNGT)) :: jwrkt
3560 Integer :: lmtna, lmtnc, irng1, irng2
3561 Integer :: nval, iind, iwrkd, iwrk, iwrkf, jinda, iinda, iindb
3563 nval = min(
SIZE(xdont),
SIZE(irngt))
3576 Do iind = 2, nval, 2
3577 If (xdont(iind-1) <= xdont(iind))
Then
3578 irngt(iind-1) = iind - 1
3581 irngt(iind-1) = iind
3582 irngt(iind) = iind - 1
3585 If (modulo(nval, 2) /= 0)
Then
3602 Do iwrkd = 0, nval - 1, 4
3603 If ((iwrkd+4) > nval)
Then
3604 If ((iwrkd+2) >= nval)
Exit
3608 If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3)))
Exit
3612 If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3)))
Then
3613 irng2 = irngt(iwrkd+2)
3614 irngt(iwrkd+2) = irngt(iwrkd+3)
3615 irngt(iwrkd+3) = irng2
3620 irng1 = irngt(iwrkd+1)
3621 irngt(iwrkd+1) = irngt(iwrkd+3)
3622 irngt(iwrkd+3) = irngt(iwrkd+2)
3623 irngt(iwrkd+2) = irng1
3630 If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) cycle
3634 If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3)))
Then
3635 irng2 = irngt(iwrkd+2)
3636 irngt(iwrkd+2) = irngt(iwrkd+3)
3637 If (xdont(irng2) <= xdont(irngt(iwrkd+4)))
Then
3639 irngt(iwrkd+3) = irng2
3642 irngt(iwrkd+3) = irngt(iwrkd+4)
3643 irngt(iwrkd+4) = irng2
3649 irng1 = irngt(iwrkd+1)
3650 irng2 = irngt(iwrkd+2)
3651 irngt(iwrkd+1) = irngt(iwrkd+3)
3652 If (xdont(irng1) <= xdont(irngt(iwrkd+4)))
Then
3653 irngt(iwrkd+2) = irng1
3654 If (xdont(irng2) <= xdont(irngt(iwrkd+4)))
Then
3656 irngt(iwrkd+3) = irng2
3659 irngt(iwrkd+3) = irngt(iwrkd+4)
3660 irngt(iwrkd+4) = irng2
3664 irngt(iwrkd+2) = irngt(iwrkd+4)
3665 irngt(iwrkd+3) = irng1
3666 irngt(iwrkd+4) = irng2
3681 If (lmtna >= nval)
Exit
3690 jinda = iwrkf + lmtna
3691 iwrkf = iwrkf + lmtnc
3692 If (iwrkf >= nval)
Then
3693 If (jinda >= nval)
Exit
3709 jwrkt(1:lmtna) = irngt(iwrkd:jinda)
3711 xvala = xdont(jwrkt(iinda))
3712 xvalb = xdont(irngt(iindb))
3719 If (xvala > xvalb)
Then
3720 irngt(iwrk) = irngt(iindb)
3722 If (iindb > iwrkf)
Then
3724 irngt(iwrk+1:iwrkf) = jwrkt(iinda:lmtna)
3727 xvalb = xdont(irngt(iindb))
3729 irngt(iwrk) = jwrkt(iinda)
3731 If (iinda > lmtna) exit
3732 xvala = xdont(jwrkt(iinda))
3745 End Subroutine i_mrgrnk
3750 type(
mg_t),
intent(inout) :: mg
3752 if (mg%n_extra_vars == 0 .and. mg%is_allocated)
then
3753 error stop
"ahelmholtz_set_methods: mg%n_extra_vars == 0"
3755 mg%n_extra_vars = max(2, mg%n_extra_vars)
3767 select case (mg%geometry_type)
3769 mg%box_op => box_ahelmh
3771 select case (mg%smoother_type)
3773 mg%box_smoother => box_gs_ahelmh
3775 error stop
"ahelmholtz_set_methods: unsupported smoother type"
3778 error stop
"ahelmholtz_set_methods: unsupported geometry"
3784 real(
dp),
intent(in) :: lambda
3787 error stop
"ahelmholtz_set_lambda: lambda < 0 not allowed"
3793 subroutine box_gs_ahelmh(mg, id, nc, redblack_cntr)
3794 type(
mg_t),
intent(inout) :: mg
3795 integer,
intent(in) :: id
3796 integer,
intent(in) :: nc
3797 integer,
intent(in) :: redblack_cntr
3798 integer :: i, j, i0, di
3800 real(
dp) :: idr2(2*2), u(2*2)
3801 real(
dp) :: a0(2*2), a(2*2), c(2*2)
3804 idr2(1:2*2:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
3805 idr2(2:2*2:2) = idr2(1:2*2:2)
3817 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
3823 i0 = 2 - iand(ieor(redblack_cntr, j), 1)
3826 a0(1:2) = cc(i, j, i_eps1)
3827 a0(3:4) = cc(i, j, i_eps2)
3828 u(1:2) = cc(i-1:i+1:2, j, n)
3829 a(1:2) = cc(i-1:i+1:2, j, i_eps1)
3830 u(3:4) = cc(i, j-1:j+1:2, n)
3831 a(3:4) = cc(i, j-1:j+1:2, i_eps2)
3832 c(:) = 2 * a0(:) * a(:) / (a0(:) + a(:)) * idr2
3835 (sum(c(:) * u(:)) - cc(i, j,
mg_irhs)) / &
3840 end subroutine box_gs_ahelmh
3843 subroutine box_ahelmh(mg, id, nc, i_out)
3844 type(
mg_t),
intent(inout) :: mg
3845 integer,
intent(in) :: id
3846 integer,
intent(in) :: nc
3847 integer,
intent(in) :: i_out
3849 real(
dp) :: idr2(2*2), a0(2*2), u0, u(2*2), a(2*2)
3852 idr2(1:2*2:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
3853 idr2(2:2*2:2) = idr2(1:2*2:2)
3855 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
3860 a0(1:2) = cc(i, j, i_eps1)
3861 a0(3:4) = cc(i, j, i_eps2)
3862 a(1:2) = cc(i-1:i+1:2, j, i_eps1)
3863 a(3:4) = cc(i, j-1:j+1:2, i_eps2)
3865 u(1:2) = cc(i-1:i+1:2, j, n)
3866 u(3:4) = cc(i, j-1:j+1:2, n)
3868 cc(i, j, i_out) = sum(2 * idr2 * &
3869 a0(:)*a(:)/(a0(:) + a(:)) * (u(:) - u0)) - &
3874 end subroutine box_ahelmh
To fill ghost cells near physical boundaries.
subroutine, public mg_comm_init(mg, comm)
Prolong from a parent to a child with index offset dix. This method could sometimes work better than ...
subroutine, public mg_fas_vcycle(mg, highest_lvl, max_res, standalone)
Perform FAS V-cycle (full approximation scheme).
subroutine, public vhelmholtz_set_methods(mg)
subroutine, public mg_ghost_cell_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
integer, parameter, public mg_num_children
subroutine, public mg_timer_end(timer)
subroutine, public mg_restrict(mg, iv)
Restrict all levels.
integer, parameter, public mg_vhelmholtz
Indicates a variable-coefficient Helmholtz equation.
integer, parameter, public mg_iveps1
Indexes of anisotropic variable coefficients.
integer, parameter, public mg_ahelmholtz
Indicates a anisotropic-coefficient Helmholtz equation.
integer, parameter, public mg_vlaplacian
Indicates a variable-coefficient Laplacian.
integer, dimension(2, 4), parameter, public mg_neighb_dix
integer, parameter, public mg_helmholtz
Indicates a constant-coefficient Helmholtz equation.
integer, parameter, public mg_physical_boundary
Special value that indicates there is a physical boundary.
subroutine, public mg_prolong_sparse(mg, p_id, dix, nc, iv, fine)
Prolong from a parent to a child with index offset dix.
subroutine, public mg_get_face_coords(box, nb, nc, x)
Get coordinates at the face of a box.
subroutine, public sort_and_transfer_buffers(mg, dsize)
integer, parameter, public dp
Type of reals.
pure integer function, public mg_highest_uniform_lvl(mg)
integer function, public mg_ix_to_ichild(ix)
Compute the child index for a box with spatial index ix. With child index we mean the index in the ch...
integer, parameter, public mg_ndim
Problem dimension.
subroutine, public mg_prolong_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
subroutine, public diffusion_solve_acoeff(mg, dt, order, max_res)
Solve a diffusion equation implicitly, assuming anisotropic diffusion coefficient which has been stor...
subroutine, public mg_fas_fmg(mg, have_guess, max_res)
Perform FAS-FMG cycle (full approximation scheme, full multigrid).
integer, parameter, public mg_ires
Index of residual.
integer, parameter, public mg_neighb_highy
subroutine, public mg_set_neighbors_lvl(mg, lvl)
integer, parameter, public mg_iveps3
subroutine, public mg_load_balance_simple(mg)
Load balance all boxes in the multigrid tree, by simply distributing the load per grid level....
subroutine, public helmholtz_set_methods(mg)
integer, parameter, public mg_bc_continuous
Value to indicate a continuous boundary condition.
subroutine, public mg_timers_show(mg)
integer, parameter, public mg_iphi
Index of solution.
subroutine, public mg_add_children(mg, id)
subroutine, public mg_prolong(mg, lvl, iv, iv_to, method, add)
Prolong variable iv from lvl to variable iv_to at lvl+1.
integer, dimension(4), parameter, public mg_neighb_dim
integer, dimension(2, 4), parameter, public mg_child_dix
integer, parameter, public mg_spherical
Spherical coordinate system.
subroutine, public mg_timer_start(timer)
integer(i8) function, public mg_number_of_unknowns(mg)
Determine total number of unknowns (on leaves)
integer, parameter, public mg_neighb_highx
integer, parameter, public mg_bc_dirichlet
Value to indicate a Dirichlet boundary condition.
subroutine, public helmholtz_set_lambda(lambda)
integer, parameter, public mg_laplacian
Indicates a standard Laplacian.
integer, dimension(4), parameter, public mg_neighb_rev
subroutine, public ahelmholtz_set_lambda(lambda)
real(dp), public, protected ahelmholtz_lambda
Module which contains multigrid procedures for a Helmholtz operator of the form: div(D grad(phi)) - l...
integer, parameter, public mg_lvl_hi
Maximum allowed grid level.
integer, dimension(4), parameter, public mg_neighb_high_pm
integer, parameter, public mg_bc_neumann
Value to indicate a Neumann boundary condition.
integer, parameter, public mg_irhs
Index of right-hand side.
logical, dimension(2, 4), parameter, public mg_child_low
integer, parameter, public mg_iveps2
elemental logical function, public mg_has_children(box)
Return .true. if a box has children.
subroutine sort_sendbuf(gc, dsize)
Sort send buffers according to the idbuf array.
subroutine, public diffusion_solve_vcoeff(mg, dt, order, max_res)
Solve a diffusion equation implicitly, assuming a variable diffusion coefficient which has been store...
integer, parameter, public mg_cylindrical
Cylindrical coordinate system.
real(dp), public, protected vhelmholtz_lambda
Module for implicitly solving diffusion equations.
subroutine, public mg_restrict_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
integer, parameter, public mg_smoother_jacobi
subroutine, public laplacian_set_methods(mg)
subroutine, public mg_deallocate_storage(mg)
Deallocate all allocatable arrays.
subroutine, public ahelmholtz_set_methods(mg)
subroutine, public mg_load_balance(mg)
Load balance all boxes in the multigrid tree. Compared to mg_load_balance_simple, this method does a ...
subroutine, public mg_fill_ghost_cells_lvl(mg, lvl, iv)
Fill ghost cells at a grid level.
subroutine, public vhelmholtz_set_lambda(lambda)
subroutine, public mg_set_next_level_ids(mg, lvl)
integer, parameter, public mg_num_neighbors
integer, parameter, public mg_cartesian
Cartesian coordinate system.
integer, dimension(2, 4), parameter, public mg_child_adj_nb
subroutine, public mg_set_refinement_boundaries(boxes, level)
Create a list of refinement boundaries (from the coarse side)
integer, parameter, public mg_neighb_lowx
integer, parameter, public mg_neighb_lowy
integer, parameter, public mg_num_vars
Number of predefined multigrid variables.
subroutine subtract_mean(mg, iv, include_ghostcells)
integer, parameter, public mg_lvl_lo
Minimum allowed grid level.
subroutine, public mg_set_leaves_parents(boxes, level)
Create a list of leaves and a list of parents for a level.
subroutine, public mg_load_balance_parents(mg)
Load balance the parents (non-leafs). Assign them to the rank that has most children.
subroutine, public mg_fill_ghost_cells(mg, iv)
Fill ghost cells at all grid levels.
subroutine, public mg_restrict_lvl(mg, iv, lvl)
Restrict from lvl to lvl-1.
subroutine, public vlaplacian_set_methods(mg)
integer, parameter, public mg_iveps
Index of the variable coefficient (at cell centers)
subroutine, public mg_allocate_storage(mg)
Allocate communication buffers and local boxes for a tree that has already been created.
real(dp) function get_sum(mg, iv)
subroutine, public mg_build_rectangle(mg, domain_size, box_size, dx, r_min, periodic, n_finer)
subroutine, public mg_apply_op(mg, i_out, op)
Apply operator to the tree and store in variable i_out.
integer, parameter, public i8
Type for 64-bit integers.
subroutine, public mg_set_methods(mg)
integer, parameter, public mg_max_num_vars
Maximum number of variables.
pure integer function, dimension(2), public mg_get_child_offset(mg, id)
Get the offset of a box with respect to its parent (e.g. in 2d, there can be a child at offset 0,...
integer, parameter, public mg_iold
Index of previous solution (used for correction)
logical, dimension(4), parameter, public mg_neighb_low
integer, parameter, public mg_max_timers
Maximum number of timers to use.
integer, parameter, public mg_smoother_gs
real(dp), public, protected helmholtz_lambda
Module for load balancing a tree (that already has been constructed). The load balancing determines w...
subroutine, public diffusion_solve(mg, dt, diffusion_coeff, order, max_res)
Solve a diffusion equation implicitly, assuming a constant diffusion coefficient. The solution at tim...
subroutine, public mg_phi_bc_store(mg)
Store boundary conditions for the solution variable, this can speed up calculations if the same bound...
integer, dimension(4, 2), parameter, public mg_child_rev
integer, parameter, public mg_no_box
Special value that indicates there is no box.
integer, parameter, public mg_smoother_gsrb
integer function, public mg_add_timer(mg, name)
Buffer type (one is used for each pair of communicating processes)
Lists of blocks per refinement level.