Merge pull request #1788 from blattms/amg_cpr_dune_2.5_rebased

Amg cpr reusing hierarchy and supporting dune 2.5

CMakeLists.txt

@@ -157,6 +157,16 @@ opm_add_test(flow
   flow/flow_ebos_oilwater_polymer.cpp
   flow/flow_ebos_oilwater_polymer_injectivity.cpp)
 
+opm_add_test(flow_blackoil_dunecpr
+  ONLY_COMPILE
+  DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF}
+  SOURCES flow/flow_blackoil_dunecpr.cpp
+  EXE_NAME flow_blackoil_dunecpr
+  DEPENDS "opmsimulators"
+  LIBRARIES "opmsimulators")
+
+
+
 if (BUILD_FLOW)
   install(TARGETS flow DESTINATION bin)
   opm_add_bash_completion(flow)

CMakeLists_files.cmake

@@ -118,6 +118,9 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/AquiferCarterTracy.hpp
   opm/autodiff/AquiferFetkovich.hpp
   opm/autodiff/BlackoilAmg.hpp
+  opm/autodiff/BlackoilAmgCpr.hpp
+  opm/autodiff/amgcpr.hh
+  opm/autodiff/twolevelmethodcpr.hh
   opm/autodiff/BlackoilDetails.hpp
   opm/autodiff/BlackoilModelParametersEbos.hpp
   opm/autodiff/BlackoilAquiferModel.hpp
@@ -129,6 +132,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/FlowMainEbos.hpp
   opm/autodiff/GraphColoring.hpp
   opm/autodiff/ISTLSolverEbos.hpp
+  opm/autodiff/ISTLSolverEbosCpr.hpp
   opm/autodiff/IterationReport.hpp
   opm/autodiff/MatrixBlock.hpp
   opm/autodiff/moduleVersion.hpp

flow/flow_blackoil_dunecpr.cpp

@@ -0,0 +1,96 @@
+/*
+  Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015, 2017 IRIS 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/>.
+*/
+#include "config.h"
+#include "flow/flow_tag.hpp"
+//#include <opm/linearsolvers/amgclsolverbackend.hh>
+#include  <opm/autodiff/ISTLSolverEbosCpr.hpp>
+//#include <ewoms/linear/superlubackend.hh>
+
+BEGIN_PROPERTIES
+NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem));
+NEW_PROP_TAG(FluidState);
+//SET_TYPE_PROP(EclBaseProblem, Problem, Ewoms::EclProblem<TypeTag>);
+SET_PROP(EclFlowProblemSimple, FluidState)
+    {
+    private:
+      typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+      typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+      enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
+      enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+      enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+      enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+      typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+      typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+      static const bool compositionSwitchEnabled = Indices::gasEnabled;
+ 
+    public:
+//typedef Opm::BlackOilFluidSystemSimple<Scalar> type;
+       typedef Opm::BlackOilFluidState<Evaluation, FluidSystem, enableTemperature, enableEnergy, compositionSwitchEnabled,  Indices::numPhases > type;
+};
+
+SET_BOOL_PROP(EclFlowProblemSimple,MatrixAddWellContributions,true);
+SET_INT_PROP(EclFlowProblemSimple,LinearSolverVerbosity,1);
+SET_SCALAR_PROP(EclFlowProblemSimple, LinearSolverReduction, 1e-4);
+SET_INT_PROP(EclFlowProblemSimple, LinearSolverMaxIter, 20);
+SET_BOOL_PROP(EclFlowProblemSimple, UseAmg, true);//probably not used
+SET_BOOL_PROP(EclFlowProblemSimple, UseCpr, true);
+SET_INT_PROP(EclFlowProblemSimple, CprMaxEllIter, 1);
+SET_INT_PROP(EclFlowProblemSimple, CprEllSolvetype, 3);
+SET_INT_PROP(EclFlowProblemSimple, CprReuseSetup, 3);
+SET_INT_PROP(EclFlowProblemSimple, CprSolverVerbose, 3);
+SET_STRING_PROP(EclFlowProblemSimple, SystemStrategy, "quasiimpes");
+END_PROPERTIES
+
+namespace Ewoms {
+  namespace Properties {
+
+    SET_PROP(EclFlowProblemSimple, FluidSystem)
+    {
+    private:
+      //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+      typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+      typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+
+    public:
+        typedef Opm::BlackOilFluidSystem<Scalar> type;
+    };
+    //NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem));
+    SET_TYPE_PROP(EclFlowProblemSimple, IntensiveQuantities, Ewoms::BlackOilIntensiveQuantities<TypeTag>);
+    //SET_TYPE_PROP(EclFlowProblemSimple, LinearSolverBackend, Opm::ISTLSolverEbos<TypeTag>);
+    //SET_TAG_PROP(EclFlowProblemSimple, LinearSolverSplice, ParallelBiCGStabLinearSolver);
+    //SET_TYPE_PROP(EclFlowProblemSimple, LinearSolverBackend, Ewoms::Linear::ParallelBiCGStabSolverBackend<TypeTag>);//not work
+    //SET_TYPE_PROP(EclFlowProblemSimple, LinearSolverBackend, Ewoms::Linear::SuperLUBackend<TypeTag>)//not work
+    //SET_TAG_PROP(EclFlowProblem, FluidState, Opm::BlackOilFluidState);
+    SET_TYPE_PROP(EclFlowProblemSimple, LinearSolverBackend, Opm::ISTLSolverEbosCpr<TypeTag>);
+    SET_BOOL_PROP(EclFlowProblemSimple, EnableStorageCache, true);
+    SET_BOOL_PROP(EclFlowProblemSimple, EnableIntensiveQuantityCache, true);
+    
+    //SET_INT_PROP(EclFlowProblemSimple, NumWellAdjoint, 1);
+    //SET_BOOL_PROP(EclFlowProblem, EnableStorageCache, true);
+    //SET_BOOL_PROP(EclFlowProblem, EnableIntensiveQuantityCache, true);
+  }
+}
+
+int main(int argc, char** argv)
+{
+  typedef TTAG(EclFlowProblemSimple) TypeTag;
+  return mainFlow<TypeTag>(argc, argv);
+}

flow/flow_tag.hpp

@@ -0,0 +1,212 @@
+/*
+  Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015, 2017 IRIS 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 FLOW_TAG_HPP
+#define FLOW_TAG_HPP
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/autodiff/FlowMainEbos.hpp>
+#include <ewoms/common/propertysystem.hh>
+#include <ewoms/common/parametersystem.hh>
+#include <opm/autodiff/MissingFeatures.hpp>
+#include <opm/common/utility/parameters/ParameterGroup.hpp>
+#include <opm/material/common/ResetLocale.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
+
+//#include <opm/material/fluidsystems/BlackOilFluidSystemSimple.hpp>
+//#include <opm/material/fluidsystems/BlackOilFluidSystemSimple.hpp> 
+#include <ewoms/models/blackoil/blackoilintensivequantities.hh>
+#include <opm/material/fluidstates/BlackOilFluidState.hpp>
+//#include <opm/material/fluidstates/BlackOilFluidStateSimple.hpp>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/misc/mpimanager.hh>
+#else
+#include <dune/common/parallel/mpihelper.hh>
+#endif
+
+
+BEGIN_PROPERTIES
+
+// this is a dummy type tag that is used to setup the parameters before the actual
+// simulator.
+NEW_TYPE_TAG(FlowEarlyBird, INHERITS_FROM(EclFlowProblem));
+
+END_PROPERTIES
+
+
+
+
+namespace Opm {
+  template <class TypeTag>
+  void flowEbosSetDeck(Deck &deck, EclipseState& eclState)
+  {
+    typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard;
+    Vanguard::setExternalDeck(&deck, &eclState);
+  }
+  
+// ----------------- Main program -----------------
+  template <class TypeTag>
+  int flowEbosMain(int argc, char** argv)
+  {
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    Opm::resetLocale();
+    
+#if HAVE_DUNE_FEM
+    Dune::Fem::MPIManager::initialize(argc, argv);
+#else
+    Dune::MPIHelper::instance(argc, argv);
+#endif
+    Opm::FlowMainEbos<TypeTag> mainfunc;
+    return mainfunc.execute(argc, argv);
+  }
+
+}
+
+
+
+
+
+namespace detail
+{
+    boost::filesystem::path simulationCaseName( const std::string& casename ) {
+        namespace fs = boost::filesystem;
+
+        const auto exists = []( const fs::path& f ) -> bool {
+            if( !fs::exists( f ) ) return false;
+
+            if( fs::is_regular_file( f ) ) return true;
+
+            return fs::is_symlink( f )
+            && fs::is_regular_file( fs::read_symlink( f ) );
+        };
+
+        auto simcase = fs::path( casename );
+
+        if( exists( simcase ) ) {
+            return simcase;
+        }
+
+        for( const auto& ext : { std::string("data"), std::string("DATA") } ) {
+            if( exists( simcase.replace_extension( ext ) ) ) {
+                return simcase;
+            }
+        }
+
+        throw std::invalid_argument( "Cannot find input case " + casename );
+    }
+}
+
+
+// ----------------- Main program -----------------
+template<class TypeTag>
+int mainFlow(int argc, char** argv)
+{
+    // MPI setup.
+#if HAVE_DUNE_FEM
+    Dune::Fem::MPIManager::initialize(argc, argv);
+    int mpiRank = Dune::Fem::MPIManager::rank();
+#else
+    // the design of the plain dune MPIHelper class is quite flawed: there is no way to
+    // get the instance without having the argc and argv parameters available and it is
+    // not possible to determine the MPI rank and size without an instance. (IOW: the
+    // rank() and size() methods are supposed to be static.)
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+    int mpiRank = mpiHelper.rank();
+#endif
+
+    // we always want to use the default locale, and thus spare us the trouble
+    // with incorrect locale settings.
+    Opm::resetLocale();
+
+    // this is a work-around for a catch 22: we do not know what code path to use without
+    // parsing the deck, but we don't know the deck without having access to the
+    // parameters and this requires to know the type tag to be used. To solve this, we
+    // use a type tag just for parsing the parameters before we instantiate the actual
+    // simulator object. (Which parses the parameters again, but since this is done in an
+    // identical manner it does not matter.)
+    typedef TTAG(FlowEarlyBird) PreTypeTag;
+    typedef GET_PROP_TYPE(PreTypeTag, Problem) PreProblem;
+
+    PreProblem::setBriefDescription("Simple Flow, an advanced reservoir simulator for ECL-decks provided by the Open Porous Media project.");
+    
+    int status = Opm::FlowMainEbos<PreTypeTag>::setupParameters_(argc, argv);
+    if (status != 0)
+        // if setupParameters_ returns a value smaller than 0, there was no error, but
+        // the program should abort. This is the case e.g. for the --help and the
+        // --print-properties parameters.
+        return (status >= 0)?status:0;
+
+    bool outputCout = false;
+    if (mpiRank == 0)
+        outputCout = EWOMS_GET_PARAM(PreTypeTag, bool, EnableTerminalOutput);
+
+    std::string deckFilename = EWOMS_GET_PARAM(PreTypeTag, std::string, EclDeckFileName);
+    typedef typename GET_PROP_TYPE(PreTypeTag, Vanguard) PreVanguard;
+    try {
+        deckFilename = PreVanguard::canonicalDeckPath(deckFilename).string();
+    }
+    catch (const std::exception& e) {
+        Ewoms::Parameters::printUsage<PreTypeTag>(PreProblem::helpPreamble(argc, const_cast<const char**>(argv)),
+                                                  e.what());
+        return 1;
+    }
+
+    // Create Deck and EclipseState.
+    try {
+      Opm::Parser parser;
+      typedef std::pair<std::string, Opm::InputError::Action> ParseModePair;
+      typedef std::vector<ParseModePair> ParseModePairs;
+      ParseModePairs tmp;
+      tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_RANDOM_SLASH, Opm::InputError::IGNORE));
+      tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_MISSING_DIMS_KEYWORD, Opm::InputError::WARN));
+      tmp.push_back(ParseModePair(Opm::ParseContext::SUMMARY_UNKNOWN_WELL, Opm::InputError::WARN));
+      tmp.push_back(ParseModePair(Opm::ParseContext::SUMMARY_UNKNOWN_GROUP, Opm::InputError::WARN));
+      Opm::ParseContext parseContext(tmp);
+      Opm::ErrorGuard errorGuard;
+
+      std::shared_ptr<Opm::Deck> deck = std::make_shared< Opm::Deck >( parser.parseFile(deckFilename , parseContext, errorGuard) );
+      if ( outputCout ) {
+          Opm::checkDeck(*deck, parser, parseContext, errorGuard);
+	Opm::MissingFeatures::checkKeywords(*deck);
+      }
+      //Opm::Runspec runspec( *deck );
+      //const auto& phases = runspec.phases();
+
+      std::shared_ptr<Opm::EclipseState> eclipseState = std::make_shared< Opm::EclipseState > ( *deck, parseContext, errorGuard );
+      Opm::flowEbosSetDeck<TypeTag>(*deck, *eclipseState);
+      return Opm::flowEbosMain<TypeTag>(argc, argv);
+    }
+  catch (const std::invalid_argument& e)
+    {
+      if (outputCout) {
+	std::cerr << "Failed to create valid EclipseState object." << std::endl;
+	std::cerr << "Exception caught: " << e.what() << std::endl;
+      }
+      throw;
+    }
+    
+    return EXIT_SUCCESS;
+}
+#endif

opm/autodiff/BlackoilAmg.hpp

@@ -23,6 +23,8 @@
 #include <opm/autodiff/ParallelOverlappingILU0.hpp>
 #include <opm/autodiff/FlowLinearSolverParameters.hpp>
 #include <opm/autodiff/CPRPreconditioner.hpp>
+#include <opm/autodiff/amgcpr.hh>
+#include <opm/autodiff/twolevelmethodcpr.hh>
 #include <dune/istl/paamg/twolevelmethod.hh>
 #include <dune/istl/paamg/aggregates.hh>
 #include <dune/istl/bvector.hh>
@@ -77,18 +79,6 @@ namespace Opm
 namespace Detail
 {
 
-/**
- * \brief Creates a MatrixAdapter as an operator
- *
- * The first argument is used to specify the return type using function overloading.
- * \param matrix The matrix to wrap.
- */
-template<class M, class X, class Y, class T>
-Dune::MatrixAdapter<M,X,Y> createOperator(const Dune::MatrixAdapter<M,X,Y>&, const M& matrix, const T&)
-{
-    return Dune::MatrixAdapter<M,X,Y>(matrix);
-}
-
 /**
  * \brief Creates a MatrixAdapter as an operator, storing it in a unique_ptr.
  *
@@ -96,22 +86,24 @@ Dune::MatrixAdapter<M,X,Y> createOperator(const Dune::MatrixAdapter<M,X,Y>&, con
  * \param matrix The matrix to wrap.
  */
 template<class M, class X, class Y, class T>
