Blocked bases, block matrices and block preconditioners

Many PDE systems have a natural block structure. Think of the Stokes equation with a velocity and a pressure part: it can be useful to access the block of the system matrix that implements the velocity Laplacian, or to precondition the velocity and the pressure block separately.

In dune-functions, the block structure of the linear algebra is determined by the index-merging strategy of the global basis. Composite and power bases can be combined with flat strategies (flatLexicographic(), flatInterleaved()) that produce a single integer index per basis function, or with blocked strategies (blockedLexicographic(), blockedInterleaved()) that produce multi-indices referring to an entry in a hierarchy of nested vectors and matrices. See The interface for functions in the dune-functions module for details about indexing and index-merging strategies.

This page explains what AMDiS currently supports with respect to blocking.

The AMDiS GlobalBasis is always flat

The basis type AMDiS::GlobalBasis (amdis/functions/GlobalBasis.hpp) converts every pre-basis into one with flat index-merging strategies when it is constructed. This is done by flatPreBasis() from amdis/functions/FlatPreBasis.hpp. blockedLexicographic is converted into flatLexicographic, and blockedInterleaved into flatInterleaved. The conversion is applied recursively to all nodes of the basis tree.

So you may pass a basis description with blocked strategies to a ProblemStat, DOFVector, or GlobalBasis. It is accepted, but the resulting basis uses flat indices:

using namespace Dune::Functions::BasisFactory;

// both bases are flat: blockedLexicographic is converted to flatLexicographic
GlobalBasis basis{gridView, power<2>(lagrange<1>(), blockedLexicographic())};
ProblemStat prob{"prob", grid, power<2>(lagrange<1>(), blockedLexicographic())};

The opposite transformation, flatBasis(basis) from amdis/functions/FlatBasis.hpp, is available as well and creates a flat copy of a given basis.

The block structure is still present in the basis tree. Sub-bases are addressed with tree paths, e.g., prob.addMatrixOperator(op, _v, _p) assembles into the velocity-pressure coupling block. The degrees of freedom of the different components are only stored in a common flat index space.

Mapping flat indices to blocked indices

To access the block structure of an assembled flat system, AMDiS provides utilities that relate the flat indices to the multi-indices of a corresponding blocked basis.

blockedBasis()

The function blockedBasis(basis) in amdis/functions/BlockedBasis.hpp creates a new Dune::Functions::DefaultGlobalBasis with the same tree as basis. In this new basis, only the outermost level uses a blocked index-merging strategy and all inner nodes are flat:

• a composite node is converted to BlockedLexicographic,
• a power node with FlatLexicographic or FlatInterleaved is converted to BlockedLexicographic or BlockedInterleaved, respectively,
• all child nodes get the flat analog of their strategy.

The type of the result is available as BlockedBasis_t<Basis>. The lower-level blockedPreBasis() and BlockedPreBasis_t from amdis/functions/BlockedPreBasis.hpp do the same on pre-basis level.

For a Taylor-Hood basis, this results in multi-indices of length two: the first index selects the velocity block (0) or the pressure block (1), and the second index is the flat index inside that block.

#include <amdis/functions/BlockedBasis.hpp>

ProblemStat<TaylorHoodBasis<Grid>> prob("stokes");
prob.initialize(INIT_ALL);

// a basis with indices [block, index-in-block]
auto bbasis = blockedBasis(*prob.globalBasis());

BasisMapping and BlockedBasisMapping

The class BasisMapping<SourceBasis, TargetBasis> (amdis/functions/BasisMapping.hpp) traverses the grid once and stores for each index of the source basis the corresponding multi-index of the target basis. Both bases must describe the same function space on the same grid and only differ in their indexing. If one of the bases is updated after grid adaptation, the mapping is recomputed automatically.

BlockedBasisMapping<RowBasis, ColBasis = RowBasis> (amdis/functions/BlockedBasisMapping.hpp) combines the two steps for a pair of flat row and column bases. It constructs the corresponding blocked bases with blockedBasis() and stores a flat-to-blocked index mapping for each of them:

#include <amdis/functions/BlockedBasisMapping.hpp>

using Basis = typename decltype(prob)::GlobalBasis;
BlockedBasisMapping<Basis> mapping{prob.globalBasis()};

// multi-index [block, index-in-block] of the flat row index i
auto const& r = mapping.row(i);
// multi-index [block, index-in-block] of the flat column index j
auto const& c = mapping.col(j);

// the blocked bases can be used to query the block sizes
auto const& rowSizeInfo = mapping.rowSizeInfo();

