StokesEquation

What's new?

• implementation of the incompressible, steady Stokes equation
• advanced: {implementation of the Stokes flow behind a grid as an application}

Prerequisites

• implementation of numerical fluxes, -> chapter NumFlux
• spatial operator, -> chapter DifferentialOperator
• implementation of the gradient operator, -> chapter PoissonAsASystem
• implementation of the SIP operator, -> chapter SIP

Problem statement

The Stokes-equation is given as

$$ -\frac{1}{Re} \Delta \vec{u} + \nabla p \ = \vec{g} $$$$ \text{div} (\vec{u}) = 0 $$

Where $Re \in \mathbb{R} $ denotes the Reynolds number. We consider two types of boundary conditions for the Stokes equation, Dirichlet (on $\Gamma_D \subset \Omega$) and Neumman (on $\Gamma_N \subset \Omega$). Those are defined as $$ \vec{u} =\vec{u}_D \ \text{ on } \Gamma_D\ \text{ (Dirichlet)}, \\ % \left( -\frac{1}{Re}\ \nabla \vec{u} + I_p \psi \right) \vec{n}\vert_{\delta \Omega} = 0 \ \text{ on } \Gamma_N\ \text{ (Neumann) } . $$

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 StokesEq.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\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();

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

using ilPSP.LinSolvers;
using BoSSS.Solution.Tecplot;
using ilPSP.Connectors.Matlab;

/// BoSSScmdSilent BoSSSexeSilent
using NUnit.Framework;

To inicate at which point of the boundary which condition is valid, i.e. wether a certain point is eiter Dirichlet or Neumann we define IsDirichletBndy which defines a mapping $$ \vec{x} \mapsto \{ \text{true}, \text{false} \}, $$ where true actually indicates a Dirichlet boundary. Since this function is defined as a global delegate, it can be altered later on. In the same manner, the function UDiri defines the Dirichlet-value for the velocity at the boundary.

static class BndyMap { 

    public static Func<double[],bool> IsDirichletBndy = null;

    public static Func<double[],double[]> UDiri = null;

}

Velocity divergence and pressure gradient

At first, we implement the velocity divergence, i.e. the continuity equation. We use the strong form, i.e. $$ b(\vec{u},v) = \oint_{\Gamma \backslash \GammaD} \bar{v} \quad \lbrack \lbrack {\vec{u}}\rbrack\rbrack \cdot \vec{n}\Gamma

 dA 
 -
 \int_{\Omega} \text{div}(\vec{u}) v ~ dV.

$$

public class Divergence : 
        BoSSS.Foundation.IEdgeForm, // edge integrals
        BoSSS.Foundation.IVolumeForm     // volume integrals
{
    /// The parameter list for the divergence is empty:
    public IList<string> ParameterOrdering { 
        get { return null; } 
    }
 
    /// We have a vector argument variable, 
    /// the velocity [ u, v ] = u
    /// (our trial function):
    public IList<String> ArgumentOrdering { 
        get { return new string[] { "u", "v" }; } 
    }
 
    public TermActivationFlags VolTerms {
        get {
            return TermActivationFlags.AllOn;
        }
    }
 
    public TermActivationFlags InnerEdgeTerms {
        get {
            return TermActivationFlags.AllOn; 
        }
    }
 
    public TermActivationFlags BoundaryEdgeTerms {
       get {
           return TermActivationFlags.AllOn;
        }
    }
 
    /// In the volume part, the integrand is div(u)*v :
    public double VolumeForm(ref CommonParamsVol cpv, 
        double[] U, double[,] GradU, 
        double V, double[] GradV) {
        double Acc = 0;
        for(int d = 0; d < cpv.D; d++) {
            Acc -= GradU[d,d]*V;
        }
        return Acc;
    }
 
    /// On interior cell boundaries, we use a velocity penalty,
    /// $\mean{v} \jump{u} \cdot n_\Gamma$:
    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 Acc = 0;
        for(int d = 0; d < inp.D; d++) {
            Acc += 0.5*(V_IN + V_OT)*(U_IN[d] - U_OT[d])*inp.Normal[d];
        }
        return Acc;
    }
 
    /// On the domain boundary, we have to distinguish between 
    /// Dirichlet- and Neumann-boundary conditions; the function
    /// \code{uDiri} defines which of the two actually applies:
    public double BoundaryEdgeForm(ref CommonParamsBnd inp, 
        double[] U_IN, double[,] GradU_IN, double V_IN, double[] GradV_OT) {
 
        double Acc = 0;
 
        if(!BndyMap.IsDirichletBndy(inp.X)) {
            // On the Neumann boundary, we do not know an outer value for the
            // velocity, so there is no penalization at all:
            Acc = 0;    
        } else {
            // On the Dirichlet boundary, the outer value for the velocity
            // is given by the function/delegate 'UDiri':
            double[] UD = BndyMap.UDiri(inp.X);
            for(int d = 0; d < inp.D; d++) {
                Acc += (U_IN[d] - UD[d])*inp.Normal[d]*V_IN;
            }
        }
        return Acc;
    }
}

