Geometric Multigrid

[lroberts36]

Solution and residual of Poisson equation after fifteen iterations of MG-BiCGStab with a step function, circular right hand side centered at zero, $D=1000$ in the L-shaped region and $D=1$ elsewhere.

PR Summary

This PR adds the infrastructure required for implementing geometric multi-grid (GMG) solvers in Parthenon. Additionally, it adds an example that solves the Poisson equation $\nabla \cdot D \nabla u = f$ (possibly with rapidly spatially varying $D$) using both raw FAS multigrid v-cycles with FAC like grids and using this same multigrid solver as a pre-conditioner for the Krylov subspace method BiCGStab. Good convergence is seen on both uniform and non-uniform grids (including flux correction) and the MG-BiCGStab method provides rapid convergence on problems with rapidly varying $D$ (which cannot be efficiently solved by either BiCGStab or MG V-cylces alone). The implementation of a GMG solver task list and its use as a preconditioner is shown in the poisson_gmg example, but I have not included generic solver routines within Parthenon itself. Generic solvers probably could be provided, since a user probably only needs to provide a method for applying there matrix to a vector and a method for smoothing (e.g. for doing Jacobi iterations).

In terms of infrastructure, multigrid generally requires a hierarchy of increasingly coarsened meshes, a way to restrict a field from a finer mesh to the next coarser mesh, a way to prolongate data from a coarser mesh to a finer mesh, and a mechanism for intra-level communication between blocks at each GMG level.

To do this in Parthenon, a hierarchy of de-refined block meshes are created from the base mesh during initialization and after each AMR step¹. The blocks in these coarsened meshes are assigned a global id (gid) by iterating through them in Morton order, but starting from the largest gid of the next finest mesh. This gives each block in the GMG hierarchy a unique identifier. Once these de-refined block lists are built, a same level neighbor block finding algorithm is performed on each level², the parent block on the next grid up in the GMG hierarchy is found, and all daughter blocks on the next grid down in the GMG hierarchy are found. All of these relationships are stored as arrays of NeighborBlocks in MeshBlock. The coarsened levels can go above the root grid of the simulation and the blocks change in size at these coarser levels where required to cover only the root grid. At least for cubic grids and other grids that have a number of root blocks and zones that are a large enough power of two, the coarsening will continue to negative grid levels and possibly to a block with a single zone.

With these GMG levels, neighbor relationships and unique ids in hand, the pre-existing same level boundary communication routines can automatically handle inter-level communication on each level of the GMG hierarchy. With a small addition to ForEachBoundary template³ and a different choice of prolongation and restriction regions, the pre-existing communication infrastructure can also be repurposed for prolongation and restriction from one GMG grid to another⁴. GMG prolongation and restriction only occurs for a variable if it has Metadata::GMGProlongate or Metadat::GMGRestrict set, respectively.

There is a zoo of different ways to implement GMG (different smoothers, different types of cycles, using it as a preconditioner, etc.), and I have only tested a reasonably small subset. Getting GMG to work on a uniform grid was fairly straightforward, although it did require being careful about how the boundary conditions are defined (mainly that the boundary condition was defined at the surface of the simulation domain and not just imposed at the cell center of the first ghost zone, since that position moves as you go to coarser levels). I eventually only focused on cell-centered fields but all of the infrastructure should work for other topological types, they would just required different ways of dealing with physical boundary conditions. Here, I also experimented with MG preconditioned CG and found that it worked fairly well. For refined meshes, I abandoned CG since it is very difficult to formulate partially refined problems as a symmetric matrix (I think this would rely on some very specialized prolongation and restriction operators).

For refined meshes in my implementation, each GMG level only includes blocks that are at the current level, which implies that at fine levels the blocks may not cover the entire mesh. Dealing with accurately setting the boundary values of these blocks is a little tricky, and I found that convergence is poor or non-exististent if the fine-coarse boundaries of a given level are not self-consistently updated (since the boundary prolongation from the coarse grid to the fine grid also depends on interior values of the fine grid that are being updated by a smoothing operation). This means that each smoothing step, boundary communication must occur between all blocks corresponding to all internal and leaf nodes at a given level in the tree and with all leaf nodes at the next coarser level which abut blocks at the current level. Therefore, the GMG block lists contain blocks at two levels, but smoothing operations etc. only occur on the subset of those blocks that are at the fine level. This is effected by using the new functionality in SparsePack that allows choosing a subset of blocks in a MeshData to include in a SparsePack. I think this general structure is how most codes implement the refined grids, with blocks only on a single level, but the Athena++ implementation differs I think. It may be that in the future we determine there is some better choice of GMG grids, which will require changes to the underlying infrastructure as written.

