Kaskade 7 > Examples > Simple Example - Poisson's equation

Simple Example - Poisson's equation

Create a script that formulates the Poisson equation and specifies the details for Kaskade 7 to find the solution.

Problem Formulation

Consider the following problem where we search for $u(x,y)$ that satisfies the 2D Poisson’s equation: \[ -\Delta u = -\nabla \cdot ( \nabla u ) = 1 \quad\text{ in $\Omega$} \] \[ u = 0 \quad\text{ on $\partial\Omega$} \] Here $\Omega$ denotes the open unit square. Define \[ F(u) := \frac{1}{2}\nabla u^T\nabla u - u \quad\text{in $\Omega$} \] \[ G(u) := \gamma\frac{1}{2}u^2 \quad\text{on $\partial\Omega$.} \] Here, $\gamma$ is a large penalty parameter for the boundary condition. Then the above Poisson’s equation can be converted into the minimization problem

\[ \min_{u \in H^1(\Omega)} J(u) = \int_\Omega F(u)\, dx + \int_{\delta\Omega} G(u) \,dS. \]

Later we will see, how we should implement the functionals $F,G$ to be used in Kaskade 7.

Makefile Preparation

In later sessions we will create a script named laplace.cpp and a header file named laplace.hh. We will recommend here a basic makefile for the purpose of this tutorial. In the terminal, we begin by creating a directory and a makefile:

mkdir poisson
cd poisson
touch Makefile

The following is a simple example of the makefile for the purpose of this tutorial. Please note that you should replace #MAKEFILE_LOCAL_LOCATION# with the file path of Makefile.Local in the installation directory of Kaskade 7.

include #MAKEFILE_LOCAL_LOCATION#
include $(KASKADE7)/Makefile.Rules

ALLLIBS = $(KASKADELIB) $(DUNELIB) $(UGLIB) $(BOOSTLIB) $(UMFPACKLIB) $(BLASLIB) $(FTNLIB) $(NUMALIB)

laplace: laplace.o $(KASKADE7)/libs/libkaskade.a
	$(CXX) $(FLAGS) $< $(KASKADE7)/libs/umfpack_solve.o $(ALLLIBS) $(LINKFLAGS) -o $@

Grid and Variable Set Initialization

We will begin by defining a grid for our domain, and the variable $u$ (function) on such domain. We start with the main compilation unit laplace.cpp in the same directory with the following codes.

#include <iostream>

#include "dune/grid/config.h"
#include "dune/grid/uggrid.hh"

#include "fem/assemble.hh"
#include "fem/istlinterface.hh"
#include "fem/functional_aux.hh"
#include "fem/gridmanager.hh"
#include "fem/lagrangespace.hh"
#include "fem/variables.hh"

#include "io/vtk.hh"
#include "linalg/direct.hh"
#include "utilities/gridGeneration.hh"

using namespace Kaskade;
#include "laplace.hh"

int main()
{
  std::cout << "Start Laplacian tutorial program..." << std::endl;

  constexpr int dim = 2;
  constexpr double side_length = 1;
  constexpr int refinements = 5;
  constexpr int order = 2;
  constexpr double penalty = 1e6;

  // define the square grid
  using Grid = Dune::UGGrid<dim>;
  GridManager<Grid> temperatureGM( createUnitSquare<Grid>(side_length) );
  temperatureGM.globalRefine(refinements);

  // construction of finite element space for the scalar solution u.
  using LeafView = Grid::LeafGridView;
  using H1Space = FEFunctionSpace<ContinuousLagrangeMapper<double,LeafView>>;
  H1Space temperatureSpace(temperatureGM, temperatureGM.grid().leafGridView(), order);

  using Spaces = boost::fusion::vector<H1Space const*>;
  Spaces temperatureSpaces(&temperatureSpace);

  using VariableDescriptions = boost::fusion::vector<Variable<SpaceIndex<0>,Components<1>>>;
  using VariableSetDesc = VariableSetDescription<Spaces,VariableDescriptions>;
  VariableSetDesc temperatureSpacesVarSetDesc(temperatureSpaces,{ "u" });
  auto u = temperatureSpacesVarSetDesc.variableSet();

  // visualize the grid with the variable u initialized as zero everywhere
  writeVTK(u,"temperature_initial",IoOptions().setOrder(order).setPrecision(7));
  std::cout << "Grid and initial variable defined" << std::endl;
}

Below is a breakdown of the important parts of the script.

