recipe-linear-system-solution-istl.cc

See explanation at Solving a linear system using ISTL

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
 
#include <dune/common/filledarray.hh>
#include <dune/common/parallel/mpihelper.hh>
 
#include <dune/grid/yaspgrid.hh>
 
#include <dune/pdelab.hh>
 
template<typename GV, typename RF>
class GenericEllipticProblem
{
  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
 
public:
  typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
 
  static constexpr bool permeabilityIsConstantPerCell()
  {
    return true;
  }
 
  typename Traits::PermTensorType
  A (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
  {
    typename Traits::PermTensorType I;
    for (std::size_t i=0; i<Traits::dimDomain; i++)
      for (std::size_t j=0; j<Traits::dimDomain; j++)
        I[i][j] = (i==j) ? 1 : 0;
    return I;
  }
 
  typename Traits::RangeType
  b (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
  {
    typename Traits::RangeType v(0.0);
    return v;
  }
 
  typename Traits::RangeFieldType
  c (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
  {
    return 0.0;
  }
 
  typename Traits::RangeFieldType
  f (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
  {
    return 0.0;
  }
 
  /* return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet for Dirichlet boundary conditions
   * return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Neumann for flux boundary conditions
   * return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Outflow for outflow boundary conditions
   */
  BCType
  bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
  {
    return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
  }
 
  typename Traits::RangeFieldType
  g (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
  {
    typename Traits::DomainType x = e.geometry().global(xlocal);
    if (x[0] < 0.5)
      return 1.0;
 
    return 0.0;
  }
 
  typename Traits::RangeFieldType
  j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
  {
    return 0.0;
  }
 
  typename Traits::RangeFieldType
  o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
  {
    return 0.0;
  }
};
 
 
int main(int argc, char **argv)
{
  // initialize MPI, finalize is done automatically on exit
  Dune::MPIHelper::instance(argc,argv);
 
  // define parameters
  typedef double NumberType;
 
  // need a grid in order to test grid functions
  constexpr unsigned int dim = 1;
  constexpr unsigned int degree = 1;
  constexpr std::size_t nonzeros = Dune::power(2*degree+1,dim);
 
  Dune::FieldVector<NumberType,dim> L(1.0);
  std::array<int,dim> N(Dune::filledArray<dim,int>(5));
 
  typedef Dune::YaspGrid<dim> Grid;
  Grid grid(L,N);
 
 
  // make grid
  typedef Dune::YaspGrid<dim> GM;
 
 
  // make problem parameters
  typedef GenericEllipticProblem<typename GM::LeafGridView,NumberType> Problem;
  Problem problem;
  typedef Dune::PDELab::ConvectionDiffusionBoundaryConditionAdapter<Problem> BCType;
  BCType bctype(grid.leafGridView(),problem);
 
 
  typedef typename GM::ctype DF;
  typedef Dune::PDELab::QkLocalFiniteElementMap<GM::LeafGridView,DF,NumberType,1> FEM;
  FEM fem(grid.leafGridView());
 
  typedef Dune::PDELab::GridFunctionSpace<GM::LeafGridView,FEM,
  Dune::PDELab::ConformingDirichletConstraints,
  Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed,1>> GFS;
  GFS gfs(grid.leafGridView(),fem);
 
  typedef typename GFS::template ConstraintsContainer<NumberType>::Type CC;
  CC cc;
  Dune::PDELab::constraints(bctype,gfs,cc);
 
  // local operator for finite elemenent problem
  typedef Dune::PDELab::ConvectionDiffusionFEM<Problem,FEM> LOP;
  LOP lop(problem);
 
  typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
  // [Grid operator]
  typedef Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,NumberType,NumberType,NumberType,CC,CC> GO;
  auto go = GO(gfs,cc,gfs,cc,lop,MBE(nonzeros));
 
  // [Make degree of freedom vector]
  typedef Dune::PDELab::Backend::Vector<GFS,NumberType> X;
  X x(gfs,0.0);
  // [Set it to match boundary conditions]
  typedef Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem> G;
  G g(grid.leafGridView(),problem);
  Dune::PDELab::interpolate(g,gfs,x);
 
  // [Assemble residual]
  X d(gfs,0.0);
  go.residual(x,d);
 
  // [Assemble matrix]
  typedef GO::Jacobian M;
  M A(go);
  go.jacobian(x,A);
 
  // [Extract ISTL types]
  typedef Dune::PDELab::Backend::Native<M> ISTLM;
  typedef Dune::PDELab::Backend::Native<X> ISTLX;
 
  // [Access to ISTL objects]
  using Dune::PDELab::Backend::native;
 
  // [Define operator]
  typedef Dune::MatrixAdapter<ISTLM,ISTLX,ISTLX> Operator;
  Operator matrixop(native(A));
 
  // [Define preconditioner and solver]
  Dune::SeqJac<ISTLM,ISTLX,ISTLX> preconditioner(native(A), 1, .5);
  Dune::CGSolver<ISTLX> solver(matrixop, preconditioner, 1e-3, 10, 2);
 
  // [Run solver]
  X v(gfs,0.0);
  Dune::InverseOperatorResult res;
  solver.apply(native(v), native(d), res);
 
  // [Subtract defect correction]
  x -= v;
 
  // [Solution output]
  Dune::printvector(std::cout, native(x), "Solution", "");
}