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

Adding functions for pressure system manipulation to well models.

Browse Source 
Also add well-aware operator and transfer policy.
This will be used for CPR with custom operators.
This commit is contained in:
hnil 2021-01-04 15:18:23 +01:00 committed by Atgeirr Flø Rasmussen
parent 0bb293aeb0
commit f3acfcde0b
15 changed files with 821 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
CMakeLists_files.cmake
View File

@ -42,12 +42,12 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.cpp
opm/simulators/flow/ValidationFunctions.cpp
opm/simulators/linalg/ExtractParallelGridInformationToISTL.cpp
opm/simulators/linalg/FlexibleSolver1.cpp
opm/simulators/linalg/FlexibleSolver2.cpp
opm/simulators/linalg/FlexibleSolver3.cpp
opm/simulators/linalg/FlexibleSolver4.cpp
opm/simulators/linalg/FlexibleSolver5.cpp
opm/simulators/linalg/FlexibleSolver6.cpp
# opm/simulators/linalg/FlexibleSolver1.cpp
# opm/simulators/linalg/FlexibleSolver2.cpp
# opm/simulators/linalg/FlexibleSolver3.cpp
# opm/simulators/linalg/FlexibleSolver4.cpp
# opm/simulators/linalg/FlexibleSolver5.cpp
# opm/simulators/linalg/FlexibleSolver6.cpp
opm/simulators/linalg/PropertyTree.cpp
opm/simulators/linalg/setupPropertyTree.cpp
opm/simulators/utils/PartiallySupportedFlowKeywords.cpp

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