In the first part we state the dimension of the grid type we will use. We create the unit square grid with the function createUnitSquare<Grid> from gridGeneration.hh, which is passed and wrapped in a grid manager. Grid managers take care of the internal complications so that grid refinement and be performed easily and consistently.

using Grid = Dune::UGGrid<dim>;
GridManager<Grid> temperatureGM( createUnitSquare<Grid>(side_length) );
temperatureGM.globalRefine(refinements);

We specified the type of grid view when we define the type of continuous Lagrange mapper. A leaf grid view (as opposed to level grid view) is adopted here. A continuous Lagrange mapper is a manager of degree of freedom for the continuous finite element function space of piecewise polynomials. Here order specifies the degree of the Ansatz polynomials considered in the function space.

using LeafView = Grid::LeafGridView;
using H1Space = FEFunctionSpace<ContinuousLagrangeMapper<double,LeafView>>;
H1Space temperatureSpace(temperatureGM, temperatureGM.grid().leafGridView(), order);

Sometimes we might consider more than one function in the solver. In general we will need to pass a Boost vector of function space when we apply the later steps in the solver.

using Spaces = boost::fusion::vector<H1Space const*>;
Spaces temperatureSpaces(&temperatureSpace);

To specify the function variable we need for the problem. We will need to define the type VariableSetDescription, which link the type of function spaces (a vector) and a vector of variables. When we define the variable vector type, called VariableDescriptions here, we need to specify each of their space indices and components. Space index defines the order we arrange the variables and should range from $0$ to $n-1$, where $n$ is the number of varialbes in consideration. Component refers to the dimension of the output of the functional variable, which should be one unless vector functions are in consideration. We also use a string “u” to name the variable, which will be important when we need to identify the values in the results.

using VariableDescriptions = boost::fusion::vector<Variable<SpaceIndex<0>,Components<1>>>;
using VariableSetDesc = VariableSetDescription<Spaces,VariableDescriptions>;
VariableSetDesc temperatureSpacesVarSetDesc(temperatureSpaces,{ "u" });
auto u = temperatureSpacesVarSetDesc.variableSet();

Visualization of the Grid Using ParaView

The script outputs the variable u as a VTK file without performing any calculation. We can visualize the grid using software like ParaView. Before we execute the script, we should create an empty header file laplace.hh (details of this header file will be provided in later session) because our CPP file contains a line that imports the content of such header file.

Run the following command in the terminal to create a D file needed for the future make commands.

make laplace.d

Then we can compile the script by running

make laplace

To execute the script, we use the command

./laplace

After which we should find the file temperature_initial.vtu created from the program. In this tutorial, we will visualize it using ParaView, where a unit square with a value of zero should be observable once the file is loaded. Initial grid

When there are more than one variables, it is possible to choose another variable in the drop down menu: Choose variable

It is also possible to view the grid with the grid lines: View grid lines

Numerical Evaluation of the Solution of Poisson’s Equation

The following piece of codes should be appended within the main() function in laplace.cpp.

  // define a functional class which contains information associated with the Poisson's equation
  using Functional = HeatFunctional<double,VariableSetDesc>;
  Functional F(penalty);

  // assemble matrices and vectors which contain the information function space discretization and the differential equation
  using Assembler = VariationalFunctionalAssembler<LinearizationAt<Functional> >;
  Assembler assembler(temperatureSpaces);
  assembler.assemble(linearization(F,u));

  using Operator = AssembledGalerkinOperator<Assembler>;
  Operator A(assembler);

  using CoefficientVectors = VariableSetDesc::CoefficientVectorRepresentation<0,1>::type;
  CoefficientVectors solution = temperatureSpacesVarSetDesc.zeroCoefficientVector();
  CoefficientVectors rhs(assembler.rhs());

  // solve the system of linear equation and obtain the coefficients which represent the solution  
  directInverseOperator(A).applyscaleadd(-1.0,rhs,solution);
  component<0>(u) += component<0>(solution);

  // visualize the solution
  writeVTK(u,"temperature_output",IoOptions().setOrder(order).setPrecision(7));

  std::cout << "End Laplacian tutorial program" << std::endl;

We start by constructing a functional variable with a penalty parameter. The class HeatFunctional (a variational functional class) is defined in laplace.hh, details of which will be explained in the next section.

using Functional = HeatFunctional<double,VariableSetDesc>;
Functional F(penalty);

Using a variational functional assembler, one can compute the entries of the matrx and vectors of the linear system of equations representing the discretized differential equation in the discretized functional space of interest.

