- using the
`XDGTimestepper`

to perform a high-order, implicit time integration on the Heat equation - performing an ad-hoc time convergence study within notebook. We point out an
**important difference**regarding the error computed against the exact solution vs. the error computed against the numerical solution with the finest timestep (aka. "experimental error")

- implementation of the Symmetrical Interior Penalty (SIP) form for the Laplace operator, chapter
**SIP**

We are investigating the instationary Heat equation $$ \partial_t u - \Delta u = 0 \text{ in } (x,y) \in (-1,1)^2, \ \ \ u = 0 \text{ on } \partial (-1,1)^2 $$ with the initial value $$ u_0(x,y) = \cos(x \pi / 2) \cos(y \pi / 2), $$ and the exact solution $$ u_{\text{ex}}(t,x,y) = ( -\pi^2 / 2 ) u_0(x,y). $$

First, we initialize the new worksheet; Note:

- This tutorial can be found in the source code repository as as
`HeatEquationSolver.ipynb`

. One can directly load this into Jupyter to interactively work with the following code examples. **In the following line, the reference to**. You must either set`BoSSSpad.dll`

is required`#r "BoSSSpad.dll"`

to something which is appropirate for your computer (e.g.`C:\Program Files (x86)\FDY\BoSSS\bin\Release\net6.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`

).

In [1]:

```
#r "BoSSSpad.dll"
using System;
using System.Collections.Generic;
using System.Linq;
using ilPSP;
using ilPSP.Utils;
using BoSSS.Platform;
using BoSSS.Platform.LinAlg;
using BoSSS.Foundation;
using BoSSS.Foundation.XDG;
using BoSSS.Foundation.Grid;
using BoSSS.Foundation.Grid.Classic;
using BoSSS.Foundation.Grid.RefElements;
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.AdvancedSolvers;
using BoSSS.Solution.Gnuplot;
using BoSSS.Application.BoSSSpad;
using BoSSS.Application.XNSE_Solver;
using BoSSS.Application.XNSFE_Solver;
using static BoSSS.Application.BoSSSpad.BoSSSshell;
Init();
```

In [2]:

```
using BoSSS.Solution.XdgTimestepping;
```

Here, 1 2D, $16 \times 16$ mesh is used:

In [3]:

```
var grid = Grid2D.Cartesian2DGrid(GenericBlas.Linspace(-1,1,17), GenericBlas.Linspace(-1,1,17));
```

Regarding the spatial part of the operator, we re-use the SIP-implementation
which is already available in the `BoSSS.Solution.NSECommon`

library.
For this implementation, whe have to specify the diffusion coefficient (here to -1.0)
as well as the Location of the Dirichlet boundary.

In [4]:

```
class Laplace : BoSSS.Solution.NSECommon.SIPLaplace {
public Laplace() : base(1.4, "u") { } // override the constuctor;
// The factor 1.4 is a constant "safty factor" multiper for the penalty parameter.
// The base implementation takes care about all other penalty factor dependencies,
// i.e. dependence of the penatly factor on local cell size, DG polynomial degree and cell shape.
// The second argument, "u", is the name of the variable for the `ArgumentOrdering`.
// Specifies whetherss a specific point (`inp.X`) is either a Dirichlet or a Neumann boundary;
// Since in our case ther, the entire boundary should be Diriclet, we always return true.
protected override bool IsDirichlet(ref CommonParamsBnd inp) {
return true;
}
// diffusion coefficient
override public double Nu(double[] x, double[] p, int jCell) {
return -1.0;
}
}
```

We compose the differential operator from the previously created
`Laplace`

implementation.
The `Laplace`

is also called an "equation component".
Individual equation components can be combined int a differential/spatial operator.
This "componentization" allows to later re-use the components in different operators.

In [5]:

```
var Op = new DifferentialOperator(1, 0, 1, QuadOrderFunc.Linear(), "u", "R1");
Op.EquationComponents["R1"].Add(new Laplace()); // adding
Op.TemporalOperator = new ConstantTemporalOperator(Op);
Op.IsLinear = true;
Op.Commit();
```

