Merge pull request #5823 from akva2/move_wellop_impl

Move implementation of WellOperators out of BlackoilWellModel
This commit is contained in:
Arne Morten Kvarving
2024-12-20 12:35:21 +01:00
committed by GitHub
4 changed files with 70 additions and 125 deletions

View File

@@ -84,14 +84,28 @@ public:
void apply(const X& x, Y& y) const override
{
OPM_TIMEBLOCK(apply);
wellMod_.apply(x, y);
for (const auto& well : this->wellMod_) {
this->applySingleWell(x, y, well, well->cells());
}
}
//! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
void applyscaleadd(field_type alpha, const X& x, Y& y) const override
{
OPM_TIMEBLOCK(applyscaleadd);
wellMod_.applyScaleAdd(alpha, x, y);
if (this->wellMod_.empty()) {
return;
}
if (scaleAddRes_.size() != y.size()) {
scaleAddRes_.resize(y.size());
}
scaleAddRes_ = 0.0;
// scaleAddRes_ = - C D^-1 B x
apply(x, scaleAddRes_);
// Ax = Ax + alpha * scaleAddRes_
y.axpy(alpha, scaleAddRes_);
}
/// Category for operator.
@@ -125,8 +139,37 @@ public:
protected:
const WellModel& wellMod_;
};
template<class WellType, class ArrayType>
void applySingleWell(const X& x, Y& y,
const WellType& well,
const ArrayType& cells) const
{
// Well equations B and C uses only the perforated cells, so need to apply on local vectors
x_local_.resize(cells.size());
Ax_local_.resize(cells.size());
for (size_t i = 0; i < cells.size(); ++i) {
x_local_[i] = x[cells[i]];
Ax_local_[i] = y[cells[i]];
}
well->apply(x_local_, Ax_local_);
for (size_t i = 0; i < cells.size(); ++i) {
// only need to update Ax
y[cells[i]] = Ax_local_[i];
}
}
// These members are used to avoid reallocation.
// Their state is not relevant between function calls, so they can
// (and must) be mutable, as the functions using them are const.
mutable X x_local_{};
mutable Y Ax_local_{};
mutable Y scaleAddRes_{};
};
template <class WellModel, class X, class Y>
class DomainWellModelAsLinearOperator : public WellModelAsLinearOperator<WellModel, X, Y>
@@ -142,13 +185,13 @@ public:
void apply(const X& x, Y& y) const override
{
OPM_TIMEBLOCK(apply);
this->wellMod_.applyDomain(x, y, domainIndex_);
}
void applyscaleadd(field_type alpha, const X& x, Y& y) const override
{
OPM_TIMEBLOCK(applyscaleadd);
this->wellMod_.applyScaleAddDomain(alpha, x, y, domainIndex_);
std::size_t well_index = 0;
for (const auto& well : this->wellMod_) {
if (this->wellMod_.well_domain().at(well->name()) == domainIndex_) {
this->applySingleWell(x, y, well, this->wellMod_.well_local_cells()[well_index]);
}
++well_index;
}
}
void addWellPressureEquations(PressureMatrix& jacobian,

View File

@@ -299,21 +299,11 @@ template<class Scalar> class WellContributions;
return this->computeWellBlockAveragePressures(this->gravity_);
}
// subtract B*inv(D)*C * x from A*x
void apply(const BVector& x, BVector& Ax) const;
void applyDomain(const BVector& x, BVector& Ax, const int domainIndex) const;
#if COMPILE_GPU_BRIDGE
// accumulate the contributions of all Wells in the WellContributions object
void getWellContributions(WellContributions<Scalar>& x) const;
#endif
// apply well model with scaling of alpha
void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
void applyScaleAddDomain(const Scalar alpha, const BVector& x, BVector& Ax, const int domainIndex) const;
// Check if well equations is converged.
ConvergenceReport getWellConvergence(const std::vector<Scalar>& B_avg, const bool checkWellGroupControls = false) const;
@@ -389,6 +379,14 @@ template<class Scalar> class WellContributions;
void setupDomains(const std::vector<Domain>& domains);
const SparseTable<int>& well_local_cells() const
{
return well_local_cells_;
}
auto begin() const { return well_container_.begin(); }
auto end() const { return well_container_.end(); }
bool empty() const { return well_container_.empty(); }
protected:
Simulator& simulator_;
@@ -411,8 +409,9 @@ template<class Scalar> class WellContributions;
createTypedWellPointer(const int wellID,
const int time_step) const;
WellInterfacePtr createWellForWellTest(const std::string& well_name, const int report_step, DeferredLogger& deferred_logger) const;
WellInterfacePtr createWellForWellTest(const std::string& well_name,
const int report_step,
DeferredLogger& deferred_logger) const;
const ModelParameters param_;
std::size_t global_num_cells_{};
@@ -430,9 +429,6 @@ template<class Scalar> class WellContributions;
// Pre-step network solve at static reservoir conditions (group and well states might be updated)
void doPreStepNetworkRebalance(DeferredLogger& deferred_logger);
// used to better efficiency of calcuation
mutable BVector scaleAddRes_{};
std::vector<Scalar> B_avg_{};
// Store the local index of the wells perforated cells in the domain, if using subdomains
@@ -574,11 +570,10 @@ template<class Scalar> class WellContributions;
BlackoilWellModel(Simulator& simulator, const PhaseUsage& pu);
// These members are used to avoid reallocation in specific functions
// (e.g., apply, applyDomain) instead of using local variables.
// instead of using local variables.
// Their state is not relevant between function calls, so they can
// (and must) be mutable, as the functions using them are const.
mutable BVector x_local_;
mutable BVector Ax_local_;
mutable BVector res_local_;
mutable GlobalEqVector linearize_res_local_;
};

