[Discussion]: Geometric Multigrid via PETSc DM

[pbauman]

I'm opening this up as a catch-all thread for the design discussion for this interface. This is work myself, @roystgnr, and @bboutkov have been pecking away at for awhile over in GRINS and we're starting to migrate what we've done into libMesh.

The idea is to interact with the DM infrastructure in PETSc. By setting up this interaction, we'll interact with PETSc from the command line to setup/use solvers in the same way Matt Knepley does with his Plex examples (unstructured grids) or PETSc DMDA (structured grids). In particular, using Geometric Multigrid (GMG) as a preconditioner or solver and potentially only on particular solution blocks. The plan is to use the mesh hierarchy from libMesh and setup the needed data structures in PETSc and then we can just do this from the command line.

I'll reference this issue (and the discussion below) in PRs. I'll try and briefly describe what's needed, where we are, current bottlenecks, etc.

What's Needed: Step 1. Setup PetscSection and PetscSF interaction. This is essentially DoF info and DoF partitioning. These alone will facilitate recursive fieldsplits at the command line (my understanding is that the current interface we use does not support the recursive fieldsplit). This will be step one of our migration.

Step 2. Setup interpolation (and possibly restriction) matrices between each of the meshes in the hierarchy and then hooks into PETSc to supply these. The first part of this (creating actual matrices from GenericProjector) was migrated by @roystgnr, based on @bboutkov's work, in #1506. We still need to provide those matrices (and function hooks) into the PETSc DM infrastructure. This will facilitate "Galerkin" multigrid.

Step 3. Setup interface to facilitate nonlinear multigrid ("Full Approximation Scheme", FAS). This will require actually setting functions for assembling residual (and possibly Jacobian) on each mesh level. We haven't really experimented this with yet.

Status: Step 1. I'm going to open a PR for this shortly.

Step 2. This will be the next step of the migration. We've experimented with this and need to clean it up; the clean up will probably be "interactive" within the PR. The support we have so far is explicitly traversing the mesh hierarchy and constructing these matrices on the fly. The code does correctly support DistributedMesh for most cases (periodic boundary conditions are the current sticking point), but the scaling sucks. But we get expected solver behavior on Laplace for 1D, 2D, and 3D, for quads, triangles, hexes, and tets.

Step 3. Haven't investigated in depth yet, but will reuse everything from Step 1 + Step 2 + add separate "injection" for solution and then see how we need to interact with PETSc since Step 2 only caches interpolation matrices and not the mesh state. For FAS, I'm guessing we'll need to cache the mesh state in each stage of the hierarchy, but that's speculation at this point.

Bottlenecks:

1. Periodic constraints and hanging node constraints. My understanding is that we need some additional work in the GenericProjector objects to support this. @roystgnr, @bboutkov can comment more.

2. Scaling: We've been digging and I'm starting to believe it's related/caused by repartitioning (which I believe is the case of #1468, but I'll open a separate issue for that once we have more concrete data and ideas). After talking with @permcody, it looks like MOOSE folks are using the splitter app for large meshes and reading those in after doing the splitting. Should we provide a capability to also do projection matrix construction for mesh hierarchies in a similar vein to the splitter app so that the matrices are read at runtime instead of constructed on the fly to make this more immediately useful? The scaling will still be awful, but is probably less awful than what we see right now.

Design Discussion: We'll need to plop these features in both PetscNonlinearSolver and PetscDiffSolver. At the moment, I'm pursuing a composition (as opposed to inheritance) design. I've created a PetscDMWrapper class that should be able to contain everything and the solver classes need only provide the System and SNES objects. The PetscDMWrapper will be setup all the data structures, including traversal of mesh hierarchy, and attach the necessary bits to the SNES object that is provided by the solver class. Then we need only drop the wrapper object in each of the solvers + appropriate function calls in the appropriate places. I'm working with PetscDiffSolver to begin because that's what I can most easily test, but I'm hoping this will be as easy for PetscNonlinearSolver to use as well. The only catch I can think, beyond discussion of controlling usage at command line and multiple System, will be how FAS usage will play out. Thoughts on this?

Usage Discussion:

