Solving a linear system within PDELab

Here we show how to solve a linear system of equations originating from a PDE using PDELab.

First, we set up a GridOperator as in Assembling a linear system from a PDE

  typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,NumberType,NumberType,NumberType,CC,CC> GO;
  auto go = GO(gfs,cc,gfs,cc,lop,MBE(nonzeros));

Next, we set up our degree of freedom vector

  typedef Dune::PDELab::Backend::Vector<GFS,NumberType> X;
  X x(gfs,0.0);

and ensure it matches the Dirichlet boundary conditions at constrained degrees of freedom. In addition to specifying Dirichlet constrained degrees of freedom, it also serves as initial guess at unconstrained ones.

  typedef Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem> G;
  G g(grid.leafGridView(),problem);
  Dune::PDELab::interpolate(g,gfs,x);

Now we choose the preconditioner and solver we want to use

  typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<GO> LS;
  LS ls(100,3);

and plug it into a StationaryLinearProblemSolver. This takes care of assembling as well as solving the system.

  typedef Dune::PDELab::StationaryLinearProblemSolver<GO,LS,X> SLP;
  SLP slp(go,ls,x,1e-10);
  slp.apply(); // here all the work is done!

Finally, let's print the result to console via

  Dune::printvector(std::cout, native(x), "Solution", "");

There is a number of alternative solvers and preconditioners available we could use instead, for example this one:

  typedef Dune::PDELab::ISTLBackend_SEQ_CG_ILU0 LS;

Full example code: recipe-linear-system-solution-pdelab.cc