[petsc-dev] post 3.1 reorganization of PETSc DMMG code

[Richard Katz]

Hi Barry,

I think this would be very useful.  In all my magma dynamics codes, there is a block of equations that is composed of Stokes, the Compaction equation (good Helmholtz), and Continuity.  Marc has shown that multigrid works very well on the Compaction equation.  It would be interesting to experiment with various combinations of preconditioners, in particular ILU + multigrid.

It also makes sense to me to merge DMMG and SNES.  I usually write my codes with the DMMG interface, but have also tried SNES, and found it clearer (probably because i don't currently use multigrid).  

I should also note that for teaching purposes (and because I am lazy), I have students solve Laplace, Poisson, etc with the DMMG interface, because I want them to learn about Newton's method (all the research problems they work on are non-linear).  So it is good to have two different ways to solve linear systems: as a LINEAR system (KSP), and as a linear SYSTEM (DMMG/SNES).

As to how to "get subphysics" I don't really have any suggestion besides keeping it user-friendly, by which I mean user-like-me-friendly, of course. :)  The idea of having special-purpose, pre-defined ISs like IS_FIELD_1, etc seems good to me. This sort of functionality would be helpful in improving the semi-Lagrangian advection implementation.

With ongoing, profound gratitude to the developers of PETSc,
Rich




On 19 Mar 2010, at 21:23, Barry Smith wrote:

> 
>  Now that the 3.1 repository has been cloned I'd like to start work on the reorganization of PETSc to pull all the functionality of DMMG into the SNES, PC and TS components. I want to do this because
> 
> 1) the DMMG "solver" are not compossible with each other and other solvers. For example, inside a PCFIELDSPLIT preconditioner I'd like to use a multigrid solver on one component and I'd like the management of that multigrid solver to be handled by something as simple as the DMMG.
> 
> 2) the DMMG solvers cannot be used by TS
> 
> 3) there is currently two ways to solve linear systems, with KSPSolve() directly or with DMMGSolve(). One of my most fundamental design principles is having two DIFFERENT ways of doing something is bad
> 
> 4) DMMG very nicely handles a hierarchy grids and sets everything in the PCMG structure automatically for you, but it does not handle, for example, setting up a PCFIELDSPLIT preconditioner
> 
> Thus I propose we completely eliminate the DMMG code and instead have each TS, SNES, PC object (for KSP objects this is pretty much irrelevant) be able to take a DM and "fill themselves up" by using appropriate information from the DM; just as the DMMG "fills up the PCMG" with the correct interpolation/restriction, coarser grid matrices etc. Thus much of the current DMMG code will end up inside the PCMG.  Similarly PCFIELDSPLIT will be able to take its DM and fill up much of itself.
> 
> To accomplish this we will need the DM to be a richer object. Currently a DM can give coarsened versions of itself (and refined versions) we need to add a way to have it give "subphysics" versions of itself. For example if we are solving a 2-component PDE and have an associated DM, then we would like to be able to get a new DM associated with a single component. For example, say we wish to use PCFIELDSPLIT with multigrid on the first component and ILU on the second. The PCFIELDSPLIT would have its two PCs and pass down to them the sub DM. Now the multigrid subPC would use the DM it is given to coarsen and create the hierarchy of grids and interpolations etc. Currently when Jed codes this he has to manage that himself or do some ugly crud where he wraps a DMMG inside that fieldsplit. If we do this right it could be extremely powerful.
> 
> Tasks
> 
> (1) merge DMMG linear into PCMG. I think this is straightforward and doesn't have any conceptual challenges at all
> (2) merge DMMG nonlinear in SNES. With (1) done this is pretty much done except for the conceptual issue of how to handle grid sequencing?
> 
> This solves the multigrid problem.
> 
> (3) introduce into DM the ability to "get subphysics" DM. Here the conceptual question is how to indicate them. Do we introduce the idea of having each DM as having N "fields" and then to be able to request "field 1 and 2" etc? PCFIELDSPLIT has this concept with PCFieldSplitSetFields(). This is easy enough to implement but I don't like it. The reason I don't like it is that we already have the concept of IS that is used to indicate subsets of things; I HATE to have two different ways of doing something. We also have PCFieldSplitSetIS() which is "generalization" of PCFieldSplitSetFields(), if I wasn't so lazy when I added PCFieldSplitSetIS() I would have immediately removed PCFieldSplitSetFields(). So I would like to have DMGetSubDM(DM,MPI_Comm,IS,DM*). For DA with suitable stride and block stride IS this is easily done and makes sense, an arbitrary IS would, of course, have to generate an error.  For DMCOMPOSITE providing sub IS's would be cumbersome and not intuitive at all. So perhaps we do need both types of getters. (Side note, if we introduce some builtins like IS_FIELD_1, IS_FIELD_2 then perhaps we could have one interface that was reasonably painless for the simple cases of fields but is general? I'll think about this.)
> 
>  Suggestions, comments improvements?
> 
>   Barry
> 
> 
> 
> 

___________________________________
Richard Foa Katz
RCUK Academic Fellow,   Univ Oxford
http://www.earth.ox.ac.uk/~richardk