OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Determine index of pressure from model used.

Browse Source 
Previously, the user had to specify it in the json file read from the
FlexibleSolver or 1 was used. Unfortunately, the index depends on the
model used and it seem rather opaque to a user what that index is.

With this commit we determine the pressure index from the model.
This commit is contained in:
Markus Blatt
2021-06-08 21:51:34 +02:00
parent ed131190eb
commit db0f19ba88
12 changed files with 160 additions and 93 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
opm/simulators/linalg/FlexibleSolver.hpp
View File

@@ -50,14 +50,16 @@ public:
/// Create a sequential solver.
FlexibleSolver(AbstractOperatorType& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()>& weightsCalculator = std::function<VectorType()>());
const std::function<VectorType()>& weightsCalculator,
std::size_t pressureIndex);
/// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
template <class Comm>
FlexibleSolver(AbstractOperatorType& op,
const Comm& comm,
const Opm::PropertyTree& prm,
const std::function<VectorType()>& weightsCalculator = std::function<VectorType()>());
const std::function<VectorType()>& 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 <class Comm>
void initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator, const Comm& comm);
const std::function<VectorType()> weightsCalculator, const Comm& comm,
std::size_t pressureIndex);
void initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator, const Dune::Amg::SequentialInformation&);
const std::function<VectorType()> 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<VectorType()> weightsCalculator);
const std::function<VectorType()> weightsCalculator,
std::size_t pressureIndex);
AbstractOperatorType* linearoperator_for_solver_;
std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;

34
opm/simulators/linalg/FlexibleSolver_impl.hpp
View File

@@ -39,9 +39,11 @@ namespace Dune
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver(AbstractOperatorType& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()>& weightsCalculator)
const std::function<VectorType()>& 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<VectorType()>& weightsCalculator)
const std::function<VectorType()>& weightsCalculator,
std::size_t pressureIndex)
{
init(op, comm, prm, weightsCalculator);
init(op, comm, prm, weightsCalculator, pressureIndex);
}
template <class MatrixType, class VectorType>
@@ -97,7 +100,8 @@ namespace Dune
initOpPrecSp(AbstractOperatorType& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
const Comm& comm)
const Comm& comm,
std::size_t pressureIndex)
{
// Parallel case.
using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
@@ -107,7 +111,8 @@ namespace Dune
preconditioner_ = Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*op_prec,
child ? *child : Opm::PropertyTree(),
weightsCalculator,
comm);
comm,
pressureIndex);
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
linearoperator_for_precond_ = op_prec;
}
@@ -118,7 +123,8 @@ namespace Dune
initOpPrecSp(AbstractOperatorType& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
const Dune::Amg::SequentialInformation&)
const Dune::Amg::SequentialInformation&,
std::size_t pressureIndex)
{
// Sequential case.
using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
@@ -127,7 +133,8 @@ namespace Dune
auto child = prm.get_child_optional("preconditioner");
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*op_prec,
child ? *child : Opm::PropertyTree(),
weightsCalculator);
weightsCalculator,
pressureIndex);
scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
linearoperator_for_precond_ = op_prec;
}
@@ -184,9 +191,10 @@ namespace Dune
init(AbstractOperatorType& op,
const Comm& comm,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator)
const std::function<VectorType()> 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<OBM<N>, BV<N>>; \
template Dune::FlexibleSolver<BM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op, \
const Comm& comm, \
const Opm::PropertyTree& prm, \
const std::function<BV<N>()>& weightsCalculator); \
const std::function<BV<N>()>& weightsCalculator, \
std::size_t); \
template Dune::FlexibleSolver<OBM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op, \
const Comm& comm, \
const Opm::PropertyTree& prm, \
const std::function<BV<N>()>& weightsCalculator);
const std::function<BV<N>()>& weightsCalculator, \
std::size_t);
#else // HAVE_MPI

