diff --git a/opm/simulators/linalg/FlexibleSolver.hpp b/opm/simulators/linalg/FlexibleSolver.hpp index 352506a92..925e01412 100644 --- a/opm/simulators/linalg/FlexibleSolver.hpp +++ b/opm/simulators/linalg/FlexibleSolver.hpp @@ -50,14 +50,16 @@ public: /// Create a sequential solver. FlexibleSolver(AbstractOperatorType& op, const Opm::PropertyTree& prm, - const std::function& weightsCalculator = std::function()); + const std::function& weightsCalculator, + std::size_t pressureIndex); /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication). template FlexibleSolver(AbstractOperatorType& op, const Comm& comm, const Opm::PropertyTree& prm, - const std::function& weightsCalculator = std::function()); + const std::function& weightsCalculator, + std::size_t pressureIndex); virtual void apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res) override; @@ -75,10 +77,12 @@ private: // Machinery for making sequential or parallel operators/preconditioners/scalar products. template void initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm, - const std::function weightsCalculator, const Comm& comm); + const std::function weightsCalculator, const Comm& comm, + std::size_t pressureIndex); void initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm, - const std::function weightsCalculator, const Dune::Amg::SequentialInformation&); + const std::function weightsCalculator, const Dune::Amg::SequentialInformation&, + std::size_t pressureIndex); void initSolver(const Opm::PropertyTree& prm, const bool is_iorank); @@ -88,7 +92,8 @@ private: void init(AbstractOperatorType& op, const Comm& comm, const Opm::PropertyTree& prm, - const std::function weightsCalculator); + const std::function weightsCalculator, + std::size_t pressureIndex); AbstractOperatorType* linearoperator_for_solver_; std::shared_ptr linearoperator_for_precond_; diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp index 69fdb1e35..e1674f834 100644 --- a/opm/simulators/linalg/FlexibleSolver_impl.hpp +++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp @@ -39,9 +39,11 @@ namespace Dune FlexibleSolver:: FlexibleSolver(AbstractOperatorType& op, const Opm::PropertyTree& prm, - const std::function& weightsCalculator) + const std::function& weightsCalculator, + std::size_t pressureIndex) { - init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator); + init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator, + pressureIndex); } /// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication). @@ -51,9 +53,10 @@ namespace Dune FlexibleSolver(AbstractOperatorType& op, const Comm& comm, const Opm::PropertyTree& prm, - const std::function& weightsCalculator) + const std::function& weightsCalculator, + std::size_t pressureIndex) { - init(op, comm, prm, weightsCalculator); + init(op, comm, prm, weightsCalculator, pressureIndex); } template @@ -97,7 +100,8 @@ namespace Dune initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm, const std::function weightsCalculator, - const Comm& comm) + const Comm& comm, + std::size_t pressureIndex) { // Parallel case. using ParOperatorType = Dune::OverlappingSchwarzOperator; @@ -107,7 +111,8 @@ namespace Dune preconditioner_ = Opm::PreconditionerFactory::create(*op_prec, child ? *child : Opm::PropertyTree(), weightsCalculator, - comm); + comm, + pressureIndex); scalarproduct_ = Dune::createScalarProduct(comm, op.category()); linearoperator_for_precond_ = op_prec; } @@ -118,7 +123,8 @@ namespace Dune initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm, const std::function weightsCalculator, - const Dune::Amg::SequentialInformation&) + const Dune::Amg::SequentialInformation&, + std::size_t pressureIndex) { // Sequential case. using SeqOperatorType = Dune::MatrixAdapter; @@ -127,7 +133,8 @@ namespace Dune auto child = prm.get_child_optional("preconditioner"); preconditioner_ = Opm::PreconditionerFactory::create(*op_prec, child ? *child : Opm::PropertyTree(), - weightsCalculator); + weightsCalculator, + pressureIndex); scalarproduct_ = std::make_shared>(); linearoperator_for_precond_ = op_prec; } @@ -184,9 +191,10 @@ namespace Dune init(AbstractOperatorType& op, const Comm& comm, const Opm::PropertyTree& prm, - const std::function weightsCalculator) + const std::function weightsCalculator, + std::size_t pressureIndex) { - initOpPrecSp(op, prm, weightsCalculator, comm); + initOpPrecSp(op, prm, weightsCalculator, comm, pressureIndex); initSolver(prm, comm.communicator().rank() == 0); } @@ -216,11 +224,13 @@ template class Dune::FlexibleSolver, BV>; \ template Dune::FlexibleSolver, BV>::FlexibleSolver(AbstractOperatorType& op, \ const Comm& comm, \ const Opm::PropertyTree& prm, \ - const std::function()>& weightsCalculator); \ + const std::function()>& weightsCalculator, \ +std::size_t); \ template Dune::FlexibleSolver, BV>::FlexibleSolver(AbstractOperatorType& op, \ const Comm& comm, \ const Opm::PropertyTree& prm, \ - const std::function()>& weightsCalculator); + const std::function()>& weightsCalculator, \ +std::size_t); #else // HAVE_MPI diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 43119c9c2..53dcd8957 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -91,6 +91,7 @@ namespace Opm using AbstractOperatorType = Dune::AssembledLinearOperator; using WellModelOperator = WellModelAsLinearOperator; using ElementMapper = GetPropType; + constexpr static std::size_t pressureIndex = GetPropType::pressureSwitchIdx; #if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA static const unsigned int block_size = Matrix::block_type::rows; @@ -360,24 +361,28 @@ namespace Opm if (useWellConn_) { using ParOperatorType = Dune::OverlappingSchwarzOperator; linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *comm_); - flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator, + pressureIndex); } else { 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); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator, + pressureIndex); } #endif } else { if (useWellConn_) { using SeqOperatorType = Dune::MatrixAdapter; linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix()); - flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator, + pressureIndex); } else { using SeqOperatorType = WellModelMatrixAdapter; wellOperator_ = std::make_unique(simulator_.problem().wellModel()); linearOperatorForFlexibleSolver_ = std::make_unique(getMatrix(), *wellOperator_); - flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator); + flexibleSolver_ = std::make_unique(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator, + pressureIndex); } } } @@ -429,15 +434,18 @@ namespace Opm if (preconditionerType == "cpr" || preconditionerType == "cprt") { const bool transpose = preconditionerType == "cprt"; const auto weightsType = prm_.get("preconditioner.weight_type"s, "quasiimpes"s); - const auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1); if (weightsType == "quasiimpes") { - // weighs will be created as default in the solver - weightsCalculator = [this, transpose, pressureIndex]() { - return Amg::getQuasiImpesWeights(this->getMatrix(), pressureIndex, transpose); + // weights will be created as default in the solver + // assignment p = pressureIndex prevent compiler warning about + // capturing variable with non-automatic storage duration + weightsCalculator = [this, transpose, p = pressureIndex]() { + return Amg::getQuasiImpesWeights(this->getMatrix(), p, transpose); }; } else if (weightsType == "trueimpes") { - weightsCalculator = [this, pressureIndex]() { - return this->getTrueImpesWeights(pressureIndex); + // assignment p = pressureIndex prevent compiler warning about + // capturing variable with non-automatic storage duration + weightsCalculator = [this, p = pressureIndex]() { + return this->getTrueImpesWeights(p); }; } else { OPM_THROW(std::invalid_argument, diff --git a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp index 9ca7221b4..67b3e2f9c 100644 --- a/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp +++ b/opm/simulators/linalg/ISTLSolverEbosFlexible.hpp @@ -83,6 +83,7 @@ class ISTLSolverEbosFlexible using ThreadManager = GetPropType; typedef typename GridView::template Codim<0>::Entity Element; using ElementContext = GetPropType; + constexpr static std::size_t pressureIndex = GetPropType::pressureSwitchIdx; public: static void registerParameters() @@ -156,7 +157,8 @@ public: if (matrixAddWellContributions_) { using ParOperatorType = Dune::OverlappingSchwarzOperator; auto op = std::make_unique(mat.istlMatrix(), *comm_); - auto sol = std::make_unique(*op, *comm_, prm_, weightsCalculator); + auto sol = std::make_unique(*op, *comm_, prm_, weightsCalculator, + pressureIndex); solver_ = std::move(sol); linear_operator_ = std::move(op); } else { @@ -167,7 +169,8 @@ public: 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); + auto sol = std::make_unique(*op, *comm_, prm_, weightsCalculator, + pressureIndex); solver_ = std::move(sol); linear_operator_ = std::move(op); well_operator_ = std::move(well_op); @@ -177,14 +180,16 @@ public: if (matrixAddWellContributions_) { using SeqOperatorType = Dune::MatrixAdapter; auto op = std::make_unique(mat.istlMatrix()); - auto sol = std::make_unique(*op, prm_, weightsCalculator); + auto sol = std::make_unique(*op, prm_, weightsCalculator, + pressureIndex); 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); + auto sol = std::make_unique(*op, prm_, weightsCalculator, + pressureIndex); solver_ = std::move(sol); linear_operator_ = std::move(op); well_operator_ = std::move(well_op); @@ -270,15 +275,14 @@ protected: if (preconditionerType == "cpr" || preconditionerType == "cprt") { const bool transpose = preconditionerType == "cprt"; const auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes"s); - 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); + weightsCalculator = [&mat, transpose, p = pressureIndex]() { + return Opm::Amg::getQuasiImpesWeights(mat, p, transpose); }; } else if (weightsType == "trueimpes") { - weightsCalculator = [this, &b, pressureIndex]() { - return this->getTrueImpesWeights(b, pressureIndex); + weightsCalculator = [this, &b, p = pressureIndex]() { + return this->getTrueImpesWeights(b, p); }; } else { OPM_THROW(std::invalid_argument, diff --git a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp index fec4493be..dd9eddcf8 100644 --- a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp +++ b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp @@ -84,15 +84,17 @@ public: using PrecFactory = Opm::PreconditionerFactory; OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm, - const std::function weightsCalculator) + const std::function weightsCalculator, + std::size_t pressureIndex) : linear_operator_(linearoperator) , finesmoother_(PrecFactory::create(linearoperator, prm.get_child_optional("finesmoother") ? - prm.get_child("finesmoother") : Opm::PropertyTree())) + prm.get_child("finesmoother") : Opm::PropertyTree(), + std::function(), pressureIndex)) , comm_(nullptr) , weightsCalculator_(weightsCalculator) , weights_(weightsCalculator()) - , levelTransferPolicy_(dummy_comm_, weights_, prm.get("pressure_var_index")) + , levelTransferPolicy_(dummy_comm_, weights_, pressureIndex) , coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : Opm::PropertyTree()) , twolevel_method_(linearoperator, finesmoother_, @@ -113,15 +115,18 @@ public: } OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm, - const std::function weightsCalculator, const Communication& comm) + const std::function weightsCalculator, + std::size_t pressureIndex, const Communication& comm) : linear_operator_(linearoperator) , finesmoother_(PrecFactory::create(linearoperator, prm.get_child_optional("finesmoother") ? - prm.get_child("finesmoother"): Opm::PropertyTree(), comm)) + prm.get_child("finesmoother"): Opm::PropertyTree(), + std::function(), + comm, pressureIndex)) , comm_(&comm) , weightsCalculator_(weightsCalculator) , weights_(weightsCalculator()) - , levelTransferPolicy_(*comm_, weights_, prm.get("pressure_var_index", 1)) + , levelTransferPolicy_(*comm_, weights_, pressureIndex) , coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : Opm::PropertyTree()) , twolevel_method_(linearoperator, finesmoother_, diff --git a/opm/simulators/linalg/PreconditionerFactory.hpp b/opm/simulators/linalg/PreconditionerFactory.hpp index 9355cd293..9b2c35a91 100644 --- a/opm/simulators/linalg/PreconditionerFactory.hpp +++ b/opm/simulators/linalg/PreconditionerFactory.hpp @@ -36,6 +36,7 @@ #include #include +#include namespace Opm { @@ -57,8 +58,10 @@ public: using PrecPtr = std::shared_ptr>; /// The type of creator functions passed to addCreator(). - using Creator = std::function&)>; - using ParCreator = std::function&, const Comm&)>; + using Creator = std::function&, std::size_t)>; + using ParCreator = std::function&, std::size_t, const Comm&)>; /// Create a new serial preconditioner and return a pointer to it. /// \param op operator to be preconditioned. @@ -66,9 +69,10 @@ public: /// \param weightsCalculator Calculator for weights used in CPR. /// \return (smart) pointer to the created preconditioner. static PrecPtr create(const Operator& op, const PropertyTree& prm, - const std::function& weightsCalculator = std::function()) + const std::function& weightsCalculator = {}, + std::size_t pressureIndex = std::numeric_limits::max()) { - return instance().doCreate(op, prm, weightsCalculator); + return instance().doCreate(op, prm, weightsCalculator, pressureIndex); } /// Create a new parallel preconditioner and return a pointer to it. @@ -78,9 +82,10 @@ public: /// \param weightsCalculator Calculator for weights used in CPR. /// \return (smart) pointer to the created preconditioner. static PrecPtr create(const Operator& op, const PropertyTree& prm, - const std::function& weightsCalculator, const Comm& comm) + const std::function& weightsCalculator, const Comm& comm, + std::size_t pressureIndex = std::numeric_limits::max()) { - return instance().doCreate(op, prm, weightsCalculator, comm); + return instance().doCreate(op, prm, weightsCalculator, pressureIndex, comm); } /// Create a new parallel preconditioner and return a pointer to it. @@ -88,9 +93,10 @@ public: /// \param prm parameters for the preconditioner, in particular its type. /// \param comm communication object (typically OwnerOverlapCopyCommunication). /// \return (smart) pointer to the created preconditioner. - static PrecPtr create(const Operator& op, const PropertyTree& prm, const Comm& comm) + static PrecPtr create(const Operator& op, const PropertyTree& prm, const Comm& comm, + std::size_t pressureIndex = std::numeric_limits::max()) { - return instance().doCreate(op, prm, std::function(), comm); + return instance().doCreate(op, prm, std::function(), pressureIndex, comm); } /// Add a creator for a serial preconditioner to the PreconditionerFactory. /// After the call, the user may obtain a preconditioner by @@ -257,32 +263,32 @@ private: using V = Vector; using P = PropertyTree; using C = Comm; - doAddCreator("ILU0", [](const O& op, const P& prm, const std::function&, const C& comm) { + doAddCreator("ILU0", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { return createParILU(op, prm, comm, 0); }); - doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function&, const C& comm) { + doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { return createParILU(op, prm, comm, prm.get("ilulevel", 0)); }); - doAddCreator("ILUn", [](const O& op, const P& prm, const std::function&, const C& comm) { + doAddCreator("ILUn", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { return createParILU(op, prm, comm, prm.get("ilulevel", 0)); }); doAddCreator("Jac", [](const O& op, const P& prm, const std::function&, - const C& comm) { + std::size_t, const C& comm) { const int n = prm.get("repeats", 1); const double w = prm.get("relaxation", 1.0); return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); }); - doAddCreator("GS", [](const O& op, const P& prm, const std::function&, const C& comm) { + doAddCreator("GS", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { const int n = prm.get("repeats", 1); const double w = prm.get("relaxation", 1.0); return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); }); - doAddCreator("SOR", [](const O& op, const P& prm, const std::function&, const C& comm) { + doAddCreator("SOR", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { const int n = prm.get("repeats", 1); const double w = prm.get("relaxation", 1.0); return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); }); - doAddCreator("SSOR", [](const O& op, const P& prm, const std::function&, const C& comm) { + doAddCreator("SSOR", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { const int n = prm.get("repeats", 1); const double w = prm.get("relaxation", 1.0); return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); @@ -293,7 +299,7 @@ private: // later, but at this point no other operators are compatible // with the AMG hierarchy construction. if constexpr (std::is_same_v>) { - doAddCreator("amg", [](const O& op, const P& prm, const std::function&, const C& comm) { + doAddCreator("amg", [](const O& op, const P& prm, const std::function&, std::size_t, const C& comm) { const std::string smoother = prm.get("smoother", "ParOverILU0"); if (smoother == "ILU0" || smoother == "ParOverILU0") { using Smoother = Opm::ParallelOverlappingILU0; @@ -306,13 +312,21 @@ private: }); } - doAddCreator("cpr", [](const O& op, const P& prm, const std::function weightsCalculator, const C& comm) { + doAddCreator("cpr", [](const O& op, const P& prm, const std::function weightsCalculator, std::size_t pressureIndex, const C& comm) { assert(weightsCalculator); - return std::make_shared>(op, prm, weightsCalculator, comm); + if (pressureIndex == std::numeric_limits::max()) + { + OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); + } + return std::make_shared>(op, prm, weightsCalculator, pressureIndex, comm); }); - doAddCreator("cprt", [](const O& op, const P& prm, const std::function weightsCalculator, const C& comm) { + doAddCreator("cprt", [](const O& op, const P& prm, const std::function weightsCalculator, std::size_t pressureIndex, const C& comm) { assert(weightsCalculator); - return std::make_shared>(op, prm, weightsCalculator, comm); + if (pressureIndex == std::numeric_limits::max()) + { + OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); + } + return std::make_shared>(op, prm, weightsCalculator, pressureIndex, comm); }); } @@ -325,39 +339,39 @@ private: using M = Matrix; using V = Vector; using P = PropertyTree; - doAddCreator("ILU0", [](const O& op, const P& prm, const std::function&) { + doAddCreator("ILU0", [](const O& op, const P& prm, const std::function&, std::size_t) { const double w = prm.get("relaxation", 1.0); return std::make_shared>( op.getmat(), 0, w, Opm::MILU_VARIANT::ILU); }); - doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function&) { + doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function&, std::size_t) { const double w = prm.get("relaxation", 1.0); const int n = prm.get("ilulevel", 0); return std::make_shared>( op.getmat(), n, w, Opm::MILU_VARIANT::ILU); }); - doAddCreator("ILUn", [](const O& op, const P& prm, const std::function&) { + doAddCreator("ILUn", [](const O& op, const P& prm, const std::function&, std::size_t) { const int n = prm.get("ilulevel", 0); const double w = prm.get("relaxation", 1.0); return std::make_shared>( op.getmat(), n, w, Opm::MILU_VARIANT::ILU); }); - doAddCreator("Jac", [](const O& op, const P& prm, const std::function&) { + doAddCreator("Jac", [](const O& op, const P& prm, const std::function&, std::size_t) { const int n = prm.get("repeats", 1); const double w = prm.get("relaxation", 1.0); return wrapPreconditioner>(op.getmat(), n, w); }); - doAddCreator("GS", [](const O& op, const P& prm, const std::function&) { + doAddCreator("GS", [](const O& op, const P& prm, const std::function&, std::size_t) { const int n = prm.get("repeats", 1); const double w = prm.get("relaxation", 1.0); return wrapPreconditioner>(op.getmat(), n, w); }); - doAddCreator("SOR", [](const O& op, const P& prm, const std::function&) { + doAddCreator("SOR", [](const O& op, const P& prm, const std::function&, std::size_t) { const int n = prm.get("repeats", 1); const double w = prm.get("relaxation", 1.0); return wrapPreconditioner>(op.getmat(), n, w); }); - doAddCreator("SSOR", [](const O& op, const P& prm, const std::function&) { + doAddCreator("SSOR", [](const O& op, const P& prm, const std::function&, std::size_t) { const int n = prm.get("repeats", 1); const double w = prm.get("relaxation", 1.0); return wrapPreconditioner>(op.getmat(), n, w); @@ -366,7 +380,7 @@ private: // Only add AMG preconditioners to the factory if the operator // is an actual matrix operator. if constexpr (std::is_same_v>) { - doAddCreator("amg", [](const O& op, const P& prm, const std::function&) { + doAddCreator("amg", [](const O& op, const P& prm, const std::function&, std::size_t) { const std::string smoother = prm.get("smoother", "ParOverILU0"); if (smoother == "ILU0" || smoother == "ParOverILU0") { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) @@ -395,7 +409,7 @@ private: OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << "."); } }); - doAddCreator("kamg", [](const O& op, const P& prm, const std::function&) { + doAddCreator("kamg", [](const O& op, const P& prm, const std::function&, std::size_t) { const std::string smoother = prm.get("smoother", "ParOverILU0"); if (smoother == "ILU0" || smoother == "ParOverILU0") { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) @@ -427,7 +441,7 @@ private: OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << "."); } }); - doAddCreator("famg", [](const O& op, const P& prm, const std::function&) { + doAddCreator("famg", [](const O& op, const P& prm, const std::function&, std::size_t) { auto crit = amgCriterion(prm); Dune::Amg::Parameters parms; parms.setNoPreSmoothSteps(1); @@ -435,11 +449,19 @@ private: return wrapPreconditioner>(op, crit, parms); }); } - doAddCreator("cpr", [](const O& op, const P& prm, const std::function& weightsCalculator) { - return std::make_shared>(op, prm, weightsCalculator); + doAddCreator("cpr", [](const O& op, const P& prm, const std::function& weightsCalculator, std::size_t pressureIndex) { + if (pressureIndex == std::numeric_limits::max()) + { + OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); + } + return std::make_shared>(op, prm, weightsCalculator, pressureIndex); }); - doAddCreator("cprt", [](const O& op, const P& prm, const std::function& weightsCalculator) { - return std::make_shared>(op, prm, weightsCalculator); + doAddCreator("cprt", [](const O& op, const P& prm, const std::function& weightsCalculator, std::size_t pressureIndex) { + if (pressureIndex == std::numeric_limits::max()) + { + OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR"); + } + return std::make_shared>(op, prm, weightsCalculator, pressureIndex); }); } @@ -461,7 +483,8 @@ private: // Actually creates the product object. PrecPtr doCreate(const Operator& op, const PropertyTree& prm, - const std::function weightsCalculator) + const std::function weightsCalculator, + std::size_t pressureIndex) { const std::string& type = prm.get("type", "ParOverILU0"); auto it = creators_.find(type); @@ -474,11 +497,12 @@ private: msg << std::endl; OPM_THROW(std::invalid_argument, msg.str()); } - return it->second(op, prm, weightsCalculator); + return it->second(op, prm, weightsCalculator, pressureIndex); } PrecPtr doCreate(const Operator& op, const PropertyTree& prm, - const std::function weightsCalculator, const Comm& comm) + const std::function weightsCalculator, + std::size_t pressureIndex, const Comm& comm) { const std::string& type = prm.get("type", "ParOverILU0"); auto it = parallel_creators_.find(type); @@ -492,7 +516,7 @@ private: msg << std::endl; OPM_THROW(std::invalid_argument, msg.str()); } - return it->second(op, prm, weightsCalculator, comm); + return it->second(op, prm, weightsCalculator, pressureIndex, comm); } // Actually adds the creator. diff --git a/opm/simulators/linalg/PressureSolverPolicy.hpp b/opm/simulators/linalg/PressureSolverPolicy.hpp index dde25169e..684ed1c77 100644 --- a/opm/simulators/linalg/PressureSolverPolicy.hpp +++ b/opm/simulators/linalg/PressureSolverPolicy.hpp @@ -63,7 +63,10 @@ namespace Amg : linsolver_() { assert(op.category() == Dune::SolverCategory::overlapping); - linsolver_ = std::make_unique(op, comm, prm, std::function()); + // Assuming that we do not use Cpr as Pressure solver and use hard + // coded pressure index that might be wrong but should be unused. + linsolver_ = std::make_unique(op, comm, prm, std::function(), + /* pressureIndex = */ 1); } #endif // HAVE_MPI @@ -73,7 +76,10 @@ namespace Amg : linsolver_() { assert(op.category() != Dune::SolverCategory::overlapping); - linsolver_ = std::make_unique(op, prm, std::function()); + // Assuming that we do not use Cpr as Pressure solver and use hard + // coded pressure index that might be wrong but should be unused. + linsolver_ = std::make_unique(op, prm, std::function(), + /* pressureIndex = */ 1); } diff --git a/opm/simulators/linalg/PressureTransferPolicy.hpp b/opm/simulators/linalg/PressureTransferPolicy.hpp index 7d6aad760..b05640acb 100644 --- a/opm/simulators/linalg/PressureTransferPolicy.hpp +++ b/opm/simulators/linalg/PressureTransferPolicy.hpp @@ -150,10 +150,14 @@ public: return *coarseLevelCommunication_; } + std::size_t getPressureIndex() const + { + return pressure_var_index_; + } private: Communication* communication_; const FineVectorType& weights_; - const int pressure_var_index_; + const std::size_t pressure_var_index_; std::shared_ptr coarseLevelCommunication_; std::shared_ptr coarseLevelMatrix_; }; diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 7dbe69e33..b48df879f 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -108,7 +108,6 @@ setupCPR(const std::string& conf, const FlowLinearSolverParameters& p) } prm.put("preconditioner.finesmoother.type", "ParOverILU0"s); prm.put("preconditioner.finesmoother.relaxation", 1.0); - prm.put("preconditioner.pressure_var_index", 1); prm.put("preconditioner.verbosity", 0); prm.put("preconditioner.coarsesolver.maxiter", 1); prm.put("preconditioner.coarsesolver.tol", 1e-1); diff --git a/tests/options_flexiblesolver.json b/tests/options_flexiblesolver.json index 435584d12..357c42bb6 100644 --- a/tests/options_flexiblesolver.json +++ b/tests/options_flexiblesolver.json @@ -25,8 +25,7 @@ "solver": "bicgstab" }, "verbosity": "11", - "weights_filename" : "weight_cpr.txt", - "pressure_var_index" : "1" + "weights_filename" : "weight_cpr.txt" }, "verbosity": "10", "solver": "bicgstab" diff --git a/tests/test_flexiblesolver.cpp b/tests/test_flexiblesolver.cpp index c37781a7d..5f0d832d2 100644 --- a/tests/test_flexiblesolver.cpp +++ b/tests/test_flexiblesolver.cpp @@ -71,12 +71,12 @@ testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, con { return Opm::Amg::getQuasiImpesWeights(matrix, - prm.get("preconditioner.pressure_var_index"), + 1, transpose); }; using SeqOperatorType = Dune::MatrixAdapter; SeqOperatorType op(matrix); - Dune::FlexibleSolver solver(op, prm, wc); + Dune::FlexibleSolver solver(op, prm, wc, 1); Vector x(rhs.size()); Dune::InverseOperatorResult res; solver.apply(x, rhs, res); diff --git a/tests/test_preconditionerfactory.cpp b/tests/test_preconditionerfactory.cpp index 25c065689..a3944da05 100644 --- a/tests/test_preconditionerfactory.cpp +++ b/tests/test_preconditionerfactory.cpp @@ -100,10 +100,10 @@ testPrec(const Opm::PropertyTree& prm, const std::string& matrix_filename, const { return Opm::Amg::getQuasiImpesWeights(matrix, - prm.get("preconditioner.pressure_var_index"), + 1, transpose); }; - auto prec = PrecFactory::create(op, prm.get_child("preconditioner"), wc); + auto prec = PrecFactory::create(op, prm.get_child("preconditioner"), wc, 1); Dune::BiCGSTABSolver solver(op, *prec, prm.get("tol"), prm.get("maxiter"), prm.get("verbosity")); Vector x(rhs.size()); Dune::InverseOperatorResult res; @@ -191,7 +191,8 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner) // Add preconditioner to factory for block size 1. - PF<1>::addCreator("nothing", [](const O<1>&, const Opm::PropertyTree&, const std::function()>&) { + PF<1>::addCreator("nothing", [](const O<1>&, const Opm::PropertyTree&, const std::function()>&, + std::size_t) { return Dune::wrapPreconditioner>>(); }); @@ -206,7 +207,8 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner) } // Add preconditioner to factory for block size 3. - PF<3>::addCreator("nothing", [](const O<3>&, const Opm::PropertyTree&, const std::function()>&) { + PF<3>::addCreator("nothing", [](const O<3>&, const Opm::PropertyTree&, const std::function()>&, + std::size_t) { return Dune::wrapPreconditioner>>(); }); @@ -298,7 +300,8 @@ testPrecRepeating(const Opm::PropertyTree& prm, const std::string& matrix_filena using PrecFactory = Opm::PreconditionerFactory; // Add no-oppreconditioner to factory for block size 1. - PrecFactory::addCreator("nothing", [](const Operator&, const Opm::PropertyTree&, const std::function&) { + PrecFactory::addCreator("nothing", [](const Operator&, const Opm::PropertyTree&, const std::function&, + std::size_t) { return Dune::wrapPreconditioner>(); });