diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 354b7e0ad..fc3fb8f2a 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -146,7 +146,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/linalg/bda/cusparseSolverBackend.hpp
opm/simulators/linalg/bda/WellContributions.hpp
opm/simulators/linalg/BlackoilAmg.hpp
- opm/simulators/linalg/BlackoilAmgCpr.hpp
opm/simulators/linalg/amgcpr.hh
opm/simulators/linalg/twolevelmethodcpr.hh
opm/simulators/linalg/CPRPreconditioner.hpp
@@ -155,7 +154,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/linalg/FlowLinearSolverParameters.hpp
opm/simulators/linalg/GraphColoring.hpp
opm/simulators/linalg/ISTLSolverEbos.hpp
- opm/simulators/linalg/ISTLSolverEbosCpr.hpp
opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
opm/simulators/linalg/MatrixBlock.hpp
opm/simulators/linalg/OwningBlockPreconditioner.hpp
diff --git a/opm/simulators/linalg/BlackoilAmgCpr.hpp b/opm/simulators/linalg/BlackoilAmgCpr.hpp
deleted file mode 100644
index a8fbe6ac0..000000000
--- a/opm/simulators/linalg/BlackoilAmgCpr.hpp
+++ /dev/null
@@ -1,173 +0,0 @@
-/*
- 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 .
-*/
-#ifndef OPM_AMGCPR_HEADER_INCLUDED
-#define OPM_AMGCPR_HEADER_INCLUDED
-
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-
-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
- class BlackoilAmgCpr
- : public Dune::Preconditioner
- {
- 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::Arguments;
-
- protected:
- using Matrix = typename Operator::matrix_type;
- using CoarseOperator = typename Detail::ScalarType::value;
- using CoarseSmoother = typename Detail::ScalarType::value;
- using FineCriterion =
- typename Detail::OneComponentCriterionType::value;
- using CoarseCriterion = typename Detail::ScalarType::value;
- using LevelTransferPolicy =
- OneComponentAggregationLevelTransferPolicy;
- using CoarseSolverPolicy =
- Detail::OneStepAMGCoarseSolverPolicy;
- using TwoLevelMethod =
- Dune::Amg::TwoLevelMethodCpr;
- public:
- Dune::SolverCategory::Category category() const override
- {
- return std::is_same::value ?
- Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
- }
-
- /**
- * \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(*scaledMatrixOperator_,
- smargs, comm)),
- levelTransferPolicy_(criterion, comm, param.cpr_pressure_aggregation_),
- coarseSolverPolicy_(¶m, 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_ = Detail::constructSmoother(*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 scaledMatrix_;
- std::unique_ptr scaledMatrixOperator_;
- std::shared_ptr smoother_;
- LevelTransferPolicy levelTransferPolicy_;
- CoarseSolverPolicy coarseSolverPolicy_;
- TwoLevelMethod twoLevelMethod_;
- };
-
-} // end namespace Opm
-#endif
diff --git a/opm/simulators/linalg/ISTLSolverEbosCpr.hpp b/opm/simulators/linalg/ISTLSolverEbosCpr.hpp
deleted file mode 100644
index 7fb7e2c9f..000000000
--- a/opm/simulators/linalg/ISTLSolverEbosCpr.hpp
+++ /dev/null
@@ -1,303 +0,0 @@
-/*
- 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 .
-*/
-
-#ifndef OPM_ISTLSOLVERCPR_EBOS_HEADER_INCLUDED
-#define OPM_ISTLSOLVERCPR_EBOS_HEADER_INCLUDED
-
-#include
-#include
-#include
-#include
-
-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 ISTLSolverEbosCpr : public ISTLSolverEbos
- {
- protected:
- // Types and indices from superclass.
- using SuperClass = ISTLSolverEbos;
- 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;
- using MatrixAdapter = Dune::MatrixAdapter;
-
- using CouplingMetric = Opm::Amg::Element;
- using CritBase = Dune::Amg::SymmetricCriterion;
- using Criterion = Dune::Amg::CoarsenCriterion;
-
- using CprSmootherFine = Opm::ParallelOverlappingILU0;
- using CprSmootherCoarse = CprSmootherFine;
- using BlackoilAmgType = BlackoilAmgCpr;
- using OperatorSerial = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false>;
-
-#if HAVE_MPI
- using POrComm = Dune::OwnerOverlapCopyCommunication;
- using ParallelMatrixAdapter = Dune::OverlappingSchwarzOperator;
- using ParallelCprSmootherFine = Opm::ParallelOverlappingILU0;
- using ParallelCprSmootherCoarse = ParallelCprSmootherFine;
- using ParallelBlackoilAmgType = BlackoilAmgCpr;
- using OperatorParallel = WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true>;
- using ParallelScalarProduct = Dune::OverlappingSchwarzScalarProduct;
-#else
- using POrComm = Dune::Amg::SequentialInformation;
- using ParallelBlackoilAmgType = BlackoilAmgType;
- using ParallelScalarProduct = Dune::SeqScalarProduct;
- using ParallelMatrixAdapter = MatrixAdapter;
- using OperatorParallel = OperatorSerial;
-#endif
-
- public:
- static void registerParameters()
- {
- FlowLinearSolverParameters::registerParameters();
- }
-
- /// 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()
- {}
-
- void prepare(const SparseMatrixAdapter& M, Vector& b)
- {
- if (oldMat != nullptr)
- std::cout << "old was "<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 =
- std::any_cast( 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_ ));
- }
- }
-
- constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::overlapping;
- auto sp = Dune::createScalarProduct(*comm_, category);
- sp_ = std::move(sp);
-
- using AMGOperator = Dune::OverlappingSchwarzOperator;
- // 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 POrCommType = Dune::Amg::SequentialInformation;
- POrCommType parallelInformation_arg;
- typedef OperatorSerial LinearOperator;
-
- constexpr Dune::SolverCategory::Category category=Dune::SolverCategory::sequential;
- auto sp = Dune::createScalarProduct(parallelInformation_arg, category);
- sp_ = std::move(sp);
-
- // 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
- 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_;
-
- // 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::value,
- BlackoilAmgType, ParallelBlackoilAmgType>::type;
- using SpType = typename std::conditional::value,
- Dune::SeqScalarProduct,
- ParallelScalarProduct >::type;
- using OperatorType = typename std::conditional::value,
- MatrixAdapter, ParallelMatrixAdapter>::type;
- typedef typename AmgType::Smoother Smoother;
- typedef typename Dune::Amg::SmootherTraits::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(*opA_);
- int newton_iteration = this->simulator_.model().newtonMethod().numIterations();
- 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(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(wellOpA, reinterpret_cast(*sp_), reinterpret_cast(*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 > 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 >;
- SPPointer sp_;
- std::shared_ptr< Dune::BiCGSTABSolver > linsolve_;
- const void* oldMat;
- std::shared_ptr comm_;
- }; // end ISTLSolver
-
-} // namespace Opm
-#endif