The gradient-operator

We use the variational formulation of the gradient operator, as it is explained in the section concerning the Poisson System

class Gradient_d :
        BoSSS.Foundation.IEdgeForm, // edge integrals
        BoSSS.Foundation.IVolumeForm     // volume integrals
{
    public Gradient_d(int _d) {
        this.d = _d;
    }
 
    /// The component index of the gradient:
    int d;
 
    /// As usual, we do not use parameters:
    public IList<string> ParameterOrdering { 
        get { return null; } 
    }
 
    /// We have one argument, the pressure $\psi$:
    public IList<String> ArgumentOrdering { 
        get { return new string[] { "psi" }; } 
    }
 
    public TermActivationFlags VolTerms {
        get { return TermActivationFlags.AllOn; }
    }
 
    public TermActivationFlags InnerEdgeTerms {
        get { return (TermActivationFlags.AllOn); }
    }
 
    public TermActivationFlags BoundaryEdgeTerms {
       get { return TermActivationFlags.AllOn; }
    }
 
    /// The volume integrand, for a vector-valued test-function $\vec{v}$
    /// would be $-\operatorname{div}{\vec{v}} \psi$. Our test function $v$
    /// is scalar-valued, so e.g. for $\code{d} = 0$ we have
    /// $\vec{v} = (v,0)$. In this case, our volume integrand reduces as 
    /// $-\operatorname{div}{\vec{v}} \psi = -\partial_x v \psi$:
    public double VolumeForm(ref CommonParamsVol cpv, 
           double[] Psi, double[,] GradPsi, 
           double V, double[] GradV) {
 
        double Acc = 0;
        Acc -= Psi[0]*GradV[d];
        return Acc;
    }        
 
    /// On interior cell edges, we simply use a central-difference flux.
    /// Again, we consider a scalar test function, so we have
    /// $ \jump{\psi} \vec{v} \cdot \vec{n} = \jump{\psi} v n_d $,
    /// where $n_d$ is the $d$--th component of $\vec{n}$:
    public double InnerEdgeForm(ref CommonParams inp, 
        double[] Psi_IN, double[] Psi_OT, 
        double[,] GradPsi_IN, double[,] GradPsi_OT, 
        double V_IN, double V_OT, double[] GradV_IN, double[] GradV_OT) {
 
        double Acc = 0;
        Acc += 0.5*(Psi_IN[0] + Psi_OT[0])*inp.Normal[this.d]*(V_IN - V_OT);
        return Acc;  
     }
 
    public double BoundaryEdgeForm(ref CommonParamsBnd inp, 
        double[] Psi_IN, double[,] GradPsi_IN, double V_IN, double[] GradV_OT) {
 
        double Acc = 0;
        if(!BndyMap.IsDirichletBndy(inp.X)) {
            // On the Neumann boundary, we want the total stress to be zero,
            // so there is no contribution from the pressure:
            Acc = 0;
        } else {
            // On the Dirichlet boundary, we do not know an outer value for 
            // the pressure, so we have to take the inner value:
            Acc += Psi_IN[0]*inp.Normal[this.d]*V_IN;
        }        
        return Acc;              
    }
}

