MPI-AMRVAC  3.1
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.
4  {^ifoned
6  }
7  {^iftwod
9  }
10  {^ifthreed
11  use m_octree_mg_3d
12  }
13 
14 
15  implicit none
16  public
17 
18  !> Data structure containing the multigrid tree.
19  type(mg_t) :: mg
20 
21  !> If defined, this routine is called after a new multigrid tree is
22  !> constructed.
23  procedure(after_new_tree), pointer :: mg_after_new_tree => null()
24 
25  interface
26  subroutine after_new_tree()
27  end subroutine after_new_tree
28  end interface
29 
30 contains
31 
32  !> Setup multigrid for usage
33  subroutine mg_setup_multigrid()
35  use mod_geometry
36 
37  if (ndim /= mg_ndim) &
38  error stop "Multigrid module was compiled for different ndim"
39 
40  select case (coordinate)
41  case (cartesian)
42  continue
43  case (cylindrical)
44  if (ndim == 3) error stop "Multigrid does not support cylindrical 3D"
45  mg%geometry_type = mg_cylindrical
46  case default
47  error stop "Multigrid does not support your geometry"
48  end select
49 
50  if (any([ block_nx^d ] /= block_nx1)) &
51  error stop "Multigrid requires all block_nx to be equal"
52 
53  call mg_comm_init(mg)
54  call mg_set_methods(mg)
56  end subroutine mg_setup_multigrid
57 
58  !> Set multigrid boundary conditions for the solution according to variable iw
59  subroutine mg_copy_boundary_conditions(mg, iw)
61  type(mg_t), intent(inout) :: mg
62  integer, intent(in) :: iw
63  character(len=std_len) :: bnd_name(mg_num_neighbors)
64  integer :: n
65 
66  do n = 1, mg_num_neighbors
67  select case (typeboundary(iw, n))
68  case (bc_symm)
69  mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
70  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
71  case (bc_asymm)
72  mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
73  mg%bc(n, mg_iphi)%bc_value = 0.0_dp
74  case (bc_cont)
75  mg%bc(n, mg_iphi)%bc_type = mg_bc_continuous
76  mg%bc(n, mg_iphi)%bc_value = 0.0_dp ! Not needed
77  case (bc_periodic)
78  ! Nothing to do here
79  case default
80  write(*,*) "mg_copy_boundary_conditions: unknown boundary type"
81  end select
82  end do
83  end subroutine mg_copy_boundary_conditions
84 
85  !> If the grid has changed, rebuild the full multigrid tree
86  subroutine mg_update_refinement(n_coarsen, n_refine)
88  integer, intent(in) :: n_coarsen
89  integer, intent(in) :: n_refine
90 
91  ! Don't build multigrid tree while doing initial refinement
92  if (.not. time_advance) return
93 
94  if (.not. mg%is_allocated) then
96  else if (n_coarsen + n_refine > 0) then
99  end if
100  end subroutine mg_update_refinement
101 
102  !> Copy a variable to the multigrid tree, including a layer of ghost cells
103  subroutine mg_copy_to_tree(iw_from, iw_to, restrict, restrict_gc, factor, state_from)
105  use mod_forest
106  integer, intent(in) :: iw_from !< Variable to use as right-hand side
107  integer, intent(in) :: iw_to !< Copy to this variable
108  logical, intent(in), optional :: restrict !< Restrict variable on multigrid tree
109  logical, intent(in), optional :: restrict_gc !< Fill ghost cells after restrict
110  real(dp), intent(in), optional :: factor !< out = factor * in
111  !> If present, use this as the input state
112  type(state), intent(in), optional, target :: state_from(max_blocks)
113  integer :: iigrid, igrid, id
114  integer :: nc, lvl, ilo^D, ihi^D
115  type(tree_node), pointer :: pnode
116  real(dp) :: fac
117  logical :: do_restrict, do_gc
118 
119  if (.not. mg%is_allocated) &
120  error stop "mg_copy_to_tree: tree not allocated yet"
121 
122  fac = 1.0_dp; if (present(factor)) fac = factor
123  do_restrict = .false.; if (present(restrict)) do_restrict = restrict
124  do_gc = .false.; if (present(restrict_gc)) do_gc = restrict_gc
125 
126  ilo^d=ixmlo^d-1;
127  ihi^d=ixmhi^d+1;
128 
129  do iigrid = 1, igridstail
130  igrid = igrids(iigrid);
131  pnode => igrid_to_node(igrid, mype)%node
132  id = pnode%id
133  lvl = mg%boxes(id)%lvl
134  nc = mg%box_size_lvl(lvl)
135 
136  ! Include one layer of ghost cells on grid leaves
137  if (present(state_from)) then
138  mg%boxes(id)%cc({0:nc+1}, iw_to) = &
139  fac * state_from(igrid)%w({ilo^d:ihi^d}, iw_from)
140  else
141  mg%boxes(id)%cc({0:nc+1}, iw_to) = fac * ps(igrid)%w({ilo^d:ihi^d}, iw_from)
142  end if
143  end do
144 
145  if (do_restrict) then
146  call mg_restrict(mg, iw_to)
147  if (do_gc) call mg_fill_ghost_cells(mg, iw_to)
148  end if
149 
150  end subroutine mg_copy_to_tree
151 
152  !> Copy a variable from the multigrid tree
153  subroutine mg_copy_from_tree(iw_from, iw_to, state_to)
155  use mod_forest
156  integer, intent(in) :: iw_from !< Variable to use as right-hand side
157  integer, intent(in) :: iw_to !< Copy to this variable
158  !> If present, use this as the output state
159  type(state), intent(inout), optional, target :: state_to(max_blocks)
160  integer :: iigrid, igrid, id
161  integer :: nc, lvl
162  type(tree_node), pointer :: pnode
163 
164  if (.not. mg%is_allocated) &
165  error stop "mg_copy_from_tree: tree not allocated yet"
166 
167  do iigrid = 1, igridstail
168  igrid = igrids(iigrid);
169  pnode => igrid_to_node(igrid, mype)%node
170  id = pnode%id
171  lvl = mg%boxes(id)%lvl
172  nc = mg%box_size_lvl(lvl)
173 
174  if (present(state_to)) then
175  state_to(igrid)%w(ixm^t, iw_to) = mg%boxes(id)%cc({1:nc}, iw_from)
176  else
177  ps(igrid)%w(ixm^t, iw_to) = mg%boxes(id)%cc({1:nc}, iw_from)
178  end if
179  end do
180  end subroutine mg_copy_from_tree
181 
182  !> Copy from multigrid tree with one layer of ghost cells. Corner ghost cells
183  !> are not used/set.
184  subroutine mg_copy_from_tree_gc(iw_from, iw_to, state_to)
186  use mod_forest
187  integer, intent(in) :: iw_from !< Variable to use as right-hand side
188  integer, intent(in) :: iw_to !< Copy to this variable
189  !> If present, use this as the output state
190  type(state), intent(inout), optional, target :: state_to(max_blocks)
191  integer :: iigrid, igrid, id
192  integer :: nc, lvl, ilo^D, ihi^D
193  type(tree_node), pointer :: pnode
194 
195  if (.not. mg%is_allocated) &
196  error stop "mg_copy_from_tree_gc: tree not allocated yet"
197 
198  ilo^d=ixmlo^d-1;
199  ihi^d=ixmhi^d+1;
200 
201  do iigrid = 1, igridstail
202  igrid = igrids(iigrid);
203  pnode => igrid_to_node(igrid, mype)%node
204  id = pnode%id
205  lvl = mg%boxes(id)%lvl
206  nc = mg%box_size_lvl(lvl)
207 
208  if (present(state_to)) then
209  state_to(igrid)%w({ilo^d:ihi^d}, iw_to) = mg%boxes(id)%cc({0:nc+1}, iw_from)
210  else
211  ps(igrid)%w({ilo^d:ihi^d}, iw_to) = mg%boxes(id)%cc({0:nc+1}, iw_from)
212  end if
213  end do
214  end subroutine mg_copy_from_tree_gc
215 
216  !> Generate a multigrid tree that includes the amrvac tree, but also contains
217  !> coarser grid levels. A number of checks has already been performed in
218  !> mg_setup_multigrid, so we don't repeat these checks here.
219  subroutine mg_tree_from_amrvac(mg)
220  use mod_forest
222  type(mg_t), intent(inout) :: mg
223  integer :: i, n, id, ix(ndim)
224  integer :: n_boxes_total, i_c, c_id, c_ix(ndim)
225  integer :: min_lvl, lvl
226  integer :: nb, nb_ix, nb_dim
227  integer :: n_finer
228  type(tree_node), pointer :: pnode, pnode_ch
229  type(tree_node_ptr), allocatable :: id_to_node(:)
230  real(dp) :: dr_coarse
231 
232  ! Estimate number of finer blocks
233  n_finer = nparents+nleafs
234 
235  call mg_build_rectangle(mg, [ domain_nx^d ], block_nx1, dx(:,1), &
236  [ xprobmin^d ], periodb, n_finer)
237 
238  mg%highest_lvl = levmax
239  n_boxes_total = mg%n_boxes + n_finer
240 
241  ! To link the two trees
242  allocate(id_to_node(n_boxes_total))
243 
244  ! Link base level
245  do i = 1, size(mg%lvls(1)%ids)
246  id = mg%lvls(1)%ids(i)
247  ix = mg%boxes(id)%ix
248 
249  pnode => tree_root({ix(^d)})%node
250  pnode%id = id
251  id_to_node(id)%node => pnode
252  mg%boxes(id)%rank = pnode%ipe
253  end do
254 
255  ! Add refinement
256  do lvl = 1, mg%highest_lvl
257  do i = 1, size(mg%lvls(lvl)%ids)
258  id = mg%lvls(lvl)%ids(i)
259  pnode => id_to_node(id)%node
260 
261  if (.not. pnode%leaf) then
262  call mg_add_children(mg, id)
263 
264  do i_c = 1, mg_num_children
265  c_id = mg%boxes(id)%children(i_c)
266  c_ix = mg_child_dix(:, i_c) + 1
267  pnode_ch => pnode%child({c_ix(^d)})%node
268  id_to_node(c_id)%node => pnode_ch
269  pnode_ch%id = c_id
270  mg%boxes(c_id)%rank = pnode_ch%ipe
271  end do
272  end if
273  end do
274 
275  call mg_set_leaves_parents(mg%boxes, mg%lvls(lvl))
276 
277  if (lvl < mg%highest_lvl) then
278  call mg_set_next_level_ids(mg, lvl)
279  call mg_set_neighbors_lvl(mg, lvl+1)
280  end if
281  end do
282 
283  ! Store boxes with refinement boundaries (from the coarse side)
284  do lvl = 1, mg%highest_lvl
285  call mg_set_refinement_boundaries(mg%boxes, mg%lvls(lvl))
286  end do
287 
288  ! Assign boxes to MPI processes
289  call mg_load_balance_parents(mg)
290 
291  ! Allocate storage for boxes owned by this process
292  call mg_allocate_storage(mg)
293 
294  if (associated(mg_after_new_tree)) then
295  call mg_after_new_tree()
296  end if
297 
298  end subroutine mg_tree_from_amrvac
299 
300 end module mod_multigrid_coupling
301 
integer, parameter, public mg_iphi
Index of solution.
integer, parameter, public mg_bc_neumann
Value to indicate a Neumann boundary condition.
integer, parameter, public mg_cylindrical
Cylindrical coordinate system.
integer, parameter, public mg_num_children
subroutine, public mg_set_methods(mg)
subroutine, public mg_fill_ghost_cells(mg, iv)
Fill ghost cells at all grid levels.
subroutine, public mg_build_rectangle(mg, domain_size, box_size, dx, r_min, periodic, n_finer)
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 mg_restrict(mg, iv)
Restrict all levels.
subroutine, public mg_add_children(mg, id)
subroutine, public mg_set_refinement_boundaries(boxes, level)
Create a list of refinement boundaries (from the coarse side)
subroutine, public mg_set_next_level_ids(mg, lvl)
subroutine, public mg_set_neighbors_lvl(mg, lvl)
integer, parameter, public mg_bc_dirichlet
Value to indicate a Dirichlet boundary condition.
subroutine, public mg_allocate_storage(mg)
Allocate communication buffers and local boxes for a tree that has already been created.
integer, parameter, public mg_ndim
Problem dimension.
subroutine, public mg_deallocate_storage(mg)
Deallocate all allocatable arrays.
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_dix
integer, parameter, public mg_bc_continuous
Value to indicate a continuous boundary condition.
Module with basic grid data structures.
Definition: mod_forest.t:2
integer, save nleafs
Number of leaf block.
Definition: mod_forest.t:76
type(tree_node_ptr), dimension(:^d &), allocatable, save tree_root
Pointers to the coarse grid.
Definition: mod_forest.t:29
integer, save nparents
Number of parent blocks.
Definition: mod_forest.t:73
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
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:7
integer, parameter cartesian
Definition: mod_geometry.t:8
integer, parameter cylindrical
Definition: mod_geometry.t:10
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
integer, parameter bc_asymm
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer mype
The rank of the current MPI task.
integer block_nx
number of cells for each dimension in grid block excluding ghostcells
integer, dimension(:), allocatable, parameter d
integer ixm
the mesh range of a physical block without ghost cells
integer, parameter bc_periodic
logical, dimension(ndim) periodb
True for dimensions with periodic boundaries.
integer, parameter bc_cont
double precision, dimension(:,:), allocatable dx
integer, parameter bc_symm
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
subroutine mg_copy_boundary_conditions(mg, iw)
Set multigrid boundary conditions for the solution according to variable iw.
subroutine mg_update_refinement(n_coarsen, n_refine)
If the grid has changed, rebuild the full multigrid tree.
subroutine mg_copy_to_tree(iw_from, iw_to, restrict, restrict_gc, factor, state_from)
Copy a variable to the multigrid tree, including a layer of ghost cells.
procedure(after_new_tree), pointer mg_after_new_tree
If defined, this routine is called after a new multigrid tree is constructed.
subroutine mg_copy_from_tree_gc(iw_from, iw_to, state_to)
Copy from multigrid tree with one layer of ghost cells. Corner ghost cells are not used/set.
subroutine mg_tree_from_amrvac(mg)
Generate a multigrid tree that includes the amrvac tree, but also contains coarser grid levels....
subroutine mg_setup_multigrid()
Setup multigrid for usage.
subroutine mg_copy_from_tree(iw_from, iw_to, state_to)
Copy a variable from the multigrid tree.
Pointer to a tree_node.
Definition: mod_forest.t:6
The data structure that contains information about a tree node/grid block.
Definition: mod_forest.t:11