MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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
238  type(mg_lvl_t) :: lvls(mg_lvl_lo:mg_lvl_hi)
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
419  public :: mg_load_balance_parents
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
432  public :: sort_and_transfer_buffers
433 
434  !! File ../src/m_ghost_cells.f90
435 
436  ! Public methods
437  public :: mg_ghost_cell_buffer_size
438  public :: mg_fill_ghost_cells
439  public :: mg_fill_ghost_cells_lvl
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
481  public :: mg_restrict_buffer_size
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
504 contains
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
641  subroutine mg_deallocate_storage(mg)
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 
889  subroutine laplacian_set_methods(mg)
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 
1026  subroutine vhelmholtz_set_methods(mg)
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)) - &
1146  vhelmholtz_lambda * 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.
1541  subroutine mg_load_balance_simple(mg)
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 
1758  subroutine vlaplacian_set_methods(mg)
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 
2881  subroutine helmholtz_set_methods(mg)
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 
3786  subroutine ahelmholtz_set_methods(mg)
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)) - &
3919  ahelmholtz_lambda * u0
3920  end do
3921  end do
3922  end do
3923  end associate
3924  end subroutine box_ahelmh
3925 
3926 end module m_octree_mg_3d
To fill ghost cells near physical 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 subtract_mean(mg, iv, include_ghostcells)
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.
real(dp) function get_sum(mg, iv)
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 sort_sendbuf(gc, dsize)
Sort send buffers according to the idbuf array.
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.