Tests on pressure gradient and velocity divergence

If our implementation is correct, we created a discretization of $$ \left[ \begin{array}{cc} 0 & \nabla \\ -\operatorname{div} & 0 \\ \end{array} \right] $$ so the matrix should have the form $$ \left[ \begin{array}{cc} 0 & B \\ B^T & 0 \\ \end{array} \right] =: M, $$ i.e. $M$ should be symmetric. We are testing this using a channel flow configuration using an equidistant grid.

• Domain: $\Omega := (0,10) \times (-1,1)$
• Neumann boundary (i.e. the outlet of the channel on the right): $\Gamma_N := \{ (x,y) | \ x = 10 \}$
• Dirichlet boundary (i.e. inlet and walls): $ \Gamma_D := \delta \Omega \setminus \Gamma_D$
• Velocity boundary value at inlet and walls: $ \vec{u}_D := (1 - y^2, 0)$

We create a grid, a DG basis for velocity and pressure and a variable mapping:

var xNodesChannel = GenericBlas.Linspace(0,10,31);// 30 cells in x-direction
var yNodesChannel = GenericBlas.Linspace(-1,1,7); // 6 cells in y-direction
var grdChannel    = Grid2D.Cartesian2DGrid(xNodesChannel,yNodesChannel);
var VelBChannel   = new Basis(grdChannel, 2);  // velocity basis
var PsiBChannel   = new Basis(grdChannel, 1);  // pressure basis 
var varMapChannel = new UnsetteledCoordinateMapping(
                       VelBChannel,VelBChannel,PsiBChannel); // variable mapping

We specify the boundary conditions as delegates:

Func<double[],bool> IsDirichletBndy_Channel 
        = (X => Math.Abs(X[0] - 10) > 1.0e-10); // its Dirichlet, if x != 10
Func<double[],double[]> UDiri_Channel 
        = (X => new double[2] { 1.0 - X[1]*X[1], 0});

Let's create the operator which contains only the pressure gradient and velocity divergence:

DifferentialOperator GradDiv = new DifferentialOperator(3,3, // 3 vars. in dom. & codom.
                           QuadOrderFunc.Linear(), // linear operator
                           "u", "v", "psi",  // names of domain variables
                           "mom_x", "mom_y", "conti"); // names of codom. vars
GradDiv.EquationComponents["mom_x"].Add(new Gradient_d(0)); 
GradDiv.EquationComponents["mom_y"].Add(new Gradient_d(1)); 
GradDiv.EquationComponents["conti"].Add(new Divergence());
GradDiv.Commit();

We create the matrix of the GradDiv-operator for the channel configuration. Before that, we have to set values for the global IsDirichletBndy and UDiri-variables.

BndyMap. IsDirichletBndy  = IsDirichletBndy_Channel;
BndyMap.UDiri             = UDiri_Channel;
var GradDivMatrix_Channel = GradDiv.ComputeMatrix(varMapChannel,
                                                  null,
                                                  varMapChannel);

Finally, we can test the symmetry of the matrix:

var ErrMtx = GradDivMatrix_Channel - GradDivMatrix_Channel.Transpose();
ErrMtx.InfNorm();

/// NUnit test (few random tests) BoSSScmdSilent
Assert.LessOrEqual(ErrMtx.InfNorm(), 1e-12);

Adding the viscous operator, forming the Stokes operator

We use the SIP-operator from chapter SIP to model the viscous terms:

