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