MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
m_octree_mg_3d.t
Go to the documentation of this file.
1
2
3! Single module generated from the octree-mg sources.
4! This file can be easier to include in existing projects.
5!
6! Notes:
7! 1. The module name is here extended by _1d, _2d or _3d
8! 2. The free space Poisson solver is not included here.
9! 3. It is best to make changes in the original repository at
10! https://github.com/jannisteunissen/octree-mg
11! 4. This file can be generated as follows:
12! cd octree-mg/single_module
13! ./to_single_module.sh
14!
15! The modules can be compiled with:
16! mpif90 -c m_octree_mg_1d.f90 [other options]
17! mpif90 -c m_octree_mg_2d.f90 [other options]
18! mpif90 -c m_octree_mg_3d.f90 [other options]
19! mpif90 -c m_octree_mg.f90 -cpp -DNDIM=<1,2,3> [other options]
20
21
23 use mpi
24 implicit none
25 private
26
27 !! File ../src/m_data_structures.f90
28
29 !> Type of reals
30 integer, parameter, public :: dp = kind(0.0d0)
31
32 !> Type for 64-bit integers
33 integer, parameter, public :: i8 = selected_int_kind(18)
34
35 !> Indicates a standard Laplacian
36 integer, parameter, public :: mg_laplacian = 1
37
38 !> Indicates a variable-coefficient Laplacian
39 integer, parameter, public :: mg_vlaplacian = 2
40
41 !> Indicates a constant-coefficient Helmholtz equation
42 integer, parameter, public :: mg_helmholtz = 3
43
44 !> Indicates a variable-coefficient Helmholtz equation
45 integer, parameter, public :: mg_vhelmholtz = 4
46
47 !> Indicates a anisotropic-coefficient Helmholtz equation
48 integer, parameter, public :: mg_ahelmholtz = 5
49
50 !> Cartesian coordinate system
51 integer, parameter, public :: mg_cartesian = 1
52 !> Cylindrical coordinate system
53 integer, parameter, public :: mg_cylindrical = 2
54 !> Spherical coordinate system
55 integer, parameter, public :: mg_spherical = 3
56
57 integer, parameter, public :: mg_smoother_gs = 1
58 integer, parameter, public :: mg_smoother_gsrb = 2
59 integer, parameter, public :: mg_smoother_jacobi = 3
60
61 !> Problem dimension
62 integer, parameter, public :: mg_ndim = 3
63
64 !> Number of predefined multigrid variables
65 integer, parameter, public :: mg_num_vars = 4
66 !> Maximum number of variables
67 integer, parameter, public :: mg_max_num_vars = 10
68 !> Index of solution
69 integer, parameter, public :: mg_iphi = 1
70 !> Index of right-hand side
71 integer, parameter, public :: mg_irhs = 2
72 !> Index of previous solution (used for correction)
73 integer, parameter, public :: mg_iold = 3
74 !> Index of residual
75 integer, parameter, public :: mg_ires = 4
76
77 !> Index of the variable coefficient (at cell centers)
78 integer, parameter, public :: mg_iveps = 5
79
80 !> Indexes of anisotropic variable coefficients
81 integer, parameter, public :: mg_iveps1 = 5
82 integer, parameter, public :: mg_iveps2 = 6
83 integer, parameter, public :: mg_iveps3 = 7
84
85 !> Minimum allowed grid level
86 integer, parameter, public :: mg_lvl_lo = -20
87 !> Maximum allowed grid level
88 integer, parameter, public :: mg_lvl_hi = 20
89
90 !> Value to indicate a Dirichlet boundary condition
91 integer, parameter, public :: mg_bc_dirichlet = -10
92
93 !> Value to indicate a Neumann boundary condition
94 integer, parameter, public :: mg_bc_neumann = -11
95
96 !> Value to indicate a continuous boundary condition
97 integer, parameter, public :: mg_bc_continuous = -12
98
99 !> Special value that indicates there is no box
100 integer, parameter, public :: mg_no_box = 0
101 !> Special value that indicates there is a physical boundary
102 integer, parameter, public :: mg_physical_boundary = -1
103
104 !> Maximum number of timers to use
105 integer, parameter, public :: mg_max_timers = 20
106
107 ! Numbering of children (same location as **corners**)
108 integer, parameter, public :: mg_num_children = 8
109
110 ! Index offset for each child
111 integer, parameter, public :: mg_child_dix(3, 8) = reshape( &
112 [0,0,0, 1,0,0, 0,1,0, 1,1,0, &
113 0,0,1, 1,0,1, 0,1,1, 1,1,1], [3,8])
114 ! Reverse child index in each direction
115 integer, parameter, public :: mg_child_rev(8, 3) = reshape( &
116 [2,1,4,3,6,5,8,7, 3,4,1,2,7,8,5,6, 5,6,7,8,1,2,3,4], [8,3])
117 ! Children adjacent to a neighbor
118 integer, parameter, public :: mg_child_adj_nb(4, 6) = reshape( &
119 [1,3,5,7, 2,4,6,8, 1,2,5,6, 3,4,7,8, 1,2,3,4, 5,6,7,8], [4,6])
120 ! Which children have a low index per dimension
121 logical, parameter, public :: mg_child_low(3, 8) = reshape([ &
122 .true., .true., .true., .false., .true., .true., &
123 .true., .false., .true., .false., .false., .true., &
124 .true., .true., .false., .false., .true., .false., &
125 .true., .false., .false., .false., .false., .false.], [3, 8])
126
127 ! Neighbor topology information
128 integer, parameter, public :: mg_num_neighbors = 6
129 integer, parameter, public :: mg_neighb_lowx = 1
130 integer, parameter, public :: mg_neighb_highx = 2
131 integer, parameter, public :: mg_neighb_lowy = 3
132 integer, parameter, public :: mg_neighb_highy = 4
133 integer, parameter, public :: mg_neighb_lowz = 5
134 integer, parameter, public :: mg_neighb_highz = 6
135 ! Index offsets of neighbors
136 integer, parameter, public :: mg_neighb_dix(3, 6) = reshape( &
137 [-1,0,0, 1,0,0, 0,-1,0, 0,1,0, 0,0,-1, 0,0,1], [3,6])
138 ! Which neighbors have a lower index
139 logical, parameter, public :: mg_neighb_low(6) = &
140 [.true., .false., .true., .false., .true., .false.]
141 ! Opposite of nb_low, but now as -1,1 integers
142 integer, parameter, public :: mg_neighb_high_pm(6) = [-1, 1, -1, 1, -1, 1]
143 ! Reverse neighbors
144 integer, parameter, public :: mg_neighb_rev(6) = [2, 1, 4, 3, 6, 5]
145 ! Direction (dimension) for a neighbor
146 integer, parameter, public :: mg_neighb_dim(6) = [1, 1, 2, 2, 3, 3]
147
148 !> Lists of blocks per refinement level
149 type, public :: mg_lvl_t
150 integer, allocatable :: leaves(:)
151 integer, allocatable :: parents(:)
152 integer, allocatable :: ref_bnds(:)
153 integer, allocatable :: ids(:)
154 integer, allocatable :: my_leaves(:)
155 integer, allocatable :: my_parents(:)
156 integer, allocatable :: my_ref_bnds(:)
157 integer, allocatable :: my_ids(:)
158 end type mg_lvl_t
159
160 !> Box data structure
161 type, public :: mg_box_t
162 integer :: rank !< Which process owns this box
163 integer :: id !< Box id (index in boxes(:) array)
164 integer :: lvl !< Refinement level
165 integer :: ix(3) !< Spatial index
166 integer :: parent !< Id of parent
167 integer :: children(2**3) !< Ids of children
168 integer :: neighbors(2*3) !< Ids of neighbors
169 real(dp) :: r_min(3) !< Minimum coordinate
170 real(dp) :: dr(3) !< Grid spacing
171 !> Cell-centered data
172 real(dp), allocatable :: cc(:, :, :, :)
173 end type mg_box_t
174
175 !> Buffer type (one is used for each pair of communicating processes)
176 type, public :: mg_buf_t
177 integer :: i_send !< Index in send array
178 integer :: i_recv
179 integer :: i_ix
180 integer, allocatable :: ix(:) !< Will be used to sort the data
181 real(dp), allocatable :: send(:)
182 real(dp), allocatable :: recv(:)
183 end type mg_buf_t
184
185 type, public :: mg_comm_t
186 integer, allocatable :: n_send(:, :)
187 integer, allocatable :: n_recv(:, :)
188 end type mg_comm_t
189
190 type, public :: mg_bc_t
191 integer :: bc_type = mg_bc_dirichlet !< Type of boundary condition
192 real(dp) :: bc_value = 0.0_dp !< Value (for e.g. Dirichlet or Neumann)
193 !> To set user-defined boundary conditions (overrides bc(:))
194 procedure(mg_subr_bc), pointer, nopass :: boundary_cond => null()
195 !> To set a user-defined refinement boundary method
196 procedure(mg_subr_rb), pointer, nopass :: refinement_bnd => null()
197 end type mg_bc_t
198
199 type, public :: mg_timer_t
200 character(len=20) :: name
201 real(dp) :: t = 0.0_dp
202 real(dp) :: t0
203 end type mg_timer_t
204
205 type, public :: mg_t
206 !> Whether the multigrid tree structure has been created
207 logical :: tree_created = .false.
208 !> Whether storage has been allocated for boxes
209 logical :: is_allocated = .false.
210 !> Number of extra cell-centered variable (e.g., for coefficients)
211 integer :: n_extra_vars = 0
212 !> MPI communicator
213 integer :: comm = -1
214 !> Number of MPI tasks
215 integer :: n_cpu = -1
216 !> MPI rank of this task
217 integer :: my_rank = -1
218 !> Size of boxes in cells, be equal for all dimensions
219 integer :: box_size = -1
220 !> Highest grid level in the tree
221 integer :: highest_lvl = -1
222 !> Lowest grid level in the tree
223 integer :: lowest_lvl = -1
224 !> First normal level of the quadtree/octree, at coarser levels parents
225 !> have only one child
226 integer :: first_normal_lvl = -1
227 !> Total number of boxes in the tree (among all processes)
228 integer :: n_boxes = 0
229 !> Size of boxes per level (differs for coarsest levels)
230 integer :: box_size_lvl(mg_lvl_lo:mg_lvl_hi)
231 !> Size of domain per level (if uniformly refined)
232 integer :: domain_size_lvl(3, mg_lvl_lo:mg_lvl_hi)
233 !> Grid spacing per level
234 real(dp) :: dr(3, mg_lvl_lo:mg_lvl_hi)
235 !> Minimum coordinates
236 real(dp) :: r_min(3)
237 !> List of all levels
239 !> Array with all boxes in the tree. Only boxes owned by this task are
240 !> actually allocated
241 type(mg_box_t), allocatable :: boxes(:)
242 !> Buffer for communication with each other task
243 type(mg_buf_t), allocatable :: buf(:)
244
245 !> Communication info for restriction
246 type(mg_comm_t) :: comm_restrict
247 !> Communication info for prolongation
248 type(mg_comm_t) :: comm_prolong
249 !> Communication info for ghost cell filling
250 type(mg_comm_t) :: comm_ghostcell
251
252 !> Whether boundary condition data has been stored for mg solution
253 logical :: phi_bc_data_stored = .false.
254
255 !> Whether a dimension is periodic
256 logical :: periodic(3) = .false.
257
258 !> To store pre-defined boundary conditions per direction per variable
260
261 !> Type of operator
262 integer :: operator_type = mg_laplacian
263 !> Type of grid geometry
264 integer :: geometry_type = mg_cartesian
265
266 !> Whether the mean has to be subtracted from the multigrid solution
267 logical :: subtract_mean = .false.
268 !> Type of multigrid smoother
269 integer :: smoother_type = mg_smoother_gs
270 !> Number of substeps for the smoother (for GSRB this is 2)
271 integer :: n_smoother_substeps = 1
272 !> Number of cycles when doing downwards relaxation
273 integer :: n_cycle_down = 2
274 !> Number of cycles when doing upwards relaxation
275 integer :: n_cycle_up = 2
276 !> Maximum number of cycles on the coarse grid
277 integer :: max_coarse_cycles = 1000
278 integer :: coarsest_grid(3) = 2
279 !> Stop coarse grid when max. residual is smaller than this
280 real(dp) :: residual_coarse_abs = 1e-8_dp
281 !> Stop coarse grid when residual has been reduced by this factor
282 real(dp) :: residual_coarse_rel = 1e-8_dp
283
284 !> Multigrid operator (e.g., Laplacian)
285 procedure(mg_box_op), pointer, nopass :: box_op => null()
286
287 !> Multigrid smoother
288 procedure(mg_box_gsrb), pointer, nopass :: box_smoother => null()
289
290 !> Multigrid prolongation method
291 procedure(mg_box_prolong), pointer, nopass :: box_prolong => null()
292
293 !> Number of timers
294 integer :: n_timers = 0
295 !> Values for the timers
296 type(mg_timer_t) :: timers(mg_max_timers)
297 end type mg_t
298
299 interface
300 !> To fill ghost cells near physical boundaries
301 subroutine mg_subr_bc(box, nc, iv, nb, bc_type, bc)
302 import
303 type(mg_box_t), intent(in) :: box
304 integer, intent(in) :: nc
305 integer, intent(in) :: iv !< Index of variable
306 integer, intent(in) :: nb !< Direction
307 integer, intent(out) :: bc_type !< Type of b.c.
308 !> Boundary values
309 real(dp), intent(out) :: bc(nc, nc)
310 end subroutine mg_subr_bc
311
312 !> To fill ghost cells near refinement boundaries
313 subroutine mg_subr_rb(box, nc, iv, nb, cgc)
314 import
315 type(mg_box_t), intent(inout) :: box
316 integer, intent(in) :: nc
317 integer, intent(in) :: iv !< Index of variable
318 integer, intent(in) :: nb !< Direction
319 !> Coarse data
320 real(dp), intent(in) :: cgc(nc, nc)
321 end subroutine mg_subr_rb
322
323 !> Subroutine that performs A * cc(..., i_in) = cc(..., i_out)
324 subroutine mg_box_op(mg, id, nc, i_out)
325 import
326 type(mg_t), intent(inout) :: mg
327 integer, intent(in) :: id
328 integer, intent(in) :: nc
329 integer, intent(in) :: i_out
330 end subroutine mg_box_op
331
332 !> Subroutine that performs Gauss-Seidel relaxation
333 subroutine mg_box_gsrb(mg, id, nc, redblack_cntr)
334 import
335 type(mg_t), intent(inout) :: mg
336 integer, intent(in) :: id
337 integer, intent(in) :: nc
338 integer, intent(in) :: redblack_cntr
339 end subroutine mg_box_gsrb
340
341 !> Subroutine that performs prolongation to a single child
342 subroutine mg_box_prolong(mg, p_id, dix, nc, iv, fine)
343 import
344 type(mg_t), intent(inout) :: mg
345 integer, intent(in) :: p_id !< Id of parent
346 integer, intent(in) :: dix(3) !< Offset of child in parent grid
347 integer, intent(in) :: nc !< Child grid size
348 integer, intent(in) :: iv !< Prolong from this variable
349 real(dp), intent(out) :: fine(nc, nc, nc) !< Prolonged values
350 end subroutine mg_box_prolong
351 end interface
352
353 ! Public methods
354 public :: mg_subr_bc
355 public :: mg_subr_rb
356 public :: mg_box_op
357 public :: mg_box_gsrb
358 public :: mg_box_prolong
359 public :: mg_has_children
360 public :: mg_ix_to_ichild
361 public :: mg_get_child_offset
362 public :: mg_highest_uniform_lvl
363 public :: mg_number_of_unknowns
364 public :: mg_get_face_coords
365 public :: mg_add_timer
366 public :: mg_timer_start
367 public :: mg_timer_end
368 public :: mg_timers_show
369
370 !! File ../src/m_allocate_storage.f90
371
372 public :: mg_allocate_storage
373 public :: mg_deallocate_storage
374
375 !! File ../src/m_diffusion.f90
376
377!> Module for implicitly solving diffusion equations
378
379 public :: diffusion_solve
380 public :: diffusion_solve_vcoeff
381 public :: diffusion_solve_acoeff
382
383 !! File ../src/m_laplacian.f90
384
385!> Module which contains multigrid procedures for a Laplacian operator
386
387 public :: laplacian_set_methods
388
389 !! File ../src/m_vhelmholtz.f90
390
391!> Module which contains multigrid procedures for a Helmholtz operator of the
392!> form: div(D grad(phi)) - lambda*phi = f, where D has a smooth spatial
393!> variation
394
395 !> The lambda used for the Helmholtz equation (should be >= 0)
396 real(dp), public, protected :: vhelmholtz_lambda = 0.0_dp
397
398 public :: vhelmholtz_set_methods
399 public :: vhelmholtz_set_lambda
400
401 !! File ../src/m_build_tree.f90
402
403 ! Public methods
404 public :: mg_build_rectangle
405 public :: mg_add_children
406 public :: mg_set_leaves_parents
407 public :: mg_set_next_level_ids
409 public :: mg_set_neighbors_lvl
410
411 !! File ../src/m_load_balance.f90
412!> Module for load balancing a tree (that already has been constructed). The
413!> load balancing determines which ranks (MPI processes) allocated physical
414!> storage for boxes. The tree structure itself is present on all processes.
415
416 ! Public methods
417 public :: mg_load_balance_simple
418 public :: mg_load_balance
420
421 !! File ../src/m_vlaplacian.f90
422
423!> Module which contains multigrid procedures for a variable-coefficient
424!> Laplacian operator, assuming the variation is smooth
425
426 public :: vlaplacian_set_methods
427
428 !! File ../src/m_communication.f90
429
430 ! Public methods
431 public :: mg_comm_init
433
434 !! File ../src/m_ghost_cells.f90
435
436 ! Public methods
438 public :: mg_fill_ghost_cells
440 public :: mg_phi_bc_store
441
442 !! File ../src/m_prolong.f90
443
444 ! Public methods
445 public :: mg_prolong
446 public :: mg_prolong_buffer_size
447 public :: mg_prolong_sparse
448
449 !! File ../src/m_helmholtz.f90
450
451!> Module which contains multigrid procedures for a Helmholtz operator of the
452!> form: laplacian(phi) - lambda*phi = f
453
454 !> The lambda used for the Helmholtz equation (should be >= 0)
455 real(dp), public, protected :: helmholtz_lambda = 0.0_dp
456
457 public :: helmholtz_set_methods
458 public :: helmholtz_set_lambda
459
460 !! File ../src/m_multigrid.f90
461
462 integer :: timer_total_vcycle = -1
463 integer :: timer_total_fmg = -1
464 integer :: timer_smoother = -1
465 integer :: timer_smoother_gc = -1
466 integer :: timer_coarse = -1
467 integer :: timer_correct = -1
468 integer :: timer_update_coarse = -1
469
470 ! Public methods
471 public :: mg_fas_vcycle
472 public :: mg_fas_fmg
473 public :: mg_set_methods
474 public :: mg_apply_op
475
476 !! File ../src/m_restrict.f90
477
478 ! Public methods
479 public :: mg_restrict
480 public :: mg_restrict_lvl
482
483 !! File ../src/m_mrgrnk.f90
484
485 Integer, Parameter :: kdp = selected_real_kind(15)
486 public :: mrgrnk
487 private :: kdp
488 private :: i_mrgrnk
489 interface mrgrnk
490 module procedure i_mrgrnk
491 end interface mrgrnk
492
493 !! File ../src/m_ahelmholtz.f90
494
495!> Module which contains multigrid procedures for a Helmholtz operator of the
496!> form: div(D grad(phi)) - lambda*phi = f, where D has a smooth spatial
497!> variation and a component in each spatial direction
498
499 !> The lambda used for the Helmholtz equation (should be >= 0)
500 real(dp), public, protected :: ahelmholtz_lambda = 0.0_dp
501
502 public :: ahelmholtz_set_methods
503 public :: ahelmholtz_set_lambda
504contains
505
506 !! File ../src/m_data_structures.f90
507
508 !> Return .true. if a box has children
509 elemental logical function mg_has_children(box)
510 type(mg_box_t), intent(in) :: box
511
512 ! Boxes are either fully refined or not, so we only need to check one of the
513 ! children
514 mg_has_children = (box%children(1) /= mg_no_box)
515 end function mg_has_children
516
517 !> Compute the child index for a box with spatial index ix. With child index
518 !> we mean the index in the children(:) array of its parent.
519 integer function mg_ix_to_ichild(ix)
520 integer, intent(in) :: ix(3) !< Spatial index of the box
521 ! The index can range from 1 (all ix odd) and 2**$D (all ix even)
522 mg_ix_to_ichild = 8 - 4 * iand(ix(3), 1) - &
523 2 * iand(ix(2), 1) - iand(ix(1), 1)
524 end function mg_ix_to_ichild
525
526 !> Get the offset of a box with respect to its parent (e.g. in 2d, there can
527 !> be a child at offset 0,0, one at n_cell/2,0, one at 0,n_cell/2 and one at
528 !> n_cell/2, n_cell/2)
529 pure function mg_get_child_offset(mg, id) result(ix_offset)
530 type(mg_t), intent(in) :: mg
531 integer, intent(in) :: id
532 integer :: ix_offset(3)
533
534 if (mg%boxes(id)%lvl <= mg%first_normal_lvl) then
535 ix_offset(:) = 0
536 else
537 ix_offset = iand(mg%boxes(id)%ix-1, 1) * &
538 ishft(mg%box_size, -1) ! * n_cell / 2
539 end if
540 end function mg_get_child_offset
541
542 pure function mg_highest_uniform_lvl(mg) result(lvl)
543 type(mg_t), intent(in) :: mg
544 integer :: lvl
545
546 do lvl = mg%first_normal_lvl, mg%highest_lvl-1
547 ! Exit if a grid is partially refined
548 if (size(mg%lvls(lvl)%leaves) /= 0 .and. &
549 size(mg%lvls(lvl)%parents) /= 0) exit
550 end do
551 ! If the loop did not exit, we get lvl equals mg%highest_lvl
552 end function mg_highest_uniform_lvl
553
554 !> Determine total number of unknowns (on leaves)
555 function mg_number_of_unknowns(mg) result(n_unknowns)
556 type(mg_t), intent(in) :: mg
557 integer :: lvl
558 integer(i8) :: n_unknowns
559
560 n_unknowns = 0
561 do lvl = mg%first_normal_lvl, mg%highest_lvl
562 n_unknowns = n_unknowns + size(mg%lvls(lvl)%leaves)
563 end do
564 n_unknowns = n_unknowns * int(mg%box_size**3, i8)
565 end function mg_number_of_unknowns
566
567 !> Get coordinates at the face of a box
568 subroutine mg_get_face_coords(box, nb, nc, x)
569 type(mg_box_t), intent(in) :: box
570 integer, intent(in) :: nb
571 integer, intent(in) :: nc
572 real(dp), intent(out) :: x(nc, nc, 3)
573 integer :: i, j, ixs(3-1)
574 integer :: nb_dim
575 real(dp) :: rmin(3)
576
577 nb_dim = mg_neighb_dim(nb)
578
579 ! Determine directions perpendicular to neighbor
580 ixs = [(i, i = 1, 3-1)]
581 ixs(nb_dim:) = ixs(nb_dim:) + 1
582
583 rmin = box%r_min
584 if (.not. mg_neighb_low(nb)) then
585 rmin(nb_dim) = rmin(nb_dim) + box%dr(nb_dim) * nc
586 end if
587
588 do j = 1, nc
589 do i = 1, nc
590 x(i, j, :) = rmin
591 x(i, j, ixs) = x(i, j, ixs) + ([i, j] - 0.5d0) * box%dr(ixs)
592 end do
593 end do
594 end subroutine mg_get_face_coords
595
596 integer function mg_add_timer(mg, name)
597 type(mg_t), intent(inout) :: mg
598 character(len=*), intent(in) :: name
599
600 mg%n_timers = mg%n_timers + 1
601 mg_add_timer = mg%n_timers
602 mg%timers(mg_add_timer)%name = name
603 end function mg_add_timer
604
605 subroutine mg_timer_start(timer)
606 use mpi
607 type(mg_timer_t), intent(inout) :: timer
608 timer%t0 = mpi_wtime()
609 end subroutine mg_timer_start
610
611 subroutine mg_timer_end(timer)
612 use mpi
613 type(mg_timer_t), intent(inout) :: timer
614 timer%t = timer%t + mpi_wtime() - timer%t0
615 end subroutine mg_timer_end
616
617 subroutine mg_timers_show(mg)
618 use mpi
619 type(mg_t), intent(in) :: mg
620 integer :: n, ierr
621 real(dp) :: tmin(mg%n_timers)
622 real(dp) :: tmax(mg%n_timers)
623
624 call mpi_reduce(mg%timers(1:mg%n_timers)%t, tmin, mg%n_timers, &
625 mpi_double, mpi_min, 0, mg%comm, ierr)
626 call mpi_reduce(mg%timers(1:mg%n_timers)%t, tmax, mg%n_timers, &
627 mpi_double, mpi_max, 0, mg%comm, ierr)
628
629 if (mg%my_rank == 0) then
630 write(*, "(A20,2A16)") "name ", "min(s)", "max(s)"
631 do n = 1, mg%n_timers
632 write(*, "(A20,2F16.6)") mg%timers(n)%name, &
633 tmin(n), tmax(n)
634 end do
635 end if
636 end subroutine mg_timers_show
637
638 !! File ../src/m_allocate_storage.f90
639
640 !> Deallocate all allocatable arrays
642 type(mg_t), intent(inout) :: mg
643 integer :: lvl
644
645 if (.not. mg%is_allocated) &
646 error stop "deallocate_storage: tree is not allocated"
647
648 deallocate(mg%boxes)
649 deallocate(mg%buf)
650
651 deallocate(mg%comm_restrict%n_send)
652 deallocate(mg%comm_restrict%n_recv)
653
654 deallocate(mg%comm_prolong%n_send)
655 deallocate(mg%comm_prolong%n_recv)
656
657 deallocate(mg%comm_ghostcell%n_send)
658 deallocate(mg%comm_ghostcell%n_recv)
659
660 do lvl = mg%lowest_lvl, mg%highest_lvl
661 deallocate(mg%lvls(lvl)%ids)
662 deallocate(mg%lvls(lvl)%leaves)
663 deallocate(mg%lvls(lvl)%parents)
664 deallocate(mg%lvls(lvl)%ref_bnds)
665 deallocate(mg%lvls(lvl)%my_ids)
666 deallocate(mg%lvls(lvl)%my_leaves)
667 deallocate(mg%lvls(lvl)%my_parents)
668 deallocate(mg%lvls(lvl)%my_ref_bnds)
669 end do
670
671 mg%is_allocated = .false.
672 mg%n_boxes = 0
673 mg%phi_bc_data_stored = .false.
674 end subroutine mg_deallocate_storage
675
676 !> Allocate communication buffers and local boxes for a tree that has already
677 !> been created
678 subroutine mg_allocate_storage(mg)
679
680 type(mg_t), intent(inout) :: mg
681 integer :: i, id, lvl, nc
682 integer :: n_send(0:mg%n_cpu-1, 3)
683 integer :: n_recv(0:mg%n_cpu-1, 3)
684 integer :: dsize(3)
685 integer :: n_in, n_out, n_id
686
687 if (.not. mg%tree_created) &
688 error stop "allocate_storage: tree is not yet created"
689
690 if (mg%is_allocated) &
691 error stop "allocate_storage: tree is already allocated"
692
693 do lvl = mg%lowest_lvl, mg%highest_lvl
694 nc = mg%box_size_lvl(lvl)
695 do i = 1, size(mg%lvls(lvl)%my_ids)
696 id = mg%lvls(lvl)%my_ids(i)
697 allocate(mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, &
698 mg_num_vars + mg%n_extra_vars))
699
700 ! Set all initial values to zero
701 mg%boxes(id)%cc(:, :, :, :) = 0.0_dp
702 end do
703 end do
704
705 allocate(mg%buf(0:mg%n_cpu-1))
706
707 call mg_ghost_cell_buffer_size(mg, n_send(:, 1), &
708 n_recv(:, 1), dsize(1))
709 call mg_restrict_buffer_size(mg, n_send(:, 2), &
710 n_recv(:, 2), dsize(2))
711 call mg_prolong_buffer_size(mg, n_send(:, 3), &
712 n_recv(:, 3), dsize(3))
713
714 do i = 0, mg%n_cpu-1
715 n_out = maxval(n_send(i, :) * dsize(:))
716 n_in = maxval(n_recv(i, :) * dsize(:))
717 n_id = maxval(n_send(i, :))
718 allocate(mg%buf(i)%send(n_out))
719 allocate(mg%buf(i)%recv(n_in))
720 allocate(mg%buf(i)%ix(n_id))
721 end do
722
723 mg%is_allocated = .true.
724 end subroutine mg_allocate_storage
725
726 !! File ../src/m_diffusion.f90
727
728 !> Solve a diffusion equation implicitly, assuming a constant diffusion
729 !> coefficient. The solution at time t should be stored in mg_iphi, which is
730 !> on output replaced by the solution at time t+dt.
731 subroutine diffusion_solve(mg, dt, diffusion_coeff, order, max_res)
732
733 type(mg_t), intent(inout) :: mg
734 real(dp), intent(in) :: dt
735 real(dp), intent(in) :: diffusion_coeff
736 integer, intent(in) :: order
737 real(dp), intent(in) :: max_res
738 integer, parameter :: max_its = 10
739 integer :: n
740 real(dp) :: res
741
742 mg%operator_type = mg_helmholtz
743 call mg_set_methods(mg)
744
745 select case (order)
746 case (1)
747 call helmholtz_set_lambda(1/(dt * diffusion_coeff))
748 call set_rhs(mg, -1/(dt * diffusion_coeff), 0.0_dp)
749 case (2)
750 call helmholtz_set_lambda(0.0d0)
751 call mg_apply_op(mg, mg_irhs)
752 call helmholtz_set_lambda(2/(dt * diffusion_coeff))
753 call set_rhs(mg, -2/(dt * diffusion_coeff), -1.0_dp)
754 case default
755 error stop "diffusion_solve: order should be 1 or 2"
756 end select
757
758 ! Start with an FMG cycle
759 call mg_fas_fmg(mg, .true., max_res=res)
760
761 ! Add V-cycles if necessary
762 do n = 1, max_its
763 if (res <= max_res) exit
764 call mg_fas_vcycle(mg, max_res=res)
765 end do
766
767 if (n == max_its + 1) then
768 if (mg%my_rank == 0) &
769 print *, "Did you specify boundary conditions correctly?"
770 error stop "diffusion_solve: no convergence"
771 end if
772 end subroutine diffusion_solve
773
774 !> Solve a diffusion equation implicitly, assuming a variable diffusion
775 !> coefficient which has been stored in mg_iveps (also on coarse grids). The
776 !> solution at time t should be stored in mg_iphi, which is on output replaced
777 !> by the solution at time t+dt.
778 subroutine diffusion_solve_vcoeff(mg, dt, order, max_res)
779
780 type(mg_t), intent(inout) :: mg
781 real(dp), intent(in) :: dt
782 integer, intent(in) :: order
783 real(dp), intent(in) :: max_res
784 integer, parameter :: max_its = 10
785 integer :: n
786 real(dp) :: res
787
788 mg%operator_type = mg_vhelmholtz
789 call mg_set_methods(mg)
790
791 select case (order)
792 case (1)
793 call vhelmholtz_set_lambda(1/dt)
794 call set_rhs(mg, -1/dt, 0.0_dp)
795 case (2)
796 call vhelmholtz_set_lambda(0.0d0)
797 call mg_apply_op(mg, mg_irhs)
798 call vhelmholtz_set_lambda(2/dt)
799 call set_rhs(mg, -2/dt, -1.0_dp)
800 case default
801 error stop "diffusion_solve: order should be 1 or 2"
802 end select
803
804 ! Start with an FMG cycle
805 call mg_fas_fmg(mg, .true., max_res=res)
806
807 ! Add V-cycles if necessary
808 do n = 1, max_its
809 if (res <= max_res) exit
810 call mg_fas_vcycle(mg, max_res=res)
811 end do
812
813 if (n == max_its + 1) then
814 if (mg%my_rank == 0) then
815 print *, "Did you specify boundary conditions correctly?"
816 print *, "Or is the variation in diffusion too large?"
817 end if
818 error stop "diffusion_solve: no convergence"
819 end if
820 end subroutine diffusion_solve_vcoeff
821
822 !> Solve a diffusion equation implicitly, assuming anisotropic diffusion
823 !> coefficient which has been stored in mg_iveps1, mg_iveps2, mg_iveps3
824 !> (also on coarse grids). The
825 !> solution at time t should be stored in mg_iphi, which is on output replaced
826 !> by the solution at time t+dt.
827 subroutine diffusion_solve_acoeff(mg, dt, order, max_res)
828
829 type(mg_t), intent(inout) :: mg
830 real(dp), intent(in) :: dt
831 integer, intent(in) :: order
832 real(dp), intent(in) :: max_res
833 integer, parameter :: max_its = 10
834 integer :: n
835 real(dp) :: res
836
837 mg%operator_type = mg_ahelmholtz
838 call mg_set_methods(mg)
839
840 select case (order)
841 case (1)
842 call ahelmholtz_set_lambda(1/dt)
843 call set_rhs(mg, -1/dt, 0.0_dp)
844 case (2)
845 call ahelmholtz_set_lambda(0.0d0)
846 call mg_apply_op(mg, mg_irhs)
847 call ahelmholtz_set_lambda(2/dt)
848 call set_rhs(mg, -2/dt, -1.0_dp)
849 case default
850 error stop "diffusion_solve: order should be 1 or 2"
851 end select
852
853 ! Start with an FMG cycle
854 call mg_fas_fmg(mg, .true., max_res=res)
855
856 ! Add V-cycles if necessary
857 do n = 1, max_its
858 if (res <= max_res) exit
859 call mg_fas_vcycle(mg, max_res=res)
860 end do
861
862 if (n == max_its + 1) then
863 if (mg%my_rank == 0) then
864 print *, "Did you specify boundary conditions correctly?"
865 print *, "Or is the variation in diffusion too large?"
866 end if
867 error stop "diffusion_solve: no convergence"
868 end if
869 end subroutine diffusion_solve_acoeff
870
871 subroutine set_rhs(mg, f1, f2)
872 type(mg_t), intent(inout) :: mg
873 real(dp), intent(in) :: f1, f2
874 integer :: lvl, i, id, nc
875
876 do lvl = 1, mg%highest_lvl
877 nc = mg%box_size_lvl(lvl)
878 do i = 1, size(mg%lvls(lvl)%my_leaves)
879 id = mg%lvls(lvl)%my_leaves(i)
880 mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, mg_irhs) = &
881 f1 * mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, mg_iphi) + &
882 f2 * mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, mg_irhs)
883 end do
884 end do
885 end subroutine set_rhs
886
887 !! File ../src/m_laplacian.f90
888
890 type(mg_t), intent(inout) :: mg
891
892 if (all(mg%periodic)) then
893 ! For a fully periodic Laplacian, remove the mean from the rhs and phi so
894 ! that a unique and periodic solution can be found
895 mg%subtract_mean = .true.
896 end if
897
898 select case (mg%geometry_type)
899 case (mg_cartesian)
900 mg%box_op => box_lpl
901
902 select case (mg%smoother_type)
903 ! case (smoother_jacobi)
904 ! mg%box_smoother => box_jacobi_lpl
906 mg%box_smoother => box_gs_lpl
907 case default
908 error stop "laplacian_set_methods: unsupported smoother type"
909 end select
910 case default
911 error stop "laplacian_set_methods: unsupported geometry"
912 end select
913
914 end subroutine laplacian_set_methods
915
916 !> Perform Gauss-Seidel relaxation on box for a Laplacian operator
917 subroutine box_gs_lpl(mg, id, nc, redblack_cntr)
918 type(mg_t), intent(inout) :: mg
919 integer, intent(in) :: id
920 integer, intent(in) :: nc
921 integer, intent(in) :: redblack_cntr !< Iteration counter
922 integer :: i, j, k, i0, di
923 real(dp) :: idr2(3), fac
924 logical :: redblack
925 real(dp), parameter :: sixth = 1/6.0_dp
926
927 idr2 = 1/mg%dr(:, mg%boxes(id)%lvl)**2
928 fac = 0.5_dp / sum(idr2)
929 i0 = 1
930
931 redblack = (mg%smoother_type == mg_smoother_gsrb)
932 if (redblack) then
933 di = 2
934 else
935 di = 1
936 end if
937
938 ! The parity of redblack_cntr determines which cells we use. If
939 ! redblack_cntr is even, we use the even cells and vice versa.
940 associate(cc => mg%boxes(id)%cc, n => mg_iphi)
941 do k = 1, nc
942 do j = 1, nc
943 if (redblack) &
944 i0 = 2 - iand(ieor(redblack_cntr, k+j), 1)
945 do i = i0, nc, di
946 cc(i, j, k, n) = fac * ( &
947 idr2(1) * (cc(i+1, j, k, n) + cc(i-1, j, k, n)) + &
948 idr2(2) * (cc(i, j+1, k, n) + cc(i, j-1, k, n)) + &
949 idr2(3) * (cc(i, j, k+1, n) + cc(i, j, k-1, n)) - &
950 cc(i, j, k, mg_irhs))
951 end do
952 end do
953 end do
954 end associate
955 end subroutine box_gs_lpl
956
957! !> Perform Jacobi relaxation on box for a Laplacian operator
958! subroutine box_jacobi_lpl(mg, id, nc, cntr)
959! type(mg_t), intent(inout) :: mg
960! integer, intent(in) :: id
961! integer, intent(in) :: nc
962! integer, intent(in) :: cntr !< Not used
963! integer :: i, j, k
964! real(dp), parameter :: w = 2.0_dp / 3
965! real(dp) :: tmp(0:nc+1, 0:nc+1, 0:nc+1)
966! real(dp) :: dr2
967! #if 3 == 3
968! real(dp), parameter :: sixth = 1/6.0_dp
969! #endif
970
971! dr2 = mg%dr(mg%boxes(id)%lvl)**2
972
973! associate (box => mg%boxes(id))
974! tmp = box%cc(:, :, :, mg_iphi)
975! do j=1, nc; do i=1, nc; do k=1, nc
976! #if 3 == 2
977! box%cc(i, j, mg_iphi) = (1-w) * box%cc(i, j, mg_iphi) + &
978! 0.25_dp * w * ( &
979! tmp(i+1, j) + tmp(i-1, j) + &
980! tmp(i, j+1) + tmp(i, j-1) - &
981! dr2 * box%cc(i, j, mg_irhs))
982! #elif 3 == 3
983! box%cc(i, j, k, mg_iphi) = (1-w) * &
984! box%cc(i, j, k, mg_iphi) + &
985! sixth * w * ( &
986! tmp(i+1, j, k) + tmp(i-1, j, k) + &
987! tmp(i, j+1, k) + tmp(i, j-1, k) + &
988! tmp(i, j, k+1) + tmp(i, j, k-1) - &
989! dr2 * box%cc(i, j, k, mg_irhs))
990! #endif
991! end do; end do; end do
992! end associate
993! end subroutine box_jacobi_lpl
994
995 !> Perform Laplacian operator on a box
996 subroutine box_lpl(mg, id, nc, i_out)
997 type(mg_t), intent(inout) :: mg
998 integer, intent(in) :: id
999 integer, intent(in) :: nc
1000 integer, intent(in) :: i_out !< Index of variable to store Laplacian in
1001 integer :: i, j, k
1002 real(dp) :: idr2(3)
1003
1004 idr2 = 1 / mg%dr(:, mg%boxes(id)%lvl)**2
1005
1006 associate(cc => mg%boxes(id)%cc, n => mg_iphi)
1007 do k = 1, nc
1008 do j = 1, nc
1009 do i = 1, nc
1010 cc(i, j, k, i_out) = &
1011 idr2(1) * (cc(i-1, j, k, n) + cc(i+1, j, k, n) &
1012 - 2 * cc(i, j, k, n)) &
1013 + idr2(2) * (cc(i, j-1, k, n) + cc(i, j+1, k, n) &
1014 - 2 * cc(i, j, k, n)) &
1015 + idr2(3) * (cc(i, j, k-1, n) + cc(i, j, k+1, n) &
1016 - 2 * cc(i, j, k, n))
1017 end do
1018 end do
1019 end do
1020 end associate
1021 end subroutine box_lpl
1022
1023
1024 !! File ../src/m_vhelmholtz.f90
1025
1027 type(mg_t), intent(inout) :: mg
1028
1029 if (mg%n_extra_vars == 0 .and. mg%is_allocated) then
1030 error stop "vhelmholtz_set_methods: mg%n_extra_vars == 0"
1031 else
1032 mg%n_extra_vars = max(1, mg%n_extra_vars)
1033 end if
1034
1035 mg%subtract_mean = .false.
1036
1037 ! Use Neumann zero boundary conditions for the variable coefficient, since
1038 ! it is needed in ghost cells.
1039 mg%bc(:, mg_iveps)%bc_type = mg_bc_neumann
1040 mg%bc(:, mg_iveps)%bc_value = 0.0_dp
1041
1042 select case (mg%geometry_type)
1043 case (mg_cartesian)
1044 mg%box_op => box_vhelmh
1045
1046 select case (mg%smoother_type)
1048 mg%box_smoother => box_gs_vhelmh
1049 case default
1050 error stop "vhelmholtz_set_methods: unsupported smoother type"
1051 end select
1052 case default
1053 error stop "vhelmholtz_set_methods: unsupported geometry"
1054 end select
1055
1056 end subroutine vhelmholtz_set_methods
1057
1058 subroutine vhelmholtz_set_lambda(lambda)
1059 real(dp), intent(in) :: lambda
1060
1061 if (lambda < 0) &
1062 error stop "vhelmholtz_set_lambda: lambda < 0 not allowed"
1063
1064 vhelmholtz_lambda = lambda
1065 end subroutine vhelmholtz_set_lambda
1066
1067 !> Perform Gauss-Seidel relaxation on box for a Helmholtz operator
1068 subroutine box_gs_vhelmh(mg, id, nc, redblack_cntr)
1069 type(mg_t), intent(inout) :: mg
1070 integer, intent(in) :: id
1071 integer, intent(in) :: nc
1072 integer, intent(in) :: redblack_cntr !< Iteration counter
1073 integer :: i, j, k, i0, di
1074 logical :: redblack
1075 real(dp) :: idr2(2*3), u(2*3)
1076 real(dp) :: a0, a(2*3), c(2*3)
1077
1078 ! Duplicate 1/dr^2 array to multiply neighbor values
1079 idr2(1:2*3:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1080 idr2(2:2*3:2) = idr2(1:2*3:2)
1081 i0 = 1
1082
1083 redblack = (mg%smoother_type == mg_smoother_gsrb)
1084 if (redblack) then
1085 di = 2
1086 else
1087 di = 1
1088 end if
1089
1090 ! The parity of redblack_cntr determines which cells we use. If
1091 ! redblack_cntr is even, we use the even cells and vice versa.
1092 associate(cc => mg%boxes(id)%cc, n => mg_iphi, &
1093 i_eps => mg_iveps)
1094 do k = 1, nc
1095 do j = 1, nc
1096 if (redblack) &
1097 i0 = 2 - iand(ieor(redblack_cntr, k+j), 1)
1098 do i = i0, nc, di
1099 a0 = cc(i, j, k, i_eps)
1100 u(1:2) = cc(i-1:i+1:2, j, k, n)
1101 a(1:2) = cc(i-1:i+1:2, j, k, i_eps)
1102 u(3:4) = cc(i, j-1:j+1:2, k, n)
1103 a(3:4) = cc(i, j-1:j+1:2, k, i_eps)
1104 u(5:6) = cc(i, j, k-1:k+1:2, n)
1105 a(5:6) = cc(i, j, k-1:k+1:2, i_eps)
1106 c(:) = 2 * a0 * a(:) / (a0 + a(:)) * idr2
1107
1108 cc(i, j, k, n) = (sum(c(:) * u(:)) - &
1109 cc(i, j, k, mg_irhs)) / &
1110 (sum(c(:)) + vhelmholtz_lambda)
1111 end do
1112 end do
1113 end do
1114 end associate
1115 end subroutine box_gs_vhelmh
1116
1117 !> Perform Helmholtz operator on a box
1118 subroutine box_vhelmh(mg, id, nc, i_out)
1119 type(mg_t), intent(inout) :: mg
1120 integer, intent(in) :: id
1121 integer, intent(in) :: nc
1122 integer, intent(in) :: i_out !< Index of variable to store Helmholtz in
1123 integer :: i, j, k
1124 real(dp) :: idr2(2*3), a0, u0, u(2*3), a(2*3)
1125
1126 ! Duplicate 1/dr^2 array to multiply neighbor values
1127 idr2(1:2*3:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1128 idr2(2:2*3:2) = idr2(1:2*3:2)
1129
1130 associate(cc => mg%boxes(id)%cc, n => mg_iphi, &
1131 i_eps => mg_iveps)
1132 do k = 1, nc
1133 do j = 1, nc
1134 do i = 1, nc
1135 u0 = cc(i, j, k, n)
1136 a0 = cc(i, j, k, i_eps)
1137 u(1:2) = cc(i-1:i+1:2, j, k, n)
1138 u(3:4) = cc(i, j-1:j+1:2, k, n)
1139 u(5:6) = cc(i, j, k-1:k+1:2, n)
1140 a(1:2) = cc(i-1:i+1:2, j, k, i_eps)
1141 a(3:4) = cc(i, j-1:j+1:2, k, i_eps)
1142 a(5:6) = cc(i, j, k-1:k+1:2, i_eps)
1143
1144 cc(i, j, k, i_out) = sum(2 * idr2 * &
1145 a0*a(:)/(a0 + a(:)) * (u(:) - u0)) - &
1147 end do
1148 end do
1149 end do
1150 end associate
1151 end subroutine box_vhelmh
1152
1153 !! File ../src/m_build_tree.f90
1154
1155 subroutine mg_build_rectangle(mg, domain_size, box_size, dx, r_min, &
1156 periodic, n_finer)
1157 type(mg_t), intent(inout) :: mg
1158 integer, intent(in) :: domain_size(3)
1159 integer, intent(in) :: box_size
1160 real(dp), intent(in) :: dx(3)
1161 real(dp), intent(in) :: r_min(3)
1162 logical, intent(in) :: periodic(3)
1163 integer, intent(in) :: n_finer
1164 integer :: i, j, k, lvl, n, id, nx(3), ib, nb(3)
1165 integer :: boxes_per_dim(3, mg_lvl_lo:1)
1166 integer :: periodic_offset(3)
1167
1168 if (modulo(box_size, 2) /= 0) &
1169 error stop "box_size should be even"
1170 if (any(modulo(domain_size, box_size) /= 0)) &
1171 error stop "box_size does not divide domain_size"
1172
1173 if (all(periodic)) then
1174 mg%subtract_mean = .true.
1175 end if
1176
1177 nx = domain_size
1178 mg%box_size = box_size
1179 mg%box_size_lvl(1) = box_size
1180 mg%domain_size_lvl(:, 1) = domain_size
1181 mg%first_normal_lvl = 1
1182 mg%dr(:, 1) = dx
1183 mg%r_min(:) = r_min
1184 mg%periodic = periodic
1185 boxes_per_dim(:, :) = 0
1186 boxes_per_dim(:, 1) = domain_size / box_size
1187
1188 do lvl = 1, mg_lvl_lo+1, -1
1189 ! For a Gauss-Seidel (non red-black) smoother, we should avoid boxes
1190 ! containing a single cell
1191 if (any(modulo(nx, 2) == 1 .or. nx == mg%coarsest_grid .or. &
1192 (mg%box_size_lvl(lvl) == mg%coarsest_grid .and. &
1193 mg%smoother_type == mg_smoother_gs))) exit
1194
1195 if (all(modulo(nx/mg%box_size_lvl(lvl), 2) == 0)) then
1196 mg%box_size_lvl(lvl-1) = mg%box_size_lvl(lvl)
1197 boxes_per_dim(:, lvl-1) = boxes_per_dim(:, lvl)/2
1198 mg%first_normal_lvl = lvl-1
1199 else
1200 mg%box_size_lvl(lvl-1) = mg%box_size_lvl(lvl)/2
1201 boxes_per_dim(:, lvl-1) = boxes_per_dim(:, lvl)
1202 end if
1203
1204 mg%dr(:, lvl-1) = mg%dr(:, lvl) * 2
1205 nx = nx / 2
1206 mg%domain_size_lvl(:, lvl-1) = nx
1207 end do
1208
1209 mg%lowest_lvl = lvl
1210 mg%highest_lvl = 1
1211
1212 do lvl = 2, mg_lvl_hi
1213 mg%dr(:, lvl) = mg%dr(:, lvl-1) * 0.5_dp
1214 mg%box_size_lvl(lvl) = box_size
1215 mg%domain_size_lvl(:, lvl) = 2 * mg%domain_size_lvl(:, lvl-1)
1216 end do
1217
1218 n = sum(product(boxes_per_dim, dim=1)) + n_finer
1219 allocate(mg%boxes(n))
1220
1221 ! Create lowest level
1222 nx = boxes_per_dim(:, mg%lowest_lvl)
1223 periodic_offset = [nx(1)-1, (nx(2)-1)*nx(1), &
1224 (nx(3)-1) * nx(2) * nx(1)]
1225
1226 do k=1,nx(3); do j=1,nx(2); do i=1,nx(1)
1227 mg%n_boxes = mg%n_boxes + 1
1228 n = mg%n_boxes
1229
1230 mg%boxes(n)%rank = 0
1231 mg%boxes(n)%id = n
1232 mg%boxes(n)%lvl = mg%lowest_lvl
1233 mg%boxes(n)%ix(:) = [i, j, k]
1234 mg%boxes(n)%r_min(:) = r_min + (mg%boxes(n)%ix(:) - 1) * &
1235 mg%box_size_lvl(mg%lowest_lvl) * mg%dr(:, mg%lowest_lvl)
1236 mg%boxes(n)%dr(:) = mg%dr(:, mg%lowest_lvl)
1237 mg%boxes(n)%parent = mg_no_box
1238 mg%boxes(n)%children(:) = mg_no_box
1239
1240 ! Set default neighbors
1241 mg%boxes(n)%neighbors(:) = [n-1, n+1, n-nx(1), n+nx(1), &
1242 n-nx(1)*nx(2), n+nx(1)*nx(2)]
1243
1244 ! Handle boundaries
1245 nb=[i, j, k]
1246 do ib=1,3
1247 if(nb(ib)==1 .and. .not.periodic(ib)) then
1248 mg%boxes(n)%neighbors(ib*2-1) = mg_physical_boundary
1249 end if
1250 if(nb(ib)==1 .and. periodic(ib)) then
1251 mg%boxes(n)%neighbors(ib*2-1) = n + periodic_offset(ib)
1252 end if
1253 if(nb(ib)==nx(ib) .and. .not.periodic(ib)) then
1254 mg%boxes(n)%neighbors(ib*2) = mg_physical_boundary
1255 end if
1256 if(nb(ib)==nx(ib) .and. periodic(ib)) then
1257 mg%boxes(n)%neighbors(ib*2) = n - periodic_offset(ib)
1258 end if
1259 end do
1260 !where ([i, j, k] == 1 .and. .not. periodic)
1261 ! mg%boxes(n)%neighbors(1:mg_num_neighbors:2) = &
1262 ! mg_physical_boundary
1263 !end where
1264 !where ([i, j, k] == 1 .and. periodic)
1265 ! mg%boxes(n)%neighbors(1:mg_num_neighbors:2) = &
1266 ! n + periodic_offset
1267 !end where
1268
1269 !where ([i, j, k] == nx .and. .not. periodic)
1270 ! mg%boxes(n)%neighbors(2:mg_num_neighbors:2) = &
1271 ! mg_physical_boundary
1272 !end where
1273 !where ([i, j, k] == nx .and. periodic)
1274 ! mg%boxes(n)%neighbors(2:mg_num_neighbors:2) = &
1275 ! n - periodic_offset
1276 !end where
1277 end do; end do; end do
1278
1279 mg%lvls(mg%lowest_lvl)%ids = [(n, n=1, mg%n_boxes)]
1280
1281 ! Add higher levels
1282 do lvl = mg%lowest_lvl, 0
1283 if (mg%box_size_lvl(lvl+1) == mg%box_size_lvl(lvl)) then
1284 do i = 1, size(mg%lvls(lvl)%ids)
1285 id = mg%lvls(lvl)%ids(i)
1286 call mg_add_children(mg, id)
1287 end do
1288
1289 call mg_set_leaves_parents(mg%boxes, mg%lvls(lvl))
1290 call mg_set_next_level_ids(mg, lvl)
1291 call mg_set_neighbors_lvl(mg, lvl+1)
1292 else
1293 do i = 1, size(mg%lvls(lvl)%ids)
1294 id = mg%lvls(lvl)%ids(i)
1295 call add_single_child(mg, id, size(mg%lvls(lvl)%ids))
1296 end do
1297
1298 call mg_set_leaves_parents(mg%boxes, mg%lvls(lvl))
1299 call mg_set_next_level_ids(mg, lvl)
1300 end if
1301 end do
1302
1303 call mg_set_leaves_parents(mg%boxes, mg%lvls(1))
1304
1305 ! No refinement boundaries
1306 do lvl = mg%lowest_lvl, 1
1307 if (allocated(mg%lvls(lvl)%ref_bnds)) &
1308 deallocate(mg%lvls(lvl)%ref_bnds)
1309 allocate(mg%lvls(lvl)%ref_bnds(0))
1310 end do
1311
1312 mg%tree_created = .true.
1313 end subroutine mg_build_rectangle
1314
1315 subroutine mg_set_neighbors_lvl(mg, lvl)
1316 type(mg_t), intent(inout) :: mg
1317 integer, intent(in) :: lvl
1318 integer :: i, id
1319
1320 do i = 1, size(mg%lvls(lvl)%ids)
1321 id = mg%lvls(lvl)%ids(i)
1322 call set_neighbs(mg%boxes, id)
1323 end do
1324 end subroutine mg_set_neighbors_lvl
1325
1326 subroutine mg_set_next_level_ids(mg, lvl)
1327 type(mg_t), intent(inout) :: mg
1328 integer, intent(in) :: lvl
1329 integer :: n, i, id
1330
1331 if (allocated(mg%lvls(lvl+1)%ids)) &
1332 deallocate(mg%lvls(lvl+1)%ids)
1333
1334 ! Set next level ids to children of this level
1335 if (mg%box_size_lvl(lvl+1) == mg%box_size_lvl(lvl)) then
1336 n = mg_num_children * size(mg%lvls(lvl)%parents)
1337 allocate(mg%lvls(lvl+1)%ids(n))
1338
1339 n = mg_num_children
1340 do i = 1, size(mg%lvls(lvl)%parents)
1341 id = mg%lvls(lvl)%parents(i)
1342 mg%lvls(lvl+1)%ids(n*(i-1)+1:n*i) = mg%boxes(id)%children
1343 end do
1344 else
1345 n = size(mg%lvls(lvl)%parents)
1346 allocate(mg%lvls(lvl+1)%ids(n))
1347
1348 n = 1
1349 do i = 1, size(mg%lvls(lvl)%parents)
1350 id = mg%lvls(lvl)%parents(i)
1351 mg%lvls(lvl+1)%ids(i) = mg%boxes(id)%children(1)
1352 end do
1353 end if
1354
1355 end subroutine mg_set_next_level_ids
1356
1357 ! Set the neighbors of id (using their parent)
1358 subroutine set_neighbs(boxes, id)
1359 type(mg_box_t), intent(inout) :: boxes(:)
1360 integer, intent(in) :: id
1361 integer :: nb, nb_id
1362
1363 do nb = 1, mg_num_neighbors
1364 if (boxes(id)%neighbors(nb) == mg_no_box) then
1365 nb_id = find_neighb(boxes, id, nb)
1366 if (nb_id > mg_no_box) then
1367 boxes(id)%neighbors(nb) = nb_id
1368 boxes(nb_id)%neighbors(mg_neighb_rev(nb)) = id
1369 end if
1370 end if
1371 end do
1372 end subroutine set_neighbs
1373
1374 !> Get the id of neighbor nb of boxes(id), through its parent
1375 function find_neighb(boxes, id, nb) result(nb_id)
1376 type(mg_box_t), intent(in) :: boxes(:) !< List with all the boxes
1377 integer, intent(in) :: id !< Box whose neighbor we are looking for
1378 integer, intent(in) :: nb !< Neighbor index
1379 integer :: nb_id, p_id, c_ix, d, old_pid
1380
1381 p_id = boxes(id)%parent
1382 old_pid = p_id
1383 c_ix = mg_ix_to_ichild(boxes(id)%ix)
1384 d = mg_neighb_dim(nb)
1385
1386 ! Check if neighbor is in same direction as ix is (low/high). If so,
1387 ! use neighbor of parent
1388 if (mg_child_low(d, c_ix) .eqv. mg_neighb_low(nb)) then
1389 p_id = boxes(p_id)%neighbors(nb)
1390 end if
1391
1392 ! The child ix of the neighbor is reversed in direction d
1393 nb_id = boxes(p_id)%children(mg_child_rev(c_ix, d))
1394 end function find_neighb
1395
1396 !> Create a list of leaves and a list of parents for a level
1397 subroutine mg_set_leaves_parents(boxes, level)
1398 type(mg_box_t), intent(in) :: boxes(:) !< List of boxes
1399 type(mg_lvl_t), intent(inout) :: level !< Level type which contains the indices of boxes
1400 integer :: i, id, i_leaf, i_parent
1401 integer :: n_parents, n_leaves
1402
1403 n_parents = count(mg_has_children(boxes(level%ids)))
1404 n_leaves = size(level%ids) - n_parents
1405
1406 if (.not. allocated(level%parents)) then
1407 allocate(level%parents(n_parents))
1408 else if (n_parents /= size(level%parents)) then
1409 deallocate(level%parents)
1410 allocate(level%parents(n_parents))
1411 end if
1412
1413 if (.not. allocated(level%leaves)) then
1414 allocate(level%leaves(n_leaves))
1415 else if (n_leaves /= size(level%leaves)) then
1416 deallocate(level%leaves)
1417 allocate(level%leaves(n_leaves))
1418 end if
1419
1420 i_leaf = 0
1421 i_parent = 0
1422 do i = 1, size(level%ids)
1423 id = level%ids(i)
1424 if (mg_has_children(boxes(id))) then
1425 i_parent = i_parent + 1
1426 level%parents(i_parent) = id
1427 else
1428 i_leaf = i_leaf + 1
1429 level%leaves(i_leaf) = id
1430 end if
1431 end do
1432 end subroutine mg_set_leaves_parents
1433
1434 !> Create a list of refinement boundaries (from the coarse side)
1435 subroutine mg_set_refinement_boundaries(boxes, level)
1436 type(mg_box_t), intent(in) :: boxes(:)
1437 type(mg_lvl_t), intent(inout) :: level
1438 integer, allocatable :: tmp(:)
1439 integer :: i, id, nb, nb_id, ix
1440
1441 if (allocated(level%ref_bnds)) deallocate(level%ref_bnds)
1442
1443 if (size(level%parents) == 0) then
1444 ! There are no refinement boundaries
1445 allocate(level%ref_bnds(0))
1446 else
1447 allocate(tmp(size(level%leaves)))
1448 ix = 0
1449 do i = 1, size(level%leaves)
1450 id = level%leaves(i)
1451
1452 do nb = 1, mg_num_neighbors
1453 nb_id = boxes(id)%neighbors(nb)
1454 if (nb_id > mg_no_box) then
1455 if (mg_has_children(boxes(nb_id))) then
1456 ix = ix + 1
1457 tmp(ix) = id
1458 exit
1459 end if
1460 end if
1461 end do
1462 end do
1463
1464 allocate(level%ref_bnds(ix))
1465 level%ref_bnds(:) = tmp(1:ix)
1466 end if
1467 end subroutine mg_set_refinement_boundaries
1468
1469 subroutine mg_add_children(mg, id)
1470 type(mg_t), intent(inout) :: mg
1471 integer, intent(in) :: id !< Id of box that gets children
1472 integer :: lvl, i, nb, child_nb(2**(3-1))
1473 integer :: c_ids(mg_num_children)
1474 integer :: c_id, c_ix_base(3)
1475
1476 if (mg%n_boxes + mg_num_children > size(mg%boxes)) then
1477 error stop "mg_add_children: not enough space"
1478 end if
1479
1480 c_ids = [(mg%n_boxes+i, i=1,mg_num_children)]
1481 mg%n_boxes = mg%n_boxes + mg_num_children
1482 mg%boxes(id)%children = c_ids
1483 c_ix_base = 2 * mg%boxes(id)%ix - 1
1484 lvl = mg%boxes(id)%lvl+1
1485
1486 do i = 1, mg_num_children
1487 c_id = c_ids(i)
1488 mg%boxes(c_id)%rank = mg%boxes(id)%rank
1489 mg%boxes(c_id)%ix = c_ix_base + mg_child_dix(:, i)
1490 mg%boxes(c_id)%lvl = lvl
1491 mg%boxes(c_id)%parent = id
1492 mg%boxes(c_id)%children = mg_no_box
1493 mg%boxes(c_id)%neighbors = mg_no_box
1494 mg%boxes(c_id)%r_min = mg%boxes(id)%r_min + &
1495 mg%dr(:, lvl) * mg_child_dix(:, i) * mg%box_size
1496 mg%boxes(c_id)%dr(:) = mg%dr(:, lvl)
1497 end do
1498
1499 ! Set boundary conditions at children
1500 do nb = 1, mg_num_neighbors
1501 if (mg%boxes(id)%neighbors(nb) < mg_no_box) then
1502 child_nb = c_ids(mg_child_adj_nb(:, nb)) ! Neighboring children
1503 mg%boxes(child_nb)%neighbors(nb) = mg%boxes(id)%neighbors(nb)
1504 end if
1505 end do
1506 end subroutine mg_add_children
1507
1508 subroutine add_single_child(mg, id, n_boxes_lvl)
1509 type(mg_t), intent(inout) :: mg
1510 integer, intent(in) :: id !< Id of box that gets children
1511 integer, intent(in) :: n_boxes_lvl
1512 integer :: lvl, c_id
1513
1514 c_id = mg%n_boxes + 1
1515 mg%n_boxes = mg%n_boxes + 1
1516 mg%boxes(id)%children(1) = c_id
1517 lvl = mg%boxes(id)%lvl+1
1518
1519 mg%boxes(c_id)%rank = mg%boxes(id)%rank
1520 mg%boxes(c_id)%ix = mg%boxes(id)%ix
1521 mg%boxes(c_id)%lvl = lvl
1522 mg%boxes(c_id)%parent = id
1523 mg%boxes(c_id)%children = mg_no_box
1524 where (mg%boxes(id)%neighbors == mg_physical_boundary)
1525 mg%boxes(c_id)%neighbors = mg%boxes(id)%neighbors
1526 elsewhere
1527 mg%boxes(c_id)%neighbors = mg%boxes(id)%neighbors + n_boxes_lvl
1528 end where
1529 mg%boxes(c_id)%r_min = mg%boxes(id)%r_min
1530 mg%boxes(c_id)%dr(:) = mg%dr(:, lvl)
1531
1532 end subroutine add_single_child
1533
1534 !! File ../src/m_load_balance.f90
1535
1536 !> Load balance all boxes in the multigrid tree, by simply distributing the
1537 !> load per grid level. This method will only work well for uniform grids.
1538 !>
1539 !> Note that in a typical application the load balancing of the leaves is
1540 !> already determined, then mg_load_balance_parents can be used.
1542 type(mg_t), intent(inout) :: mg
1543 integer :: i, id, lvl, single_cpu_lvl
1544 integer :: work_left, my_work, i_cpu
1545
1546 ! Up to this level, all boxes have to be on a single processor because they
1547 ! have a different size and the communication routines do not support this
1548 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1549
1550 do lvl = mg%lowest_lvl, single_cpu_lvl
1551 do i = 1, size(mg%lvls(lvl)%ids)
1552 id = mg%lvls(lvl)%ids(i)
1553 mg%boxes(id)%rank = 0
1554 end do
1555 end do
1556
1557 ! Distribute the boxes equally. Due to the way the mesh is constructed, the
1558 ! mg%lvls(lvl)%ids array already contains a Morton-like ordering.
1559 do lvl = single_cpu_lvl+1, mg%highest_lvl
1560 work_left = size(mg%lvls(lvl)%ids)
1561 my_work = 0
1562 i_cpu = 0
1563
1564 do i = 1, size(mg%lvls(lvl)%ids)
1565 if ((mg%n_cpu - i_cpu - 1) * my_work >= work_left) then
1566 i_cpu = i_cpu + 1
1567 my_work = 0
1568 end if
1569
1570 my_work = my_work + 1
1571 work_left = work_left - 1
1572
1573 id = mg%lvls(lvl)%ids(i)
1574 mg%boxes(id)%rank = i_cpu
1575 end do
1576 end do
1577
1578 do lvl = mg%lowest_lvl, mg%highest_lvl
1579 call update_lvl_info(mg, mg%lvls(lvl))
1580 end do
1581
1582 end subroutine mg_load_balance_simple
1583
1584 !> Load balance all boxes in the multigrid tree. Compared to
1585 !> mg_load_balance_simple, this method does a better job of setting the ranks
1586 !> of parent boxes
1587 !>
1588 !> Note that in a typical application the load balancing of the leaves is
1589 !> already determined, then mg_load_balance_parents can be used.
1590 subroutine mg_load_balance(mg)
1591 type(mg_t), intent(inout) :: mg
1592 integer :: i, id, lvl, single_cpu_lvl
1593 integer :: work_left, my_work(0:mg%n_cpu), i_cpu
1594 integer :: c_ids(mg_num_children)
1595 integer :: c_ranks(mg_num_children)
1596 integer :: coarse_rank
1597
1598 ! Up to this level, all boxes have to be on a single processor because they
1599 ! have a different size and the communication routines do not support this
1600 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1601
1602 ! Distribute the boxes equally. Due to the way the mesh is constructed, the
1603 ! mg%lvls(lvl)%ids array already contains a Morton-like ordering.
1604 do lvl = mg%highest_lvl, single_cpu_lvl+1, -1
1605 ! For parents determine the rank based on their child ranks
1606 my_work(:) = 0
1607
1608 do i = 1, size(mg%lvls(lvl)%parents)
1609 id = mg%lvls(lvl)%parents(i)
1610
1611 c_ids = mg%boxes(id)%children
1612 c_ranks = mg%boxes(c_ids)%rank
1613 i_cpu = most_popular(c_ranks, my_work, mg%n_cpu)
1614 mg%boxes(id)%rank = i_cpu
1615 my_work(i_cpu) = my_work(i_cpu) + 1
1616 end do
1617
1618 work_left = size(mg%lvls(lvl)%leaves)
1619 i_cpu = 0
1620
1621 do i = 1, size(mg%lvls(lvl)%leaves)
1622 ! Skip this CPU if it already has enough work
1623 if ((mg%n_cpu - i_cpu - 1) * my_work(i_cpu) >= &
1624 work_left + sum(my_work(i_cpu+1:))) then
1625 i_cpu = i_cpu + 1
1626 end if
1627
1628 my_work(i_cpu) = my_work(i_cpu) + 1
1629 work_left = work_left - 1
1630
1631 id = mg%lvls(lvl)%leaves(i)
1632 mg%boxes(id)%rank = i_cpu
1633 end do
1634 end do
1635
1636 ! Determine most popular CPU for coarse grids
1637 if (single_cpu_lvl < mg%highest_lvl) then
1638 coarse_rank = most_popular(mg%boxes(&
1639 mg%lvls(single_cpu_lvl+1)%ids)%rank, my_work, mg%n_cpu)
1640 else
1641 coarse_rank = 0
1642 end if
1643
1644 do lvl = mg%lowest_lvl, single_cpu_lvl
1645 do i = 1, size(mg%lvls(lvl)%ids)
1646 id = mg%lvls(lvl)%ids(i)
1647 mg%boxes(id)%rank = coarse_rank
1648 end do
1649 end do
1650
1651 do lvl = mg%lowest_lvl, mg%highest_lvl
1652 call update_lvl_info(mg, mg%lvls(lvl))
1653 end do
1654
1655 end subroutine mg_load_balance
1656
1657 !> Load balance the parents (non-leafs). Assign them to the rank that has most
1658 !> children.
1660 type(mg_t), intent(inout) :: mg
1661 integer :: i, id, lvl
1662 integer :: c_ids(mg_num_children)
1663 integer :: c_ranks(mg_num_children)
1664 integer :: single_cpu_lvl, coarse_rank
1665 integer :: my_work(0:mg%n_cpu), i_cpu
1666
1667 ! Up to this level, all boxes have to be on a single processor because they
1668 ! have a different size and the communication routines do not support this
1669 single_cpu_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
1670
1671 do lvl = mg%highest_lvl-1, single_cpu_lvl+1, -1
1672 my_work(:) = 0
1673
1674 ! Determine amount of work for the leaves
1675 do i = 1, size(mg%lvls(lvl)%leaves)
1676 id = mg%lvls(lvl)%leaves(i)
1677 i_cpu = mg%boxes(id)%rank
1678 my_work(i_cpu) = my_work(i_cpu) + 1
1679 end do
1680
1681 do i = 1, size(mg%lvls(lvl)%parents)
1682 id = mg%lvls(lvl)%parents(i)
1683
1684 c_ids = mg%boxes(id)%children
1685 c_ranks = mg%boxes(c_ids)%rank
1686 i_cpu = most_popular(c_ranks, my_work, mg%n_cpu)
1687 mg%boxes(id)%rank = i_cpu
1688 my_work(i_cpu) = my_work(i_cpu) + 1
1689 end do
1690
1691 end do
1692
1693 ! Determine most popular CPU for coarse grids
1694 if (single_cpu_lvl < mg%highest_lvl) then
1695 coarse_rank = most_popular(mg%boxes(&
1696 mg%lvls(single_cpu_lvl+1)%ids)%rank, my_work, mg%n_cpu)
1697 else
1698 coarse_rank = 0
1699 end if
1700
1701 do lvl = mg%lowest_lvl, single_cpu_lvl
1702 do i = 1, size(mg%lvls(lvl)%ids)
1703 id = mg%lvls(lvl)%ids(i)
1704 mg%boxes(id)%rank = coarse_rank
1705 end do
1706 end do
1707
1708 do lvl = mg%lowest_lvl, mg%highest_lvl
1709 call update_lvl_info(mg, mg%lvls(lvl))
1710 end do
1711
1712 end subroutine mg_load_balance_parents
1713
1714 !> Determine most popular rank in the list. In case of ties, assign the rank
1715 !> with the least work.
1716 pure integer function most_popular(list, work, n_cpu)
1717 integer, intent(in) :: list(:) !< List of MPI ranks
1718 integer, intent(in) :: n_cpu
1719 integer, intent(in) :: work(0:n_cpu-1) !< Existing work per rank
1720 integer :: i, best_count, current_count
1721 integer :: my_work, best_work
1722
1723 best_count = 0
1724 best_work = 0
1725 most_popular = -1
1726
1727 do i = 1, size(list)
1728 current_count = count(list == list(i))
1729 my_work = work(list(i))
1730
1731 ! In case of ties, select task with lowest work
1732 if (current_count > best_count .or. &
1733 current_count == best_count .and. my_work < best_work) then
1734 best_count = current_count
1735 best_work = my_work
1736 most_popular = list(i)
1737 end if
1738 end do
1739
1740 end function most_popular
1741
1742 subroutine update_lvl_info(mg, lvl)
1743 type(mg_t), intent(inout) :: mg
1744 type(mg_lvl_t), intent(inout) :: lvl
1745
1746 lvl%my_ids = pack(lvl%ids, &
1747 mg%boxes(lvl%ids)%rank == mg%my_rank)
1748 lvl%my_leaves = pack(lvl%leaves, &
1749 mg%boxes(lvl%leaves)%rank == mg%my_rank)
1750 lvl%my_parents = pack(lvl%parents, &
1751 mg%boxes(lvl%parents)%rank == mg%my_rank)
1752 lvl%my_ref_bnds = pack(lvl%ref_bnds, &
1753 mg%boxes(lvl%ref_bnds)%rank == mg%my_rank)
1754 end subroutine update_lvl_info
1755
1756 !! File ../src/m_vlaplacian.f90
1757
1759 type(mg_t), intent(inout) :: mg
1760
1761 if (mg%n_extra_vars == 0 .and. mg%is_allocated) then
1762 error stop "vlaplacian_set_methods: mg%n_extra_vars == 0"
1763 else
1764 mg%n_extra_vars = max(1, mg%n_extra_vars)
1765 end if
1766
1767 mg%subtract_mean = .false.
1768
1769 ! Use Neumann zero boundary conditions for the variable coefficient, since
1770 ! it is needed in ghost cells.
1771 mg%bc(:, mg_iveps)%bc_type = mg_bc_neumann
1772 mg%bc(:, mg_iveps)%bc_value = 0.0_dp
1773
1774 select case (mg%geometry_type)
1775 case (mg_cartesian)
1776 mg%box_op => box_vlpl
1777
1778 ! TODO: test whether a custom prolongation works better
1779 ! mg%box_prolong => vlpl_prolong
1780
1781 select case (mg%smoother_type)
1783 mg%box_smoother => box_gs_vlpl
1784 case default
1785 error stop "vlaplacian_set_methods: unsupported smoother type"
1786 end select
1787
1788 case default
1789 error stop "vlaplacian_set_methods: unsupported geometry"
1790 end select
1791
1792 end subroutine vlaplacian_set_methods
1793
1794 !> Perform Gauss-Seidel relaxation on box for a Laplacian operator
1795 subroutine box_gs_vlpl(mg, id, nc, redblack_cntr)
1796 type(mg_t), intent(inout) :: mg
1797 integer, intent(in) :: id
1798 integer, intent(in) :: nc
1799 integer, intent(in) :: redblack_cntr !< Iteration counter
1800 integer :: i, j, k, i0, di
1801 logical :: redblack
1802 real(dp) :: idr2(2*3), u(2*3)
1803 real(dp) :: a0, a(2*3), c(2*3)
1804
1805 ! Duplicate 1/dr^2 array to multiply neighbor values
1806 idr2(1:2*3:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1807 idr2(2:2*3:2) = idr2(1:2*3:2)
1808 i0 = 1
1809
1810 redblack = (mg%smoother_type == mg_smoother_gsrb)
1811 if (redblack) then
1812 di = 2
1813 else
1814 di = 1
1815 end if
1816
1817 ! The parity of redblack_cntr determines which cells we use. If
1818 ! redblack_cntr is even, we use the even cells and vice versa.
1819 associate(cc => mg%boxes(id)%cc, n => mg_iphi, &
1820 i_eps => mg_iveps)
1821 do k = 1, nc
1822 do j = 1, nc
1823 if (redblack) &
1824 i0 = 2 - iand(ieor(redblack_cntr, k+j), 1)
1825 do i = i0, nc, di
1826 a0 = cc(i, j, k, i_eps)
1827 u(1:2) = cc(i-1:i+1:2, j, k, n)
1828 a(1:2) = cc(i-1:i+1:2, j, k, i_eps)
1829 u(3:4) = cc(i, j-1:j+1:2, k, n)
1830 a(3:4) = cc(i, j-1:j+1:2, k, i_eps)
1831 u(5:6) = cc(i, j, k-1:k+1:2, n)
1832 a(5:6) = cc(i, j, k-1:k+1:2, i_eps)
1833 c(:) = 2 * a0 * a(:) / (a0 + a(:)) * idr2
1834
1835 cc(i, j, k, n) = (sum(c(:) * u(:)) - &
1836 cc(i, j, k, mg_irhs)) / sum(c(:))
1837 end do
1838 end do
1839 end do
1840 end associate
1841 end subroutine box_gs_vlpl
1842
1843 !> Perform Laplacian operator on a box
1844 subroutine box_vlpl(mg, id, nc, i_out)
1845 type(mg_t), intent(inout) :: mg
1846 integer, intent(in) :: id
1847 integer, intent(in) :: nc
1848 integer, intent(in) :: i_out !< Index of variable to store Laplacian in
1849 integer :: i, j, k
1850 real(dp) :: idr2(2*3), a0, u0, u(2*3), a(2*3)
1851
1852 ! Duplicate 1/dr^2 array to multiply neighbor values
1853 idr2(1:2*3:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
1854 idr2(2:2*3:2) = idr2(1:2*3:2)
1855
1856 associate(cc => mg%boxes(id)%cc, n => mg_iphi, &
1857 i_eps => mg_iveps)
1858 do k = 1, nc
1859 do j = 1, nc
1860 do i = 1, nc
1861 u0 = cc(i, j, k, n)
1862 a0 = cc(i, j, k, i_eps)
1863 u(1:2) = cc(i-1:i+1:2, j, k, n)
1864 u(3:4) = cc(i, j-1:j+1:2, k, n)
1865 u(5:6) = cc(i, j, k-1:k+1:2, n)
1866 a(1:2) = cc(i-1:i+1:2, j, k, i_eps)
1867 a(3:4) = cc(i, j-1:j+1:2, k, i_eps)
1868 a(5:6) = cc(i, j, k-1:k+1:2, i_eps)
1869
1870 cc(i, j, k, i_out) = sum(2 * idr2 * &
1871 a0*a(:)/(a0 + a(:)) * (u(:) - u0))
1872 end do
1873 end do
1874 end do
1875 end associate
1876 end subroutine box_vlpl
1877
1878 !> Prolong from a parent to a child with index offset dix. This method could
1879 !> sometimes work better than the default prolongation, which does not take
1880 !> the variation in epsilon into account
1881! subroutine vlpl_prolong(mg, p_id, dix, nc, iv, fine)
1882! type(mg_t), intent(inout) :: mg
1883! integer, intent(in) :: p_id !< Id of parent
1884! integer, intent(in) :: dix(3) !< Offset of child in parent grid
1885! integer, intent(in) :: nc !< Child grid size
1886! integer, intent(in) :: iv !< Prolong from this variable
1887! real(dp), intent(out) :: fine(nc, nc, nc) !< Prolonged values
1888
1889! integer :: i, j, k, hnc
1890! #if 3 == 2
1891! integer :: ic, jc
1892! real(dp) :: f0, flx, fhx, fly, fhy
1893! real(dp) :: c0, clx, chx, cly, chy
1894! #elif 3 == 3
1895! integer :: ic, jc, kc
1896! real(dp) :: f0, flx, fhx, fly, fhy, flz, fhz
1897! real(dp) :: c0, clx, chx, cly, chy, clz, chz
1898! #endif
1899
1900! hnc = nc/2
1901
1902! associate (crs => mg%boxes(p_id)%cc, &
1903! i_eps => mg_iveps)
1904! #if 3 == 2
1905! do j = 1, hnc
1906! jc = j + dix(2)
1907! do i = 1, hnc
1908! ic = i + dix(1)
1909
1910! f0 = crs(ic, jc, iv)
1911! flx = crs(ic-1, jc, iv)
1912! fhx = crs(ic+1, jc, iv)
1913! fly = crs(ic, jc-1, iv)
1914! fhy = crs(ic, jc+1, iv)
1915
1916! c0 = 2 * crs(ic, jc, i_eps)
1917! clx = crs(ic-1, jc, i_eps)
1918! chx = crs(ic+1, jc, i_eps)
1919! cly = crs(ic, jc-1, i_eps)
1920! chy = crs(ic, jc+1, i_eps)
1921
1922! fine(2*i-1, 2*j-1) = (f0*c0 + flx*clx + fly*cly) / &
1923! (c0 + clx + cly)
1924! fine(2*i , 2*j-1) = (f0*c0 + fhx*chx + fly*cly) / &
1925! (c0 + chx + cly)
1926! fine(2*i-1, 2*j) = (f0*c0 + flx*clx + fhy*chy) / &
1927! (c0 + clx + chy)
1928! fine(2*i , 2*j) = (f0*c0 + fhx*chx + fhy*chy) / &
1929! (c0 + chx + chy)
1930! end do
1931! end do
1932! #elif 3 == 3
1933! do k = 1, hnc
1934! kc = k + dix(3)
1935! do j = 1, hnc
1936! jc = j + dix(2)
1937! do i = 1, hnc
1938! ic = i + dix(1)
1939
1940! f0 = crs(ic, jc, kc, iv)
1941! flx = crs(ic-1, jc, kc, iv)
1942! fhx = crs(ic+1, jc, kc, iv)
1943! fly = crs(ic, jc-1, kc, iv)
1944! fhy = crs(ic, jc+1, kc, iv)
1945! flz = crs(ic, jc, kc-1, iv)
1946! fhz = crs(ic, jc, kc+1, iv)
1947
1948! c0 = crs(ic, jc, kc, i_eps)
1949! clx = crs(ic-1, jc, kc, i_eps)
1950! chx = crs(ic+1, jc, kc, i_eps)
1951! cly = crs(ic, jc-1, kc, i_eps)
1952! chy = crs(ic, jc+1, kc, i_eps)
1953! clz = crs(ic, jc, kc-1, i_eps)
1954! chz = crs(ic, jc, kc+1, i_eps)
1955
1956! fine(2*i-1, 2*j-1, 2*k-1) = (f0*c0 + flx*clx + fly*cly + flz*clz) / &
1957! (c0 + clx + cly + clz)
1958! fine(2*i, 2*j-1, 2*k-1) = (f0*c0 + fhx*chx + fly*cly + flz*clz) / &
1959! (c0 + chx + cly + clz)
1960! fine(2*i-1, 2*j, 2*k-1) = (f0*c0 + flx*clx + fhy*chy + flz*clz) / &
1961! (c0 + clx + chy + clz)
1962! fine(2*i, 2*j, 2*k-1) = (f0*c0 + fhx*chx + fhy*chy + flz*clz) / &
1963! (c0 + chx + chy + clz)
1964! fine(2*i-1, 2*j-1, 2*k) = (f0*c0 + flx*clx + fly*cly + fhz*chz) / &
1965! (c0 + clx + cly + chz)
1966! fine(2*i, 2*j-1, 2*k) = (f0*c0 + fhx*chx + fly*cly + fhz*chz) / &
1967! (c0 + chx + cly + chz)
1968! fine(2*i-1, 2*j, 2*k) = (f0*c0 + flx*clx + fhy*chy + fhz*chz) / &
1969! (c0 + clx + chy + chz)
1970! fine(2*i, 2*j, 2*k) = (f0*c0 + fhx*chx + fhy*chy + fhz*chz) / &
1971! (c0 + chx + chy + chz)
1972! end do
1973! end do
1974! end do
1975! #endif
1976! end associate
1977! end subroutine vlpl_prolong
1978
1979 !! File ../src/m_communication.f90
1980
1981 !> Initialize MPI if needed, and store MPI information
1982 subroutine mg_comm_init(mg, comm)
1983 use mpi
1984 type(mg_t), intent(inout) :: mg
1985 !> MPI communicator (default: MPI_COMM_WORLD)
1986 integer, intent(in), optional :: comm
1987 integer :: ierr
1988 logical :: initialized
1989
1990 call mpi_initialized(initialized, ierr)
1991 if (.not. initialized) then
1992 call mpi_init(ierr)
1993 end if
1994
1995 if (present(comm)) then
1996 mg%comm = comm
1997 else
1998 mg%comm = mpi_comm_world
1999 end if
2000
2001 call mpi_comm_rank(mg%comm, mg%my_rank, ierr)
2002 call mpi_comm_size(mg%comm, mg%n_cpu, ierr)
2003 end subroutine mg_comm_init
2004
2005 subroutine sort_and_transfer_buffers(mg, dsize)
2006 use mpi
2007 type(mg_t), intent(inout) :: mg
2008 integer, intent(in) :: dsize
2009 integer :: i, n_send, n_recv
2010 integer :: send_req(mg%n_cpu)
2011 integer :: recv_req(mg%n_cpu)
2012 integer :: ierr
2013
2014 n_send = 0
2015 n_recv = 0
2016
2017 do i = 0, mg%n_cpu - 1
2018 if (mg%buf(i)%i_send > 0) then
2019 n_send = n_send + 1
2020 call sort_sendbuf(mg%buf(i), dsize)
2021 call mpi_isend(mg%buf(i)%send, mg%buf(i)%i_send, mpi_double, &
2022 i, 0, mg%comm, send_req(n_send), ierr)
2023 end if
2024 if (mg%buf(i)%i_recv > 0) then
2025 n_recv = n_recv + 1
2026 call mpi_irecv(mg%buf(i)%recv, mg%buf(i)%i_recv, mpi_double, &
2027 i, 0, mg%comm, recv_req(n_recv), ierr)
2028 end if
2029 end do
2030
2031 call mpi_waitall(n_recv, recv_req(1:n_recv), mpi_statuses_ignore, ierr)
2032 call mpi_waitall(n_send, send_req(1:n_send), mpi_statuses_ignore, ierr)
2033
2034 end subroutine sort_and_transfer_buffers
2035
2036 !> Sort send buffers according to the idbuf array
2037 subroutine sort_sendbuf(gc, dsize)
2038
2039 type(mg_buf_t), intent(inout) :: gc
2040 integer, intent(in) :: dsize !< Size of send buffer elements
2041 integer :: ix_sort(gc%i_ix)
2042 real(dp) :: buf_cpy(gc%i_send)
2043 integer :: i, j, n
2044
2045 call mrgrnk(gc%ix(1:gc%i_ix), ix_sort)
2046
2047 buf_cpy = gc%send(1:gc%i_send)
2048
2049 do n = 1, gc%i_ix
2050 i = (n-1) * dsize
2051 j = (ix_sort(n)-1) * dsize
2052 gc%send(i+1:i+dsize) = buf_cpy(j+1:j+dsize)
2053 end do
2054 gc%ix(1:gc%i_ix) = gc%ix(ix_sort)
2055
2056 end subroutine sort_sendbuf
2057
2058 !! File ../src/m_ghost_cells.f90
2059
2060 !> Specify minimum buffer size (per process) for communication
2061 subroutine mg_ghost_cell_buffer_size(mg, n_send, n_recv, dsize)
2062 type(mg_t), intent(inout) :: mg
2063 integer, intent(out) :: n_send(0:mg%n_cpu-1)
2064 integer, intent(out) :: n_recv(0:mg%n_cpu-1)
2065 integer, intent(out) :: dsize
2066 integer :: i, id, lvl, nc
2067
2068 allocate(mg%comm_ghostcell%n_send(0:mg%n_cpu-1, &
2069 mg%first_normal_lvl:mg%highest_lvl))
2070 allocate(mg%comm_ghostcell%n_recv(0:mg%n_cpu-1, &
2071 mg%first_normal_lvl:mg%highest_lvl))
2072
2073 dsize = mg%box_size**(3-1)
2074
2075 do lvl = mg%first_normal_lvl, mg%highest_lvl
2076 nc = mg%box_size_lvl(lvl)
2077 mg%buf(:)%i_send = 0
2078 mg%buf(:)%i_recv = 0
2079 mg%buf(:)%i_ix = 0
2080
2081 do i = 1, size(mg%lvls(lvl)%my_ids)
2082 id = mg%lvls(lvl)%my_ids(i)
2083 call buffer_ghost_cells(mg, id, nc, 1, dry_run=.true.)
2084 end do
2085
2086 if (lvl > 1) then
2087 do i = 1, size(mg%lvls(lvl-1)%my_ref_bnds)
2088 id = mg%lvls(lvl-1)%my_ref_bnds(i)
2089 call buffer_refinement_boundaries(mg, id, nc, 1, dry_run=.true.)
2090 end do
2091 end if
2092
2093 ! Set ghost cells to received data
2094 mg%buf(:)%i_recv = 0
2095 do i = 1, size(mg%lvls(lvl)%my_ids)
2096 id = mg%lvls(lvl)%my_ids(i)
2097 call set_ghost_cells(mg, id, nc, 1, dry_run=.true.)
2098 end do
2099
2100 mg%comm_ghostcell%n_send(:, lvl) = mg%buf(:)%i_send/dsize
2101 mg%comm_ghostcell%n_recv(:, lvl) = mg%buf(:)%i_recv/dsize
2102 end do
2103
2104 n_send = maxval(mg%comm_ghostcell%n_send, dim=2)
2105 n_recv = maxval(mg%comm_ghostcell%n_recv, dim=2)
2106 end subroutine mg_ghost_cell_buffer_size
2107
2108 !> Store boundary conditions for the solution variable, this can speed up
2109 !> calculations if the same boundary conditions are re-used.
2110 subroutine mg_phi_bc_store(mg)
2111 type(mg_t), intent(inout) :: mg
2112 integer :: lvl, nc
2113
2114 nc = mg%box_size
2115
2116 do lvl = mg%lowest_lvl, mg%highest_lvl
2117 nc = mg%box_size_lvl(lvl)
2118 call mg_phi_bc_store_lvl(mg, lvl, nc)
2119 end do
2120
2121 mg%phi_bc_data_stored = .true.
2122 end subroutine mg_phi_bc_store
2123
2124 subroutine mg_phi_bc_store_lvl(mg, lvl, nc)
2125 type(mg_t), intent(inout) :: mg
2126 integer, intent(in) :: lvl
2127 integer, intent(in) :: nc
2128 real(dp) :: bc(nc, nc)
2129 integer :: i, id, nb, nb_id, bc_type
2130
2131 do i = 1, size(mg%lvls(lvl)%my_ids)
2132 id = mg%lvls(lvl)%my_ids(i)
2133 do nb = 1, mg_num_neighbors
2134 nb_id = mg%boxes(id)%neighbors(nb)
2135 if (nb_id < mg_no_box) then
2136 ! Physical boundary
2137 if (associated(mg%bc(nb, mg_iphi)%boundary_cond)) then
2138 call mg%bc(nb, mg_iphi)%boundary_cond(mg%boxes(id), nc, &
2139 mg_iphi, nb, bc_type, bc)
2140 else
2141 bc_type = mg%bc(nb, mg_iphi)%bc_type
2142 bc = mg%bc(nb, mg_iphi)%bc_value
2143 end if
2144
2145 ! Store the boundary condition type. This is not globally set in
2146 ! the tree, but all negative values are treated the same in
2147 ! other parts of the code
2148 mg%boxes(id)%neighbors(nb) = bc_type
2149
2150 ! Store ghost cell data in the right-hand side
2151 call box_set_gc(mg%boxes(id), nb, nc, mg_irhs, bc)
2152 end if
2153 end do
2154 end do
2155 end subroutine mg_phi_bc_store_lvl
2156
2157 !> Fill ghost cells at all grid levels
2158 subroutine mg_fill_ghost_cells(mg, iv)
2159 type(mg_t) :: mg
2160 integer, intent(in) :: iv !< Index of variable
2161 integer :: lvl
2162
2163 do lvl = mg%lowest_lvl, mg%highest_lvl
2164 call mg_fill_ghost_cells_lvl(mg, lvl, iv)
2165 end do
2166 end subroutine mg_fill_ghost_cells
2167
2168 !> Fill ghost cells at a grid level
2169 subroutine mg_fill_ghost_cells_lvl(mg, lvl, iv)
2170
2171 type(mg_t) :: mg
2172 integer, intent(in) :: lvl
2173 integer, intent(in) :: iv !< Index of variable
2174 integer :: i, id, dsize, nc
2175
2176 if (lvl < mg%lowest_lvl) &
2177 error stop "fill_ghost_cells_lvl: lvl < lowest_lvl"
2178 if (lvl > mg%highest_lvl) &
2179 error stop "fill_ghost_cells_lvl: lvl > highest_lvl"
2180
2181 nc = mg%box_size_lvl(lvl)
2182
2183 if (lvl >= mg%first_normal_lvl) then
2184 dsize = nc**(3-1)
2185 mg%buf(:)%i_send = 0
2186 mg%buf(:)%i_recv = 0
2187 mg%buf(:)%i_ix = 0
2188
2189 do i = 1, size(mg%lvls(lvl)%my_ids)
2190 id = mg%lvls(lvl)%my_ids(i)
2191 call buffer_ghost_cells(mg, id, nc, iv, .false.)
2192 end do
2193
2194 if (lvl > 1) then
2195 do i = 1, size(mg%lvls(lvl-1)%my_ref_bnds)
2196 id = mg%lvls(lvl-1)%my_ref_bnds(i)
2197 call buffer_refinement_boundaries(mg, id, nc, iv, .false.)
2198 end do
2199 end if
2200
2201 ! Transfer data between processes
2202 mg%buf(:)%i_recv = mg%comm_ghostcell%n_recv(:, lvl) * dsize
2203 call sort_and_transfer_buffers(mg, dsize)
2204
2205 ! Set ghost cells to received data
2206 mg%buf(:)%i_recv = 0
2207 end if
2208
2209 do i = 1, size(mg%lvls(lvl)%my_ids)
2210 id = mg%lvls(lvl)%my_ids(i)
2211 call set_ghost_cells(mg, id, nc, iv, .false.)
2212 end do
2213 end subroutine mg_fill_ghost_cells_lvl
2214
2215 subroutine buffer_ghost_cells(mg, id, nc, iv, dry_run)
2216 type(mg_t), intent(inout) :: mg
2217 integer, intent(in) :: id
2218 integer, intent(in) :: nc
2219 integer, intent(in) :: iv
2220 logical, intent(in) :: dry_run
2221 integer :: nb, nb_id, nb_rank
2222
2223 do nb = 1, mg_num_neighbors
2224 nb_id = mg%boxes(id)%neighbors(nb)
2225
2226 if (nb_id > mg_no_box) then
2227 ! There is a neighbor
2228 nb_rank = mg%boxes(nb_id)%rank
2229
2230 if (nb_rank /= mg%my_rank) then
2231 call buffer_for_nb(mg, mg%boxes(id), nc, iv, nb_id, nb_rank, &
2232 nb, dry_run)
2233 end if
2234 end if
2235 end do
2236 end subroutine buffer_ghost_cells
2237
2238 subroutine buffer_refinement_boundaries(mg, id, nc, iv, dry_run)
2239 type(mg_t), intent(inout) :: mg
2240 integer, intent(in) :: id
2241 integer, intent(in) :: nc
2242 integer, intent(in) :: iv
2243 logical, intent(in) :: dry_run
2244 integer :: nb, nb_id, c_ids(2**(3-1))
2245 integer :: n, c_id, c_rank
2246
2247 do nb = 1, mg_num_neighbors
2248 nb_id = mg%boxes(id)%neighbors(nb)
2249 if (nb_id > mg_no_box) then
2250 if (mg_has_children(mg%boxes(nb_id))) then
2251 c_ids = mg%boxes(nb_id)%children(&
2253
2254 do n = 1, mg_num_children/2
2255 c_id = c_ids(n)
2256 c_rank = mg%boxes(c_id)%rank
2257
2258 if (c_rank /= mg%my_rank) then
2259 ! Send all coarse ghost cells
2260 call buffer_for_fine_nb(mg, mg%boxes(id), nc, iv, c_id, &
2261 c_rank, nb, dry_run)
2262 end if
2263 end do
2264 end if
2265 end if
2266 end do
2267 end subroutine buffer_refinement_boundaries
2268
2269 !> The routine that actually fills the ghost cells
2270 subroutine set_ghost_cells(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 real(dp) :: bc(nc, nc)
2277 integer :: nb, nb_id, nb_rank, bc_type
2278
2279 do nb = 1, mg_num_neighbors
2280 nb_id = mg%boxes(id)%neighbors(nb)
2281
2282 if (nb_id > mg_no_box) then
2283 ! There is a neighbor
2284 nb_rank = mg%boxes(nb_id)%rank
2285
2286 if (nb_rank /= mg%my_rank) then
2287 call fill_buffered_nb(mg, mg%boxes(id), nb_rank, &
2288 nb, nc, iv, dry_run)
2289 else if (.not. dry_run) then
2290 call copy_from_nb(mg%boxes(id), mg%boxes(nb_id), &
2291 nb, nc, iv)
2292 end if
2293 else if (nb_id == mg_no_box) then
2294 ! Refinement boundary
2295 call fill_refinement_bnd(mg, id, nb, nc, iv, dry_run)
2296 else if (.not. dry_run) then
2297 ! Physical boundary
2298 if (mg%phi_bc_data_stored .and. iv == mg_iphi) then
2299 ! Copy the boundary conditions stored in the ghost cells of the
2300 ! right-hand side
2301 call box_get_gc(mg%boxes(id), nb, nc, mg_irhs, bc)
2302 bc_type = nb_id
2303 else
2304 if (associated(mg%bc(nb, iv)%boundary_cond)) then
2305 call mg%bc(nb, iv)%boundary_cond(mg%boxes(id), nc, iv, &
2306 nb, bc_type, bc)
2307 else
2308 bc_type = mg%bc(nb, iv)%bc_type
2309 bc = mg%bc(nb, iv)%bc_value
2310 end if
2311 end if
2312
2313 call box_set_gc(mg%boxes(id), nb, nc, iv, bc)
2314 call bc_to_gc(mg, id, nc, iv, nb, bc_type)
2315 end if
2316 end do
2317 end subroutine set_ghost_cells
2318
2319 subroutine fill_refinement_bnd(mg, id, nb, nc, iv, dry_run)
2320 type(mg_t), intent(inout) :: mg
2321 integer, intent(in) :: id
2322 integer, intent(in) :: nc
2323 integer, intent(in) :: iv
2324 integer, intent(in) :: nb
2325 logical, intent(in) :: dry_run
2326 real(dp) :: gc(nc, nc)
2327 integer :: p_id, p_nb_id, ix_offset(3)
2328 integer :: i, dsize, p_nb_rank
2329
2330 dsize = nc**(3-1)
2331 p_id = mg%boxes(id)%parent
2332 p_nb_id = mg%boxes(p_id)%neighbors(nb)
2333 p_nb_rank = mg%boxes(p_nb_id)%rank
2334
2335 if (p_nb_rank /= mg%my_rank) then
2336 i = mg%buf(p_nb_rank)%i_recv
2337 if (.not. dry_run) then
2338 gc = reshape(mg%buf(p_nb_rank)%recv(i+1:i+dsize), shape(gc))
2339 end if
2340 mg%buf(p_nb_rank)%i_recv = mg%buf(p_nb_rank)%i_recv + dsize
2341 else if (.not. dry_run) then
2342 ix_offset = mg_get_child_offset(mg, id)
2343 call box_gc_for_fine_neighbor(mg%boxes(p_nb_id), mg_neighb_rev(nb), &
2344 ix_offset, nc, iv, gc)
2345 end if
2346
2347 if (.not. dry_run) then
2348 if (associated(mg%bc(nb, iv)%refinement_bnd)) then
2349 call mg%bc(nb, iv)%refinement_bnd(mg%boxes(id), nc, iv, nb, gc)
2350 else
2351 call sides_rb(mg%boxes(id), nc, iv, nb, gc)
2352 end if
2353 end if
2354 end subroutine fill_refinement_bnd
2355
2356 subroutine copy_from_nb(box, box_nb, nb, nc, iv)
2357 type(mg_box_t), intent(inout) :: box
2358 type(mg_box_t), intent(in) :: box_nb
2359 integer, intent(in) :: nb
2360 integer, intent(in) :: nc
2361 integer, intent(in) :: iv
2362 real(dp) :: gc(nc, nc)
2363
2364 call box_gc_for_neighbor(box_nb, mg_neighb_rev(nb), nc, iv, gc)
2365 call box_set_gc(box, nb, nc, iv, gc)
2366 end subroutine copy_from_nb
2367
2368 subroutine buffer_for_nb(mg, box, nc, iv, nb_id, nb_rank, nb, dry_run)
2369 type(mg_t), intent(inout) :: mg
2370 type(mg_box_t), intent(inout) :: box
2371 integer, intent(in) :: nc
2372 integer, intent(in) :: iv
2373 integer, intent(in) :: nb_id
2374 integer, intent(in) :: nb_rank
2375 integer, intent(in) :: nb
2376 logical, intent(in) :: dry_run
2377 integer :: i, dsize
2378 real(dp) :: gc(nc, nc)
2379
2380 i = mg%buf(nb_rank)%i_send
2381 dsize = nc**(3-1)
2382
2383 if (.not. dry_run) then
2384 call box_gc_for_neighbor(box, nb, nc, iv, gc)
2385 mg%buf(nb_rank)%send(i+1:i+dsize) = pack(gc, .true.)
2386 end if
2387
2388 ! Later the buffer is sorted, using the fact that loops go from low to high
2389 ! box id, and we fill ghost cells according to the neighbor order
2390 i = mg%buf(nb_rank)%i_ix
2391 if (.not. dry_run) then
2392 mg%buf(nb_rank)%ix(i+1) = mg_num_neighbors * nb_id + mg_neighb_rev(nb)
2393 end if
2394
2395 mg%buf(nb_rank)%i_send = mg%buf(nb_rank)%i_send + dsize
2396 mg%buf(nb_rank)%i_ix = mg%buf(nb_rank)%i_ix + 1
2397 end subroutine buffer_for_nb
2398
2399 subroutine buffer_for_fine_nb(mg, box, nc, iv, fine_id, fine_rank, nb, dry_run)
2400 type(mg_t), intent(inout) :: mg
2401 type(mg_box_t), intent(inout) :: box
2402 integer, intent(in) :: nc
2403 integer, intent(in) :: iv
2404 integer, intent(in) :: fine_id
2405 integer, intent(in) :: fine_rank
2406 integer, intent(in) :: nb
2407 logical, intent(in) :: dry_run
2408 integer :: i, dsize, ix_offset(3)
2409 real(dp) :: gc(nc, nc)
2410
2411 i = mg%buf(fine_rank)%i_send
2412 dsize = nc**(3-1)
2413
2414 if (.not. dry_run) then
2415 ix_offset = mg_get_child_offset(mg, fine_id)
2416 call box_gc_for_fine_neighbor(box, nb, ix_offset, nc, iv, gc)
2417 mg%buf(fine_rank)%send(i+1:i+dsize) = pack(gc, .true.)
2418 end if
2419
2420 ! Later the buffer is sorted, using the fact that loops go from low to high
2421 ! box id, and we fill ghost cells according to the neighbor order
2422 i = mg%buf(fine_rank)%i_ix
2423 if (.not. dry_run) then
2424 mg%buf(fine_rank)%ix(i+1) = mg_num_neighbors * fine_id + &
2425 mg_neighb_rev(nb)
2426 end if
2427
2428 mg%buf(fine_rank)%i_send = mg%buf(fine_rank)%i_send + dsize
2429 mg%buf(fine_rank)%i_ix = mg%buf(fine_rank)%i_ix + 1
2430 end subroutine buffer_for_fine_nb
2431
2432 subroutine fill_buffered_nb(mg, box, nb_rank, nb, nc, iv, dry_run)
2433 type(mg_t), intent(inout) :: mg
2434 type(mg_box_t), intent(inout) :: box
2435 integer, intent(in) :: nb_rank
2436 integer, intent(in) :: nb
2437 integer, intent(in) :: nc
2438 integer, intent(in) :: iv
2439 logical, intent(in) :: dry_run
2440 integer :: i, dsize
2441 real(dp) :: gc(nc, nc)
2442
2443 i = mg%buf(nb_rank)%i_recv
2444 dsize = nc**(3-1)
2445
2446 if (.not. dry_run) then
2447 gc = reshape(mg%buf(nb_rank)%recv(i+1:i+dsize), shape(gc))
2448 call box_set_gc(box, nb, nc, iv, gc)
2449 end if
2450 mg%buf(nb_rank)%i_recv = mg%buf(nb_rank)%i_recv + dsize
2451
2452 end subroutine fill_buffered_nb
2453
2454 subroutine box_gc_for_neighbor(box, nb, nc, iv, gc)
2455 type(mg_box_t), intent(in) :: box
2456 integer, intent(in) :: nb, nc, iv
2457 real(dp), intent(out) :: gc(nc, nc)
2458
2459 select case (nb)
2460 case (mg_neighb_lowx)
2461 gc = box%cc(1, 1:nc, 1:nc, iv)
2462 case (mg_neighb_highx)
2463 gc = box%cc(nc, 1:nc, 1:nc, iv)
2464 case (mg_neighb_lowy)
2465 gc = box%cc(1:nc, 1, 1:nc, iv)
2466 case (mg_neighb_highy)
2467 gc = box%cc(1:nc, nc, 1:nc, iv)
2468 case (mg_neighb_lowz)
2469 gc = box%cc(1:nc, 1:nc, 1, iv)
2470 case (mg_neighb_highz)
2471 gc = box%cc(1:nc, 1:nc, nc, iv)
2472 end select
2473 end subroutine box_gc_for_neighbor
2474
2475 !> Get ghost cells for a fine neighbor
2476 subroutine box_gc_for_fine_neighbor(box, nb, di, nc, iv, gc)
2477 type(mg_box_t), intent(in) :: box
2478 integer, intent(in) :: nb !< Direction of fine neighbor
2479 integer, intent(in) :: di(3) !< Index offset of fine neighbor
2480 integer, intent(in) :: nc, iv
2481 real(dp), intent(out) :: gc(nc, nc)
2482 real(dp) :: tmp(0:nc/2+1, 0:nc/2+1), grad(3-1)
2483 integer :: i, j, hnc
2484
2485 hnc = nc/2
2486
2487 ! First fill a temporary array with data next to the fine grid
2488 select case (nb)
2489 case (mg_neighb_lowx)
2490 tmp = box%cc(1, di(2):di(2)+hnc+1, di(3):di(3)+hnc+1, iv)
2491 case (mg_neighb_highx)
2492 tmp = box%cc(nc, di(2):di(2)+hnc+1, di(3):di(3)+hnc+1, iv)
2493 case (mg_neighb_lowy)
2494 tmp = box%cc(di(1):di(1)+hnc+1, 1, di(3):di(3)+hnc+1, iv)
2495 case (mg_neighb_highy)
2496 tmp = box%cc(di(1):di(1)+hnc+1, nc, di(3):di(3)+hnc+1, iv)
2497 case (mg_neighb_lowz)
2498 tmp = box%cc(di(1):di(1)+hnc+1, di(2):di(2)+hnc+1, 1, iv)
2499 case (mg_neighb_highz)
2500 tmp = box%cc(di(1):di(1)+hnc+1, di(2):di(2)+hnc+1, nc, iv)
2501 case default
2502 error stop
2503 end select
2504
2505 ! Now interpolate the coarse grid data to obtain values 'straight' next to
2506 ! the fine grid points
2507 do j = 1, hnc
2508 do i = 1, hnc
2509 grad(1) = 0.125_dp * (tmp(i+1, j) - tmp(i-1, j))
2510 grad(2) = 0.125_dp * (tmp(i, j+1) - tmp(i, j-1))
2511 gc(2*i-1, 2*j-1) = tmp(i, j) - grad(1) - grad(2)
2512 gc(2*i, 2*j-1) = tmp(i, j) + grad(1) - grad(2)
2513 gc(2*i-1, 2*j) = tmp(i, j) - grad(1) + grad(2)
2514 gc(2*i, 2*j) = tmp(i, j) + grad(1) + grad(2)
2515 end do
2516 end do
2517 end subroutine box_gc_for_fine_neighbor
2518
2519 subroutine box_get_gc(box, nb, nc, iv, gc)
2520 type(mg_box_t), intent(in) :: box
2521 integer, intent(in) :: nb, nc, iv
2522 real(dp), intent(out) :: gc(nc, nc)
2523
2524 select case (nb)
2525 case (mg_neighb_lowx)
2526 gc = box%cc(0, 1:nc, 1:nc, iv)
2527 case (mg_neighb_highx)
2528 gc = box%cc(nc+1, 1:nc, 1:nc, iv)
2529 case (mg_neighb_lowy)
2530 gc = box%cc(1:nc, 0, 1:nc, iv)
2531 case (mg_neighb_highy)
2532 gc = box%cc(1:nc, nc+1, 1:nc, iv)
2533 case (mg_neighb_lowz)
2534 gc = box%cc(1:nc, 1:nc, 0, iv)
2535 case (mg_neighb_highz)
2536 gc = box%cc(1:nc, 1:nc, nc+1, iv)
2537 end select
2538 end subroutine box_get_gc
2539
2540 subroutine box_set_gc(box, nb, nc, iv, gc)
2541 type(mg_box_t), intent(inout) :: box
2542 integer, intent(in) :: nb, nc, iv
2543 real(dp), intent(in) :: gc(nc, nc)
2544
2545 select case (nb)
2546 case (mg_neighb_lowx)
2547 box%cc(0, 1:nc, 1:nc, iv) = gc
2548 case (mg_neighb_highx)
2549 box%cc(nc+1, 1:nc, 1:nc, iv) = gc
2550 case (mg_neighb_lowy)
2551 box%cc(1:nc, 0, 1:nc, iv) = gc
2552 case (mg_neighb_highy)
2553 box%cc(1:nc, nc+1, 1:nc, iv) = gc
2554 case (mg_neighb_lowz)
2555 box%cc(1:nc, 1:nc, 0, iv) = gc
2556 case (mg_neighb_highz)
2557 box%cc(1:nc, 1:nc, nc+1, iv) = gc
2558 end select
2559 end subroutine box_set_gc
2560
2561 subroutine bc_to_gc(mg, id, nc, iv, nb, bc_type)
2562 type(mg_t), intent(inout) :: mg
2563 integer, intent(in) :: id
2564 integer, intent(in) :: nc
2565 integer, intent(in) :: iv
2566 integer, intent(in) :: nb !< Neighbor direction
2567 integer, intent(in) :: bc_type !< Type of b.c.
2568 real(dp) :: c0, c1, c2, dr
2569
2570 ! If we call the interior point x1, x2 and the ghost point x0, then a
2571 ! Dirichlet boundary value b can be imposed as:
2572 ! x0 = -x1 + 2*b
2573 ! A Neumann b.c. can be imposed as:
2574 ! x0 = x1 +/- dx * b
2575 ! A continuous boundary (same slope) as:
2576 ! x0 = 2 * x1 - x2
2577 ! Below, we set coefficients to handle these cases
2578 select case (bc_type)
2579 case (mg_bc_dirichlet)
2580 c0 = 2
2581 c1 = -1
2582 c2 = 0
2583 case (mg_bc_neumann)
2584 dr = mg%dr(mg_neighb_dim(nb), mg%boxes(id)%lvl)
2585 c0 = dr * mg_neighb_high_pm(nb) ! This gives a + or - sign
2586 c1 = 1
2587 c2 = 0
2588 case (mg_bc_continuous)
2589 c0 = 0
2590 c1 = 2
2591 c2 = -1
2592 case default
2593 error stop "bc_to_gc: unknown boundary condition"
2594 end select
2595
2596 select case (nb)
2597 case (mg_neighb_lowx)
2598 mg%boxes(id)%cc(0, 1:nc, 1:nc, iv) = &
2599 c0 * mg%boxes(id)%cc(0, 1:nc, 1:nc, iv) + &
2600 c1 * mg%boxes(id)%cc(1, 1:nc, 1:nc, iv) + &
2601 c2 * mg%boxes(id)%cc(2, 1:nc, 1:nc, iv)
2602 case (mg_neighb_highx)
2603 mg%boxes(id)%cc(nc+1, 1:nc, 1:nc, iv) = &
2604 c0 * mg%boxes(id)%cc(nc+1, 1:nc, 1:nc, iv) + &
2605 c1 * mg%boxes(id)%cc(nc, 1:nc, 1:nc, iv) + &
2606 c2 * mg%boxes(id)%cc(nc-1, 1:nc, 1:nc, iv)
2607 case (mg_neighb_lowy)
2608 mg%boxes(id)%cc(1:nc, 0, 1:nc, iv) = &
2609 c0 * mg%boxes(id)%cc(1:nc, 0, 1:nc, iv) + &
2610 c1 * mg%boxes(id)%cc(1:nc, 1, 1:nc, iv) + &
2611 c2 * mg%boxes(id)%cc(1:nc, 2, 1:nc, iv)
2612 case (mg_neighb_highy)
2613 mg%boxes(id)%cc(1:nc, nc+1, 1:nc, iv) = &
2614 c0 * mg%boxes(id)%cc(1:nc, nc+1, 1:nc, iv) + &
2615 c1 * mg%boxes(id)%cc(1:nc, nc, 1:nc, iv) + &
2616 c2 * mg%boxes(id)%cc(1:nc, nc-1, 1:nc, iv)
2617 case (mg_neighb_lowz)
2618 mg%boxes(id)%cc(1:nc, 1:nc, 0, iv) = &
2619 c0 * mg%boxes(id)%cc(1:nc, 1:nc, 0, iv) + &
2620 c1 * mg%boxes(id)%cc(1:nc, 1:nc, 1, iv) + &
2621 c2 * mg%boxes(id)%cc(1:nc, 1:nc, 2, iv)
2622 case (mg_neighb_highz)
2623 mg%boxes(id)%cc(1:nc, 1:nc, nc+1, iv) = &
2624 c0 * mg%boxes(id)%cc(1:nc, 1:nc, nc+1, iv) + &
2625 c1 * mg%boxes(id)%cc(1:nc, 1:nc, nc, iv) + &
2626 c2 * mg%boxes(id)%cc(1:nc, 1:nc, nc-1, iv)
2627 end select
2628 end subroutine bc_to_gc
2629
2630 !> Fill ghost cells near refinement boundaries which preserves diffusive fluxes.
2631 subroutine sides_rb(box, nc, iv, nb, gc)
2632 type(mg_box_t), intent(inout) :: box
2633 integer, intent(in) :: nc
2634 integer, intent(in) :: iv
2635 integer, intent(in) :: nb !< Ghost cell direction
2636 !> Interpolated coarse grid ghost cell data (but not yet in the nb direction)
2637 real(dp), intent(in) :: gc(nc, nc)
2638 integer :: di, dj, dk
2639 integer :: i, j, k, ix, dix
2640
2641 if (mg_neighb_low(nb)) then
2642 ix = 1
2643 dix = 1
2644 else
2645 ix = nc
2646 dix = -1
2647 end if
2648
2649 select case (mg_neighb_dim(nb))
2650 case (1)
2651 i = ix
2652 di = dix
2653 do k = 1, nc
2654 dk = -1 + 2 * iand(k, 1)
2655 do j = 1, nc
2656 dj = -1 + 2 * iand(j, 1)
2657 box%cc(i-di, j, k, iv) = 0.5_dp * gc(j, k) &
2658 + 0.75_dp * box%cc(i, j, k, iv) &
2659 - 0.25_dp * box%cc(i+di, j, k, iv)
2660 end do
2661 end do
2662 case (2)
2663 j = ix
2664 dj = dix
2665 do k = 1, nc
2666 dk = -1 + 2 * iand(k, 1)
2667 do i = 1, nc
2668 di = -1 + 2 * iand(i, 1)
2669 box%cc(i, j-dj, k, iv) = 0.5_dp * gc(i, k) &
2670 + 0.75_dp * box%cc(i, j, k, iv) &
2671 - 0.25_dp * box%cc(i, j+dj, k, iv)
2672 end do
2673 end do
2674 case (3)
2675 k = ix
2676 dk = dix
2677 do j = 1, nc
2678 dj = -1 + 2 * iand(j, 1)
2679 do i = 1, nc
2680 di = -1 + 2 * iand(i, 1)
2681 box%cc(i, j, k-dk, iv) = 0.5_dp * gc(i, j) &
2682 + 0.75_dp * box%cc(i, j, k, iv) &
2683 - 0.25_dp * box%cc(i, j, k+dk, iv)
2684 end do
2685 end do
2686 end select
2687
2688 end subroutine sides_rb
2689
2690 !! File ../src/m_prolong.f90
2691
2692 !> Specify minimum buffer size (per process) for communication
2693 subroutine mg_prolong_buffer_size(mg, n_send, n_recv, dsize)
2694 type(mg_t), intent(inout) :: mg
2695 integer, intent(out) :: n_send(0:mg%n_cpu-1)
2696 integer, intent(out) :: n_recv(0:mg%n_cpu-1)
2697 integer, intent(out) :: dsize
2698 integer :: lvl, min_lvl
2699
2700 if (.not. allocated(mg%comm_restrict%n_send)) then
2701 error stop "Call restrict_buffer_size before prolong_buffer_size"
2702 end if
2703
2704 min_lvl = max(mg%first_normal_lvl-1, mg%lowest_lvl)
2705 allocate(mg%comm_prolong%n_send(0:mg%n_cpu-1, &
2706 min_lvl:mg%highest_lvl))
2707 allocate(mg%comm_prolong%n_recv(0:mg%n_cpu-1, &
2708 min_lvl:mg%highest_lvl))
2709
2710 mg%comm_prolong%n_recv(:, mg%highest_lvl) = 0
2711 mg%comm_prolong%n_send(:, mg%highest_lvl) = 0
2712
2713 do lvl = min_lvl, mg%highest_lvl-1
2714 mg%comm_prolong%n_recv(:, lvl) = &
2715 mg%comm_restrict%n_send(:, lvl+1)
2716 mg%comm_prolong%n_send(:, lvl) = &
2717 mg%comm_restrict%n_recv(:, lvl+1)
2718 end do
2719
2720 ! Send fine grid points, because this is more flexible than sending coarse
2721 ! grid points (e.g., when multiple variables are used for interpolation)
2722 dsize = (mg%box_size)**3
2723 n_send = maxval(mg%comm_prolong%n_send, dim=2)
2724 n_recv = maxval(mg%comm_prolong%n_recv, dim=2)
2725 end subroutine mg_prolong_buffer_size
2726
2727 !> Prolong variable iv from lvl to variable iv_to at lvl+1
2728 subroutine mg_prolong(mg, lvl, iv, iv_to, method, add)
2729
2730 type(mg_t), intent(inout) :: mg
2731 integer, intent(in) :: lvl !< Level to prolong from
2732 integer, intent(in) :: iv !< Source variable
2733 integer, intent(in) :: iv_to !< Target variable
2734 procedure(mg_box_prolong) :: method !< Prolongation method
2735 logical, intent(in) :: add !< If true, add to current values
2736 integer :: i, id, dsize, nc
2737
2738 if (lvl == mg%highest_lvl) error stop "cannot prolong highest level"
2739 if (lvl < mg%lowest_lvl) error stop "cannot prolong below lowest level"
2740
2741 ! Below the first normal level, all boxes are on the same CPU
2742 if (lvl >= mg%first_normal_lvl-1) then
2743 dsize = mg%box_size**3
2744 mg%buf(:)%i_send = 0
2745 mg%buf(:)%i_ix = 0
2746
2747 do i = 1, size(mg%lvls(lvl)%my_ids)
2748 id = mg%lvls(lvl)%my_ids(i)
2749 call prolong_set_buffer(mg, id, mg%box_size, iv, method)
2750 end do
2751
2752 mg%buf(:)%i_recv = mg%comm_prolong%n_recv(:, lvl) * dsize
2753 call sort_and_transfer_buffers(mg, dsize)
2754 mg%buf(:)%i_recv = 0
2755 end if
2756
2757 nc = mg%box_size_lvl(lvl+1)
2758 do i = 1, size(mg%lvls(lvl+1)%my_ids)
2759 id = mg%lvls(lvl+1)%my_ids(i)
2760 call prolong_onto(mg, id, nc, iv, iv_to, add, method)
2761 end do
2762 end subroutine mg_prolong
2763
2764 !> In case the fine grid is on a different CPU, perform the prolongation and
2765 !> store the fine-grid values in the send buffer.
2766 !>
2767 subroutine prolong_set_buffer(mg, id, nc, iv, method)
2768 type(mg_t), intent(inout) :: mg
2769 integer, intent(in) :: id
2770 integer, intent(in) :: nc
2771 integer, intent(in) :: iv
2772 procedure(mg_box_prolong) :: method
2773 integer :: i, dix(3)
2774 integer :: i_c, c_id, c_rank, dsize
2775 real(dp) :: tmp(nc, nc, nc)
2776
2777 dsize = nc**3
2778
2779 do i_c = 1, mg_num_children
2780 c_id = mg%boxes(id)%children(i_c)
2781 if (c_id > mg_no_box) then
2782 c_rank = mg%boxes(c_id)%rank
2783 if (c_rank /= mg%my_rank) then
2784 dix = mg_get_child_offset(mg, c_id)
2785 call method(mg, id, dix, nc, iv, tmp)
2786
2787 i = mg%buf(c_rank)%i_send
2788 mg%buf(c_rank)%send(i+1:i+dsize) = pack(tmp, .true.)
2789 mg%buf(c_rank)%i_send = mg%buf(c_rank)%i_send + dsize
2790
2791 i = mg%buf(c_rank)%i_ix
2792 mg%buf(c_rank)%ix(i+1) = c_id
2793 mg%buf(c_rank)%i_ix = mg%buf(c_rank)%i_ix + 1
2794 end if
2795 end if
2796 end do
2797 end subroutine prolong_set_buffer
2798
2799 !> Prolong onto a child box
2800 subroutine prolong_onto(mg, id, nc, iv, iv_to, add, method)
2801 type(mg_t), intent(inout) :: mg
2802 integer, intent(in) :: id
2803 integer, intent(in) :: nc
2804 integer, intent(in) :: iv !< Prolong from this variable
2805 integer, intent(in) :: iv_to !< Prolong to this variable
2806 logical, intent(in) :: add !< If true, add to current values
2807 procedure(mg_box_prolong) :: method
2808 integer :: hnc, p_id, p_rank, i, dix(3), dsize
2809 real(dp) :: tmp(nc, nc, nc)
2810
2811 hnc = nc/2
2812 p_id = mg%boxes(id)%parent
2813 p_rank = mg%boxes(p_id)%rank
2814
2815 if (p_rank == mg%my_rank) then
2816 dix = mg_get_child_offset(mg, id)
2817 call method(mg, p_id, dix, nc, iv, tmp)
2818 else
2819 dsize = nc**3
2820 i = mg%buf(p_rank)%i_recv
2821 tmp = reshape(mg%buf(p_rank)%recv(i+1:i+dsize), [nc, nc, nc])
2822 mg%buf(p_rank)%i_recv = mg%buf(p_rank)%i_recv + dsize
2823 end if
2824
2825 if (add) then
2826 mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, iv_to) = &
2827 mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, iv_to) + tmp
2828 else
2829 mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, iv_to) = tmp
2830 end if
2831
2832 end subroutine prolong_onto
2833
2834 !> Prolong from a parent to a child with index offset dix
2835 subroutine mg_prolong_sparse(mg, p_id, dix, nc, iv, fine)
2836 type(mg_t), intent(inout) :: mg
2837 integer, intent(in) :: p_id !< Id of parent
2838 integer, intent(in) :: dix(3) !< Offset of child in parent grid
2839 integer, intent(in) :: nc !< Child grid size
2840 integer, intent(in) :: iv !< Prolong from this variable
2841 real(dp), intent(out) :: fine(nc, nc, nc) !< Prolonged values
2842
2843 integer :: i, j, k, hnc
2844 integer :: ic, jc, kc
2845 real(dp) :: f0, flx, fhx, fly, fhy, flz, fhz
2846
2847 hnc = nc/2
2848
2849 associate(crs => mg%boxes(p_id)%cc)
2850 do k = 1, hnc
2851 kc = k + dix(3)
2852 do j = 1, hnc
2853 jc = j + dix(2)
2854 do i = 1, hnc
2855 ic = i + dix(1)
2856
2857 f0 = 0.25_dp * crs(ic, jc, kc, iv)
2858 flx = 0.25_dp * crs(ic-1, jc, kc, iv)
2859 fhx = 0.25_dp * crs(ic+1, jc, kc, iv)
2860 fly = 0.25_dp * crs(ic, jc-1, kc, iv)
2861 fhy = 0.25_dp * crs(ic, jc+1, kc, iv)
2862 flz = 0.25_dp * crs(ic, jc, kc-1, iv)
2863 fhz = 0.25_dp * crs(ic, jc, kc+1, iv)
2864
2865 fine(2*i-1, 2*j-1, 2*k-1) = f0 + flx + fly + flz
2866 fine(2*i, 2*j-1, 2*k-1) = f0 + fhx + fly + flz
2867 fine(2*i-1, 2*j, 2*k-1) = f0 + flx + fhy + flz
2868 fine(2*i, 2*j, 2*k-1) = f0 + fhx + fhy + flz
2869 fine(2*i-1, 2*j-1, 2*k) = f0 + flx + fly + fhz
2870 fine(2*i, 2*j-1, 2*k) = f0 + fhx + fly + fhz
2871 fine(2*i-1, 2*j, 2*k) = f0 + flx + fhy + fhz
2872 fine(2*i, 2*j, 2*k) = f0 + fhx + fhy + fhz
2873 end do
2874 end do
2875 end do
2876 end associate
2877 end subroutine mg_prolong_sparse
2878
2879 !! File ../src/m_helmholtz.f90
2880
2882 type(mg_t), intent(inout) :: mg
2883
2884 mg%subtract_mean = .false.
2885
2886 select case (mg%geometry_type)
2887 case (mg_cartesian)
2888 mg%box_op => box_helmh
2889
2890 select case (mg%smoother_type)
2892 mg%box_smoother => box_gs_helmh
2893 case default
2894 error stop "helmholtz_set_methods: unsupported smoother type"
2895 end select
2896 case default
2897 error stop "helmholtz_set_methods: unsupported geometry"
2898 end select
2899
2900 end subroutine helmholtz_set_methods
2901
2902 subroutine helmholtz_set_lambda(lambda)
2903 real(dp), intent(in) :: lambda
2904
2905 if (lambda < 0) &
2906 error stop "helmholtz_set_lambda: lambda < 0 not allowed"
2907
2908 helmholtz_lambda = lambda
2909 end subroutine helmholtz_set_lambda
2910
2911 !> Perform Gauss-Seidel relaxation on box for a Helmholtz operator
2912 subroutine box_gs_helmh(mg, id, nc, redblack_cntr)
2913 type(mg_t), intent(inout) :: mg
2914 integer, intent(in) :: id
2915 integer, intent(in) :: nc
2916 integer, intent(in) :: redblack_cntr !< Iteration counter
2917 integer :: i, j, k, i0, di
2918 real(dp) :: idr2(3), fac
2919 logical :: redblack
2920 real(dp), parameter :: sixth = 1/6.0_dp
2921
2922 idr2 = 1/mg%dr(:, mg%boxes(id)%lvl)**2
2923 fac = 1.0_dp / (2 * sum(idr2) + helmholtz_lambda)
2924 i0 = 1
2925
2926 redblack = (mg%smoother_type == mg_smoother_gsrb)
2927 if (redblack) then
2928 di = 2
2929 else
2930 di = 1
2931 end if
2932
2933 ! The parity of redblack_cntr determines which cells we use. If
2934 ! redblack_cntr is even, we use the even cells and vice versa.
2935 associate(cc => mg%boxes(id)%cc, n => mg_iphi)
2936 do k = 1, nc
2937 do j = 1, nc
2938 if (redblack) &
2939 i0 = 2 - iand(ieor(redblack_cntr, k+j), 1)
2940 do i = i0, nc, di
2941 cc(i, j, k, n) = fac * ( &
2942 idr2(1) * (cc(i+1, j, k, n) + cc(i-1, j, k, n)) + &
2943 idr2(2) * (cc(i, j+1, k, n) + cc(i, j-1, k, n)) + &
2944 idr2(3) * (cc(i, j, k+1, n) + cc(i, j, k-1, n)) - &
2945 cc(i, j, k, mg_irhs))
2946 end do
2947 end do
2948 end do
2949 end associate
2950 end subroutine box_gs_helmh
2951
2952 !> Perform Helmholtz operator on a box
2953 subroutine box_helmh(mg, id, nc, i_out)
2954 type(mg_t), intent(inout) :: mg
2955 integer, intent(in) :: id
2956 integer, intent(in) :: nc
2957 integer, intent(in) :: i_out !< Index of variable to store Helmholtz in
2958 integer :: i, j, k
2959 real(dp) :: idr2(3)
2960
2961 idr2 = 1 / mg%dr(:, mg%boxes(id)%lvl)**2
2962
2963 associate(cc => mg%boxes(id)%cc, n => mg_iphi)
2964 do k = 1, nc
2965 do j = 1, nc
2966 do i = 1, nc
2967 cc(i, j, k, i_out) = &
2968 idr2(1) * (cc(i-1, j, k, n) + cc(i+1, j, k, n) &
2969 - 2 * cc(i, j, k, n)) &
2970 + idr2(2) * (cc(i, j-1, k, n) + cc(i, j+1, k, n) &
2971 - 2 * cc(i, j, k, n)) &
2972 + idr2(3) * (cc(i, j, k-1, n) + cc(i, j, k+1, n) &
2973 - 2 * cc(i, j, k, n)) - &
2974 helmholtz_lambda * cc(i, j, k, n)
2975 end do
2976 end do
2977 end do
2978 end associate
2979 end subroutine box_helmh
2980
2981 !! File ../src/m_multigrid.f90
2982
2983 subroutine mg_set_methods(mg)
2984
2985 type(mg_t), intent(inout) :: mg
2986
2987 ! Set default prolongation method (routines below can override this)
2988 mg%box_prolong => mg_prolong_sparse
2989
2990 select case (mg%operator_type)
2991 case (mg_laplacian)
2992 call laplacian_set_methods(mg)
2993 case (mg_vlaplacian)
2994 call vlaplacian_set_methods(mg)
2995 case (mg_helmholtz)
2996 call helmholtz_set_methods(mg)
2997 case (mg_vhelmholtz)
2998 call vhelmholtz_set_methods(mg)
2999 case (mg_ahelmholtz)
3000 call ahelmholtz_set_methods(mg)
3001 case default
3002 error stop "mg_set_methods: unknown operator"
3003 end select
3004
3005 ! For red-black, perform two smoothing sub-steps so that all unknowns are
3006 ! updated per cycle
3007 if (mg%smoother_type == mg_smoother_gsrb) then
3008 mg%n_smoother_substeps = 2
3009 else
3010 mg%n_smoother_substeps = 1
3011 end if
3012 end subroutine mg_set_methods
3013
3014 subroutine check_methods(mg)
3015 type(mg_t), intent(inout) :: mg
3016
3017 if (.not. associated(mg%box_op) .or. &
3018 .not. associated(mg%box_smoother)) then
3019 call mg_set_methods(mg)
3020 end if
3021
3022 end subroutine check_methods
3023
3024 subroutine mg_add_timers(mg)
3025 type(mg_t), intent(inout) :: mg
3026 timer_total_vcycle = mg_add_timer(mg, "mg total V-cycle")
3027 timer_total_fmg = mg_add_timer(mg, "mg total FMG cycle")
3028 timer_smoother = mg_add_timer(mg, "mg smoother")
3029 timer_smoother_gc = mg_add_timer(mg, "mg smoother g.c.")
3030 timer_coarse = mg_add_timer(mg, "mg coarse")
3031 timer_correct = mg_add_timer(mg, "mg correct")
3032 timer_update_coarse = mg_add_timer(mg, "mg update coarse")
3033 end subroutine mg_add_timers
3034
3035 !> Perform FAS-FMG cycle (full approximation scheme, full multigrid).
3036 subroutine mg_fas_fmg(mg, have_guess, max_res)
3037 type(mg_t), intent(inout) :: mg
3038 logical, intent(in) :: have_guess !< If false, start from phi = 0
3039 real(dp), intent(out), optional :: max_res !< Store max(abs(residual))
3040 integer :: lvl, i, id
3041
3042 call check_methods(mg)
3043 if (timer_smoother == -1) call mg_add_timers(mg)
3044
3045 call mg_timer_start(mg%timers(timer_total_fmg))
3046
3047 if (.not. have_guess) then
3048 do lvl = mg%highest_lvl, mg%lowest_lvl, -1
3049 do i = 1, size(mg%lvls(lvl)%my_ids)
3050 id = mg%lvls(lvl)%my_ids(i)
3051 mg%boxes(id)%cc(:, :, :, mg_iphi) = 0.0_dp
3052 end do
3053 end do
3054 end if
3055
3056 ! Ensure ghost cells are filled correctly
3057 call mg_fill_ghost_cells_lvl(mg, mg%highest_lvl, mg_iphi)
3058
3059 do lvl = mg%highest_lvl, mg%lowest_lvl+1, -1
3060 ! Set rhs on coarse grid and restrict phi
3061 call mg_timer_start(mg%timers(timer_update_coarse))
3062 call update_coarse(mg, lvl)
3063 call mg_timer_end(mg%timers(timer_update_coarse))
3064 end do
3065
3066 if (mg%subtract_mean) then
3067 ! For fully periodic solutions, the mean source term has to be zero
3068 call subtract_mean(mg, mg_irhs, .false.)
3069 end if
3070
3071 do lvl = mg%lowest_lvl, mg%highest_lvl
3072 ! Store phi_old
3073 do i = 1, size(mg%lvls(lvl)%my_ids)
3074 id = mg%lvls(lvl)%my_ids(i)
3075 mg%boxes(id)%cc(:, :, :, mg_iold) = &
3076 mg%boxes(id)%cc(:, :, :, mg_iphi)
3077 end do
3078
3079 if (lvl > mg%lowest_lvl) then
3080 ! Correct solution at this lvl using lvl-1 data
3081 ! phi = phi + prolong(phi_coarse - phi_old_coarse)
3082 call mg_timer_start(mg%timers(timer_correct))
3083 call correct_children(mg, lvl-1)
3084 call mg_timer_end(mg%timers(timer_correct))
3085
3086 ! Update ghost cells
3087 call mg_fill_ghost_cells_lvl(mg, lvl, mg_iphi)
3088 end if
3089
3090 ! Perform V-cycle, possibly set residual on last iteration
3091 if (lvl == mg%highest_lvl) then
3092 call mg_fas_vcycle(mg, lvl, max_res, standalone=.false.)
3093 else
3094 call mg_fas_vcycle(mg, lvl, standalone=.false.)
3095 end if
3096 end do
3097
3098 call mg_timer_end(mg%timers(timer_total_fmg))
3099 end subroutine mg_fas_fmg
3100
3101 !> Perform FAS V-cycle (full approximation scheme).
3102 subroutine mg_fas_vcycle(mg, highest_lvl, max_res, standalone)
3103 use mpi
3104 type(mg_t), intent(inout) :: mg
3105 integer, intent(in), optional :: highest_lvl !< Maximum level for V-cycle
3106 real(dp), intent(out), optional :: max_res !< Store max(abs(residual))
3107 !> Whether the V-cycle is called by itself (default: true)
3108 logical, intent(in), optional :: standalone
3109 integer :: lvl, min_lvl, i, max_lvl, ierr
3110 real(dp) :: init_res, res
3111 logical :: is_standalone
3112
3113 is_standalone = .true.
3114 if (present(standalone)) is_standalone = standalone
3115
3116 call check_methods(mg)
3117 if (timer_smoother == -1) call mg_add_timers(mg)
3118
3119 if (is_standalone) &
3120 call mg_timer_start(mg%timers(timer_total_vcycle))
3121
3122 if (mg%subtract_mean .and. .not. present(highest_lvl)) then
3123 ! Assume that this is a stand-alone call. For fully periodic solutions,
3124 ! ensure the mean source term is zero.
3125 call subtract_mean(mg, mg_irhs, .false.)
3126 end if
3127
3128 min_lvl = mg%lowest_lvl
3129 max_lvl = mg%highest_lvl
3130 if (present(highest_lvl)) max_lvl = highest_lvl
3131
3132 ! Ensure ghost cells are filled correctly
3133 if (is_standalone) then
3134 call mg_fill_ghost_cells_lvl(mg, max_lvl, mg_iphi)
3135 end if
3136
3137 do lvl = max_lvl, min_lvl+1, -1
3138 ! Downwards relaxation
3139 call smooth_boxes(mg, lvl, mg%n_cycle_down)
3140
3141 ! Set rhs on coarse grid, restrict phi, and copy mg_iphi to mg_iold for the
3142 ! correction later
3143 call mg_timer_start(mg%timers(timer_update_coarse))
3144 call update_coarse(mg, lvl)
3145 call mg_timer_end(mg%timers(timer_update_coarse))
3146 end do
3147
3148 call mg_timer_start(mg%timers(timer_coarse))
3149 if (.not. all(mg%boxes(mg%lvls(min_lvl)%ids)%rank == &
3150 mg%boxes(mg%lvls(min_lvl)%ids(1))%rank)) then
3151 error stop "Multiple CPUs for coarse grid (not implemented yet)"
3152 end if
3153
3154 init_res = max_residual_lvl(mg, min_lvl)
3155 do i = 1, mg%max_coarse_cycles
3156 call smooth_boxes(mg, min_lvl, mg%n_cycle_up+mg%n_cycle_down)
3157 res = max_residual_lvl(mg, min_lvl)
3158 if (res < mg%residual_coarse_rel * init_res .or. &
3159 res < mg%residual_coarse_abs) exit
3160 end do
3161 call mg_timer_end(mg%timers(timer_coarse))
3162
3163 ! Do the upwards part of the v-cycle in the tree
3164 do lvl = min_lvl+1, max_lvl
3165 ! Correct solution at this lvl using lvl-1 data
3166 ! phi = phi + prolong(phi_coarse - phi_old_coarse)
3167 call mg_timer_start(mg%timers(timer_correct))
3168 call correct_children(mg, lvl-1)
3169
3170 ! Have to fill ghost cells after correction
3171 call mg_fill_ghost_cells_lvl(mg, lvl, mg_iphi)
3172 call mg_timer_end(mg%timers(timer_correct))
3173
3174 ! Upwards relaxation
3175 call smooth_boxes(mg, lvl, mg%n_cycle_up)
3176 end do
3177
3178 if (present(max_res)) then
3179 init_res = 0.0_dp
3180 do lvl = min_lvl, max_lvl
3181 res = max_residual_lvl(mg, lvl)
3182 init_res = max(res, init_res)
3183 end do
3184 call mpi_allreduce(init_res, max_res, 1, &
3185 mpi_double, mpi_max, mg%comm, ierr)
3186 end if
3187
3188 ! Subtract mean(phi) from phi
3189 if (mg%subtract_mean) then
3190 call subtract_mean(mg, mg_iphi, .true.)
3191 end if
3192
3193 if (is_standalone) &
3194 call mg_timer_end(mg%timers(timer_total_vcycle))
3195 end subroutine mg_fas_vcycle
3196
3197 subroutine subtract_mean(mg, iv, include_ghostcells)
3198 use mpi
3199 type(mg_t), intent(inout) :: mg
3200 integer, intent(in) :: iv
3201 logical, intent(in) :: include_ghostcells
3202 integer :: i, id, lvl, nc, ierr
3203 real(dp) :: sum_iv, mean_iv, volume
3204
3205 nc = mg%box_size
3206 sum_iv = get_sum(mg, iv)
3207 call mpi_allreduce(sum_iv, mean_iv, 1, &
3208 mpi_double, mpi_sum, mg%comm, ierr)
3209
3210 ! Divide by total grid volume to get mean
3211 volume = nc**3 * product(mg%dr(:, 1)) * size(mg%lvls(1)%ids)
3212 mean_iv = mean_iv / volume
3213
3214 do lvl = mg%lowest_lvl, mg%highest_lvl
3215 nc = mg%box_size_lvl(lvl)
3216
3217 do i = 1, size(mg%lvls(lvl)%my_ids)
3218 id = mg%lvls(lvl)%my_ids(i)
3219 if (include_ghostcells) then
3220 mg%boxes(id)%cc(:, :, :, iv) = &
3221 mg%boxes(id)%cc(:, :, :, iv) - mean_iv
3222 else
3223 mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, iv) = &
3224 mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, iv) - mean_iv
3225 end if
3226 end do
3227 end do
3228 end subroutine subtract_mean
3229
3230 real(dp) function get_sum(mg, iv)
3231 type(mg_t), intent(in) :: mg
3232 integer, intent(in) :: iv
3233 integer :: lvl, i, id, nc
3234 real(dp) :: w
3235
3236 get_sum = 0.0_dp
3237 do lvl = 1, mg%highest_lvl
3238 nc = mg%box_size_lvl(lvl)
3239 w = product(mg%dr(:, lvl)) ! Adjust for non-Cartesian cases
3240 do i = 1, size(mg%lvls(lvl)%my_leaves)
3241 id = mg%lvls(lvl)%my_leaves(i)
3242 get_sum = get_sum + w * &
3243 sum(mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, iv))
3244 end do
3245 end do
3246 end function get_sum
3247
3248 real(dp) function max_residual_lvl(mg, lvl)
3249 type(mg_t), intent(inout) :: mg
3250 integer, intent(in) :: lvl
3251 integer :: i, id, nc
3252 real(dp) :: res
3253
3254 nc = mg%box_size_lvl(lvl)
3255 max_residual_lvl = 0.0_dp
3256
3257 do i = 1, size(mg%lvls(lvl)%my_ids)
3258 id = mg%lvls(lvl)%my_ids(i)
3259 call residual_box(mg, id, nc)
3260 res = maxval(abs(mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, mg_ires)))
3261 max_residual_lvl = max(max_residual_lvl, res)
3262 end do
3263 end function max_residual_lvl
3264
3265 ! subroutine solve_coarse_grid(mg)
3266 ! use m_fishpack
3267 ! type(mg_t), intent(inout) :: mg
3268
3269 ! real(dp) :: rhs(mg%box_size, mg%box_size, mg%box_size)
3270 ! real(dp) :: rmin(3), rmax(3)
3271 ! integer :: nc, nx(3), my_boxes, total_boxes
3272
3273 ! my_boxes = size(mg%lvls(1)%my_ids)
3274 ! total_boxes = size(mg%lvls(1)%ids)
3275 ! nc = mg%box_size
3276
3277 ! if (my_boxes == total_boxes) then
3278 ! nx(:) = nc
3279 ! rmin = [0.0_dp, 0.0_dp, 0.0_dp]
3280 ! rmax = mg%dr(1) * [nc, nc, nc]
3281 ! rhs = mg%boxes(1)%cc(1:nc, 1:nc, 1:nc, mg_irhs)
3282
3283 ! #if 3 == 2
3284 ! call fishpack_2d(nx, rhs, mg%bc, rmin, rmax)
3285 ! #elif 3 == 3
3286 ! call fishpack_3d(nx, rhs, mg%bc, rmin, rmax)
3287 ! #endif
3288
3289 ! mg%boxes(1)%cc(1:nc, 1:nc, 1:nc, mg_iphi) = rhs
3290 ! else if (my_boxes > 0) then
3291 ! error stop "Boxes at level 1 at different processors"
3292 ! end if
3293
3294 ! call fill_ghost_cells_lvl(mg, 1)
3295 ! end subroutine solve_coarse_grid
3296
3297 ! Set rhs on coarse grid, restrict phi, and copy mg_iphi to mg_iold for the
3298 ! correction later
3299 subroutine update_coarse(mg, lvl)
3300 type(mg_t), intent(inout) :: mg !< Tree containing full grid
3301 integer, intent(in) :: lvl !< Update coarse values at lvl-1
3302 integer :: i, id, nc, nc_c
3303
3304 nc = mg%box_size_lvl(lvl)
3305 nc_c = mg%box_size_lvl(lvl-1)
3306
3307 ! Compute residual
3308 do i = 1, size(mg%lvls(lvl)%my_ids)
3309 id = mg%lvls(lvl)%my_ids(i)
3310 call residual_box(mg, id, nc)
3311 end do
3312
3313 ! Restrict phi and the residual
3314 call mg_restrict_lvl(mg, mg_iphi, lvl)
3315 call mg_restrict_lvl(mg, mg_ires, lvl)
3316
3317 call mg_fill_ghost_cells_lvl(mg, lvl-1, mg_iphi)
3318
3319 ! Set rhs_c = laplacian(phi_c) + restrict(res) where it is refined, and
3320 ! store current coarse phi in old.
3321 do i = 1, size(mg%lvls(lvl-1)%my_parents)
3322 id = mg%lvls(lvl-1)%my_parents(i)
3323
3324 ! Set rhs = L phi
3325 call mg%box_op(mg, id, nc_c, mg_irhs)
3326
3327 ! Add the fine grid residual to rhs
3328 mg%boxes(id)%cc(1:nc_c, 1:nc_c, 1:nc_c, mg_irhs) = &
3329 mg%boxes(id)%cc(1:nc_c, 1:nc_c, 1:nc_c, mg_irhs) + &
3330 mg%boxes(id)%cc(1:nc_c, 1:nc_c, 1:nc_c, mg_ires)
3331
3332 ! Story a copy of phi
3333 mg%boxes(id)%cc(:, :, :, mg_iold) = &
3334 mg%boxes(id)%cc(:, :, :, mg_iphi)
3335 end do
3336 end subroutine update_coarse
3337
3338 ! Sets phi = phi + prolong(phi_coarse - phi_old_coarse)
3339 subroutine correct_children(mg, lvl)
3340 type(mg_t), intent(inout) :: mg
3341 integer, intent(in) :: lvl
3342 integer :: i, id
3343
3344 do i = 1, size(mg%lvls(lvl)%my_parents)
3345 id = mg%lvls(lvl)%my_parents(i)
3346
3347 ! Store the correction in mg_ires
3348 mg%boxes(id)%cc(:, :, :, mg_ires) = &
3349 mg%boxes(id)%cc(:, :, :, mg_iphi) - &
3350 mg%boxes(id)%cc(:, :, :, mg_iold)
3351 end do
3352
3353 call mg_prolong(mg, lvl, mg_ires, mg_iphi, mg%box_prolong, add=.true.)
3354 end subroutine correct_children
3355
3356 subroutine smooth_boxes(mg, lvl, n_cycle)
3357 type(mg_t), intent(inout) :: mg
3358 integer, intent(in) :: lvl
3359 integer, intent(in) :: n_cycle !< Number of cycles to perform
3360 integer :: n, i, id, nc
3361
3362 nc = mg%box_size_lvl(lvl)
3363
3364 do n = 1, n_cycle * mg%n_smoother_substeps
3365 call mg_timer_start(mg%timers(timer_smoother))
3366 do i = 1, size(mg%lvls(lvl)%my_ids)
3367 id = mg%lvls(lvl)%my_ids(i)
3368 call mg%box_smoother(mg, id, nc, n)
3369 end do
3370 call mg_timer_end(mg%timers(timer_smoother))
3371
3372 call mg_timer_start(mg%timers(timer_smoother_gc))
3373 call mg_fill_ghost_cells_lvl(mg, lvl, mg_iphi)
3374 call mg_timer_end(mg%timers(timer_smoother_gc))
3375 end do
3376 end subroutine smooth_boxes
3377
3378 subroutine residual_box(mg, id, nc)
3379 type(mg_t), intent(inout) :: mg
3380 integer, intent(in) :: id
3381 integer, intent(in) :: nc
3382
3383 call mg%box_op(mg, id, nc, mg_ires)
3384
3385 mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, mg_ires) = &
3386 mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, mg_irhs) &
3387 - mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, mg_ires)
3388 end subroutine residual_box
3389
3390 !> Apply operator to the tree and store in variable i_out
3391 subroutine mg_apply_op(mg, i_out, op)
3392 type(mg_t), intent(inout) :: mg
3393 integer, intent(in) :: i_out
3394 procedure(mg_box_op), optional :: op
3395 integer :: lvl, i, id, nc
3396
3397 do lvl = mg%lowest_lvl, mg%highest_lvl
3398 nc = mg%box_size_lvl(lvl)
3399 do i = 1, size(mg%lvls(lvl)%my_ids)
3400 id = mg%lvls(lvl)%my_ids(i)
3401 if (present(op)) then
3402 call op(mg, id, nc, i_out)
3403 else
3404 call mg%box_op(mg, id, nc, i_out)
3405 end if
3406 end do
3407 end do
3408 end subroutine mg_apply_op
3409
3410 !! File ../src/m_restrict.f90
3411
3412 !> Specify minimum buffer size (per process) for communication
3413 subroutine mg_restrict_buffer_size(mg, n_send, n_recv, dsize)
3414 type(mg_t), intent(inout) :: mg
3415 integer, intent(out) :: n_send(0:mg%n_cpu-1)
3416 integer, intent(out) :: n_recv(0:mg%n_cpu-1)
3417 integer, intent(out) :: dsize
3418 integer :: n_out(0:mg%n_cpu-1, mg%first_normal_lvl:mg%highest_lvl)
3419 integer :: n_in(0:mg%n_cpu-1, mg%first_normal_lvl:mg%highest_lvl)
3420 integer :: lvl, i, id, p_id, p_rank
3421 integer :: i_c, c_id, c_rank, min_lvl
3422
3423 n_out(:, :) = 0
3424 n_in(:, :) = 0
3425 min_lvl = max(mg%lowest_lvl+1, mg%first_normal_lvl)
3426
3427 do lvl = min_lvl, mg%highest_lvl
3428 ! Number of messages to receive (at lvl-1)
3429 do i = 1, size(mg%lvls(lvl-1)%my_parents)
3430 id = mg%lvls(lvl-1)%my_parents(i)
3431 do i_c = 1, mg_num_children
3432 c_id = mg%boxes(id)%children(i_c)
3433
3434 if (c_id > mg_no_box) then
3435 c_rank = mg%boxes(c_id)%rank
3436 if (c_rank /= mg%my_rank) then
3437 n_in(c_rank, lvl) = n_in(c_rank, lvl) + 1
3438 end if
3439 end if
3440 end do
3441 end do
3442
3443 ! Number of messages to send
3444 do i = 1, size(mg%lvls(lvl)%my_ids)
3445 id = mg%lvls(lvl)%my_ids(i)
3446
3447 p_id = mg%boxes(id)%parent
3448 p_rank = mg%boxes(p_id)%rank
3449 if (p_rank /= mg%my_rank) then
3450 n_out(p_rank, lvl) = n_out(p_rank, lvl) + 1
3451 end if
3452
3453 end do
3454 end do
3455
3456 allocate(mg%comm_restrict%n_send(0:mg%n_cpu-1, &
3457 mg%first_normal_lvl:mg%highest_lvl))
3458 allocate(mg%comm_restrict%n_recv(0:mg%n_cpu-1, &
3459 mg%first_normal_lvl:mg%highest_lvl))
3460 mg%comm_restrict%n_send = n_out
3461 mg%comm_restrict%n_recv = n_in
3462
3463 dsize = (mg%box_size/2)**3
3464 n_send = maxval(n_out, dim=2)
3465 n_recv = maxval(n_in, dim=2)
3466 end subroutine mg_restrict_buffer_size
3467
3468 !> Restrict all levels
3469 subroutine mg_restrict(mg, iv)
3470 type(mg_t), intent(inout) :: mg
3471 integer, intent(in) :: iv
3472 integer :: lvl
3473
3474 do lvl = mg%highest_lvl, mg%lowest_lvl+1, -1
3475 call mg_restrict_lvl(mg, iv, lvl)
3476 end do
3477 end subroutine mg_restrict
3478
3479 !> Restrict from lvl to lvl-1
3480 subroutine mg_restrict_lvl(mg, iv, lvl)
3481
3482 type(mg_t), intent(inout) :: mg
3483 integer, intent(in) :: iv
3484 integer, intent(in) :: lvl
3485 integer :: i, id, dsize, nc
3486
3487 if (lvl <= mg%lowest_lvl) error stop "cannot restrict lvl <= lowest_lvl"
3488
3489 nc = mg%box_size_lvl(lvl)
3490
3491 if (lvl >= mg%first_normal_lvl) then
3492 dsize = (nc/2)**3
3493
3494 mg%buf(:)%i_send = 0
3495 mg%buf(:)%i_ix = 0
3496
3497 do i = 1, size(mg%lvls(lvl)%my_ids)
3498 id = mg%lvls(lvl)%my_ids(i)
3499 call restrict_set_buffer(mg, id, iv)
3500 end do
3501
3502 mg%buf(:)%i_recv = mg%comm_restrict%n_recv(:, lvl) * dsize
3503 call sort_and_transfer_buffers(mg, dsize)
3504 mg%buf(:)%i_recv = 0
3505 end if
3506
3507 do i = 1, size(mg%lvls(lvl-1)%my_parents)
3508 id = mg%lvls(lvl-1)%my_parents(i)
3509 call restrict_onto(mg, id, nc, iv)
3510 end do
3511 end subroutine mg_restrict_lvl
3512
3513 subroutine restrict_set_buffer(mg, id, iv)
3514 type(mg_t), intent(inout) :: mg
3515 integer, intent(in) :: id
3516 integer, intent(in) :: iv
3517 integer :: i, j, k, n, hnc, p_id, p_rank
3518 real(dp) :: tmp(mg%box_size/2, mg%box_size/2, mg%box_size/2)
3519
3520 hnc = mg%box_size/2
3521 p_id = mg%boxes(id)%parent
3522 p_rank = mg%boxes(p_id)%rank
3523
3524 if (p_rank /= mg%my_rank) then
3525 do k = 1, hnc
3526 do j = 1, hnc
3527 do i = 1, hnc
3528 tmp(i, j, k) = 0.125_dp * sum(mg%boxes(id)%cc(2*i-1:2*i, &
3529 2*j-1:2*j, 2*k-1:2*k, iv))
3530 end do
3531 end do
3532 end do
3533
3534 ! Buffer
3535 n = size(tmp)
3536 i = mg%buf(p_rank)%i_send
3537 mg%buf(p_rank)%send(i+1:i+n) = pack(tmp, .true.)
3538 mg%buf(p_rank)%i_send = mg%buf(p_rank)%i_send + n
3539
3540 ! To later sort the send buffer according to parent order
3541 i = mg%buf(p_rank)%i_ix
3542 n = mg_ix_to_ichild(mg%boxes(id)%ix)
3543 mg%buf(p_rank)%ix(i+1) = mg_num_children * p_id + n
3544 mg%buf(p_rank)%i_ix = mg%buf(p_rank)%i_ix + 1
3545 end if
3546 end subroutine restrict_set_buffer
3547
3548 subroutine restrict_onto(mg, id, nc, iv)
3549 type(mg_t), intent(inout) :: mg
3550 integer, intent(in) :: id
3551 integer, intent(in) :: nc
3552 integer, intent(in) :: iv
3553 integer :: i, j, k, hnc, dsize, i_c, c_id
3554 integer :: c_rank, dix(3)
3555
3556 hnc = nc/2
3557 dsize = hnc**3
3558
3559 do i_c = 1, mg_num_children
3560 c_id = mg%boxes(id)%children(i_c)
3561 if (c_id == mg_no_box) cycle ! For coarsened grid
3562 c_rank = mg%boxes(c_id)%rank
3563 dix = mg_get_child_offset(mg, c_id)
3564
3565 if (c_rank == mg%my_rank) then
3566 do j=1, hnc; do i=1, hnc; do k=1, hnc
3567 mg%boxes(id)%cc(dix(1)+i, dix(2)+j, dix(3)+k, iv) = &
3568 0.125_dp * sum(mg%boxes(c_id)%cc(2*i-1:2*i, &
3569 2*j-1:2*j, 2*k-1:2*k, iv))
3570 end do; end do; end do
3571 else
3572 i = mg%buf(c_rank)%i_recv
3573 mg%boxes(id)%cc(dix(1)+1:dix(1)+hnc, &
3574 dix(2)+1:dix(2)+hnc, dix(3)+1:dix(3)+hnc, iv) = &
3575 reshape(mg%buf(c_rank)%recv(i+1:i+dsize), [hnc, hnc, hnc])
3576 mg%buf(c_rank)%i_recv = mg%buf(c_rank)%i_recv + dsize
3577 end if
3578 end do
3579
3580 end subroutine restrict_onto
3581
3582 !! File ../src/m_mrgrnk.f90
3583
3584 Subroutine i_mrgrnk (XDONT, IRNGT)
3585 ! __________________________________________________________
3586 ! MRGRNK = Merge-sort ranking of an array
3587 ! For performance reasons, the first 2 passes are taken
3588 ! out of the standard loop, and use dedicated coding.
3589 ! __________________________________________________________
3590 ! __________________________________________________________
3591 Integer, Dimension (:), Intent (In) :: XDONT
3592 Integer, Dimension (:), Intent (Out) :: IRNGT
3593 ! __________________________________________________________
3594 Integer :: XVALA, XVALB
3595 !
3596 Integer, Dimension (SIZE(IRNGT)) :: JWRKT
3597 Integer :: LMTNA, LMTNC, IRNG1, IRNG2
3598 Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
3599 !
3600 nval = min(SIZE(xdont), SIZE(irngt))
3601 Select Case (nval)
3602 Case (:0)
3603 Return
3604 Case (1)
3605 irngt(1) = 1
3606 Return
3607 Case Default
3608 Continue
3609 End Select
3610 !
3611 ! Fill-in the index array, creating ordered couples
3612 !
3613 Do iind = 2, nval, 2
3614 If (xdont(iind-1) <= xdont(iind)) Then
3615 irngt(iind-1) = iind - 1
3616 irngt(iind) = iind
3617 Else
3618 irngt(iind-1) = iind
3619 irngt(iind) = iind - 1
3620 End If
3621 End Do
3622 If (modulo(nval, 2) /= 0) Then
3623 irngt(nval) = nval
3624 End If
3625 !
3626 ! We will now have ordered subsets A - B - A - B - ...
3627 ! and merge A and B couples into C - C - ...
3628 !
3629 lmtna = 2
3630 lmtnc = 4
3631 !
3632 ! First iteration. The length of the ordered subsets goes from 2 to 4
3633 !
3634 Do
3635 If (nval <= 2) Exit
3636 !
3637 ! Loop on merges of A and B into C
3638 !
3639 Do iwrkd = 0, nval - 1, 4
3640 If ((iwrkd+4) > nval) Then
3641 If ((iwrkd+2) >= nval) Exit
3642 !
3643 ! 1 2 3
3644 !
3645 If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) Exit
3646 !
3647 ! 1 3 2
3648 !
3649 If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3))) Then
3650 irng2 = irngt(iwrkd+2)
3651 irngt(iwrkd+2) = irngt(iwrkd+3)
3652 irngt(iwrkd+3) = irng2
3653 !
3654 ! 3 1 2
3655 !
3656 Else
3657 irng1 = irngt(iwrkd+1)
3658 irngt(iwrkd+1) = irngt(iwrkd+3)
3659 irngt(iwrkd+3) = irngt(iwrkd+2)
3660 irngt(iwrkd+2) = irng1
3661 End If
3662 Exit
3663 End If
3664 !
3665 ! 1 2 3 4
3666 !
3667 If (xdont(irngt(iwrkd+2)) <= xdont(irngt(iwrkd+3))) cycle
3668 !
3669 ! 1 3 x x
3670 !
3671 If (xdont(irngt(iwrkd+1)) <= xdont(irngt(iwrkd+3))) Then
3672 irng2 = irngt(iwrkd+2)
3673 irngt(iwrkd+2) = irngt(iwrkd+3)
3674 If (xdont(irng2) <= xdont(irngt(iwrkd+4))) Then
3675 ! 1 3 2 4
3676 irngt(iwrkd+3) = irng2
3677 Else
3678 ! 1 3 4 2
3679 irngt(iwrkd+3) = irngt(iwrkd+4)
3680 irngt(iwrkd+4) = irng2
3681 End If
3682 !
3683 ! 3 x x x
3684 !
3685 Else
3686 irng1 = irngt(iwrkd+1)
3687 irng2 = irngt(iwrkd+2)
3688 irngt(iwrkd+1) = irngt(iwrkd+3)
3689 If (xdont(irng1) <= xdont(irngt(iwrkd+4))) Then
3690 irngt(iwrkd+2) = irng1
3691 If (xdont(irng2) <= xdont(irngt(iwrkd+4))) Then
3692 ! 3 1 2 4
3693 irngt(iwrkd+3) = irng2
3694 Else
3695 ! 3 1 4 2
3696 irngt(iwrkd+3) = irngt(iwrkd+4)
3697 irngt(iwrkd+4) = irng2
3698 End If
3699 Else
3700 ! 3 4 1 2
3701 irngt(iwrkd+2) = irngt(iwrkd+4)
3702 irngt(iwrkd+3) = irng1
3703 irngt(iwrkd+4) = irng2
3704 End If
3705 End If
3706 End Do
3707 !
3708 ! The Cs become As and Bs
3709 !
3710 lmtna = 4
3711 Exit
3712 End Do
3713 !
3714 ! Iteration loop. Each time, the length of the ordered subsets
3715 ! is doubled.
3716 !
3717 Do
3718 If (lmtna >= nval) Exit
3719 iwrkf = 0
3720 lmtnc = 2 * lmtnc
3721 !
3722 ! Loop on merges of A and B into C
3723 !
3724 Do
3725 iwrk = iwrkf
3726 iwrkd = iwrkf + 1
3727 jinda = iwrkf + lmtna
3728 iwrkf = iwrkf + lmtnc
3729 If (iwrkf >= nval) Then
3730 If (jinda >= nval) Exit
3731 iwrkf = nval
3732 End If
3733 iinda = 1
3734 iindb = jinda + 1
3735 !
3736 ! Shortcut for the case when the max of A is smaller
3737 ! than the min of B. This line may be activated when the
3738 ! initial set is already close to sorted.
3739 !
3740 ! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
3741 !
3742 ! One steps in the C subset, that we build in the final rank array
3743 !
3744 ! Make a copy of the rank array for the merge iteration
3745 !
3746 jwrkt(1:lmtna) = irngt(iwrkd:jinda)
3747 !
3748 xvala = xdont(jwrkt(iinda))
3749 xvalb = xdont(irngt(iindb))
3750 !
3751 Do
3752 iwrk = iwrk + 1
3753 !
3754 ! We still have unprocessed values in both A and B
3755 !
3756 If (xvala > xvalb) Then
3757 irngt(iwrk) = irngt(iindb)
3758 iindb = iindb + 1
3759 If (iindb > iwrkf) Then
3760 ! Only A still with unprocessed values
3761 irngt(iwrk+1:iwrkf) = jwrkt(iinda:lmtna)
3762 Exit
3763 End If
3764 xvalb = xdont(irngt(iindb))
3765 Else
3766 irngt(iwrk) = jwrkt(iinda)
3767 iinda = iinda + 1
3768 If (iinda > lmtna) exit! Only B still with unprocessed values
3769 xvala = xdont(jwrkt(iinda))
3770 End If
3771 !
3772 End Do
3773 End Do
3774 !
3775 ! The Cs become As and Bs
3776 !
3777 lmtna = 2 * lmtna
3778 End Do
3779 !
3780 Return
3781 !
3782 End Subroutine i_mrgrnk
3783
3784 !! File ../src/m_ahelmholtz.f90
3785
3787 type(mg_t), intent(inout) :: mg
3788
3789 if (mg%n_extra_vars == 0 .and. mg%is_allocated) then
3790 error stop "ahelmholtz_set_methods: mg%n_extra_vars == 0"
3791 else
3792 mg%n_extra_vars = max(3, mg%n_extra_vars)
3793 end if
3794
3795 ! Use Neumann zero boundary conditions for the variable coefficient, since
3796 ! it is needed in ghost cells.
3797 mg%bc(:, mg_iveps1)%bc_type = mg_bc_neumann
3798 mg%bc(:, mg_iveps1)%bc_value = 0.0_dp
3799
3800 mg%bc(:, mg_iveps2)%bc_type = mg_bc_neumann
3801 mg%bc(:, mg_iveps2)%bc_value = 0.0_dp
3802
3803 mg%bc(:, mg_iveps3)%bc_type = mg_bc_neumann
3804 mg%bc(:, mg_iveps3)%bc_value = 0.0_dp
3805
3806 select case (mg%geometry_type)
3807 case (mg_cartesian)
3808 mg%box_op => box_ahelmh
3809
3810 select case (mg%smoother_type)
3812 mg%box_smoother => box_gs_ahelmh
3813 case default
3814 error stop "ahelmholtz_set_methods: unsupported smoother type"
3815 end select
3816 case default
3817 error stop "ahelmholtz_set_methods: unsupported geometry"
3818 end select
3819
3820 end subroutine ahelmholtz_set_methods
3821
3822 subroutine ahelmholtz_set_lambda(lambda)
3823 real(dp), intent(in) :: lambda
3824
3825 if (lambda < 0) &
3826 error stop "ahelmholtz_set_lambda: lambda < 0 not allowed"
3827
3828 ahelmholtz_lambda = lambda
3829 end subroutine ahelmholtz_set_lambda
3830
3831 !> Perform Gauss-Seidel relaxation on box for a Helmholtz operator
3832 subroutine box_gs_ahelmh(mg, id, nc, redblack_cntr)
3833 type(mg_t), intent(inout) :: mg
3834 integer, intent(in) :: id
3835 integer, intent(in) :: nc
3836 integer, intent(in) :: redblack_cntr !< Iteration counter
3837 integer :: i, j, k, i0, di
3838 logical :: redblack
3839 real(dp) :: idr2(2*3), u(2*3)
3840 real(dp) :: a0(2*3), a(2*3), c(2*3)
3841
3842 ! Duplicate 1/dr^2 array to multiply neighbor values
3843 idr2(1:2*3:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
3844 idr2(2:2*3:2) = idr2(1:2*3:2)
3845 i0 = 1
3846
3847 redblack = (mg%smoother_type == mg_smoother_gsrb)
3848 if (redblack) then
3849 di = 2
3850 else
3851 di = 1
3852 end if
3853
3854 ! The parity of redblack_cntr determines which cells we use. If
3855 ! redblack_cntr is even, we use the even cells and vice versa.
3856 associate(cc => mg%boxes(id)%cc, n => mg_iphi, &
3857 i_eps1 => mg_iveps1, &
3858 i_eps2 => mg_iveps2, &
3859 i_eps3 => mg_iveps3)
3860
3861 do k = 1, nc
3862 do j = 1, nc
3863 if (redblack) &
3864 i0 = 2 - iand(ieor(redblack_cntr, k+j), 1)
3865 do i = i0, nc, di
3866 a0(1:2) = cc(i, j, k, i_eps1)
3867 a0(3:4) = cc(i, j, k, i_eps2)
3868 a0(4:5) = cc(i, j, k, i_eps3)
3869 u(1:2) = cc(i-1:i+1:2, j, k, n)
3870 a(1:2) = cc(i-1:i+1:2, j, k, i_eps1)
3871 u(3:4) = cc(i, j-1:j+1:2, k, n)
3872 a(3:4) = cc(i, j-1:j+1:2, k, i_eps2)
3873 u(5:6) = cc(i, j, k-1:k+1:2, n)
3874 a(5:6) = cc(i, j, k-1:k+1:2, i_eps3)
3875 c(:) = 2 * a0(:) * a(:) / (a0(:) + a(:)) * idr2
3876
3877 cc(i, j, k, n) = (sum(c(:) * u(:)) - &
3878 cc(i, j, k, mg_irhs)) / &
3879 (sum(c(:)) + ahelmholtz_lambda)
3880 end do
3881 end do
3882 end do
3883 end associate
3884 end subroutine box_gs_ahelmh
3885
3886 !> Perform Helmholtz operator on a box
3887 subroutine box_ahelmh(mg, id, nc, i_out)
3888 type(mg_t), intent(inout) :: mg
3889 integer, intent(in) :: id
3890 integer, intent(in) :: nc
3891 integer, intent(in) :: i_out !< Index of variable to store Helmholtz in
3892 integer :: i, j, k
3893 real(dp) :: idr2(2*3), a0(2*3), u0, u(2*3), a(2*3)
3894
3895 ! Duplicate 1/dr^2 array to multiply neighbor values
3896 idr2(1:2*3:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
3897 idr2(2:2*3:2) = idr2(1:2*3:2)
3898
3899 associate(cc => mg%boxes(id)%cc, n => mg_iphi, &
3900 i_eps1 => mg_iveps1, &
3901 i_eps2 => mg_iveps2, &
3902 i_eps3 => mg_iveps3)
3903 do k = 1, nc
3904 do j = 1, nc
3905 do i = 1, nc
3906 u0 = cc(i, j, k, n)
3907 a0(1:2) = cc(i, j, k, i_eps1)
3908 a0(3:4) = cc(i, j, k, i_eps2)
3909 a0(5:6) = cc(i, j, k, i_eps3)
3910 u(1:2) = cc(i-1:i+1:2, j, k, n)
3911 u(3:4) = cc(i, j-1:j+1:2, k, n)
3912 u(5:6) = cc(i, j, k-1:k+1:2, n)
3913 a(1:2) = cc(i-1:i+1:2, j, k, i_eps1)
3914 a(3:4) = cc(i, j-1:j+1:2, k, i_eps2)
3915 a(5:6) = cc(i, j, k-1:k+1:2, i_eps3)
3916
3917 cc(i, j, k, i_out) = sum(2 * idr2 * &
3918 a0(:)*a(:)/(a0(:) + a(:)) * (u(:) - u0)) - &
3920 end do
3921 end do
3922 end do
3923 end associate
3924 end subroutine box_ahelmh
3925
3926end module m_octree_mg_3d
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.
subroutine, public helmholtz_set_lambda(lambda)
subroutine, public mg_set_methods(mg)
subroutine, public mg_restrict_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
integer, parameter, public mg_irhs
Index of right-hand side.
integer, parameter, public mg_vhelmholtz
Indicates a variable-coefficient Helmholtz equation.
integer, dimension(4, 6), parameter, public mg_child_adj_nb
integer, parameter, public mg_ahelmholtz
Indicates a anisotropic-coefficient Helmholtz equation.
integer, dimension(3, 6), parameter, public mg_neighb_dix
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...
subroutine, public mg_fas_fmg(mg, have_guess, max_res)
Perform FAS-FMG cycle (full approximation scheme, full multigrid).
integer, parameter, public mg_smoother_gsrb
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_num_children
subroutine, public mg_timers_show(mg)
integer, parameter, public mg_max_timers
Maximum number of timers to use.
subroutine, public mg_fas_vcycle(mg, highest_lvl, max_res, standalone)
Perform FAS V-cycle (full approximation scheme).
subroutine, public vhelmholtz_set_methods(mg)
integer, parameter, public mg_bc_dirichlet
Value to indicate a Dirichlet boundary condition.
integer, parameter, public mg_iveps
Index of the variable coefficient (at cell centers)
subroutine, public mg_load_balance_parents(mg)
Load balance the parents (non-leafs). Assign them to the rank that has most children.
integer, dimension(3, 8), parameter, public mg_child_dix
subroutine, public mg_phi_bc_store(mg)
Store boundary conditions for the solution variable, this can speed up calculations if the same bound...
integer(i8) function, public mg_number_of_unknowns(mg)
Determine total number of unknowns (on leaves)
integer, parameter, public mg_vlaplacian
Indicates a variable-coefficient Laplacian.
integer, parameter, public mg_cylindrical
Cylindrical coordinate system.
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_build_rectangle(mg, domain_size, box_size, dx, r_min, periodic, n_finer)
subroutine, public vhelmholtz_set_lambda(lambda)
integer, parameter, public mg_ires
Index of residual.
integer, parameter, public mg_smoother_gs
subroutine, public laplacian_set_methods(mg)
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_ghost_cell_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
subroutine, public mg_set_neighbors_lvl(mg, lvl)
integer, parameter, public mg_smoother_jacobi
integer, parameter, public mg_neighb_lowy
integer, parameter, public mg_bc_continuous
Value to indicate a continuous boundary condition.
subroutine, public mg_set_leaves_parents(boxes, level)
Create a list of leaves and a list of parents for a level.
subroutine, public helmholtz_set_methods(mg)
subroutine, public mg_add_children(mg, id)
integer, parameter, public mg_lvl_lo
Minimum allowed grid level.
integer, parameter, public mg_num_neighbors
integer, parameter, public mg_neighb_highy
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_iold
Index of previous solution (used for correction)
integer, parameter, public mg_bc_neumann
Value to indicate a Neumann boundary condition.
integer, dimension(6), parameter, public mg_neighb_high_pm
integer, parameter, public mg_ndim
Problem dimension.
integer, parameter, public mg_iphi
Index of solution.
integer, parameter, public mg_helmholtz
Indicates a constant-coefficient Helmholtz equation.
subroutine, public mg_allocate_storage(mg)
Allocate communication buffers and local boxes for a tree that has already been created.
subroutine, public mg_timer_end(timer)
integer, parameter, public mg_num_vars
Number of predefined multigrid variables.
subroutine, public ahelmholtz_set_methods(mg)
real(dp), public, protected vhelmholtz_lambda
Module for implicitly solving diffusion equations.
subroutine, public mg_restrict(mg, iv)
Restrict all levels.
integer, parameter, public mg_neighb_lowx
subroutine, public mg_set_next_level_ids(mg, lvl)
integer, parameter, public mg_physical_boundary
Special value that indicates there is a physical boundary.
integer, dimension(6), parameter, public mg_neighb_rev
integer, parameter, public mg_spherical
Spherical coordinate system.
integer, parameter, public mg_laplacian
Indicates a standard Laplacian.
subroutine, public mg_get_face_coords(box, nb, nc, x)
Get coordinates at the face of a box.
subroutine, public mg_timer_start(timer)
subroutine, public mg_prolong_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
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...
integer, parameter, public mg_neighb_highx
integer, parameter, public mg_cartesian
Cartesian coordinate system.
integer, parameter, public mg_iveps2
integer, parameter, public i8
Type for 64-bit integers.
integer, parameter, public mg_iveps1
Indexes of anisotropic variable coefficients.
subroutine, public mg_fill_ghost_cells_lvl(mg, lvl, iv)
Fill ghost cells at a grid level.
integer, parameter, public mg_neighb_lowz
subroutine, public diffusion_solve_vcoeff(mg, dt, order, max_res)
Solve a diffusion equation implicitly, assuming a variable diffusion coefficient which has been store...
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 ...
integer, parameter, public mg_lvl_hi
Maximum allowed grid level.
subroutine, public vlaplacian_set_methods(mg)
logical, dimension(6), parameter, public mg_neighb_low
integer, parameter, public mg_iveps3
elemental logical function, public mg_has_children(box)
Return .true. if a box has children.
subroutine, public mg_deallocate_storage(mg)
Deallocate all allocatable arrays.
real(dp), public, protected helmholtz_lambda
Module for load balancing a tree (that already has been constructed). The load balancing determines w...
subroutine, public mg_set_refinement_boundaries(boxes, level)
Create a list of refinement boundaries (from the coarse side)
pure integer function, dimension(3), 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_no_box
Special value that indicates there is no box.
subroutine, public mg_restrict_lvl(mg, iv, lvl)
Restrict from lvl to lvl-1.
subroutine, public diffusion_solve_acoeff(mg, dt, order, max_res)
Solve a diffusion equation implicitly, assuming anisotropic diffusion coefficient which has been stor...
pure integer function, public mg_highest_uniform_lvl(mg)
integer, dimension(6), parameter, public mg_neighb_dim
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 ahelmholtz_set_lambda(lambda)
subroutine, public mg_fill_ghost_cells(mg, iv)
Fill ghost cells at all grid levels.
integer function, public mg_add_timer(mg, name)
subroutine, public mg_apply_op(mg, i_out, op)
Apply operator to the tree and store in variable i_out.
integer, dimension(8, 3), parameter, public mg_child_rev
logical, dimension(3, 8), parameter, public mg_child_low
integer, parameter, public mg_max_num_vars
Maximum number of variables.
integer, parameter, public dp
Type of reals.
subroutine, public sort_and_transfer_buffers(mg, dsize)
integer, parameter, public mg_neighb_highz
Box data structure.
Buffer type (one is used for each pair of communicating processes)
Lists of blocks per refinement level.