LinearSolver.hpp

#pragma once

#include <dune/istl/solver.hh>

#include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/linearalgebra/istl/ISTLSolverCreator.hpp>
#include <amdis/linearalgebra/istl/Preconditioners.hpp>
#include <amdis/linearalgebra/istl/Solvers.hpp>
#include <amdis/linearalgebra/istl/Traits.hpp>

namespace AMDiS
{
  template <class Mat, class VecX, class VecY = VecX>
  class LinearSolver
      : public LinearSolverInterface<Mat,VecX,VecY>
  {
    using Traits = SolverTraits<Mat,VecX,VecY>;
    using Map = CreatorMap<tag::solver<Traits>>;
    using Creator = ISTLSolverCreatorBase<Traits>;

  public:
    LinearSolver(std::string const& name, std::string const& prefix)
      : solverCreator_(dynamic_cast<Creator*>(Map::get(name, prefix)))
    {
      assert(solverCreator_ != nullptr);
      solverCreator_->init(prefix);
    }

    void init(Mat const& A) override
    {
      assert(solverCreator_ != nullptr);
      solver_ = solverCreator_->createSolver(A.matrix(), A.comm());
    }

    void finish() override
    {
      solver_.reset();
    }

    void apply(VecX& x, VecY const& b, Dune::InverseOperatorResult& stat) override
    {
      assert(!!solver_);
      auto rhs = b.vector();
      solver_->apply(x.vector(), rhs, stat);
    }

  private:
    ISTLSolverCreatorBase<Traits>* solverCreator_ = nullptr;
    std::shared_ptr<typename Traits::Solver> solver_;
  };

} // end namespace AMDiS