1. How to trigger? If we're not doing GMG or fieldsplit, we shouldn't require the user to go through this code path. I was thinking we look for --use_petsc_dm on the command line to trigger this construction. Once we know to follow this code path, then the other relevant PETSc options will be picked up through PETSc and all will be fine. Thoughts?

2. Old style fieldsplit: This will be add features we can't use in the current fieldsplit paradigm, so I think we should prefer this, but in the name of backward compability, we should preserve the old (petsc_auto_fieldsplit). I'm thinking we just go the new codepath if we're using the DM, and the old code path if not. Should deprecate the old code path? The catch is petsc_auto_fieldsplit I believe will allow us to reference variable names in the command line functions, whereas the DM path currently does not, you have to reference by variable index. I think we can prod PETSc folks to fix it upstream (or we can send a patch), but in the near term, the syntax will not be as nice as I'd like.

3. I haven't thought about the multiple System case at all. But it looks like PetscNonlinearSolver only understands one System anyway. Could someone from MOOSE provide a little insight how they would like the multiple System case to play out? Would it be enough to prefix all the options with the system name (which I believe is what is done now)? I believe we could do this with all the GMG and Fieldsplit options.

Other points of discussion and questions are also welcome here.

CC'ing @fdkong since he's expressed interest in the past.

[lindsayad]

After #2047, what remains for this ticket?

[roystgnr]

I pointed out one failing to @lindsayad by phone just now: until I can figure out the right way to fix #1853 and thereby #1832, GMG isn't going to be reliable on adaptively-refined meshes.

Two other concerns for which I'm not sure of the answers, though; hopefully @bboutkov can tell me if I'm right or wrong:

1. How hard will it be to get this working with our other PETSc solver interfaces? It looks to me like adding PetscNonlinearSolver support at this point will be a 10-line patch: just add a PetscDMWrapper to the class and call init_and_attach_petscdm() at the appropriate spot. And it looks like adding PetscLinearSolver support would be nearly as simple, except that first we'll have to refactor init_and_attach_petscdm to allow another overload which uses KSPSetDM instead of SNESSetDM. Am I missing anything?

2. What restrictions do we have on snes_type options? It looks to me like the interface you wrote ought to be enough for PETSc to do multigrid at the nonlinear level, not just after linearization within a Newton-Krylov solve; is that right?

[bboutkov]

I think this list is fairly complete, but there's some additional quips/items I'll mention.

until I can figure out the right way to fix #1853 and thereby #1832, GMG isn't going to be reliable on adaptively-refined meshes.

This still stands AFAIK.

1.  It looks to me like adding PetscNonlinearSolver support at this point will be a 10-line patch: just add a PetscDMWrapper to the class and call `init_and_attach_petscdm()` at the appropriate spot. 

Also pretty much what I expected. I just wasn't aware of an easy place for me to test this as I was mainly using GRINS and was under the impression that PetscNonlinearSolver is what MOOSE uses so wanted to make sure this was working properly and didn't want to guess at the implementation. I think you would need to also properly hook in to call PetscDMWrapper::clear() when appropriate.

And it looks like adding PetscLinearSolver support would be nearly as simple, except that first we'll have to refactor init_and_attach_petscdm to allow another overload which uses KSPSetDM instead of SNESSetDM. Am I missing anything?

Seems right to me although I only really worked through the SNES interface, but I expect that such a change should just work.

2. What restrictions do we have on snes_type options?  It looks to me like the interface you wrote ought to be enough for PETSc to do multigrid at the nonlinear level, not just after linearization within a Newton-Krylov solve; is that right?

Mostly. I did some preliminary testing with snes_type while experimenting with some multiphysics problems, I think some will work, but for example snes_type=fas needs further changes. I think that at the least in order to get FAS to work we would need to implement DMShellSetCreateInjection which handles state transfer between levels which would need to fit into the PetscDMWrapper while initializing the infrastructure. This in turn might require us to store the entire mesh at each multigrid level which could be nontrivial.

3. Periodic BC support is still lacking. Our Projection matrices for GMG are sparse, while DofMap::constrain_element_matrix only has a dense variant and since we make calls to system.reinit_constraints() during init_and_attach_petscdm this will cause problems.

Think that's it for the major items as far as I recall.