From c94eec872fd2f3eaab79f92835a3e2e2faedfef8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 19 Jun 2020 16:45:34 +0200 Subject: [PATCH] Allow well operators with FlexibleSolver. --- opm/simulators/linalg/FlexibleSolver.hpp | 34 ++-- opm/simulators/linalg/FlexibleSolver_impl.hpp | 103 +++++----- opm/simulators/linalg/ISTLSolverEbos.hpp | 34 ++-- .../linalg/ISTLSolverEbosFlexible.hpp | 178 +++++++++++------- .../linalg/PressureSolverPolicy.hpp | 18 +- tests/test_flexiblesolver.cpp | 4 +- 6 files changed, 225 insertions(+), 146 deletions(-) diff --git a/opm/simulators/linalg/FlexibleSolver.hpp b/opm/simulators/linalg/FlexibleSolver.hpp index 826611302..d6383d414 100644 --- a/opm/simulators/linalg/FlexibleSolver.hpp +++ b/opm/simulators/linalg/FlexibleSolver.hpp @@ -43,54 +43,56 @@ public: using MatrixType = MatrixTypeT; using VectorType = VectorTypeT; + /// Base class type of the operator passed to the solver. + using AbstractOperatorType = Dune::AssembledLinearOperator; + /// Base class type of the contained preconditioner. + using AbstractPrecondType = Dune::PreconditionerWithUpdate; + /// Create a sequential solver. - FlexibleSolver(const MatrixType& matrix, + FlexibleSolver(AbstractOperatorType& op, const boost::property_tree::ptree& prm, - const std::function& weightsCalculator = std::function()); + const std::function& weightsCalculator = std::function()); /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication). template - FlexibleSolver(const MatrixType& matrix, + FlexibleSolver(AbstractOperatorType& op, const Comm& comm, const boost::property_tree::ptree& prm, - const std::function& weightsCalculator = std::function()); + const std::function& weightsCalculator = std::function()); virtual void apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res) override; virtual void apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res) override; - /// Type of the contained preconditioner. - using AbstractPrecondType = Dune::PreconditionerWithUpdate; - /// Access the contained preconditioner. AbstractPrecondType& preconditioner(); virtual Dune::SolverCategory::Category category() const override; private: - using AbstractOperatorType = Dune::AssembledLinearOperator; using AbstractScalarProductType = Dune::ScalarProduct; using AbstractSolverType = Dune::InverseOperator; // Machinery for making sequential or parallel operators/preconditioners/scalar products. template - void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, - const std::function weightsCalculator, const Comm& comm); + void initOpPrecSp(AbstractOperatorType& op, const boost::property_tree::ptree& prm, + const std::function weightsCalculator, const Comm& comm); - void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, - const std::function weightsCalculator, const Dune::Amg::SequentialInformation&); + void initOpPrecSp(AbstractOperatorType& op, const boost::property_tree::ptree& prm, + const std::function weightsCalculator, const Dune::Amg::SequentialInformation&); - void initSolver(const boost::property_tree::ptree& prm, bool isMaster); + void initSolver(const boost::property_tree::ptree& prm, const bool is_iorank); // Main initialization routine. // Call with Comm == Dune::Amg::SequentialInformation to get a serial solver. template - void init(const MatrixType& matrix, + void init(AbstractOperatorType& op, const Comm& comm, const boost::property_tree::ptree& prm, - const std::function weightsCalculator); + const std::function weightsCalculator); - std::shared_ptr linearoperator_; + AbstractOperatorType* linearoperator_for_solver_; + std::shared_ptr linearoperator_for_precond_; std::shared_ptr preconditioner_; std::shared_ptr scalarproduct_; std::shared_ptr linsolver_; diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp index 8394bd055..054830540 100644 --- a/opm/simulators/linalg/FlexibleSolver_impl.hpp +++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp @@ -39,23 +39,23 @@ namespace Dune /// Create a sequential solver. template FlexibleSolver:: - FlexibleSolver(const MatrixType& matrix, + FlexibleSolver(AbstractOperatorType& op, const boost::property_tree::ptree& prm, const std::function& weightsCalculator) { - init(matrix, Dune::Amg::SequentialInformation(), prm, weightsCalculator); + init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator); } /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication). template template FlexibleSolver:: - FlexibleSolver(const MatrixType& matrix, + FlexibleSolver(AbstractOperatorType& op, const Comm& comm, const boost::property_tree::ptree& prm, const std::function& weightsCalculator) { - init(matrix, comm, prm, weightsCalculator); + init(op, comm, prm, weightsCalculator); } template @@ -88,7 +88,7 @@ namespace Dune FlexibleSolver:: category() const { - return linearoperator_->category(); + return linearoperator_for_solver_->category(); } // Machinery for making sequential or parallel operators/preconditioners/scalar products. @@ -96,56 +96,64 @@ namespace Dune template void FlexibleSolver:: - initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, - const std::function weightsCalculator, const Comm& comm) + initOpPrecSp(AbstractOperatorType& op, + const boost::property_tree::ptree& prm, + const std::function weightsCalculator, + const Comm& comm) { // Parallel case. - using ParOperatorType = Dune::OverlappingSchwarzOperator; using pt = const boost::property_tree::ptree; - auto linop = std::make_shared(matrix, comm); - linearoperator_ = linop; + using ParOperatorType = Dune::OverlappingSchwarzOperator; + linearoperator_for_solver_ = &op; + auto op_prec = std::make_shared(op.getmat(), comm); auto child = prm.get_child_optional("preconditioner"); - preconditioner_ - = Opm::PreconditionerFactory::create(*linop, child? *child : pt(), - weightsCalculator, comm); - scalarproduct_ = Dune::createScalarProduct(comm, linearoperator_->category()); + preconditioner_ = Opm::PreconditionerFactory::create(*op_prec, + child ? *child : pt(), + weightsCalculator, + comm); + scalarproduct_ = Dune::createScalarProduct(comm, op.category()); + linearoperator_for_precond_ = op_prec; } template void FlexibleSolver:: - initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, - const std::function weightsCalculator, const Dune::Amg::SequentialInformation&) + initOpPrecSp(AbstractOperatorType& op, + const boost::property_tree::ptree& prm, + const std::function weightsCalculator, + const Dune::Amg::SequentialInformation&) { // Sequential case. - using SeqOperatorType = Dune::MatrixAdapter; using pt = const boost::property_tree::ptree; - auto linop = std::make_shared(matrix); - linearoperator_ = linop; + using SeqOperatorType = Dune::MatrixAdapter; + linearoperator_for_solver_ = &op; + auto op_prec = std::make_shared(op.getmat()); auto child = prm.get_child_optional("preconditioner"); - preconditioner_ = Opm::PreconditionerFactory::create(*linop, child? *child : pt(), + preconditioner_ = Opm::PreconditionerFactory::create(*op_prec, + child ? *child : pt(), weightsCalculator); scalarproduct_ = std::make_shared>(); + linearoperator_for_precond_ = op_prec; } template void FlexibleSolver:: - initSolver(const boost::property_tree::ptree& prm, bool isMaster) + initSolver(const boost::property_tree::ptree& prm, const bool is_iorank) { const double tol = prm.get("tol", 1e-2); const int maxiter = prm.get("maxiter", 200); - const int verbosity = isMaster? prm.get("verbosity", 0) : 0; + const int verbosity = is_iorank ? prm.get("verbosity", 0) : 0; const std::string solver_type = prm.get("solver", "bicgstab"); if (solver_type == "bicgstab") { - linsolver_.reset(new Dune::BiCGSTABSolver(*linearoperator_, + linsolver_.reset(new Dune::BiCGSTABSolver(*linearoperator_for_solver_, *scalarproduct_, *preconditioner_, tol, // desired residual reduction factor maxiter, // maximum number of iterations verbosity)); } else if (solver_type == "loopsolver") { - linsolver_.reset(new Dune::LoopSolver(*linearoperator_, + linsolver_.reset(new Dune::LoopSolver(*linearoperator_for_solver_, *scalarproduct_, *preconditioner_, tol, // desired residual reduction factor @@ -153,7 +161,7 @@ namespace Dune verbosity)); } else if (solver_type == "gmres") { int restart = prm.get("restart", 15); - linsolver_.reset(new Dune::RestartedGMResSolver(*linearoperator_, + linsolver_.reset(new Dune::RestartedGMResSolver(*linearoperator_for_solver_, *scalarproduct_, *preconditioner_, tol, @@ -163,7 +171,7 @@ namespace Dune #if HAVE_SUITESPARSE_UMFPACK } else if (solver_type == "umfpack") { bool dummy = false; - linsolver_.reset(new Dune::UMFPack(linearoperator_->getmat(), verbosity, dummy)); + linsolver_.reset(new Dune::UMFPack(linearoperator_for_solver_->getmat(), verbosity, dummy)); #endif } else { OPM_THROW(std::invalid_argument, "Properties: Solver " << solver_type << " not known."); @@ -177,13 +185,13 @@ namespace Dune template void FlexibleSolver:: - init(const MatrixType& matrix, + init(AbstractOperatorType& op, const Comm& comm, const boost::property_tree::ptree& prm, const std::function weightsCalculator) { - initOpPrecSp(matrix, prm, weightsCalculator, comm); - initSolver(prm, comm.communicator().rank()==0); + initOpPrecSp(op, prm, weightsCalculator, comm); + initSolver(prm, comm.communicator().rank() == 0); } } // namespace Dune @@ -198,28 +206,33 @@ using BM = Dune::BCRSMatrix>; template using OBM = Dune::BCRSMatrix>; -// INSTANTIATE_CONSTRUCTOR instantiates the constructor that is a template, -// this is only needed in the MPI case, since otherwise the Comm type is -// not a template argument but always SequentialInformation. #if HAVE_MPI + using Comm = Dune::OwnerOverlapCopyCommunication; -#define INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n) \ -template Dune::FlexibleSolver, BV>::FlexibleSolver(const MatrixType& matrix, \ + +// Note: we must instantiate the constructor that is a template. +// This is only needed in the parallel case, since otherwise the Comm type is +// not a template argument but always SequentialInformation. + +#define INSTANTIATE_FLEXIBLESOLVER(N) \ +template class Dune::FlexibleSolver, BV>; \ +template class Dune::FlexibleSolver, BV>; \ +template Dune::FlexibleSolver, BV>::FlexibleSolver(AbstractOperatorType& op, \ + const Comm& comm, \ + const boost::property_tree::ptree& prm, \ + const std::function()>& weightsCalculator); \ +template Dune::FlexibleSolver, BV>::FlexibleSolver(AbstractOperatorType& op, \ const Comm& comm, \ const boost::property_tree::ptree& prm, \ - const std::function()>& weightsCalculator); -#else -#define INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n) -#endif + const std::function()>& weightsCalculator); -// INSTANTIATE instantiates the class including any templated constructors if necessary. -#define INSTANTIATE_FLEXIBLESOLVER(n) \ -/* Variants using Dune::FieldMatrix blocks. */ \ -template class Dune::FlexibleSolver, BV>; \ -/* Variants using Opm::MatrixBlock blocks. */ \ -template class Dune::FlexibleSolver, BV>; \ -INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n) +#else // HAVE_MPI +#define INSTANTIATE_FLEXIBLESOLVER(N) \ +template class Dune::FlexibleSolver, BV>; \ +template class Dune::FlexibleSolver, BV>; + +#endif // HAVE_MPI #endif // OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 75d9ea855..615654ad2 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -113,6 +113,8 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) typedef typename GridView::template Codim<0>::Entity Element; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; using FlexibleSolverType = Dune::FlexibleSolver; + using AbstractOperatorType = Dune::AssembledLinearOperator; + using WellModelOperator = WellModelAsLinearOperator; // Due to miscibility oil <-> gas the water eqn is the one we can replace with a pressure equation. static const bool waterEnabled = Indices::waterEnabled; static const int pindex = (waterEnabled) ? BlackOilDefaultIndexTraits::waterCompIdx : BlackOilDefaultIndexTraits::oilCompIdx; @@ -176,17 +178,10 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) #endif extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); - - if (!useWellConn_ && useFlexible_) - { - OPM_THROW(std::logic_error, "Flexible solvers and CPR need the well contribution in the matrix. Please run with" - " --matrix-add-well-contributions=true"); - } - ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst); interiorCellNum_ = detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_); - if ( isParallel() && (!ownersFirst_ || parameters_.linear_solver_use_amg_ || useFlexible_ ) ) { + if ( isParallel() && (!ownersFirst_ || parameters_.linear_solver_use_amg_) ) { detail::setWellConnections(gridForConn, simulator_.vanguard().schedule().getWellsatEnd(), useWellConn_, wellConnectionsGraph_); // For some reason simulator_.model().elementMapper() is not initialized at this stage // Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work. @@ -741,11 +736,26 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) if (recreate_solver || !flexibleSolver_) { if (isParallel()) { #if HAVE_MPI - assert(noGhostMat_); - flexibleSolver_.reset(new FlexibleSolverType(getMatrix(), *comm_, prm_, weightsCalculator)); + if (useWellConn_) { + assert(noGhostMat_); + using ParOperatorType = Dune::OverlappingSchwarzOperator; + linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *comm_); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator); + } else { + if (!ownersFirst_) { + OPM_THROW(std::runtime_error, "In parallel, the flexible solver requires " + "--owner-cells-first=true when --matrix-add-well-contributions=false is used."); + } + using ParOperatorType = WellModelGhostLastMatrixAdapter; + wellOperator_ = std::make_unique(simulator_.problem().wellModel()); + linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *wellOperator_, interiorCellNum_); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator); + } #endif } else { - flexibleSolver_.reset(new FlexibleSolverType(getMatrix(), prm_, weightsCalculator)); + using SeqLinearOperator = Dune::MatrixAdapter; + linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix()); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator); } } else @@ -1004,6 +1014,8 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) Vector *rhs_; std::unique_ptr flexibleSolver_; + std::unique_ptr linearOperatorForFlexibleSolver_; + std::unique_ptr> wellOperator_; std::vector overlapRows_; std::vector interiorRows_; std::vector> wellConnectionsGraph_; diff --git a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp index 5403a9744..6361db03d 100644 --- a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp +++ b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp @@ -61,18 +61,19 @@ class ISTLSolverEbosFlexible using Simulator = typename GET_PROP_TYPE(TypeTag, Simulator); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using MatrixType = typename SparseMatrixAdapter::IstlMatrix; + using WellModel = typename GET_PROP_TYPE(TypeTag, EclWellModel); #if HAVE_MPI using Communication = Dune::OwnerOverlapCopyCommunication; #else using Communication = int; // Dummy type. #endif + using AbstractOperatorType = Dune::AssembledLinearOperator; + using WellModelOpType = WellModelAsLinearOperator; using SolverType = Dune::FlexibleSolver; // for quasiImpesWeights typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - //typedef typename GET_PROP_TYPE(TypeTag, EclWellModel) WellModel; - //typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename SparseMatrixAdapter::IstlMatrix Matrix; typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType; typedef typename Vector::block_type BlockVector; @@ -90,6 +91,9 @@ public: explicit ISTLSolverEbosFlexible(const Simulator& simulator) : simulator_(simulator) + , ownersFirst_(EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst)) + , matrixAddWellContributions_(EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions)) + , interiorCellNum_(detail::numMatrixRowsToUseInSolver(simulator_.vanguard().grid(), ownersFirst_)) { parameters_.template init(); prm_ = setupPropertyTree(parameters_); @@ -134,76 +138,53 @@ public: parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1); firstcall = false; } - makeOverlapRowsInvalid(mat.istlMatrix()); + if (isParallel() && matrixAddWellContributions_) { + makeOverlapRowsInvalid(mat.istlMatrix()); + } #endif - // Decide if we should recreate the solver or just do - // a minimal preconditioner update. - const int newton_iteration = this->simulator_.model().newtonMethod().numIterations(); - bool recreate_solver = false; - if (this->parameters_.cpr_reuse_setup_ == 0) { - // Always recreate solver. - recreate_solver = true; - } else if (this->parameters_.cpr_reuse_setup_ == 1) { - // Recreate solver on the first iteration of every timestep. - if (newton_iteration == 0) { - recreate_solver = true; - } - } else if (this->parameters_.cpr_reuse_setup_ == 2) { - // Recreate solver if the last solve used more than 10 iterations. - if (this->iterations() > 10) { - recreate_solver = true; - } - } else { - assert(this->parameters_.cpr_reuse_setup_ == 3); - assert(recreate_solver == false); - // Never recreate solver. - } - std::function weightsCalculator; + matrix_ = &mat.istlMatrix(); // Store pointer for output if needed. + std::function weightsCalculator = getWeightsCalculator(mat.istlMatrix(), b); - auto preconditionerType = prm_.get("preconditioner.type", "cpr"); - if( preconditionerType == "cpr" || - preconditionerType == "cprt" - ) - { - bool transpose = false; - if(preconditionerType == "cprt"){ - transpose = true; - } - - auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes"); - auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1); - if(weightsType == "quasiimpes") { - // weighs will be created as default in the solver - weightsCalculator = - [&mat, transpose, pressureIndex](){ - return Opm::Amg::getQuasiImpesWeights( - mat.istlMatrix(), - pressureIndex, - transpose); - }; - - }else if(weightsType == "trueimpes" ){ - weightsCalculator = - [this, &b, pressureIndex](){ - return this->getTrueImpesWeights(b, pressureIndex); - }; - }else{ - OPM_THROW(std::invalid_argument, "Weights type " << weightsType << "not implemented for cpr." - << " Please use quasiimpes or trueimpes."); - } - } - - if (recreate_solver || !solver_) { + if (shouldCreateSolver()) { if (isParallel()) { #if HAVE_MPI - matrix_ = &mat.istlMatrix(); - solver_.reset(new SolverType(mat.istlMatrix(), *comm_, prm_, weightsCalculator)); -#endif + if (matrixAddWellContributions_) { + using ParOperatorType = Dune::OverlappingSchwarzOperator; + auto op = std::make_unique(mat.istlMatrix(), *comm_); + auto sol = std::make_unique(*op, *comm_, prm_, weightsCalculator); + solver_ = std::move(sol); + linear_operator_ = std::move(op); + } else { + if (!ownersFirst_) { + OPM_THROW(std::runtime_error, "In parallel, the flexible solver requires " + "--owner-cells-first=true when --matrix-add-well-contributions=false is used."); + } + using ParOperatorType = WellModelGhostLastMatrixAdapter; + auto well_op = std::make_unique(simulator_.problem().wellModel()); + auto op = std::make_unique(mat.istlMatrix(), *well_op, interiorCellNum_); + auto sol = std::make_unique(*op, *comm_, prm_, weightsCalculator); + solver_ = std::move(sol); + linear_operator_ = std::move(op); + well_operator_ = std::move(well_op); + } +#endif // HAVE_MPI } else { - matrix_ = &mat.istlMatrix(); - solver_.reset(new SolverType(mat.istlMatrix(), prm_, weightsCalculator)); + if (matrixAddWellContributions_) { + using SeqOperatorType = Dune::MatrixAdapter; + auto op = std::make_unique(mat.istlMatrix()); + auto sol = std::make_unique(*op, prm_, weightsCalculator); + solver_ = std::move(sol); + linear_operator_ = std::move(op); + } else { + using SeqOperatorType = WellModelMatrixAdapter; + auto well_op = std::make_unique(simulator_.problem().wellModel()); + auto op = std::make_unique(mat.istlMatrix(), *well_op); + auto sol = std::make_unique(*op, prm_, weightsCalculator); + solver_ = std::move(sol); + linear_operator_ = std::move(op); + well_operator_ = std::move(well_op); + } } rhs_ = b; } else { @@ -245,6 +226,63 @@ public: protected: + bool shouldCreateSolver() const + { + // Decide if we should recreate the solver or just do + // a minimal preconditioner update. + if (!solver_) { + return true; + } + const int newton_iteration = this->simulator_.model().newtonMethod().numIterations(); + bool recreate_solver = false; + if (this->parameters_.cpr_reuse_setup_ == 0) { + // Always recreate solver. + recreate_solver = true; + } else if (this->parameters_.cpr_reuse_setup_ == 1) { + // Recreate solver on the first iteration of every timestep. + if (newton_iteration == 0) { + recreate_solver = true; + } + } else if (this->parameters_.cpr_reuse_setup_ == 2) { + // Recreate solver if the last solve used more than 10 iterations. + if (this->iterations() > 10) { + recreate_solver = true; + } + } else { + assert(this->parameters_.cpr_reuse_setup_ == 3); + assert(recreate_solver == false); + // Never recreate solver. + } + return recreate_solver; + } + + std::function getWeightsCalculator(const MatrixType& mat, const VectorType& b) const + { + std::function weightsCalculator; + + auto preconditionerType = prm_.get("preconditioner.type", "cpr"); + if (preconditionerType == "cpr" || preconditionerType == "cprt") { + const bool transpose = preconditionerType == "cprt"; + const auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes"); + const auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1); + if (weightsType == "quasiimpes") { + // weighs will be created as default in the solver + weightsCalculator = [&mat, transpose, pressureIndex]() { + return Opm::Amg::getQuasiImpesWeights(mat, pressureIndex, transpose); + }; + } else if (weightsType == "trueimpes") { + weightsCalculator = [this, &b, pressureIndex]() { + return this->getTrueImpesWeights(b, pressureIndex); + }; + } else { + OPM_THROW(std::invalid_argument, + "Weights type " << weightsType << "not implemented for cpr." + << " Please use quasiimpes or trueimpes."); + } + } + return weightsCalculator; + } + /// Zero out off-diagonal blocks on rows corresponding to overlap cells /// Diagonal blocks on ovelap rows are set to diag(1.0). void makeOverlapRowsInvalid(MatrixType& matrix) const @@ -267,7 +305,7 @@ protected: } } - VectorType getTrueImpesWeights(const VectorType& b,const int pressureVarIndex) + VectorType getTrueImpesWeights(const VectorType& b, const int pressureVarIndex) const { VectorType weights(b.size()); ElementContext elemCtx(simulator_); @@ -289,14 +327,20 @@ protected: } } + const Simulator& simulator_; MatrixType* matrix_; + std::unique_ptr well_operator_; + std::unique_ptr linear_operator_; std::unique_ptr solver_; FlowLinearSolverParameters parameters_; boost::property_tree::ptree prm_; VectorType rhs_; Dune::InverseOperatorResult res_; std::any parallelInformation_; + bool ownersFirst_; + bool matrixAddWellContributions_; + int interiorCellNum_; std::unique_ptr comm_; std::vector overlapRows_; std::vector interiorRows_; diff --git a/opm/simulators/linalg/PressureSolverPolicy.hpp b/opm/simulators/linalg/PressureSolverPolicy.hpp index 2831b69b8..699ac1611 100644 --- a/opm/simulators/linalg/PressureSolverPolicy.hpp +++ b/opm/simulators/linalg/PressureSolverPolicy.hpp @@ -25,6 +25,7 @@ #include #include +#include namespace Dune { @@ -55,19 +56,24 @@ namespace Amg * The operator will use one step of AMG to approximately solve * the coarse level system. */ - struct PressureInverseOperator : public Dune::InverseOperator { - template - PressureInverseOperator(Operator& op, const boost::property_tree::ptree& prm, const Comm& comm) + struct PressureInverseOperator : public Dune::InverseOperator + { + template + PressureInverseOperator(Operator& op, + const boost::property_tree::ptree& prm, + const Dune::OwnerOverlapCopyCommunication& comm) : linsolver_() { assert(op.category() == Dune::SolverCategory::overlapping); - linsolver_ = std::make_unique(op.getmat(), comm, prm, std::function()); + linsolver_ = std::make_unique(op, comm, prm, std::function()); } - PressureInverseOperator(Operator& op, const boost::property_tree::ptree& prm, const SequentialInformation&) + PressureInverseOperator(Operator& op, + const boost::property_tree::ptree& prm, + const SequentialInformation&) : linsolver_() { assert(op.category() != Dune::SolverCategory::overlapping); - linsolver_ = std::make_unique(op.getmat(), prm, std::function()); + linsolver_ = std::make_unique(op, prm, std::function()); } diff --git a/tests/test_flexiblesolver.cpp b/tests/test_flexiblesolver.cpp index 215795daf..f02e14f4e 100644 --- a/tests/test_flexiblesolver.cpp +++ b/tests/test_flexiblesolver.cpp @@ -75,7 +75,9 @@ testSolver(const boost::property_tree::ptree& prm, const std::string& matrix_fil prm.get("preconditioner.pressure_var_index"), transpose); }; - Dune::FlexibleSolver solver(matrix, prm, wc); + using SeqOperatorType = Dune::MatrixAdapter; + SeqOperatorType op(matrix); + Dune::FlexibleSolver solver(op, prm, wc); Vector x(rhs.size()); Dune::InverseOperatorResult res; solver.apply(x, rhs, res);