**Note on historical development and inconsistent naming conventions:**
Until 2020-API changes, (API-layer 2, aka. `BoSSS.Foundation`

) BoSSS concentrated on spatial operators.
Thus, differential operators were called "spatial operators".
Time discretization was taken car of in higher levels and/or left to the user.
Later, as an attempt to unify the temporal integration, the `ITemporalOperator`

was created,
which was added to the spatial operator, rendering its naming inconsistent.
At some later point (mid 2023), the spatial operator was renamed to "differential operator", but the oly name
might still appear in multiple places, especially in the documentation.

In [6]:

```
SinglePhaseField TimeIntegrate(SinglePhaseField u0, double EndTime, int NoOfTimeteps) {
var u1 = u0.CloneAs();
var Timestepper = new XdgTimestepping(Op, new DGField[] { u1 }, new DGField[] { new SinglePhaseField(u0.Basis)},
TimeSteppingScheme.SDIRK_54);
double dt = EndTime/NoOfTimeteps;
for(int i = 0; i < NoOfTimeteps; i++)
Timestepper.Solve(dt*i, dt);
return u1;
}
```

Specify the initial value; For the initial field we also specify the DG polynomial degree 4, which also sets the DG polynomial degree for the numerical solution.

In [7]:

```
var u0 = new SinglePhaseField(new Basis(grid, 4), "u");
u0.ProjectField((double[] X) => Math.Cos(X[0]*Math.PI*0.5)*Math.Cos(X[1]*Math.PI*0.5));
```

In [8]:

```
Tecplot("u0", u0);
```

In [9]:

```
//for(int i = 0; i < 10; i++) {
// Console.WriteLine(" -------------- timestep: " + (i + 1));
// var u1 = TimeIntegrate(u0, 0.1, 1);
// Tecplot("u" + (i+1), u1);
// u0 = u1;
//}
```

In [10]:

```
var NoOfTimestepS = new int[] { 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000/*, 2000, 5000*/ };
var solutions = new List<SinglePhaseField>();
for(int i = 0; i < NoOfTimestepS.Length; i++) {
Console.WriteLine(" -------------- computing no of timeteps: " + NoOfTimestepS[i]);
var u1 = TimeIntegrate(u0, 0.1, NoOfTimestepS[i]);
//Tecplot("u1." + i, u1);
solutions.Add(u1);
}
```

In [11]:

```
for(int i = 0; i < solutions.Count; i++)
Tecplot("u1." + i, solutions[i]);
```

Computing the error against the exact solution:

In [12]:

```
double uExEq(double[] X) {
double Lambda = -Math.PI.Pow2()*0.5;
double tend = 0.1;
double u0 = Math.Cos(X[0]*Math.PI*0.5)*Math.Cos(X[1]*Math.PI*0.5);
return Math.Exp(Lambda*tend)*u0;
}
```

In [13]:

```
double[] ErrEx = new double[solutions.Count];
for(int i = 0; i < ErrEx.Length; i++) {
ErrEx[i] = solutions[i].L2Error(uExEq);
}
```

In [14]:

```
ErrEx
```

Out[14]:

index | value |
---|---|

0 | 1.5494514264296767E-05 |

1 | 9.55948408019166E-07 |

2 | 5.799024708378291E-08 |

3 | 5.26873013066484E-08 |

4 | 5.266573529027633E-08 |

5 | 5.2665658359614885E-08 |

6 | 5.2665658506708205E-08 |

7 | 5.266565850813002E-08 |

8 | 5.266565850719167E-08 |

9 | 5.266565851339356E-08 |

Computing the error against the finest solution:

In [15]:

```
var uFinest = solutions.Last();
```

In [16]:

```
var Errs = solutions.Take(solutions.Count - 1).Select(u => u.L2Error(uFinest)).ToArray();
```

In [17]:

```
Errs
```

Out[17]:

index | value |
---|---|

0 | 1.549442980822616E-05 |

1 | 9.545014425038409E-07 |

2 | 2.4278077172585493E-08 |

