Heat Equation Solver using most recent API

What's new?

• 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")

Prerequisites

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

Problem statement

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). $$

Solution within the BoSSS framework

First, we initialize the new worksheet; Note:

1. 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.
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\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).

#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();

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

using BoSSS.Solution.XdgTimestepping;

Creation of Grid and Differential Operator

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

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.

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.

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.

Performing Temporal Integration

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.

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));

Tecplot("u0", u0);

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u0...

done.

//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;
//}

Performing a temporal convergence study

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);
}

   --------------  computing no of timeteps: 1

Warning: Multigrid sequence is requested, but not setup for this mesh; performing a zero-aggregation with one level; consider configuring multigrid using the `AggregationGridData[] MultigridSequence(...)` method.

   --------------  computing no of timeteps: 2

   --------------  computing no of timeteps: 5

   --------------  computing no of timeteps: 10

   --------------  computing no of timeteps: 20

   --------------  computing no of timeteps: 50

   --------------  computing no of timeteps: 100

   --------------  computing no of timeteps: 200

   --------------  computing no of timeteps: 500

   --------------  computing no of timeteps: 1000

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

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.0...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.1...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.2...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.3...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.4...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.5...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.6...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.7...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.8...

done.

Supersampling 4 requested, but limiting to 3 because higher values would very likely exceed this computers memory.

Note: Higher supersampling values are supported by external plot application, or by using e.g. the Tecplot class directly.

Writing output file C:\Gitlab-Runner\builds\PYSUsL3w\1\kummer\bosss-public\src\L4-application\PublicTestRunner\bin\Release\net6.0\u1.9...

done.

Computing the error against the exact solution:

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;
}

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

ErrEx

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:

var uFinest = solutions.Last();

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

Errs

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:

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

diffs

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$.

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();
}

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

HiDeg

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

Plotting different errors

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.!

Confirming the order of the time integration:

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

-4.001274177790234

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.