What's new?

• implementation of a saddle point problem into BoSSS (e.g. the Poisson equation as a system)
• central-difference-form and strong form
• comparison of both formulations
• advanced: algebraic reduction of the poisson problem

Prerequisites

• implementation of numerical fluxes -> chapter NumFlux
• spatial operator -> chapter DifferentialOperator

1 Problem statement

Within this exercise, we are going to investigate the discretization of a Poisson equation as a system. Obviously, it is possible to discretize the Poisson equation as a system of first-order-PDE's, introducing a vector field $\vec{\sigma}$: $$ \begin{align} \vec{\sigma} + \nabla u & = 0, & & \text{ in } \Omega \ \operatorname{div}(\vec{\sigma}) & = g_{\Omega}, & & \text{ in } \Omega \ u & = g_D, & & \text{ on } \Gamma_D \

• \vec{\sigma} \cdot \vec{n}_{\partial \Omega} & = g_N, & & \text{ on } \GammaN \end{align} $$ resp. in matrix-notation: $$ \begin{align*} \begin{bmatrix} 1 & \nabla \\ \operatorname{div} & 0 \\ \end{bmatrix}\cdot \begin{bmatrix} \vec{\sigma}\\ u \end{bmatrix}= \begin{bmatrix} 0 \ g{\Omega} \end{bmatrix} \end{align*} $$ This exercise, together with the previous one, will form the foundation for an incompressible Stokes- resp. Navier-Stokes solver.

2 Solution within the BoSSS framework

#r "BoSSSpad.dll"
using System;
using System.Collections.Generic;
using System.Linq;
using ilPSP;
using ilPSP.Utils;
using BoSSS.Platform;
using BoSSS.Foundation;
using BoSSS.Foundation.Grid;
using BoSSS.Foundation.Grid.Classic;
using BoSSS.Foundation.IO;
using BoSSS.Solution;
using BoSSS.Solution.Control;
using BoSSS.Solution.GridImport;
using BoSSS.Solution.Statistic;
using BoSSS.Solution.Utils;
using BoSSS.Solution.Gnuplot;
using BoSSS.Application.BoSSSpad;
using BoSSS.Application.XNSE_Solver;
using static BoSSS.Application.BoSSSpad.BoSSSshell;
Init();

The below script needs to be able to find the current output cell; this is an easy method to get it.

using ilPSP.LinSolvers; 
using BoSSS.Solution.Tecplot; 
using ilPSP.Connectors.Matlab;

2.1 Tests on the divergence

Common base-class for $\text{div}$-implementations

We are going to implement two different formulations of the divergence-operator for which going to show equivalence. We implement a common base-class for both formulations:

abstract public class BaseDivergence :  
        BoSSS.Foundation.IEdgeForm, // edge integrals 
        BoSSS.Foundation.IVolumeForm     // volume integrals 
{ 
    /// We don't use parameters (e.g. variable viscosity, ...)
    /// at this point: so the parameter list can be null, resp. empty:
    public IList<string> ParameterOrdering {  
        get { return null; }  
    } 
 
    /// But we have a vector argument variable, 
    /// $ [ \sigma_1, \sigma_2 ] = \vec{\sigma} $
    /// (our trial function):
    public IList<String> ArgumentOrdering {  
        get { return new string[] { "sigma1", "sigma2" }; }  
    } 
 
    public TermActivationFlags VolTerms { 
        get { 
            return TermActivationFlags.AllOn; 
        } 
    } 
 
    public TermActivationFlags InnerEdgeTerms { 
        get { 
            return (TermActivationFlags.AllOn);  
        } 
    } 
 
    public TermActivationFlags BoundaryEdgeTerms { 
       get { 
           return TermActivationFlags.AllOn; 
        } 
    } 
 
    /// The following functions cover the actual math.
    /// For any discretization of the divergence-operator, we have to specify:
    /// \begin{itemize}
    ///    \item a volume integrand,
    ///    \item an edge integrand for inner edges, i.e. on $ \Gamma_i$,
    ///    \item an edge integrand for boundary edges, 
    ///          i.e. on $\partial \Omega$.
    /// \end{itemize}
    /// These functions are declared as \code{abstract}, meaning that one has 
    /// to specify them in classes derived from \code{BaseLaplace}.
 
    abstract public double VolumeForm(ref CommonParamsVol cpv,  
           double[] U, double[,] GradU,  
           double V, double[] GradV);         
 
    abstract public double InnerEdgeForm(ref CommonParams inp,  
        double[] U_IN, double[] U_OT, double[,] GradU_IN, double[,] GradU_OT,  
        double V_IN, double V_OT, double[] GradV_IN, double[] GradV_OT); 
 
    abstract public double BoundaryEdgeForm(ref CommonParamsBnd inp,  
        double[] U_IN, double[,] GradU_IN, double V_IN, double[] GradV_OT); 
}

