Step 3: Solve and VTK output

This step solves the linear system and writes the nodal solution for visualization.

Learning goals

After this page, you should be able to:

• solve with CGSolver and SeqSSOR preconditioning,
• read basic convergence information from InverseOperatorResult,
• write a .vtu file with VTKWriter.

Iterative solve

We solve Ax=b using a preconditioned conjugate gradient method. The solver works fine, since the the matrix A is symmetric positive definite. Also the preconditioner, a symmetric SSOR method, is compatible with this matrix structure.

We use dune-istl as linear algebra backend. There, the iterative solvers expect a linear operator that acts on the resiual vectors. A Matrix can be transformed into a linear operator using the MatrixAdapter wrapper, where the application of the operator is simply the matrix-vector product.

using Operator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
Operator op(A);
Dune::SeqSSOR<Matrix, Vector, Vector> preconditioner(A, 1, 1.0);
Dune::CGSolver<Vector> solver(op, preconditioner, 1e-12, 5000, 0);

Dune::InverseOperatorResult res;
solver.apply(x, b, res);

Optional direct solver with UMFPack

The module dune-istl provides a large collection of iterative solvers, but also interfaces to external solver libraries, which provide sparse direct solvers, like the SuiteSparse library. This can be used, after activating the corresponding package for the code. In CMake, call add_dune_suitesparse_flags(example09) to set necessary compile flags and link libraries. Then you can switch the solver implementation to Dune::UMFPack<Matrix>.

VTK export

We copy the nodal vector to a plain std::vector<double> and write:

Dune::VTKWriter<GridView> vtk(gv, Dune::VTK::conforming);
vtk.addVertexData(nodalSolution, "u_h");
vtk.write("poisson_p1");

This produces poisson_p1.vtu.

Runnable example: examples/09_poisson_core/poisson_core.cc

#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

#include <dune/common/fvector.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/localfunctions/lagrange/lagrangecube.hh>

namespace {

constexpr int dim = 2;
using Grid = Dune::YaspGrid<dim>;
using GridView = Grid::LeafGridView;
using Matrix = Dune::BCRSMatrix<double>;
using Vector = Dune::BlockVector<double>;
using LocalFiniteElement = Dune::LagrangeCubeLocalFiniteElement<double, double, dim, 1>;

Grid makeGrid()
{
  return Grid({1.0, 1.0}, {32, 32});
}

std::array<std::size_t, 4> elementDofs(const GridView& gv, const GridView::Codim<0>::Entity& e)
{
  const auto& indexSet = gv.indexSet();
  std::array<std::size_t, 4> v{};
  for (int i = 0; i < 4; ++i)
    v[i] = indexSet.subIndex(e, i, dim);
  return v;
}

Matrix setupMatrixPattern(const GridView& gv, std::size_t numDofs)
{
  Dune::MatrixIndexSet pattern(numDofs, numDofs);
  for (const auto& e : Dune::elements(gv)) {
    const auto v = elementDofs(gv, e);
    for (int i = 0; i < 4; ++i)
      for (int j = 0; j < 4; ++j)
        pattern.add(v[i], v[j]);
  }

  Matrix A;
  pattern.exportIdx(A);
  A = 0.0;
  return A;
}

void assembleSystem(const GridView& gv, Matrix& A, Vector& b)
{
  LocalFiniteElement lfe;
  const auto& basis = lfe.localBasis();
  std::vector<LocalFiniteElement::Traits::LocalBasisType::Traits::RangeType> phi;
  std::vector<LocalFiniteElement::Traits::LocalBasisType::Traits::JacobianType> jacobian;

  for (const auto& e : Dune::elements(gv)) {
    const auto v = elementDofs(gv, e);
    const auto geo = e.geometry();
    const auto quadrature = Dune::QuadratureRules<double, dim>::rule(geo.type(), 2);

    for (const auto& qp : quadrature) {
      const auto local = qp.position();
      basis.evaluateFunction(local, phi);
      basis.evaluateJacobian(local, jacobian);

      const auto jacInvT = geo.jacobianInverseTransposed(local);
      const double dx = qp.weight() * geo.integrationElement(local);

      std::array<Dune::FieldVector<double, dim>, 4> grad{};
      for (int i = 0; i < 4; ++i)
        jacInvT.mv(jacobian[i][0], grad[i]);

      for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j)
          A[v[i]][v[j]] += dx * (grad[i] * grad[j]);
        b[v[i]] += dx * phi[i][0]; // f(x)=1
      }
    }
  }
}