-std::unique_ptr< Dune::MatrixAdapter<M,X,Y> > createOperatorPtr(const Dune::MatrixAdapter<M,X,Y>&, const M& matrix, const T&)
+std::unique_ptr< Dune::MatrixAdapter<M,X,Y> >
+createOperator(const Dune::MatrixAdapter<M,X,Y>&, const M& matrix, const T&)
 {
     return std::make_unique< Dune::MatrixAdapter<M,X,Y> >(matrix);
 }
+
 /**
- * \brief Creates an OverlappingSchwarzOperator as an operator.
+ * \brief Creates an OverlappingSchwarzOperator as an operator, storing it in a unique_ptr.
  *
  * The first argument is used to specify the return type using function overloading.
  * \param matrix The matrix to wrap.
  * \param comm The object encapsulating the parallelization information.
  */
 template<class M, class X, class Y, class T>
-Dune::OverlappingSchwarzOperator<M,X,Y,T> createOperator(const Dune::OverlappingSchwarzOperator<M,X,Y,T>&,
+std::unique_ptr< Dune::OverlappingSchwarzOperator<M,X,Y,T> >
-                                                         const M& matrix, const T& comm)
+createOperator(const Dune::OverlappingSchwarzOperator<M,X,Y,T>&, const M& matrix, const T& comm)
 {
-    return Dune::OverlappingSchwarzOperator<M,X,Y,T>(matrix, comm);
+    return std::make_unique< Dune::OverlappingSchwarzOperator<M,X,Y,T> >(matrix, comm);
 }
 
 //! \brief Applies diagonal scaling to the discretization Matrix (Scheichl, 2003)
@@ -122,10 +114,9 @@ Dune::OverlappingSchwarzOperator<M,X,Y,T> createOperator(const Dune::Overlapping
 //! \param comm The communication objecte describing the data distribution.
 //! \param pressureEqnIndex The index of the pressure in the matrix block
 //! \retun A pair of the scaled matrix and the associated operator-
-template<class Operator, class Communication, class Vector>
+template<class Operator, class Vector>
-std::tuple<std::unique_ptr<typename Operator::matrix_type>, Operator>
+std::unique_ptr<typename Operator::matrix_type>
-scaleMatrixDRS(const Operator& op, const Communication& comm,
+scaleMatrixDRS(const Operator& op, std::size_t pressureEqnIndex, const Vector& weights, const Opm::CPRParameter& param)
-               std::size_t pressureEqnIndex, const Vector& weights, const Opm::CPRParameter& param)
 {
     using Matrix = typename Operator::matrix_type;
     using Block = typename Matrix::block_type;
@@ -144,7 +135,7 @@ scaleMatrixDRS(const Operator& op, const Communication& comm,
             }
         }
     }
-    return std::make_tuple(std::move(matrix), createOperator(op, *matrix, comm));
+    return matrix;
 }
 
 //! \brief Applies diagonal scaling to the discretization Matrix (Scheichl, 2003)
@@ -356,7 +347,7 @@ public:
     /** @brief The type of the arguments used for constructing the smoother. */
     typedef typename Dune::Amg::SmootherTraits<S>::Arguments SmootherArgs;
     /** @brief The type of the AMG construct on the coarse level.*/
-    typedef Dune::Amg::AMG<Operator,X,Smoother,Communication> AMGType;
+    typedef Dune::Amg::AMGCPR<Operator,X,Smoother,Communication> AMGType;
     /**
      * @brief Constructs the coarse solver policy.
      * @param args The arguments used for constructing the smoother.
@@ -384,7 +375,7 @@ private:
                            const Criterion& crit,
                            const typename AMGType::SmootherArgs& args,
                            const Communication& comm)
-            : param_(param), amg_(), smoother_(), op_(op), comm_(comm)
+            : param_(param), amg_(), smoother_(), op_(op), crit_(crit), comm_(comm)
         {
             if ( param_->cpr_use_amg_ )
             {
@@ -400,6 +391,11 @@ private:
             }
         }
 
+        void updateAmgPreconditioner(typename AMGType::Operator& op)
+        {
+            amg_->updateSolver(crit_, op, comm_);
+        }
+
 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
         Dune::SolverCategory::Category category() const override
         {
@@ -408,7 +404,7 @@ private:
         }
 #endif
 
-        void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res)
+        void apply(X& x, X& b, double reduction, Dune::InverseOperatorResult& res) override
         {
             DUNE_UNUSED_PARAMETER(reduction);
             DUNE_UNUSED_PARAMETER(res);
@@ -518,7 +514,7 @@ private:
 #endif
         }
 
-        void apply(X& x, X& b, Dune::InverseOperatorResult& res)
+        void apply(X& x, X& b, Dune::InverseOperatorResult& res) override
         {
             return apply(x,b,1e-8,res);
         }
@@ -535,6 +531,7 @@ private:
         std::unique_ptr<AMGType> amg_;
         std::unique_ptr<Smoother> smoother_;
         const typename AMGType::Operator& op_;
+        Criterion crit_;
         const Communication& comm_;
     };
 
@@ -560,8 +557,13 @@ public:
                                                          smootherArgs_,
                                                          transfer.getCoarseLevelCommunication());
 
-        return inv; //std::shared_ptr<InverseOperator<X,X> >(inv);
+        return inv;
+    }
 
+    template<class LTP>
+    void setCoarseOperator(LTP& transferPolicy)
+    {
+        coarseOperator_= transferPolicy.getCoarseLevelOperator();
     }
 
 private:
@@ -706,7 +708,7 @@ void buildCoarseSparseMatrix(M& coarseMatrix, G& fineGraph, const V& visitedMap,
  */
 template<class Operator, class Criterion, class Communication, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
 class OneComponentAggregationLevelTransferPolicy
-    : public Dune::Amg::LevelTransferPolicy<Operator, typename Detail::ScalarType<Operator>::value>
+    : public Dune::Amg::LevelTransferPolicyCpr<Operator, typename Detail::ScalarType<Operator>::value>
 {
     typedef Dune::Amg::AggregatesMap<typename Operator::matrix_type::size_type> AggregatesMap;
 public:
@@ -789,7 +791,7 @@ public:
             {
                 coarseLevelCommunication_->freeGlobalLookup();
             }
-            calculateCoarseEntries(fineOperator.getmat());
+            calculateCoarseEntriesWithAggregatesMap(fineOperator);
         }
         else
         {
@@ -829,10 +831,10 @@ public:
         this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
     }
 
-    template<class M>
+    void calculateCoarseEntriesWithAggregatesMap(const Operator& fineOperator)
-    void calculateCoarseEntries(const M& fineMatrix)
     {
-      *coarseLevelMatrix_ = 0;
+        const auto& fineMatrix = fineOperator.getmat();
+        *coarseLevelMatrix_ = 0;
         for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end();
             row != rowEnd; ++row)
         {
@@ -852,6 +854,23 @@ public:
         }
     }
 
+    virtual void calculateCoarseEntries(const Operator& fineOperator)
+    {
+        const auto& fineMatrix = fineOperator.getmat();
+        *coarseLevelMatrix_ = 0;
+        for(auto row = fineMatrix.begin(), rowEnd = fineMatrix.end();
+            row != rowEnd; ++row)
+            {
+                const auto& i = row.index();
+                for(auto entry = row->begin(), entryEnd = row->end();
+                    entry != entryEnd; ++entry)
+                    {
+                        const auto& j = entry.index();
+                        (*coarseLevelMatrix_)[i][j] += (*entry)[COMPONENT_INDEX][VARIABLE_INDEX];
+                    }
+            }
+    }
+
     void moveToCoarseLevel(const typename FatherType::FineRangeType& fine)
     {
         // Set coarse vector to zero
@@ -980,9 +999,9 @@ protected:
                                              CoarseCriterion,
                                              LevelTransferPolicy>;
     using TwoLevelMethod =
-        Dune::Amg::TwoLevelMethod<Operator,
+        Dune::Amg::TwoLevelMethodCpr<Operator,
-                                  CoarseSolverPolicy,
+                                     CoarseSolverPolicy,
-                                  Smoother>;
+                                     Smoother>;
 public:
 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
     Dune::SolverCategory::Category category() const override
@@ -1011,29 +1030,29 @@ public:
                 const SmootherArgs& smargs, const Communication& comm)
         : param_(param),
           weights_(weights),
-          scaledMatrixOperator_(Detail::scaleMatrixDRS(fineOperator, comm,
+          scaledMatrix_(Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights, param)),
-                                                       COMPONENT_INDEX, weights, param)),
+          scaledMatrixOperator_(Detail::createOperator(fineOperator, *scaledMatrix_, comm)),
-          smoother_(Detail::constructSmoother<Smoother>(std::get<1>(scaledMatrixOperator_),
+          smoother_(Detail::constructSmoother<Smoother>(*scaledMatrixOperator_,
                                                         smargs, comm)),
           levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_),
           coarseSolverPolicy_(&param, smargs, criterion),
-          twoLevelMethod_(std::get<1>(scaledMatrixOperator_), smoother_,
+          twoLevelMethod_(*scaledMatrixOperator_, smoother_,
                           levelTransferPolicy_, coarseSolverPolicy_, 0, 1)
     {}
 
     void pre(typename TwoLevelMethod::FineDomainType& x,
-             typename TwoLevelMethod::FineRangeType& b)
+             typename TwoLevelMethod::FineRangeType& b) override
     {
         twoLevelMethod_.pre(x,b);
     }
 
-    void post(typename TwoLevelMethod::FineDomainType& x)
+    void post(typename TwoLevelMethod::FineDomainType& x) override
     {
         twoLevelMethod_.post(x);
     }
 
     void apply(typename TwoLevelMethod::FineDomainType& v,
-               const typename TwoLevelMethod::FineRangeType& d)
+               const typename TwoLevelMethod::FineRangeType& d) override
     {
         auto scaledD = d;
         Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_);
@@ -1042,7 +1061,8 @@ public:
 private:
     const CPRParameter& param_;
     const typename TwoLevelMethod::FineDomainType& weights_;
-    std::tuple<std::unique_ptr<Matrix>, Operator> scaledMatrixOperator_;
+    std::unique_ptr<Matrix> scaledMatrix_;
+    std::unique_ptr<Operator> scaledMatrixOperator_;
     std::shared_ptr<Smoother> smoother_;
     LevelTransferPolicy levelTransferPolicy_;
     CoarseSolverPolicy coarseSolverPolicy_;

