OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-28 02:00:59 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2694 from atgeirr/flexible-wellop

Browse Source 
Allow well operators with FlexibleSolver.
This commit is contained in:
Markus Blatt 2020-07-01 20:30:44 +02:00 committed by GitHub
commit c01e6d14f7
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 225 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -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<MatrixType, VectorType, VectorType>;
/// Base class type of the contained preconditioner.
using AbstractPrecondType = Dune::PreconditionerWithUpdate<VectorType, VectorType>;
/// Create a sequential solver.
FlexibleSolver(const MatrixType& matrix,
FlexibleSolver(AbstractOperatorType& op,
const boost::property_tree::ptree& prm,
const std::function<VectorTypeT()>& weightsCalculator = std::function<VectorTypeT()>());
const std::function<VectorType()>& weightsCalculator = std::function<VectorType()>());
/// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
template <class Comm>
FlexibleSolver(const MatrixType& matrix,
FlexibleSolver(AbstractOperatorType& op,
const Comm& comm,
const boost::property_tree::ptree& prm,
const std::function<VectorTypeT()>& weightsCalculator = std::function<VectorTypeT()>());
const std::function<VectorType()>& weightsCalculator = std::function<VectorType()>());
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<VectorType, VectorType>;
/// Access the contained preconditioner.
AbstractPrecondType& preconditioner();
virtual Dune::SolverCategory::Category category() const override;
private:
using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
using AbstractScalarProductType = Dune::ScalarProduct<VectorType>;
using AbstractSolverType = Dune::InverseOperator<VectorType, VectorType>;
// Machinery for making sequential or parallel operators/preconditioners/scalar products.
template <class Comm>
void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
const std::function<VectorTypeT()> weightsCalculator, const Comm& comm);
void initOpPrecSp(AbstractOperatorType& op, const boost::property_tree::ptree& prm,
const std::function<VectorType()> weightsCalculator, const Comm& comm);
void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
const std::function<VectorTypeT()> weightsCalculator, const Dune::Amg::SequentialInformation&);
void initOpPrecSp(AbstractOperatorType& op, const boost::property_tree::ptree& prm,
const std::function<VectorType()> 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 <class Comm>
void init(const MatrixType& matrix,
void init(AbstractOperatorType& op,
const Comm& comm,
const boost::property_tree::ptree& prm,
const std::function<VectorTypeT()> weightsCalculator);
const std::function<VectorType()> weightsCalculator);
std::shared_ptr<AbstractOperatorType> linearoperator_;
AbstractOperatorType* linearoperator_for_solver_;
std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;
std::shared_ptr<AbstractPrecondType> preconditioner_;
std::shared_ptr<AbstractScalarProductType> scalarproduct_;
std::shared_ptr<AbstractSolverType> linsolver_;

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

@ -39,23 +39,23 @@ namespace Dune
/// Create a sequential solver.
template <class MatrixType, class VectorType>
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver(const MatrixType& matrix,
FlexibleSolver(AbstractOperatorType& op,
const boost::property_tree::ptree& prm,
const std::function<VectorType()>& 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 <class MatrixType, class VectorType>
template <class Comm>
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver(const MatrixType& matrix,
FlexibleSolver(AbstractOperatorType& op,
const Comm& comm,
const boost::property_tree::ptree& prm,
const std::function<VectorType()>& weightsCalculator)
{
init(matrix, comm, prm, weightsCalculator);
init(op, comm, prm, weightsCalculator);
}
template <class MatrixType, class VectorType>
@ -88,7 +88,7 @@ namespace Dune
FlexibleSolver<MatrixType, VectorType>::
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 <class Comm>
void
FlexibleSolver<MatrixType, VectorType>::
initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
const std::function<VectorType()> weightsCalculator, const Comm& comm)
initOpPrecSp(AbstractOperatorType& op,
const boost::property_tree::ptree& prm,
const std::function<VectorType()> weightsCalculator,
const Comm& comm)
{
// Parallel case.
using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
using pt = const boost::property_tree::ptree;
auto linop = std::make_shared<ParOperatorType>(matrix, comm);
linearoperator_ = linop;
using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
linearoperator_for_solver_ = &op;
auto op_prec = std::make_shared<ParOperatorType>(op.getmat(), comm);
auto child = prm.get_child_optional("preconditioner");
preconditioner_
= Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, child? *child : pt(),
weightsCalculator, comm);
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, linearoperator_->category());
preconditioner_ = Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*op_prec,
child ? *child : pt(),
weightsCalculator,
comm);
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
linearoperator_for_precond_ = op_prec;
}
template <class MatrixType, class VectorType>
void
FlexibleSolver<MatrixType, VectorType>::
initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
const std::function<VectorType()> weightsCalculator, const Dune::Amg::SequentialInformation&)
initOpPrecSp(AbstractOperatorType& op,
const boost::property_tree::ptree& prm,
const std::function<VectorType()> weightsCalculator,
const Dune::Amg::SequentialInformation&)
{
// Sequential case.
using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
using pt = const boost::property_tree::ptree;
auto linop = std::make_shared<SeqOperatorType>(matrix);
linearoperator_ = linop;
using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
linearoperator_for_solver_ = &op;
auto op_prec = std::make_shared<SeqOperatorType>(op.getmat());
auto child = prm.get_child_optional("preconditioner");
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*linop, child? *child : pt(),
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*op_prec,
child ? *child : pt(),
weightsCalculator);
scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
linearoperator_for_precond_ = op_prec;
}
template <class MatrixType, class VectorType>
void
FlexibleSolver<MatrixType, VectorType>::
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<double>("tol", 1e-2);
const int maxiter = prm.get<int>("maxiter", 200);
const int verbosity = isMaster? prm.get<int>("verbosity", 0) : 0;
const int verbosity = is_iorank ? prm.get<int>("verbosity", 0) : 0;
const std::string solver_type = prm.get<std::string>("solver", "bicgstab");
if (solver_type == "bicgstab") {
linsolver_.reset(new Dune::BiCGSTABSolver<VectorType>(*linearoperator_,
linsolver_.reset(new Dune::BiCGSTABSolver<VectorType>(*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<VectorType>(*linearoperator_,
linsolver_.reset(new Dune::LoopSolver<VectorType>(*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<int>("restart", 15);
linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*linearoperator_,
linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*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<MatrixType>(linearoperator_->getmat(), verbosity, dummy));
linsolver_.reset(new Dune::UMFPack<MatrixType>(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 <class Comm>
void
FlexibleSolver<MatrixType, VectorType>::
init(const MatrixType& matrix,
init(AbstractOperatorType& op,
const Comm& comm,
const boost::property_tree::ptree& prm,
const std::function<VectorType()> 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<Dune::FieldMatrix<double, N, N>>;
template <int N>
using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<double, N, N>>;
// 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<int, int>;
#define INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n) \
template Dune::FlexibleSolver<OBM<n>, BV<n>>::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<BM<N>, BV<N>>; \
template class Dune::FlexibleSolver<OBM<N>, BV<N>>; \
template Dune::FlexibleSolver<BM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op, \
const Comm& comm, \
const boost::property_tree::ptree& prm, \
const std::function<BV<N>()>& weightsCalculator); \
template Dune::FlexibleSolver<OBM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op, \
const Comm& comm, \
const boost::property_tree::ptree& prm, \
const std::function<BV<n>()>& weightsCalculator);
#else
#define INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n)
#endif
const std::function<BV<N>()>& weightsCalculator);
// INSTANTIATE instantiates the class including any templated constructors if necessary.
#define INSTANTIATE_FLEXIBLESOLVER(n) \
/* Variants using Dune::FieldMatrix blocks. */ \
template class Dune::FlexibleSolver<BM<n>, BV<n>>; \
/* Variants using Opm::MatrixBlock blocks. */ \
template class Dune::FlexibleSolver<OBM<n>, BV<n>>; \
INSTANTIATE_FLEXIBLESOLVER_CONSTRUCTOR(n)
#else // HAVE_MPI
#define INSTANTIATE_FLEXIBLESOLVER(N) \
template class Dune::FlexibleSolver<BM<N>, BV<N>>; \
template class Dune::FlexibleSolver<OBM<N>, BV<N>>;
#endif // HAVE_MPI
#endif // OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED

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

@ -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<Matrix, Vector>;
using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
using WellModelOperator = WellModelAsLinearOperator<WellModel, Vector, Vector>;
// 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<Matrix, Vector, Vector, Comm>;
linearOperatorForFlexibleSolver_ = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*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<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);
}
#endif
} else {
flexibleSolver_.reset(new FlexibleSolverType(getMatrix(), prm_, weightsCalculator));
using SeqLinearOperator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
linearOperatorForFlexibleSolver_ = std::make_unique<SeqLinearOperator>(getMatrix());
flexibleSolver_ = std::make_unique<FlexibleSolverType>(*linearOperatorForFlexibleSolver_, prm_, weightsCalculator);
}
}
else
@ -1004,6 +1014,8 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
Vector *rhs_;
std::unique_ptr<FlexibleSolverType> flexibleSolver_;
std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
std::vector<int> overlapRows_;
std::vector<int> interiorRows_;
std::vector<std::set<int>> wellConnectionsGraph_;

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

@ -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<int, int>;
#else
using Communication = int; // Dummy type.
#endif
using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
using WellModelOpType = WellModelAsLinearOperator<WellModel, VectorType, VectorType>;
using SolverType = Dune::FlexibleSolver<MatrixType, VectorType>;
// 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<TypeTag>();
prm_ = setupPropertyTree<TypeTag>(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<VectorType()> weightsCalculator;
matrix_ = &mat.istlMatrix(); // Store pointer for output if needed.
std::function<VectorType()> 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<MatrixType,
VectorType>(
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<MatrixType, VectorType, VectorType, Communication>;
auto op = std::make_unique<ParOperatorType>(mat.istlMatrix(), *comm_);
auto sol = std::make_unique<SolverType>(*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<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);
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<MatrixType, VectorType, VectorType>;
auto op = std::make_unique<SeqOperatorType>(mat.istlMatrix());
auto sol = std::make_unique<SolverType>(*op, prm_, weightsCalculator);
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);
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<VectorType()> getWeightsCalculator(const MatrixType& mat, const VectorType& b) const
{
std::function<VectorType()> 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<MatrixType, VectorType>(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<WellModelOpType> well_operator_;
std::unique_ptr<AbstractOperatorType> linear_operator_;
std::unique_ptr<SolverType> 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<Communication> comm_;
std::vector<int> overlapRows_;
std::vector<int> interiorRows_;

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

@ -25,6 +25,7 @@
#include <boost/property_tree/ptree.hpp>
#include <dune/istl/solver.hh>
#include <dune/istl/owneroverlapcopy.hh>
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<X, X> {
template <class Comm>
PressureInverseOperator(Operator& op, const boost::property_tree::ptree& prm, const Comm& comm)
struct PressureInverseOperator : public Dune::InverseOperator<X, X>
{
template <typename GlobalIndex, typename LocalIndex>
PressureInverseOperator(Operator& op,
const boost::property_tree::ptree& prm,
const Dune::OwnerOverlapCopyCommunication<GlobalIndex, LocalIndex>& comm)
: linsolver_()
{
assert(op.category() == Dune::SolverCategory::overlapping);
linsolver_ = std::make_unique<Solver>(op.getmat(), comm, prm, std::function<X()>());
linsolver_ = std::make_unique<Solver>(op, comm, prm, std::function<X()>());
}
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<Solver>(op.getmat(), prm, std::function<X()>());
linsolver_ = std::make_unique<Solver>(op, prm, std::function<X()>());
}

4
tests/test_flexiblesolver.cpp
View File

@ -75,7 +75,9 @@ testSolver(const boost::property_tree::ptree& prm, const std::string& matrix_fil
prm.get<int>("preconditioner.pressure_var_index"),
transpose);
};
Dune::FlexibleSolver<Matrix, Vector> solver(matrix, prm, wc);
using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
SeqOperatorType op(matrix);
Dune::FlexibleSolver<Matrix, Vector> solver(op, prm, wc);
Vector x(rhs.size());
Dune::InverseOperatorResult res;
solver.apply(x, rhs, res);