Blocked linear system

issue 198

First step: move to Thyra objects

A first step towards the block linear algebra is to switch to Thyra linear algebra structures, without introducing blocked structures yet. This will change/add some interfaces, mostly in AbstractDiscretization (added) and Application/Workset (changed), plus some internal changes in the evaluators (only gather/scatter ones, in theory). Delaying the introduction of block structures to a second step will allow for a quick(er) merge into master, so that new features can be developed with the new interfaces.

Here's a tentative todo list for the first step. If someone completes them, mark them as done.

• (a) Add getters in AbstractDiscretization for Thyra vector spaces. These can simply create the VS on the fly, wrapping the corresponding Tpetra_Map.
• (b) Replace Tpetra/Epetra structures in Workset with Thyra ones. As a half-step, we can simply add Thyra structures,
• (c) Change Application/ModelEvaluatorT, to work directly on Thyra objects, rather than Tpetra/Epetra. In particular, ModelEvaluatorT can simply forward InArgs contents to the Application.
• (d) Change gather/scatter evaluators, so that they can work with a Workset storing only Thyra objects.

Steps (b)-(d) can be performed incrementally on the objects stored in Workset. For instance, we can start from the residual evaluation only. All other structures will be left as Tpetra/Epetra.

Note: ModelEvaluator (epetra version) does not inherit from Thyra's ME, which complicates things. There is a Thyra adapter for EpetraExt ME, which would allow to use Thyra for both Epetra and Tpetra builds, but it's not used in Albany. I see two options: (1) make ModelEvaluator inherit from the Thyra adapter to EpetraExt's ME; (2) remove Epetra capabilities altogether. The first option is the most backward compatible. But the fact that Epetra-only builds are not possible makes me wonder why shouldn't we go for (2)...

Update 8/6/18

The thyra refactor is ongoing. So far, all inputs in the Workset (x,xdot,xdotdot,Vx,Vxdot,Vxdotdot,Vp) have been replaced with Thyra structures. The refactor for outputs (f,jac,g,gx,gxdot,gxdotdot,gp,fp,...) is ongoing. When this is completed (i.e., all structures in the workset that were Tpetra/Epetra are converted to thyra), there will be a cleaning phase, where all conversions thyra<->*petra are pushed down to the stack, to the very place where the knowledge of the backend linear algebra package is absolutely necessary (read/write/load/save/import/export).

At that point, the thyra refactor will be pretty much over, and we can start thinking about blocks. This will include a serious thinking about how to perform assembly: by block (1 sweep per block), by dimension of mesh objects (1 sweep per dimension), a hybrid approach?

Summary of block discretization meeting (8/21/18)

Participants:

• Luca Bertagna
• Mauro Perego
• Irina Tezaur
• Jerry Watkins
• Dan Ibanez
• Eric Cyr

Discussion:

• DofManager: talk to Roger, find out if we can make panzer's adapter-stk independent from disc-fe. If yes, then aim at pulling DofManager into Albany (offer to help).
• assembly: perform assembly by geometric entity. All n-dimensional assembly goes in the same field manager. This means different blocks on same geometry, are gathered, assembled and scattered together. Also, if you have a 3d-2d coupling, have 1 FM for the 3d stuff, and 1 FM for the bc of the 3d problem AND the 'volume' 2d assembly.
• even if using block, support the possibility of assembling monolithic linear algebra structures. Note: it's tempting to say 'if you want monolithic, place all dofs in the same block', but that may not be possible if the dofs have different discretizations.
• equation sets: can we implement 'block problems'? Perhaps have a 'BlockProblem', or 'AggregateProblem', which delegates the construction of the evaluators blocks to its sub-problems...
• observer: how to save a block solution? Each block to a different discretization? Or blocks share mesh? Ideally, support both.
• evaluators: must support test BF different from trial BF. Question: should we have a specialized impl for when test=trial? It would be faster.
• Similar problem for interpolation evaluators: we could fuse all cell/side evaluators in a single evaluator. However, for HGrad, the cell interpolation does not require any metric (while the side interpolation does), so we could have a special case for HGrad (which i 100% of Albany, right now), and make it faster.
• We should probably make every problem a block problem, possibly with 1 block. However, we should support an input file that has no block information, for backward compatibility (and for easy problems). In case no block info is passed, assume 1 block, and handle it.
• In case of a single block, we should allow the user to not have teko enabled. This helps to keep dependencies of Albany under control. As a matter of fact, as of today, teko has a required dependency on all the Epetra stack... (see Trilinos#3328)

Updates (11/14/18)

@bartgol and @ikalash