DUNE PDELab (git)

Solving a linear system using ISTL

Here we show how to solve a linear system of equations originating from a PDE directly using ISTL instead of PDELab's abstractions.

Note that generally, we recommend going the all-PDELab route as explained in Solving a linear system within PDELab.

First, we assemble our residual as in Assembling a linear system from a PDE

X d(gfs,0.0);
go.residual(x,d);

as well as the system matrix we want to solve:

typedef GO::Jacobian M;
M A(go);
go.jacobian(x,A);

So far, these are data types provided by PDELab. We can access the underlying ISTL data types via

typedef Dune::PDELab::Backend::Native<M> ISTLM;
typedef Dune::PDELab::Backend::Native<X> ISTLX;

and the ISTL vectors or matrices behind a PDELab object via:

using Dune::PDELab::Backend::native;

Now we can wrap our matrix in a LinearOperator object

Operator matrixop(native(A));
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:135

and pass that into our desired solver along with a preconditioner of our choice.

Dune::SeqJac<ISTLM,ISTLX,ISTLX> preconditioner(native(A), 1, .5);
Dune::CGSolver<ISTLX> solver(matrixop, preconditioner, 1e-3, 10, 2);
conjugate gradient method
Definition: solvers.hh:193
The sequential jacobian preconditioner.
Definition: preconditioners.hh:413

Applying the solver now returns the solution of the defect problem, i.e. A*v=d.

X v(gfs,0.0);
solver.apply(native(v), native(d), res);
Statistics about the application of an inverse operator.
Definition: solver.hh:50

Finally, we can subtract that from our initial guess to obtain the solution of A*x=b:

x -= v;

In the end, let's print the solution to console.

Dune::printvector(std::cout, native(x), "Solution", "");
void printvector(std::ostream &s, const V &v, std::string title, std::string rowtext, int columns=1, int width=10, int precision=2)
Print an ISTL vector.
Definition: io.hh:89

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

Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)