opm/autodiff/BlackoilAmgCpr.hpp

@@ -0,0 +1,181 @@
+/*
+  Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
+
+  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_AMGCPR_HEADER_INCLUDED
+#define OPM_AMGCPR_HEADER_INCLUDED
+
+#include <opm/autodiff/twolevelmethodcpr.hh>
+#include <ewoms/linear/matrixblock.hh>
+#include <opm/autodiff/ParallelOverlappingILU0.hpp>
+#include <opm/autodiff/FlowLinearSolverParameters.hpp>
+#include <opm/autodiff/CPRPreconditioner.hpp>
+#include <opm/autodiff/amgcpr.hh>
+#include <dune/istl/paamg/twolevelmethod.hh>
+#include <dune/istl/paamg/aggregates.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/schwarz.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/scalarproducts.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+namespace Opm
+{
+
+    /**
+     * \brief An algebraic twolevel or multigrid approach for solving blackoil (supports CPR with and without AMG)
+     *
+     * This preconditioner first decouples the component used for coarsening using a simple scaling
+     * approach (e.g. Scheichl, Masson 2013,\see scaleMatrixDRS). Then it constructs the
+     * coarse level system. The coupling is defined by the weights corresponding to the element located at
+     * (COMPONENT_INDEX, VARIABLE_INDEX) in the block matrix. Then the coarse level system is constructed
+     * either by extracting these elements, or by applying aggregation to them directly. This coarse level
+     * can be solved either by AMG or by ILU. The preconditioner is configured using CPRParameter.
+     * \tparam O The type of the operator (encapsulating a BCRSMatrix).
+     * \tparam S The type of the smoother.
+     * \tparam C The type of coarsening criterion to use.
+     * \tparam P The type of the class describing the parallelization.
+     * \tparam COMPONENT_INDEX The index of the component to use for coarsening (usually water).
+     * \tparam VARIABLE_INDEX The index of the variable to use for coarsening (usually pressure).
+     */
+    template<typename O, typename S, typename SC, typename C,
+             typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
+    class BlackoilAmgCpr
+        : public Dune::Preconditioner<typename O::domain_type, typename O::range_type>
+    {
+    public:
+        /** \brief The type of the operator (encapsulating a BCRSMatrix). */
+        using Operator = O;
+        /** \brief The type of coarsening criterion to use. */
+        using Criterion = C;
+        /** \brief The type of the class describing the parallelization. */
+        using Communication = P;
+        /** \brief The type of the smoother. */
+        using Smoother = S;
+        /** \brief The type of the smoother arguments for construction. */
+        using SmootherArgs   = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
+
+    protected:
+        using Matrix = typename Operator::matrix_type;
+        using CoarseOperator = typename Detail::ScalarType<Operator>::value;
+        using CoarseSmoother = typename Detail::ScalarType<SC>::value;
+        using FineCriterion  =
+            typename Detail::OneComponentCriterionType<Criterion,COMPONENT_INDEX, VARIABLE_INDEX>::value;
+        using CoarseCriterion =  typename Detail::ScalarType<Criterion>::value;
+        using LevelTransferPolicy =
+            OneComponentAggregationLevelTransferPolicy<Operator,
+                                                       FineCriterion,
+                                                       Communication,
+                                                       COMPONENT_INDEX,
+                                                       VARIABLE_INDEX>;
+        using CoarseSolverPolicy   =
+            Detail::OneStepAMGCoarseSolverPolicy<CoarseOperator,
+                                                 CoarseSmoother,
+                                                 CoarseCriterion,
+                                                 LevelTransferPolicy>;
+        using TwoLevelMethod =
+            Dune::Amg::TwoLevelMethodCpr<Operator,
+                                         CoarseSolverPolicy,
+                                         Smoother>;
+    public:
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+        Dune::SolverCategory::Category category() const override
+        {
+            return std::is_same<Communication, Dune::Amg::SequentialInformation>::value ?
+                Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
+        }
+#else
+        // define the category
+        enum {
+            //! \brief The category the precondtioner is part of.
+            category = Operator::category
+        };
+#endif
+
+        /**
+         * \brief Constructor.
+         * \param param The parameters used for configuring the solver.
+         * \param fineOperator The operator of the fine level.
+         * \param criterion The criterion describing the coarsening approach.
+         * \param smargs The arguments for constructing the smoother.
+         * \param comm The information about the parallelization.
+         */
+        BlackoilAmgCpr(const CPRParameter& param,
+                       const typename TwoLevelMethod::FineDomainType& weights,
+                       const Operator& fineOperator, const Criterion& criterion,
+                       const SmootherArgs& smargs, const Communication& comm)
+            : param_(param),
+              weights_(weights),
+              scaledMatrix_(Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights_, param)),
+              scaledMatrixOperator_(Detail::createOperator(fineOperator, *scaledMatrix_, comm)),
+              smoother_(Detail::constructSmoother<Smoother>(*scaledMatrixOperator_,
+                                                            smargs, comm)),
+              levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_),
+              coarseSolverPolicy_(&param, smargs, criterion),
+              twoLevelMethod_(*scaledMatrixOperator_,
+                              smoother_,
+                              levelTransferPolicy_,
+                              coarseSolverPolicy_, 0, 1)
+        {
+        }
+
+        void updatePreconditioner(const Operator& fineOperator,
+                                  const SmootherArgs& smargs,
+                                  const Communication& comm)
+        {
+            *scaledMatrix_ = *Detail::scaleMatrixDRS(fineOperator, COMPONENT_INDEX, weights_, param_);
+            smoother_.reset(Detail::constructSmoother<Smoother>(*scaledMatrixOperator_, smargs, comm));
+            twoLevelMethod_.updatePreconditioner(*scaledMatrixOperator_,
+                                                 smoother_,
+                                                 coarseSolverPolicy_);
+        }
+    
+        void pre(typename TwoLevelMethod::FineDomainType& x,
+                 typename TwoLevelMethod::FineRangeType& b) override
+        {
+            twoLevelMethod_.pre(x,b);
+        }
+
+        void post(typename TwoLevelMethod::FineDomainType& x) override
+        {
+            twoLevelMethod_.post(x);
+        }
+
+        void apply(typename TwoLevelMethod::FineDomainType& v,
+                   const typename TwoLevelMethod::FineRangeType& d) override
+        {
+            auto scaledD = d;
+            Detail::scaleVectorDRS(scaledD, COMPONENT_INDEX, param_, weights_);
+            twoLevelMethod_.apply(v, scaledD);
+        }
+
+    private:
+        const CPRParameter& param_;
+        const typename TwoLevelMethod::FineDomainType& weights_;
+        std::unique_ptr<Matrix> scaledMatrix_;
+        std::unique_ptr<Operator> scaledMatrixOperator_;
+        std::shared_ptr<Smoother> smoother_;
+        LevelTransferPolicy levelTransferPolicy_;
+        CoarseSolverPolicy coarseSolverPolicy_;
+        TwoLevelMethod twoLevelMethod_;
+    };
+
+} // end namespace Opm
+#endif

opm/autodiff/BlackoilModelEbos.hpp

