30 integer,
parameter,
public ::
dp = kind(0.0d0)
33 integer,
parameter,
public ::
i8 = selected_int_kind(18)
111 integer,
parameter,
public ::
mg_child_dix(1, 2) = reshape([0,1], [1,2])
113 integer,
parameter,
public ::
mg_child_rev(2, 1) = reshape([2,1], [2,1])
117 logical,
parameter,
public ::
mg_child_low(1, 2) = reshape([.true., .false.], [1, 2])
138 integer,
allocatable :: leaves(:)
139 integer,
allocatable :: parents(:)
140 integer,
allocatable :: ref_bnds(:)
141 integer,
allocatable :: ids(:)
142 integer,
allocatable :: my_leaves(:)
143 integer,
allocatable :: my_parents(:)
144 integer,
allocatable :: my_ref_bnds(:)
145 integer,
allocatable :: my_ids(:)
155 integer :: children(2**1)
156 integer :: neighbors(2*1)
160 real(
dp),
allocatable :: cc(:, :)
168 integer,
allocatable :: ix(:)
169 real(
dp),
allocatable :: send(:)
170 real(
dp),
allocatable :: recv(:)
174 integer,
allocatable :: n_send(:, :)
175 integer,
allocatable :: n_recv(:, :)
180 real(
dp) :: bc_value = 0.0_dp
182 procedure(
mg_subr_bc),
pointer,
nopass :: boundary_cond => null()
184 procedure(
mg_subr_rb),
pointer,
nopass :: refinement_bnd => null()
188 character(len=20) :: name
189 real(
dp) :: t = 0.0_dp
195 logical :: tree_created = .false.
197 logical :: is_allocated = .false.
199 integer :: n_extra_vars = 0
203 integer :: n_cpu = -1
205 integer :: my_rank = -1
207 integer :: box_size = -1
209 integer :: highest_lvl = -1
211 integer :: lowest_lvl = -1
214 integer :: first_normal_lvl = -1
216 integer :: n_boxes = 0
241 logical :: phi_bc_data_stored = .false.
244 logical :: periodic(1) = .false.
255 logical :: subtract_mean = .false.
259 integer :: n_smoother_substeps = 1
261 integer :: n_cycle_down = 2
263 integer :: n_cycle_up = 2
265 integer :: max_coarse_cycles = 1000
266 integer :: coarsest_grid(1) = 2
268 real(
dp) :: residual_coarse_abs = 1e-8_dp
270 real(
dp) :: residual_coarse_rel = 1e-8_dp
273 procedure(
mg_box_op),
pointer,
nopass :: box_op => null()
282 integer :: n_timers = 0
292 integer,
intent(in) :: nc
293 integer,
intent(in) :: iv
294 integer,
intent(in) :: nb
295 integer,
intent(out) :: bc_type
297 real(dp),
intent(out) :: bc(1)
303 type(
mg_box_t),
intent(inout) :: box
304 integer,
intent(in) :: nc
305 integer,
intent(in) :: iv
306 integer,
intent(in) :: nb
308 real(dp),
intent(in) :: cgc(1)
314 type(
mg_t),
intent(inout) :: mg
315 integer,
intent(in) :: id
316 integer,
intent(in) :: nc
317 integer,
intent(in) :: i_out
323 type(
mg_t),
intent(inout) :: mg
324 integer,
intent(in) :: id
325 integer,
intent(in) :: nc
326 integer,
intent(in) :: redblack_cntr
332 type(
mg_t),
intent(inout) :: mg
333 integer,
intent(in) :: p_id
334 integer,
intent(in) :: dix(1)
335 integer,
intent(in) :: nc
336 integer,
intent(in) :: iv
337 real(dp),
intent(out) :: fine(nc)
450 integer :: timer_total_vcycle = -1
451 integer :: timer_total_fmg = -1
452 integer :: timer_smoother = -1
453 integer :: timer_smoother_gc = -1
454 integer :: timer_coarse = -1
455 integer :: timer_correct = -1
456 integer :: timer_update_coarse = -1
473 Integer,
Parameter :: kdp = selected_real_kind(15)
478 module procedure i_mrgrnk
508 integer,
intent(in) :: ix(1)
517 type(
mg_t),
intent(in) :: mg
518 integer,
intent(in) :: id
519 integer :: ix_offset(1)
521 if (mg%boxes(id)%lvl <= mg%first_normal_lvl)
then
524 ix_offset = iand(mg%boxes(id)%ix-1, 1) * &
525 ishft(mg%box_size, -1)
530 type(
mg_t),
intent(in) :: mg
533 do lvl = mg%first_normal_lvl, mg%highest_lvl-1
535 if (
size(mg%lvls(lvl)%leaves) /= 0 .and. &
536 size(mg%lvls(lvl)%parents) /= 0)
exit
543 type(
mg_t),
intent(in) :: mg
545 integer(i8) :: n_unknowns
548 do lvl = mg%first_normal_lvl, mg%highest_lvl
549 n_unknowns = n_unknowns +
size(mg%lvls(lvl)%leaves)
551 n_unknowns = n_unknowns * int(mg%box_size**3,
i8)
557 integer,
intent(in) :: nb
558 integer,
intent(in) :: nc
559 real(
dp),
intent(out) :: x(2)
568 rmin(nb_dim) = rmin(nb_dim) + box%dr(nb_dim) * nc
572 x(2) = rmin(1) + box%dr(1) * nc
576 type(
mg_t),
intent(inout) :: mg
577 character(len=*),
intent(in) :: name
579 mg%n_timers = mg%n_timers + 1
587 timer%t0 = mpi_wtime()
593 timer%t = timer%t + mpi_wtime() - timer%t0
598 type(
mg_t),
intent(in) :: mg
600 real(
dp) :: tmin(mg%n_timers)
601 real(
dp) :: tmax(mg%n_timers)
602 real(
dp),
allocatable :: t_array(:)
604 allocate(t_array(mg%n_timers))
605 do n = 1, mg%n_timers
606 t_array(n) = mg%timers(n)%t
608 call mpi_reduce(t_array, tmin, mg%n_timers, mpi_double, mpi_min, 0, mg%comm, ierr)
609 call mpi_reduce(t_array, tmax, mg%n_timers, mpi_double, mpi_max, 0, mg%comm, ierr)
611 if (mg%my_rank == 0)
then
612 write(*,
"(A20,2A16)")
"name ",
"min(s)",
"max(s)"
613 do n = 1, mg%n_timers
614 write(*,
"(A20,2F16.6)") mg%timers(n)%name, &
624 type(
mg_t),
intent(inout) :: mg
627 if (.not. mg%is_allocated) &
628 error stop
"deallocate_storage: tree is not allocated"
633 deallocate(mg%comm_restrict%n_send)
634 deallocate(mg%comm_restrict%n_recv)
636 deallocate(mg%comm_prolong%n_send)
637 deallocate(mg%comm_prolong%n_recv)
639 deallocate(mg%comm_ghostcell%n_send)
640 deallocate(mg%comm_ghostcell%n_recv)
642 do lvl = mg%lowest_lvl, mg%highest_lvl
643 deallocate(mg%lvls(lvl)%ids)
644 deallocate(mg%lvls(lvl)%leaves)
645 deallocate(mg%lvls(lvl)%parents)
646 deallocate(mg%lvls(lvl)%ref_bnds)
647 deallocate(mg%lvls(lvl)%my_ids)
648 deallocate(mg%lvls(lvl)%my_leaves)
649 deallocate(mg%lvls(lvl)%my_parents)
650 deallocate(mg%lvls(lvl)%my_ref_bnds)
653 mg%is_allocated = .false.
655 mg%phi_bc_data_stored = .false.
662 type(
mg_t),
intent(inout) :: mg
663 integer :: i, id, lvl, nc
664 integer :: n_send(0:mg%n_cpu-1, 3)
665 integer :: n_recv(0:mg%n_cpu-1, 3)
667 integer :: n_in, n_out, n_id
669 if (.not. mg%tree_created) &
670 error stop
"allocate_storage: tree is not yet created"
672 if (mg%is_allocated) &
673 error stop
"allocate_storage: tree is already allocated"
675 do lvl = mg%lowest_lvl, mg%highest_lvl
676 nc = mg%box_size_lvl(lvl)
677 do i = 1,
size(mg%lvls(lvl)%my_ids)
678 id = mg%lvls(lvl)%my_ids(i)
679 allocate(mg%boxes(id)%cc(0:nc+1, &
683 mg%boxes(id)%cc(:, :) = 0.0_dp
687 allocate(mg%buf(0:mg%n_cpu-1))
690 n_recv(:, 1), dsize(1))
692 n_recv(:, 2), dsize(2))
694 n_recv(:, 3), dsize(3))
697 n_out = maxval(n_send(i, :) * dsize(:))
698 n_in = maxval(n_recv(i, :) * dsize(:))
699 n_id = maxval(n_send(i, :))
700 allocate(mg%buf(i)%send(n_out))
701 allocate(mg%buf(i)%recv(n_in))
702 allocate(mg%buf(i)%ix(n_id))
705 mg%is_allocated = .true.
715 type(
mg_t),
intent(inout) :: mg
716 real(
dp),
intent(in) :: dt
717 real(
dp),
intent(in) :: diffusion_coeff
718 integer,
intent(in) :: order
719 real(
dp),
intent(in) :: max_res
720 integer,
parameter :: max_its = 10
730 call set_rhs(mg, -1/(dt * diffusion_coeff), 0.0_dp)
735 call set_rhs(mg, -2/(dt * diffusion_coeff), -1.0_dp)
737 error stop
"diffusion_solve: order should be 1 or 2"
745 if (res <= max_res)
exit
749 if (n == max_its + 1)
then
750 if (mg%my_rank == 0) &
751 print *,
"Did you specify boundary conditions correctly?"
752 error stop
"diffusion_solve: no convergence"
762 type(
mg_t),
intent(inout) :: mg
763 real(
dp),
intent(in) :: dt
764 integer,
intent(in) :: order
765 real(
dp),
intent(in) :: max_res
766 integer,
parameter :: max_its = 10
776 call set_rhs(mg, -1/dt, 0.0_dp)
781 call set_rhs(mg, -2/dt, -1.0_dp)
783 error stop
"diffusion_solve: order should be 1 or 2"
791 if (res <= max_res)
exit
795 if (n == max_its + 1)
then
796 if (mg%my_rank == 0)
then
797 print *,
"Did you specify boundary conditions correctly?"
798 print *,
"Or is the variation in diffusion too large?"
800 error stop
"diffusion_solve: no convergence"
811 type(
mg_t),
intent(inout) :: mg
812 real(
dp),
intent(in) :: dt
813 integer,
intent(in) :: order
814 real(
dp),
intent(in) :: max_res
815 integer,
parameter :: max_its = 10
825 call set_rhs(mg, -1/dt, 0.0_dp)
830 call set_rhs(mg, -2/dt, -1.0_dp)
832 error stop
"diffusion_solve: order should be 1 or 2"
840 if (res <= max_res)
exit
844 if (n == max_its + 1)
then
845 if (mg%my_rank == 0)
then
846 print *,
"Did you specify boundary conditions correctly?"
847 print *,
"Or is the variation in diffusion too large?"
849 error stop
"diffusion_solve: no convergence"
853 subroutine set_rhs(mg, f1, f2)
854 type(
mg_t),
intent(inout) :: mg
855 real(
dp),
intent(in) :: f1, f2
856 integer :: lvl, i, id, nc
858 do lvl = 1, mg%highest_lvl
859 nc = mg%box_size_lvl(lvl)
860 do i = 1,
size(mg%lvls(lvl)%my_leaves)
861 id = mg%lvls(lvl)%my_leaves(i)
862 mg%boxes(id)%cc(1:nc,
mg_irhs) = &
863 f1 * mg%boxes(id)%cc(1:nc,
mg_iphi) + &
864 f2 * mg%boxes(id)%cc(1:nc,
mg_irhs)
867 end subroutine set_rhs
872 type(
mg_t),
intent(inout) :: mg
874 if (all(mg%periodic))
then
877 mg%subtract_mean = .true.
880 select case (mg%geometry_type)
884 select case (mg%smoother_type)
888 mg%box_smoother => box_gs_lpl
890 error stop
"laplacian_set_methods: unsupported smoother type"
893 error stop
"laplacian_set_methods: unsupported geometry"
899 subroutine box_gs_lpl(mg, id, nc, redblack_cntr)
900 type(
mg_t),
intent(inout) :: mg
901 integer,
intent(in) :: id
902 integer,
intent(in) :: nc
903 integer,
intent(in) :: redblack_cntr
905 real(
dp) :: idr2(1), fac
908 idr2 = 1/mg%dr(:, mg%boxes(id)%lvl)**2
909 fac = 0.5_dp / sum(idr2)
921 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
922 if (redblack) i0 = 2 - iand(redblack_cntr, 1)
926 idr2(1) * (cc(i+1, n) + cc(i-1, n)) - &
930 end subroutine box_gs_lpl
971 subroutine box_lpl(mg, id, nc, i_out)
972 type(
mg_t),
intent(inout) :: mg
973 integer,
intent(in) :: id
974 integer,
intent(in) :: nc
975 integer,
intent(in) :: i_out
979 idr2 = 1 / mg%dr(:, mg%boxes(id)%lvl)**2
981 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
984 idr2(1) * (cc(i-1, n) + cc(i+1, n) - 2 * cc(i, n))
987 end subroutine box_lpl
993 type(
mg_t),
intent(inout) :: mg
995 if (mg%n_extra_vars == 0 .and. mg%is_allocated)
then
996 error stop
"vhelmholtz_set_methods: mg%n_extra_vars == 0"
998 mg%n_extra_vars = max(1, mg%n_extra_vars)
1001 mg%subtract_mean = .false.
1006 mg%bc(:,
mg_iveps)%bc_value = 0.0_dp
1008 select case (mg%geometry_type)
1010 mg%box_op => box_vhelmh
1012 select case (mg%smoother_type)
1014 mg%box_smoother => box_gs_vhelmh
1016 error stop
"vhelmholtz_set_methods: unsupported smoother type"
1019 error stop
"vhelmholtz_set_methods: unsupported geometry"
1025 real(
dp),
intent(in) :: lambda
1028 error stop
"vhelmholtz_set_lambda: lambda < 0 not allowed"
1034 subroutine box_gs_vhelmh(mg, id, nc, redblack_cntr)
1035 type(
mg_t),
intent(inout) :: mg
1036 integer,
intent(in) :: id
1037 integer,
intent(in) :: nc
1038 integer,
intent(in) :: redblack_cntr
1039 integer :: i, i0, di
1041 real(
dp) :: idr2(2*1), u(2*1)
1042 real(
dp) :: a0, a(2*1), c(2*1)
1045 idr2(1:2*1:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1046 idr2(2:2*1:2) = idr2(1:2*1:2)
1058 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
1060 if (redblack) i0 = 2 - iand(redblack_cntr, 1)
1064 u(1:2) = cc(i-1:i+1:2, n)
1065 a(1:2) = cc(i-1:i+1:2, i_eps)
1066 c(:) = 2 * a0 * a(:) / (a0 + a(:)) * idr2
1069 (sum(c(:) * u(:)) - cc(i,
mg_irhs)) / &
1073 end subroutine box_gs_vhelmh
1076 subroutine box_vhelmh(mg, id, nc, i_out)
1077 type(
mg_t),
intent(inout) :: mg
1078 integer,
intent(in) :: id
1079 integer,
intent(in) :: nc
1080 integer,
intent(in) :: i_out
1082 real(
dp) :: idr2(2*1), a0, u0, u(2*1), a(2*1)
1085 idr2(1:2*1:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1086 idr2(2:2*1:2) = idr2(1:2*1:2)
1088 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
1092 a(1:2) = cc(i-1:i+1:2, i_eps)
1094 u(1:2) = cc(i-1:i+1:2, n)
1096 cc(i, i_out) = sum(2 * idr2 * &
1097 a0*a(:)/(a0 + a(:)) * (u(:) - u0)) - &
1101 end subroutine box_vhelmh
1107 type(
mg_t),
intent(inout) :: mg
1108 integer,
intent(in) :: domain_size(1)
1109 integer,
intent(in) :: box_size
1110 real(
dp),
intent(in) :: dx(1)
1111 real(
dp),
intent(in) :: r_min(1)
1112 logical,
intent(in) :: periodic(1)
1113 integer,
intent(in) :: n_finer
1114 integer :: i, lvl, n, id, nx(1)
1115 integer :: boxes_per_dim(1,
mg_lvl_lo:1)
1116 integer :: periodic_offset(1)
1118 if (modulo(box_size, 2) /= 0) &
1119 error stop
"box_size should be even"
1120 if (any(modulo(domain_size, box_size) /= 0)) &
1121 error stop
"box_size does not divide domain_size"
1123 if (all(periodic))
then
1124 mg%subtract_mean = .true.
1128 mg%box_size = box_size
1129 mg%box_size_lvl(1) = box_size
1130 mg%domain_size_lvl(:, 1) = domain_size
1131 mg%first_normal_lvl = 1
1134 mg%periodic = periodic
1135 boxes_per_dim(:, :) = 0
1136 boxes_per_dim(:, 1) = domain_size / box_size
1141 if (any(modulo(nx, 2) == 1 .or. nx == mg%coarsest_grid .or. &
1142 (mg%box_size_lvl(lvl) == mg%coarsest_grid .and. &
1145 if (all(modulo(nx/mg%box_size_lvl(lvl), 2) == 0))
then
1146 mg%box_size_lvl(lvl-1) = mg%box_size_lvl(lvl)
1147 boxes_per_dim(:, lvl-1) = boxes_per_dim(:, lvl)/2
1148 mg%first_normal_lvl = lvl-1
1150 mg%box_size_lvl(lvl-1) = mg%box_size_lvl(lvl)/2
1151 boxes_per_dim(:, lvl-1) = boxes_per_dim(:, lvl)
1154 mg%dr(:, lvl-1) = mg%dr(:, lvl) * 2
1156 mg%domain_size_lvl(:, lvl-1) = nx
1163 mg%dr(:, lvl) = mg%dr(:, lvl-1) * 0.5_dp
1164 mg%box_size_lvl(lvl) = box_size
1165 mg%domain_size_lvl(:, lvl) = 2 * mg%domain_size_lvl(:, lvl-1)
1168 n = sum(product(boxes_per_dim, dim=1)) + n_finer
1169 allocate(mg%boxes(n))
1172 nx = boxes_per_dim(:, mg%lowest_lvl)
1173 periodic_offset = [nx(1)-1]
1176 mg%n_boxes = mg%n_boxes + 1
1179 mg%boxes(n)%rank = 0
1181 mg%boxes(n)%lvl = mg%lowest_lvl
1182 mg%boxes(n)%ix(:) = [i]
1183 mg%boxes(n)%r_min(:) = r_min + (mg%boxes(n)%ix(:) - 1) * &
1184 mg%box_size_lvl(mg%lowest_lvl) * mg%dr(:, mg%lowest_lvl)
1185 mg%boxes(n)%dr(:) = mg%dr(:, mg%lowest_lvl)
1190 mg%boxes(n)%neighbors(:) = [n-1, n+1]
1193 if(i==1 .and. .not.periodic(1))
then
1196 if(i==1 .and. periodic(1))
then
1197 mg%boxes(n)%neighbors(1) = n + periodic_offset(1)
1199 if(i==nx(1) .and. .not.periodic(1))
then
1202 if(i==nx(1) .and. periodic(1))
then
1203 mg%boxes(n)%neighbors(2) = n - periodic_offset(1)
1224 mg%lvls(mg%lowest_lvl)%ids = [(n, n=1, mg%n_boxes)]
1227 do lvl = mg%lowest_lvl, 0
1228 if (mg%box_size_lvl(lvl+1) == mg%box_size_lvl(lvl))
then
1229 do i = 1,
size(mg%lvls(lvl)%ids)
1230 id = mg%lvls(lvl)%ids(i)
1238 do i = 1,
size(mg%lvls(lvl)%ids)
1239 id = mg%lvls(lvl)%ids(i)
1240 call add_single_child(mg, id,
size(mg%lvls(lvl)%ids))
1251 do lvl = mg%lowest_lvl, 1
1252 if (
allocated(mg%lvls(lvl)%ref_bnds)) &
1253 deallocate(mg%lvls(lvl)%ref_bnds)
1254 allocate(mg%lvls(lvl)%ref_bnds(0))
1257 mg%tree_created = .true.
1261 type(
mg_t),
intent(inout) :: mg
1262 integer,
intent(in) :: lvl
1265 do i = 1,
size(mg%lvls(lvl)%ids)
1266 id = mg%lvls(lvl)%ids(i)
1267 call set_neighbs(mg%boxes, id)
1272 type(
mg_t),
intent(inout) :: mg
1273 integer,
intent(in) :: lvl
1276 if (
allocated(mg%lvls(lvl+1)%ids)) &
1277 deallocate(mg%lvls(lvl+1)%ids)
1280 if (mg%box_size_lvl(lvl+1) == mg%box_size_lvl(lvl))
then
1282 allocate(mg%lvls(lvl+1)%ids(n))
1285 do i = 1,
size(mg%lvls(lvl)%parents)
1286 id = mg%lvls(lvl)%parents(i)
1287 mg%lvls(lvl+1)%ids(n*(i-1)+1:n*i) = mg%boxes(id)%children
1290 n =
size(mg%lvls(lvl)%parents)
1291 allocate(mg%lvls(lvl+1)%ids(n))
1294 do i = 1,
size(mg%lvls(lvl)%parents)
1295 id = mg%lvls(lvl)%parents(i)
1296 mg%lvls(lvl+1)%ids(i) = mg%boxes(id)%children(1)
1303 subroutine set_neighbs(boxes, id)
1304 type(
mg_box_t),
intent(inout) :: boxes(:)
1305 integer,
intent(in) :: id
1306 integer :: nb, nb_id
1309 if (boxes(id)%neighbors(nb) ==
mg_no_box)
then
1310 nb_id = find_neighb(boxes, id, nb)
1312 boxes(id)%neighbors(nb) = nb_id
1317 end subroutine set_neighbs
1320 function find_neighb(boxes, id, nb)
result(nb_id)
1321 type(
mg_box_t),
intent(in) :: boxes(:)
1322 integer,
intent(in) :: id
1323 integer,
intent(in) :: nb
1324 integer :: nb_id, p_id, c_ix, d, old_pid
1326 p_id = boxes(id)%parent
1334 p_id = boxes(p_id)%neighbors(nb)
1339 end function find_neighb
1343 type(
mg_box_t),
intent(in) :: boxes(:)
1344 type(
mg_lvl_t),
intent(inout) :: level
1345 integer :: i, id, i_leaf, i_parent
1346 integer :: n_parents, n_leaves
1349 n_leaves =
size(level%ids) - n_parents
1351 if (.not.
allocated(level%parents))
then
1352 allocate(level%parents(n_parents))
1353 else if (n_parents /=
size(level%parents))
then
1354 deallocate(level%parents)
1355 allocate(level%parents(n_parents))
1358 if (.not.
allocated(level%leaves))
then
1359 allocate(level%leaves(n_leaves))
1360 else if (n_leaves /=
size(level%leaves))
then
1361 deallocate(level%leaves)
1362 allocate(level%leaves(n_leaves))
1367 do i = 1,
size(level%ids)
1370 i_parent = i_parent + 1
1371 level%parents(i_parent) = id
1374 level%leaves(i_leaf) = id
1381 type(
mg_box_t),
intent(in) :: boxes(:)
1382 type(
mg_lvl_t),
intent(inout) :: level
1383 integer,
allocatable :: tmp(:)
1384 integer :: i, id, nb, nb_id, ix
1386 if (
allocated(level%ref_bnds))
deallocate(level%ref_bnds)
1388 if (
size(level%parents) == 0)
then
1390 allocate(level%ref_bnds(0))
1392 allocate(tmp(
size(level%leaves)))
1394 do i = 1,
size(level%leaves)
1395 id = level%leaves(i)
1398 nb_id = boxes(id)%neighbors(nb)
1409 allocate(level%ref_bnds(ix))
1410 level%ref_bnds(:) = tmp(1:ix)
1415 type(
mg_t),
intent(inout) :: mg
1416 integer,
intent(in) :: id
1417 integer :: lvl, i, nb, child_nb(2**(1-1))
1419 integer :: c_id, c_ix_base(1)
1422 error stop
"mg_add_children: not enough space"
1427 mg%boxes(id)%children = c_ids
1428 c_ix_base = 2 * mg%boxes(id)%ix - 1
1429 lvl = mg%boxes(id)%lvl+1
1433 mg%boxes(c_id)%rank = mg%boxes(id)%rank
1435 mg%boxes(c_id)%lvl = lvl
1436 mg%boxes(c_id)%parent = id
1439 mg%boxes(c_id)%r_min = mg%boxes(id)%r_min + &
1441 mg%boxes(c_id)%dr(:) = mg%dr(:, lvl)
1446 if (mg%boxes(id)%neighbors(nb) <
mg_no_box)
then
1448 mg%boxes(child_nb)%neighbors(nb) = mg%boxes(id)%neighbors(nb)
1453 subroutine add_single_child(mg, id, n_boxes_lvl)
1454 type(
mg_t),
intent(inout) :: mg
1455 integer,
intent(in) :: id
1456 integer,
intent(in) :: n_boxes_lvl
1457 integer :: lvl, c_id
1459 c_id = mg%n_boxes + 1
1460 mg%n_boxes = mg%n_boxes + 1
1461 mg%boxes(id)%children(1) = c_id
1462 lvl = mg%boxes(id)%lvl+1
1464 mg%boxes(c_id)%rank = mg%boxes(id)%rank
1465 mg%boxes(c_id)%ix = mg%boxes(id)%ix
1466 mg%boxes(c_id)%lvl = lvl
1467 mg%boxes(c_id)%parent = id
1470 mg%boxes(c_id)%neighbors = mg%boxes(id)%neighbors
1472 mg%boxes(c_id)%neighbors = mg%boxes(id)%neighbors + n_boxes_lvl
1474 mg%boxes(c_id)%r_min = mg%boxes(id)%r_min
1475 mg%boxes(c_id)%dr(:) = mg%dr(:, lvl)
1477 end subroutine add_single_child
1487 type(
mg_t),
intent(inout) :: mg
1488 integer :: i, id, lvl, single_cpu_lvl
1489 integer :: work_left, my_work, i_cpu
1493 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1495 do lvl = mg%lowest_lvl, single_cpu_lvl
1496 do i = 1,
size(mg%lvls(lvl)%ids)
1497 id = mg%lvls(lvl)%ids(i)
1498 mg%boxes(id)%rank = 0
1504 do lvl = single_cpu_lvl+1, mg%highest_lvl
1505 work_left =
size(mg%lvls(lvl)%ids)
1509 do i = 1,
size(mg%lvls(lvl)%ids)
1510 if ((mg%n_cpu - i_cpu - 1) * my_work >= work_left)
then
1515 my_work = my_work + 1
1516 work_left = work_left - 1
1518 id = mg%lvls(lvl)%ids(i)
1519 mg%boxes(id)%rank = i_cpu
1523 do lvl = mg%lowest_lvl, mg%highest_lvl
1524 call update_lvl_info(mg, mg%lvls(lvl))
1536 type(
mg_t),
intent(inout) :: mg
1537 integer :: i, id, lvl, single_cpu_lvl
1538 integer :: work_left, my_work(0:mg%n_cpu), i_cpu
1541 integer :: coarse_rank
1545 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1549 do lvl = mg%highest_lvl, single_cpu_lvl+1, -1
1553 do i = 1,
size(mg%lvls(lvl)%parents)
1554 id = mg%lvls(lvl)%parents(i)
1556 c_ids = mg%boxes(id)%children
1557 c_ranks = mg%boxes(c_ids)%rank
1558 i_cpu = most_popular(c_ranks, my_work, mg%n_cpu)
1559 mg%boxes(id)%rank = i_cpu
1560 my_work(i_cpu) = my_work(i_cpu) + 1
1563 work_left =
size(mg%lvls(lvl)%leaves)
1566 do i = 1,
size(mg%lvls(lvl)%leaves)
1568 if ((mg%n_cpu - i_cpu - 1) * my_work(i_cpu) >= &
1569 work_left + sum(my_work(i_cpu+1:)))
then
1573 my_work(i_cpu) = my_work(i_cpu) + 1
1574 work_left = work_left - 1
1576 id = mg%lvls(lvl)%leaves(i)
1577 mg%boxes(id)%rank = i_cpu
1582 if (single_cpu_lvl < mg%highest_lvl)
then
1583 coarse_rank = most_popular(mg%boxes(&
1584 mg%lvls(single_cpu_lvl+1)%ids)%rank, my_work, mg%n_cpu)
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 = coarse_rank
1596 do lvl = mg%lowest_lvl, mg%highest_lvl
1597 call update_lvl_info(mg, mg%lvls(lvl))
1605 type(
mg_t),
intent(inout) :: mg
1606 integer :: i, id, lvl, j, jd
1609 integer :: single_cpu_lvl, coarse_rank
1610 integer :: my_work(0:mg%n_cpu), i_cpu
1611 integer,
allocatable :: rank_array(:)
1615 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1617 do lvl = mg%highest_lvl-1, single_cpu_lvl+1, -1
1621 do i = 1,
size(mg%lvls(lvl)%leaves)
1622 id = mg%lvls(lvl)%leaves(i)
1623 i_cpu = mg%boxes(id)%rank
1624 my_work(i_cpu) = my_work(i_cpu) + 1
1627 do i = 1,
size(mg%lvls(lvl)%parents)
1628 id = mg%lvls(lvl)%parents(i)
1630 c_ids = mg%boxes(id)%children
1631 c_ranks = mg%boxes(c_ids)%rank
1632 i_cpu = most_popular(c_ranks, my_work, mg%n_cpu)
1633 mg%boxes(id)%rank = i_cpu
1634 my_work(i_cpu) = my_work(i_cpu) + 1
1649 if(single_cpu_lvl < mg%highest_lvl)
then
1650 allocate(rank_array(
size(mg%lvls(single_cpu_lvl+1)%ids)))
1651 do j = 1,
size(mg%lvls(single_cpu_lvl+1)%ids)
1652 jd = mg%lvls(single_cpu_lvl+1)%ids(j)
1653 rank_array(j) = mg%boxes(jd)%rank
1655 coarse_rank = most_popular(rank_array, my_work, mg%n_cpu)
1656 deallocate(rank_array)
1661 do lvl = mg%lowest_lvl, single_cpu_lvl
1662 do i = 1,
size(mg%lvls(lvl)%ids)
1663 id = mg%lvls(lvl)%ids(i)
1664 mg%boxes(id)%rank = coarse_rank
1668 do lvl = mg%lowest_lvl, mg%highest_lvl
1669 call update_lvl_info(mg, mg%lvls(lvl))
1676 pure integer function most_popular(list, work, n_cpu)
1677 integer,
intent(in) :: list(:)
1678 integer,
intent(in) :: n_cpu
1679 integer,
intent(in) :: work(0:n_cpu-1)
1680 integer :: i, best_count, current_count
1681 integer :: my_work, best_work
1687 do i = 1,
size(list)
1688 current_count = count(list == list(i))
1689 my_work = work(list(i))
1692 if (current_count > best_count .or. &
1693 current_count == best_count .and. my_work < best_work)
then
1694 best_count = current_count
1696 most_popular = list(i)
1700 end function most_popular
1702 subroutine update_lvl_info(mg, lvl)
1703 type(
mg_t),
intent(inout) :: mg
1704 type(
mg_lvl_t),
intent(inout) :: lvl
1706 lvl%my_ids = pack(lvl%ids, &
1707 mg%boxes(lvl%ids)%rank == mg%my_rank)
1708 lvl%my_leaves = pack(lvl%leaves, &
1709 mg%boxes(lvl%leaves)%rank == mg%my_rank)
1710 lvl%my_parents = pack(lvl%parents, &
1711 mg%boxes(lvl%parents)%rank == mg%my_rank)
1712 lvl%my_ref_bnds = pack(lvl%ref_bnds, &
1713 mg%boxes(lvl%ref_bnds)%rank == mg%my_rank)
1714 end subroutine update_lvl_info
1719 type(
mg_t),
intent(inout) :: mg
1721 if (mg%n_extra_vars == 0 .and. mg%is_allocated)
then
1722 error stop
"vlaplacian_set_methods: mg%n_extra_vars == 0"
1724 mg%n_extra_vars = max(1, mg%n_extra_vars)
1727 mg%subtract_mean = .false.
1732 mg%bc(:,
mg_iveps)%bc_value = 0.0_dp
1734 select case (mg%geometry_type)
1736 mg%box_op => box_vlpl
1741 select case (mg%smoother_type)
1743 mg%box_smoother => box_gs_vlpl
1745 error stop
"vlaplacian_set_methods: unsupported smoother type"
1749 error stop
"vlaplacian_set_methods: unsupported geometry"
1755 subroutine box_gs_vlpl(mg, id, nc, redblack_cntr)
1756 type(
mg_t),
intent(inout) :: mg
1757 integer,
intent(in) :: id
1758 integer,
intent(in) :: nc
1759 integer,
intent(in) :: redblack_cntr
1760 integer :: i, i0, di
1762 real(
dp) :: idr2(2*1), u(2*1)
1763 real(
dp) :: a0, a(2*1), c(2*1)
1766 idr2(1:2*1:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1767 idr2(2:2*1:2) = idr2(1:2*1:2)
1779 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
1781 if (redblack) i0 = 2 - iand(redblack_cntr, 1)
1785 u(1:2) = cc(i-1:i+1:2, n)
1786 a(1:2) = cc(i-1:i+1:2, i_eps)
1787 c(:) = 2 * a0 * a(:) / (a0 + a(:)) * idr2
1790 (sum(c(:) * u(:)) - cc(i,
mg_irhs)) / sum(c(:))
1793 end subroutine box_gs_vlpl
1796 subroutine box_vlpl(mg, id, nc, i_out)
1797 type(
mg_t),
intent(inout) :: mg
1798 integer,
intent(in) :: id
1799 integer,
intent(in) :: nc
1800 integer,
intent(in) :: i_out
1802 real(
dp) :: idr2(2*1), a0, u0, u(2*1), a(2*1)
1805 idr2(1:2*1:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1806 idr2(2:2*1:2) = idr2(1:2*1:2)
1808 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
1812 a(1:2) = cc(i-1:i+1:2, i_eps)
1814 u(1:2) = cc(i-1:i+1:2, n)
1816 cc(i, i_out) = sum(2 * idr2 * &
1817 a0*a(:)/(a0 + a(:)) * (u(:) - u0))
1820 end subroutine box_vlpl
1928 type(
mg_t),
intent(inout) :: mg
1930 integer,
intent(in),
optional :: comm
1932 logical :: initialized
1934 call mpi_initialized(initialized, ierr)
1935 if (.not. initialized)
then
1939 if (
present(comm))
then
1942 mg%comm = mpi_comm_world
1945 call mpi_comm_rank(mg%comm, mg%my_rank, ierr)
1946 call mpi_comm_size(mg%comm, mg%n_cpu, ierr)
1951 type(
mg_t),
intent(inout) :: mg
1952 integer,
intent(in) :: dsize
1953 integer :: i, n_send, n_recv
1954 integer :: send_req(mg%n_cpu)
1955 integer :: recv_req(mg%n_cpu)
1961 do i = 0, mg%n_cpu - 1
1962 if (mg%buf(i)%i_send > 0)
then
1964 call sort_sendbuf(mg%buf(i), dsize)
1965 call mpi_isend(mg%buf(i)%send, mg%buf(i)%i_send, mpi_double, &
1966 i, 0, mg%comm, send_req(n_send), ierr)
1968 if (mg%buf(i)%i_recv > 0)
then
1970 call mpi_irecv(mg%buf(i)%recv, mg%buf(i)%i_recv, mpi_double, &
1971 i, 0, mg%comm, recv_req(n_recv), ierr)
1975 call mpi_waitall(n_recv, recv_req(1:n_recv), mpi_statuses_ignore, ierr)
1976 call mpi_waitall(n_send, send_req(1:n_send), mpi_statuses_ignore, ierr)
1981 subroutine sort_sendbuf(gc, dsize)
1983 type(
mg_buf_t),
intent(inout) :: gc
1984 integer,
intent(in) :: dsize
1985 integer :: ix_sort(gc%i_ix)
1986 real(
dp) :: buf_cpy(gc%i_send)
1989 call mrgrnk(gc%ix(1:gc%i_ix), ix_sort)
1991 buf_cpy = gc%send(1:gc%i_send)
1995 j = (ix_sort(n)-1) * dsize
1996 gc%send(i+1:i+dsize) = buf_cpy(j+1:j+dsize)
1998 gc%ix(1:gc%i_ix) = gc%ix(ix_sort)
2000 end subroutine sort_sendbuf
2006 type(
mg_t),
intent(inout) :: mg
2007 integer,
intent(out) :: n_send(0:mg%n_cpu-1)
2008 integer,
intent(out) :: n_recv(0:mg%n_cpu-1)
2009 integer,
intent(out) :: dsize
2010 integer :: i, id, lvl, nc
2012 allocate(mg%comm_ghostcell%n_send(0:mg%n_cpu-1, &
2013 mg%first_normal_lvl:mg%highest_lvl))
2014 allocate(mg%comm_ghostcell%n_recv(0:mg%n_cpu-1, &
2015 mg%first_normal_lvl:mg%highest_lvl))
2017 dsize = mg%box_size**(1-1)
2019 do lvl = mg%first_normal_lvl, mg%highest_lvl
2020 nc = mg%box_size_lvl(lvl)
2021 mg%buf(:)%i_send = 0
2022 mg%buf(:)%i_recv = 0
2025 do i = 1,
size(mg%lvls(lvl)%my_ids)
2026 id = mg%lvls(lvl)%my_ids(i)
2027 call buffer_ghost_cells(mg, id, nc, 1, dry_run=.true.)
2031 do i = 1,
size(mg%lvls(lvl-1)%my_ref_bnds)
2032 id = mg%lvls(lvl-1)%my_ref_bnds(i)
2033 call buffer_refinement_boundaries(mg, id, nc, 1, dry_run=.true.)
2038 mg%buf(:)%i_recv = 0
2039 do i = 1,
size(mg%lvls(lvl)%my_ids)
2040 id = mg%lvls(lvl)%my_ids(i)
2041 call set_ghost_cells(mg, id, nc, 1, dry_run=.true.)
2044 mg%comm_ghostcell%n_send(:, lvl) = mg%buf(:)%i_send/dsize
2045 mg%comm_ghostcell%n_recv(:, lvl) = mg%buf(:)%i_recv/dsize
2048 n_send = maxval(mg%comm_ghostcell%n_send, dim=2)
2049 n_recv = maxval(mg%comm_ghostcell%n_recv, dim=2)
2055 type(
mg_t),
intent(inout) :: mg
2060 do lvl = mg%lowest_lvl, mg%highest_lvl
2061 nc = mg%box_size_lvl(lvl)
2062 call mg_phi_bc_store_lvl(mg, lvl, nc)
2065 mg%phi_bc_data_stored = .true.
2068 subroutine mg_phi_bc_store_lvl(mg, lvl, nc)
2069 type(
mg_t),
intent(inout) :: mg
2070 integer,
intent(in) :: lvl
2071 integer,
intent(in) :: nc
2073 integer :: i, id, nb, nb_id, bc_type
2075 do i = 1,
size(mg%lvls(lvl)%my_ids)
2076 id = mg%lvls(lvl)%my_ids(i)
2078 nb_id = mg%boxes(id)%neighbors(nb)
2081 if (
associated(mg%bc(nb,
mg_iphi)%boundary_cond))
then
2082 call mg%bc(nb,
mg_iphi)%boundary_cond(mg%boxes(id), nc, &
2085 bc_type = mg%bc(nb,
mg_iphi)%bc_type
2086 bc = mg%bc(nb,
mg_iphi)%bc_value
2092 mg%boxes(id)%neighbors(nb) = bc_type
2095 call box_set_gc(mg%boxes(id), nb, nc,
mg_irhs, bc)
2099 end subroutine mg_phi_bc_store_lvl
2104 integer,
intent(in) :: iv
2107 do lvl = mg%lowest_lvl, mg%highest_lvl
2116 integer,
intent(in) :: lvl
2117 integer,
intent(in) :: iv
2118 integer :: i, id, dsize, nc
2120 if (lvl < mg%lowest_lvl) &
2121 error stop
"fill_ghost_cells_lvl: lvl < lowest_lvl"
2122 if (lvl > mg%highest_lvl) &
2123 error stop
"fill_ghost_cells_lvl: lvl > highest_lvl"
2125 nc = mg%box_size_lvl(lvl)
2127 if (lvl >= mg%first_normal_lvl)
then
2129 mg%buf(:)%i_send = 0
2130 mg%buf(:)%i_recv = 0
2133 do i = 1,
size(mg%lvls(lvl)%my_ids)
2134 id = mg%lvls(lvl)%my_ids(i)
2135 call buffer_ghost_cells(mg, id, nc, iv, .false.)
2139 do i = 1,
size(mg%lvls(lvl-1)%my_ref_bnds)
2140 id = mg%lvls(lvl-1)%my_ref_bnds(i)
2141 call buffer_refinement_boundaries(mg, id, nc, iv, .false.)
2146 mg%buf(:)%i_recv = mg%comm_ghostcell%n_recv(:, lvl) * dsize
2150 mg%buf(:)%i_recv = 0
2153 do i = 1,
size(mg%lvls(lvl)%my_ids)
2154 id = mg%lvls(lvl)%my_ids(i)
2155 call set_ghost_cells(mg, id, nc, iv, .false.)
2159 subroutine buffer_ghost_cells(mg, id, nc, iv, dry_run)
2160 type(
mg_t),
intent(inout) :: mg
2161 integer,
intent(in) :: id
2162 integer,
intent(in) :: nc
2163 integer,
intent(in) :: iv
2164 logical,
intent(in) :: dry_run
2165 integer :: nb, nb_id, nb_rank
2168 nb_id = mg%boxes(id)%neighbors(nb)
2172 nb_rank = mg%boxes(nb_id)%rank
2174 if (nb_rank /= mg%my_rank)
then
2175 call buffer_for_nb(mg, mg%boxes(id), nc, iv, nb_id, nb_rank, &
2180 end subroutine buffer_ghost_cells
2182 subroutine buffer_refinement_boundaries(mg, id, nc, iv, dry_run)
2183 type(
mg_t),
intent(inout) :: mg
2184 integer,
intent(in) :: id
2185 integer,
intent(in) :: nc
2186 integer,
intent(in) :: iv
2187 logical,
intent(in) :: dry_run
2188 integer :: nb, nb_id, c_ids(2**(1-1))
2189 integer :: n, c_id, c_rank
2192 nb_id = mg%boxes(id)%neighbors(nb)
2195 c_ids = mg%boxes(nb_id)%children(&
2200 c_rank = mg%boxes(c_id)%rank
2202 if (c_rank /= mg%my_rank)
then
2204 call buffer_for_fine_nb(mg, mg%boxes(id), nc, iv, c_id, &
2205 c_rank, nb, dry_run)
2211 end subroutine buffer_refinement_boundaries
2214 subroutine set_ghost_cells(mg, id, nc, iv, dry_run)
2215 type(
mg_t),
intent(inout) :: mg
2216 integer,
intent(in) :: id
2217 integer,
intent(in) :: nc
2218 integer,
intent(in) :: iv
2219 logical,
intent(in) :: dry_run
2221 integer :: nb, nb_id, nb_rank, bc_type
2224 nb_id = mg%boxes(id)%neighbors(nb)
2228 nb_rank = mg%boxes(nb_id)%rank
2230 if (nb_rank /= mg%my_rank)
then
2231 call fill_buffered_nb(mg, mg%boxes(id), nb_rank, &
2232 nb, nc, iv, dry_run)
2233 else if (.not. dry_run)
then
2234 call copy_from_nb(mg%boxes(id), mg%boxes(nb_id), &
2239 call fill_refinement_bnd(mg, id, nb, nc, iv, dry_run)
2240 else if (.not. dry_run)
then
2242 if (mg%phi_bc_data_stored .and. iv ==
mg_iphi)
then
2245 call box_get_gc(mg%boxes(id), nb, nc,
mg_irhs, bc)
2248 if (
associated(mg%bc(nb, iv)%boundary_cond))
then
2249 call mg%bc(nb, iv)%boundary_cond(mg%boxes(id), nc, iv, &
2252 bc_type = mg%bc(nb, iv)%bc_type
2253 bc = mg%bc(nb, iv)%bc_value
2257 call box_set_gc(mg%boxes(id), nb, nc, iv, bc)
2258 call bc_to_gc(mg, id, nc, iv, nb, bc_type)
2261 end subroutine set_ghost_cells
2263 subroutine fill_refinement_bnd(mg, id, nb, nc, iv, dry_run)
2264 type(
mg_t),
intent(inout) :: mg
2265 integer,
intent(in) :: id
2266 integer,
intent(in) :: nc
2267 integer,
intent(in) :: iv
2268 integer,
intent(in) :: nb
2269 logical,
intent(in) :: dry_run
2271 integer :: p_id, p_nb_id, ix_offset(1)
2272 integer :: i, dsize, p_nb_rank
2275 p_id = mg%boxes(id)%parent
2276 p_nb_id = mg%boxes(p_id)%neighbors(nb)
2277 p_nb_rank = mg%boxes(p_nb_id)%rank
2279 if (p_nb_rank /= mg%my_rank)
then
2280 i = mg%buf(p_nb_rank)%i_recv
2281 if (.not. dry_run)
then
2282 gc = reshape(mg%buf(p_nb_rank)%recv(i+1:i+dsize), shape(gc))
2284 mg%buf(p_nb_rank)%i_recv = mg%buf(p_nb_rank)%i_recv + dsize
2285 else if (.not. dry_run)
then
2287 call box_gc_for_fine_neighbor(mg%boxes(p_nb_id),
mg_neighb_rev(nb), &
2288 ix_offset, nc, iv, gc)
2291 if (.not. dry_run)
then
2292 if (
associated(mg%bc(nb, iv)%refinement_bnd))
then
2293 call mg%bc(nb, iv)%refinement_bnd(mg%boxes(id), nc, iv, nb, gc)
2295 call sides_rb(mg%boxes(id), nc, iv, nb, gc)
2298 end subroutine fill_refinement_bnd
2300 subroutine copy_from_nb(box, box_nb, nb, nc, iv)
2301 type(
mg_box_t),
intent(inout) :: box
2302 type(
mg_box_t),
intent(in) :: box_nb
2303 integer,
intent(in) :: nb
2304 integer,
intent(in) :: nc
2305 integer,
intent(in) :: iv
2308 call box_gc_for_neighbor(box_nb,
mg_neighb_rev(nb), nc, iv, gc)
2309 call box_set_gc(box, nb, nc, iv, gc)
2310 end subroutine copy_from_nb
2312 subroutine buffer_for_nb(mg, box, nc, iv, nb_id, nb_rank, nb, dry_run)
2313 type(
mg_t),
intent(inout) :: mg
2314 type(
mg_box_t),
intent(inout) :: box
2315 integer,
intent(in) :: nc
2316 integer,
intent(in) :: iv
2317 integer,
intent(in) :: nb_id
2318 integer,
intent(in) :: nb_rank
2319 integer,
intent(in) :: nb
2320 logical,
intent(in) :: dry_run
2324 i = mg%buf(nb_rank)%i_send
2327 if (.not. dry_run)
then
2328 call box_gc_for_neighbor(box, nb, nc, iv, gc)
2329 mg%buf(nb_rank)%send(i+1:i+dsize) = pack(gc, .true.)
2334 i = mg%buf(nb_rank)%i_ix
2335 if (.not. dry_run)
then
2339 mg%buf(nb_rank)%i_send = mg%buf(nb_rank)%i_send + dsize
2340 mg%buf(nb_rank)%i_ix = mg%buf(nb_rank)%i_ix + 1
2341 end subroutine buffer_for_nb
2343 subroutine buffer_for_fine_nb(mg, box, nc, iv, fine_id, fine_rank, nb, dry_run)
2344 type(
mg_t),
intent(inout) :: mg
2345 type(
mg_box_t),
intent(inout) :: box
2346 integer,
intent(in) :: nc
2347 integer,
intent(in) :: iv
2348 integer,
intent(in) :: fine_id
2349 integer,
intent(in) :: fine_rank
2350 integer,
intent(in) :: nb
2351 logical,
intent(in) :: dry_run
2352 integer :: i, dsize, ix_offset(1)
2355 i = mg%buf(fine_rank)%i_send
2358 if (.not. dry_run)
then
2360 call box_gc_for_fine_neighbor(box, nb, ix_offset, nc, iv, gc)
2361 mg%buf(fine_rank)%send(i+1:i+dsize) = pack(gc, .true.)
2366 i = mg%buf(fine_rank)%i_ix
2367 if (.not. dry_run)
then
2372 mg%buf(fine_rank)%i_send = mg%buf(fine_rank)%i_send + dsize
2373 mg%buf(fine_rank)%i_ix = mg%buf(fine_rank)%i_ix + 1
2374 end subroutine buffer_for_fine_nb
2376 subroutine fill_buffered_nb(mg, box, nb_rank, nb, nc, iv, dry_run)
2377 type(
mg_t),
intent(inout) :: mg
2378 type(
mg_box_t),
intent(inout) :: box
2379 integer,
intent(in) :: nb_rank
2380 integer,
intent(in) :: nb
2381 integer,
intent(in) :: nc
2382 integer,
intent(in) :: iv
2383 logical,
intent(in) :: dry_run
2387 i = mg%buf(nb_rank)%i_recv
2390 if (.not. dry_run)
then
2391 gc = mg%buf(nb_rank)%recv(i+1)
2392 call box_set_gc(box, nb, nc, iv, gc)
2394 mg%buf(nb_rank)%i_recv = mg%buf(nb_rank)%i_recv + dsize
2396 end subroutine fill_buffered_nb
2398 subroutine box_gc_for_neighbor(box, nb, nc, iv, gc)
2400 integer,
intent(in) :: nb, nc, iv
2401 real(
dp),
intent(out) :: gc(1)
2409 end subroutine box_gc_for_neighbor
2412 subroutine box_gc_for_fine_neighbor(box, nb, di, nc, iv, gc)
2414 integer,
intent(in) :: nb
2415 integer,
intent(in) :: di(1)
2416 integer,
intent(in) :: nc, iv
2417 real(
dp),
intent(out) :: gc(1)
2428 tmp = box%cc(nc, iv)
2436 end subroutine box_gc_for_fine_neighbor
2438 subroutine box_get_gc(box, nb, nc, iv, gc)
2440 integer,
intent(in) :: nb, nc, iv
2441 real(
dp),
intent(out) :: gc(1)
2447 gc = box%cc(nc+1, iv)
2449 end subroutine box_get_gc
2451 subroutine box_set_gc(box, nb, nc, iv, gc)
2452 type(
mg_box_t),
intent(inout) :: box
2453 integer,
intent(in) :: nb, nc, iv
2454 real(
dp),
intent(in) :: gc(1)
2458 box%cc(0, iv) = gc(1)
2460 box%cc(nc+1, iv) = gc(1)
2462 end subroutine box_set_gc
2464 subroutine bc_to_gc(mg, id, nc, iv, nb, bc_type)
2465 type(
mg_t),
intent(inout) :: mg
2466 integer,
intent(in) :: id
2467 integer,
intent(in) :: nc
2468 integer,
intent(in) :: iv
2469 integer,
intent(in) :: nb
2470 integer,
intent(in) :: bc_type
2471 real(
dp) :: c0, c1, c2, dr
2481 select case (bc_type)
2496 error stop
"bc_to_gc: unknown boundary condition"
2501 mg%boxes(id)%cc(0, iv) = &
2502 c0 * mg%boxes(id)%cc(0, iv) + &
2503 c1 * mg%boxes(id)%cc(1, iv) + &
2504 c2 * mg%boxes(id)%cc(2, iv)
2506 mg%boxes(id)%cc(nc+1, iv) = &
2507 c0 * mg%boxes(id)%cc(nc+1, iv) + &
2508 c1 * mg%boxes(id)%cc(nc, iv) + &
2509 c2 * mg%boxes(id)%cc(nc-1, iv)
2511 end subroutine bc_to_gc
2514 subroutine sides_rb(box, nc, iv, nb, gc)
2515 type(
mg_box_t),
intent(inout) :: box
2516 integer,
intent(in) :: nc
2517 integer,
intent(in) :: iv
2518 integer,
intent(in) :: nb
2520 real(
dp),
intent(in) :: gc(1)
2522 integer :: i, ix, dix
2536 box%cc(i-di, iv) = (2 * gc(1) + box%cc(i, iv))/3.0_dp
2539 end subroutine sides_rb
2545 type(
mg_t),
intent(inout) :: mg
2546 integer,
intent(out) :: n_send(0:mg%n_cpu-1)
2547 integer,
intent(out) :: n_recv(0:mg%n_cpu-1)
2548 integer,
intent(out) :: dsize
2549 integer :: lvl, min_lvl
2551 if (.not.
allocated(mg%comm_restrict%n_send))
then
2552 error stop
"Call restrict_buffer_size before prolong_buffer_size"
2555 min_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
2556 allocate(mg%comm_prolong%n_send(0:mg%n_cpu-1, &
2557 min_lvl:mg%highest_lvl))
2558 allocate(mg%comm_prolong%n_recv(0:mg%n_cpu-1, &
2559 min_lvl:mg%highest_lvl))
2561 mg%comm_prolong%n_recv(:, mg%highest_lvl) = 0
2562 mg%comm_prolong%n_send(:, mg%highest_lvl) = 0
2564 do lvl = min_lvl, mg%highest_lvl-1
2565 mg%comm_prolong%n_recv(:, lvl) = &
2566 mg%comm_restrict%n_send(:, lvl+1)
2567 mg%comm_prolong%n_send(:, lvl) = &
2568 mg%comm_restrict%n_recv(:, lvl+1)
2573 dsize = (mg%box_size)**1
2574 n_send = maxval(mg%comm_prolong%n_send, dim=2)
2575 n_recv = maxval(mg%comm_prolong%n_recv, dim=2)
2581 type(
mg_t),
intent(inout) :: mg
2582 integer,
intent(in) :: lvl
2583 integer,
intent(in) :: iv
2584 integer,
intent(in) :: iv_to
2586 logical,
intent(in) :: add
2587 integer :: i, id, dsize, nc
2589 if (lvl == mg%highest_lvl) error stop
"cannot prolong highest level"
2590 if (lvl < mg%lowest_lvl) error stop
"cannot prolong below lowest level"
2593 if (lvl >= mg%first_normal_lvl-1)
then
2594 dsize = mg%box_size**1
2595 mg%buf(:)%i_send = 0
2598 do i = 1,
size(mg%lvls(lvl)%my_ids)
2599 id = mg%lvls(lvl)%my_ids(i)
2600 call prolong_set_buffer(mg, id, mg%box_size, iv, method)
2603 mg%buf(:)%i_recv = mg%comm_prolong%n_recv(:, lvl) * dsize
2605 mg%buf(:)%i_recv = 0
2608 nc = mg%box_size_lvl(lvl+1)
2609 do i = 1,
size(mg%lvls(lvl+1)%my_ids)
2610 id = mg%lvls(lvl+1)%my_ids(i)
2611 call prolong_onto(mg, id, nc, iv, iv_to, add, method)
2618 subroutine prolong_set_buffer(mg, id, nc, iv, method)
2619 type(
mg_t),
intent(inout) :: mg
2620 integer,
intent(in) :: id
2621 integer,
intent(in) :: nc
2622 integer,
intent(in) :: iv
2624 integer :: i, dix(1)
2625 integer :: i_c, c_id, c_rank, dsize
2631 c_id = mg%boxes(id)%children(i_c)
2633 c_rank = mg%boxes(c_id)%rank
2634 if (c_rank /= mg%my_rank)
then
2636 call method(mg, id, dix, nc, iv, tmp)
2638 i = mg%buf(c_rank)%i_send
2639 mg%buf(c_rank)%send(i+1:i+dsize) = pack(tmp, .true.)
2640 mg%buf(c_rank)%i_send = mg%buf(c_rank)%i_send + dsize
2642 i = mg%buf(c_rank)%i_ix
2643 mg%buf(c_rank)%ix(i+1) = c_id
2644 mg%buf(c_rank)%i_ix = mg%buf(c_rank)%i_ix + 1
2648 end subroutine prolong_set_buffer
2651 subroutine prolong_onto(mg, id, nc, iv, iv_to, add, method)
2652 type(
mg_t),
intent(inout) :: mg
2653 integer,
intent(in) :: id
2654 integer,
intent(in) :: nc
2655 integer,
intent(in) :: iv
2656 integer,
intent(in) :: iv_to
2657 logical,
intent(in) :: add
2659 integer :: hnc, p_id, p_rank, i, dix(1), dsize
2663 p_id = mg%boxes(id)%parent
2664 p_rank = mg%boxes(p_id)%rank
2666 if (p_rank == mg%my_rank)
then
2668 call method(mg, p_id, dix, nc, iv, tmp)
2671 i = mg%buf(p_rank)%i_recv
2672 tmp = reshape(mg%buf(p_rank)%recv(i+1:i+dsize), [nc])
2673 mg%buf(p_rank)%i_recv = mg%buf(p_rank)%i_recv + dsize
2677 mg%boxes(id)%cc(1:nc, iv_to) = &
2678 mg%boxes(id)%cc(1:nc, iv_to) + tmp
2680 mg%boxes(id)%cc(1:nc, iv_to) = tmp
2683 end subroutine prolong_onto
2687 type(
mg_t),
intent(inout) :: mg
2688 integer,
intent(in) :: p_id
2689 integer,
intent(in) :: dix(1)
2690 integer,
intent(in) :: nc
2691 integer,
intent(in) :: iv
2692 real(
dp),
intent(out) :: fine(nc)
2696 real(
dp) :: f0, flx, fhx
2700 associate(crs => mg%boxes(p_id)%cc)
2704 f0 = 0.75_dp * crs(ic, iv)
2705 flx = 0.25_dp * crs(ic-1, iv)
2706 fhx = 0.25_dp * crs(ic+1, iv)
2708 fine(2*i-1) = f0 + flx
2709 fine(2*i ) = f0 + fhx
2717 type(
mg_t),
intent(inout) :: mg
2719 mg%subtract_mean = .false.
2721 select case (mg%geometry_type)
2723 mg%box_op => box_helmh
2725 select case (mg%smoother_type)
2727 mg%box_smoother => box_gs_helmh
2729 error stop
"helmholtz_set_methods: unsupported smoother type"
2732 error stop
"helmholtz_set_methods: unsupported geometry"
2738 real(
dp),
intent(in) :: lambda
2741 error stop
"helmholtz_set_lambda: lambda < 0 not allowed"
2747 subroutine box_gs_helmh(mg, id, nc, redblack_cntr)
2748 type(
mg_t),
intent(inout) :: mg
2749 integer,
intent(in) :: id
2750 integer,
intent(in) :: nc
2751 integer,
intent(in) :: redblack_cntr
2752 integer :: i, i0, di
2753 real(
dp) :: idr2(1), fac
2756 idr2 = 1/mg%dr(:, mg%boxes(id)%lvl)**2
2769 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
2770 if (redblack) i0 = 2 - iand(redblack_cntr, 1)
2773 cc(i, n) = fac * ( &
2774 idr2(1) * (cc(i+1, n) + cc(i-1, n)) - &
2778 end subroutine box_gs_helmh
2781 subroutine box_helmh(mg, id, nc, i_out)
2782 type(
mg_t),
intent(inout) :: mg
2783 integer,
intent(in) :: id
2784 integer,
intent(in) :: nc
2785 integer,
intent(in) :: i_out
2789 idr2 = 1 / mg%dr(:, mg%boxes(id)%lvl)**2
2791 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi)
2794 idr2(1) * (cc(i-1, n) + cc(i+1, n) - 2 * cc(i, n)) - &
2798 end subroutine box_helmh
2804 type(
mg_t),
intent(inout) :: mg
2809 select case (mg%operator_type)
2821 error stop
"mg_set_methods: unknown operator"
2827 mg%n_smoother_substeps = 2
2829 mg%n_smoother_substeps = 1
2833 subroutine check_methods(mg)
2834 type(
mg_t),
intent(inout) :: mg
2836 if (.not.
associated(mg%box_op) .or. &
2837 .not.
associated(mg%box_smoother))
then
2841 end subroutine check_methods
2843 subroutine mg_add_timers(mg)
2844 type(
mg_t),
intent(inout) :: mg
2845 timer_total_vcycle =
mg_add_timer(mg,
"mg total V-cycle")
2846 timer_total_fmg =
mg_add_timer(mg,
"mg total FMG cycle")
2848 timer_smoother_gc =
mg_add_timer(mg,
"mg smoother g.c.")
2851 timer_update_coarse =
mg_add_timer(mg,
"mg update coarse")
2852 end subroutine mg_add_timers
2856 type(
mg_t),
intent(inout) :: mg
2857 logical,
intent(in) :: have_guess
2858 real(
dp),
intent(out),
optional :: max_res
2859 integer :: lvl, i, id
2861 call check_methods(mg)
2862 if (timer_smoother == -1)
call mg_add_timers(mg)
2866 if (.not. have_guess)
then
2867 do lvl = mg%highest_lvl, mg%lowest_lvl, -1
2868 do i = 1,
size(mg%lvls(lvl)%my_ids)
2869 id = mg%lvls(lvl)%my_ids(i)
2870 mg%boxes(id)%cc(:,
mg_iphi) = 0.0_dp
2878 do lvl = mg%highest_lvl, mg%lowest_lvl+1, -1
2881 call update_coarse(mg, lvl)
2885 if (mg%subtract_mean)
then
2887 call subtract_mean(mg,
mg_irhs, .false.)
2890 do lvl = mg%lowest_lvl, mg%highest_lvl
2892 do i = 1,
size(mg%lvls(lvl)%my_ids)
2893 id = mg%lvls(lvl)%my_ids(i)
2894 mg%boxes(id)%cc(:,
mg_iold) = &
2898 if (lvl > mg%lowest_lvl)
then
2902 call correct_children(mg, lvl-1)
2910 if (lvl == mg%highest_lvl)
then
2923 type(
mg_t),
intent(inout) :: mg
2924 integer,
intent(in),
optional :: highest_lvl
2925 real(
dp),
intent(out),
optional :: max_res
2927 logical,
intent(in),
optional :: standalone
2928 integer :: lvl, min_lvl, i, max_lvl, ierr
2929 real(
dp) :: init_res, res
2930 logical :: is_standalone
2932 is_standalone = .true.
2933 if (
present(standalone)) is_standalone = standalone
2935 call check_methods(mg)
2936 if (timer_smoother == -1)
call mg_add_timers(mg)
2938 if (is_standalone) &
2941 if (mg%subtract_mean .and. .not.
present(highest_lvl))
then
2944 call subtract_mean(mg,
mg_irhs, .false.)
2947 min_lvl = mg%lowest_lvl
2948 max_lvl = mg%highest_lvl
2949 if (
present(highest_lvl)) max_lvl = highest_lvl
2952 if (is_standalone)
then
2956 do lvl = max_lvl, min_lvl+1, -1
2958 call smooth_boxes(mg, lvl, mg%n_cycle_down)
2963 call update_coarse(mg, lvl)
2968 if (.not. all(mg%boxes(mg%lvls(min_lvl)%ids)%rank == &
2969 mg%boxes(mg%lvls(min_lvl)%ids(1))%rank))
then
2970 error stop
"Multiple CPUs for coarse grid (not implemented yet)"
2973 init_res = max_residual_lvl(mg, min_lvl)
2974 do i = 1, mg%max_coarse_cycles
2975 call smooth_boxes(mg, min_lvl, mg%n_cycle_up+mg%n_cycle_down)
2976 res = max_residual_lvl(mg, min_lvl)
2977 if (res < mg%residual_coarse_rel * init_res .or. &
2978 res < mg%residual_coarse_abs)
exit
2983 do lvl = min_lvl+1, max_lvl
2987 call correct_children(mg, lvl-1)
2994 call smooth_boxes(mg, lvl, mg%n_cycle_up)
2997 if (
present(max_res))
then
2999 do lvl = min_lvl, max_lvl
3000 res = max_residual_lvl(mg, lvl)
3001 init_res = max(res, init_res)
3003 call mpi_allreduce(init_res, max_res, 1, &
3004 mpi_double, mpi_max, mg%comm, ierr)
3008 if (mg%subtract_mean)
then
3009 call subtract_mean(mg,
mg_iphi, .true.)
3012 if (is_standalone) &
3016 subroutine subtract_mean(mg, iv, include_ghostcells)
3018 type(
mg_t),
intent(inout) :: mg
3019 integer,
intent(in) :: iv
3020 logical,
intent(in) :: include_ghostcells
3021 integer :: i, id, lvl, nc, ierr
3022 real(
dp) :: sum_iv, mean_iv, volume
3025 sum_iv = get_sum(mg, iv)
3026 call mpi_allreduce(sum_iv, mean_iv, 1, &
3027 mpi_double, mpi_sum, mg%comm, ierr)
3030 volume = nc**1 * product(mg%dr(:, 1)) *
size(mg%lvls(1)%ids)
3031 mean_iv = mean_iv / volume
3033 do lvl = mg%lowest_lvl, mg%highest_lvl
3034 nc = mg%box_size_lvl(lvl)
3036 do i = 1,
size(mg%lvls(lvl)%my_ids)
3037 id = mg%lvls(lvl)%my_ids(i)
3038 if (include_ghostcells)
then
3039 mg%boxes(id)%cc(:, iv) = &
3040 mg%boxes(id)%cc(:, iv) - mean_iv
3042 mg%boxes(id)%cc(1:nc, iv) = &
3043 mg%boxes(id)%cc(1:nc, iv) - mean_iv
3047 end subroutine subtract_mean
3049 real(
dp) function get_sum(mg, iv)
3050 type(
mg_t),
intent(in) :: mg
3051 integer,
intent(in) :: iv
3052 integer :: lvl, i, id, nc
3056 do lvl = 1, mg%highest_lvl
3057 nc = mg%box_size_lvl(lvl)
3058 w = product(mg%dr(:, lvl))
3059 do i = 1,
size(mg%lvls(lvl)%my_leaves)
3060 id = mg%lvls(lvl)%my_leaves(i)
3061 get_sum = get_sum + w * &
3062 sum(mg%boxes(id)%cc(1:nc, iv))
3065 end function get_sum
3067 real(
dp) function max_residual_lvl(mg, lvl)
3068 type(
mg_t),
intent(inout) :: mg
3069 integer,
intent(in) :: lvl
3070 integer :: i, id, nc
3073 nc = mg%box_size_lvl(lvl)
3074 max_residual_lvl = 0.0_dp
3076 do i = 1,
size(mg%lvls(lvl)%my_ids)
3077 id = mg%lvls(lvl)%my_ids(i)
3078 call residual_box(mg, id, nc)
3079 res = maxval(abs(mg%boxes(id)%cc(1:nc,
mg_ires)))
3080 max_residual_lvl = max(max_residual_lvl, res)
3082 end function max_residual_lvl
3118 subroutine update_coarse(mg, lvl)
3119 type(
mg_t),
intent(inout) :: mg
3120 integer,
intent(in) :: lvl
3121 integer :: i, id, nc, nc_c
3123 nc = mg%box_size_lvl(lvl)
3124 nc_c = mg%box_size_lvl(lvl-1)
3127 do i = 1,
size(mg%lvls(lvl)%my_ids)
3128 id = mg%lvls(lvl)%my_ids(i)
3129 call residual_box(mg, id, nc)
3140 do i = 1,
size(mg%lvls(lvl-1)%my_parents)
3141 id = mg%lvls(lvl-1)%my_parents(i)
3144 call mg%box_op(mg, id, nc_c,
mg_irhs)
3147 mg%boxes(id)%cc(1:nc_c,
mg_irhs) = &
3148 mg%boxes(id)%cc(1:nc_c,
mg_irhs) + &
3149 mg%boxes(id)%cc(1:nc_c,
mg_ires)
3152 mg%boxes(id)%cc(:,
mg_iold) = &
3155 end subroutine update_coarse
3158 subroutine correct_children(mg, lvl)
3159 type(
mg_t),
intent(inout) :: mg
3160 integer,
intent(in) :: lvl
3163 do i = 1,
size(mg%lvls(lvl)%my_parents)
3164 id = mg%lvls(lvl)%my_parents(i)
3167 mg%boxes(id)%cc(:,
mg_ires) = &
3168 mg%boxes(id)%cc(:,
mg_iphi) - &
3173 end subroutine correct_children
3175 subroutine smooth_boxes(mg, lvl, n_cycle)
3176 type(
mg_t),
intent(inout) :: mg
3177 integer,
intent(in) :: lvl
3178 integer,
intent(in) :: n_cycle
3179 integer :: n, i, id, nc
3181 nc = mg%box_size_lvl(lvl)
3183 do n = 1, n_cycle * mg%n_smoother_substeps
3185 do i = 1,
size(mg%lvls(lvl)%my_ids)
3186 id = mg%lvls(lvl)%my_ids(i)
3187 call mg%box_smoother(mg, id, nc, n)
3195 end subroutine smooth_boxes
3197 subroutine residual_box(mg, id, nc)
3198 type(
mg_t),
intent(inout) :: mg
3199 integer,
intent(in) :: id
3200 integer,
intent(in) :: nc
3202 call mg%box_op(mg, id, nc,
mg_ires)
3204 mg%boxes(id)%cc(1:nc,
mg_ires) = &
3205 mg%boxes(id)%cc(1:nc,
mg_irhs) &
3206 - mg%boxes(id)%cc(1:nc,
mg_ires)
3207 end subroutine residual_box
3211 type(
mg_t),
intent(inout) :: mg
3212 integer,
intent(in) :: i_out
3214 integer :: lvl, i, id, nc
3216 do lvl = mg%lowest_lvl, mg%highest_lvl
3217 nc = mg%box_size_lvl(lvl)
3218 do i = 1,
size(mg%lvls(lvl)%my_ids)
3219 id = mg%lvls(lvl)%my_ids(i)
3220 if (
present(op))
then
3221 call op(mg, id, nc, i_out)
3223 call mg%box_op(mg, id, nc, i_out)
3233 type(
mg_t),
intent(inout) :: mg
3234 integer,
intent(out) :: n_send(0:mg%n_cpu-1)
3235 integer,
intent(out) :: n_recv(0:mg%n_cpu-1)
3236 integer,
intent(out) :: dsize
3237 integer :: n_out(0:mg%n_cpu-1, mg%first_normal_lvl:mg%highest_lvl)
3238 integer :: n_in(0:mg%n_cpu-1, mg%first_normal_lvl:mg%highest_lvl)
3239 integer :: lvl, i, id, p_id, p_rank
3240 integer :: i_c, c_id, c_rank, min_lvl
3244 min_lvl = max(mg%lowest_lvl+1, mg%first_normal_lvl)
3246 do lvl = min_lvl, mg%highest_lvl
3248 do i = 1,
size(mg%lvls(lvl-1)%my_parents)
3249 id = mg%lvls(lvl-1)%my_parents(i)
3251 c_id = mg%boxes(id)%children(i_c)
3254 c_rank = mg%boxes(c_id)%rank
3255 if (c_rank /= mg%my_rank)
then
3256 n_in(c_rank, lvl) = n_in(c_rank, lvl) + 1
3263 do i = 1,
size(mg%lvls(lvl)%my_ids)
3264 id = mg%lvls(lvl)%my_ids(i)
3266 p_id = mg%boxes(id)%parent
3267 p_rank = mg%boxes(p_id)%rank
3268 if (p_rank /= mg%my_rank)
then
3269 n_out(p_rank, lvl) = n_out(p_rank, lvl) + 1
3275 allocate(mg%comm_restrict%n_send(0:mg%n_cpu-1, &
3276 mg%first_normal_lvl:mg%highest_lvl))
3277 allocate(mg%comm_restrict%n_recv(0:mg%n_cpu-1, &
3278 mg%first_normal_lvl:mg%highest_lvl))
3279 mg%comm_restrict%n_send = n_out
3280 mg%comm_restrict%n_recv = n_in
3282 dsize = (mg%box_size/2)**1
3283 n_send = maxval(n_out, dim=2)
3284 n_recv = maxval(n_in, dim=2)
3289 type(
mg_t),
intent(inout) :: mg
3290 integer,
intent(in) :: iv
3293 do lvl = mg%highest_lvl, mg%lowest_lvl+1, -1
3301 type(
mg_t),
intent(inout) :: mg
3302 integer,
intent(in) :: iv
3303 integer,
intent(in) :: lvl
3304 integer :: i, id, dsize, nc
3306 if (lvl <= mg%lowest_lvl) error stop
"cannot restrict lvl <= lowest_lvl"
3308 nc = mg%box_size_lvl(lvl)
3310 if (lvl >= mg%first_normal_lvl)
then
3313 mg%buf(:)%i_send = 0
3316 do i = 1,
size(mg%lvls(lvl)%my_ids)
3317 id = mg%lvls(lvl)%my_ids(i)
3318 call restrict_set_buffer(mg, id, iv)
3321 mg%buf(:)%i_recv = mg%comm_restrict%n_recv(:, lvl) * dsize
3323 mg%buf(:)%i_recv = 0
3326 do i = 1,
size(mg%lvls(lvl-1)%my_parents)
3327 id = mg%lvls(lvl-1)%my_parents(i)
3328 call restrict_onto(mg, id, nc, iv)
3332 subroutine restrict_set_buffer(mg, id, iv)
3333 type(
mg_t),
intent(inout) :: mg
3334 integer,
intent(in) :: id
3335 integer,
intent(in) :: iv
3336 integer :: i, n, hnc, p_id, p_rank
3337 real(
dp) :: tmp(mg%box_size/2)
3340 p_id = mg%boxes(id)%parent
3341 p_rank = mg%boxes(p_id)%rank
3343 if (p_rank /= mg%my_rank)
then
3346 sum(mg%boxes(id)%cc(2*i-1:2*i, iv))
3351 i = mg%buf(p_rank)%i_send
3352 mg%buf(p_rank)%send(i+1:i+n) = pack(tmp, .true.)
3353 mg%buf(p_rank)%i_send = mg%buf(p_rank)%i_send + n
3356 i = mg%buf(p_rank)%i_ix
3359 mg%buf(p_rank)%i_ix = mg%buf(p_rank)%i_ix + 1
3361 end subroutine restrict_set_buffer
3363 subroutine restrict_onto(mg, id, nc, iv)
3364 type(
mg_t),
intent(inout) :: mg
3365 integer,
intent(in) :: id
3366 integer,
intent(in) :: nc
3367 integer,
intent(in) :: iv
3368 integer :: i, hnc, dsize, i_c, c_id
3369 integer :: c_rank, dix(1)
3375 c_id = mg%boxes(id)%children(i_c)
3377 c_rank = mg%boxes(c_id)%rank
3380 if (c_rank == mg%my_rank)
then
3382 mg%boxes(id)%cc(dix(1)+i, iv) = 0.5_dp * &
3383 sum(mg%boxes(c_id)%cc(2*i-1:2*i, iv))
3386 i = mg%buf(c_rank)%i_recv
3387 mg%boxes(id)%cc(dix(1)+1:dix(1)+hnc, iv) = &
3388 reshape(mg%buf(c_rank)%recv(i+1:i+dsize), [hnc])
3389 mg%buf(c_rank)%i_recv = mg%buf(c_rank)%i_recv + dsize
3393 end subroutine restrict_onto
3397 Subroutine i_mrgrnk (XDONT, IRNGT)
3404 Integer,
Dimension (:),
Intent (In) :: XDONT
3405 Integer,
Dimension (:),
Intent (Out) :: IRNGT
3407 Integer :: XVALA, XVALB
3409 Integer,
Dimension (SIZE(IRNGT)) :: JWRKT
3410 Integer :: LMTNA, LMTNC, IRNG1, IRNG2
3411 Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
3413 nval = min(
SIZE(xdont),
SIZE(irngt))
3426 Do iind = 2, nval, 2
3427 If (xdont(iind-1) <= xdont(iind))
Then
3428 irngt(iind-1) = iind - 1
3431 irngt(iind-1) = iind
3432 irngt(iind) = iind - 1
3435 If (modulo(nval, 2) /= 0)
Then
3452 Do iwrkd = 0, nval - 1, 4
3453 If ((iwrkd+4) > nval)
Then
3454 If ((iwrkd+2) >= nval)
Exit
3458 If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3)))
Exit
3462 If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3)))
Then
3463 irng2 = irngt(iwrkd+2)
3464 irngt(iwrkd+2) = irngt(iwrkd+3)
3465 irngt(iwrkd+3) = irng2
3470 irng1 = irngt(iwrkd+1)
3471 irngt(iwrkd+1) = irngt(iwrkd+3)
3472 irngt(iwrkd+3) = irngt(iwrkd+2)
3473 irngt(iwrkd+2) = irng1
3480 If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) cycle
3484 If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3)))
Then
3485 irng2 = irngt(iwrkd+2)
3486 irngt(iwrkd+2) = irngt(iwrkd+3)
3487 If (xdont(irng2) <= xdont(irngt(iwrkd+4)))
Then
3489 irngt(iwrkd+3) = irng2
3492 irngt(iwrkd+3) = irngt(iwrkd+4)
3493 irngt(iwrkd+4) = irng2
3499 irng1 = irngt(iwrkd+1)
3500 irng2 = irngt(iwrkd+2)
3501 irngt(iwrkd+1) = irngt(iwrkd+3)
3502 If (xdont(irng1) <= xdont(irngt(iwrkd+4)))
Then
3503 irngt(iwrkd+2) = irng1
3504 If (xdont(irng2) <= xdont(irngt(iwrkd+4)))
Then
3506 irngt(iwrkd+3) = irng2
3509 irngt(iwrkd+3) = irngt(iwrkd+4)
3510 irngt(iwrkd+4) = irng2
3514 irngt(iwrkd+2) = irngt(iwrkd+4)
3515 irngt(iwrkd+3) = irng1
3516 irngt(iwrkd+4) = irng2
3531 If (lmtna >= nval)
Exit
3540 jinda = iwrkf + lmtna
3541 iwrkf = iwrkf + lmtnc
3542 If (iwrkf >= nval)
Then
3543 If (jinda >= nval)
Exit
3559 jwrkt(1:lmtna) = irngt(iwrkd:jinda)
3561 xvala = xdont(jwrkt(iinda))
3562 xvalb = xdont(irngt(iindb))
3569 If (xvala > xvalb)
Then
3570 irngt(iwrk) = irngt(iindb)
3572 If (iindb > iwrkf)
Then
3574 irngt(iwrk+1:iwrkf) = jwrkt(iinda:lmtna)
3577 xvalb = xdont(irngt(iindb))
3579 irngt(iwrk) = jwrkt(iinda)
3581 If (iinda > lmtna) exit
3582 xvala = xdont(jwrkt(iinda))
3595 End Subroutine i_mrgrnk
3600 type(
mg_t),
intent(inout) :: mg
3602 if (mg%n_extra_vars == 0 .and. mg%is_allocated)
then
3603 error stop
"ahelmholtz_set_methods: mg%n_extra_vars == 0"
3605 mg%n_extra_vars = max(1, mg%n_extra_vars)
3615 select case (mg%geometry_type)
3617 mg%box_op => box_ahelmh
3619 select case (mg%smoother_type)
3621 mg%box_smoother => box_gs_ahelmh
3623 error stop
"ahelmholtz_set_methods: unsupported smoother type"
3626 error stop
"ahelmholtz_set_methods: unsupported geometry"
3632 real(
dp),
intent(in) :: lambda
3635 error stop
"ahelmholtz_set_lambda: lambda < 0 not allowed"
3641 subroutine box_gs_ahelmh(mg, id, nc, redblack_cntr)
3642 type(
mg_t),
intent(inout) :: mg
3643 integer,
intent(in) :: id
3644 integer,
intent(in) :: nc
3645 integer,
intent(in) :: redblack_cntr
3646 integer :: i, i0, di
3648 real(
dp) :: idr2(2*1), u(2*1)
3649 real(
dp) :: a0(2*1), a(2*1), c(2*1)
3652 idr2(1:2*1:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
3653 idr2(2:2*1:2) = idr2(1:2*1:2)
3665 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
3669 if (redblack) i0 = 2 - iand(redblack_cntr, 1)
3672 a0(1:2) = cc(i, i_eps1)
3673 u(1:2) = cc(i-1:i+1:2, n)
3674 a(1:2) = cc(i-1:i+1:2, i_eps1)
3675 c(:) = 2 * a0(:) * a(:) / (a0(:) + a(:)) * idr2
3678 (sum(c(:) * u(:)) - cc(i,
mg_irhs)) / &
3682 end subroutine box_gs_ahelmh
3685 subroutine box_ahelmh(mg, id, nc, i_out)
3686 type(
mg_t),
intent(inout) :: mg
3687 integer,
intent(in) :: id
3688 integer,
intent(in) :: nc
3689 integer,
intent(in) :: i_out
3691 real(
dp) :: idr2(2*1), a0(2*1), u0, u(2*1), a(2*1)
3694 idr2(1:2*1:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
3695 idr2(2:2*1:2) = idr2(1:2*1:2)
3697 associate(cc => mg%boxes(id)%cc, n =>
mg_iphi, &
3701 a0(1:2) = cc(i, i_eps1)
3702 a(1:2) = cc(i-1:i+1:2, i_eps1)
3704 u(1:2) = cc(i-1:i+1:2, n)
3706 cc(i, i_out) = sum(2 * idr2 * &
3707 a0(:)*a(:)/(a0(:) + a(:)) * (u(:) - u0)) - &
3711 end subroutine box_ahelmh
Subroutine that performs Gauss-Seidel relaxation.
Subroutine that performs A * cc(..., i_in) = cc(..., i_out)
Subroutine that performs prolongation to a single child.
To fill ghost cells near physical boundaries.
To fill ghost cells near refinement boundaries.
integer(i8) function, public mg_number_of_unknowns(mg)
Determine total number of unknowns (on leaves)
integer, parameter, public mg_iphi
Index of solution.
integer, parameter, public i8
Type for 64-bit integers.
subroutine, public mg_prolong(mg, lvl, iv, iv_to, method, add)
Prolong variable iv from lvl to variable iv_to at lvl+1.
subroutine, public mg_fill_ghost_cells_lvl(mg, lvl, iv)
Fill ghost cells at a grid level.
integer, parameter, public mg_iveps3
integer, dimension(2, 1), parameter, public mg_child_rev
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_bc_neumann
Value to indicate a Neumann boundary condition.
logical, dimension(1, 2), parameter, public mg_child_low
integer, parameter, public mg_cylindrical
Cylindrical coordinate system.
integer, parameter, public mg_num_children
subroutine, public mg_restrict_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
integer, parameter, public mg_lvl_hi
Maximum allowed grid level.
integer, parameter, public mg_laplacian
Indicates a standard Laplacian.
integer, parameter, public mg_cartesian
Cartesian coordinate system.
subroutine, public mg_apply_op(mg, i_out, op)
Apply operator to the tree and store in variable i_out.
integer, parameter, public mg_physical_boundary
Special value that indicates there is a physical boundary.
elemental logical function, public mg_has_children(box)
Return .true. if a box has children.
subroutine, public mg_timer_end(timer)
integer, parameter, public mg_num_vars
Number of predefined multigrid variables.
integer, parameter, public mg_smoother_gsrb
integer, parameter, public mg_iveps
Index of the variable coefficient (at cell centers)
subroutine, public mg_set_methods(mg)
integer, parameter, public mg_ahelmholtz
Indicates a anisotropic-coefficient Helmholtz equation.
integer, dimension(1, 2), parameter, public mg_neighb_dix
subroutine, public vhelmholtz_set_methods(mg)
real(dp), public, protected vhelmholtz_lambda
Module for implicitly solving diffusion equations.
integer, parameter, public mg_no_box
Special value that indicates there is no box.
integer, parameter, public mg_lvl_lo
Minimum allowed grid level.
subroutine, public mg_restrict_lvl(mg, iv, lvl)
Restrict from lvl to lvl-1.
subroutine, public helmholtz_set_methods(mg)
subroutine, public mg_fill_ghost_cells(mg, iv)
Fill ghost cells at all grid levels.
integer, parameter, public mg_smoother_jacobi
subroutine, public mg_build_rectangle(mg, domain_size, box_size, dx, r_min, periodic, n_finer)
pure integer function, dimension(1), 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,...
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 vhelmholtz_set_lambda(lambda)
subroutine, public mg_ghost_cell_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
subroutine, public mg_timers_show(mg)
integer, parameter, public mg_max_num_vars
Maximum number of variables.
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...
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_restrict(mg, iv)
Restrict all levels.
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 mg_fas_fmg(mg, have_guess, max_res)
Perform FAS-FMG cycle (full approximation scheme, full multigrid).
subroutine, public helmholtz_set_lambda(lambda)
integer, parameter, public mg_iveps1
Indexes of anisotropic variable coefficients.
subroutine, public mg_add_children(mg, id)
subroutine, public laplacian_set_methods(mg)
subroutine, public mg_set_refinement_boundaries(boxes, level)
Create a list of refinement boundaries (from the coarse side)
integer, dimension(2), parameter, public mg_neighb_high_pm
subroutine, public mg_prolong_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
subroutine, public mg_set_next_level_ids(mg, lvl)
subroutine, public mg_set_neighbors_lvl(mg, lvl)
integer, dimension(2), parameter, public mg_neighb_rev
integer, parameter, public dp
Type of reals.
integer, parameter, public mg_bc_dirichlet
Value to indicate a Dirichlet boundary condition.
integer, parameter, public mg_iveps2
integer, parameter, public mg_vlaplacian
Indicates a variable-coefficient Laplacian.
integer, parameter, public mg_helmholtz
Indicates a constant-coefficient Helmholtz equation.
pure integer function, public mg_highest_uniform_lvl(mg)
subroutine, public mg_allocate_storage(mg)
Allocate communication buffers and local boxes for a tree that has already been created.
subroutine, public diffusion_solve_acoeff(mg, dt, order, max_res)
Solve a diffusion equation implicitly, assuming anisotropic diffusion coefficient which has been stor...
integer, parameter, public mg_ndim
Problem dimension.
logical, dimension(2), parameter, public mg_neighb_low
integer, parameter, public mg_smoother_gs
subroutine, public mg_timer_start(timer)
subroutine, public vlaplacian_set_methods(mg)
integer, parameter, public mg_neighb_highx
integer, parameter, public mg_irhs
Index of right-hand side.
subroutine, public mg_deallocate_storage(mg)
Deallocate all allocatable arrays.
subroutine, public mg_get_face_coords(box, nb, nc, x)
Get coordinates at the face of a box.
integer, parameter, public mg_neighb_lowx
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_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.
integer, dimension(1, 2), parameter, public mg_child_adj_nb
subroutine, public ahelmholtz_set_lambda(lambda)
integer, dimension(2), parameter, public mg_neighb_dim
subroutine, public mg_load_balance(mg)
Load balance all boxes in the multigrid tree. Compared to mg_load_balance_simple, this method does a ...
integer, parameter, public mg_spherical
Spherical coordinate system.
integer function, public mg_add_timer(mg, name)
subroutine, public mg_fas_vcycle(mg, highest_lvl, max_res, standalone)
Perform FAS V-cycle (full approximation scheme).
integer, dimension(1, 2), parameter, public mg_child_dix
integer, parameter, public mg_bc_continuous
Value to indicate a continuous boundary condition.
integer, parameter, public mg_max_timers
Maximum number of timers to use.
integer, parameter, public mg_num_neighbors
subroutine, public ahelmholtz_set_methods(mg)
integer, parameter, public mg_ires
Index of residual.
integer, parameter, public mg_iold
Index of previous solution (used for correction)
subroutine, public sort_and_transfer_buffers(mg, dsize)
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_vhelmholtz
Indicates a variable-coefficient Helmholtz equation.
subroutine, public mg_phi_bc_store(mg)
Store boundary conditions for the solution variable, this can speed up calculations if the same bound...
Buffer type (one is used for each pair of communicating processes)
Lists of blocks per refinement level.