@ -35,27 +35,25 @@ namespace Dune
/// (operator, scalar product, iterative solver and preconditioner) and sets
/// them up based on runtime parameters, using the PreconditionerFactory for
/// setting up preconditioners.
template <class MatrixTypeT, class VectorTypeT>
class FlexibleSolver : public Dune::InverseOperator<VectorTypeT, VectorTypeT>
template <class Operator>
class FlexibleSolver : public Dune::InverseOperator<typename Operator::domain_type,
typename Operator::range_type>
{
public:
using MatrixType = MatrixTypeT;
using VectorType = VectorTypeT;
using VectorType = typename Operator::domain_type; // Assuming symmetry: domain == range
/// 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(AbstractOperatorType& op,
FlexibleSolver(Operator& op,
const Opm::PropertyTree& prm,
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,
FlexibleSolver(Operator& op,
const Comm& comm,
const Opm::PropertyTree& prm,
const std::function<VectorType()>& weightsCalculator,
@ -76,11 +74,11 @@ private:
// Machinery for making sequential or parallel operators/preconditioners/scalar products.
template <class Comm>
void initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm,
void initOpPrecSp(Operator& op, const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator, const Comm& comm,
std::size_t pressureIndex);
void initOpPrecSp(AbstractOperatorType& op, const Opm::PropertyTree& prm,
void initOpPrecSp(Operator& op, const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator, const Dune::Amg::SequentialInformation&,
std::size_t pressureIndex);
@ -89,14 +87,13 @@ private:
// Main initialization routine.
// Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
template <class Comm>
void init(AbstractOperatorType& op,
void init(Operator& op,
const Comm& comm,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
std::size_t pressureIndex);
AbstractOperatorType* linearoperator_for_solver_;
std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;
Operator* linearoperator_for_solver_;
std::shared_ptr<AbstractPrecondType> preconditioner_;
std::shared_ptr<AbstractScalarProductType> scalarproduct_;
std::shared_ptr<AbstractSolverType> linsolver_;
@ -104,6 +101,7 @@ private:
} // namespace Dune
#include <opm/simulators/linalg/FlexibleSolver_impl.hpp>
#endif // OPM_FLEXIBLE_SOLVER_HEADER_INCLUDED

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

@ -36,9 +36,9 @@
namespace Dune
{
/// Create a sequential solver.
template <class MatrixType, class VectorType>
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver(AbstractOperatorType& op,
template <class Operator>
FlexibleSolver<Operator>::
FlexibleSolver(Operator& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()>& weightsCalculator,
std::size_t pressureIndex)
@ -48,10 +48,10 @@ namespace Dune
}
/// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
template <class MatrixType, class VectorType>
template <class Operator>
template <class Comm>
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver(AbstractOperatorType& op,
FlexibleSolver<Operator>::
FlexibleSolver(Operator& op,
const Comm& comm,
const Opm::PropertyTree& prm,
const std::function<VectorType()>& weightsCalculator,
@ -60,89 +60,83 @@ namespace Dune
init(op, comm, prm, weightsCalculator, pressureIndex);
}
template <class MatrixType, class VectorType>
template <class Operator>
void
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver<Operator>::
apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
{
linsolver_->apply(x, rhs, res);
}
template <class MatrixType, class VectorType>
template <class Operator>
void
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver<Operator>::
apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res)
{
linsolver_->apply(x, rhs, reduction, res);
}
/// Access the contained preconditioner.
template <class MatrixType, class VectorType>
template <class Operator>
auto
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver<Operator>::
preconditioner() -> AbstractPrecondType&
{
return *preconditioner_;
}
template <class MatrixType, class VectorType>
template <class Operator>
Dune::SolverCategory::Category
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver<Operator>::
category() const
{
return linearoperator_for_solver_->category();
}
// Machinery for making sequential or parallel operators/preconditioners/scalar products.
template <class MatrixType, class VectorType>
template <class Operator>
template <class Comm>
void
FlexibleSolver<MatrixType, VectorType>::
initOpPrecSp(AbstractOperatorType& op,
FlexibleSolver<Operator>::
initOpPrecSp(Operator& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
const Comm& comm,
std::size_t pressureIndex)
{
// Parallel case.
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(*op_prec,
child ? *child : Opm::PropertyTree(),
weightsCalculator,
comm,
pressureIndex);
preconditioner_ = Opm::PreconditionerFactory<Operator, Comm>::create(op,
child ? *child : Opm::PropertyTree(),
weightsCalculator,
comm,
pressureIndex);
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
linearoperator_for_precond_ = op_prec;
}
template <class MatrixType, class VectorType>
template <class Operator>
void
FlexibleSolver<MatrixType, VectorType>::
initOpPrecSp(AbstractOperatorType& op,
FlexibleSolver<Operator>::
initOpPrecSp(Operator& op,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
const Dune::Amg::SequentialInformation&,
std::size_t pressureIndex)
{
// Sequential case.
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(*op_prec,
child ? *child : Opm::PropertyTree(),
weightsCalculator,
pressureIndex);
preconditioner_ = Opm::PreconditionerFactory<Operator>::create(op,
child ? *child : Opm::PropertyTree(),
weightsCalculator,
pressureIndex);
scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
linearoperator_for_precond_ = op_prec;
}
template <class MatrixType, class VectorType>
template <class Operator>
void
FlexibleSolver<MatrixType, VectorType>::
FlexibleSolver<Operator>::
initSolver(const Opm::PropertyTree& prm, const bool is_iorank)
{
const double tol = prm.get<double>("tol", 1e-2);
@ -175,6 +169,7 @@ namespace Dune
#if HAVE_SUITESPARSE_UMFPACK
} else if (solver_type == "umfpack") {
bool dummy = false;
using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
linsolver_.reset(new Dune::UMFPack<MatrixType>(linearoperator_for_solver_->getmat(), verbosity, dummy));
#endif
} else {
@ -185,11 +180,11 @@ namespace Dune
// Main initialization routine.
// Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
template <class MatrixType, class VectorType>
template <class Operator>
template <class Comm>
void
FlexibleSolver<MatrixType, VectorType>::
init(AbstractOperatorType& op,
FlexibleSolver<Operator>::
init(Operator& op,
const Comm& comm,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
@ -203,7 +198,7 @@ namespace Dune
// Macros to simplify explicit instantiation of FlexibleSolver for various block sizes.
/*
template <int N>
using BV = Dune::BlockVector<Dune::FieldVector<double, N>>;
template <int N>
@ -240,6 +235,8 @@ template class Dune::FlexibleSolver<BM<N>, BV<N>>; \
template class Dune::FlexibleSolver<OBM<N>, BV<N>>;
#endif // HAVE_MPI
*/
#define INSTANTIATE_FLEXIBLESOLVER(N) do(){}
#endif // OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED

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

@ -87,9 +87,11 @@ namespace Opm
using Matrix = typename SparseMatrixAdapter::IstlMatrix;
using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FlexibleSolverType = Dune::FlexibleSolver<Matrix, Vector>;
using AbstractSolverType = Dune::InverseOperator<Vector, Vector>;
using AbstractOperatorType = Dune::AssembledLinearOperator<Matrix, Vector, Vector>;
using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<Vector, Vector>;
using WellModelOperator = WellModelAsLinearOperator<WellModel, Vector, Vector>;
using ParOperatorType = WellModelGhostLastMatrixAdapter<Matrix, Vector, Vector, true>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
@ -381,35 +383,47 @@ namespace Opm
#if HAVE_MPI
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,
pressureIndex);
auto op = std::make_unique<ParOperatorType>(getMatrix(), *comm_);
using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*op, *comm_, prm_, weightsCalculator, pressureIndex);
preconditionerForFlexibleSolver_ = &(sol->preconditioner());
linearOperatorForFlexibleSolver_ = std::move(op);
flexibleSolver_ = std::move(sol);
} 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,
pressureIndex);
auto op = std::make_unique<ParOperatorType>(getMatrix(), *wellOperator_, interiorCellNum_);
using FlexibleSolverType = Dune::FlexibleSolver<ParOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*op, *comm_, prm_, weightsCalculator, pressureIndex);
preconditionerForFlexibleSolver_ = &(sol->preconditioner());
linearOperatorForFlexibleSolver_ = std::move(op);
flexibleSolver_ = std::move(sol);
}
#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,
pressureIndex);
auto op = std::make_unique<SeqOperatorType>(getMatrix());
using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*op, prm_, weightsCalculator, pressureIndex);
preconditionerForFlexibleSolver_ = &(sol->preconditioner());
linearOperatorForFlexibleSolver_ = std::move(op);
flexibleSolver_ = std::move(sol);
} 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,
pressureIndex);
auto op = std::make_unique<SeqOperatorType>(getMatrix(), *wellOperator_);
using FlexibleSolverType = Dune::FlexibleSolver<SeqOperatorType>;
auto sol = std::make_unique<FlexibleSolverType>(*op, prm_, weightsCalculator, pressureIndex);
preconditionerForFlexibleSolver_ = &(sol->preconditioner());
linearOperatorForFlexibleSolver_ = std::move(op);
flexibleSolver_ = std::move(sol);
}
}
}
else
{
flexibleSolver_->preconditioner().update();
preconditionerForFlexibleSolver_->update();
}
}
@ -452,9 +466,11 @@ namespace Opm
using namespace std::string_literals;
auto preconditionerType = prm_.get("preconditioner.type"s, "cpr"s);
if (preconditionerType == "cpr" || preconditionerType == "cprt") {
const bool transpose = preconditionerType == "cprt";
if (preconditionerType == "cpr" || preconditionerType == "cprt"
|| preconditionerType == "cprw" || preconditionerType == "cprwt") {
const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
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") {
// weights will be created as default in the solver
// assignment p = pressureIndex prevent compiler warning about
@ -598,8 +614,9 @@ namespace Opm
std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
std::unique_ptr<FlexibleSolverType> flexibleSolver_;
std::unique_ptr<AbstractSolverType> flexibleSolver_;
std::unique_ptr<AbstractOperatorType> linearOperatorForFlexibleSolver_;
AbstractPreconditionerType* preconditionerForFlexibleSolver_;
std::unique_ptr<WellModelAsLinearOperator<WellModel, Vector, Vector>> wellOperator_;
std::vector<int> overlapRows_;
std::vector<int> interiorRows_;

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

@ -50,7 +50,7 @@ namespace Dune
// Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
template <class MatrixTypeT, class VectorTypeT>
template <class Operator>
class FlexibleSolver;
template <typename T, typename A, int i>
@ -183,7 +183,7 @@ private:
ParCoarseOperatorType>;
using LevelTransferPolicy = Opm::PressureTransferPolicy<OperatorType, CoarseOperatorType, Communication, transpose>;
using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperatorType,
FlexibleSolver<PressureMatrixType, PressureVectorType>,
FlexibleSolver<CoarseOperatorType>,
LevelTransferPolicy>;
using TwoLevelMethod
= Dune::Amg::TwoLevelMethodCpr<OperatorType, CoarseSolverPolicy, Dune::Preconditioner<VectorType, VectorType>>;

231
opm/simulators/linalg/OwningTwoLevelPreconditionerWell.hpp Normal file
View File

@ -0,0 +1,231 @@
/*
Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/simulators/linalg/PressureSolverPolicy.hpp>
#include <opm/simulators/linalg/PressureBhpTransferPolicy.hpp>
#include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
#include <opm/simulators/linalg/twolevelmethodcpr.hh>
#include <opm/common/ErrorMacros.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/paamg/amg.hh>
#include <boost/property_tree/ptree.hpp>
#include <fstream>
#include <type_traits>
namespace Opm
{
// Circular dependency between PreconditionerFactory [which can make an OwningTwoLevelPreconditioner]
// and OwningTwoLevelPreconditioner [which uses PreconditionerFactory to choose the fine-level smoother]
// must be broken, accomplished by forward-declaration here.
//template <class Operator, class Comm = Dune::Amg::SequentialInformation>
//class PreconditionerFactory;
}
namespace Dune
{
// Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
template <class Operator>
class FlexibleSolver;
// template <typename T, typename A, int i>
// std::ostream& operator<<(std::ostream& out,
// const BlockVector< FieldVector< T, i >, A >& vector)
// {
// Dune::writeMatrixMarket(vector, out);
// return out;
// }
// template <typename T, typename A, int i>
// std::istream& operator>>(std::istream& input,
// BlockVector< FieldVector< T, i >, A >& vector)
// {
// Dune::readMatrixMarket(vector, input);
// return input;
// }
/// A version of the two-level preconditioner that is:
/// - Self-contained, because it owns its policy components.
/// - Flexible, because it uses the runtime-flexible solver
/// and preconditioner factory.
template <class OperatorType,
class VectorType,
bool transpose = false,
class Communication = Dune::Amg::SequentialInformation>
class OwningTwoLevelPreconditionerWell : public Dune::PreconditionerWithUpdate<VectorType, VectorType>
{
public:
using pt = boost::property_tree::ptree;
using MatrixType = typename OperatorType::matrix_type;
using PrecFactory = Opm::PreconditionerFactory<OperatorType, Communication>;
using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
OwningTwoLevelPreconditionerWell(const OperatorType& linearoperator, const pt& prm,
const std::function<VectorType()> weightsCalculator)
: linear_operator_(linearoperator)
, finesmoother_(PrecFactory::create(linearoperator,
prm.get_child_optional("finesmoother")?
prm.get_child("finesmoother"): pt()))
, comm_(nullptr)
, weightsCalculator_(weightsCalculator)
, weights_(weightsCalculator())
, levelTransferPolicy_(dummy_comm_, weights_, prm)//.get<int>("pressure_var_index"))
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : pt())
, twolevel_method_(linearoperator,
finesmoother_,
levelTransferPolicy_,
coarseSolverPolicy_,
prm.get<int>("pre_smooth", transpose? 1 : 0),
prm.get<int>("post_smooth", transpose? 0 : 1))
, prm_(prm)
{
if (prm.get<int>("verbosity", 0) > 10) {
std::string filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
std::ofstream outfile(filename);
if (!outfile) {
OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
}
Dune::writeMatrixMarket(weights_, outfile);
}
}
OwningTwoLevelPreconditionerWell(const OperatorType& linearoperator,
const pt& prm,
const std::function<VectorType()> weightsCalculator, const Communication& comm)
: linear_operator_(linearoperator)
, finesmoother_(PrecFactory::create(linearoperator,
prm.get_child_optional("finesmoother")?
prm.get_child("finesmoother"): pt(), comm))
, comm_(&comm)
, weightsCalculator_(weightsCalculator)
, weights_(weightsCalculator())
, levelTransferPolicy_(*comm_, weights_, prm)//.get<int>("pressure_var_index", 1))
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : pt())
, twolevel_method_(linearoperator,
finesmoother_,
levelTransferPolicy_,
coarseSolverPolicy_,
prm.get<int>("pre_smooth", transpose? 1 : 0),
prm.get<int>("post_smooth", transpose? 0 : 1))
, prm_(prm)
{
if (prm.get<int>("verbosity", 0) > 10 && comm.communicator().rank() == 0) {
auto filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
std::ofstream outfile(filename);
if (!outfile) {
OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
}
Dune::writeMatrixMarket(weights_, outfile);
}
}
virtual void pre(VectorType& x, VectorType& b) override
{
twolevel_method_.pre(x, b);
}
virtual void apply(VectorType& v, const VectorType& d) override
{
twolevel_method_.apply(v, d);
}
virtual void post(VectorType& x) override
{
twolevel_method_.post(x);
}
virtual void update() override
{
weights_ = weightsCalculator_();
updateImpl(comm_);
}
virtual Dune::SolverCategory::Category category() const override
{
return linear_operator_.category();
}
private:
using PressureMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
using ParCoarseOperatorType
= Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Communication>;
using CoarseOperatorType = std::conditional_t<std::is_same<Communication, Dune::Amg::SequentialInformation>::value,
SeqCoarseOperatorType,
ParCoarseOperatorType>;
using LevelTransferPolicy = Opm::PressureBhpTransferPolicy<OperatorType, CoarseOperatorType, Communication, transpose>;
using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperatorType,
FlexibleSolver<CoarseOperatorType>,
LevelTransferPolicy>;
using TwoLevelMethod
= Dune::Amg::TwoLevelMethodCpr<OperatorType, CoarseSolverPolicy, Dune::Preconditioner<VectorType, VectorType>>;
// Handling parallel vs serial instantiation of preconditioner factory.
template <class Comm>
void updateImpl(const Comm*)
{
// Parallel case.
// using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
// auto op_prec = std::make_shared<ParOperatorType>(linear_operator_.getmat(), *comm_);
auto child = prm_.get_child_optional("finesmoother");
finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : pt(), *comm_);
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
// linearoperator_for_precond_ = op_prec;
}
void updateImpl(const Dune::Amg::SequentialInformation*)
{
// Serial case.
// using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
// auto op_prec = std::make_shared<SeqOperatorType>(linear_operator_.getmat());
auto child = prm_.get_child_optional("finesmoother");
finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : pt());
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
// linearoperator_for_precond_ = op_prec;
}
const OperatorType& linear_operator_;
std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> finesmoother_;
const Communication* comm_;
std::function<VectorType()> weightsCalculator_;
VectorType weights_;
LevelTransferPolicy levelTransferPolicy_;
CoarseSolverPolicy coarseSolverPolicy_;
TwoLevelMethod twolevel_method_;
boost::property_tree::ptree prm_;
Communication dummy_comm_;
//std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;
};
} // namespace Dune

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

@ -24,10 +24,12 @@
#include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
#include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
#include <opm/simulators/linalg/OwningTwoLevelPreconditionerWell.hpp>
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <opm/simulators/linalg/amgcpr.hh>
#include <opm/simulators/linalg/WellOperators.hpp>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
@ -328,6 +330,27 @@ private:
}
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, weightsCalculator, pressureIndex, comm);
});
if constexpr (std::is_same_v<O, WellModelGhostLastMatrixAdapter<M, V, V, true>>) {
doAddCreator("cprw",
[](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, std::size_t pressureIndex, const C& comm) {
assert(weightsCalculator);
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<OwningTwoLevelPreconditionerWell<O, V, false, Comm>>(
op, prm, weightsCalculator, comm);
});
doAddCreator("cprwt",
[](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, std::size_t pressureIndex, const C& comm) {
assert(weightsCalculator);
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<OwningTwoLevelPreconditionerWell<O, V, true, Comm>>(
op, prm, weightsCalculator, comm);
});
}
}
// Add a useful default set of preconditioners to the factory.
@ -449,6 +472,21 @@ private:
return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
});
}
if constexpr (std::is_same_v<O, WellModelMatrixAdapter<M, V, V, false>>) {
doAddCreator("cprw", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
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<OwningTwoLevelPreconditionerWell<O, V, false>>(op, prm, weightsCalculator);
});
doAddCreator("cprwt", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
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<OwningTwoLevelPreconditionerWell<O, V, true>>(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())
{

254
opm/simulators/linalg/PressureBhpTransferPolicy.hpp Normal file
View File

@ -0,0 +1,254 @@
/*
Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2019 Dr. Blatt - HPC-Simulation-Software & Services.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include <opm/simulators/linalg/twolevelmethodcpr.hh>
#include <dune/istl/paamg/pinfo.hh>
namespace Opm
{
template <class Communication>
void extendCommunicatorWithWells(const Communication& comm,
std::shared_ptr<Communication>& commRW,
const int nw)
{
// used for extending the coarse communicator pattern
using IndexSet = typename Communication::ParallelIndexSet;
using LocalIndex = typename IndexSet::LocalIndex;
const IndexSet& indset = comm.indexSet();
IndexSet& indset_rw = commRW->indexSet();
const int max_nw = comm.communicator().max(nw);
const int num_proc = comm.communicator().size();
int glo_max = 0;
size_t loc_max = 0;
indset_rw.beginResize();
for (auto ind = indset.begin(), indend = indset.end(); ind != indend; ++ind) {
indset_rw.add(ind->global(), LocalIndex(ind->local(), ind->local().attribute(), true));
const int glo = ind->global();
const size_t loc = ind->local().local();
glo_max = std::max(glo_max, glo);
loc_max = std::max(loc_max, loc);
}
const int global_max = comm.communicator().max(glo_max);
// used append the welldofs at the end
assert(loc_max + 1 == indset.size());
size_t local_ind = loc_max + 1;
for (int i = 0; i < nw; ++i) {
// need to set unique global number
const size_t v = global_max + max_nw * num_proc + i + 1;
// set to true to not have problems with higher levels if growing of domains is used
indset_rw.add(v, LocalIndex(local_ind, Dune::OwnerOverlapCopyAttributeSet::owner, true));
++local_ind;
}
indset_rw.endResize();
assert(indset_rw.size() == indset.size() + nw);
// assume same communication pattern
commRW->remoteIndices().setNeighbours(comm.remoteIndices().getNeighbours());
commRW->remoteIndices().template rebuild<true>();
}
template <class FineOperator, class CoarseOperator, class Communication, bool transpose = false>
class PressureBhpTransferPolicy : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, CoarseOperator>
{
public:
typedef Dune::Amg::LevelTransferPolicyCpr<FineOperator, CoarseOperator> ParentType;
typedef Communication ParallelInformation;
typedef typename FineOperator::domain_type FineVectorType;
public:
PressureBhpTransferPolicy(const Communication& comm,
const FineVectorType& weights,
const boost::property_tree::ptree& prm)
: communication_(&const_cast<Communication&>(comm))
, weights_(weights)
, prm_(prm)
, pressure_var_index_(prm_.get<int>("pressure_var_index"))
{
}
virtual void createCoarseLevelSystem(const FineOperator& fineOperator) override
{
using CoarseMatrix = typename CoarseOperator::matrix_type;
const auto& fineLevelMatrix = fineOperator.getmat();
const auto& nw = fineOperator.getNumberOfExtraEquations();
if (prm_.get<bool>("add_wells")) {
const size_t average_elements_per_row
= static_cast<size_t>(std::ceil(fineLevelMatrix.nonzeroes() / fineLevelMatrix.N()));
const double overflow_fraction = 1.2;
coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N() + nw,
fineLevelMatrix.M() + nw,
average_elements_per_row,
overflow_fraction,
CoarseMatrix::implicit));
int rownum = 0;
for (const auto& row : fineLevelMatrix) {
for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
coarseLevelMatrix_->entry(rownum, col.index()) = 0.0;
}
++rownum;
}
} else {
coarseLevelMatrix_.reset(
new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise));
auto createIter = coarseLevelMatrix_->createbegin();
for (const auto& row : fineLevelMatrix) {
for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
createIter.insert(col.index());
}
++createIter;
}
}
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
if constexpr (std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
coarseLevelCommunication_ = std::make_shared<Communication>();
} else {
coarseLevelCommunication_ = std::make_shared<Communication>(
communication_->communicator(), communication_->category(), false);
}
#else
if constexpr (std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
coarseLevelCommunication_ = std::make_shared<Communication>();
} else {
coarseLevelCommunication_ = std::make_shared<Communication>(
communication_->communicator(), communication_->getSolverCategory(), false);
}
#endif
if (prm_.get<bool>("add_wells")) {
fineOperator.addWellPressureEquationsStruct(*coarseLevelMatrix_);
coarseLevelMatrix_->compress(); // all elemenst should be set
if constexpr (!std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
extendCommunicatorWithWells(*communication_, coarseLevelCommunication_, nw);
}
}
calculateCoarseEntries(fineOperator);
this->lhs_.resize(this->coarseLevelMatrix_->M());
this->rhs_.resize(this->coarseLevelMatrix_->N());
using OperatorArgs = typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_);
this->operator_ = Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs);
#else
OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_);
this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
#endif
}
virtual void calculateCoarseEntries(const FineOperator& fineOperator) override
{
const auto& fineMatrix = fineOperator.getmat();
*coarseLevelMatrix_ = 0;
auto rowCoarse = coarseLevelMatrix_->begin();
for (auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); row != rowEnd; ++row, ++rowCoarse) {
assert(row.index() == rowCoarse.index());
auto entryCoarse = rowCoarse->begin();
for (auto entry = row->begin(), entryEnd = row->end(); entry != entryEnd; ++entry, ++entryCoarse) {
assert(entry.index() == entryCoarse.index());
double matrix_el = 0;
if (transpose) {
const auto& bw = weights_[entry.index()];
for (size_t i = 0; i < bw.size(); ++i) {
matrix_el += (*entry)[pressure_var_index_][i] * bw[i];
}
} else {
const auto& bw = weights_[row.index()];
for (size_t i = 0; i < bw.size(); ++i) {
matrix_el += (*entry)[i][pressure_var_index_] * bw[i];
}
}
(*entryCoarse) = matrix_el;
}
}
if (prm_.get<bool>("add_wells")) {
assert(transpose == false); // not implemented
fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_);
}
#ifndef NDEBUG
std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
assert(rowCoarse == coarseLevelMatrix_->end());
#endif
}
virtual void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override
{
//NB we iterate over fine assumming welldofs is at the end
// Set coarse vector to zero
this->rhs_ = 0;
auto end = fine.end(), begin = fine.begin();
for (auto block = begin; block != end; ++block) {
const auto& bw = weights_[block.index()];
double rhs_el = 0.0;
if (transpose) {
rhs_el = (*block)[pressure_var_index_];
} else {
for (size_t i = 0; i < block->size(); ++i) {
rhs_el += (*block)[i] * bw[i];
}
}
this->rhs_[block - begin] = rhs_el;
}
this->lhs_ = 0;
}
virtual void moveToFineLevel(typename ParentType::FineDomainType& fine) override
{
//NB we iterate over fine assumming welldofs is at the end
auto end = fine.end(), begin = fine.begin();
for (auto block = begin; block != end; ++block) {
if (transpose) {
const auto& bw = weights_[block.index()];
for (size_t i = 0; i < block->size(); ++i) {
(*block)[i] = this->lhs_[block - begin] * bw[i];
}
} else {
(*block)[pressure_var_index_] = this->lhs_[block - begin];
}
}
}
virtual PressureBhpTransferPolicy* clone() const override
{
return new PressureBhpTransferPolicy(*this);
}
const Communication& getCoarseLevelCommunication() const
{
return *coarseLevelCommunication_;
}
private:
Communication* communication_;
const FineVectorType& weights_;
boost::property_tree::ptree prm_;
const int pressure_var_index_;
std::shared_ptr<Communication> coarseLevelCommunication_;
std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
};
} // namespace Opm

79
opm/simulators/linalg/WellOperators.hpp
View File

@ -45,30 +45,40 @@ namespace Opm
/// WellModelMatrixAdapter separately for each TypeTag, as it will only
/// depend on the matrix and vector types involved, which typically are
/// just one for each block size with block sizes 1-4.
template <class WellModel, class X, class Y>
class WellModelAsLinearOperator : public Dune::LinearOperator<X, Y>
template <class X, class Y>
class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
{
public:
using Base = Dune::LinearOperator<X, Y>;
using field_type = typename Base::field_type;
using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
virtual void addWellPressureEquations(PressureMatrix& jacobian, const X& weights) const = 0;
virtual void addWellPressureEquationsStruct(PressureMatrix& jacobian) const = 0;
virtual int getNumberOfExtraEquations() const = 0;
};
template <class WellModel, class X, class Y>
class WellModelAsLinearOperator : public Opm::LinearOperatorExtra<X, Y>
{
public:
using Base = Opm::LinearOperatorExtra<X, Y>;
using field_type = typename Base::field_type;
using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
explicit WellModelAsLinearOperator(const WellModel& wm)
: wellMod_(wm)
{
}
/*! \brief apply operator to x: \f$ y = A(x) \f$
The input vector is consistent and the output must also be
The input vector is consistent and the output must also be
consistent on the interior+border partition.
*/
void apply (const X& x, Y& y) const override
void apply(const X& x, Y& y) const override
{
wellMod_.apply(x, y );
wellMod_.apply(x, y);
}
//! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
virtual void applyscaleadd(field_type alpha, const X& x, Y& y) const override
{
wellMod_.applyScaleAdd( alpha, x, y );
wellMod_.applyScaleAdd(alpha, x, y);
}
/// Category for operator.
@ -80,6 +90,19 @@ public:
{
return Dune::SolverCategory::sequential;
}
void addWellPressureEquations(PressureMatrix& jacobian, const X& weights) const override
{
wellMod_.addWellPressureEquations(jacobian, weights);
}
void addWellPressureEquationsStruct(PressureMatrix& jacobian) const override
{
wellMod_.addWellPressureEquationsStruct(jacobian);
}
int getNumberOfExtraEquations() const override
{
return wellMod_.numLocalWells();
}
private:
const WellModel& wellMod_;
};
@ -103,7 +126,7 @@ public:
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
#else
@ -118,7 +141,7 @@ public:
//! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A,
const Dune::LinearOperator<X, Y>& wellOper,
const Opm::LinearOperatorExtra<X, Y>& wellOper,
const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >())
: A_( A ), wellOper_( wellOper ), comm_(comm)
{}
@ -153,9 +176,22 @@ public:
virtual const matrix_type& getmat() const override { return A_; }
void addWellPressureEquations(PressureMatrix& jacobian, const X& weights) const
{
wellOper_.addWellPressureEquations(jacobian, weights);
}
void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
{
wellOper_.addWellPressureEquationsStruct(jacobian);
}
int getNumberOfExtraEquations() const
{
return wellOper_.getNumberOfExtraEquations();
}
protected:
const matrix_type& A_ ;
const Dune::LinearOperator<X, Y>& wellOper_;
const Opm::LinearOperatorExtra<X, Y>& wellOper_;
std::shared_ptr< communication_type > comm_;
};
@ -176,7 +212,7 @@ public:
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
#else
@ -192,7 +228,7 @@ public:
//! constructor: just store a reference to a matrix
WellModelGhostLastMatrixAdapter (const M& A,
const Dune::LinearOperator<X, Y>& wellOper,
const Opm::LinearOperatorExtra<X, Y>& wellOper,
const size_t interiorSize )
: A_( A ), wellOper_( wellOper ), interiorSize_(interiorSize)
{}
@ -230,6 +266,19 @@ public:
virtual const matrix_type& getmat() const override { return A_; }
void addWellPressureEquations(PressureMatrix& jacobian, const X& weights) const
{
wellOper_.addWellPressureEquations(jacobian, weights);
}
void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
{
wellOper_.addWellPressureEquationsStruct(jacobian);
}
int getNumberOfExtraEquations() const
{
return wellOper_.getNumberOfExtraEquations();
}
protected:
void ghostLastProject(Y& y) const
{
@ -239,7 +288,7 @@ protected:
}
const matrix_type& A_ ;
const Dune::LinearOperator<X, Y>& wellOper_;
const Opm::LinearOperatorExtra< X, Y>& wellOper_;
size_t interiorSize_;
};

