MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
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
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
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
421
422 !! File ../src/m_ghost_cells.f90
423
424 ! Public methods
426 public :: mg_fill_ghost_cells
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
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
492contains
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
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
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
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)) - &
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.
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
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
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
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)) - &
3691 end do
3692 end associate
3693 end subroutine box_ahelmh
3694
3695end module m_octree_mg_1d
Subroutine that performs Gauss-Seidel relaxation.
Subroutine that performs A * cc(..., i_in) = cc(..., i_out)
Subroutine that performs prolongation to a single child.
To fill ghost cells near physical boundaries.
To fill ghost cells near refinement boundaries.
integer(i8) function, public mg_number_of_unknowns(mg)
Determine total number of unknowns (on leaves)
integer, parameter, public mg_iphi
Index of solution.
integer, parameter, public i8
Type for 64-bit integers.
subroutine, public mg_prolong(mg, lvl, iv, iv_to, method, add)
Prolong variable iv from lvl to variable iv_to at lvl+1.
subroutine, public mg_fill_ghost_cells_lvl(mg, lvl, iv)
Fill ghost cells at a grid level.
integer, parameter, public mg_iveps3
integer, dimension(2, 1), parameter, public mg_child_rev
real(dp), public, protected ahelmholtz_lambda
Module which contains multigrid procedures for a Helmholtz operator of the form: div(D grad(phi)) - l...
integer, parameter, public mg_bc_neumann
Value to indicate a Neumann boundary condition.
logical, dimension(1, 2), parameter, public mg_child_low
integer, parameter, public mg_cylindrical
Cylindrical coordinate system.
integer, parameter, public mg_num_children
subroutine, public mg_restrict_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
integer, parameter, public mg_lvl_hi
Maximum allowed grid level.
integer, parameter, public mg_laplacian
Indicates a standard Laplacian.
integer, parameter, public mg_cartesian
Cartesian coordinate system.
subroutine, public mg_apply_op(mg, i_out, op)
Apply operator to the tree and store in variable i_out.
integer, parameter, public mg_physical_boundary
Special value that indicates there is a physical boundary.
elemental logical function, public mg_has_children(box)
Return .true. if a box has children.
subroutine, public mg_timer_end(timer)
integer, parameter, public mg_num_vars
Number of predefined multigrid variables.
integer, parameter, public mg_smoother_gsrb
integer, parameter, public mg_iveps
Index of the variable coefficient (at cell centers)
subroutine, public mg_set_methods(mg)
integer, parameter, public mg_ahelmholtz
Indicates a anisotropic-coefficient Helmholtz equation.
integer, dimension(1, 2), parameter, public mg_neighb_dix
subroutine, public vhelmholtz_set_methods(mg)
real(dp), public, protected vhelmholtz_lambda
Module for implicitly solving diffusion equations.
integer, parameter, public mg_no_box
Special value that indicates there is no box.
integer, parameter, public mg_lvl_lo
Minimum allowed grid level.
subroutine, public mg_restrict_lvl(mg, iv, lvl)
Restrict from lvl to lvl-1.
subroutine, public helmholtz_set_methods(mg)
subroutine, public mg_fill_ghost_cells(mg, iv)
Fill ghost cells at all grid levels.
integer, parameter, public mg_smoother_jacobi
subroutine, public mg_build_rectangle(mg, domain_size, box_size, dx, r_min, periodic, n_finer)
pure integer function, dimension(1), public mg_get_child_offset(mg, id)
Get the offset of a box with respect to its parent (e.g. in 2d, there can be a child at offset 0,...
subroutine, public mg_comm_init(mg, comm)
Prolong from a parent to a child with index offset dix. This method could sometimes work better than ...
subroutine, public vhelmholtz_set_lambda(lambda)
subroutine, public mg_ghost_cell_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
subroutine, public mg_timers_show(mg)
integer, parameter, public mg_max_num_vars
Maximum number of variables.
integer function, public mg_ix_to_ichild(ix)
Compute the child index for a box with spatial index ix. With child index we mean the index in the ch...
real(dp), public, protected helmholtz_lambda
Module for load balancing a tree (that already has been constructed). The load balancing determines w...
subroutine, public diffusion_solve(mg, dt, diffusion_coeff, order, max_res)
Solve a diffusion equation implicitly, assuming a constant diffusion coefficient. The solution at tim...
subroutine, public mg_restrict(mg, iv)
Restrict all levels.
subroutine, public mg_load_balance_simple(mg)
Load balance all boxes in the multigrid tree, by simply distributing the load per grid level....
subroutine, public mg_fas_fmg(mg, have_guess, max_res)
Perform FAS-FMG cycle (full approximation scheme, full multigrid).
subroutine, public helmholtz_set_lambda(lambda)
integer, parameter, public mg_iveps1
Indexes of anisotropic variable coefficients.
subroutine, public mg_add_children(mg, id)
subroutine, public laplacian_set_methods(mg)
subroutine, public mg_set_refinement_boundaries(boxes, level)
Create a list of refinement boundaries (from the coarse side)
integer, dimension(2), parameter, public mg_neighb_high_pm
subroutine, public mg_prolong_buffer_size(mg, n_send, n_recv, dsize)
Specify minimum buffer size (per process) for communication.
subroutine, public mg_set_next_level_ids(mg, lvl)
subroutine, public mg_set_neighbors_lvl(mg, lvl)
integer, dimension(2), parameter, public mg_neighb_rev
integer, parameter, public dp
Type of reals.
integer, parameter, public mg_bc_dirichlet
Value to indicate a Dirichlet boundary condition.
integer, parameter, public mg_iveps2
integer, parameter, public mg_vlaplacian
Indicates a variable-coefficient Laplacian.
integer, parameter, public mg_helmholtz
Indicates a constant-coefficient Helmholtz equation.
pure integer function, public mg_highest_uniform_lvl(mg)
subroutine, public mg_allocate_storage(mg)
Allocate communication buffers and local boxes for a tree that has already been created.
subroutine, public diffusion_solve_acoeff(mg, dt, order, max_res)
Solve a diffusion equation implicitly, assuming anisotropic diffusion coefficient which has been stor...
integer, parameter, public mg_ndim
Problem dimension.
logical, dimension(2), parameter, public mg_neighb_low
integer, parameter, public mg_smoother_gs
subroutine, public mg_timer_start(timer)
subroutine, public vlaplacian_set_methods(mg)
integer, parameter, public mg_neighb_highx
integer, parameter, public mg_irhs
Index of right-hand side.
subroutine, public mg_deallocate_storage(mg)
Deallocate all allocatable arrays.
subroutine, public mg_get_face_coords(box, nb, nc, x)
Get coordinates at the face of a box.
integer, parameter, public mg_neighb_lowx
subroutine, public mg_prolong_sparse(mg, p_id, dix, nc, iv, fine)
Prolong from a parent to a child with index offset dix.
subroutine, public mg_set_leaves_parents(boxes, level)
Create a list of leaves and a list of parents for a level.
subroutine, public mg_load_balance_parents(mg)
Load balance the parents (non-leafs). Assign them to the rank that has most children.
integer, dimension(1, 2), parameter, public mg_child_adj_nb
subroutine, public ahelmholtz_set_lambda(lambda)
integer, dimension(2), parameter, public mg_neighb_dim
subroutine, public mg_load_balance(mg)
Load balance all boxes in the multigrid tree. Compared to mg_load_balance_simple, this method does a ...
integer, parameter, public mg_spherical
Spherical coordinate system.
integer function, public mg_add_timer(mg, name)
subroutine, public mg_fas_vcycle(mg, highest_lvl, max_res, standalone)
Perform FAS V-cycle (full approximation scheme).
integer, dimension(1, 2), parameter, public mg_child_dix
integer, parameter, public mg_bc_continuous
Value to indicate a continuous boundary condition.
integer, parameter, public mg_max_timers
Maximum number of timers to use.
integer, parameter, public mg_num_neighbors
subroutine, public ahelmholtz_set_methods(mg)
integer, parameter, public mg_ires
Index of residual.
integer, parameter, public mg_iold
Index of previous solution (used for correction)
subroutine, public sort_and_transfer_buffers(mg, dsize)
subroutine, public diffusion_solve_vcoeff(mg, dt, order, max_res)
Solve a diffusion equation implicitly, assuming a variable diffusion coefficient which has been store...
integer, parameter, public mg_vhelmholtz
Indicates a variable-coefficient Helmholtz equation.
subroutine, public mg_phi_bc_store(mg)
Store boundary conditions for the solution variable, this can speed up calculations if the same bound...
Box data structure.
Buffer type (one is used for each pair of communicating processes)
Lists of blocks per refinement level.