MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_multigrid_coupling.t
Go to the documentation of this file.
1 !> Module to couple the octree-mg library to AMRVAC. This file uses the VACPP
2 !> preprocessor, but its use is kept to a minimum.
3 {^ifoned
4 !> In 1D, only provide a dummy module.
6  implicit none
7  public
8 
9  type mg_t
10  end type mg_t
11 
12 contains
13 
14  !> Setup multigrid for usage
15  subroutine mg_setup_multigrid()
16  error stop "Multigrid not available in 1D"
17  end subroutine mg_setup_multigrid
18 
19  !> If the grid has changed, rebuild the full multigrid tree
20  subroutine mg_update_refinement(n_coarsen, n_refine)
21  integer, intent(in) :: n_coarsen
22  integer, intent(in) :: n_refine
23  error stop "Multigrid not available in 1D"
24  end subroutine mg_update_refinement
25 
26  !> Copy a variable to the multigrid tree, including a layer of ghost cells
27  subroutine mg_copy_to_tree(iw_from, iw_to, restrict, restrict_gc)
28  integer, intent(in) :: iw_from !< Variable to use as right-hand side
29  integer, intent(in) :: iw_to !< Copy to this variable
30  logical, intent(in) :: restrict !< Restrict variable on multigrid tree
31  logical, intent(in) :: restrict_gc !< Fill ghost cells after restrict
32 
33  error stop "Multigrid not available in 1D"
34  end subroutine mg_copy_to_tree
35 
36 end module mod_multigrid_coupling
37 }
38 {^nooned
40  {^iftwod
41  use m_octree_mg_2d
42  }
43  {^ifthreed
44  use m_octree_mg_3d
45  }
46 
47 
48  implicit none
49  public
50 
51  !> Data structure containing the multigrid tree.
52  type(mg_t) :: mg
53 
54  !> If defined, this routine is called after a new multigrid tree is
55  !> constructed.
56  procedure(after_new_tree), pointer :: mg_after_new_tree => null()
57 
58  interface
59  subroutine after_new_tree()
60  end subroutine after_new_tree
61  end interface
62 
63 contains
64 
65  !> Setup multigrid for usage
66  subroutine mg_setup_multigrid()
68  use mod_geometry
69 
70  if (ndim == 1) &
71  error stop "Multigrid not available in 1D"
72 
73  if (ndim /= mg_ndim) &
74  error stop "Multigrid module was compiled for different ndim"
75 
76  select case (coordinate)
77  case (cartesian)
78  if (ndim == 1) error stop "Multigrid only support 2D, 3D"
79  case (cylindrical)
80  if (ndim == 3) error stop "Multigrid does not support cylindrical 3D"
81  mg%geometry_type = mg_cylindrical
82  case default
83  error stop "Multigrid does not support your geometry"
84  end select
85 
86  if (any([ block_nx^d ] /= block_nx1)) &
87  error stop "Multigrid requires all block_nx to be equal"
88 
89  call mg_comm_init(mg)
90  call mg_set_methods(mg)
92  end subroutine mg_setup_multigrid
93 
94  !> Set multigrid boundary conditions for the solution according to variable iw
95  subroutine mg_copy_boundary_conditions(mg, iw)
97  type(mg_t), intent(inout) :: mg
98  integer, intent(in) :: iw
99  character(len=std_len) :: bnd_name(mg_num_neighbors)
100  integer :: n
101 
102  do n = 1, mg_num_neighbors
103  select case (typeboundary(iw, n))
104  case ('symm')
105  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
106  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
107  case ('asymm')
108  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
109  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
110  case ('cont')
111  mg%bc(n, mg_iphi)%bc_type = mg_bc_continuous
112  mg%bc(n, mg_iphi)%bc_value = 0.0_dp ! Not needed
113  case ('periodic')
114  ! Nothing to do here
115  case default
116  print *, "Not a standard: ", trim(typeboundary(iw, n))
117  error stop "You have to set a user-defined boundary method"
118  end select
119  end do
120  end subroutine mg_copy_boundary_conditions
121 
122  !> If the grid has changed, rebuild the full multigrid tree
123  subroutine mg_update_refinement(n_coarsen, n_refine)
125  integer, intent(in) :: n_coarsen
126  integer, intent(in) :: n_refine
127 
128  ! Don't build multigrid tree while doing initial refinement
129  if (.not. time_advance) return
130 
131  if (.not. mg%is_allocated) then
132  call mg_tree_from_amrvac(mg)
133  else if (n_coarsen + n_refine > 0) then
135  call mg_tree_from_amrvac(mg)
136  end if
137  end subroutine mg_update_refinement
138 
139  !> Copy a variable to the multigrid tree, including a layer of ghost cells
140  subroutine mg_copy_to_tree(iw_from, iw_to, restrict, restrict_gc, factor)
142  use mod_forest
143  integer, intent(in) :: iw_from !< Variable to use as right-hand side
144  integer, intent(in) :: iw_to !< Copy to this variable
145  logical, intent(in), optional :: restrict !< Restrict variable on multigrid tree
146  logical, intent(in), optional :: restrict_gc !< Fill ghost cells after restrict
147  real(dp), intent(in), optional :: factor !< out = factor * in
148  integer :: iigrid, igrid, id
149  integer :: nc, lvl
150  type(tree_node), pointer :: pnode
151  real(dp) :: fac
152  logical :: do_restrict, do_gc
153 
154  if (.not. mg%is_allocated) &
155  error stop "mg_copy_to_tree: tree not allocated yet"
156 
157  fac = 1.0_dp; if (present(factor)) fac = factor
158  do_restrict = .false.; if (present(restrict)) do_restrict = restrict
159  do_gc = .false.; if (present(restrict_gc)) do_gc = restrict_gc
160 
161  do iigrid = 1, igridstail
162  igrid = igrids(iigrid);
163  pnode => igrid_to_node(igrid, mype)%node
164  id = pnode%id
165  lvl = mg%boxes(id)%lvl
166  nc = mg%box_size_lvl(lvl)
167 
168  ! Include one layer of ghost cells on grid leaves
169  {^iftwod
170  mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = fac * &
171  ps(igrid)%w(ixmlo1-1:ixmhi1+1, ixmlo2-1:ixmhi2+1, iw_from)
172  }
173  {^ifthreed
174  mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
175  ps(igrid)%w(ixmlo1-1:ixmhi1+1, ixmlo2-1:ixmhi2+1, &
176  ixmlo3-1:ixmhi3+1, iw_from)
177  }
178  end do
179 
180  if (do_restrict) then
181  call mg_restrict(mg, iw_to)
182  if (do_gc) call mg_fill_ghost_cells(mg, iw_to)
183  end if
184 
185  end subroutine mg_copy_to_tree
186 
187 ! !> Copy a variable to the multigrid tree with a single layer of ghost cells
188 ! subroutine mg_copy_to_tree_gc(iw_from, iw_to)
189 ! use mod_global_parameters
190 ! use mod_forest
191 ! integer, intent(in) :: iw_from !< Variable to use as right-hand side
192 ! integer, intent(in) :: iw_to !< Copy to this variable
193 ! integer :: iigrid, igrid, id
194 ! integer :: nc, lvl
195 ! type(tree_node), pointer :: pnode
196 
197 ! if (.not. mg%is_allocated) &
198 ! error stop "mg_copy_to_tree: tree not allocated yet"
199 
200 ! do iigrid = 1, igridstail
201 ! igrid = igrids(iigrid);
202 ! pnode => igrid_to_node(igrid, mype)%node
203 ! id = pnode%id
204 ! lvl = mg%boxes(id)%lvl
205 ! nc = mg%box_size_lvl(lvl)
206 
207 ! {^IFTWOD
208 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = &
209 ! ps(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_from)
210 ! }
211 ! {^IFTHREED
212 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = &
213 ! ps(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
214 ! ixMlo3-1:ixMhi3+1, iw_from)
215 ! }
216 ! end do
217 ! end subroutine mg_copy_to_tree_gc
218 
219  !> Copy a variable from the multigrid tree
220  subroutine mg_copy_from_tree(iw_from, iw_to)
222  use mod_forest
223  integer, intent(in) :: iw_from !< Variable to use as right-hand side
224  integer, intent(in) :: iw_to !< Copy to this variable
225  integer :: iigrid, igrid, id
226  integer :: nc, lvl
227  type(tree_node), pointer :: pnode
228 
229  if (.not. mg%is_allocated) &
230  error stop "mg_copy_from_tree: tree not allocated yet"
231 
232  do iigrid = 1, igridstail
233  igrid = igrids(iigrid);
234  pnode => igrid_to_node(igrid, mype)%node
235  id = pnode%id
236  lvl = mg%boxes(id)%lvl
237  nc = mg%box_size_lvl(lvl)
238 
239 {^iftwod
240  ps(igrid)%w(ixmlo1:ixmhi1, ixmlo2:ixmhi2, iw_to) = &
241  mg%boxes(id)%cc(1:nc, 1:nc, iw_from)
242 }
243 {^ifthreed
244  ps(igrid)%w(ixmlo1:ixmhi1, ixmlo2:ixmhi2, ixmlo3:ixmhi3, iw_to) = &
245  mg%boxes(id)%cc(1:nc, 1:nc, 1:nc, iw_from)
246 }
247  end do
248  end subroutine mg_copy_from_tree
249 
250  !> Copy from multigrid tree with one layer of ghost cells. Corner ghost cells
251  !> are not used/set.
252  subroutine mg_copy_from_tree_gc(iw_from, iw_to)
254  use mod_forest
255  integer, intent(in) :: iw_from !< Variable to use as right-hand side
256  integer, intent(in) :: iw_to !< Copy to this variable
257  integer :: iigrid, igrid, id
258  integer :: nc, lvl
259  type(tree_node), pointer :: pnode
260 
261  if (.not. mg%is_allocated) &
262  error stop "mg_copy_from_tree_gc: tree not allocated yet"
263 
264  do iigrid = 1, igridstail
265  igrid = igrids(iigrid);
266  pnode => igrid_to_node(igrid, mype)%node
267  id = pnode%id
268  lvl = mg%boxes(id)%lvl
269  nc = mg%box_size_lvl(lvl)
270 
271 {^iftwod
272  ps(igrid)%w(ixmlo1-1:ixmhi1+1, ixmlo2-1:ixmhi2+1, iw_to) = &
273  mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_from)
274 }
275 {^ifthreed
276  ps(igrid)%w(ixmlo1-1:ixmhi1+1, ixmlo2-1:ixmhi2+1, &
277  ixmlo3-1:ixmhi3+1, iw_to) = &
278  mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_from)
279 }
280  end do
281  end subroutine mg_copy_from_tree_gc
282 
283  !> Generate a multigrid tree that includes the amrvac tree, but also contains
284  !> coarser grid levels. A number of checks has already been performed in
285  !> mg_setup_multigrid, so we don't repeat these checks here.
286  subroutine mg_tree_from_amrvac(mg)
289  type(mg_t), intent(inout) :: mg
290  integer :: i, n, id, ix(ndim)
291  integer :: n_boxes_total, i_c, c_id, c_ix(ndim)
292  integer :: min_lvl, lvl
293  integer :: nb, nb_ix, nb_dim
294  integer :: n_finer
295  type(tree_node), pointer :: pnode, pnode_ch
296  type(tree_node_ptr), allocatable :: id_to_node(:)
297  real(dp) :: dr_coarse
298 
299  ! Estimate number of finer blocks
300  n_finer = nparents+nleafs
301 
302  call mg_build_rectangle(mg, [ domain_nx^d ], block_nx1, dx(:,1), &
303  [ xprobmin^d ], periodb, n_finer)
304 
305  mg%highest_lvl = levmax
306  n_boxes_total = mg%n_boxes + n_finer
307 
308  ! To link the two trees
309  allocate(id_to_node(n_boxes_total))
310 
311  ! Link base level
312  do i = 1, size(mg%lvls(1)%ids)
313  id = mg%lvls(1)%ids(i)
314  ix = mg%boxes(id)%ix
315 
316  pnode => tree_root({ix(^d)})%node
317  pnode%id = id
318  id_to_node(id)%node => pnode
319  mg%boxes(id)%rank = pnode%ipe
320  end do
321 
322  ! Add refinement
323  do lvl = 1, mg%highest_lvl
324  do i = 1, size(mg%lvls(lvl)%ids)
325  id = mg%lvls(lvl)%ids(i)
326  pnode => id_to_node(id)%node
327 
328  if (.not. pnode%leaf) then
329  call mg_add_children(mg, id)
330 
331  do i_c = 1, mg_num_children
332  c_id = mg%boxes(id)%children(i_c)
333  c_ix = mg_child_dix(:, i_c) + 1
334  pnode_ch => pnode%child({c_ix(^d)})%node
335  id_to_node(c_id)%node => pnode_ch
336  pnode_ch%id = c_id
337  mg%boxes(c_id)%rank = pnode_ch%ipe
338  end do
339  end if
340  end do
341 
342  call mg_set_leaves_parents(mg%boxes, mg%lvls(lvl))
343 
344  if (lvl < mg%highest_lvl) then
345  call mg_set_next_level_ids(mg, lvl)
346  call mg_set_neighbors_lvl(mg, lvl+1)
347  end if
348  end do
349 
350  ! Store boxes with refinement boundaries (from the coarse side)
351  do lvl = 1, mg%highest_lvl
352  call mg_set_refinement_boundaries(mg%boxes, mg%lvls(lvl))
353  end do
354 
355  ! Assign boxes to MPI processes
356  call mg_load_balance_parents(mg)
357 
358  ! Allocate storage for boxes owned by this process
359  call mg_allocate_storage(mg)
360 
361  if (associated(mg_after_new_tree)) then
362  call mg_after_new_tree()
363  end if
364 
365  end subroutine mg_tree_from_amrvac
366 
367 end module mod_multigrid_coupling
368 }
integer, parameter cylindrical
Definition: mod_geometry.t:9
subroutine mg_update_refinement(n_coarsen, n_refine)
If the grid has changed, rebuild the full multigrid tree.
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer domain_nx
number of cells for each dimension in level-one mesh
The data structure that contains information about a tree node/grid block.
Definition: mod_forest.t:11
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition: mod_forest.t:32
subroutine, public mg_load_balance_parents(mg)
Load balance the parents (non-leafs). Assign them to the rank that has most children.
subroutine mg_tree_from_amrvac(mg)
Generate a multigrid tree that includes the amrvac tree, but also contains coarser grid levels...
subroutine, public mg_fill_ghost_cells(mg, iv)
Fill ghost cells at all grid levels.
subroutine, public mg_set_refinement_boundaries(boxes, level)
Create a list of refinement boundaries (from the coarse side)
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine, public mg_set_methods(mg)
integer coordinate
Definition: mod_geometry.t:6
subroutine, public mg_allocate_storage(mg)
Allocate communication buffers and local boxes for a tree that has already been created.
Module with basic grid data structures.
Definition: mod_forest.t:2
subroutine mg_copy_from_tree(iw_from, iw_to)
Copy a variable from the multigrid tree.
integer, parameter, public mg_bc_dirichlet
Value to indicate a Dirichlet boundary condition.
integer, parameter, public mg_num_children
type(mg_t) mg
Data structure containing the multigrid tree.
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor, but its use is kept to a minimum.
integer, parameter, public mg_iphi
Index of solution.
Pointer to a tree_node.
Definition: mod_forest.t:6
subroutine, public mg_set_leaves_parents(boxes, level)
Create a list of leaves and a list of parents for a level.
subroutine mg_copy_boundary_conditions(mg, iw)
Set multigrid boundary conditions for the solution according to variable iw.
integer, parameter, public mg_num_neighbors
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 ...
type(tree_node_ptr), dimension(:^d &), allocatable, save tree_root
Pointers to the coarse grid.
Definition: mod_forest.t:29
integer, parameter cartesian
Definition: mod_geometry.t:7
character(len=std_len), dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
subroutine, public mg_deallocate_storage(mg)
Deallocate all allocatable arrays.
subroutine, public mg_add_children(mg, id)
double precision, dimension(:,:), allocatable dx
integer, dimension(:), allocatable, parameter d
integer mype
The rank of the current MPI task.
procedure(after_new_tree), pointer mg_after_new_tree
If defined, this routine is called after a new multigrid tree is constructed.
integer, parameter, public mg_cylindrical
Cylindrical coordinate system.
integer block_nx
number of cells for each dimension in grid block excluding ghostcells
subroutine, public mg_restrict(mg, iv)
Restrict all levels.
subroutine mg_copy_to_tree(iw_from, iw_to, restrict, restrict_gc)
Copy a variable to the multigrid tree, including a layer of ghost cells.
subroutine mg_setup_multigrid()
Setup multigrid for usage.
integer, parameter, public mg_ndim
Problem dimension.
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine mg_copy_from_tree_gc(iw_from, iw_to)
Copy from multigrid tree with one layer of ghost cells. Corner ghost cells are not used/set...
subroutine, public mg_set_next_level_ids(mg, lvl)
subroutine, public mg_build_rectangle(mg, domain_size, box_size, dx, r_min, periodic, n_finer)
integer, parameter, public mg_bc_continuous
Value to indicate a continuous boundary condition.
integer, parameter, public mg_bc_neumann
Value to indicate a Neumann boundary condition.
integer, save nparents
Number of parent blocks.
Definition: mod_forest.t:73
logical, dimension(ndim) periodb
True for dimensions with periodic boundaries.
integer, save nleafs
Number of leaf block.
Definition: mod_forest.t:76
subroutine, public mg_set_neighbors_lvl(mg, lvl)
integer, dimension(2, 4), parameter, public mg_child_dix