30
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -294,6 +294,29 @@ namespace Opm {
WellInterfacePtr getWell(const std::string& well_name) const;
bool hasWell(const std::string& well_name) const;
// The number of components in the model.
int numComponents() const;
int numLocalWells() const;
int numPhases() const;
using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights) const
{
for (const auto& well : well_container_) {
well->addWellPressureEquations(jacobian, weights);
}
}
void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
{
for (const auto& well : well_container_) {
well->addWellPressureEquationsStruct(jacobian);
}
}
void initGliftEclWellMap(GLiftEclWells &ecl_well_map);
/// \brief Get list of local nonshut wells
@ -386,12 +409,7 @@ namespace Opm {
void calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
DeferredLogger& deferred_logger);
// The number of components in the model.
int numComponents() const;
int reportStepIndex() const;
void assembleWellEq(const double dt, DeferredLogger& deferred_logger);
void assembleWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger);
bool maybeDoGasLiftOptimize(DeferredLogger& deferred_logger);

7
opm/simulators/wells/StandardWell.hpp
View File

@ -94,6 +94,7 @@ namespace Opm
using PolymerModule = BlackOilPolymerModule<TypeTag>;
using FoamModule = BlackOilFoamModule<TypeTag>;
using BrineModule = BlackOilBrineModule<TypeTag>;
using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
// number of the conservation equations
static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
@ -178,7 +179,11 @@ namespace Opm
WellState& well_state,
DeferredLogger& deferred_logger) const override;
virtual void addWellContributions(SparseMatrixAdapter& mat) const override;
virtual void addWellContributions(SparseMatrixAdapter& mat) const override;
virtual void addWellPressureEquationsStruct(PressureMatrix& mat) const override;
virtual void addWellPressureEquations(PressureMatrix& mat, const BVector& x) const override;
// iterate well equations with the specified control until converged
bool iterateWellEqWithControl(const Simulator& ebosSimulator,

90
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -2170,6 +2170,96 @@ namespace Opm
template <typename TypeTag>
void
StandardWell<TypeTag>::addWellPressureEquationsStruct(PressureMatrix& jacobian) const
{
// sustem is the pressur variant of
// We need to change matrx A as follows
// A CT
// B D
// we need to add the elemenst of CT
// then we need to ad the quasiimpes type well equation for B D if the well is not
// BHP contolled
const int welldof_ind = duneC_.M() + index_of_well_;
for (auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC) {
const auto row_index = colC.index();
double matel = 0;
jacobian.entry(row_index, welldof_ind) = matel;
}
jacobian.entry(welldof_ind, welldof_ind) = 0.0;
// set the matrix elements for well reservoir coupling
for (auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
double matel = 0;
jacobian.entry(welldof_ind, col_index) = matel;
}
}
template <typename TypeTag>
void
StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights) const
{
// sustem is the pressur variant of
// We need to change matrx A as follows
// A CT
// B D
// we need to add the elemenst of CT
// then we need to ad the quasiimpes type well equation for B D if the well is not
// BHP contolled
int bhp_var_index = Bhp;
assert(duneC_.M() == weights.size());
const int welldof_ind = duneC_.M() + index_of_well_;
for (auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC) {
const auto row_ind = colC.index();
const auto& bw = weights[row_ind];
double matel = 0;
assert((*colC).M() == bw.size());
for (size_t i = 0; i < bw.size(); ++i) {
matel += (*colC)[bhp_var_index][i] * bw[i];
}
jacobian[row_ind][welldof_ind] = matel;
}
// make quasipes weights for bhp it should be trival
using VectorBlockType = typename BVector::block_type;
using MatrixBlockType = DiagMatrixBlockWellType;
VectorBlockType bweights(0.0);
double diagElem = 0;
{
// const DiagMatrixBlockWellType& invA = invDuneD_[0][0];
VectorBlockType rhs(0.0);
rhs[bhp_var_index] = 1.0;
MatrixBlockType inv_diag_block = invDuneD_[0][0];
auto inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
// diag_block_transpose.solve(bweights, rhs);
inv_diag_block_transpose.mv(rhs, bweights);
// NB how to scale to make it most symmetric
// double abs_max = *std::max_element(
// bweights.begin(), bweights.end(), [](double a, double b) { return std::fabs(a) < std::fabs(b); });
// bweights /= std::fabs(abs_max);
}
//
jacobian[welldof_ind][welldof_ind] = 1.0;
// set the matrix elements for well reservoir coupling
for (auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
const auto& bw = bweights;
double matel = 0;
for (size_t i = 0; i < bw.size(); ++i) {
matel += (*colB)[bhp_var_index][i] * bw[i];
}
jacobian[welldof_ind][col_index] = matel;
}
}
template<typename TypeTag>
typename StandardWell<TypeTag>::EvalWell
StandardWell<TypeTag>::