3 | 1.5146627871626243E-09 |

4 | 9.468514161427916E-11 |

5 | 2.4732738808243437E-12 |

6 | 1.841265450576875E-13 |

7 | 1.0779650298681105E-13 |

8 | 1.4613824767462205E-13 |

Computing the error between each solution and the finer one:

In [18]:

```
double[] diffs = new double[solutions.Count - 1];
for(int i = 0; i < diffs.Length; i++) {
diffs[i] = solutions[i].L2Error(solutions[i + 1]);
}
```

In [19]:

```
diffs
```

Out[19]:

index | value |
---|---|

0 | 1.4539928365852596E-05 |

1 | 9.302233653390543E-07 |

2 | 2.2763414385634946E-08 |

3 | 1.419977648613582E-09 |

4 | 9.221200246126197E-11 |

5 | 2.2901664486873485E-12 |

6 | 7.904200985429604E-14 |

7 | 4.00552668286762E-14 |

8 | 1.4613824767462205E-13 |

Computing the experimental error of the DG discretization with DG polynomial degree $k$ against degree $k-1$.

In [20]:

```
double ErrHi(SinglePhaseField u) {
SinglePhaseField uLo = new SinglePhaseField(new Basis(u.Basis.GridDat, u.Basis.Degree - 1));
uLo.AccLaidBack(1.0, u);
var uHi = u.CloneAs();
uHi.AccLaidBack(-1.0, uLo);
return uHi.L2Norm();
}
```

In [21]:

```
double[] HiDeg = new double[solutions.Count];
for(int i = 0; i < HiDeg.Length; i++) {
HiDeg[i] = ErrHi(solutions[i]);
}
```

In [22]:

```
HiDeg
```

Out[22]:

index | value |
---|---|

0 | 1.7335106741759476E-06 |

1 | 1.7335276040134437E-06 |

2 | 1.7335271671465917E-06 |

3 | 1.73352709675771E-06 |

4 | 1.7335270927023246E-06 |

5 | 1.7335270924588835E-06 |

6 | 1.733527092442236E-06 |

7 | 1.73352709245955E-06 |

8 | 1.733527092459966E-06 |

9 | 1.733527092461772E-06 |

In [23]:

```
var gp = new Gnuplot();
double[] ts = NoOfTimestepS.Take(NoOfTimestepS.Length).Select(x => (double)x).ToArray();
double[] ts1 = ts.Take(NoOfTimestepS.Length - 1).ToArray();
gp.PlotLogXLogY(ts1, Errs, title:"Exp Error", format:new PlotFormat("-xk"));
gp.PlotLogXLogY(ts1, diffs, title:"Diffs", format:new PlotFormat("-ob"));
gp.PlotLogXLogY(ts, ErrEx, title:"Exact Error", format:new PlotFormat("-+r"));
gp.PlotLogXLogY(ts, HiDeg, title:"High Degree Norm", format:new PlotFormat(":ok"));
gp.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.!

Out[23]:

Confirming the order of the time integration:

In [24]:

```
double slope = ts.Take(6).LogLogRegressionSlope(Errs.Take(6));
slope
```

Out[24]:

-4.001274177790234

In [25]:

```
NUnit.Framework.Assert.LessOrEqual(slope, -3.5, "time convergence of the DIRK_54 method is not reached.")
```

Note that the lower limit (lower plateau) of the error plots looks different, depending

- whether we compare against the finest solution (aka. experimental error, "Exp Error" in plot),
- or we compare against the "real" exact solution ("Exact Error" in plot).

Fro finer time-steps, the experimental error is reducing further and further,
the lower plateau is the floating-point accuracy limit of the `double`

data type.
The exact error hits its lower plateau at a much higher level,
where the accuracy limit by the DG polynomial degree and the spatial accuracy is reached.
I.e., using a finer resolution and/or higher DG degree will bring the lower limit of the red curve further down.
The blue and red curve are not affected by the spatial resolution - one would have to
use increased floating-point accuracy (not feasible for BoSSS) to get this limit lower.