dune-pdelab 2.7-git
Loading...
Searching...
No Matches
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

typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper.
Definition: backend/interface.hh:176

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

std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192

Now we can wrap our matrix in a LinearOperator object

typedef Dune::MatrixAdapter<ISTLM,ISTLX,ISTLX> Operator;
Operator matrixop(native(A));

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);

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

X v(gfs,0.0);
Dune::InverseOperatorResult res;
solver.apply(native(v), native(d), res);

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", "");

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