Within this tutorial the Lax-Friedrichs flux will be implemented as an alternative to the Upwind flux (see chapter DifferentialOperator). For both fluxes a convergence study will be performed and the "experimental order of convergence" (EOC) will be computed. For the implementation of the Lax-Friedrichs flux and the Upwind flux it is recommended to work through chapter DifferentialOperator first as the definition of fluxes is explained there in more detail. For the second part of the tutorial, chapter GridInstantiation has already taught the basics of doing a convergence study.
As an examplary problem, we consider the scalar transport equation
$$ \frac{\delta c}{\delta t} + \nabla \cdot (\vec{u} c) = 0 $$in the domain $\Omega = [-1, 1] \times [-1, 1]$, where the concentration $c = c(x,y,t) \in \mathbb{R}$ is unknown and the velocity field is given by $$ \vec{u} = \begin{pmatrix} y\\-x \end{pmatrix} $$ The analytical solution for the problem above is given by $$ c_\text{Exact}(x,y,t) = \cos(\cos(t) x - \sin(t) y) \quad \text{ for } (x,y) \in \Omega, $$ which can be used to describe the initial and boundary conditions.
#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();
First, we define functions for the given velocity field and the exact solution
public static class MyGlobals {
public static double u(double[] X) => X[1];
public static double v(double[] X) => -X[0];
public static double cExact(double[] X, double t) => Math.Cos(Math.Cos(t)*X[0] - Math.Sin(t)*X[1]);
}
Before implementing the ScalarTransportFlux as done in Tutorial 4, we define both fluxes as functions in terms of the parameters Uin, Uout, the normal vector n and the velocityVector. Recalling the Upwind flux, the corresponding flux function is defined as
using BoSSS.Platform.LinAlg;
Func<double, double, Vector, Vector, double> upwindFlux =
delegate(double Uin, double Uout, Vector n, Vector velocityVector) {
if (velocityVector * n > 0) {
return (velocityVector * Uin) * n;
} else {
return (velocityVector * Uout) * n;
}
};
The second flux we are considering now, is the Lax-Friedrichs flux $$ \hat{f}(c_h^-, c_h^+, \vec{u} \cdot \vec{n}) = \overline{\vec{u} \, c_h} \cdot \vec{n} + \frac{C}{2} \lbrack \lbrack {c_h} \rbrack \rbrack. $$ The constant $C \in \mathbb{R}^+$ has to be sufficiently large in order to guarantee the stability of the numerical scheme. We choose $$ C = \max \vert{\vec{u} \cdot \vec{n}}\vert = 1 $$ for the given problem.
double C = 0.3;
Func<double, double, Vector, Vector, double> laxFriedrichsFlux =
delegate(double Uin, double Uout, Vector n, Vector velocityVector) {
return 0.5 * (Uin + Uout) * velocityVector * n - C * (Uout - Uin);
};
For the implementation of the ScalarTransportFlux we want to generate instances by the definition of the numericalFlux.
Therefore a new constructor is implemented and the implementation of InnerEdgeFlux is rewritten in terms of the numericalFlux.
class ScalarTransportFlux : NonlinearFlux {
private Func<double, double, Vector, Vector, double> numericalFlux;
// Provides instances of this class with a specific flux implementation
public ScalarTransportFlux(
Func<double, double, Vector, Vector, double> numericalFlux) {
this.numericalFlux = numericalFlux;
}
public override IList<string> ArgumentOrdering {
get { return new string[] { "ch" }; }
}
protected override void Flux(
double time, double[] x, double[] U, double[] output) {
output[0] = MyGlobals.u(x) * U[0];
output[1] = MyGlobals.v(x) * U[0];
}
// Makes use of the flux implementation supplied in the constructor
protected override double InnerEdgeFlux(
double time, double[] x, double[] normal,
double[] Uin, double[] Uout, int jEdge) {
Vector n = new Vector(normal);
Vector velocityVector = new Vector(MyGlobals.u(x), MyGlobals.v(x));
return numericalFlux(Uin[0], Uout[0], n, velocityVector);
}
protected override double BorderEdgeFlux(
double time, double[] x, double[] normal, byte EdgeTag,
double[] Uin, int jEdge) {
double[] Uout = new double[] { MyGlobals.cExact(x, time) };
return InnerEdgeFlux(time, x, normal, Uin, Uout, jEdge);
}
}
For both numerical fluxes we define a DifferentialOperator that uses the corresponding flux to compute div(...)
var upwindOperator = new DifferentialOperator(
new string[] { "ch" },
new string[] { "div" },
QuadOrderFunc.NonLinear(2));
upwindOperator.EquationComponents["div"].Add(
new ScalarTransportFlux(upwindFlux));
upwindOperator.Commit();
var laxFriedrichsOperator = new DifferentialOperator(
new string[] { "ch" },
new string[] { "div" },
QuadOrderFunc.NonLinear(2));
laxFriedrichsOperator.EquationComponents["div"].Add(
new ScalarTransportFlux(laxFriedrichsFlux));
laxFriedrichsOperator.Commit();
In the following a convergence study for the discretization error will be performed on grids with $2\times2$, $4\times4$, $8\times8$ and $16\times16$ cells.
The error at $t = \pi /4$ will be investigated for the polynomial degrees $p = 0, \ldots, 3$ of the ansatzfunctions.
The study will be done for both numerical fluxes.
We start by defining a series of nested grids for the convergence study
int[] numbersOfCells = new int[] { 2, 4, 8, 16 };
GridData[] grids = new GridData[numbersOfCells.Length];
for (int i = 0; i < numbersOfCells.Length; i++) {
double[] nodes = GenericBlas.Linspace(-1.0, 1.0, numbersOfCells[i] + 1);
GridCommons grid = Grid2D.Cartesian2DGrid(nodes, nodes);
grids[i] = new GridData(grid);
}
Then, a SinglePhaseField is defined for each polynomial degree on each grid. The initial value is projected on each field.
int[] dgDegrees = new int[] { 0, 1, 2, 3 };
SinglePhaseField[,] fields =
new SinglePhaseField[dgDegrees.Length, numbersOfCells.Length];
for (int i = 0; i < dgDegrees.Length; i++) {
for (int j = 0; j < numbersOfCells.Length; j++) {
Basis basis = new Basis(grids[j], dgDegrees[i]);
fields[i, j] = new SinglePhaseField(
basis, "ch_" + dgDegrees[i] + "_" + grids[j]);
fields[i, j].ProjectField(X => MyGlobals.cExact(X, 0.0));
}
}
Since the error at a fourth revolution is considered, the initial concentration has to be simulated until $\textbf{endTime} = \pi / 4$. For the timestepping we are using the classical fourth order Runge-Kutta scheme.
The timestep size should be chosen small enough in order to reduce the temporal error, so that the spatial error dominates for the convergence study.
using BoSSS.Solution.Timestepping;
double endTime = Math.PI / 4.0;
int numberOfTimesteps = 100;
double dt = endTime / numberOfTimesteps;
A MultidimensionalArray is a special double array that enables shallow resizing and data extraction. For example, sub-arrays can be extracted without needing to copy the entries of the MultidimensionalArray. BoSSS makes extensive use of this class because this is crucial for the performance. Here, we create a three-dimensional array, where the first index corresponds to the flux function, the second to the DG degree and the third to the number of cells.
MultidimensionalArray errors =
MultidimensionalArray.Create(2, dgDegrees.Length, grids.Length);
for (int i = 0; i < dgDegrees.Length; i++) {
for (int j = 0; j < numbersOfCells.Length; j++) {
// Clones an object and casts it to the type of the original object
// at the same time.
SinglePhaseField c = fields[i, j].CloneAs();
// Instantiate the timestepper (classical Runge-Kutta scheme)
// for the evolution of the concentration c with Upwind flux
RungeKutta timeStepper = new RungeKutta(
RungeKutta.RungeKuttaSchemes.RungeKutta1901,
upwindOperator, c);
// Integrate in time for each timestep
for (int k = 0; k < numberOfTimesteps; k++) {
timeStepper.Perform(dt);
}
// Compute the error w.r.t. analytical solution
errors[0, i, j] = c.L2Error(X => MyGlobals.cExact(X, endTime));
// Simulate with the Lax-Friedrichs flux
c = fields[i, j].CloneAs();
timeStepper = new RungeKutta(
RungeKutta.RungeKuttaSchemes.RungeKutta1901,
laxFriedrichsOperator, c);
for (int k = 0; k < numberOfTimesteps; k++) {
timeStepper.Perform(dt);
}
errors[1, i, j] = c.L2Error(X => MyGlobals.cExact(X, endTime));
}
}
For the convergence plots the error is plotted over the mesh size, where the coarsest grid size is defined as size = 1.
var sizes = numbersOfCells.Select(s => 2.0 / s).ToArray();
First, we are looking at the convergence plot of the Upwind flux. Therefore, an instance of Gnuplot has to be generated.
var gpUpwind = new Gnuplot();
gpUpwind.SetYRange(Math.Pow(10,-7), Math.Pow(10,0));
Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe
for (int i = 0; i < dgDegrees.Length; i++) {
/// Here, we use the shallow extraction features of
/// \code{MultidimensionalArray} mentioned above.
/// The command \code{ExtractSubArrayShallow(0, i, -1)}
/// is roughly equivalent to writing "errors[0, i, :]" in Matlab
var errorsForDegree = errors.ExtractSubArrayShallow(0, i, -1).To1DArray();
// some formatting of the plot data for the convergence study
PlotFormat format = new PlotFormat(lineColor: (LineColors)(i+1),
Style: Styles.LinesPoints,
pointType: PointTypes.Diamond,
pointSize: 2.0);
// Create log-log plot mesh size vs. error for dgDegrees[i]
gpUpwind.PlotLogXLogY(sizes, errorsForDegree,
"Upwind, order " + dgDegrees[i], format);
}
gpUpwind.PlotNow() // perform the plotting
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.!
The convergence plot is also done for the Lax-Friedrichs flux.
var gpLaxFriedrichs = new Gnuplot();
gpLaxFriedrichs.SetYRange(Math.Pow(10,-7), Math.Pow(10,0));
Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe
for (int i = 0; i < dgDegrees.Length; i++) {
var errorsForDegree = errors.ExtractSubArrayShallow(1, i, -1).To1DArray();
PlotFormat format = new PlotFormat(lineColor: (LineColors)(i+1),
Style: Styles.LinesPoints,
pointType: PointTypes.Diamond,
pointSize: 2.0);
gpLaxFriedrichs.PlotLogXLogY(sizes, errorsForDegree,
"LaxFriedrichs, order " + dgDegrees[i],
format);
}
gpLaxFriedrichs.PlotNow()
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.!
For the determination of the experimental order of convergence (EOC) the linear regression is used to compute the slope of the line in the log-log plot. In the following, the ordinary least squares method is implemented to estimate the regression coefficient of the slope for given sets of x- and y-values.
Func<double[], double[], double> slope =
delegate(double[] xValues, double[] yValues) {
if (xValues.Length != yValues.Length) {
throw new ArgumentException();
}
xValues = xValues.Select(s => Math.Log10(s)).ToArray();
yValues = yValues.Select(s => Math.Log10(s)).ToArray();
double xAverage = xValues.Sum() / xValues.Length;
double yAverage = yValues.Sum() / yValues.Length;
double v1 = 0.0;
double v2 = 0.0;
// Computation of the regression coefficient for the slope
for (int i = 0; i < yValues.Length; i++) {
v1 += (xValues[i] - xAverage) * (yValues[i] - yAverage);
v2 += Math.Pow(xValues[i] - xAverage, 2);
}
return v1 / v2;
};
Verification that the slope is computed correctly
Math.Abs(slope(sizes, new double[] { 64.0, 16.0, 4.0, 1.0 }) - 2.0) < 1e-14
The experimental orders of convergence (EOC) are computed for both fluxes for each polynomial degree of the ansatzfunctions
for (int i = 0; i < dgDegrees.Length; i++) {
Console.WriteLine();
var errorsForDegree = errors.ExtractSubArrayShallow(0, i, -1).To1DArray();
Console.WriteLine(
"Upwind flux, order {0}:\t\t {1:F2}", i,
slope(sizes, errorsForDegree));
errorsForDegree = errors.ExtractSubArrayShallow(1, i, -1).To1DArray();
Console.WriteLine(
"Lax-Friedrichs flux, order {0}:\t {1:F2}", i,
slope(sizes, errorsForDegree));
}
Upwind flux, order 0: 0.79 Lax-Friedrichs flux, order 0: 0.73 Upwind flux, order 1: 1.98 Lax-Friedrichs flux, order 1: 2.00 Upwind flux, order 2: 2.96 Lax-Friedrichs flux, order 2: 2.92 Upwind flux, order 3: 3.95 Lax-Friedrichs flux, order 3: 3.97