using Assembler = VariationalFunctionalAssembler<LinearizationAt<Functional> >;
Assembler assembler(temperatureSpaces);
assembler.assemble(linearization(F,u));

using Operator = AssembledGalerkinOperator<Assembler>;
Operator A(assembler);

using CoefficientVectors = VariableSetDesc::CoefficientVectorRepresentation<0,1>::type;
CoefficientVectors solution = temperatureSpacesVarSetDesc.zeroCoefficientVector();
CoefficientVectors rhs(assembler.rhs());

The solution can be computed by solving the system of equations $Au=rhs$. The variable $u$ already contains information of the discretized function space considered for the output of the values of the varialbe on the grid from the coefficients.

directInverseOperator(A).applyscaleadd(-1.0,rhs,solution);
component<0>(u) += component<0>(solution);

Defining the Variational Functional Class

A variational functional class named HeatFunctional is written in the header file laplace.hh. It should contain the following variables/structures:

Variables/Function	Purpose
`DomainCache`	Evaluation of F (in the domain) and its derivatives
`BoundaryCache`	Evaluation of G (on the boundary) and its derivatives
`D1`	Specification of block structures related to the first derivatives
`D2`	Specification of block structures related to the second derivatives
`integrationOrder()`	A required specification of the integration order required by the base class

In particular DomainCache and BoundaryCache are important class objects with the following members:

Variables/Functions	Purpose (for DomainCache)
`moveTo()`	Fix the value of class variables that are constant within an entity (such as a face)
`evaluateAt()`	Fix the value of class variables that are constant at a point in the domain
`d0()`	Evaluate the value of F at the specified point
`d1()`	Evaluate the first derivative F'
`d2()`	Evaluate the second derivative F''

Efficient computation is achieved by isolating the retrieval of certain values associated with a certain point. Therefore, the computation in d0, d1 and d2 will make use of the values computed in evaluteAt().

The following is the laplace.hh we use for this tutorial.

#include <algorithm>
#include "fem/diffops/boundaryConditions.hh"
#include "fem/functional_aux.hh"
#include "fem/fixdune.hh"
#include "fem/variables.hh"
#include "utilities/linalg/scalarproducts.hh"

// Simple Poisson equation -\Delta u = f with homogeneous Dirichlet boundary conditions
// on \Omega (which will be the unit square in our case).
// The corresponding variational functional on H^1_0(\Omega) is J = \int_\Omega |\nabla u|^2 - f*u dx.
// Let us call the integrand F(x,u,\nabla u) = |\nabla u|^2 - f*u.
template <class RType, class VarSet>
class HeatFunctional : public Kaskade::FunctionalBase<Kaskade::VariationalFunctional>
{
public:
  using Scalar = RType;
  using OriginVars = VarSet;
  using AnsatzVars = VarSet;
  using TestVars = VarSet;

  static constexpr int dim = AnsatzVars::Grid::dimension;
  static constexpr int uIdx = 0;
  static constexpr int uSpaceIdx = spaceIndex<AnsatzVars,uIdx>;

  // The domain cache defines F as well as its first and second directional derivative,
  // evaluated at the current iterate u.
  class DomainCache : public Kaskade::CacheBase<HeatFunctional,DomainCache>
  {
  public:
    DomainCache(HeatFunctional const&,
                typename AnsatzVars::VariableSet const& vars_,
                int flags=7):
      vars(vars_)
    {}


    template <class Position, class Evaluators>
    void evaluateAt(Position const& x, Evaluators const& evaluators)
    {
      u  = vars.template value<uIdx>(evaluators);
      du = vars.template derivative<uIdx>(evaluators);
      f = 1.0;
    }

    Scalar
    d0() const
    {
      return sp(du,du)/2 - f*u;
    }

    template<int row>
    D1Result<TestVars,row> d1(Kaskade::VariationalArg<Scalar,dim,TestVars::template Components<row>::m> const& arg) const
    {
      return sp(du,arg.derivative) - f*arg.value;
    }

    template<int row, int col>
    Scalar d2_impl(Kaskade::VariationalArg<Scalar,dim,TestVars::template Components<row>::m> const &arg1,
                   Kaskade::VariationalArg<Scalar,dim,AnsatzVars::template Components<row>::m> const &arg2) const
    {
      return sp(arg1.derivative,arg2.derivative);
    }

  private:
    typename AnsatzVars::VariableSet const& vars;
    Dune::FieldVector<Scalar,1> u, f;
    Dune::FieldVector<Scalar,dim> du;
    Kaskade::LinAlg::EuclideanScalarProduct sp;
  };