/// We are going to use both, Dirichlet- and Neumann-boundary conditions
/// in this exercise; the function \code{IsDirichletBndy} is used to
/// specify the type of boundary condition at point \code{X}:
static class BndyMap {
static public Func<double[],bool> IsDirichletBndy = delegate(double[] X) { 
    double x = X[0]; 
    double y = X[1]; 
    if(Math.Abs(x - (-1.0)) < 1.0e-8) 
        return true;     
    if(Math.Abs(y - (-1.0)) < 1.0e-8) 
        return true;     
    return false; 
};
}

Formulation (i): Central-difference-form of $\text{div}$

The implementation of the central-difference form is as follows:

class Divergence_cendiff : BaseDivergence { 
 
 
    /// The volume form is equal to 
    /// $ -\vec{\sigma} \cdot \nabla v$:
    override public double VolumeForm(ref CommonParamsVol cpv,  
        double[] Sigma, double[,] GradSigma,  
        double V, double[] GradV) { 
        double Acc = 0; 
        for(int d = 0; d < cpv.D; d++) { 
            Acc -= Sigma[d]*GradV[d]; 
        } 
        return Acc; 
    } 
 
    /// At the cell boundaries, we use a central-difference-flux,
    /// i.e. $\mean{\vec{\sigma}} \cdot \vec{n}_{\Gamma} \jump{v}$:
    override public double InnerEdgeForm(ref CommonParams inp,  
        double[] Sigma_IN, double[] Sigma_OT, double[,] GradSigma_IN, double[,] GradSigma_OT,  
        double V_IN, double V_OT, double[] GradV_IN, double[] GradV_OT) { 
 
        double Acc = 0; 
        for(int d = 0; d < inp.D; d++) { 
            Acc += 0.5*(Sigma_IN[d] + Sigma_OT[d])*inp.Normal[d]*(V_IN - V_OT); 
        } 
        return Acc; 
    } 
 
    override public double BoundaryEdgeForm(ref CommonParamsBnd inp,  
        double[] Sigma_IN, double[,] GradSigma_IN, double V_IN, double[] GradV_OT) { 
 
        double Acc = 0; 
 
        if(BndyMap.IsDirichletBndy(inp.X)) { 
            /// Dirichlet-boundary: by taking the inner value of $\vec{\sigma}$, 
            /// this is a free boundary with respect to $\vec{\sigma}$.
            for(int d = 0; d < inp.D; d++) { 
                Acc += Sigma_IN[d]*inp.Normal[d]*V_IN; 
            } 
        } else { 
            /// Neumann-boundary
            double gNeu = 0.0; 
            Acc += gNeu*V_IN; 
        } 
        return Acc; 
    } 
}

Formulation (ii): 'Strong' form of $\text{div}$:

Here, we use the form $$ b(\vec{\sigma},v) = \oint_{\Gamma \backslash \GammaD} M(v) J(\vec{\sigma}) \cdot \vec{n}\Gamma

