A starting point for a polymorphic solver-based TDycore interface.

[jeff-cohere]

This small PR lays part of the foundation for a move to a more solver-based interface. Our objective here, in the spirit of #197, is to move from our current approach of determining the behavior of the dycore (if tests in several functions scattered throughout the library, to match the dycore against its type of discretization, etc) to a new approach based on function pointers. The closest handy analogy for most is C++'s implementation of polymorphism using virtual tables. But PETSc itself uses this approach heavily.

In include/tdycore.h, line 162, you can see that I've written that certain functions using the old approach will likely be replaced by a smaller set of new polymorphic functions. Specifically, we are adding two new functions:

• TDySetSNES(TDy dycore, SNES solver) - takes an SNES nonlinear solver object and sets it up to solve the system relevant to the dycore.
• TDySetTS(TDy dycore, TS solver) - takes a TS ODE solver object and sets it up to solve the system relevant to the dycore.

These functions are simple wrappers and are implemented in tdycore.c. You can see all they do is pass the call along to the appropriate function pointer, or complain with an error if no such function pointer exists. This allows any given implementation of the dycore to simply supply the right function pointers instead of being added to a long chain of if tests in several functions throughout the code base.

Additionally (if people like the direction of this PR), we'll be retrofitting other functions such as TDyComputeErrorNorms to be polymorphic. Any function that contains if tests that try to figure out "which algorithm we're using in the dycore" will be replaced with an interface function that calls the appropriate implementation-specific function pointer.

The idea here is that we can use the solver interface to set specific solver parameters, and to step forward in time. We also talked about a higher level TDycore interface that does some of this stuff in solver-neutral language. That's very easy to add on top of this kind of interface.

I know the PETSc folks are very familiar with this approach--I expect them mostly to correct my "spelling"of polymorphic C code to their tastes. For others, I'm happy to chat about function pointers, polymorphism in general, and/or how this will look in specific cases.

[codecov-commenter]

Codecov Report

Merging #211 (83e1177) into master (adc0c5e) will not change coverage. The diff coverage is n/a.

:exclamation: Current head 83e1177 differs from pull request most recent head 90ae49d. Consider uploading reports for the commit 90ae49d to get more accurate results

@@           Coverage Diff           @@
##           master     #211   +/-   ##
=======================================
  Coverage   49.19%   49.19%           
=======================================
  Files           4        4           
  Lines         683      683           
=======================================
  Hits          336      336           
  Misses        347      347           

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update adc0c5e...90ae49d. Read the comment docs.

[jedbrown]

Thanks for taking this up. One thought on TDycore being used as a library: all the SNESSetX and TSSetX interfaces related to a model actually call corresponding DMSNESSetX and DMTSSetX. The primary thing TDycore could do is to wire up a DM with all the model stuff needed, and drivers could TSSetDM() and solve as they see fit (e.g., perhaps within optimization, etc.). So TDycore doesn't need to take responsibility for operating the solver, it can just set up the DM with the information the solver will use.

[jeff-cohere]

Interesting. Thanks for pointing this out, Jed. If everything happens through the DM, what you suggest does seem like the way to go. I'll try to figure out a simple interface to do things this way.

[jedbrown]

We'll be doing something like this for a libCEED+PETSc mechanics solver (tentatively https://gitlab.com/micromorph/ratel), which is currently just a fork of the mini-app that lives in the libCEED repo, but is intended to evolve into a full-featured library+app.

[jeff-cohere]

To make sure I understand: can Jacobian and residual/RHS functions be supplied to the solvers through the DM? One of the benefits of the TDySetSNES and TDySetTS approach above is that TDycore can set all this stuff on the solver and then the solver can be tweaked and used from there. Most of the PETSc examples I see surrounding e.g. SNESSetDM still involve some problem-specific solver setup.

In other words, the interface I proposed looks like this:

Approach 1

TDy dycore;
TDyCreate(comm, &dycore); // or whatever
TDySetFromOptions(dycore);
// set up the dycore for a specific problem, to the degree not specified by options
...
// now create a solver
SNES solver;
SNESCreate(comm, &solver);
SNESSetFromOptions(solver);
// set up the solver to do all the TDycore things
TDySetSNES(dycore, solver);

// forget about the dycore for a bit and use the solver itself
while (! finished) {
  SNESSolve(...);
}

// refer back to the dycore as needed for postprocessing...

So we're both moving toward a workflow that sets up a solver using the dycore, and then uses that solver. It seems to me that you'd still need to set Jacobian/residual functions for the solver after setting the DM. That's fine--it just means we have to make polymorphic SNES/TS Jacobian/residual/RHS functions that are passed to the solver:

