changed: move implementation of WellOperators out of BlackoilWellModel

opm/simulators/linalg/WellOperators.hpp

@@ -84,14 +84,43 @@ public:
     void apply(const X& x, Y& y) const override
     {
         OPM_TIMEBLOCK(apply);
-        wellMod_.apply(x, y);
+        for (const auto& well : this->wellMod_) {
+            // 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] = 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];
+            }
+        }
     }
 
     //! 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 +154,14 @@ public:
 
 protected:
     const WellModel& wellMod_;
-};
 
+    // 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 +177,47 @@ public:
     void apply(const X& x, Y& y) const override
     {
         OPM_TIMEBLOCK(apply);
-        this->wellMod_.applyDomain(x, y, domainIndex_);
+        std::size_t well_index = 0;
+        for (const auto& well : this->wellMod_) {
+            if (this->wellMod_.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 = this->wellMod_.well_local_cells()[well_index];
+                this->x_local_.resize(local_cells.size());
+                this->Ax_local_.resize(local_cells.size());
+
+                for (size_t i = 0; i < local_cells.size(); ++i) {
+                    this->x_local_[i] = x[local_cells[i]];
+                    this->Ax_local_[i] = y[local_cells[i]];
+                }
+
+                well->apply(this->x_local_, this->Ax_local_);
+
+                for (size_t i = 0; i < local_cells.size(); ++i) {
+                    // only need to update Ax
+                    y[local_cells[i]] = this->Ax_local_[i];
+                }
+            }
+            ++well_index;
+        }
     }
 
     void applyscaleadd(field_type alpha, const X& x, Y& y) const override
     {
         OPM_TIMEBLOCK(applyscaleadd);
-        this->wellMod_.applyScaleAddDomain(alpha, x, y, domainIndex_);
+        if (this->wellMod_.empty()) {
+            return;
+        }
+
+        if (this->scaleAddRes_.size() != y.size()) {
+            this->scaleAddRes_.resize(y.size());
+        }
+
+        this->scaleAddRes_ = 0.0;
+        // scaleAddRes_  = - C D^-1 B x
+        this->apply(x, this->scaleAddRes_);
+        // Ax = Ax + alpha * scaleAddRes_
+        y.axpy(alpha, this->scaleAddRes_);
     }
 
     void addWellPressureEquations(PressureMatrix& jacobian,

opm/simulators/wells/BlackoilWellModel.hpp

@@ -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;
 
@@ -439,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
@@ -583,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_;
         };

opm/simulators/wells/BlackoilWellModel_impl.hpp

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