ParallelOverlappingILU0.hpp

/*
  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
  Copyright 2015 Statoil AS

  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
#define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
#include <opm/common/TimingMacros.hpp>
#include <opm/simulators/linalg/MILU.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <dune/istl/paamg/smoother.hh>

#include <cstddef>
#include <vector>
#include <type_traits>

namespace Opm
{

//template<class M, class X, class Y, class C>
//class ParallelOverlappingILU0;
template<class Matrix, class Domain, class Range, class ParallelInfo>
class ParallelOverlappingILU0;

template<class F>
class ParallelOverlappingILU0Args
    : public Dune::Amg::DefaultSmootherArgs<F>
{
 public:
    explicit ParallelOverlappingILU0Args(MILU_VARIANT milu = MILU_VARIANT::ILU )
        : milu_(milu), n_(0)
    {}
    void setMilu(MILU_VARIANT milu)
    {
        milu_ = milu;
    }
    MILU_VARIANT getMilu() const
    {
        return milu_;
    }
    void setN(int n)
    {
        n_ = n;
    }
    int getN() const
    {
        return n_;
    }
 private:
    MILU_VARIANT milu_;
    int n_;
};
} // end namespace Opm

namespace Dune
{

namespace Amg
{


template<class M, class X, class Y, class C>
struct SmootherTraits<Opm::ParallelOverlappingILU0<M,X,Y,C> >
{
    using Arguments = Opm::ParallelOverlappingILU0Args<typename M::field_type>;
};

template<class Matrix, class Domain, class Range, class ParallelInfo>
struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
{
    using T = Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>;
    using Arguments = DefaultParallelConstructionArgs<T,ParallelInfo>;

    using ParallelOverlappingILU0Pointer = std::shared_ptr<T>;

    static inline ParallelOverlappingILU0Pointer construct(Arguments& args)
    {
        return ParallelOverlappingILU0Pointer(
                new T(args.getMatrix(),
                      args.getComm(),
                      args.getArgs().getN(),
                      args.getArgs().relaxationFactor,
                      args.getArgs().getMilu()) );
    }
};

} // end namespace Amg

} // end namespace Dune

namespace Opm
{
template<class Matrix, class Domain, class Range, class ParallelInfoT>
class ParallelOverlappingILU0
    : public Dune::PreconditionerWithUpdate<Domain,Range>
{
    using ParallelInfo = ParallelInfoT;

public:
    using matrix_type = typename std::remove_const<Matrix>::type;
    using domain_type = Domain;
    using range_type = Range;
    using field_type = typename Domain::field_type;

    using block_type = typename matrix_type::block_type;
    using size_type = typename matrix_type::size_type;

    virtual bool hasPerfectUpdate() const override {
        return true;
    }

protected:
    struct CRS
    {
      CRS() : nRows_( 0 ) {}

      size_type rows() const { return nRows_; }

      size_type nonZeros() const
      {
        assert( rows_[ rows() ] != size_type(-1) );
        return rows_[ rows() ];
      }

      void resize( const size_type nRows )
      {
          if( nRows_ != nRows )
          {
            nRows_ = nRows ;
            rows_.resize( nRows_+1, size_type(-1) );
          }
      }

      void reserveAdditional( const size_type nonZeros )
      {
          const size_type needed = values_.size() + nonZeros ;
          if( values_.capacity() < needed )
          {
              const size_type estimate = needed * 1.1;
              values_.reserve( estimate );
              cols_.reserve( estimate );
          }
      }

      void push_back( const block_type& value, const size_type index )
      {
          values_.push_back( value );
          cols_.push_back( index );
      }

      void clear()
      {
          rows_.clear();
          values_.clear();
          cols_.clear();
          nRows_= 0;
      }

      std::vector<size_type> rows_;
      std::vector<block_type> values_;
      std::vector<size_type> cols_;
      size_type nRows_;
    };

public:
    Dune::SolverCategory::Category category() const override;

    ParallelOverlappingILU0 (const Matrix& A,
                             const int n, const field_type w,
                             MILU_VARIANT milu, bool redblack = false,
                             bool reorder_sphere = true);

    ParallelOverlappingILU0 (const Matrix& A,
                             const ParallelInfo& comm, const int n, const field_type w,
                             MILU_VARIANT milu, bool redblack = false,
                             bool reorder_sphere = true);

    ParallelOverlappingILU0 (const Matrix& A,
                             const field_type w, MILU_VARIANT milu,
                             bool redblack = false,
                             bool reorder_sphere = true);

    ParallelOverlappingILU0 (const Matrix& A,
                             const ParallelInfo& comm, const field_type w,
                             MILU_VARIANT milu, bool redblack = false,
                             bool reorder_sphere = true);

    ParallelOverlappingILU0 (const Matrix& A,
                             const ParallelInfo& comm,
                             const field_type w, MILU_VARIANT milu,
                             size_type interiorSize, bool redblack = false,
                             bool reorder_sphere = true);

    void pre (Domain&, Range&) override
    {}

    void apply (Domain& v, const Range& d) override;

    template <class V>
    void copyOwnerToAll( V& v ) const;

    void post (Range&) override
    {}

    void update() override;

protected:
    Range& reorderD(const Range& d);

    Domain& reorderV(Domain& v);

    void reorderBack(const Range& reorderedV, Range& v);

    std::unique_ptr<Matrix> ILU_;
    CRS lower_;
    CRS upper_;
    std::vector< block_type > inv_;
    std::vector< std::size_t > ordering_;
    Range reorderedD_;
    Domain reorderedV_;

    const ParallelInfo* comm_;
    const field_type w_;
    const bool relaxation_;
    size_type interiorSize_;
    const Matrix* A_;
    int iluIteration_;
    MILU_VARIANT milu_;
    bool redBlack_;
    bool reorderSphere_;
};

} // end namespace Opm
#endif