What's new

• structure of stiffness matrix
• stability regions of explicit time integration schemes
• determination and significance of most critical Eigenvalue
• calculation of stable time step $dt_0$
• influence of grid resolution and polynomial degree on time step size

Prerequisites

• implementation of fluxes, chapter DifferentialOperator and NumFlux
• explicit time integration, chapter DifferentialOperator
• basics of convergence studies, chapters GridInstantiation and NumFlux
• grid generation, chapter GridInstantiation

In this tutorial, you will learn the basics of time discretization. Using the Upwind, Central and two Lax-Friedrichs fluxes the tutorial first visualizes the structure of the stiffness matrix and its Eigenvalues. Then the stability regions of explicit time integration schemes, e.g. the explicit Euler method, are considered. Using the fourth order Runge-Kutta scheme the Eigenvalues are studied for different polynomial orders and grid resolutions in order to find the most critical Eigenvalue. Finally, the influence of grid resolution and polynomial degree on the time step size is examined and the stable time step $dt_0$ determined in order to obtain a stable and converging solution.

Problem statement

To keep it simple, the one-dimensional scalar transport equation $$ \frac{\delta g}{\delta t} + \nabla \cdot \left(uc\right) = 0 $$ will be used, where $g=g(x,t) \in \mathbb{R}$ is the unknown concentration and $u=1$ the velocity field in a periodic domain~$\Omega = \left[-\pi,\pi\right]$.

The linear PDE can be written in its semi-discrete form $$ \frac{\delta \vec{g}}{\delta t} + M \, \vec{g} + **\vec b** = 0 $$ with the unknown DG coefficients~$\vec{g}$, the system or stiffness matrix~$M$ and the affine part~$\vec{b}$. The affine part is given by source terms and boundary conditions. In this case, $\vec{b}= 0$.

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.