dA

\int_{\Omega} \text{div}(\vec{\sigma}) \cdot v dV $$ where M,J denote the mean and jump operator, respectively. This is actually the negative divergence, which will be more useful later on.

class Divergence_strong : BaseDivergence { 
 
    /// We have to implement \code{VolumeForm},
    /// \emph{InnerEdgeForm} and \code{BoundaryEdgeForm}:
    override public double VolumeForm(ref CommonParamsVol cpv,  
        double[] Sigma, double[,] GradSigma,  
        double V, double[] GradV) { 
        double Acc = 0; 
        for(int d = 0; d < cpv.D; d++) { 
            Acc -= GradSigma[d,d]*V; 
        } 
        return Acc; 
    } 
 
    override public double InnerEdgeForm(ref CommonParams inp,  
        double[] Sigma_IN, double[] Sigma_OT, double[,] GradSigma_IN, double[,] GradSigma_OT,  
        double V_IN, double V_OT, double[] GradV_IN, double[] GradV_OT) { 
 
        double Acc = 0; 
        for(int d = 0; d < inp.D; d++) { 
            Acc += 0.5*(V_IN + V_OT)*(Sigma_IN[d] - Sigma_OT[d])*inp.Normal[d]; 
        } 
        return Acc; 
    } 
 
    override public double BoundaryEdgeForm(ref CommonParamsBnd inp,  
        double[] Sigma_IN, double[,] GradSigma_IN, double V_IN, double[] GradV_OT) { 
 
        double Acc = 0; 
 
        if(BndyMap.IsDirichletBndy(inp.X)) { 
            Acc = 0;
        } else { 
            double gNeu = 0.0; 
            for(int d = 0; d < inp.D; d++) { 
                Acc += Sigma_IN[d]*inp.Normal[d]*V_IN; 
            } 
            Acc -= gNeu*V_IN; 
        } 
        return Acc; 
    } 
}

3 Equality test

We are going to test the equivalence of both formulationt on a 2D grid, using a DG basis of degree 1:

var grd2D               = Grid2D.Cartesian2DGrid(GenericBlas.Linspace(-1,1,6),                                                    GenericBlas.Linspace(-1,1,7)); 
 var b                   = new Basis(grd2D, 1); 
 SinglePhaseField sigma1 = new SinglePhaseField(b,"sigma1"); 
 SinglePhaseField sigma2 = new SinglePhaseField(b,"sigma2"); 
 CoordinateVector sigma  = new CoordinateVector(sigma1,sigma2); 
 var TrialMapping        = sigma.Mapping; 
 var TestMapping         = new UnsetteledCoordinateMapping(b);

/// We create the matrix of the central-difference formulation:
 var OpDiv_cendiff = (new Divergence_cendiff()).Operator(); 
 var MtxDiv_cendiff = OpDiv_cendiff.ComputeMatrix(TrialMapping,  
                                                  null,  
                                                  TestMapping);

We create the matrix of the strong formulation and show that the matrices of both formulations are equal.

We use the \code{InfNorm(...)}-method to identify whether a matrix is (approximately) zero or not.

var OpDiv_strong  = (new Divergence_strong()).Operator(); 
 var MtxDiv_strong = OpDiv_strong.ComputeMatrix(TrialMapping, null, TestMapping); 
 var TestP         = MtxDiv_cendiff + MtxDiv_strong; 
 TestP.InfNorm();

4 The gradient-operator

For the variational formulation of the gradient operator, a vector-valued test-function is required. Unfourtunately, this is not supported by BoSSS. Therefore we have to discretize the gradent component-wise, i.e. as $\partial_{x}$ and $\partial_y$.

A single derivative can obviously be expressed as a divergence by the identity $ \partial_{x_d} = \text{div}( \vec{e}_d u ) $.