View File

@@ -208,6 +208,11 @@ public:
const std::map<std::string, double>& wellOpenTimes() const { return well_open_times_; }
const std::map<std::string, double>& wellCloseTimes() const { return well_close_times_; }
const std::map<std::string, int>& well_domain() const
{
return well_domain_;
}
bool reportStepStarts() const { return report_step_starts_; }
bool shouldBalanceNetwork(const int reportStepIndex,

View File

@@ -1850,62 +1850,6 @@ namespace Opm {
}
// Ax = A x - C D^-1 B x
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
apply(const BVector& x, BVector& Ax) const
{
for (auto& well : well_container_) {
// Well equations B and C uses only the perforated cells, so need to apply on local vectors
const auto& cells = well->cells();
x_local_.resize(cells.size());
Ax_local_.resize(cells.size());
for (size_t i = 0; i < cells.size(); ++i) {
x_local_[i] = x[cells[i]];
Ax_local_[i] = Ax[cells[i]];
}
well->apply(x_local_, Ax_local_);
for (size_t i = 0; i < cells.size(); ++i) {
// only need to update Ax
Ax[cells[i]] = Ax_local_[i];
}
}
}
// Ax = A x - C D^-1 B x
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
applyDomain(const BVector& x, BVector& Ax, const int domainIndex) const
{
for (size_t well_index = 0; well_index < well_container_.size(); ++well_index) {
auto& well = well_container_[well_index];
if (this->well_domain_.at(well->name()) == domainIndex) {
// Well equations B and C uses only the perforated cells, so need to apply on local vectors
// transfer global cells index to local subdomain cells index
const auto& local_cells = well_local_cells_[well_index];
x_local_.resize(local_cells.size());
Ax_local_.resize(local_cells.size());
for (size_t i = 0; i < local_cells.size(); ++i) {
x_local_[i] = x[local_cells[i]];
Ax_local_[i] = Ax[local_cells[i]];
}
well->apply(x_local_, Ax_local_);
for (size_t i = 0; i < local_cells.size(); ++i) {
// only need to update Ax
Ax[local_cells[i]] = Ax_local_[i];
}
}
}
}
#if COMPILE_GPU_BRIDGE
template<typename TypeTag>
void
@@ -1944,48 +1888,6 @@ namespace Opm {
}
#endif
// Ax = Ax - alpha * C D^-1 B x
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const
{
if (this->well_container_.empty()) {
return;
}
if( scaleAddRes_.size() != Ax.size() ) {
scaleAddRes_.resize( Ax.size() );
}
scaleAddRes_ = 0.0;
// scaleAddRes_ = - C D^-1 B x
apply( x, scaleAddRes_ );
// Ax = Ax + alpha * scaleAddRes_
Ax.axpy( alpha, scaleAddRes_ );
}
// Ax = Ax - alpha * C D^-1 B x
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
applyScaleAddDomain(const Scalar alpha, const BVector& x, BVector& Ax, const int domainIndex) const
{
if (this->well_container_.empty()) {
return;
}
if( scaleAddRes_.size() != Ax.size() ) {
scaleAddRes_.resize( Ax.size() );
}
scaleAddRes_ = 0.0;
// scaleAddRes_ = - C D^-1 B x
applyDomain(x, scaleAddRes_, domainIndex);
// Ax = Ax + alpha * scaleAddRes_
Ax.axpy( alpha, scaleAddRes_ );
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::