Updated for upstream changes.

This commit is contained in:
Atgeirr Flø Rasmussen
2022-06-08 17:02:02 +02:00
parent f3acfcde0b
commit fef06a77af
6 changed files with 43 additions and 43 deletions
@@ -80,22 +80,23 @@ template <class OperatorType,
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)
OwningTwoLevelPreconditionerWell(const OperatorType& linearoperator,
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
const std::size_t pressureIndex)
: linear_operator_(linearoperator)
, finesmoother_(PrecFactory::create(linearoperator,
prm.get_child_optional("finesmoother")?
prm.get_child("finesmoother"): pt()))
prm.get_child("finesmoother") : Opm::PropertyTree()))
, 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())
, levelTransferPolicy_(dummy_comm_, weights_, prm, pressureIndex)
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver") ? prm.get_child("coarsesolver") : Opm::PropertyTree())
, twolevel_method_(linearoperator,
finesmoother_,
levelTransferPolicy_,
@@ -115,17 +116,19 @@ public:
}
OwningTwoLevelPreconditionerWell(const OperatorType& linearoperator,
const pt& prm,
const std::function<VectorType()> weightsCalculator, const Communication& comm)
const Opm::PropertyTree& prm,
const std::function<VectorType()> weightsCalculator,
const std::size_t pressureIndex,
const Communication& comm)
: linear_operator_(linearoperator)
, finesmoother_(PrecFactory::create(linearoperator,
prm.get_child_optional("finesmoother")?
prm.get_child("finesmoother"): pt(), comm))
prm.get_child("finesmoother") : Opm::PropertyTree(), 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())
, levelTransferPolicy_(*comm_, weights_, prm, pressureIndex)
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver") ? prm.get_child("coarsesolver") : Opm::PropertyTree())
, twolevel_method_(linearoperator,
finesmoother_,
levelTransferPolicy_,
@@ -194,7 +197,7 @@ private:
// 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_);
finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : Opm::PropertyTree(), *comm_);
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
// linearoperator_for_precond_ = op_prec;
}
@@ -205,7 +208,7 @@ private:
// 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());
finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : Opm::PropertyTree());
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
// linearoperator_for_precond_ = op_prec;
}
@@ -218,7 +221,7 @@ private:
LevelTransferPolicy levelTransferPolicy_;
CoarseSolverPolicy coarseSolverPolicy_;
TwoLevelMethod twolevel_method_;
boost::property_tree::ptree prm_;
Opm::PropertyTree prm_;
Communication dummy_comm_;
//std::shared_ptr<AbstractOperatorType> linearoperator_for_precond_;
};
@@ -339,7 +339,7 @@ private:
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);
op, prm, weightsCalculator, pressureIndex, comm);
});
doAddCreator("cprwt",
[](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, std::size_t pressureIndex, const C& comm) {
@@ -348,7 +348,7 @@ private:
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);
op, prm, weightsCalculator, pressureIndex, comm);
});
}
}
@@ -473,17 +473,17 @@ private:
});
}
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) {
doAddCreator("cprw", [](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<OwningTwoLevelPreconditionerWell<O, V, false>>(op, prm, weightsCalculator);
return std::make_shared<OwningTwoLevelPreconditionerWell<O, V, false>>(op, prm, weightsCalculator, pressureIndex);
});
doAddCreator("cprwt", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
doAddCreator("cprwt", [](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<OwningTwoLevelPreconditionerWell<O, V, true>>(op, prm, weightsCalculator);
return std::make_shared<OwningTwoLevelPreconditionerWell<O, V, true>>(op, prm, weightsCalculator, pressureIndex);
});
}
@@ -78,11 +78,12 @@ namespace Opm
public:
PressureBhpTransferPolicy(const Communication& comm,
const FineVectorType& weights,
const boost::property_tree::ptree& prm)
const Opm::PropertyTree& prm,
const std::size_t pressureIndex)
: communication_(&const_cast<Communication&>(comm))
, weights_(weights)
, prm_(prm)
, pressure_var_index_(prm_.get<int>("pressure_var_index"))
, pressure_var_index_(pressureIndex)
{
}
@@ -243,11 +244,10 @@ namespace Opm
private:
Communication* communication_;
const FineVectorType& weights_;
boost::property_tree::ptree prm_;
PropertyTree prm_;
const int pressure_var_index_;
std::shared_ptr<Communication> coarseLevelCommunication_;
std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
};
} // namespace Opm
+6 -8
View File
@@ -294,13 +294,6 @@ 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
@@ -409,7 +402,12 @@ namespace Opm {
void calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
DeferredLogger& deferred_logger);
void assembleWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger);
// The number of components in the model.
int numComponents() const;
int reportStepIndex() const;
void assembleWellEq(const double dt, Opm::DeferredLogger& deferred_logger);
bool maybeDoGasLiftOptimize(DeferredLogger& deferred_logger);
+8 -9
View File
@@ -2181,8 +2181,8 @@ namespace Opm
// 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 int welldof_ind = this->duneC_.M() + this->index_of_well_;
for (auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
const auto row_index = colC.index();
double matel = 0;
jacobian.entry(row_index, welldof_ind) = matel;
@@ -2191,7 +2191,7 @@ namespace Opm
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) {
for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
double matel = 0;
jacobian.entry(welldof_ind, col_index) = matel;
@@ -2213,9 +2213,9 @@ namespace Opm
// 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) {
assert(this->duneC_.M() == weights.size());
const int welldof_ind = this->duneC_.M() + this->index_of_well_;
for (auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
const auto row_ind = colC.index();
const auto& bw = weights[row_ind];
double matel = 0;
@@ -2227,14 +2227,13 @@ namespace Opm
}
// 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 = this->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);
@@ -2246,7 +2245,7 @@ namespace Opm
//
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) {
for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
const auto& bw = bweights;
double matel = 0;