|
Teko Version of the Day
|
This page collects the mechanisms you reach for once the built-in preconditioners and their options are not quite enough: supplying operators through the RequestHandler, reordering blocks, reusing expensive inverses across nonlinear iterations, and registering your own preconditioner factory so it can be named from XML.
Some preconditioners need operators that cannot be derived from the assembled system — for example a velocity mass matrix (SIMPLE/LSC), a pressure Laplacian, or a probing graph. Rather than complicate every factory's constructor, Teko uses a request/callback pattern:
Teko::RequestMesg (a named request, e.g. "Velocity Mass Matrix") when it needs the operator.Teko::RequestCallback<DataT> with a RequestHandler that answers matching messages.The callback base (src/Teko_RequestCallback.hpp) has three methods:
A minimal callback that supplies a velocity mass matrix:
Register it and attach the handler to whatever builds the preconditioner:
| Message | Issued by | Enabled when |
|---|---|---|
"Velocity Mass Matrix" | LSC Basic Inverse, SIMPLE | "Use Mass Scaling" = true |
"W-Scaling Vector" | LSC Basic Inverse | "Use W-Scaling" = true |
"Pressure Laplace Operator" | LSC Pressure Laplace, PCD | always (for those strategies) |
"Velocity Mass Operator" | LSC Pressure Laplace | always |
"PCD Operator" | PCD strategy | always |
"Probing Graph" | Probing Preconditioner | "User Will Set Probing Graph" = true |
The "Required Parameters" sublist of a Stratimikos-preconditioner library entry is the parameter-driven counterpart: its contents are removed and re-supplied through the RequestHandler at build time.
Beyond the top-level "Reorder Type" string (Configuration Model), you can reorder a block operator programmatically. Teko::blockedReorderFromString parses the nested-bracket grammar into a BlockReorderManager, and Teko::buildReorderedLinearOp applies it (src/Teko_BlockedReordering.hpp):
The bracket string mirrors the block tree: bare integers are leaf block indices, and [ ... ] groups its contents into a nested block. "[2 [0 1]]" builds a 2×2 outer operator whose (0,0) entry is block 2 and whose (1,1) entry is the nested 2×2 of blocks 0 and 1.
In a Newton / nonlinear loop the operator changes every iteration but its sparsity and the preconditioner structure often do not. Rebuild the inverse in place instead of reconstructing it, which lets factories reuse symbolic factorizations, prolongation operators, and cached explicit matrices when the chosen backend supports reuse:
A BlockPreconditionerFactory may build preconditioners for many different blocked operators. Do not store operator-specific objects, such as cached inverses, directly in the factory unless the factory is intentionally single-use. Use the BlockPreconditionerState& state argument to buildPreconditionerOperator instead; Teko pairs that state object with the specific operator instance being built or rebuilt.
For cached inverse operators, a common pattern is:
state.getModifiableOp("invP") returns a reference to a stored RCP. Assigning through that reference updates the object kept by the state. The next rebuild for the same preconditioner instance retrieves the same handle and can refresh it in place.
If the default state is not enough, derive from Teko::BlockPreconditionerState and override your factory's buildPreconditionerState() method to return the richer state object. Inside buildPreconditionerOperator, dynamically cast the incoming state to that derived type before using the extra fields.
The explicit utilities have overloads that take a destination ModifiableLinearOp, for example explicitAdd(A, B, dest) and explicitMultiply(A, B, dest). Use these when a preconditioner forms the same auxiliary operator every nonlinear iteration and only the values change:
If P is null, Teko allocates it. If it is compatible with the new inputs, Teko can refill and reuse it; otherwise it falls back to constructing a new operator. This is especially useful for Schur-complement approximations and mass-scaled operators that would otherwise be reallocated every rebuild.
Once you have written a BlockPreconditionerFactory (see the step1 example), register it so it can be selected by a "Type" string from the Inverse Library:
After this, "Type" = "My Block Prec" works anywhere a Teko block type is expected, and your factory's initializeFromParameterList receives the entry's remaining keys. Override initializeFromParameterList to read your own parameters, and (optionally) provide getRequestHandler()-based lookups for any operators you need from the application.
After implementing a new preconditioner factory, first test its action on known vectors before using it inside a Krylov solve. A good unit test builds a small block operator, applies the preconditioner to a predetermined vector, and compares the result with an independently computed reference result (for example from a small script or a direct hand calculation). Then add an integration test that uses the preconditioner inside Belos on a representative matrix.
Helpful diagnostics include:
"Diagnostic Inverse"** (reference) wraps any inner inverse to time it and optionally print its residual — drop it in temporarily to find the expensive sub-solve."Test Block Operator"** and **"Write Block Operator"** (top-level parameters) verify and dump the segregated blocks when you suspect the strided blocking or reordering is wrong.StridedTpetraOperator::testAgainstFullOperator and BlockedTpetraOperator::testAgainstFullOperator compare the wrapped block operator with the original monolithic Tpetra::Operator on random vectors.Teuchos::describe(*op, Teuchos::VERB_EXTREME) prints the structure and, for some concrete operators, values of a Teko::LinearOp.Teko's core operator and vector types are aliases for Teuchos::RCP smart-pointer handles around Thyra objects:
A BlockedLinearOp handle can be upcast to a LinearOp handle, and a BlockedMultiVector handle can be upcast to a MultiVector handle. The helper casts in Teko_Utilities.hpp make those relationships explicit:
Because these types are Teuchos::RCP handles, assignment is a shallow copy: it copies the handle, not the matrix or vector data. Use the explicit deepcopy helpers when you need an independent vector object:
A few small utilities are often useful when constructing custom operators:
Wrap native matrices exactly once (Thyra::tpetraLinearOp), then stay in Teko/Thyra types. Block assembly (block2x2, zeroBlockedOp/setBlock/endBlockFill) produces Teko::BlockedLinearOp, which the factories consume directly. The TpetraBlockPreconditioner wrapper exists so a finished Teko preconditioner can be handed back to a solver that only speaks the Tpetra Operator interface.
1.9.8