Block matrices and block preconditioners (MTL backend)

With the MTL backend (-DBACKEND=MTL), the flat-to-blocked mapping is used to split an assembled flat matrix into blocks. This is the basis for the block preconditioners described below.

BlockMatrix

BlockMatrix<Matrix> (amdis/linearalgebra/mtl/BlockMatrix.hpp) is a dense mtl::dense2D matrix of sparse sub-matrices of type Matrix. It is filled from a flat matrix and a BlockedBasisMapping:

#include <amdis/linearalgebra/mtl/BlockMatrix.hpp>

using Matrix = typename MTLSparseMatrix<double>::BaseMatrix;
Matrix const& A = prob.systemMatrix()->impl().matrix();

BlockMatrix<Matrix> B(A, mapping);       // copy all blocks
// BlockMatrix<Matrix> B(A, mapping, true); // copy only the diagonal blocks

std::size_t numRowBlocks = B.rows();     // number of row blocks
std::size_t blockSize = B.rows(0);       // number of rows in block 0
Matrix const& A00 = B[0][0];             // the sub-matrix (0,0)

The block matrix is a copy and does not change when the flat matrix changes. Call B.copyFrom(A, mapping) to update it.

Block preconditioners

The preconditioner base class BlockPreconditioner<M,V,RowBasis,ColBasis> (amdis/linearalgebra/mtl/BlockPreconditioner.hpp) stores a BlockedBasisMapping and a BlockMatrix. It provides access to the sub-matrices with subMatrix(i,j) and to sub-vector views with subVector(vec, rowIndices(i)). The following preconditioners derive from it:

BlockDiagonalPreconditioner applies one preconditioner per diagonal block, $P = \operatorname{diag}(P_0, P_1, \ldots)$. The StokesPreconditioner is a block-diagonal preconditioner for Taylor-Hood discretizations of the Stokes equation. It replaces the pressure block by the pressure mass matrix, which approximates the Schur complement. The SIMPLE and SIMPLER preconditioners are classic block preconditioners for the (Navier-)Stokes equations.

These preconditioners need the global basis to compute the block structure. So they must be registered with the basis before the linear solver is created. Split the initialization of the problem into two steps:

#include <amdis/linearalgebra/mtl/precons/StokesPreconditioner.hpp>

ProblemStat<TaylorHoodBasis<Grid>> prob("stokes");
prob.initialize(INIT_ALL - INIT_SOLVER);

// register the preconditioner under the name "stokes"
registerStokesPreconditioner(prob.globalBasis());

// now create the solver, which reads the preconditioner from the initfile
prob.initialize(INIT_SOLVER);

The preconditioner and the preconditioners of the diagonal blocks are selected in the initfile. Each block preconditioner is configured with the prefix [solver-prefix]->precon->sub precon[i]. The block preconditioner can be an inner iterative or direct solver (key solver) or any other registered preconditioner:

stokes->solver: gcr
stokes->solver->relative tolerance: 1e-8
stokes->solver->precon: stokes

% velocity block: inner direct solver
stokes->solver->precon->sub precon[0]: solver
stokes->solver->precon->sub precon[0]->solver: direct

% pressure (mass matrix) block
stokes->solver->precon->sub precon[1]: solver
stokes->solver->precon->sub precon[1]->solver: direct

The complete examples are in examples/stokes_mtl.cc and examples/navier_stokes_mtl.cc, with the parameter files examples/init/stokes_mtl.dat.2d and examples/init/navier_stokes_mtl.dat.2d. They are built only with the MTL backend.

Discrete functions on blocked dune-functions bases

A DiscreteFunction (amdis/gridfunctions/DiscreteFunction.hpp) can also wrap a coefficient container that you manage yourself, together with a plain dune-functions basis. In this case the basis is used as given. A blocked basis is combined with a nested container:

using namespace Dune::Functions::BasisFactory;
auto basis = makeBasis(gridView, power<2>(lagrange<2>(), blockedInterleaved()));

Dune::BlockVector<Dune::FieldVector<double,2>> coeffs(basis.size());
DiscreteFunction u{coeffs, basis};

This is useful to evaluate or post-process data that is stored in blocked containers, e.g., data from other Dune modules. Interpolation into a DiscreteFunction (operator<<, interpolate) requires an AMDiS vector type (a VectorFacade, e.g., the coefficients of a DOFVector), so it is not available for such plain containers. These functions are not connected to the linear system of a ProblemStat.