Definition of the operator matrix} \label{sec:scalarConvection_matrix

using BoSSS.Platform.LinAlg;

using ilPSP.LinSolvers;

The function~u is defined for later usage.

static Func<double[], double> u = (X => 1.0);

The function~GetOperatorMatrix constructs the operator matrix~$M$ corresponding to the semi-discrete ODE system.

Func<int, int, LinearFlux, MsrMatrix> GetOperatorMatrix =       
    delegate(int numberOfCells, int dgDegree, LinearFlux flux) {       
        double[] nodes = GenericBlas.Linspace(-Math.PI, Math.PI,       
            numberOfCells + 1);       
        GridData gridData = new GridData(Grid1D.LineGrid(       
            nodes, periodic: true));       
 
        Basis basis = new Basis(gridData, dgDegree);       
        SinglePhaseField g = new SinglePhaseField(basis);       
 
        var op = new DifferentialOperator(       
            new string[] { "g" },         // domain variable       
            new string[] { "div" },     // co-domain variable      
            QuadOrderFunc.Linear());       
        op.EquationComponents["div"].Add(flux);       
 
        op.Commit();       
 
        /// The \code{MsrMatrix} is a sparse matrix, where \emph{MSR} stands for \emph{Mutable Sparse Row}. That is, the
        /// matrix can be changed, and the entries are stored in a compressed format where we only have
        /// a small number of column entries per row. This is crucial to be able to handle larger
        /// systems and is required by most linear solvers.
        MsrMatrix operatorMatrix = new MsrMatrix(g.Mapping);       
        double[] affineOffset = new double[g.Mapping.LocalLength];       
        /// Compute the matrix~\code{operatorMatrix} and the affine part~\code{affineOffset} such that an operator
        /// evaluation~$op(g)$ can be expressed as \code{operatorMatrix} $\times$ \code{g} + \code{affineOffset}.
        op.ComputeMatrix(g.Mapping, null, g.Mapping, operatorMatrix,       
            affineOffset, false);       
 
        /// Since we do not enforce any boundary conditions, \code{affineOffset} should be zero.
        if (affineOffset.Any(d => d != 0.0)) {       
            throw new Exception(       
                "We have only periodic boundary conditions."       
                + " Affine part should be zero!");       
        }       
 
        return operatorMatrix;       
    };

Definition of numerical fluxes

As in the previous tutorials, the Upwind flux is used.

Func<double, double, double, double, double> upwindFlux =       
    delegate (double Uin, double Uout, double n, double velocity) {       
        if (velocity * n > 0) {       
            return (velocity * Uin) * n;       
        } else {       
            return (velocity * Uout) * n;       
        }       
    };

Moreover, the Lax-Friedrichs flux with a constant C (here, different values are used) is implemented.

If C equals zero, the Lax-Friedrichs flux is equivalent to the central flux.

Func<double, Func<double, double, double, double, double>>       
  laxFriedrichsFlux =       
    C => delegate (double Uin, double Uout, double n, double velocity) {       
             return 0.5 * (Uin + Uout) * velocity * n - C * (Uout - Uin);       
         };

In general, the flux implementation is similar to the one in the previous tutorials, but in this case, it is one-dimensional and uses LinearFlux as a base class. This class helps to create the entries of the operator matrix discussed above.

class LinearTransportFlux : LinearFlux {       
 
    private Func<double, double, double, double, double> numericalFlux;       
 
    public LinearTransportFlux(Func<double, double, double, double,       
        double> numericalFlux) {       
            this.numericalFlux = numericalFlux;       
    }       
 
    public override IList<string> ArgumentOrdering {       
        get {       
            return new string[] { "g" };       
        }       
    }       
 
    /// The volume term is very similar to what we have done before.
    protected override void Flux(ref CommonParamsVol inp, double[] U,       
        double[] output) {       
        output[0] = u(inp.Xglobal) * U[0];       
    }       
 
    /// The \code{InnerEdgeFlux} implementation is the same as in the last tutorials, with a slightly different signature.
    protected override double InnerEdgeFlux(ref CommonParams inp,       
        double[] Uin, double[] Uout) {       
            return numericalFlux(Uin[0], Uout[0], inp.Normal[0], u(inp.X));       
    }       
 
    /// We are working on periodic grids only, so the \code{BorderEdgeFlux} is irrelevant.
    protected override double BorderEdgeFlux(ref CommonParamsBnd inp,       
        double[] Uin) {       
            throw new NotImplementedException("Should never be called!");       
    }       
}

Definition of operators and construction of system matrices

A dictionary is created to access the different fluxes easily.

var fluxes = new Dictionary<string, Func<double, double, double, double,       
    double>>() {       
        { "Upwind",                   upwindFlux },       
        { "Central",                  laxFriedrichsFlux(0.0) },       
        { "Lax-Friedrichs (C = 0.1)", laxFriedrichsFlux(0.1) },       
        { "Lax-Friedrichs (C = 0.3)", laxFriedrichsFlux(0.3) }       
};

Some initial configurations (polynomial degree and numberOfCells) are done, as well as defining the system matrices for the above-defined fluxes.

int degree = 2;       
int numberOfCells = 10;       
MsrMatrix[] matrices = new MsrMatrix[fluxes.Count];       
for (int i = 0; i < fluxes.Count; i++) {       
    LinearFlux flux = new LinearTransportFlux(fluxes.ElementAt(i).Value);       
    matrices[i] = GetOperatorMatrix(numberOfCells, degree, flux);       
}

Investigation of operator matrices

The style of the plots is defined.

PlotFormat format = new PlotFormat(       
    Style: Styles.Points,       
    pointType: PointTypes.Asterisk,       
    pointSize: 0.1);

SetMultiplot enables multiple plots per graphic.

Gnuplot gp1 = new Gnuplot(baseLineFormat: format);       
gp1.SetMultiplot(2, fluxes.Count / 2);       
gp1.SetXLabel("Matrix column");       
gp1.SetYLabel("Matrix row");       
// Invert y-axis       
gp1.SetYRange(numberOfCells * (degree + 1), 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 < fluxes.Count; i++) {       
    // Select the active sub-plot       
    gp1.SetSubPlot(i / 2, i % 2, title: fluxes.ElementAt(i).Key);       
    /// \code{SetSubPlot} visualizes the non-zero entries of the system matrix. The block size defines the number
    /// of entries per cell. In our case, we have \code{degree} + 1 polynomials per cell such that the
    /// matrix-block associated with a single cell is a (\code{degree} + 1) $\times$ (\code{degree} + 1) sub-matrix.
    gp1.PlotMatrixStructure(matrices[i].ToFullMatrixOnProc0(),       
    blockSize: degree + 1);       
    // Finalize the sub-plot       
    gp1.WriteDeferredPlotCommands();       
}

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

Now, we can see purple and green blocks in the matrix plots. The purple blocks are the diagonal blocks that couple the DOF inside in a cell.

The off-diagonal green blocks define the coupling between neighboring cells.

Moreover, the differences in the matrix structure for the different fluxes,

especially when comparing the upwind and the central flux to the Lax-Friedrichs variants, is investigate

• Upwind flux: Using the upwind flux, the coupling between cells is asymmetric because the flow

     direction is taken into account. As a result, the DOF for a given cell depend on
     the DOF of its left neighbor (green blocks above the diagonal), but not on the DOF of
     its right neighbor (no green blocks below the diagonal in case of Upwind flux).

• Central flux: It can be shown (not proven in this tutorial) that the first entry of each

     diagonal block vanishes. That is because of the symmetry of the basis functions, since
     the corresponding term equals zero as well. Please note, that this is a special case for the
     basis function used here.

Analysis of the spectrum of common explicit ODE integrators

The discrete spectrum of a linear operator in a periodic domain is given by the eigenvalues of the corresponding system matrix. Here, we will use Matlab to analyze the spectrum for the different fluxes. For the stability analysis, we have to use the standard form:

$$\frac{\delta \vec u}{\delta t} = -M \, \vec u $$

using ilPSP.Connectors.Matlab;

The BatchmodeConnector initializes an interface to Matlab:

BatchmodeConnector connector = new BatchmodeConnector();

We will use the array~eigenValues to store the eigenvalues of the matrix for each flux. The first index corresponds to the flux, the second index to the eigenvalue and the third index to the part of the eigenvalue (0: real part; 1: imaginary part).

MultidimensionalArray eigenvalues = MultidimensionalArray.Create(       
    fluxes.Count, (int)matrices[0].NoOfCols, 2);

for (int i = 0; i < fluxes.Count; i++) {       
    // Transfer sparse matrix to Matlab and name the Matlab       
    // variables "M0", "M1", ...       
    connector.PutSparseMatrix(matrices[i], "M" + i);       
    // Compute \emph{all} eigenvalues of the individual matrices       
    // and sort them by magnitude       
    connector.Cmd("eigenvalues{0} = sort(eig(full(M{0})))", i);       
 
    // Separate real and imaginary part; negating the real part       
    // accounts for the minus sign in the standard form discussed above       
    connector.Cmd("result{0} = [-real(eigenvalues{0}) " +       
        "imag(eigenvalues{0})]", i);       
 
    // Retrieve results       
    connector.GetMatrix(eigenvalues.ExtractSubArrayShallow(i, -1, -1),       
        "result" + i);       
}

Calling Matlab takes some time on most machines, even though the individual calculations are quite fast (suppressing the output). This is due to the slow startup of Matlab...

connector.Execute(PrintOutput: false);

Now, the eigenvalues are plotted in the complex plane.

Plot2Ddata[,] specMulPlot = new Plot2Ddata[2,2];

for (int i = 0; i < fluxes.Count; i++) {       
    string key   = fluxes.ElementAt(i).Key;
    var specPlot = new Plot2Ddata();
    specMulPlot[i / 2, i % 2] = specPlot;
    specPlot.Title     = key;  
    specPlot.XrangeMin = -15;
    specPlot.XrangeMax = +1;
    specPlot.YrangeMin = -15;
    specPlot.YrangeMax = +15;
    specPlot.Xlabel    = "Real";
    specPlot.Ylabel    = "Imaginary";
 
    var p = specPlot.AddDataGroup(
        eigenvalues.ExtractSubArrayShallow(i, -1, 0).To1DArray(),       
        eigenvalues.ExtractSubArrayShallow(i, -1, 1).To1DArray());       
    p.Format.Style     = Styles.Points;
    p.Format.PointSize = 2;
}

specMulPlot.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 will compare the spectra of the different operator matrices. The spectrum for the central flux case does not contain any negative real parts. Consequently, the solution of the ODE system is unstable in a sense that disturbances (e.g. due to round-off or projection errors) will never be damped. Next, let us have a closer look at the influence of the spectrum of the system matrix on the maximum admissible time step size. The stability can be analyzed by using the so-called stability function~$f(z)$ of an ODE solver. That is, the system is stable if all eigenvalues of $dt \, M$ are located within the stability region $|f(z)| < 1$ of the solver. In the following, we will nvestigate the stability for the case of the Upwind flux.,

using System.Numerics;

The stability function of Runge-Kutta schemes for the orders 1 to 4 are stored in a Dictionary.

var stabilityFunctions = new Dictionary<string, Func<Complex, Complex>>() {       
    { "Euler (order 1)",       z => 1 + z },       
    { "Heun (order 2)",        z => 1.0 + z + 0.5*z*z },       
    { "Kutta (order 3)",       z => 1.0 + z + 0.5*z*z + 1.0/6.0*z*z*z },       
    { "Runge-Kutta (order 4)", z => 1.0 + z + 0.5*z*z + 1.0/6.0*z*z*z +        
        1.0/24.0*z*z*z*z }       
};

We create some sampling nodes for the following plots.

double[] xNodes = GenericBlas.Linspace(-4.0, 1.0, 50).ToArray();       
double[] yNodes = GenericBlas.Linspace(-4.0, 4.0, 100).ToArray();

We plot the boundary of the stability regions for the different schemes and orders, respectively.

Gnuplot gp3 = new Gnuplot(baseLineFormat: format);       
gp3.SetMultiplot(2, stabilityFunctions.Count / 2);       
gp3.SetXLabel("Real");       
gp3.SetYLabel("Imaginary");       
gp3.Cmd("set grid");       
gp3.Cmd("set size square");       
gp3.SetXRange(-4.0, 4.0);       
gp3.SetYRange(-4.0, 4.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 < stabilityFunctions.Count; i++) {       
    string key = stabilityFunctions.ElementAt(i).Key;       
    gp3.SetSubPlot(i / 2, i % 2, title: key);       
 
    /// We plot the iso-contour where $|\code{stabilityFunction}| = 1.0$.
    gp3.PlotContour(       
        xNodes,       
        yNodes,       
        (x, y) => stabilityFunctions[key](new Complex(x, y)).Magnitude,       
        new double[] { 1.0 },       
        title: key);       
 
    gp3.WriteDeferredPlotCommands();       
}

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