  // The boundary cache implements the boundary conditions. We do not implement Dirichlet boundary conditions
  // by eliminating degrees of freedom or equivalently modify the stiffness matrix and right hand side, but
  // by penalizing the deviation of u from the Dirichlet value ud (zero in our simple example), in the form of
  // \int_{\partial \Omega} g(x,u,\nabla u) ds.
  //
  // The boundary cache implements g and its first and second directional derivatives. Different penalizations
  // exist. The simplest is the quadratic Dirichlet penalty g = \gamma/2 * (u-ud)^2. For convergence on mesh
  // refinement, i.e. h->0, the penalty factor gamma must be scaled with the face diameter. This is provided
  // for convenience by the class DirichletPenaltyBoundary. A more sophisticated penalty formulation is
  // Nitsche's method, provided by the class DirichletNitscheBoundary. Here we use Nitsche's method on the
  // right edge of the unit square and the Dirichlet penalty method on the other sides.
  class BoundaryCache
  {
  public:
    BoundaryCache(HeatFunctional<RType,AnsatzVars> const& functional,
                  typename AnsatzVars::VariableSet const& vars_,
                  int flags=7):
      vars(vars_), penalty(functional.penalty), uDirichletBoundaryValue(0.0)
    {}

    template <class FaceIterator>
    void moveTo(FaceIterator const& fi)
    {
      useNitsche = fi->centerUnitOuterNormal()[0] > 0.5;
      if (useNitsche)
        nitsche.moveTo(fi);
      else
        dirichlet.moveTo(fi);
    }

    template <class Position, class Evaluators>
    void evaluateAt(Position const& xi, Evaluators const& evaluators)
    {
      auto u = vars.template value<uIdx>(evaluators);
      if (useNitsche)
      {
        auto ux = vars.template derivative<uIdx>(evaluators);
        nitsche.setBoundaryData(xi,penalty,u,uDirichletBoundaryValue,ux,unitMatrix<Scalar,dim>());
      }
      else
        dirichlet.setBoundaryData(penalty,u,uDirichletBoundaryValue);
    }

    Scalar d0() const
    {
      if (useNitsche)
        return nitsche.d0();
      else
        return dirichlet.d0();
    }

    template<int row>
    Dune::FieldVector<Scalar,1> d1(Kaskade::VariationalArg<Scalar,dim> const& argT) const
    {
      if (useNitsche)
        return nitsche.d1(argT);
      else
        return dirichlet.d1(argT);
    }

    template<int row, int col>
    Dune::FieldMatrix<Scalar,1,1> d2(Kaskade::VariationalArg<Scalar,dim> const &argT,
                                     Kaskade::VariationalArg<Scalar,dim> const &argA) const
    {
      if (useNitsche)
        return nitsche.d2(argT,argA);
      else
        return dirichlet.d2(argT,argA);
    }

  private:
    typename AnsatzVars::VariableSet const& vars;
    Scalar penalty;
    Dune::FieldVector<Scalar,1> uDirichletBoundaryValue;
    DirichletPenaltyBoundary<typename AnsatzVars::GridView,1> dirichlet;
    DirichletNitscheBoundary<typename AnsatzVars::GridView,1> nitsche;
    bool useNitsche;
  };

  template <int row>
  struct D1: public Kaskade::FunctionalBase<Kaskade::VariationalFunctional>::D1<row>
  {
    static bool const present   = true;
    static bool const constant  = false;

  };

  template <int row, int col>
  struct D2: public Kaskade::FunctionalBase<Kaskade::VariationalFunctional>::D2<row,col>
  {
    static bool const present = true;
    static bool const symmetric = true;
    static bool const lumped = false;
  };

  HeatFunctional(double gamma=1e6)
  : penalty(gamma)
  {}

  template <class Cell>
  int integrationOrder(Cell const& /* cell */, int shapeFunctionOrder, bool boundary) const
  {
    if (boundary)
      return 2*shapeFunctionOrder;
    else
    {
      int stiffnessMatrixIntegrationOrder = 2*(shapeFunctionOrder-1);
      int sourceTermIntegrationOrder = shapeFunctionOrder;        // as rhs f is constant, i.e. of order 0

      return std::max(stiffnessMatrixIntegrationOrder,sourceTermIntegrationOrder);
    }
  }

private:
  double penalty;
};

Page last modified: 2023-04-25 15:16:35 +0100 +0100