class Gradient_d : 
        BoSSS.Foundation.IEdgeForm, // edge integrals 
        BoSSS.Foundation.IVolumeForm     // volume integrals 
{ 
    public Gradient_d(int _d) { 
        this.d = _d; 
    } 
 
    /// The component index of the gradient:
    int d; 
 
    /// As ususal, we do not use parameters:
    public IList<string> ParameterOrdering {  
        get { return null; }  
    } 
 
    /// We have one argument $u$:
    public IList<String> ArgumentOrdering {  
        get { return new string[] { "u" }; }  
    } 
 
    public TermActivationFlags VolTerms { 
        get { return TermActivationFlags.AllOn; } 
    } 
 
    public TermActivationFlags InnerEdgeTerms { 
        get { return (TermActivationFlags.AllOn); } 
    } 
 
    public TermActivationFlags BoundaryEdgeTerms { 
       get { return TermActivationFlags.AllOn; } 
    } 
 
    /// Now, we implement 
    /// \begin{itemize}
    ///    \item the volume form $u \vec{e}_d \cdot \nabla v$
    ///    \item the boundary form 
    ///       $\mean{u \ \vec{e}_d} \cdot \vec{n}_\Gamma \jump{v}$
    /// \end{itemize}
    public double VolumeForm(ref CommonParamsVol cpv,  
           double[] U, double[,] GradU,  
           double V, double[] GradV) { 
 
        double Acc = 0; 
        Acc -= U[0]*GradV[this.d]; 
        return Acc; 
    }         
 
    public double InnerEdgeForm(ref CommonParams inp,  
        double[] U_IN, double[] U_OT, double[,] GradU_IN, double[,] GradU_OT,  
        double V_IN, double V_OT, double[] GradV_IN, double[] GradV_OT) { 
 
        double Acc = 0; 
        Acc += 0.5*(U_IN[0] + U_OT[0])*inp.Normal[this.d]*(V_IN - V_OT); 
        return Acc;   
     } 
 
    public double BoundaryEdgeForm(ref CommonParamsBnd inp,  
        double[] U_IN, double[,] GradU_IN, double V_IN, double[] GradV_OT) { 
 
        double Acc = 0; 
        if(BndyMap.IsDirichletBndy(inp.X)) { 
            double u_Diri = 0.0; 
            Acc += u_Diri*inp.Normal[this.d]*V_IN; 
        } else { 
            Acc += U_IN[0]*inp.Normal[this.d]*V_IN; 
        } 
        return Acc;               
    } 
}

Now, we are ready to assemble the full $\nabla$ operator as $\left[ \begin{array}{c} \partial_x \\ \partial_y \end{array} \right]$.

var OpGrad = new DifferentialOperator(1,2,QuadOrderFunc.Linear(),"u","c1","c2"); 
OpGrad.EquationComponents["c1"].Add(new Gradient_d(0)); 
OpGrad.EquationComponents["c2"].Add(new Gradient_d(1)); 
OpGrad.Commit();

As an additional test, we create the gradient-matrix and verify that its transpose is equal to the negative MtxDiv-matrix:

var MtxGrad = OpGrad.ComputeMatrix(TestMapping, null, TrialMapping); 
var Test2   = MtxGrad.Transpose() - MtxDiv_strong; 
Test2.InfNorm();

5 The complete Poisson-system

Assembly of the system

We also need the identity-matrix in the top-left corner of the Poisson-system:

public class Identity :  
        BoSSS.Foundation.IVolumeForm  
{ 
    public IList<string> ParameterOrdering {  
        get { return new string[0]; }  
    } 
 
    public string component;  
 
    public IList<String> ArgumentOrdering {  
        get { return new string[] { component }; }  
    } 
 
    public TermActivationFlags VolTerms { 
        get { 
            return TermActivationFlags.AllOn; 
        } 
    } 
 
    public double VolumeForm(ref CommonParamsVol cpv,  
           double[] U, double[,] GradU,  
           double V, double[] GradV) { 
        return U[0]*V;            
    } 
}