public class Viscous : 
        IEdgeForm,   // edge integrals
        IVolumeForm, // volume integrals
        IEquationComponentCoefficient // update of coefficients required for penalty parameters 
{
    /// The velocity component:
    int d;
 
    public Viscous(int _d) {
        this.d = _d;    
    }
 
    /// We implement Reynolds number and the polynomial degree, 
    /// as well as the cell-wise length scales (required for 
    /// the computation of the penalty factor) as global, static variables.
    public static double Re;
   

 
    /// We do not use parameters:
    public IList<string> ParameterOrdering { 
        get { return new string[0]; } 
    }
 
    /// Depending on \code{d}, the argument variable
    /// should be either $u$ or $v$:
    public IList<String> ArgumentOrdering { 
        get { 
            switch(d) {
                case 0  : return new string[] { "u" }; 
                case 1  : return new string[] { "v" }; 
                default : throw new Exception();
            }
        } 
    }
 
    /// The \code{TermActivationFlags}, as usual:
    public TermActivationFlags VolTerms {
        get {
            return TermActivationFlags.GradUxGradV;
        }
    }
 
    public TermActivationFlags InnerEdgeTerms {
        get {
            return TermActivationFlags.AllOn;
        }
    }
 
    public TermActivationFlags BoundaryEdgeTerms {
       get {
           return TermActivationFlags.AllOn;
        }
    }
 
    /// The integrand for the volume integral:
    public double VolumeForm(ref CommonParamsVol cpv, 
           double[] U, double[,] GradU,
           double V, double[] GradV) {               
        double acc = 0;
        for(int d = 0; d < cpv.D; d++)
            acc += GradU[0, d] * GradV[d];
        return (1/Re)*acc;
    }
 
 
    /// The integrand for the integral on the inner edges:
    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];
 
            // the symmetry term -({{ \/v }} [[ u ]])*Normal
            Acc -= 0.5 * (GradV_IN[d] + GradV_OT[d])*(U_IN[0] - U_OT[0])
                       * inp.Normal[d];;
        }
 
        // the penalty term eta*[[u]]*[[v]]
        Acc += eta*(U_IN[0] - U_OT[0])*(V_IN - V_OT);
        return (1/Re)*Acc;
 
    }
 
    /// The integrand on boundary edges, i.e. on $\partial \Omega$:
    public double BoundaryEdgeForm(ref CommonParamsBnd inp, 
        double[] U_IN, double[,] GradU_IN, double V_IN, double[] GradV_IN) {
 
 
        double Acc = 0.0;
 
        if(!BndyMap.IsDirichletBndy(inp.X)) {
            // Neumann boundary conditions, i.e. zero-stress:
            Acc = 0;
        } else {
            // Dirichlet boundary conditions
            double uBnd = BndyMap.UDiri(inp.X)[d];
 
            for(int d = 0; d < inp.D; d++) { // loop over vector components 
                // consistency term:
                Acc -= (GradU_IN[0, d])*(V_IN) * inp.Normal[d];
                // symmetry term:
                Acc -= (GradV_IN[d])*(U_IN[0]- uBnd) * inp.Normal[d];
            }
 
            // penalty term
            double eta = PenaltyFactor(inp.jCellIn, -1);
            Acc += eta*(U_IN[0] - uBnd)*(V_IN);
        }
 
        return (1/Re)*Acc;
    }
            
    MultidimensionalArray cj;
    double penalty_base;
            
    double PenaltyFactor(int jCellIn, int jCellOut) {
        double PenaltySafety = 2;
        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;
    }
            
            
    /// Update of penalty length scales.
    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

        cj = ((GridData)(cs.GrdDat)).Cells.cj;
    }
}

(16,26): warning CS0649: Field 'Viscous.Re' is never assigned to, and will always have its default value 0

Finally, we are ready to implement the Stokes operator:

DifferentialOperator Stokes = new DifferentialOperator(3,3, // 3 vars. in dom. & codom.
                         QuadOrderFunc.Linear(), // linear operator
                         "u", "v", "psi",  // names of domain variables
                         "mom_x", "mom_y", "conti"); // names of codom. vars
Stokes.EquationComponents["mom_x"].Add(new Gradient_d(0)); 
Stokes.EquationComponents["mom_x"].Add(new Viscous(0)); 
Stokes.EquationComponents["mom_y"].Add(new Gradient_d(1)); 
Stokes.EquationComponents["mom_y"].Add(new Viscous(1));
Stokes.EquationComponents["conti"].Add(new Divergence());
Stokes.Commit();

Again, we create the matrix (now, for the Stokes operator) and check its symmetry;