@@ -298,12 +298,16 @@ namespace Opm {
                 wellModel().linearize(ebosSimulator().model().linearizer().jacobian(),
                                       ebosSimulator().model().linearizer().residual());
 
+                // Solve the linear system.
+                linear_solve_setup_time_ = 0.0;
                 try {
                     solveJacobianSystem(x);
+                    report.linear_solve_setup_time += linear_solve_setup_time_;
                     report.linear_solve_time += perfTimer.stop();
                     report.total_linear_iterations += linearIterationsLastSolve();
                 }
                 catch (...) {
+                    report.linear_solve_setup_time += linear_solve_setup_time_;
                     report.linear_solve_time += perfTimer.stop();
                     report.total_linear_iterations += linearIterationsLastSolve();
 
@@ -473,7 +477,10 @@ namespace Opm {
             x = 0.0;
 
             auto& ebosSolver = ebosSimulator_.model().newtonMethod().linearSolver();
+            Dune::Timer perfTimer;
+            perfTimer.start();
             ebosSolver.prepare(ebosJac, ebosResid);
+            linear_solve_setup_time_ = perfTimer.stop();
             ebosSolver.setResidual(ebosResid);
             // actually, the error needs to be calculated after setResidual in order to
             // account for parallelization properly. since the residual of ECFV
@@ -906,7 +913,7 @@ namespace Opm {
         double dsMax() const { return param_.ds_max_; }
         double drMaxRel() const { return param_.dr_max_rel_; }
         double maxResidualAllowed() const { return param_.max_residual_allowed_; }
-
+        double linear_solve_setup_time_;
     public:
         std::vector<bool> wasSwitched_;
     };

opm/autodiff/ISTLSolverEbos.hpp

@@ -76,6 +76,7 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
 // Implementation for ISTL-matrix based operator
 //=====================================================================
 
+
 /*!
    \brief Adapter to turn a matrix into a linear operator.
 
@@ -129,8 +130,16 @@ public:
     }
 #endif
   }
+  WellModelMatrixAdapter (const M& A,
+                          const M& A_for_precond,
+                          const WellModel& wellMod,
+                          std::shared_ptr<communication_type> comm )
+      : A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_(comm)
+  {
+  }
 
-  virtual void apply( const X& x, Y& y ) const
+
+  virtual void apply( const X& x, Y& y ) const override
   {
     A_.mv( x, y );
 
@@ -144,7 +153,7 @@ public:
   }
 
   // y += \alpha * A * x
-  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
   {
     A_.usmv(alpha,x,y);
 
@@ -157,18 +166,18 @@ public:
 #endif
   }
 
-  virtual const matrix_type& getmat() const { return A_for_precond_; }
+  virtual const matrix_type& getmat() const override { return A_for_precond_; }
 
-  communication_type* comm()
+    std::shared_ptr<communication_type> comm()
   {
-      return comm_.operator->();
+      return comm_;
   }
 
 protected:
   const matrix_type& A_ ;
   const matrix_type& A_for_precond_ ;
-  const WellModel& wellMod_;
+  const WellModel& wellMod_;    
-  std::unique_ptr< communication_type > comm_;
+  std::shared_ptr< communication_type > comm_;
 };
 
     /// This class solves the fully implicit black-oil system by
@@ -178,6 +187,7 @@ protected:
     template <class TypeTag>
     class ISTLSolverEbos
     {
+    protected:
         typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
         typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
         typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;

opm/autodiff/ISTLSolverEbosCpr.hpp

@@ -0,0 +1,319 @@
+/*
+  Copyright 2016 IRIS 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_ISTLSOLVERCPR_EBOS_HEADER_INCLUDED
+#define OPM_ISTLSOLVERCPR_EBOS_HEADER_INCLUDED
+
+#include <opm/autodiff/ISTLSolverEbos.hpp>
+#include <opm/autodiff/BlackoilAmgCpr.hpp>
+#include <utility>
+#include <memory>
+
+namespace Opm
+{
+//=====================================================================
+// Implementation for ISTL-matrix based operator
+//=====================================================================
+
+
+    /// This class solves the fully implicit black-oil system by
+    /// solving the reduced system (after eliminating well variables)
+    /// as a block-structured matrix (one block for all cell variables) for a fixed
+    /// number of cell variables np .
+    /// \tparam MatrixBlockType The type of the matrix block used.
+    /// \tparam VectorBlockType The type of the vector block used.
+    /// \tparam pressureIndex The index of the pressure component in the vector
+    ///                       vector block. It is used to guide the AMG coarsening.
+    ///                       Default is zero.
+    template <class TypeTag>
+    class ISTLSolverEbosCpr : public ISTLSolverEbos<TypeTag>
+    {
+    protected:
+        // Types and indices from superclass.
+        using SuperClass = ISTLSolverEbos<TypeTag>;
+        using Matrix = typename SuperClass::Matrix;
+        using Vector = typename SuperClass::Vector;
+        using WellModel = typename SuperClass::WellModel;
+        using Simulator = typename SuperClass::Simulator;
+        using SparseMatrixAdapter = typename SuperClass::SparseMatrixAdapter;
+        enum { pressureEqnIndex = SuperClass::pressureEqnIndex };
+        enum { pressureVarIndex = SuperClass::pressureVarIndex };
+
+        // New properties in this subclass.
+        using Preconditioner            = Dune::Preconditioner<Vector, Vector>;
+        using MatrixAdapter             = Dune::MatrixAdapter<Matrix,Vector, Vector>;
+
+        using CouplingMetric            = Opm::Amg::Element<pressureEqnIndex,pressureVarIndex>;
+        using CritBase                  = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
+        using Criterion                 = Dune::Amg::CoarsenCriterion<CritBase>;
+
+        using ParallelMatrixAdapter     = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, Dune::OwnerOverlapCopyCommunication<int,int> >;
+        using CprSmootherFine           = Opm::ParallelOverlappingILU0<Matrix, Vector, Vector, Dune::Amg::SequentialInformation>;
+        using CprSmootherCoarse         = CprSmootherFine;
+        using BlackoilAmgType           = BlackoilAmgCpr<MatrixAdapter,CprSmootherFine, CprSmootherCoarse, Criterion, Dune::Amg::SequentialInformation,
+                                                 pressureEqnIndex, pressureVarIndex>;
+        using ParallelCprSmootherFine   = Opm::ParallelOverlappingILU0<Matrix, Vector, Vector, Dune::OwnerOverlapCopyCommunication<int,int> >;
+        using ParallelCprSmootherCoarse = ParallelCprSmootherFine;
+        using ParallelBlackoilAmgType   = BlackoilAmgCpr<ParallelMatrixAdapter, ParallelCprSmootherFine, ParallelCprSmootherCoarse, Criterion,
+                                                 Dune::OwnerOverlapCopyCommunication<int,int>, pressureEqnIndex, pressureVarIndex>;
+
+        using OperatorSerial = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false>;
+        using OperatorParallel = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true>;
+
+    public:
+        static void registerParameters()
+        {
+            FlowLinearSolverParameters::registerParameters<TypeTag>();
+        }
+
+        /// Construct a system solver.
+        /// \param[in] parallelInformation In the case of a parallel run
+        ///                                with dune-istl the information about the parallelization.
+        explicit ISTLSolverEbosCpr(const Simulator& simulator)
+            : SuperClass(simulator), oldMat()
+        {
+            extractParallelGridInformationToISTL(this->simulator_.vanguard().grid(), this->parallelInformation_);
+            detail::findOverlapRowsAndColumns(this->simulator_.vanguard().grid(), this->overlapRowAndColumns_);
+        }
+
+        void prepare(const SparseMatrixAdapter& M, Vector& b)
+        {
+            if (oldMat != nullptr)
+                std::cout << "old was "<<oldMat<<" new is "<<&M.istlMatrix()<<std::endl;
+            oldMat = &M.istlMatrix();
+            int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
+            if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
+                SuperClass::matrix_.reset(new Matrix(M.istlMatrix()));
+            } else {
+                *SuperClass::matrix_ = M.istlMatrix();
+            }
+            SuperClass::rhs_ = &b;
+            SuperClass::scaleSystem();
+            const WellModel& wellModel = this->simulator_.problem().wellModel();
+
+#if HAVE_MPI
+            if( this->isParallel() ) {
+
+                //remove ghost rows in local matrix without doing a copy.
+                this->makeOverlapRowsInvalid(*(this->matrix_));
+
+                if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
+                    //Not sure what actual_mat_for_prec is, so put ebosJacIgnoreOverlap as both variables
+                    //to be certain that correct matrix is used for preconditioning.
+                    if( ! comm_ )
+                    {
+                        opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel,
+                                                                this->parallelInformation_ ));
+                        comm_ = opAParallel_->comm();
+                        assert(comm_->indexSet().size()==0);
+                        const size_t size = opAParallel_->getmat().N();
+
+                        const ParallelISTLInformation& info =
+                            boost::any_cast<const ParallelISTLInformation&>( this->parallelInformation_);
+
+                        // As we use a dune-istl with block size np the number of components
+                        // per parallel is only one.
+                        info.copyValuesTo(comm_->indexSet(), comm_->remoteIndices(),
+                                          size, 1);
+                    }
+                    else
+                    {
+                        opAParallel_.reset(new OperatorParallel(*(this->matrix_), *(this->matrix_), wellModel,
+                                                                comm_ ));
+                    }
+                }
+
+                using POrComm =  Dune::OwnerOverlapCopyCommunication<int,int>;
+
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+                constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::overlapping;
+                auto sp = Dune::createScalarProduct<Vector,POrComm>(*comm_, category);
+                sp_ = std::move(sp);
+#else
+                constexpr int  category = Dune::SolverCategory::overlapping;
+                typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
+                typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
+                SPPointer sp(ScalarProductChooser::construct(*comm_));
+                sp_ = std::move(sp);
+#endif
+
+                using AMGOperator = Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector, POrComm>;
+                // If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>;
+                if( ! std::is_same< OperatorParallel, AMGOperator > :: value &&
+                    ( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) ) {
+                    // create new operator in case linear operator and matrix operator differ
+                    opA_.reset( new AMGOperator( opAParallel_->getmat(), *comm_ ));
+                }
+
+                prepareSolver(*opAParallel_, *comm_);
+
+            } else
+#endif
+            {
+
+                if (newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_)) {
+                    opASerial_.reset(new OperatorSerial(*(this->matrix_), *(this->matrix_), wellModel));
+                }
+
+                using POrComm = Dune::Amg::SequentialInformation;
+                POrComm parallelInformation_arg;
+                typedef  OperatorSerial LinearOperator;
+
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+                constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::sequential;
+                auto sp = Dune::createScalarProduct<Vector,POrComm>(parallelInformation_arg, category);
+                sp_ = std::move(sp);
+#else
+                constexpr int  category = Dune::SolverCategory::sequential;
+                typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
+                typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
+                SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
+                sp_ = std::move(sp);
+#endif
+
+                // If clause is always execute as as Linearoperator is WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false|true>;
+                if( ! std::is_same< LinearOperator, MatrixAdapter > :: value &&
+                    ( newton_iteration < 1 or not(this->parameters_.cpr_reuse_setup_) ) ) {
+                    // create new operator in case linear operator and matrix operator differ
+                    opA_.reset( new MatrixAdapter( opASerial_->getmat()));//, parallelInformation_arg ) );
+                }
+
+                prepareSolver(*opASerial_, parallelInformation_arg);
+            }
+        }
+
+        template<typename Operator, typename Comm>
+        void prepareSolver(Operator& wellOpA, Comm& comm)
+        {
+
+            Vector& istlb = *(this->rhs_);
+            comm.copyOwnerToAll(istlb, istlb);
+
+            const double relax = this->parameters_.ilu_relaxation_;
+            const MILU_VARIANT ilu_milu  = this->parameters_.ilu_milu_;
+            using Matrix         = typename MatrixAdapter::matrix_type;
+            const int verbosity    = ( this->parameters_.cpr_solver_verbose_ &&
+                                       comm.communicator().rank()==0 ) ? 1 : 0;
+
+            // TODO: revise choice of parameters
+            // int coarsenTarget = 4000;
+            int coarsenTarget = 1200;
+            Criterion criterion(15, coarsenTarget);
+            criterion.setDebugLevel( this->parameters_.cpr_solver_verbose_ ); // no debug information, 1 for printing hierarchy information
+            criterion.setDefaultValuesIsotropic(2);
+            criterion.setNoPostSmoothSteps( 1 );
+            criterion.setNoPreSmoothSteps( 1 );
+            //new guesses by hmbn
+            //criterion.setAlpha(0.01); // criterion for connection strong 1/3 is default
+            //criterion.setMaxLevel(2); //
+            //criterion.setGamma(1); //  //1 V cycle 2 WW
+
+            // Since DUNE 2.2 we also need to pass the smoother args instead of steps directly
+            using AmgType           = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
+                                                                BlackoilAmgType, ParallelBlackoilAmgType>::type;
+            using SpType            = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
+                                                                Dune::SeqScalarProduct<Vector>,
+                                                                Dune::OverlappingSchwarzScalarProduct<Vector, Comm> >::type;
+            using OperatorType      = typename std::conditional<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
+                                                                MatrixAdapter, ParallelMatrixAdapter>::type;
+            typedef typename AmgType::Smoother Smoother;
+            typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments  SmootherArgs;
+            SmootherArgs  smootherArgs;
+            smootherArgs.iterations = 1;
+            smootherArgs.relaxationFactor = relax;
+            const Opm::CPRParameter& params(this->parameters_); // strange conversion
+            ISTLUtility::setILUParameters(smootherArgs, ilu_milu);
+
+            auto& opARef = reinterpret_cast<OperatorType&>(*opA_);
+            int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
+            double dt = this->simulator_.timeStepSize();
+            bool update_preconditioner = false;
+
+            if (this->parameters_.cpr_reuse_setup_ < 1) {
+                update_preconditioner = true;
+            }
+            if (this->parameters_.cpr_reuse_setup_ < 2) {
+                if (newton_iteration < 1) {
+                    update_preconditioner = true;
+                }
+            }
+            if (this->parameters_.cpr_reuse_setup_ < 3) {
+                if (this->iterations() > 10) {
+                    update_preconditioner = true;
+                }
+            }
+
+            if ( update_preconditioner or (amg_== 0) ) {
+                amg_.reset( new AmgType( params, this->weights_, opARef, criterion, smootherArgs, comm ) );
+            } else {
+                if (this->parameters_.cpr_solver_verbose_) {
+                    std::cout << " Only update amg solver " << std::endl;
+                }
+                reinterpret_cast<AmgType*>(amg_.get())->updatePreconditioner(opARef, smootherArgs, comm);
+            }
+            // Solve.
+            //SuperClass::solve(linearOperator, x, istlb, *sp, *amg, result);
+            //references seems to do something els than refering
+
+            int verbosity_linsolve = 0;
+            if (comm.communicator().rank() == 0) {
+                verbosity_linsolve = this->parameters_.linear_solver_verbosity_;
+            }
+
+            linsolve_.reset(new Dune::BiCGSTABSolver<Vector>(wellOpA, reinterpret_cast<SpType&>(*sp_), reinterpret_cast<AmgType&>(*amg_),
+                                                             this->parameters_.linear_solver_reduction_,
+                                                             this->parameters_.linear_solver_maxiter_,
+                                                             verbosity_linsolve));
+        }
+
+        bool solve(Vector& x)
+        {
+            // Solve system.
+            Dune::InverseOperatorResult result;
+            Vector& istlb = *(this->rhs_);
+            linsolve_->apply(x, istlb, result);
+            SuperClass::checkConvergence(result);
+            if (this->parameters_.scale_linear_system_) {
+                this->scaleSolution(x);
+            }
+	    return this->converged_;
+        }
+
+
+    protected:
+
+      ///! \brief The dune-istl operator (either serial or parallel
+      std::unique_ptr< Dune::LinearOperator<Vector, Vector> > opA_;
+      ///! \brief Serial well matrix adapter
+      std::unique_ptr< OperatorSerial > opASerial_;
+      ///! \brief Parallel well matrix adapter
+      std::unique_ptr< OperatorParallel > opAParallel_;
+      ///! \brief The preconditoner to use (either serial or parallel CPR with AMG)
+      std::unique_ptr< Preconditioner > amg_;
+        
+      using SPPointer = std::shared_ptr< Dune::ScalarProduct<Vector> >;
+      SPPointer sp_;
+      std::shared_ptr< Dune::BiCGSTABSolver<Vector> > linsolve_;
+      const void* oldMat;
+      using POrComm =  Dune::OwnerOverlapCopyCommunication<int,int>;
+      std::shared_ptr<POrComm> comm_;
+    }; // end ISTLSolver
+
+} // namespace Opm
+#endif

opm/autodiff/ParallelOverlappingILU0.hpp

@@ -724,7 +724,7 @@ public:
 
       \copydoc Preconditioner::pre(X&,Y&)
     */
-    virtual void pre (Domain& x, Range& b)
+    virtual void pre (Domain& x, Range& b) override
     {
         DUNE_UNUSED_PARAMETER(x);
         DUNE_UNUSED_PARAMETER(b);
@@ -735,7 +735,7 @@ public:
 
       \copydoc Preconditioner::apply(X&,const Y&)
     */
-    virtual void apply (Domain& v, const Range& d)
+    virtual void apply (Domain& v, const Range& d) override
     {
         Range& md = reorderD(d);
         Domain& mv = reorderV(v);
@@ -806,7 +806,7 @@ public:
 
       \copydoc Preconditioner::post(X&)
     */
-    virtual void post (Range& x)
+    virtual void post (Range& x) override
     {
         DUNE_UNUSED_PARAMETER(x);
     }

opm/autodiff/amgcpr.hh

@@ -0,0 +1,938 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_AMG_AMG_CPR_HH
+#define DUNE_AMG_AMG_CPR_HH
+
+// NOTE: This file is a modified version of dune/istl/paamg/amg.hh from
+// dune-istl release 2.6.0. Modifications have been kept as minimal as possible.
+
+#include <memory>
+#include <dune/common/exceptions.hh>
+#include <dune/istl/paamg/smoother.hh>
+#include <dune/istl/paamg/transfer.hh>
+#include <dune/istl/paamg/hierarchy.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/scalarproducts.hh>
+#include <dune/istl/superlu.hh>
+#include <dune/istl/umfpack.hh>
+#include <dune/istl/solvertype.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/common/exceptions.hh>
+
+namespace Dune
+{
+  namespace Amg
+  {
+    /**
+     * @defgroup ISTL_PAAMG Parallel Algebraic Multigrid
+     * @ingroup ISTL_Prec
+     * @brief A Parallel Algebraic Multigrid based on Agglomeration.
+     */
+
+    /**
+     * @addtogroup ISTL_PAAMG
+     *
+     * @{
+     */
+
+    /** @file
+     * @author Markus Blatt
+     * @brief The AMG preconditioner.
+     */
+
+    template<class M, class X, class S, class P, class K, class A>
+    class KAMG;
+
+    template<class T>
+    class KAmgTwoGrid;
+
+    /**
+     * @brief Parallel algebraic multigrid based on agglomeration.
+     *
+     * \tparam M The matrix type
+     * \tparam X The vector type
+     * \tparam S The smoother type
+     * \tparam A An allocator for X
+     *
+     * \todo drop the smoother template parameter and replace with dynamic construction
+     */
+    template<class M, class X, class S, class PI=SequentialInformation,
+        class A=std::allocator<X> >
+    class AMGCPR : public Preconditioner<X,X>
+    {
+      template<class M1, class X1, class S1, class P1, class K1, class A1>
+      friend class KAMG;
+
+      friend class KAmgTwoGrid<AMGCPR>;
+
+    public:
+      /** @brief The matrix operator type. */
+      typedef M Operator;
+      /**
+       * @brief The type of the parallel information.
+       * Either OwnerOverlapCommunication or another type
+       * describing the parallel data distribution and
+       * providing communication methods.
+       */
+      typedef PI ParallelInformation;
+      /** @brief The operator hierarchy type. */
+      typedef MatrixHierarchy<M, ParallelInformation, A> OperatorHierarchy;
+      /** @brief The parallal data distribution hierarchy type. */
+      typedef typename OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy;
+
+      /** @brief The domain type. */
+      typedef X Domain;
+      /** @brief The range type. */
+      typedef X Range;
+      /** @brief the type of the coarse solver. */
+      typedef InverseOperator<X,X> CoarseSolver;
+      /**
+       * @brief The type of the smoother.
+       *
+       * One of the preconditioners implementing the Preconditioner interface.
+       * Note that the smoother has to fit the ParallelInformation.*/
+      typedef S Smoother;
+
+      /** @brief The argument type for the construction of the smoother. */
+      typedef typename SmootherTraits<Smoother>::Arguments SmootherArgs;
+
+      /**
+       * @brief Construct a new amg with a specific coarse solver.
+       * @param matrices The already set up matix hierarchy.
+       * @param coarseSolver The set up solver to use on the coarse
+       * grid, must match the coarse matrix in the matrix hierarchy.
+       * @param smootherArgs The  arguments needed for thesmoother to use
+       * for pre and post smoothing.
+       * @param parms The parameters for the AMG.
+       */
+      AMGCPR(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
+          const SmootherArgs& smootherArgs, const Parameters& parms);
+
+      /**
+       * @brief Construct an AMG with an inexact coarse solver based on the smoother.
+       *
+       * As coarse solver a preconditioned CG method with the smoother as preconditioner
+       * will be used. The matrix hierarchy is built automatically.
+       * @param fineOperator The operator on the fine level.
+       * @param criterion The criterion describing the coarsening strategy. E. g. SymmetricCriterion
+       * or UnsymmetricCriterion, and providing the parameters.
+       * @param smootherArgs The arguments for constructing the smoothers.
+       * @param pinfo The information about the parallel distribution of the data.
+       */
+      template<class C>
+      AMGCPR(const Operator& fineOperator, const C& criterion,
+          const SmootherArgs& smootherArgs=SmootherArgs(),
+          const ParallelInformation& pinfo=ParallelInformation());
+
+      /**
+       * @brief Copy constructor.
+       */
+      AMGCPR(const AMGCPR& amg);
+
+      ~AMGCPR();
+
+      /** \copydoc Preconditioner::pre */
+      void pre(Domain& x, Range& b);
+
+      /** \copydoc Preconditioner::apply */
+      void apply(Domain& v, const Range& d);
+
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+      //! Category of the preconditioner (see SolverCategory::Category)
+      virtual SolverCategory::Category category() const
+      {
+        return category_;
+      }
+#else
+      enum {
+            //! \brief The category the preconditioner is part of.
+            category = std::is_same<PI,Dune::Amg::SequentialInformation>::value?
+            Dune::SolverCategory::sequential:Dune::SolverCategory::overlapping
+        };
+#endif
+
+      /** \copydoc Preconditioner::post */
+      void post(Domain& x);
+
+      /**
+       * @brief Get the aggregate number of each unknown on the coarsest level.
+       * @param cont The random access container to store the numbers in.
+       */
+      template<class A1>
+      void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
+
+      std::size_t levels();
+
+      std::size_t maxlevels();
+
+      /**
+       * @brief Recalculate the matrix hierarchy.
+       *
+       * It is assumed that the coarsening for the changed fine level
+       * matrix would yield the same aggregates. In this case it suffices
+       * to recalculate all the Galerkin products for the matrices of the
+       * coarser levels.
+       */
+      void recalculateHierarchy()
+      {
+        matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
+      }
+
+      /**
+       * @brief Update the coarse solver and the hierarchies.
+       */
+      template<class C>  
+      void updateSolver(C& criterion, Operator& /* matrix */, const PI& pinfo);
+
+      /**
+       * @brief Check whether the coarse solver used is a direct solver.
+       * @return True if the coarse level solver is a direct solver.
+       */
+      bool usesDirectCoarseLevelSolver() const;
+
+    private:
+      /**
+       * @brief Create matrix and smoother hierarchies.
+       * @param criterion The coarsening criterion.
+       * @param matrix The fine level matrix operator.
+       * @param pinfo The fine level parallel information.
+       */
+      template<class C>
+      void createHierarchies(C& criterion, Operator& matrix,
+                             const PI& pinfo);
+
+      void setupCoarseSolver();
+
+      /**
+       * @brief A struct that holds the context of the current level.
+       *
+       * These are the iterators to the smoother, matrix, parallel information,
+       * and so on needed for the computations on the current level.
+       */
+      struct LevelContext
+      {
+        typedef Smoother SmootherType;
+        /**
+         * @brief The iterator over the smoothers.
+         */
+        typename Hierarchy<Smoother,A>::Iterator smoother;
+        /**
+         * @brief The iterator over the matrices.
+         */
+        typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix;
+        /**
+         * @brief The iterator over the parallel information.
+         */
+        typename ParallelInformationHierarchy::Iterator pinfo;
+        /**
+         * @brief The iterator over the redistribution information.
+         */
+        typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
+        /**
+         * @brief The iterator over the aggregates maps.
+         */
+        typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
+        /**
+         * @brief The iterator over the left hand side.
+         */
+        typename Hierarchy<Domain,A>::Iterator lhs;
+        /**
+         * @brief The iterator over the updates.
+         */
+        typename Hierarchy<Domain,A>::Iterator update;
+        /**
+         * @brief The iterator over the right hand sided.
+         */
+        typename Hierarchy<Range,A>::Iterator rhs;
+        /**
+         * @brief The level index.
+         */
+        std::size_t level;
+      };
+
+
+      /**
+       * @brief Multigrid cycle on a level.
+       * @param levelContext the iterators of the current level.
+       */
+      void mgc(LevelContext& levelContext);
+
+      void additiveMgc();
+
+      /**
+       * @brief Move the iterators to the finer level
+       * @param levelContext the iterators of the current level.
+       * @param processedFineLevel Whether the process computed on
+       *         fine level or not.
+       */
+      void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
+
+      /**
+       * @brief Move the iterators to the coarser level.
+       * @param levelContext the iterators of the current level
+       */
+      bool moveToCoarseLevel(LevelContext& levelContext);
+
+      /**
+       * @brief Initialize iterators over levels with fine level.
+       * @param levelContext the iterators of the current level
+       */
+      void initIteratorsWithFineLevel(LevelContext& levelContext);
+
+      /**  @brief The matrix we solve. */
+      std::shared_ptr<OperatorHierarchy> matrices_;
+      /** @brief The arguments to construct the smoother */
+      SmootherArgs smootherArgs_;
+      /** @brief The hierarchy of the smoothers. */
+      std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
+      /** @brief The solver of the coarsest level. */
+      std::shared_ptr<CoarseSolver> solver_;
+      /** @brief The right hand side of our problem. */
+      Hierarchy<Range,A>* rhs_;
+      /** @brief The left approximate solution of our problem. */
+      Hierarchy<Domain,A>* lhs_;
+      /** @brief The total update for the outer solver. */
+      Hierarchy<Domain,A>* update_;
+      /** @brief The type of the scalar product for the coarse solver. */
+      using ScalarProduct = Dune::ScalarProduct<X>;
+      /** @brief Scalar product on the coarse level. */
+      std::shared_ptr<ScalarProduct> scalarProduct_;
+      /** @brief Gamma, 1 for V-cycle and 2 for W-cycle. */
+      std::size_t gamma_;
+      /** @brief The number of pre and postsmoothing steps. */
+      std::size_t preSteps_;
+      /** @brief The number of postsmoothing steps. */
+      std::size_t postSteps_;
+      bool buildHierarchy_;
+      bool additive;
+      bool coarsesolverconverged;
+      std::shared_ptr<Smoother> coarseSmoother_;
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+      /** @brief The solver category. */
+      SolverCategory::Category category_;
+#endif
+      /** @brief The verbosity level. */
+      std::size_t verbosity_;
+    };
+
+    template<class M, class X, class S, class PI, class A>
+    inline AMGCPR<M,X,S,PI,A>::AMGCPR(const AMGCPR& amg)
+    : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
+      smoothers_(amg.smoothers_), solver_(amg.solver_),
+      rhs_(), lhs_(), update_(),
+      scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
+      preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
+      buildHierarchy_(amg.buildHierarchy_),
+      additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
+      coarseSmoother_(amg.coarseSmoother_),
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+      category_(amg.category_),
+#endif
+      verbosity_(amg.verbosity_)
+    {
+      if(amg.rhs_)
+        rhs_=new Hierarchy<Range,A>(*amg.rhs_);
+      if(amg.lhs_)
+        lhs_=new Hierarchy<Domain,A>(*amg.lhs_);
+      if(amg.update_)
+        update_=new Hierarchy<Domain,A>(*amg.update_);
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    AMGCPR<M,X,S,PI,A>::AMGCPR(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
+                         const SmootherArgs& smootherArgs,
+                         const Parameters& parms)
+      : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
+        smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
+        rhs_(), lhs_(), update_(), scalarProduct_(0),
+        gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
+        postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
+        additive(parms.getAdditive()), coarsesolverconverged(true),
+        coarseSmoother_(),
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+// #warning should category be retrieved from matrices?
+        category_(SolverCategory::category(*smoothers_->coarsest())),
+#endif
+        verbosity_(parms.debugLevel())
+    {
+      assert(matrices_->isBuilt());
+
+      // build the necessary smoother hierarchies
+      matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    template<class C>
+    AMGCPR<M,X,S,PI,A>::AMGCPR(const Operator& matrix,
+                         const C& criterion,
+                         const SmootherArgs& smootherArgs,
+                         const PI& pinfo)
+      : smootherArgs_(smootherArgs),
+        smoothers_(new Hierarchy<Smoother,A>), solver_(),
+        rhs_(), lhs_(), update_(), scalarProduct_(),
+        gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
+        postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
+        additive(criterion.getAdditive()), coarsesolverconverged(true),
+        coarseSmoother_(),
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+        category_(SolverCategory::category(pinfo)),
+#endif
+        verbosity_(criterion.debugLevel())
+    {
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+      if(SolverCategory::category(matrix) != SolverCategory::category(pinfo))
+        DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!");
+#endif
+      createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
+    }
+
+
+    template<class M, class X, class S, class PI, class A>
+    AMGCPR<M,X,S,PI,A>::~AMGCPR()
+    {
+      if(buildHierarchy_) {
+        if(solver_)
+          solver_.reset();
+        if(coarseSmoother_)
+          coarseSmoother_.reset();
+      }
+      if(lhs_)
+        delete lhs_;
+      lhs_=nullptr;
+      if(update_)
+        delete update_;
+      update_=nullptr;
+      if(rhs_)
+        delete rhs_;
+      rhs_=nullptr;
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    template<class C>
+    void AMGCPR<M,X,S,PI,A>::updateSolver(C& /* criterion */, Operator& /* matrix */, const PI& /* pinfo */)
+    {
+      Timer watch;
+      smoothers_.reset(new Hierarchy<Smoother,A>);
+      solver_.reset();
+      coarseSmoother_.reset();
+      scalarProduct_.reset();
+      buildHierarchy_= true;
+      coarsesolverconverged = true;
+      smoothers_.reset(new Hierarchy<Smoother,A>);
+      matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
+      matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
+      setupCoarseSolver();
+      if (verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) {
+        std::cout << "Recalculating galerkin and coarse somothers "<< matrices_->maxlevels() << " levels "
+                  << watch.elapsed() << " seconds." << std::endl;
+      }
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    template<class C>
+    void AMGCPR<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix,
+                                            const PI& pinfo)
+    {
+      Timer watch;
+      matrices_.reset(new OperatorHierarchy(matrix, pinfo));
+
+      matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
+
+      // build the necessary smoother hierarchies
+      matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
+      setupCoarseSolver();
+      if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
+        std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
+                 <<"(inclusive coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::setupCoarseSolver()
+    {
+      // test whether we should solve on the coarse level. That is the case if we
+      // have that level and if there was a redistribution on this level then our
+      // communicator has to be valid (size()>0) as the smoother might try to communicate
+      // in the constructor.
+      if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
+         && ( ! matrices_->redistributeInformation().back().isSetup() ||
+              matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
+      {
+        // We have the carsest level. Create the coarse Solver
+        SmootherArgs sargs(smootherArgs_);
+        sargs.iterations = 1;
+
+        typename ConstructionTraits<Smoother>::Arguments cargs;
+        cargs.setArgs(sargs);
+        if(matrices_->redistributeInformation().back().isSetup()) {
+          // Solve on the redistributed partitioning
+          cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
+          cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
+        }else{
+          cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
+          cargs.setComm(*matrices_->parallelInformation().coarsest());
+        }
+
+        coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
+
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+        scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
+#else
+        typedef Dune::ScalarProductChooser<X,ParallelInformation,category>
+          ScalarProductChooser;
+        // the scalar product.
+        scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
+#endif
+
+
+        typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
+
+        // Use superlu if we are purely sequential or with only one processor on the coarsest level.
+        if( SolverSelector::isDirectSolver &&
+            (std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
+           || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
+           || (matrices_->parallelInformation().coarsest().isRedistributed()
+               && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
+               && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
+          )
+        { // redistribute and 1 proc
+          if(matrices_->parallelInformation().coarsest().isRedistributed())
+          {
+            if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
+            {
+              // We are still participating on this level
+              solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
+            }
+            else
+              solver_.reset();
+          }
+          else
+          {
+            solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
+          }
+          if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
+            std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
+        }
+        else
+        {
+          if(matrices_->parallelInformation().coarsest().isRedistributed())
+          {
+            if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
+              // We are still participating on this level
+
+              // we have to allocate these types using the rebound allocator
+              // in order to ensure that we fulfill the alignement requirements
+              solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
+                                                  // Cast needed for Dune <=2.5
+                                                  reinterpret_cast<typename
+                                                                   std::conditional<std::is_same<PI,SequentialInformation>::value,
+                                                                                    Dune::SeqScalarProduct<X>,
+                                                                                    Dune::OverlappingSchwarzScalarProduct<X,PI> >::type&>(*scalarProduct_),
+                                                  *coarseSmoother_, 1E-2, 1000, 0));
+            else
+              solver_.reset();
+          }else
+          {
+              solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
+                                                  // Cast needed for Dune <=2.5
+                                                  reinterpret_cast<typename
+                                                                   std::conditional<std::is_same<PI,SequentialInformation>::value,
+                                                                                    Dune::SeqScalarProduct<X>,
+                                                                                    Dune::OverlappingSchwarzScalarProduct<X,PI> >::type&>(*scalarProduct_),
+                                                  *coarseSmoother_, 1E-2, 1000, 0));
+            // // we have to allocate these types using the rebound allocator
+            // // in order to ensure that we fulfill the alignement requirements
+            // using Alloc = typename A::template rebind<BiCGSTABSolver<X>>::other;
+            // Alloc alloc;
+            // auto p = alloc.allocate(1);
+            // alloc.construct(p,
+            //   const_cast<M&>(*matrices_->matrices().coarsest()),
+            //   *scalarProduct_,
+            //   *coarseSmoother_, 1E-2, 1000, 0);
+            // solver_.reset(p,[](BiCGSTABSolver<X>* p){
+            //     Alloc alloc;
+            //     alloc.destroy(p);
+            //     alloc.deallocate(p,1);
+            //   });
+          }
+        }
+      }
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::pre(Domain& x, Range& b)
+    {
+      // Detect Matrix rows where all offdiagonal entries are
+      // zero and set x such that  A_dd*x_d=b_d
+      // Thus users can be more careless when setting up their linear
+      // systems.
+      typedef typename M::matrix_type Matrix;
+      typedef typename Matrix::ConstRowIterator RowIter;
+      typedef typename Matrix::ConstColIterator ColIter;
+      typedef typename Matrix::block_type Block;
+      Block zero;
+      zero=typename Matrix::field_type();
+
+      const Matrix& mat=matrices_->matrices().finest()->getmat();
+      for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
+        bool isDirichlet = true;
+        bool hasDiagonal = false;
+        Block diagonal;
+        for(ColIter col=row->begin(); col!=row->end(); ++col) {
+          if(row.index()==col.index()) {
+            diagonal = *col;
+            hasDiagonal = false;
+          }else{
+            if(*col!=zero)
+              isDirichlet = false;
+          }
+        }
+        if(isDirichlet && hasDiagonal)
+          diagonal.solve(x[row.index()], b[row.index()]);
+      }
+
+      if(smoothers_->levels()>0)
+        smoothers_->finest()->pre(x,b);
+      else
+        // No smoother to make x consistent! Do it by hand
+        matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
+      Range* copy = new Range(b);
+      if(rhs_)
+        delete rhs_;
+      rhs_ = new Hierarchy<Range,A>(copy);
+      Domain* dcopy = new Domain(x);
+      if(lhs_)
+        delete lhs_;
+      lhs_ = new Hierarchy<Domain,A>(dcopy);
+      dcopy = new Domain(x);
+      if(update_)
+        delete update_;
+      update_ = new Hierarchy<Domain,A>(dcopy);
+      matrices_->coarsenVector(*rhs_);
+      matrices_->coarsenVector(*lhs_);
+      matrices_->coarsenVector(*update_);
+
+      // Preprocess all smoothers
+      typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
+      typedef typename Hierarchy<Range,A>::Iterator RIterator;
+      typedef typename Hierarchy<Domain,A>::Iterator DIterator;
+      Iterator coarsest = smoothers_->coarsest();
+      Iterator smoother = smoothers_->finest();
+      RIterator rhs = rhs_->finest();
+      DIterator lhs = lhs_->finest();
+      if(smoothers_->levels()>0) {
+
+        assert(lhs_->levels()==rhs_->levels());
+        assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
+        assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
+
+        if(smoother!=coarsest)
+          for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
+            smoother->pre(*lhs,*rhs);
+        smoother->pre(*lhs,*rhs);
+      }
+
+
+      // The preconditioner might change x and b. So we have to
+      // copy the changes to the original vectors.
+      x = *lhs_->finest();
+      b = *rhs_->finest();
+
+    }
+    template<class M, class X, class S, class PI, class A>
+    std::size_t AMGCPR<M,X,S,PI,A>::levels()
+    {
+      return matrices_->levels();
+    }
+    template<class M, class X, class S, class PI, class A>
+    std::size_t AMGCPR<M,X,S,PI,A>::maxlevels()
+    {
+      return matrices_->maxlevels();
+    }
+
+    /** \copydoc Preconditioner::apply */
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::apply(Domain& v, const Range& d)
+    {
+      LevelContext levelContext;
+
+      if(additive) {
+        *(rhs_->finest())=d;
+        additiveMgc();
+        v=*lhs_->finest();
+      }else{
+        // Init all iterators for the current level
+        initIteratorsWithFineLevel(levelContext);
+
+
+        *levelContext.lhs = v;
+        *levelContext.rhs = d;
+        *levelContext.update=0;
+        levelContext.level=0;
+
+        mgc(levelContext);
+
+        if(postSteps_==0||matrices_->maxlevels()==1)
+          levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
+
+        v=*levelContext.update;
+      }
+
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
+    {
+      levelContext.smoother = smoothers_->finest();
+      levelContext.matrix = matrices_->matrices().finest();
+      levelContext.pinfo = matrices_->parallelInformation().finest();
+      levelContext.redist =
+        matrices_->redistributeInformation().begin();
+      levelContext.aggregates = matrices_->aggregatesMaps().begin();
+      levelContext.lhs = lhs_->finest();
+      levelContext.update = update_->finest();
+      levelContext.rhs = rhs_->finest();
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    bool AMGCPR<M,X,S,PI,A>
+    ::moveToCoarseLevel(LevelContext& levelContext)
+    {
+
+      bool processNextLevel=true;
+
+      if(levelContext.redist->isSetup()) {
+        levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
+                             levelContext.rhs.getRedistributed());
+        processNextLevel = levelContext.rhs.getRedistributed().size()>0;
+        if(processNextLevel) {
+          //restrict defect to coarse level right hand side.
+          typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
+          ++levelContext.pinfo;
+          Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+          ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
+                           static_cast<const Range&>(fineRhs.getRedistributed()),
+                           *levelContext.pinfo);
+        }
+      }else{
+        //restrict defect to coarse level right hand side.
+        typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
+        ++levelContext.pinfo;
+        Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+        ::restrictVector(*(*levelContext.aggregates),
+                         *levelContext.rhs, static_cast<const Range&>(*fineRhs),
+                         *levelContext.pinfo);
+      }
+
+      if(processNextLevel) {
+        // prepare coarse system
+        ++levelContext.lhs;
+        ++levelContext.update;
+        ++levelContext.matrix;
+        ++levelContext.level;
+        ++levelContext.redist;
+
+        if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
+          // next level is not the globally coarsest one
+          ++levelContext.smoother;
+          ++levelContext.aggregates;
+        }
+        // prepare the update on the next level
+        *levelContext.update=0;
+      }
+      return processNextLevel;
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>
+    ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
+    {
+      if(processNextLevel) {
+        if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
+          // previous level is not the globally coarsest one
+          --levelContext.smoother;
+          --levelContext.aggregates;
+        }
+        --levelContext.redist;
+        --levelContext.level;
+        //prolongate and add the correction (update is in coarse left hand side)
+        --levelContext.matrix;
+
+        //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
+        --levelContext.lhs;
+        --levelContext.pinfo;
+      }
+      if(levelContext.redist->isSetup()) {
+        // Need to redistribute during prolongateVector
+        levelContext.lhs.getRedistributed()=0;
+        Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+        ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
+                           levelContext.lhs.getRedistributed(),
+                           matrices_->getProlongationDampingFactor(),
+                           *levelContext.pinfo, *levelContext.redist);
+      }else{
+        *levelContext.lhs=0;
+        Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+        ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
+                           matrices_->getProlongationDampingFactor(),
+                           *levelContext.pinfo);
+      }
+
+
+      if(processNextLevel) {
+        --levelContext.update;
+        --levelContext.rhs;
+      }
+
+      *levelContext.update += *levelContext.lhs;
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    bool AMGCPR<M,X,S,PI,A>::usesDirectCoarseLevelSolver() const
+    {
+      return IsDirectSolver< CoarseSolver>::value;
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::mgc(LevelContext& levelContext){
+      if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
+        // Solve directly
+        InverseOperatorResult res;
+        res.converged=true; // If we do not compute this flag will not get updated
+        if(levelContext.redist->isSetup()) {
+          levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
+          if(levelContext.rhs.getRedistributed().size()>0) {
+            // We are still participating in the computation
+            levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
+                                                    levelContext.rhs.getRedistributed());
+            solver_->apply(levelContext.update.getRedistributed(),
+                           levelContext.rhs.getRedistributed(), res);
+          }
+          levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
+          levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
+        }else{
+          levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
+          solver_->apply(*levelContext.update, *levelContext.rhs, res);
+        }
+
+        if (!res.converged)
+          coarsesolverconverged = false;
+      }else{
+        // presmoothing
+        presmooth(levelContext, preSteps_);
+
+#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
+        bool processNextLevel = moveToCoarseLevel(levelContext);
+
+        if(processNextLevel) {
+          // next level
+          for(std::size_t i=0; i<gamma_; i++)
+            mgc(levelContext);
+        }
+
+        moveToFineLevel(levelContext, processNextLevel);
+#else
+        *lhs=0;
+#endif
+
+        if(levelContext.matrix == matrices_->matrices().finest()) {
+          coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
+          if(!coarsesolverconverged){
+            //DUNE_THROW(MathError, "Coarse solver did not converge");
+          }
+        }
+        // postsmoothing
+        postsmooth(levelContext, postSteps_);
+
+      }
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::additiveMgc(){
+
+      // restrict residual to all levels
+      typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
+      typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
+      typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
+      typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
+
+      for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
+        ++pinfo;
+        Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+        ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
+      }
+
+      // pinfo is invalid, set to coarsest level
+      //pinfo = matrices_->parallelInformation().coarsest
+      // calculate correction for all levels
+      lhs = lhs_->finest();
+      typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
+
+      for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
+        // presmoothing
+        *lhs=0;
+        smoother->apply(*lhs, *rhs);
+      }
+
+      // Coarse level solve
+#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
+      InverseOperatorResult res;
+      pinfo->copyOwnerToAll(*rhs, *rhs);
+      solver_->apply(*lhs, *rhs, res);
+
+      if(!res.converged)
+        DUNE_THROW(MathError, "Coarse solver did not converge");
+#else
+      *lhs=0;
+#endif
+      // Prologate and add up corrections from all levels
+      --pinfo;
+      --aggregates;
+
+      for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
+        Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+        ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
+      }
+    }
+
+
+    /** \copydoc Preconditioner::post */
+    template<class M, class X, class S, class PI, class A>
+    void AMGCPR<M,X,S,PI,A>::post(Domain& x)
+    {
+      DUNE_UNUSED_PARAMETER(x);
+      // Postprocess all smoothers
+      typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
+      typedef typename Hierarchy<Domain,A>::Iterator DIterator;
+      Iterator coarsest = smoothers_->coarsest();
+      Iterator smoother = smoothers_->finest();
+      DIterator lhs = lhs_->finest();
+      if(smoothers_->levels()>0) {
+        if(smoother != coarsest  || matrices_->levels()<matrices_->maxlevels())
+          smoother->post(*lhs);
+        if(smoother!=coarsest)
+          for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
+            smoother->post(*lhs);
+        smoother->post(*lhs);
+      }
+      //delete &(*lhs_->finest());
+      delete lhs_;
+      lhs_=nullptr;
+      //delete &(*update_->finest());
+      delete update_;
+      update_=nullptr;
+      //delete &(*rhs_->finest());
+      delete rhs_;
+      rhs_=nullptr;
+    }
+
+    template<class M, class X, class S, class PI, class A>
+    template<class A1>
+    void AMGCPR<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
+    {
+      matrices_->getCoarsestAggregatesOnFinest(cont);
+    }
+
+  } // end namespace Amg
+} // end namespace Dune
+
+#endif

opm/autodiff/twolevelmethodcpr.hh

@@ -0,0 +1,562 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_ISTL_TWOLEVELMETHODCPR_HH
+#define DUNE_ISTL_TWOLEVELMETHODCPR_HH
+
+// NOTE: This file is a modified version of dune/istl/paamg/twolevelmethod.hh from
+// dune-istl release 2.6.0. Modifications have been kept as minimal as possible.
+
+#include <tuple>
+
+#include<dune/istl/operators.hh>
+//#include "amg.hh"
+//#include"galerkin.hh"
+#include<dune/istl/paamg/amg.hh>
+#include<dune/istl/paamg/galerkin.hh>
+#include<dune/istl/solver.hh>
+
+#include<dune/common/unused.hh>
+
+/**
+ * @addtogroup ISTL_PAAMG
+ * @{
+ * @file
+ * @author Markus Blatt
+ * @brief Algebraic twolevel methods.
+ */
+namespace Dune
+{
+namespace Amg
+{
+
+/**
+ * @brief Abstract base class for transfer between levels and creation
+ * of the coarse level system.
+ *
+ * @tparam FO The type of the linear operator of the finel level system. Has to be
+ * derived from AssembledLinearOperator.
+ * @tparam CO The type of the linear operator of the coarse level system. Has to be
+ * derived from AssembledLinearOperator.
+ */
+template<class FO, class CO>
+class LevelTransferPolicyCpr
+{
+public:
+  /**
+   * @brief The linear operator of the finel level system. Has to be
+   * derived from AssembledLinearOperator.
+   */
+  typedef FO FineOperatorType;
+  /**
+   * @brief The type of the range of the fine level operator.
+   */
+  typedef typename FineOperatorType::range_type FineRangeType;
+  /**
+   * @brief The type of the domain of the fine level operator.
+   */
+  typedef typename FineOperatorType::domain_type FineDomainType;
+  /**
+   * @brief The linear operator of the finel level system. Has to be
+   * derived from AssembledLinearOperator.
+   */
+  typedef CO CoarseOperatorType;
+  /**
+   * @brief The type of the range of the coarse level operator.
+   */
+  typedef typename CoarseOperatorType::range_type CoarseRangeType;
+  /**
+   * @brief The type of the domain of the coarse level operator.
+   */
+  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
+  /**
+   * @brief Get the coarse level operator.
+   * @return A shared pointer to the coarse level system.
+   */
+  std::shared_ptr<CoarseOperatorType>& getCoarseLevelOperator()
+  {
+    return operator_;
+  }
+  /**
+   * @brief Get the coarse level right hand side.
+   * @return The coarse level right hand side.
+   */
+  CoarseRangeType& getCoarseLevelRhs()
+  {
+    return rhs_;
+  }
+
+  /**
+   * @brief Get the coarse level left hand side.
+   * @return The coarse level leftt hand side.
+   */
+  CoarseDomainType& getCoarseLevelLhs()
+  {
+    return lhs_;
+  }
+  /**
+   * @brief Transfers the data to the coarse level.
+   *
+   * Restricts the residual to the right hand side of the
+   * coarse level system and initialies the left hand side
+   * of the coarse level system. These can afterwards be accessed
+   * usinf getCoarseLevelRhs() and getCoarseLevelLhs().
+   * @param fineDefect The current residual of the fine level system.
+   */
+  virtual void moveToCoarseLevel(const FineRangeType& fineRhs)=0;
+  /**
+   * @brief Updates the fine level linear system after the correction
+   * of the coarse levels system.
+   *
+   * After returning from this function the coarse level correction
+   * will have been added to fine level system.
+   * @param[inout] fineLhs The left hand side of the fine level to update
+   * with the coarse level correction.
+   */
+  virtual void moveToFineLevel(FineDomainType& fineLhs)=0;
+  /**
+   * @brief Algebraically creates the coarse level system.
+   *
+   * After returning from this function the coarse level operator
+   * can be accessed using getCoarseLevelOperator().
+   * @param fineOperator The operator of the fine level system.
+   */
+  virtual void createCoarseLevelSystem(const FineOperatorType& fineOperator)=0;
+
+  /**
+   * @brief ???.
+   */
+  virtual void calculateCoarseEntries(const FineOperatorType& fineOperator) = 0;
+
+  /** @brief Clone the current object. */
+  virtual LevelTransferPolicyCpr* clone() const =0;
+
+  /** @brief Destructor. */
+  virtual ~LevelTransferPolicyCpr(){}
+
+  protected:
+  /** @brief The coarse level rhs. */
+  CoarseRangeType rhs_;
+  /** @brief The coarse level lhs. */
+  CoarseDomainType lhs_;
+  /** @brief the coarse level linear operator. */
+  std::shared_ptr<CoarseOperatorType> operator_;
+};
+
+/**
+ * @brief A LeveTransferPolicy that used aggregation to construct the coarse level system.
+ * @tparam O The type of the fine and coarse level operator.
+ * @tparam C The criterion that describes the aggregation procedure.
+ */
+template<class O, class C>
+class AggregationLevelTransferPolicyCpr
+  : public LevelTransferPolicyCpr<O,O>
+{
+  typedef Dune::Amg::AggregatesMap<typename O::matrix_type::size_type> AggregatesMap;
+public:
+  typedef LevelTransferPolicyCpr<O,O> FatherType;
+  typedef C Criterion;
+  typedef SequentialInformation ParallelInformation;
+
+  AggregationLevelTransferPolicyCpr(const Criterion& crit)
+  : criterion_(crit)
+  {}
+
+  void createCoarseLevelSystem(const O& fineOperator)
+  {
+    prolongDamp_ = criterion_.getProlongationDampingFactor();
+    GalerkinProduct<ParallelInformation> productBuilder;
+    typedef typename Dune::Amg::MatrixGraph<const typename O::matrix_type> MatrixGraph;
+    typedef typename Dune::Amg::PropertiesGraph<MatrixGraph,Dune::Amg::VertexProperties,
+      Dune::Amg::EdgeProperties,Dune::IdentityMap,Dune::IdentityMap> PropertiesGraph;
+    MatrixGraph mg(fineOperator.getmat());
+    PropertiesGraph pg(mg,Dune::IdentityMap(),Dune::IdentityMap());
+    typedef NegateSet<typename ParallelInformation::OwnerSet> OverlapFlags;
+
+    aggregatesMap_.reset(new AggregatesMap(pg.maxVertex()+1));
+
+    int noAggregates, isoAggregates, oneAggregates, skippedAggregates;
+
+    std::tie(noAggregates, isoAggregates, oneAggregates, skippedAggregates) =
+       aggregatesMap_->buildAggregates(fineOperator.getmat(), pg, criterion_, true);
+     std::cout<<"no aggregates="<<noAggregates<<" iso="<<isoAggregates<<" one="<<oneAggregates<<" skipped="<<skippedAggregates<<std::endl;
+    // misuse coarsener to renumber aggregates
+    Dune::Amg::IndicesCoarsener<Dune::Amg::SequentialInformation,int> renumberer;
+    typedef std::vector<bool>::iterator Iterator;
+    typedef Dune::IteratorPropertyMap<Iterator, Dune::IdentityMap> VisitedMap;
+    std::vector<bool> excluded(fineOperator.getmat().N(), false);
+    VisitedMap vm(excluded.begin(), Dune::IdentityMap());
+    ParallelInformation pinfo;
+    std::size_t aggregates = renumberer.coarsen(pinfo, pg, vm,
+                                                *aggregatesMap_, pinfo,
+                                                noAggregates);
+    std::vector<bool>& visited=excluded;
+
+    typedef std::vector<bool>::iterator Iterator;
+
+    for(Iterator iter= visited.begin(), end=visited.end();
+        iter != end; ++iter)
+          *iter=false;
+    matrix_.reset(productBuilder.build(mg, vm,
+                                       SequentialInformation(),
+                                       *aggregatesMap_,
+                                       aggregates,
+                                       OverlapFlags()));
+    productBuilder.calculate(fineOperator.getmat(), *aggregatesMap_, *matrix_, pinfo, OverlapFlags());
+    this->lhs_.resize(this->matrix_->M());
+    this->rhs_.resize(this->matrix_->N());
+    this->operator_.reset(new O(*matrix_));
+  }
+
+  void moveToCoarseLevel(const typename FatherType::FineRangeType& fineRhs)
+  {
+    Transfer<std::size_t,typename FatherType::FineRangeType,ParallelInformation>
+      ::restrictVector(*aggregatesMap_, this->rhs_, fineRhs, ParallelInformation());
+    this->lhs_=0;
+  }
+
+  void moveToFineLevel(typename FatherType::FineDomainType& fineLhs)
+  {
+    Transfer<std::size_t,typename FatherType::FineRangeType,ParallelInformation>
+      ::prolongateVector(*aggregatesMap_, this->lhs_, fineLhs,
+                         prolongDamp_, ParallelInformation());
+  }
+
+  AggregationLevelTransferPolicyCpr* clone() const
+  {
+    return new AggregationLevelTransferPolicyCpr(*this);
+  }
+
+private:
+  typename O::matrix_type::field_type prolongDamp_;
+  std::shared_ptr<AggregatesMap> aggregatesMap_;
+  Criterion criterion_;
+  std::shared_ptr<typename O::matrix_type> matrix_;
+};
+
+/**
+ * @brief A policy class for solving the coarse level system using one step of AMG.
+ * @tparam O The type of the linear operator used.
+ * @tparam S The type of the smoother used in AMG.
+ * @tparam C The type of the crition used for the aggregation within AMG.
+ */
+template<class O, class S, class C>
+class OneStepAMGCoarseSolverPolicyCpr
+{
+public:
+  /** @brief The type of the linear operator used. */
+  typedef O Operator;
+  /** @brief The type of the range and domain of the operator. */
+  typedef typename O::range_type X;
+  /** @brief The type of the crition used for the aggregation within AMG.*/
+  typedef C Criterion;
+  /** @brief The type of the smoother used in AMG. */
+  typedef S Smoother;
+  /** @brief The type of the arguments used for constructing the smoother. */
+  typedef typename Dune::Amg::SmootherTraits<S>::Arguments SmootherArgs;
+  /** @brief The type of the AMG construct on the coarse level.*/
+  typedef AMG<Operator,X,Smoother> AMGType;
+  /**
+   * @brief Constructs the coarse solver policy.
+   * @param args The arguments used for constructing the smoother.
+   * @param c The crition used for the aggregation within AMG.
+   */
+  OneStepAMGCoarseSolverPolicyCpr(const SmootherArgs& args, const Criterion& c)
+    : smootherArgs_(args), criterion_(c)
+  {}
+  /** @brief Copy constructor. */
+  OneStepAMGCoarseSolverPolicyCpr(const OneStepAMGCoarseSolverPolicyCpr& other)
+  : coarseOperator_(other.coarseOperator_), smootherArgs_(other.smootherArgs_),
+    criterion_(other.criterion_)
+  {}
+private:
+  /**
+   * @brief A wrapper that makes an inverse operator out of AMG.
+   *
+   * The operator will use one step of AMG to approximately solve
+   * the coarse level system.
+   */
+  struct AMGInverseOperator : public InverseOperator<X,X>
+  {
+    AMGInverseOperator(const typename AMGType::Operator& op,
+                       const Criterion& crit,
+                       const typename AMGType::SmootherArgs& args)
+      : amg_(op, crit,args), first_(true)
+    {}
+
+    void apply(X& x, X& b, double reduction, InverseOperatorResult& res)
+    {
+      DUNE_UNUSED_PARAMETER(reduction);
+      DUNE_UNUSED_PARAMETER(res);
+      if(first_)
+      {
+        amg_.pre(x,b);
+        first_=false;
+        x_=x;
+      }
+      amg_.apply(x,b);
+    }
+
+    void apply(X& x, X& b, InverseOperatorResult& res)
+    {
+      return apply(x,b,1e-8,res);
+    }
+
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+    virtual SolverCategory::Category category() const
+    {
+      return amg_.category();
+    }
+#endif
+    
+    ~AMGInverseOperator()
+    {
+      if(!first_)
+        amg_.post(x_);
+    }
+    AMGInverseOperator(const AMGInverseOperator& other)
+    : x_(other.x_), amg_(other.amg_), first_(other.first_)
+    {
+    }
+  private:
+    X x_;
+    AMGType amg_;
+    bool first_;
+  };
+
+public:
+  /** @brief The type of solver constructed for the coarse level. */
+  typedef AMGInverseOperator CoarseLevelSolver;
+
+  /**
+    * @brief Constructs a coarse level solver.
+    *
+    * @param transferPolicy The policy describing the transfer between levels.
+    * @return A pointer to the constructed coarse level solver.
+    * @tparam P The type of the level transfer policy.
+    */
+  template<class P>
+  CoarseLevelSolver* createCoarseLevelSolver(P& transferPolicy)
+  {
+    coarseOperator_=transferPolicy.getCoarseLevelOperator();
+    AMGInverseOperator* inv = new AMGInverseOperator(*coarseOperator_,
+                                                     criterion_,
+                                                     smootherArgs_);
+
+    return inv; //std::shared_ptr<InverseOperator<X,X> >(inv);
+
+  }
+
+private:
+  /** @brief The coarse level operator. */
+  std::shared_ptr<Operator> coarseOperator_;
+  /** @brief The arguments used to construct the smoother. */
+  SmootherArgs smootherArgs_;
+  /** @brief The coarsening criterion. */
+  Criterion criterion_;
+};
+
+/**
+ * @tparam FO The type of the fine level linear operator.
+ * @tparam CSP The type of the coarse level solver policy.
+ * @tparam S The type of the fine level smoother used.
+ */
+template<class FO, class CSP, class S>
+class TwoLevelMethodCpr :
+    public Preconditioner<typename FO::domain_type, typename FO::range_type>
+{
+public:
+  /** @brief The type of the policy for constructing the coarse level solver. */
+  typedef CSP CoarseLevelSolverPolicy;
+  /** @brief The type of the coarse level solver. */
+  typedef typename CoarseLevelSolverPolicy::CoarseLevelSolver CoarseLevelSolver;
+  /**
+   * @brief The linear operator of the finel level system. Has to be
+   * derived from AssembledLinearOperator.
+   */
+  typedef FO FineOperatorType;
+  /**
+   * @brief The type of the range of the fine level operator.
+   */
+  typedef typename FineOperatorType::range_type FineRangeType;
+  /**
+   * @brief The type of the domain of the fine level operator.
+   */
+  typedef typename FineOperatorType::domain_type FineDomainType;
+  /**
+   * @brief The linear operator of the finel level system. Has to be
+   * derived from AssembledLinearOperator.
+   */
+  typedef typename CSP::Operator CoarseOperatorType;
+  /**
+   * @brief The type of the range of the coarse level operator.
+   */
+  typedef typename CoarseOperatorType::range_type CoarseRangeType;
+  /**
+   * @brief The type of the domain of the coarse level operator.
+   */
+  typedef typename CoarseOperatorType::domain_type CoarseDomainType;
+  /**
+   * @brief The type of the fine level smoother.
+   */
+  typedef S SmootherType;
+  
+  // define the category
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
+  
+#else  
+  enum {
+    //! \brief The category the preconditioner is part of.
+    category=SolverCategory::sequential
+  };
+#endif
+
+  /**
+   * @brief Constructs a two level method.
+   *
+   * @tparam CoarseSolverPolicy The policy for constructing the coarse
+   * solver, e.g. OneStepAMGCoarseSolverPolicy
+   * @param op The fine level operator.
+   * @param smoother The fine level smoother.
+   * @param policy The level transfer policy.
+   * @param coarsePolicy The policy for constructing the coarse level solver.
+   * @param preSteps The number of smoothing steps to apply before the coarse
+   * level correction.
+   * @param preSteps The number of smoothing steps to apply after the coarse
+   * level correction.
+   */
+  TwoLevelMethodCpr(const FineOperatorType& op,
+                    std::shared_ptr<SmootherType> smoother,
+                    const LevelTransferPolicyCpr<FineOperatorType,
+                                                 CoarseOperatorType>& policy,
+                    CoarseLevelSolverPolicy& coarsePolicy,
+                    std::size_t preSteps=1, std::size_t postSteps=1)
+    : operator_(&op), smoother_(smoother),
+      preSteps_(preSteps), postSteps_(postSteps)
+  {
+    policy_ = policy.clone();
+    policy_->createCoarseLevelSystem(*operator_);
+    coarseSolver_=coarsePolicy.createCoarseLevelSolver(*policy_);
+  }
+
+  TwoLevelMethodCpr(const TwoLevelMethodCpr& other)
+  : operator_(other.operator_), coarseSolver_(new CoarseLevelSolver(*other.coarseSolver_)),
+    smoother_(other.smoother_), policy_(other.policy_->clone()),
+    preSteps_(other.preSteps_), postSteps_(other.postSteps_)
+  {}
+
+  ~TwoLevelMethodCpr()
+  {
+    // Each instance has its own policy.
+    delete policy_;
+    delete coarseSolver_;
+  }
+
+  void updatePreconditioner(FineOperatorType& /* op */,
+                            std::shared_ptr<SmootherType> smoother,
+                            CoarseLevelSolverPolicy& coarsePolicy)
+  {
+    //assume new matrix is not reallocated the new precondition should anyway be made
+    smoother_ = smoother;
+    if (coarseSolver_) {
+      policy_->calculateCoarseEntries(*operator_);
+      coarsePolicy.setCoarseOperator(*policy_);
+      coarseSolver_->updateAmgPreconditioner(*(policy_->getCoarseLevelOperator()));
+    } else {
+      // we should probably not be heere
+      policy_->createCoarseLevelSystem(*operator_);
+      coarseSolver_ = coarsePolicy.createCoarseLevelSolver(*policy_);
+    }
+  }
+
+  void pre(FineDomainType& x, FineRangeType& b)
+  {
+    smoother_->pre(x,b);
+  }
+
+  void post(FineDomainType& x)
+  {
+    DUNE_UNUSED_PARAMETER(x);
+  }
+
+  void apply(FineDomainType& v, const FineRangeType& d)
+  {
+    FineDomainType u(v);
+    FineRangeType rhs(d);
+    LevelContext context;
+    SequentialInformation info;
+    context.pinfo=&info;
+    context.lhs=&u;
+    context.update=&v;
+    context.smoother=smoother_;
+    context.rhs=&rhs;
+    context.matrix=operator_;
+    // Presmoothing
+    presmooth(context, preSteps_);
+    //Coarse grid correction
+    policy_->moveToCoarseLevel(*context.rhs);
+    InverseOperatorResult res;
+    coarseSolver_->apply(policy_->getCoarseLevelLhs(), policy_->getCoarseLevelRhs(), res);
+    *context.lhs=0;
+    policy_->moveToFineLevel(*context.lhs);
+    *context.update += *context.lhs;
+    // Postsmoothing
+    postsmooth(context, postSteps_);
+
+  }
+#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)  
+//   //! Category of the preconditioner (see SolverCategory::Category)
+    virtual SolverCategory::Category category() const
+    {
+      return SolverCategory::sequential;
+    }
+#endif
+private:
+  /**
+   * @brief Struct containing the level information.
+   */
+  struct LevelContext
+  {
+    /** @brief The type of the smoother used. */
+    typedef S SmootherType;
+    /** @brief A pointer to the smoother. */
+    std::shared_ptr<SmootherType> smoother;
+    /** @brief The left hand side passed to the and returned by the smoother. */
+    FineDomainType* lhs;
+    /*
+     * @brief The right hand side holding the current residual.
+     *
+     * This is passed to the smoother as the right hand side.
+     */
+    FineRangeType* rhs;
+    /**
+     * @brief The total update calculated by the preconditioner.
+     *
+     * I.e. all update from smoothing and coarse grid correction summed up.
+     */
+    FineDomainType* update;
+    /** @parallel information */
+    SequentialInformation* pinfo;
+    /**
+     * @brief The matrix that we are solving.
+     *
+     * Needed to update the residual.
+     */
+    const FineOperatorType* matrix;
+  };
+  const FineOperatorType* operator_;
+  /** @brief The coarse level solver. */
+  CoarseLevelSolver* coarseSolver_;
+  /** @brief The fine level smoother. */
+  std::shared_ptr<S> smoother_;
+  /** @brief Policy for prolongation, restriction, and coarse level system creation. */
+  LevelTransferPolicyCpr<FO,typename CSP::Operator>* policy_;
+  /** @brief The number of presmoothing steps to apply. */
+  std::size_t preSteps_;
+  /** @brief The number of postsmoothing steps to apply. */
+  std::size_t postSteps_;
+};
+}// end namespace Amg
+}// end namespace Dune
+
+/** @} */
+#endif

opm/core/simulator/SimulatorReport.cpp

@@ -33,6 +33,7 @@ namespace Opm
           total_time(0.0),
           solver_time(0.0),
           assemble_time(0.0),
+	  linear_solve_setup_time(0.0),
           linear_solve_time(0.0),
           update_time(0.0),
           output_write_time(0.0),
@@ -49,6 +50,7 @@ namespace Opm
     {
         pressure_time += sr.pressure_time;
         transport_time += sr.transport_time;
+	linear_solve_setup_time += sr.linear_solve_setup_time;
         linear_solve_time += sr.linear_solve_time;
         solver_time += sr.solver_time;
         assemble_time += sr.assemble_time;
@@ -119,6 +121,14 @@ namespace Opm
                 }
                 os << std::endl;
 
+		t = linear_solve_setup_time + (failureReport ? failureReport->linear_solve_setup_time : 0.0);
+                os << " Linear solve setup time (seconds): " << t;
+                if (failureReport) {
+		  os << " (Failed: " << failureReport->linear_solve_setup_time << "; "
+		     << 100*failureReport->linear_solve_setup_time/t << "%)";
+                }
+                os << std::endl;
+		
                 t = update_time + (failureReport ? failureReport->update_time : 0.0);
                 os << " Update time (seconds):       " << t;
                 if (failureReport) {

opm/core/simulator/SimulatorReport.hpp

@@ -33,6 +33,7 @@ namespace Opm
         double total_time;
         double solver_time;
         double assemble_time;
+        double linear_solve_setup_time;
         double linear_solve_time;
         double update_time;
         double output_write_time;