(9,19): warning CS0649: Field 'Identity.component' is never assigned to, and will always have its default value null

We are going to implement the linear Poisson-operator $$ \left[ \begin{array}{ccc} 1 & 0 & \partial_x \\ 0 & 1 & \partial_y \\ -\partial_x & -\partial_y & 0 \end{array} \right] \cdot

\left[ \begin{array}{c} \sigma_0 \\ \sigma_1 \\ u \end{array} \right]

\left[ \begin{array}{c} c_0 \\ c_1 \\ c_2 \end{array} \right] $$ The variables $c_0$, $c_1$ and $c_2$, which correspond to the test functions are also called co-domain variables of the operator. We are using the negative divergence, since this will lead to a symmetric matrix, instead of a anti-symmetric one. By doing so, we can e.g. use a Cholesky-factorization to determine whether the system is definite or not.

var OpPoisson = new DifferentialOperator(3, 3,  
                      QuadOrderFunc.Linear(), 
                      "sigma1", "sigma2", "u", // the domain-variables 
                      "c1", "c2", "c3"); //       the co-domain variables  
/// Now we add all required components to \code{OpPoisson}:
OpPoisson.EquationComponents["c1"].Add(new Gradient_d(0)); 
OpPoisson.EquationComponents["c1"].Add(new Identity() { component = "sigma1" }); 
OpPoisson.EquationComponents["c2"].Add(new Gradient_d(1)); 
OpPoisson.EquationComponents["c2"].Add(new Identity() { component = "sigma2" }); 
OpPoisson.EquationComponents["c3"].Add(new Divergence_strong()); 
OpPoisson.Commit();

We create mappings $[\sigma_1, \sigma_2, u ]$: three different combinations of DG orders will be investigated:

• equal order: the same polynomial degree for $u$ and $\vec{\sigma}$
• mixed order: the degree of $u$ is lower than the degree

  of $\vec{\sigma}$.

