LinearSolver.hpp

#pragma once
 
#include <petscpc.h>
#include <petscksp.h>
 
#include <amdis/Initfile.hpp>
#include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/linearalgebra/petsc/Interface.hpp>
#include <amdis/linearalgebra/petsc/SolverConfig.hpp>
 
 
namespace AMDiS
{
 
  template <class Matrix, class VectorX, class VectorY = VectorX>
  class LinearSolver
      : public LinearSolverInterface<Matrix,VectorX,VectorY>
  {
    using Vector = VectorX;
 
  public:
 
    LinearSolver(std::string const& /*name*/, std::string const& prefix)
      : prefix_(prefix)
    {
      Parameters::get(prefix + "->info", info_);
    }
 
    ~LinearSolver()
    {
      finish();
    }
 
    void init(Matrix const& A) override
    {
      finish();
      KSPCreate(A.comm(), &ksp_);
 
      PETSc::KSPSetOperators(ksp_, A.matrix(), A.matrix());
      initKSP(ksp_, prefix_);
      initialized_ = true;
    }
 
    void finish() override
    {
      if (initialized_)
        KSPDestroy(&ksp_);
      initialized_ = false;
    }
 
    void apply(VectorX& x, VectorY const& b, Dune::InverseOperatorResult& stat) override
    {
      KSPSolve(ksp_, b.vector(), x.vector());
 
      if (info_ > 0)
        PETSc::KSPConvergedReasonView(ksp_, PETSC_VIEWER_STDOUT_WORLD);
 
      KSPConvergedReason reason;
      KSPGetConvergedReason(ksp_, &reason);
 
      stat.converged = (reason > 0);
    }
 
    KSP ksp() { return ksp_; }
 
  protected:
    // initialize the KSP solver from the initfile
    virtual void initKSP(KSP ksp, std::string prefix) const
    {
      PETSc::configKSP(ksp, prefix);
 
      // show details of ksp object
      if (info_ >= 10)
        KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD);
    }
 
  protected:
    std::string prefix_;
    int info_ = 0;
 
    KSP ksp_;
    bool initialized_ = false;
  };
 
} // end namespace AMDiS