This PR will close #858.

Description of changes:

• NeighborBlock: Added NeighborBlock constructor (to be used in place of NeighborBlock::SetNeighbor, although the old method still gets called some places in the code). Also added a RegionSize member variable to hold information about the size of the neighbor block, since in the GMG grid hierarchy all blocks are not the same size. NeighborBlock instances can refer to intra-level neighbors (i.e. MeshBlock::neighbors, MeshBlock::gmg_same_neighbors, and MeshBlock::gmg_composite_finer_neighbors) or to refer inter-level neighbors (i.e. parents or daughters on other GMG composite grids, MeshBlock::gmg_coarser_neighbors or MeshBlock::gmg_finer_neighbors.
• CalcIndices (which is used by BndInfo and ProResInfo for calculating buffer packing indices): Refactored a little to clarify the different types of index ranges it can provide, which is determined by the new enumerator IndexRangeType which is passed to CalcIndices. Additionally, CalcIndices is extended to work for prolongation and restriction of the interiors of blocks and to be able to handle when inter-grid transfers involve different block sizes. For inter-grid prolongation, the index range of the entire coarse buffer of the receiving zone is filled so that prolongation on the receiver can occur without a communication step. For this to work correctly, the ghosts of the coarse sending block must be filled before the inter-grid communication is called.
• ProResInfo::Get* and BndInfo:Get*: Minor refactoring to clarifly index range types.
• Boundary Communication: Five new boundary types been added to BoundaryType: gmg_same which encapsulates intra-grid communication on a GMG composite grid, gmg_restrict_send and gmg_restrict_recv which are responsible for inter-GMG grid restriction, and gmg_prolongate_send and gmg_prolongate_recv which are responsible or inter-GMG grid prolongation. The ForEachBoundary template loop has been updated to loop over these intergrid boundaries correctly, so communication routines for these new boundary types can be called, for example SendBoundBufs<gmg_restrict_send> can be called on a finer GMG level and ReceiveBoundBufs<gmg_restrict_recv> can be called on MeshData for a coarser level to restrict from fine to coarse. Small changes were required for choosing which version of RebuildBufferCache to call so that the correct prolongation range is used for inter-gmg communication. Small changes to BuildBoundaryBuffers were also required to make it work for the new BoundaryTypes. The TagMap object was also updated to be able to include more than one type of boundary (since we need MPI tags for all of the different types of boundary communicators).
• DataCollection: Added overloads so that MeshData corresponding to the blocks in a particular GMG level can be requested and have the its grid type identified.
• MeshData: Added a member object that stores the type of grid that the blocks in a MeshData are associated with.
• Metadata: Added flags GMGRestrict and GMGProlongate to denote fields which should be prolongated or restricted during GMG inter-level communication.
• SparsePack: Added a method to get the IndexShape of a variable on device from a ParArray3D, since now it is possible that the size of blocks in a MeshData object, and hence a SparsePack, may not all be the same size. This feature doesn't really end up being used though currently (it was used during testing).
• LogicalLocation: Extended to deal with "negative" levels. A lot of functionality moved to a cpp file instead of a header. A couple of bugs relating to periodic grids were fixed (although I honestly don't remember the details). There may be more efficient ways to perform some of the tasks done by LogicalLocation. If there are performance issues related to LogicalLocation, there is likely room to speed things up.
• Mesh::SetSameLevelNeighbors: Provides the same functionality as SearchAndSetNeighbors but using maps of LogicalLocation instead of the old MeshBlockTree object (I have tested fairly extensively that it gives the same results as the old code). It also extends SearchAndSetNeighbors to deal with two-level GMG composite grids where only communication into the boundaries of blocks on the finer of the two levels is required. Currently, the code finds neighbors using both SearchAndSetNeighbors and SetSameLevelNeighbors for the main grid. SearchAndSetNeighbors sets boundary information contained in MeshBlock::pbval while SetSameLevelNeighbors sets the neighbors contained in MeshBlock::neighbors, the latter is what is used now by all of the boundary communication routines for fields but I think swarms rely on pbval still. Nevertheless, things are pretty close to the point that we can remove pbval and a lot of the code in bvals_base.cpp.
• Mesh::BuildGMGHierarchy, Mesh::gmg_grid_locs, and Mesh::gmg_block_lists: Build and store the logical locations and block lists associated with each GMG composite-grid. Mesh::BuildGMGHierarchy also finds and sets the inter-grid neighbors for each block.
• Mesh::multigrid: Added flag for turning on and off multigrid parthenon/mesh/multigrid=true/false. If false, none of the internal nodes of the tree will be populated with blocks and there should be no change in how the code behaves from previous versions (although there was refactoring of internals, e.g. in Mesh and using SetSameLevelNeighbors).
• MeshGenerator: Removed unused code.
• Refactored code for doing static refinement in Mesh constructor.
• Updated Mesh::SetBlockSizeAndBoundaries and Mesh::GetBlockSize to deal with negative levels where the block size changes.
• MeshBlock: Updated index shape to correctly include ghosts for non-symmetry directions when the block size has gone to one (i.e. the coarsest GMG grid). Added vectors of neighbors that are used in communication (see ForEachBoundary).
• Added a linear prolongation operation since the default prolongation operation is non-linear and therefore (at least formally) can't be used for variables used in a linear solver.
• Created a new utilities header to store bit hacks, moved some there and added one for finding the number of binary trailing zeros (which is useful for determining how coarse a particular grid can be made).

PR Checklist

• Code passes cpplint
• New features are documented.
• Adds tests for new features.
• Code is formatted
• Changes are summarized in CHANGELOG.md
• CI has been triggered on Darwin for performance regression tests.
• Docs build
• (@lanl.gov employees) Update copyright on changed files

Footnotes

1. Rather than storing these GMG meshes in their own octree like the base mesh, they are stored in a std::map from LogicalLocation to a (gid, rank) pair in the field Mesh::gmg_grid_locs. The different level block lists are stored in Mesh::gmg_block_lists and then MeshData associated with each level of the GMG hierarchy are stored in Mesh::gmg_mesh_data. We may want to think about if we want to have a separate data collection for each grid, or store them all in the already existing Mesh::mesh_data. The map is used because it is default ordered by Morton number, so it is easy to iterate through and assign gids. That being said, it may be more performant to switch to a std::unordered_map and just put the logical locations into a vector when gids need to be assigned, sort the vector, then assign gids in the std::unordered_map. ↩

2. Neighbor finding is currently provided by BoundaryBase::SearchAndSetNeighbros, but a refactored version of the algorithm using the new LogicalLocation infrastructure (SetSameLevelNeighbors, which might be a bad name since the level here refers to GMG level). A lot of the trickiness in this PR was updating the LogicalLocation functionality to return things like block offsets when periodic boundary conditions are present. ↩

3. The current ForEachBoundary just iterates over every block in a MeshData object, then every variable present on that block, then every neighbor block of that block (skipping some of these depending on the template parameters provided). Now, if GMG boundary template parameters are passed to ForEachBoundary it does the same set of iterations but replacing the list of same grid neighbor blocks with either the list of parent or daughter blocks in the GMG hierarchy. ↩

4. Communication buffers are stored in a map with a key made up of a sending gid, a receiving gid, and an array encoding the relative offset of the blocks. This naturally extends to storing multiple grids, but does mean that new communication buffers are created for each level. In particular, the GMG prolongation and restriction buffers are the same size as the coarse buffer for each block, which feels a little redundant and maybe has an impact on the memory footprint of simulations. ↩

[pgrete]

Thanks for addressing all the comments. I'm happy to see this merged (and documented in the other PR ;) )

[lroberts36]

Thanks @pgrete, I just enabled auto-merge. Documentation definitely needs to be added about logical location, although I thought there was a good deal of documentation on this PR relative to my baseline ;). It does look like check PR dependencies is failing now, but maybe other PRs are seeing this too? (@bprather, @Yurlungur)

[Yurlungur]

It does look like check PR dependencies is failing now, but maybe other PRs are seeing this too? (@bprather, @Yurlungur)

I think the check dependencies machinery is broken

[lroberts36]

@Yurlungur, ok thanks, I will just force merge once everything else passes.

[bprather]

The dependencies check has been fixed. I re-ran it and it passes, but it looks like some others are still failing.

[bprather]

Looks like the auto-merge didn't take so I merged develop and tried it again. Speak now etc etc

EDIT Oops saw the note on the fix. Thanks @lroberts36

[lroberts36]

I took off the auto-merge because #966 needs to be merged into this first.

[lroberts36]

Thanks @bprather, I had forgotten to finish this up. I just merged in the bugfix, I hadn't realized it was approved. I am turning back on auto-merge.