(Maybe) Understanding the agglomeration procedure for eXtended Discontinuous Galerkin

Prerequisites

Loading the BoSSS-library

Note:

  1. This tutorial can be found in the source code repository as as XDGagglomeration.ipynb. One can directly load this into Jupyter to interactively work with the following code examples.
  2. In the following line, the reference to BoSSSpad.dll is required. You must either set #r "BoSSSpad.dll" to something which is appropirate for your computer (e.g. C:\Program Files (x86)\FDY\BoSSS\bin\Release\net5.0\BoSSSpad.dll if you installed the binary distribution), or, if you are working with the source code, you must compile BoSSSpad and put it side-by-side to this worksheet file (from the original location in the repository, you can use the scripts getbossspad.sh, resp. getbossspad.bat).

Implementation of a Scalar Flux

In this worksheet, Scalar Advection will be used as an example; Therefore it is required to implement the respective equation components (aka. flux) for the bulk and the interface.

Below, an upwind formulation for the equation $$ \text{div}(\vec{V} c) = 0 \text{ in } \Omega = (0,2) \times (0,1) $$ with the boundary condition $$ c = c_{\text{in}} \text{ for } \vec{V} \cdot \vec{n}_{\partial \Omega} < 0 $$ for the inflow part of the boundary is implemented:

Initialize Grid and Fields

We discretize the domain $\Omega$ with two cells $$ K_1 =(0,1) \times (0,1),~K_2 =(1,2) \times (0,1) $$ and initialize a field of degree 1 for the LevelSet and a XDGField of degree 0.

The LevelSet is set as such that the 0-isocontour of the LevelSet leads to a very small cutCell.

$$ \varphi (x,t)= -x+0.9+1.1t \quad $$

We set the XDGField to

$$c_0(c,t)=1 $$

Here we plot the field we just initalized.

Cell Agglomeration

An XDG method typically produces very small cut cells; Indeed, depending on the position of the level set, cut cells may become arbitrarily small. This typically leads to unstable methods, for both, implicit and explicit solvers.

Therefore, stabilization techniques are required; One stabilization technique is cell-agglomeration, where small cut cells are agglomerated with their largest neighbor (in the same species). Then, the XDG-problem is not solved in the original XDG-space (defined through the cut-cell mesh), but in the agglomerated XDG-space (defined through the agglomerated cut-cell-mesh).

(An alternative stabilization technique, which is not supported in BoSSS are e.g. ghost penalties.)

Here We initialize the CellAglomerator. After being initialized we can invesitgate the Volumes of the CutCells. A CutCell $K_{i,s}$ will be agglomerated (merged with another cell) if $$ \frac{\vert K_{i,s} \vert}{\vert K_{i} \vert} < agglo_{trsh}$$ In other words, if its volume divided by the Volume of the correspoding Background Cell is less than the agglomeration treshhold.

By calculating we we have:

which we can also obtain from our agglomerator.

Note that the cell volumes (as well as other cell metrics) provided by the agglomerator are valid for the agglomerated mesh. These metrics will/should be used e.g. in the penalty parameters used for Symmetric Interior Penalty (SIP).

It is, however, also possible to access the non-agglomerated metrics:

so we see that in our setting we would want to agglomerate Cut Cell $K_{2A}$ as the fraction is less then our treshhold (0.1)

Sidenote:

the standard way of handing length scales to an equation component involves two steps:

  1. Set the XEvaluatorBase.CellLengthScales to the appropriate length scale before evaluation of the operator, resp. matrix assembly
  2. Implementing the IEquationComponentCoefficient interfcae for bulk components and the ILevelSetEquationComponentCoefficient interface for level-set components. The length scales specified (see below) are passed to the CoefficientUpdate(...) method.

The Mass Matrix

For XDG, the mass matrix is quite expensive to compute. Furthermore, it also depends on the chosen quadrature order: E.g. a mass matrix for degree 2, computed with 4-th order rules is typically slightly different from a the one computed with a 6-th order rule. In order to compute it we need an Object called MassMatrixFactory.

Note, that in the simple setting presented here, the entries of the mass matrix are equal to the (non-agglomerated) cut-cell volumes.This is because:

The mass-matrix factory can also be used to obtain the inverse of the mass-matrix; these are also cached, since computing the inverse blocks is, for higher DG-orders quite expensive.

Implicit solution of an equation in the agglomerated XDG space

Definition of the operator:

Evaluation of the Residual w.r.t. the original XDG-Basis

Since the input field c is constant and complies with the boundary condition, the expected residual is zero.

Solution in the original XDG space

Here we assemble the Operator Matrix and the RHS

Then we can solve the implicit system:

$$Op(c)*c=-r(c)\approx \text{RHS}$$

Solution in the agglomerated XDG space

Next we will solve the same system but in the Aggregated space. Hence, we first need to Transform our Matrix and the RHS. This is done using the method ManipulateMatrixAndRHS(...)