• `strange' order: the degree of $u$ is higher than the degree of

  $\vec{\sigma}$.

var b3         = new Basis(grd2D, 3); 
var b2         = new Basis(grd2D, 2); 
var b4         = new Basis(grd2D, 4); 
var EqualOrder = new UnsetteledCoordinateMapping(b3,b3,b3); 
var MixedOrder = new UnsetteledCoordinateMapping(b4,b4,b3); 
var StrngOrder = new UnsetteledCoordinateMapping(b2,b2,b3);

var MtxPoisson_Equal = OpPoisson.ComputeMatrix(EqualOrder, null, EqualOrder); 
var MtxPoisson_Mixed = OpPoisson.ComputeMatrix(MixedOrder, null, MixedOrder); 
var MtxPoisson_Strng = OpPoisson.ComputeMatrix(StrngOrder, null, StrngOrder);

We show that the matrices are symmetric (use e.g. SymmetryDeviation(...)), but indefinite (use e.g. IsDefinite(...)).

double symDev_Equal = MtxPoisson_Equal.SymmetryDeviation();
symDev_Equal

2.0920765120280294E-14

double symDev_Mixed = MtxPoisson_Mixed.SymmetryDeviation();
symDev_Mixed

3.755329380794592E-14

double symDev_Strng = MtxPoisson_Strng.SymmetryDeviation();
symDev_Strng

2.090688733247248E-14

MtxPoisson_Equal.IsDefinite()

False

MtxPoisson_Mixed.IsDefinite()

False

MtxPoisson_Strng.IsDefinite()

False

/// BoSSScmdSilent BoSSSexeSilent
NUnit.Framework.Assert.LessOrEqual(symDev_Equal, 1.0e-8);
NUnit.Framework.Assert.LessOrEqual(symDev_Mixed, 1.0e-8);
NUnit.Framework.Assert.LessOrEqual(symDev_Strng, 1.0e-8);

6 Advanced topics

Algebraic reduction

Since the top-left corner of our matrix $$ \left[ \begin{array}{cc} 1 & B \\ B^T & 0 \end{array} \right] $$ is actually very easy to eliminate the variable $\vec{\sigma}$ from our system algebraically. The matrix of the reduces system is obviously $B^T \cdot B$.

Extraction of sub-matrices and elimination

From the mapping, we can actually obtain index-lists for each variable, which can then be used to extract sub-matrices from MtxPoisson_Equal, MtxPoisson_Mixed, resp. MtxPoisson_Strng.

long[] SigmaIdx_Equal = EqualOrder.GetSubvectorIndices(true, 0,1); 
long[] uIdx_Equal     = EqualOrder.GetSubvectorIndices(true, 2); 
long[] SigmaIdx_Mixed = MixedOrder.GetSubvectorIndices(true, 0,1); 
long[] uIdx_Mixed     = MixedOrder.GetSubvectorIndices(true, 2); 
long[] SigmaIdx_Strng = StrngOrder.GetSubvectorIndices(true, 0,1); 
long[] uIdx_Strng     = StrngOrder.GetSubvectorIndices(true, 2);

The extraction of the sub-matrix and the elimination, for the equal order:

var MtxPoissonRed_Equal =  
      MtxPoisson_Equal.GetSubMatrix(uIdx_Equal, SigmaIdx_Equal)  // -Divergence 
    * MtxPoisson_Equal.GetSubMatrix(SigmaIdx_Equal, uIdx_Equal); // Gradient

Finally, we also create the reduced system for the mixed and the strange order, test for the definiteness of the reduced system.

Equal and mixed order are positive definite, while the strange order is indefinite - a clear indication that something ist wrong:

var MtxPoissonRed_Mixed =  
      MtxPoisson_Mixed.GetSubMatrix(uIdx_Mixed, SigmaIdx_Mixed)  // -Divergence 
    * MtxPoisson_Mixed.GetSubMatrix(SigmaIdx_Mixed, uIdx_Mixed); // Gradient   
var MtxPoissonRed_Strng =  
      MtxPoisson_Strng.GetSubMatrix(uIdx_Strng, SigmaIdx_Strng)  // -Divergence 
    * MtxPoisson_Strng.GetSubMatrix(SigmaIdx_Strng, uIdx_Strng); // Gradient

bool isdef_red_Equal = MtxPoissonRed_Equal.IsDefinite();
isdef_red_Equal

True

bool isdef_red_Mixed = MtxPoissonRed_Mixed.IsDefinite();
isdef_red_Mixed

True

bool isdef_red_Strng = MtxPoissonRed_Strng.IsDefinite();
isdef_red_Strng

False

/// BoSSScmdSilent BoSSSexeSilent
NUnit.Framework.Assert.IsTrue(isdef_red_Equal);
NUnit.Framework.Assert.IsTrue(isdef_red_Mixed);
NUnit.Framework.Assert.IsFalse(isdef_red_Strng);

We compute the condition number of all three matrices; we observe that the mixed as well as the equal-order discretization result give rather moderate condition numbers.

For the strange orders, the condition number of the system is far to high:

double condest_Mixed = MtxPoissonRed_Mixed.condest();
condest_Mixed

7016.828053419772

double condest_Equal = MtxPoissonRed_Equal.condest();
condest_Equal

3577.101349907232

double condest_Strng = MtxPoissonRed_Strng.condest();
condest_Strng

1.224410614249654E+19

/// BoSSScmdSilent BoSSSexeSilent
NUnit.Framework.Assert.LessOrEqual(condest_Mixed, 1e5);
NUnit.Framework.Assert.LessOrEqual(condest_Equal, 1e5);
NUnit.Framework.Assert.Greater(condest_Strng, 1e10);