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