13
opm/simulators/wells/WellHelpers.hpp
View File

@ -233,8 +233,19 @@ namespace Opm {
}
} // namespace wellhelpers
template <class DenseMatrix>
DenseMatrix transposeDenseDynMatrix(const DenseMatrix& M)
{
DenseMatrix tmp{M.cols(), M.rows()};
for (int i = 0; i < M.rows(); ++i) {
for (int j = 0; j < M.cols(); ++j) {
tmp[j][i] = M[i][j];
}
}
return tmp;
}
} // namespace wellhelpers
}
#endif

10
opm/simulators/wells/WellInterface.hpp
View File

@ -95,6 +95,7 @@ public:
using MatrixBlockType = Dune::FieldMatrix<Scalar, Indices::numEq, Indices::numEq>;
using BVector = Dune::BlockVector<VectorBlockType>;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
using PressureMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
using RateConverterType =
typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType;
@ -224,6 +225,15 @@ public:
// Add well contributions to matrix
virtual void addWellContributions(SparseMatrixAdapter&) const = 0;
virtual void addWellPressureEquationsStruct(PressureMatrix&) const
{
}
virtual void addWellPressureEquations(PressureMatrix&,
const BVector& x) const
{
}
void addCellRates(RateVector& rates, int cellIdx) const;
Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const;

2
tests/test_flexiblesolver.cpp
View File

@ -74,7 +74,7 @@ testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, con
using SeqOperatorType = Dune::MatrixAdapter<Matrix, Vector, Vector>;
SeqOperatorType op(matrix);
Dune::FlexibleSolver<Matrix, Vector> solver(op, prm, wc, 1);
Dune::FlexibleSolver<SeqOperatorType> solver(op, prm, wc, 1);
Vector x(rhs.size());
Dune::InverseOperatorResult res;
solver.apply(x, rhs, res);