28
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@@ -91,6 +91,7 @@ namespace Opm
using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
using WellModelOperator = WellModelAsLinearOperator<WellModel, Vector, Vector>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::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<Matrix, Vector, Vector, Comm>;
linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator);
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
pressureIndex);
} else {
using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator);
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, *comm_, prm_, weightsCalculator,
pressureIndex);
}
#endif
} else {
if (useWellConn_) {
using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix());
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator);
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator,
pressureIndex);
} else {
using SeqOperatorType = WellModelMatrixAdapter<Matrix, Vector, Vector, false>;
wellOperator_ = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
linearOperatorForFlexibleSolver_ = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator);
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*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<Matrix, Vector>(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<Matrix, Vector>(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,

22
opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
View File

@@ -83,6 +83,7 @@ class ISTLSolverEbosFlexible
using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
typedef typename GridView::template Codim<0>::Entity Element;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
public:
static void registerParameters()
@@ -156,7 +157,8 @@ public:
if (matrixAddWellContributions_) {
using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Communication>;
auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *comm_);
auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator);
auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator,
pressureIndex);
solver_ = std::move(sol);
linear_operator_ = std::move(op);
} else {
@@ -167,7 +169,8 @@ public:
using ParOperatorType = WellModelGhostLastMatrixAdapter<MatrixType, VectorType, VectorType, true>;
auto well_op = std::make_unique<WellModelOpType>(simulator_.problem().wellModel());
auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *well_op, interiorCellNum_);
auto sol = std::make_unique<SolverType>(*op, *comm_, prm_, weightsCalculator);
auto sol = std::make_unique<SolverType>(*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<MatrixType, VectorType, VectorType>;
auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix());
auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator);
auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator,
pressureIndex);
solver_ = std::move(sol);
linear_operator_ = std::move(op);
} else {
using SeqOperatorType = WellModelMatrixAdapter<MatrixType, VectorType, VectorType, false>;
auto well_op = std::make_unique<WellModelOpType>(simulator_.problem().wellModel());
auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix(), *well_op);
auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator);
auto sol = std::make_unique<SolverType>(*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<MatrixType, VectorType>(mat, pressureIndex, transpose);
weightsCalculator = [&mat, transpose, p = pressureIndex]() {
return Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(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,

17
opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp
View File

@@ -84,15 +84,17 @@ public:
using PrecFactory = Opm::PreconditionerFactory<OperatorType, Communication>;
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator)
const std::function<VectorType()> 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<VectorType()>(), pressureIndex))
, comm_(nullptr)
, weightsCalculator_(weightsCalculator)
, weights_(weightsCalculator())
, levelTransferPolicy_(dummy_comm_, weights_, prm.get<int>("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<VectorType()> weightsCalculator, const Communication& comm)
const std::function<VectorType()> 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<VectorType()>(),
comm, pressureIndex))
, comm_(&comm)
, weightsCalculator_(weightsCalculator)
, weights_(weightsCalculator())
, levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index", 1))
, levelTransferPolicy_(*comm_, weights_, pressureIndex)
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : Opm::PropertyTree())
, twolevel_method_(linearoperator,
finesmoother_,

100
opm/simulators/linalg/PreconditionerFactory.hpp
View File

@@ -36,6 +36,7 @@
#include <map>
#include <memory>
#include <limits>
namespace Opm
{
@@ -57,8 +58,10 @@ public:
using PrecPtr = std::shared_ptr<Dune::PreconditionerWithUpdate<Vector, Vector>>;
/// The type of creator functions passed to addCreator().
using Creator = std::function<PrecPtr(const Operator&, const PropertyTree&, const std::function<Vector()>&)>;
using ParCreator = std::function<PrecPtr(const Operator&, const PropertyTree&, const std::function<Vector()>&, const Comm&)>;
using Creator = std::function<PrecPtr(const Operator&, const PropertyTree&,
const std::function<Vector()>&, std::size_t)>;
using ParCreator = std::function<PrecPtr(const Operator&, const PropertyTree&,
const std::function<Vector()>&, 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<Vector()>& weightsCalculator = std::function<Vector()>())
const std::function<Vector()>& weightsCalculator = {},
std::size_t pressureIndex = std::numeric_limits<std::size_t>::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<Vector()>& weightsCalculator, const Comm& comm)
const std::function<Vector()>& weightsCalculator, const Comm& comm,
std::size_t pressureIndex = std::numeric_limits<std::size_t>::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<std::size_t>::max())
{
return instance().doCreate(op, prm, std::function<Vector()>(), comm);
return instance().doCreate(op, prm, std::function<Vector()>(), 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<Vector()>&, const C& comm) {
doAddCreator("ILU0", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t, const C& comm) {
return createParILU(op, prm, comm, 0);
});
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t, const C& comm) {
return createParILU(op, prm, comm, prm.get<int>("ilulevel", 0));
});
doAddCreator("ILUn", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
doAddCreator("ILUn", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t, const C& comm) {
return createParILU(op, prm, comm, prm.get<int>("ilulevel", 0));
});
doAddCreator("Jac", [](const O& op, const P& prm, const std::function<Vector()>&,
const C& comm) {
std::size_t, const C& comm) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqJac<M, V, V>>>(comm, op.getmat(), n, w);
});
doAddCreator("GS", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
doAddCreator("GS", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t, const C& comm) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqGS<M, V, V>>>(comm, op.getmat(), n, w);
});
doAddCreator("SOR", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
doAddCreator("SOR", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t, const C& comm) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSOR<M, V, V>>>(comm, op.getmat(), n, w);
});
doAddCreator("SSOR", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
doAddCreator("SSOR", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t, const C& comm) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSSOR<M, V, V>>>(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<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t, const C& comm) {
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
if (smoother == "ILU0" || smoother == "ParOverILU0") {
using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
@@ -306,13 +312,21 @@ private:
});
}
doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, const C& comm) {
doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, std::size_t pressureIndex, const C& comm) {
assert(weightsCalculator);
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, weightsCalculator, comm);
if (pressureIndex == std::numeric_limits<std::size_t>::max())
{
OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
}
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, weightsCalculator, pressureIndex, comm);
});
doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, const C& comm) {
doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, std::size_t pressureIndex, const C& comm) {
assert(weightsCalculator);
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, weightsCalculator, comm);
if (pressureIndex == std::numeric_limits<std::size_t>::max())
{
OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
}
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(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<Vector()>&) {
doAddCreator("ILU0", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const double w = prm.get<double>("relaxation", 1.0);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<Vector()>&) {
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const double w = prm.get<double>("relaxation", 1.0);
const int n = prm.get<int>("ilulevel", 0);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("ILUn", [](const O& op, const P& prm, const std::function<Vector()>&) {
doAddCreator("ILUn", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const int n = prm.get<int>("ilulevel", 0);
const double w = prm.get<double>("relaxation", 1.0);
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
});
doAddCreator("Jac", [](const O& op, const P& prm, const std::function<Vector()>&) {
doAddCreator("Jac", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqJac<M, V, V>>(op.getmat(), n, w);
});
doAddCreator("GS", [](const O& op, const P& prm, const std::function<Vector()>&) {
doAddCreator("GS", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqGS<M, V, V>>(op.getmat(), n, w);
});
doAddCreator("SOR", [](const O& op, const P& prm, const std::function<Vector()>&) {
doAddCreator("SOR", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqSOR<M, V, V>>(op.getmat(), n, w);
});
doAddCreator("SSOR", [](const O& op, const P& prm, const std::function<Vector()>&) {
doAddCreator("SSOR", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqSSOR<M, V, V>>(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<O, Dune::MatrixAdapter<M, V, V>>) {
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&) {
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const std::string smoother = prm.get<std::string>("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<Vector()>&) {
doAddCreator("kamg", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
const std::string smoother = prm.get<std::string>("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<Vector()>&) {
doAddCreator("famg", [](const O& op, const P& prm, const std::function<Vector()>&, std::size_t) {
auto crit = amgCriterion(prm);
Dune::Amg::Parameters parms;
parms.setNoPreSmoothSteps(1);
@@ -435,11 +449,19 @@ private:
return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
});
}
doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm, weightsCalculator);
doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator, std::size_t pressureIndex) {
if (pressureIndex == std::numeric_limits<std::size_t>::max())
{
OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
}
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm, weightsCalculator, pressureIndex);
});
doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(op, prm, weightsCalculator);
doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator, std::size_t pressureIndex) {
if (pressureIndex == std::numeric_limits<std::size_t>::max())
{
OPM_THROW(std::logic_error, "Pressure index out of bounds. It needs to specified for CPR");
}
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(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<Vector()> weightsCalculator)
const std::function<Vector()> weightsCalculator,
std::size_t pressureIndex)
{
const std::string& type = prm.get<std::string>("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<Vector()> weightsCalculator, const Comm& comm)
const std::function<Vector()> weightsCalculator,
std::size_t pressureIndex, const Comm& comm)
{
const std::string& type = prm.get<std::string>("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.

10
opm/simulators/linalg/PressureSolverPolicy.hpp
View File

@@ -63,7 +63,10 @@ namespace Amg
: linsolver_()
{
assert(op.category() == Dune::SolverCategory::overlapping);
linsolver_ = std::make_unique<Solver>(op, comm, prm, std::function<X()>());
// 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<Solver>(op, comm, prm, std::function<X()>(),
/* pressureIndex = */ 1);
}
#endif // HAVE_MPI
@@ -73,7 +76,10 @@ namespace Amg
: linsolver_()
{
assert(op.category() != Dune::SolverCategory::overlapping);
linsolver_ = std::make_unique<Solver>(op, prm, std::function<X()>());
// 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<Solver>(op, prm, std::function<X()>(),
/* pressureIndex = */ 1);
}

6
opm/simulators/linalg/PressureTransferPolicy.hpp
View File

@@ -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<Communication> coarseLevelCommunication_;
std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
};

1
opm/simulators/linalg/setupPropertyTree.cpp
View File

@@ -106,7 +106,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);

3
tests/options_flexiblesolver.json
View File

@@ -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"

4
tests/test_flexiblesolver.cpp
View File

@@ -71,12 +71,12 @@ testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, con
{
return Opm::Amg::getQuasiImpesWeights<Matrix,
Vector>(matrix,
prm.get<int>("preconditioner.pressure_var_index"),
1,
transpose);
};
using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
SeqOperatorType op(matrix);
Dune::FlexibleSolver<Matrix, Vector> solver(op, prm, wc);
Dune::FlexibleSolver<Matrix, Vector> solver(op, prm, wc, 1);
Vector x(rhs.size());
Dune::InverseOperatorResult res;
solver.apply(x, rhs, res);

13
tests/test_preconditionerfactory.cpp
View File

@@ -100,10 +100,10 @@ testPrec(const Opm::PropertyTree& prm, const std::string& matrix_filename, const
{
return Opm::Amg::getQuasiImpesWeights<Matrix,
Vector>(matrix,
prm.get<int>("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<Vector> solver(op, *prec, prm.get<double>("tol"), prm.get<int>("maxiter"), prm.get<int>("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<V<1>()>&) {
PF<1>::addCreator("nothing", [](const O<1>&, const Opm::PropertyTree&, const std::function<V<1>()>&,
std::size_t) {
return Dune::wrapPreconditioner<NothingPreconditioner<V<1>>>();
});
@@ -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<V<3>()>&) {
PF<3>::addCreator("nothing", [](const O<3>&, const Opm::PropertyTree&, const std::function<V<3>()>&,
std::size_t) {
return Dune::wrapPreconditioner<NothingPreconditioner<V<3>>>();
});
@@ -298,7 +300,8 @@ testPrecRepeating(const Opm::PropertyTree& prm, const std::string& matrix_filena
using PrecFactory = Opm::PreconditionerFactory<Operator>;
// Add no-oppreconditioner to factory for block size 1.
PrecFactory::addCreator("nothing", [](const Operator&, const Opm::PropertyTree&, const std::function<Vector()>&) {
PrecFactory::addCreator("nothing", [](const Operator&, const Opm::PropertyTree&, const std::function<Vector()>&,
std::size_t) {
return Dune::wrapPreconditioner<NothingPreconditioner<Vector>>();
});