ISTLSolParams.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _ISTL_SOLPARAMS_H
#define _ISTL_SOLPARAMS_H
 
#include "ISTLSupport.h"
#include <dune/common/version.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/solvercategory.hh>
 
#include <iostream>
#include <set>
#include <string>
#include <vector>
 
 
class DomainDecomposition;
class LinSolParams;
class ProcessAdm;
 
 
namespace ISTL {
 
class BlockPreconditioner : public Preconditioner {
public:
  Dune::SolverCategory::Category category() const
  { return Dune::SolverCategory::sequential; }
 
  BlockPreconditioner(const ISTL::Mat& A, const DomainDecomposition& dd_,
                      const std::string& schurType);
 
  ~BlockPreconditioner() override
  {}
 
  void pre(ISTL::Vec& x, ISTL::Vec& b) override;
 
  void apply(ISTL::Vec& v, const ISTL::Vec& d) override;
 
  void post(ISTL::Vec& x) override;
 
  std::unique_ptr<ISTL::Preconditioner>& getBlockPre(size_t block)
  {
    return blockPre[block];
  }
 
  ISTL::Mat& getBlock(size_t block)
  {
    return blocks[block];
  }
 
  ISTL::Operator& getBlockOp(size_t block)
  {
    if (!blockOp[block])
      blockOp[block].reset(new ISTL::Operator(blocks[block]));
 
    return *blockOp[block];
  }
 
protected:
  static void extractBlock(ISTL::Mat& B, const ISTL::Mat& A,
                           const std::set<int>& eqs_row,
                           const std::set<int>& eqs_col,
                           const std::map<int,int>& eqs_col_g2l);
 
  static void subtractMatrices(ISTL::Mat& A,
                               const ISTL::Mat& B,
                               const ISTL::Mat& C);
 
  std::vector<std::unique_ptr<ISTL::Preconditioner>> blockPre; 
  std::vector<std::unique_ptr<ISTL::Operator>> blockOp; 
  std::vector<ISTL::Mat> blocks; 
  const DomainDecomposition& dd; 
};
 
#ifdef HAVE_MPI
template<class M, class X, class Y, class C>
class OverlappingSchwarzOperator : public Dune::AssembledLinearOperator<M,X,Y>
{
public:
  using matrix_type = M;
  using domain_type = X;
  using range_type = Y;
  using field_type = typename X::field_type;
  using communication_type = C;
 
  Dune::SolverCategory::Category category() const
  { return Dune::SolverCategory::sequential; }
 
  OverlappingSchwarzOperator (const matrix_type& A, const communication_type& com)
    : _A_(A), communication(com)
  {}
 
  void apply (const X& x, Y& y) const override
  {
//    Y y2(y.size());
    _A_.umv(x,y);
    communication.addOwnerOverlapToAll(y,y);
  }
 
  void applyscaleadd (field_type alpha, const X& x, Y& y) const override
  {
    Y y2(y.size());
    _A_.usmv(alpha, x, y2);
    y = 0;
    communication.addOwnerOverlapToAll(y2,y);
  }
 
  const matrix_type& getmat () const override
  {
    return _A_;
  }
 
private:
  const matrix_type& _A_;
  const communication_type& communication;
};
 
using ParMatrixAdapter =
  OverlappingSchwarzOperator<ISTL::Mat,
                             ISTL::Vec,
                             ISTL::Vec,
                             Dune::OwnerOverlapCopyCommunication<int,int>>;
#endif
 
}
 
 
class ISTLSolParams
{
public:
  ISTLSolParams(const LinSolParams& spar, const ProcessAdm& padm) :
    solParams(spar), adm(padm)
  {}
 
 
  std::tuple<std::unique_ptr<ISTL::InverseOperator>,
             std::unique_ptr<ISTL::Preconditioner>,
             std::unique_ptr<ISTL::Operator>>
    setupPC(ISTL::Mat& A);
 
  const LinSolParams& get() const { return solParams; }
 
protected:
  ISTL::Preconditioner* setupPCInternal(ISTL::Mat& A,
                                        ISTL::Operator& op,
                                        size_t block,
                                        std::unique_ptr<ISTL::InverseOperator>* solver);
 
  const LinSolParams& solParams; 
  const ProcessAdm& adm;      
};
 
#endif