As you can see BoSSS keeps the original size of the matrix. So as the cells being agglmerated vanish, we obtain zero-rows and zero-columns for these respective cells. In order to reobtain a uniquely solvable system we therefore need to artificially add a 1 on the diagonal, for the agglomerated cells.

Next, the solution in the agglomerated space can be extrapolated to the original space using the Extrapolate(...)-method .

During the extrapolation, the agglomerated coordinate will be overwritten; we illustrate this by setting the respective entry to NaN:

The mass matrix of the agglomerated space can also be obtained by passing the original mass matrix to the ManipulateMatrixAndRHS(...) method:

Notes on Extrapolation an Restriction

One can verify the extrapolation by asserting that the jump norm on agglomerated edges for e.g. a random field is zero.

Next, we obtain the agglomerated edges from the Cell Agglomerator for a specific species.

First we get the cell agglomerator for the "A" species using GetAgglomerator(LsTrk.GetSpeciesId("A")). Then we acces to AgglomerationEdges by AggInfo.AgglomerationEdges.

Then we can compute the JumpNorm on those edges by

The matrix of the restiction operation is can be accessed by:

One can also observe that the Extrapolation (aka. Prolongation, aka. Injection) is the transpose of this:

However, even when ignoring the zero row and column, which correspond to an un-used DOF in the agglomerated domain, one can observe that the Solution Extrapolation RHS-Restriction/Projection

The L2-projection to the agglomerated XDG space

In the next section we want to focus on the projection of a solution $u(x)$ from the original space to the agglomerated. While the injection is canonoical (i.e. any member of the aggregate XDG space must also be a member of the original XDG space, due to the sub-space structure), the restriction is not.

The $L^2$-Projection can be constructed in the following way; Let $Q^T = \text{RestrMtx}$ and $Q = \text{ExpolMtx}$; for the theoretical consideration here, we ommit the zero-rows and columns.

The solution is given as a coordinate vector $\tilde{u} $ such that $u(x)= \underline{\phi} \cdot \tilde{u}.$ We now want to find a coordinatevector corresponding to projection onto the agglomerated space that is $\underline{\phi}^A \cdot \tilde{u}^A$.

The agglomerate XDG basis $ \underline{\phi}^A $ is related to the original basis $\underline{\phi}$ through the relation $$ \underline{\phi}^A = \underline{\phi} Q . $$ Furthermore, let $M = \langle \underline{\phi}^T , \underline{\phi} \rangle $ be the mass matrix of the original XDG basis.

Then the mass matrix of $\underline{\phi}^A$ is given as $$ M^A := \langle {\underline{\phi}^{A}}^T , \underline{\phi}^{A} \rangle = \langle Q^T \underline{\phi}^T , \underline{\phi} Q \rangle = Q^T \langle \underline{\phi}^T , \underline{\phi} \rangle Q = Q^T M Q. $$

Now, the member $\underline{\phi}^A \cdot \tilde{u}^A$ of the agglomerated XDG space should represent the $L^2$-Projection of $\underline{\phi} \cdot \tilde{u}$ from the original XDG space. The ansatz for the projection is $$ \langle {\underline{\phi}^{A}}^T, \underline{\phi}^A \cdot \tilde{u}^A - \underline{\phi} \cdot \tilde{u} \rangle = 0, $$ i.e. we enforce orthogonality of the approximation residual $\underline{\phi}^A \cdot \tilde{u}^A - \underline{\phi} \cdot \tilde{u}$ to all elements of $ {\underline{\phi}^{A}}^T$.

We infer: $$

\langle {\underline{\phi}^{A}}^T, \underline{\phi}^{A} \rangle \tilde{u}^A

\langle {\underline{\phi}^{A}}^T, \underline{\phi} \rangle \tilde{u} $$ and further $$ M^A \tilde{u}^A = Q^T M \tilde{u}. $$ Therefore, the matrix of the $L^2$-Projection in the respective bases is given as: $$ \tilde{u}^A = ( {M^A}^{-1} Q^T M ) \tilde{u}. $$ One can easily verify that the composition of extrapolation (matrix $Q$) and $L^2$-Projection is the identity: $$ ( {M^A}^{-1} Q^T M ) Q = {M^A}^{-1} {M^A} = I. $$

Right now the MultiPhaseAgglomerator does not support this operation but the MultiGridOperator must be used instead

Solution using the multigrid operator

The multigrid operator performes a change of basis for the DG representation; Using e.g. the IdMass configuration ensures an orthonormal XDG basis.

Note: the following MultigridOperator.ChangeOfBasisConfig and will, at some point, be merged into the IDifferentialOperator interface.

Then, the mass matrix is the identity:

The operator matrix can be accesd via:

Note that the MultigridOperator also gets rid of the unused DOF's after agglomeration; Therefore, the matrices in the example above are only 3x3 after transformation.

An easy-to-use driver function UniSolver.Solve(...) allows to use the multigrid operator in order to solve the steady-state equation opsed by the IDifferentialOperator;

To solve an equation manualy, using the multigrid operator, the following recipe can be used: