We consider the 2D Poisson problem: $$ \Delta u = f(x,y) $$
where $f(x,y) \neq 0$ is an arbitrary function of $x$ and/or $y$. Within this exercise, we are going to investigate the Symmetric Interior Penalty discretization method (SIP) for the Laplace operator
$$a_{\text{sip}}(u,v) = \int_{\Omega} \underbrace{\nabla u \cdot \nabla v}_{\text{Volume\ term}}dV - \oint_{\Gamma \setminus \Gamma_{N }} \underbrace{ M(\nabla u) \cdot n_{\Gamma}J(v) }_{\text{consistency term}} + \underbrace{ M(\nabla v) \cdot \vec{n}_{\Gamma} J(u) }_{\text{symmetry term}} dA + \oint_{\Gamma \setminus \Gamma_{N}} \underbrace{ \eta J(u)J(v) }_{\text{penalty term}} dA$$Where $M$ shall denote the Mean and $J$ the Jump operator. The use of these fluxes including a penalty term stabilizes the DG-discretization for the Laplace operator.
First, we initialize the new worksheet; Note:
sip.ipynb.
One can directly load this into Jupyter to interactively work with the following code examples.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).#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();
We need the following packages:
using ilPSP.LinSolvers;
using ilPSP.Connectors.Matlab;
BoSSScmdSilent BoSSSexeSilent
using NUnit.Framework;
We are going to implement the SIP-form $$ a{\text{sip}}(u,v) = \int{\Omega} \underbrace{\nabla u \cdot \nabla v}_{\text{volume\ term}}dV
M {\nabla u} \cdot n_{\Gamma}J(v)
M{\nabla v} \cdot \vec{n}_{\Gamma} J(u)
public class SipLaplace :
BoSSS.Foundation.IEdgeForm, // edge integrals
BoSSS.Foundation.IVolumeForm, // volume integrals
IEquationComponentCoefficient // update of coefficients (e.g. length scales) required for penalty parameters
{
/// We do not use parameters (e.g. variable viscosity, ...)
/// at this point: so this can be null
public IList<string> ParameterOrdering {
get { return new string[0]; }
}
/// but we have one argument variable, $u$ (our trial function)
public IList<String> ArgumentOrdering {
get { return new string[] { "u" }; }
}
/// The \code{TermActivationFlags} tell \BoSSS which kind of terms,
/// i.e. products of u, v, \nabla u, and \nabla v
/// the VolumeForm(...) actually contains.
/// This additional information helps to improve the performance.
public TermActivationFlags VolTerms {
get {
return TermActivationFlags.GradUxGradV;
}
}
/// activation flags for the 'InnerEdgeForm(...)'
public TermActivationFlags InnerEdgeTerms {
get {
return (TermActivationFlags.AllOn);
// if we do not care about performance, we can activate all terms.
}
}
public TermActivationFlags BoundaryEdgeTerms {
get {
return TermActivationFlags.AllOn;
}
}
/// For the computation of the penalty factor $\eta$,
/// we require
/// some length scale for each cell and
/// the polynomial degree of the DG approximation.
MultidimensionalArray cj;
double penalty_base; // base factor must be scaled by polynomial degree
/// The additional scaling of the penalty by polynomial degree
/// and in depencence of geometry can be obtained through
/// impmenting the \code{IEquationComponentCoefficient} interface:
public void CoefficientUpdate(CoefficientSet cs, int[] DomainDGdeg, int TestDGdeg) {
int D = cs.GrdDat.SpatialDimension;
double _D = D;
double _p = DomainDGdeg.Max();
double penalty_deg_tri = (_p + 1) * (_p + _D) / _D; // formula for triangles/tetras
double penalty_deg_sqr = (_p + 1.0) * (_p + 1.0); // formula for squares/cubes
penalty_base = Math.Max(penalty_deg_tri, penalty_deg_sqr); // the conservative choice
//Console.WriteLine("Setting penalty base factor for deg " + _p + " to " + penalty_base);
cj = ((GridData)(cs.GrdDat)).Cells.cj;
}
/// The safety factor for the penalty factor should be in the order of 1.
/// A very large penalty factor increases the condition number of the
/// system, but without affecting stability.
/// A very small penalty factor yields to an unstable discretization.
public double PenaltySafety = 2.2;
/// The actual computation of the penalty factor, which should be
/// used in the \code{InnerEdgeForm} and \code{BoundaryEdgeForm} functions.
/// Hint: for the parameters \code{jCellIn}, \code{jCellOut} and \code{g},
/// take a look at
/// \code{CommonParams} and \code{CommonParamsBnd}.
double PenaltyFactor(int jCellIn, int jCellOut) {
double cj_in = cj[jCellIn];
double eta = penalty_base * cj_in * PenaltySafety;
if(jCellOut >= 0) {
double cj_out = cj[jCellOut];
eta = Math.Max(eta, penalty_base * cj_out * PenaltySafety);
}
return eta;
}
/// The following functions cover the actual math.
/// For any discretization of the Laplace operator, we have to specify:
///
/// - a volume integrand,
/// - an edge integrand for inner edges, i.e. on $\Gamma_i$,
/// - an edge integrand for boundary edges, i.e. on $\partial \Omega$.
///
/// The integrand for the volume integral:
public double VolumeForm(ref CommonParamsVol cpv,
double[] U, double[,] GradU, // the trial-function u
// (i.e. the function we search for) and its gradient
double V, double[] GradV // the test function;
) {
double acc = 0;
for(int d = 0; d < cpv.D; d++)
acc += GradU[0, d] * GradV[d];
return acc;
}
/// The integrand for the integral on the inner edges,
///
/// -( M{\nabla u} J{v}) \cdot \vec{n}_{\Gamma}
/// -( M{\nabla v} J{u}) \cdot \vec{n}_{\Gamma}
/// + \eta J{u} J{v} :
///
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 eta = PenaltyFactor(inp.jCellIn, inp.jCellOut);
double Acc = 0.0;
for(int d = 0; d < inp.D; d++) { // loop over vector components
// consistency term: -({{ \/u }} [[ v ]])*Normal
// index d: spatial direction
Acc -= 0.5 * (GradU_IN[0, d] + GradU_OT[0, d])*(V_IN - V_OT)
* inp.Normal[d];
// symmetry term: -({{ \/v }} [[ u ]])*Normal
Acc -= 0.5 * (GradV_IN[d] + GradV_OT[d])*(U_IN[0] - U_OT[0])
* inp.Normal[d];
}
// penalty term: eta*[[u]]*[[v]]
Acc += eta*(U_IN[0] - U_OT[0])*(V_IN - V_OT);
return Acc;
}
/// The integrand on boundary edges, i.e. on $\partial \Omega$, is
///
/// -( M{\nabla u} J{v}) \cdot \vec{n}_{\Gamma}
/// -( M{\nabla v} J{u}) \cdot \vec{n}_{\Gamma}
/// + \eta J{u} J{v} .
///
/// For the boundary we have to consider the special definition for
/// the mean-value operator $M{-}$ and the jump operator
/// $J{-}$ on the boundary.
public double BoundaryEdgeForm(ref CommonParamsBnd inp,
double[] U_IN, double[,] GradU_IN, double V_IN, double[] GradV_IN) {
double eta = PenaltyFactor(inp.jCellIn, -1);
double Acc = 0.0;
for(int d = 0; d < inp.D; d++) { // loop over vector components
// consistency term: -({{ \/u }} [[ v ]])*Normale
// index d: spatial direction
Acc -= (GradU_IN[0, d])*(V_IN) * inp.Normal[d];
// symmetry term: -({{ \/v }} [[ u ]])*Normale
Acc -= (GradV_IN[d])*(U_IN[0]) * inp.Normal[d];
}
// penalty term: eta*[[u]]*[[v]]
Acc += eta*(U_IN[0])*(V_IN);
return Acc;
}
}
/*
// An alternative implementation which derives from the build-in, pre-defined SIP implementation:
public class SipLaplace : BoSSS.Solution.NSECommon.SIPLaplace {
public SipLaplace() : base(2.2, "u") { }
public double PenaltySafety {
get { return base.m_penalty_base; }
set { base.m_penalty_base = value; }
}
protected override bool IsDirichlet(ref CommonParamsBnd inp) {
if(Math.Abs(Math.Abs(inp.X.y) - 1.0) < 1.0e-8)
return false; // Neumann b.c. @ y == +1,-1
else
return true; // Dirichlet b.c. @ x = +1,-1
}
override public double Nu(double[] x, double[] p, int jCell) {
return -1.0;
}
}*/
We consider the following problem: $$ \Delta u = 2,\quad -1<x<1,\quad u(-1)=u(1)=0. $$ The solution is $u_{ex}(x) = 1 - x^2$. Since this is quadratic, we can represent it exactly in a DG space of order 2. As usual, we have to set up a grid, a basis and a right-hand-side.:
var grd1D = Grid1D.LineGrid(GenericBlas.Linspace(-1,1,10));
var DGBasisOn1D = new Basis(grd1D, 2);
var RHS = new SinglePhaseField(DGBasisOn1D, "RHS");
RHS.ProjectField((double x) => 2);
var i_SipLaplace = new SipLaplace();
var Operator_SipLaplace = i_SipLaplace.Operator();
We now want to calculate the residual after inserting the exact solution as well as a wrong solution. The implementation of the exact solution:
var u_ex = new SinglePhaseField(DGBasisOn1D, "$u_{ex}$");
u_ex.ProjectField((double x) => 1.0 - x*x);
The implementation of a spurious, i.e. a wrong solution; we take the exact solution and add random values in each cell:
var u_wrong = new SinglePhaseField(DGBasisOn1D, "$u_{wrong}$");
u_wrong.ProjectField((double x) => 1.0 - x*x);
Random R = new Random();
for(int j = 0; j < grd1D.GridData.Cells.NoOfLocalUpdatedCells; j++){
double ujMean = u_wrong.GetMeanValue(j);
ujMean += R.NextDouble();
u_wrong.SetMeanValue(j, ujMean);
}
Evaluating the Laplace operator using the different solutions:
var Residual = new SinglePhaseField(DGBasisOn1D,"Resi1");
var ResidualNorm = new List<double>();
foreach(var u in new DGField[] {u_ex, u_wrong}) {
Residual.Clear();
Operator_SipLaplace.Evaluate(u, Residual); // evaluate
Residual.Acc(-1.0, RHS);
double ResiNorm = Residual.L2Norm();
ResidualNorm.Add(ResiNorm);
Console.WriteLine("Residual for " + u.Identification + " = " + ResiNorm);
}
Residual for $u_{ex}$ = 1.3694083221572641E-12
Residual for $u_{wrong}$ = 780.3769046218825
/// tests BoSSScmdSilent
Assert.LessOrEqual(ResidualNorm[0], 1e-10);
Assert.GreaterOrEqual(ResidualNorm[1], 1e-1);
If we do not know the exact solution, we have to solve a linear system. Therefore, we not only need to evaluate the operator, but we need its matrix. The Mapping controls which degree-of-freedom of the DG approximation is mapped to which row, resp. column of the matrix.
var Mapping = new UnsetteledCoordinateMapping(DGBasisOn1D);
var Matrix_SipLaplace = Operator_SipLaplace.ComputeMatrix(Mapping,null,Mapping);
Matrix_SipLaplace.NoOfCols
Matrix_SipLaplace.NoOfRows
We see that the matrix has 27 rows and columns.
Matrix_SipLaplace: Use the functions rank and det to analyze the matrix (warning: this can get costly for larger matrices!). Interpret the results:
double rank = Matrix_SipLaplace.rank();
Console.WriteLine("Matrix rank = " + rank);
double det = Matrix_SipLaplace.det();
Console.WriteLine("Determinante = " + det);
Matrix rank = 27 Determinante = 1.405008381452841E+75
So the matrix of the SIP discretization has a unique solution.
/// tests BoSSScmdSilent
Assert.AreEqual(rank, Matrix_SipLaplace.NoOfCols);
Assert.Greater(det, 1.0);
We define a two-dimensional grid:
var grd2D = Grid2D.Cartesian2DGrid(GenericBlas.Linspace(-1,1,21),
GenericBlas.Linspace(-1,1,16));
var DGBasisOn2D = new Basis(grd2D, 5);
var Mapping2D = new UnsetteledCoordinateMapping(DGBasisOn2D);
We are going to choose the PenaltySafety for the SipLaplace from the following list
double[] SFs = new double[]
{0.001, 0.002, 0.01, 0.02, 0.1, 0.2, 1, 2, 10, 20, 100};
and compute the condition number as well as the determinate. We consider the example $$ -\Delta u = \pi^2 (a_x^2 + a_y^2)/4 \cos(a_x \pi x/2) \cos(a_y \pi y/2) \text{ with } (x,y) \in (-1,1)^2 $$ and $u = 0$ on the boundary. The exact solution is $u_{Ex}(x,y) = \cos(a_x \pi x/2) \cos(a_x \pi y/2)$, where $a_x$ and $a_y$ must be odd numbers to comply with homogeneous bounary condition.
double ax = 1.0; // must be an even number to comply with homogeneous Dirichlet boundary condition
double ay = 3.0; // must be an odd number to comply with homogeneous Dirichlet boundary condition
Func<double[], double> exSol =
(X => Math.Cos(X[0]*ax*Math.PI*0.5)*Math.Cos(X[1]*ay*Math.PI*0.5));
Func<double[], double> exRhs =
(X => ((ax.Pow2() + ay.Pow2())/4.0)*Math.PI.Pow2()
*Math.Cos( X[0]*ax*Math.PI*0.5 )*Math.Cos( X[1]*ay*Math.PI*0.5 )); // == - /\ exSol
SinglePhaseField RHS = new SinglePhaseField(DGBasisOn2D, "RHS");
RHS.ProjectField(exRhs);
double[] RHSvec = RHS.CoordinateVector.ToArray();
We check our discretization once more in 2D; the residual should be low, but not exactly (resp. up to $10^{-12}$) since the solution is not polynomial and cannot be fulfilled exactly.
SinglePhaseField u = new SinglePhaseField(DGBasisOn2D,"u");
u.ProjectField(exSol);
var Matrix_SIP_sf = Operator_SipLaplace.ComputeMatrix(Mapping2D,
null,
Mapping2D);
SinglePhaseField Residual = new SinglePhaseField(DGBasisOn2D,"Residual");
Residual.Acc(1.0, RHS);
Matrix_SIP_sf.SpMV(-1.0, u.CoordinateVector, 1.0, Residual.CoordinateVector);
Console.WriteLine("Residual L2 norm: " + Residual.L2Norm());
Residual L2 norm: 0.02104929706815022
We also check that the matrix is symmetric:
var checkMatrix = Matrix_SIP_sf - Matrix_SIP_sf.Transpose();
checkMatrix.InfNorm()
/// tests BoSSScmdSilent
Assert.LessOrEqual(checkMatrix.InfNorm(), 1e-8);
Now, we assemble the matrix of the SIP for different PenaltySafety-factors. We also try to solve the linear system using an iterative method.
As Matlab is called multiple times during this command, it can take some minutes until it is done.
int cnt = 0;
var Results = new List<(double safetyFactor, double condNumber, int NoOfIterations, double L2errror, bool isDefinite)>();
foreach(double sf in SFs) {
cnt++;
i_SipLaplace.PenaltySafety = sf;
var Matrix_SIP_sf = Operator_SipLaplace.ComputeMatrix(
Mapping2D, null, Mapping2D);
double condNo1 = Matrix_SIP_sf.condest();
bool definite = Matrix_SIP_sf.IsDefinite();
/// We solve the system
///
/// Matrix\_SIP\_sf \cdot u = RHS
///
/// using a an iterative solver, the so-called
/// conjugate gradient (CG) method.
/// CG requires a positive definite matrix.
/// The function \code{Solve\_CG} returns the number of iterations.
SinglePhaseField u = new SinglePhaseField(DGBasisOn2D,"u");
u.InitRandom();
int NoOfIter = Matrix_SIP_sf.Solve_CG(u.CoordinateVector, RHSvec);
SinglePhaseField Error = new SinglePhaseField(DGBasisOn2D,"Error");
Error.ProjectField(exSol);
Error.Acc(-1.0, u);
double L2err = u.L2Error(exSol);
Console.WriteLine(sf + "\t" + condNo1.ToString("0.#E-00")
+ "\t" + NoOfIter
+ "\t" + L2err.ToString("0.#E-00")
+ "\t" + definite);
Results.Add((sf, condNo1, NoOfIter, L2err, definite));
}
0.001 8.1E04 12803 3.3E-07 False 0.002 8.1E04 14802 2.8E-07 False 0.01 7.9E04 19510 2.8E-07 False 0.02 7.9E04 24047 3.1E-07 False 0.1 8.5E05 49480 4.1E-06 False 0.2 4.5E04 4759 6.2E-07 False 1 3.5E05 1447 1.7E-06 True 2 7.9E05 2061 2.9E-06 True 10 4.3E06 4512 1.4E-05 True 20 8.6E06 5680 2.5E-05 True 100 4.3E07 12304 1.1E-04 True
/// tests BoSSScmdSilent
foreach(var r in Results) {
if(r.safetyFactor >= 1 && r.safetyFactor <= 20) {
Assert.LessOrEqual(r.condNumber, 1e7); // cond No.
Assert.LessOrEqual(r.NoOfIterations, 7000); // iter
Assert.LessOrEqual(r.L2errror, 1e-4); // L2 err
Assert.IsTrue(r.isDefinite); // definite
}
if(r.Item1 <= 0.1) {
Assert.IsFalse(r.isDefinite); // indefinite
}
}
Plot the number of conjugate gradient iterations versus the PenaltySafety.
var xValues = Results.Select(r => r.safetyFactor).ToArray();
var yValues = Results.Select(r => ((double)(r.NoOfIterations))).ToArray();
var plt = new Plot2Ddata();
plt.AddDataGroup(xValues, yValues);
/// A logarithmic scale is used for the horizontal axis.
plt.LogX = true;
/// Set Format
plt.dataGroups[0].Format = new PlotFormat(lineColor: LineColors.Blue,
pointSize: 2,
dashType: DashTypes.DotDashed,
Style: Styles.LinesPoints,
pointType:PointTypes.OpenCircle);
// Show!
plt.PlotNow()
Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!
We are going to solve the SIP-system for different grid resolutions, comparing an insufficient penalty to a penalty which is large enough.
double[] Resolution = new double[] { 2, 4, 8, 16, 32, 64 };
List<double> L2Error_indef = new List<double>();
List<double> L2Error_posdef = new List<double>();
int cnt = 0;
foreach(int Res in Resolution) {
cnt++;
//var grd2D = Grid2D.UnstructuredTriangleGrid(GenericBlas.Linspace(-1,1,(int)Res + 1),
// GenericBlas.Linspace(-1,1,(int)Res + 1));
var grd2D = Grid2D.Cartesian2DGrid(GenericBlas.Linspace(-1,1,(int)Res + 1),
GenericBlas.Linspace(-1,1,(int)Res + 1));
var gdata2D = new GridData(grd2D);
var DGBasisOn2D = new Basis(gdata2D, 2);
var Mapping2D = new UnsetteledCoordinateMapping(DGBasisOn2D);
SinglePhaseField RHS = new SinglePhaseField(DGBasisOn2D, "RHS");
RHS.ProjectField(exRhs);
SinglePhaseField uEx = new SinglePhaseField(
new Basis(gdata2D, DGBasisOn2D.Degree*2),
"Error");
uEx.ProjectField(exSol);
i_SipLaplace.PenaltySafety = 0.01;
var Matrix_SIP_indef = Operator_SipLaplace.ComputeMatrix(
Mapping2D,null,Mapping2D);
SinglePhaseField u_indef = new SinglePhaseField(DGBasisOn2D,"u_indef");
Matrix_SIP_indef.Solve_Direct(u_indef.CoordinateVector,
RHS.CoordinateVector);
var Error_indef = uEx.CloneAs();
Error_indef.AccLaidBack(-1.0, u_indef);
L2Error_indef.Add(Error_indef.L2Norm());
/// In order to have a positive definite system, we are
/// using PenaltySafety = 2!
i_SipLaplace.PenaltySafety = 2.0;
var Matrix_SIP_posdef = Operator_SipLaplace.ComputeMatrix(
Mapping2D, null, Mapping2D);
SinglePhaseField u_posdef = new SinglePhaseField(DGBasisOn2D,"u_posdef");
Matrix_SIP_posdef.Solve_Direct(u_posdef.CoordinateVector,
RHS.CoordinateVector);
var Error_posdef = uEx.CloneAs();
Error_posdef.AccLaidBack(-1.0, u_posdef);
L2Error_posdef.Add(Error_posdef.L2Norm());
//Tecplot("ConvStudy-" + cnt, uEx, u_posdef, u_indef); // activate this line for plotting!
Console.WriteLine(L2Error_indef.Last().ToString("0.#E-00")
+ "\t" + L2Error_posdef.Last().ToString("0.#E-00"));
}
1.5E01 7.2E-01 4.1E-01 1.9E-01 3.1E-02 1.7E-02 5.8E-03 1.6E-03 8.4E-04 1.7E-04
Solve_Direct fail in fallback seq (#0): Linear Solver: High residual from direct solver ilPSP.LinSolvers.PARDISO.PARDISOSolver.
L2 Norm of RHS: 24.674010999677737
L2 Norm of Solution: 0.9999956185473319
L2 Norm of Residual: 0.001995397603712859
Relative Residual norm: 1.685339351028455E-08
Matrix Inf norm: 118397.37810045174
Solve_Direct fail in fallback seq (#1): Linear Solver: High residual from direct solver ilPSP.LinSolvers.PARDISO.PARDISOSolver.
L2 Norm of RHS: 24.674010999677737
L2 Norm of Solution: 0.9999956237380531
L2 Norm of Residual: 0.0024287222272182667
Relative Residual norm: 2.0513310904213342E-08
Matrix Inf norm: 118397.37810045174
1.1E-04 2.1E-05
The convergence plot should unveil that there is something wrong if the
penalty factor is set too low.
Unfortunately, it does not, so this is some kind of anti-example;
It is in this tutorial anyway to illustrate the difficulties of numerical testing.
Interested readers migth check out the source code and
try to modify the test BoSSS.Application.SipPoisson.Tests.TestProgram.TestOperatorConvergence3D(2)
so that it fails.
The reason why the indefinite matrix still gives a solution convergence is very likely that the solver which is used in BoSSS is also (sometimes) capable of solving singular or close-to-singular systems, i.e. systems without a unique solution. In those cases, it selects a solution with a minimal solution norm. Since BoSSS uses an orthonormal basis the $L^2$ norm of the DG-Field is identical to the $l_2$ norm of the coordinate vector (Parseval's identity). Therefore, the solver by chance adds additional stability which is not part of the (instable) discretization.
While the solution of the indefinite system may look right at the first
glance, we see that we do not obtain grid convergence for
Error_indef.
The error of the positive definite system, Error_posdef, where the penalty is chosen sufficiently large converges with the expected rate.
var plt = new Plot2Ddata();
plt.AddDataGroup("indef mtx", Resolution, L2Error_indef);
plt.AddDataGroup("pos def mtx", Resolution, L2Error_posdef);
/// A double-logarithmic scale is used:
plt.LogX = true;
plt.LogY = true;
/// Set Format
plt.dataGroups[0].Format = new PlotFormat(lineColor: LineColors.Red,
pointSize: 2,
dashType: DashTypes.DotDashed,
Style: Styles.LinesPoints,
pointType:PointTypes.OpenCircle);
plt.dataGroups[1].Format = new PlotFormat("Blue-.o"); // altenatively, using MATLAB-like format strings
plt.dataGroups[1].Format.PointSize = 2;
// Show!
plt.PlotNow()
Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!
Finally, we are going to compute the convergence rate of the SIP discretization.
We compute the slope of the log-log plot:
double dk = Resolution.LogLogRegressionSlope(L2Error_posdef);
dk
/// tests BoSSScmdSilent
Assert.LessOrEqual(dk, -2.9);
An alternative way of identifying problems with the discretization is the investigation of minimal Eigenvalues (in absolute value) and the respective Eigenvectors/Eigensolutions.
This very often unveils stability issues on rather coarse meshes.
var grd2D = Grid2D.Cartesian2DGrid(GenericBlas.Linspace(-1,1,11),
GenericBlas.Linspace(-1,1,11));
var DGBasisOn2D = new Basis(grd2D, 2);
var Mapping2D = new UnsetteledCoordinateMapping(DGBasisOn2D);
i_SipLaplace.PenaltySafety = 0.01;
var Matrix_SIP_indef = Operator_SipLaplace.ComputeMatrix(Mapping2D,null,Mapping2D);
i_SipLaplace.PenaltySafety = 4.0;
var Matrix_SIP_posdef = Operator_SipLaplace.ComputeMatrix(Mapping2D,null,Mapping2D);
For both matrixes (indefinite and positive definite) one can determine the respective Eigenvectors (for the minimal Eigenvalue). Since Eigenvectors represent solution to the matrix for a RHS that is a multiple of the Eigenvector itself, it makes sense to interpret the Eigenvecor as a DG field.
(var lmin_indef, var EvectMin_indef) = Matrix_SIP_indef.MinimalEigen();
DGField dg_EvectMin_indef = new SinglePhaseField(DGBasisOn2D, "MinimaEvect-indef");
dg_EvectMin_indef.CoordinateVector.SetV(EvectMin_indef);
lmin_indef
0 -- Change norm is: 1 1 -- Change norm is: 0.9344394389137047 2 -- Change norm is: 0.005646713765192463 3 -- Change norm is: 0.03148789017389806 4 -- Change norm is: 0.05830203780337603 5 -- Change norm is: 0.10403948087978505 6 -- Change norm is: 0.1857369997615229 7 -- Change norm is: 0.3388216764622676 8 -- Change norm is: 0.6834512341491299 9 -- Change norm is: 1.666508160771979 10 -- Change norm is: 0.6214468477952105 11 -- Change norm is: 0.34292160606455924 12 -- Change norm is: 0.20092883155702643 13 -- Change norm is: 0.11750435072745946 14 -- Change norm is: 0.06780111957340085 15 -- Change norm is: 0.03865592259005921 16 -- Change norm is: 0.021857837751188398 17 -- Change norm is: 0.012296043660706371 18 -- Change norm is: 0.00689607093503553 19 -- Change norm is: 0.0038607865758993024 20 -- Change norm is: 0.0021593176408484876 21 -- Change norm is: 0.0012070154811334051 22 -- Change norm is: 0.0006744843365777928 23 -- Change norm is: 0.00037683738569684733 24 -- Change norm is: 0.00021051985398619957 25 -- Change norm is: 0.00011760020435330092 26 -- Change norm is: 6.56915705841577E-05 27 -- Change norm is: 3.6694730341451044E-05 28 -- Change norm is: 2.049715371626279E-05 29 -- Change norm is: 1.144935629991898E-05 30 -- Change norm is: 6.395393534771174E-06 31 -- Change norm is: 3.5723396553695677E-06 32 -- Change norm is: 1.9954360160836667E-06 33 -- Change norm is: 1.1146092442523889E-06 34 -- Change norm is: 6.22597462206174E-07
(var lmin_posdef, var EvectMin_posdef) = Matrix_SIP_posdef.MinimalEigen();
DGField dg_EvectMin_posdef = new SinglePhaseField(DGBasisOn2D, "MinimaEvect-posdef");
dg_EvectMin_posdef.CoordinateVector.SetV(EvectMin_posdef);
lmin_posdef
0 -- Change norm is: 1 1 -- Change norm is: 0.9112089352281744 2 -- Change norm is: 0.0052857964459833845 3 -- Change norm is: 0.00019371278300347547 4 -- Change norm is: 9.36688068916085E-06 5 -- Change norm is: 5.959769587267434E-07
In many cases, the visualization of Eigenvalues unveils spurious solutions, resp. instabilities which are hidden in the operator discretization. Compare the Visualization of both Eigenvectors, e.g. in Tecplot or Visit:
Tecplot("Eigenvectors", dg_EvectMin_indef, dg_EvectMin_posdef);
Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\Eigenvectors... done.