Approach 2 (Corrected to reflect Matt and Jed's comments)

TDy dycore;
TDyCreate(comm, &dycore); // or whatever
TDySetFromOptions(dycore);
// set up the dycore for a specific problem, to the degree not specified by options
// (within the library, the dycore configures a DM for its problem of interest)
...
// now create a solver
SNES solver;
SNESCreate(comm, &solver);
SNESSetFromOptions(solver);
// extract a DM from the dycore and use it to set up the solver
DM dm;
TDyGetDM(dycore, &dm);
SNESSetDM(solver, dm);

// forget about the dycore for a bit and use the solver itself
while (! finished) {
  SNESSolve(...);
}

// refer back to the dycore as needed for postprocessing...

These two approaches are very similar. The first one is more of a "black box" that sets the DM/function/Jacobian on the solver in one command ("TDycore, just please do whatever you need to the solver"), and the second one reflects how it's done more explicitly in a PETSc app.

[knepley]

You would replace the SNESSetFunction() with something like DMSNESSetFunctionLocal(), and similarly for the Jacobian. I also put in a callback for imposing Dirichlet conditions.

[jeff-cohere]

Maybe TDySetUpSNES is a better name than TDySetSNES in Approach 1, since TDySetSNES might imply that the SNES belongs to the dycore, which it doesn't.

[bishtgautam]

Unlike PETSc, I was not envisioning a demo to specify a user-defined residual/jacobian/ifunction/ijacobian etc. Those functions will be internal to the TDycore lib and get automatically set depending on the choice of the mode (=Richards/TH/etc) + Discretization (FE/FV/etc). I'm most likely missing the point in the above discussion of why a demo would be setting user-defined functions. Can someone explain the reasoning behind it?

[jeff-cohere]

Hey @bishtgautam ,

Approach 1 is probably a bit more friendly as a "black box" way to set up a solver. This is what I had originally envisioned, though I probably wasn't very clear about the dycore's role in a demo (both approaches use TDycore to set up a solver and then use the solver itself instead of the dycore to run a simulation). Aside from the need to assign these functions, the two approaches are basically the same.

Maybe Jed and Matt have thoughts about this, though.

[jedbrown]

I would go with 2, but DMSNESSetFunctionLocal gets called internally so it's really just SNESSetDM(snes, dm) and off to the solves. If you're solving a transient problem, it's exactly the same except TSSetDM(ts, dm).

How do initial conditions and initial guesses for Newton get created? If they require a solve, what does it look like and what data needs to be read?

[jeff-cohere]

Okay, I get it now. We do the calls to set the functions for J and R internally as part of setting up the DM, and then the dycore turns over the DM to the solver, communicating those functions. The two approaches are identical, except in the case of Approach 2, the dycore never has to talk directly with the solver. Neat. I like it.

Let's go with Approach 2. In this case, it's simplest if the dycore creates the DM itself, since it's going to be configuring these kinds of things on it. I recommend we get rid of TDySetDM to clarify the "life cycle" of the DM, and just use options (and possibly the DM itself, obtained from TDyGetDM to fiddle with it. If we keep TDySetDM, we'll have to carefully document it to make sure the caller knows "you can call this after X and not after Y" and "you might get something different than you imagine in some cases."

Good questions about initial conditions and Newton guesses. I would guess the answers depend on the range of complexity involved. If we're doing "spin up", we're going to need solves there.

[jeff-cohere]

If we end up solving a different set of equations for initial conditions, maybe it's worth thinking about another DM constructed for that purpose, to be passed to a solver that can be fiddled with, same as we're discussing for solving the transient problem?

[jedbrown]

If the equations are different, then we want the same function space (but perhaps different boundary conditions). I wish this were explicit by making the "function space" part of DM distinct from the "equation set" part of DM (a dispute I lost to Barry a decade ago).

If the initial condition problem is formulated with different essential boundary conditions, then we wouldn't use the same DM for the forward model, but we will need an affine projection between the spaces. Maybe if we can get a specific example from the science team, @knepley and I can make sure it's supported and outline a process.

[knepley]

@jedbrown I did this separation in "stealth" mode with PetscWeakForm.

Currently, I do this kind of affine projection through the back door with DMPlexSetBoundaryValues() and Toby has complained about this, so it would be good to have a proper interface for it.

[jeff-cohere]

Thanks, guys. This sounds like a productive direction for model initialization. Maybe @bishtgautam and I (and others interested) can come up with an example that illustrates what we need.

[jeff-cohere]

I'm merging this PR so we can continue working in the direction of polymorphism. I opened #213 to track our model initialization discussion.