we also have to set the Reynolds number and the polynomial degree before calling ComputeMatrix (since we are doing a rather dirty trick by using global variables).

BndyMap.IsDirichletBndy  = IsDirichletBndy_Channel;
BndyMap.UDiri            = UDiri_Channel;
Viscous.Re               = 20.0;
var StokesMatrix_Channel = Stokes.ComputeMatrix(varMapChannel,
                                                null,
                                                varMapChannel);

Testing the symmetry:

var ErrMtx1 = StokesMatrix_Channel - StokesMatrix_Channel.Transpose();
ErrMtx1.InfNorm()

1.8784973576657649E-13

/// NUnit test (few random tests) BoSSScmdSilent
Assert.LessOrEqual(ErrMtx1.InfNorm(), 1e-12);

We also verify that our Stokes-matrix has full rank, i.e. we show that matrix size and rank are equal:

StokesMatrix_Channel.NoOfRows

2700

StokesMatrix_Channel.rank()

2700

/// NUnit test (few random tests) BoSSScmdSilent
Assert.AreEqual(StokesMatrix_Channel.rank(), StokesMatrix_Channel.NoOfRows);

Solving the Stokes equation in the channel

We set the parameters and see whether we actually obtain the correct solution; the exact solution of our problem is $$ \vec{u}_{\text{ex}} = (1 - y^2, 0 )^T, \quad p_{\text{ex}} = \frac{200}{\text{Re}} - x \frac{2}{\text{Re}} $$ and since it is polynomial we should be able to obtain it \emph{exactly} in our velocity-pressure-space of degrees $(2,1)$.

BndyMap.IsDirichletBndy  = IsDirichletBndy_Channel;
BndyMap.UDiri            = UDiri_Channel;
Viscous.Re               = 20.0;
var StokesMatrix_Channel = Stokes.ComputeMatrix(varMapChannel,
                                                null,
                                                varMapChannel);
var StokesAffine_Channel = Stokes.ComputeAffine(varMapChannel,
                                                null,
                                                varMapChannel);

Now, we are ready to solve the stokes equation. \BoSSS\ provides us with a system $$ \texttt{StokesMatrix Channel} \cdot (u,v,\psi)

• \texttt{StokesAffine Channel} = 0, $$ so we have to multiply `StokesAffine_Channel` with $-1$ to get a right-hand-side.

double[] RHS = StokesAffine_Channel.CloneAs();
RHS.ScaleV(-1.0);

In order to store our solution, we have to create DG fields:

SinglePhaseField u               = new SinglePhaseField(VelBChannel,"u");
SinglePhaseField v               = new SinglePhaseField(VelBChannel,"v");
SinglePhaseField psi             = new SinglePhaseField(PsiBChannel,"psi");
CoordinateVector SolutionChannel = new CoordinateVector(u,v,psi);

Solve the linear system using a direct method:

StokesMatrix_Channel.Solve_Direct(SolutionChannel, RHS);

We export the solution to a Tecplot file, use Visit (or any other visualization software) to inspect the solution:

Tecplot("Box", SolutionChannel.Fields.ToArray());

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

done.

Advanced topics

Stokes flow behind a grid

We use the following setting:

• Domain: $ \Omega := (0,5) \times (-2,2) $
• Neumann boundary: $\Gamma_N := \{ (x,y) | \ x = 5 \} $
• Diriclet boundary: $ \Gamma_D := \partial \Omega \setminus \Gamma_D $
• Velocity at Dirichlet boundary: $ \vec{u}_D := (1 - (2 (y - \operatorname{floor}(y)) - 1)^2, 0) $

So, the boundary functions are:

Func<double[],bool> IsDirichletBndy_GridFlow 
        = (X => Math.Abs(X[0] - 5) > 1.0e-10); 
Func<double[],double[]> UDiri_GridFlow 
        = (X => new double[2] { 1.0 - (2*(X[1] - Math.Floor(X[1])) - 1).Pow2(),
                                0});

TODO:

the rest is for you! One hint: in $y$-direction, use some spacing so that you have cell boundaries at (least at) $y \in \{ -1, 0, 1 \}$.