std::vector<unsigned char> markDirichletDofs(const GridView& gv, std::size_t numDofs)
{
  std::vector<unsigned char> isDirichlet(numDofs, false);
  const auto& indexSet = gv.indexSet();
  const double eps = 1e-12;

  for (const auto& vertex : Dune::vertices(gv)) {
    const auto x = vertex.geometry().center();
    const auto i = indexSet.index(vertex);
    if (x[0] < eps || x[0] > 1.0 - eps || x[1] < eps || x[1] > 1.0 - eps)
      isDirichlet[i] = true;
  }

  return isDirichlet;
}

void applyHomogeneousDirichlet(const std::vector<unsigned char>& isDirichlet, Matrix& A, Vector& b)
{
  const std::size_t n = b.size();
  for (std::size_t i = 0; i < n; ++i) {
    for (auto entry = A[i].begin(); entry != A[i].end(); ++entry) {
      const auto j = entry.index();
      if (isDirichlet[i] || isDirichlet[j])
        *entry = (isDirichlet[i] && i == j) ? 1.0 : 0.0;
    }
    if (isDirichlet[i])
      b[i] = 0.0;
  }
}

Dune::InverseOperatorResult solveSystem(const Matrix& A, Vector& b, Vector& x)
{
  using Operator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
  Operator op(A);
  Dune::SeqSSOR<Matrix, Vector, Vector> preconditioner(A, 1, 1.0);
  Dune::CGSolver<Vector> solver(op, preconditioner, 1e-12, 5000, 0);

  Dune::InverseOperatorResult res;
  solver.apply(x, b, res);
  return res;
}

std::size_t centerDof(const GridView& gv)
{
  const auto& indexSet = gv.indexSet();
  std::size_t center = 0;
  double best = std::numeric_limits<double>::max();
  for (const auto& vertex : Dune::vertices(gv)) {
    const auto i = indexSet.index(vertex);
    const auto p = vertex.geometry().center();
    const double dist = std::abs(p[0] - 0.5) + std::abs(p[1] - 0.5);
    if (dist < best) {
      best = dist;
      center = i;
    }
  }
  return center;
}

void writeVtu(const GridView& gv, const Vector& x)
{
  const auto& indexSet = gv.indexSet();
  std::vector<double> nodalSolution(x.size(), 0.0);
  for (const auto& vertex : Dune::vertices(gv)) {
    const auto i = indexSet.index(vertex);
    nodalSolution[i] = x[i];
  }

  Dune::VTKWriter<GridView> vtk(gv, Dune::VTK::conforming);
  vtk.addVertexData(nodalSolution, "u_h");
  vtk.write("poisson_p1");
}

} // namespace

int main()
{
  Grid grid = makeGrid();
  auto gv = grid.leafGridView();
  const auto numDofs = gv.indexSet().size(dim);

  Matrix A = setupMatrixPattern(gv, numDofs);
  Vector b(numDofs);
  b = 0.0;

  assembleSystem(gv, A, b);
  const auto isDirichlet = markDirichletDofs(gv, numDofs);
  applyHomogeneousDirichlet(isDirichlet, A, b);

  Vector x(numDofs);
  x = 0.0;
  const auto res = solveSystem(A, b, x);

  writeVtu(gv, x);

  std::size_t boundaryCount = 0;
  for (const auto flag : isDirichlet)
    boundaryCount += static_cast<std::size_t>(flag);

  std::cout << "elements: " << gv.indexSet().size(0) << "\n";
  std::cout << "vertices/dofs: " << numDofs << "\n";
  std::cout << "boundary dofs: " << boundaryCount << "\n";
  std::cout << std::fixed << std::setprecision(10);
  std::cout << "u_h(0.5,0.5) approx: " << x[centerDof(gv)] << "\n";
  std::cout << "solver converged: " << std::boolalpha << res.converged << "\n";
  std::cout << "solver iterations: " << res.iterations << "\n";
  std::cout << "solver reduction: " << res.reduction << "\n";
  std::cout << "wrote: poisson_p1.vtu\n";

  return 0;
}

Build and run:

source _env/activate.sh
cd examples/09_poisson_core
dunecontrol --current all
./build/debug/example09

Output:

elements: 1024
vertices/dofs: 1089
boundary dofs: 128
u_h(0.5,0.5) approx: 0.0737281169
solver converged: true
solver iterations: 41
solver reduction: 0.0000000000
wrote: poisson_p1.vtu

Summary and next steps

You now have a complete low-level Poisson workflow with DUNE core modules and dune-localfunctions basis evaluation. Return to Core modules or extend the example with nonzero boundary data and mesh refinement studies.