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