Add well operator for NLDD domains

opm/simulators/flow/BlackoilModelNldd.hpp

@@ -169,6 +169,7 @@ public:
             loc_param.linear_solver_print_json_definition_ = false;
             const bool force_serial = true;
             domain_linsolvers_.emplace_back(model_.simulator(), loc_param, force_serial);
+            domain_linsolvers_.back().setDomainIndex(index);
         }
 
         assert(int(domains_.size()) == num_domains);

opm/simulators/linalg/ISTLSolver.hpp

@@ -352,6 +352,8 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
                 // Outch! We need to be able to scale the linear system! Hence const_cast
                 matrix_ = const_cast<Matrix*>(&M);
 
+                is_nldd_local_solver_ = parameters_[activeSolverNum_].is_nldd_local_solver_;
+
                 useWellConn_ = Parameters::Get<Parameters::MatrixAddWellContributions>();
                 // setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
             } else {
@@ -450,6 +452,11 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
 
         const CommunicationType* comm() const { return comm_.get(); }
 
+        void setDomainIndex(int index)
+        {
+            domainIndex_ = index;
+        }
+
     protected:
 #if HAVE_MPI
         using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
@@ -490,8 +497,15 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
             OPM_TIMEBLOCK(flexibleSolverPrepare);
             if (shouldCreateSolver()) {
                 if (!useWellConn_) {
-                    auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
-                    flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
+                    if (is_nldd_local_solver_) {
+                        auto wellOp = std::make_unique<DomainWellModelAsLinearOperator<WellModel, Vector, Vector>>(simulator_.problem().wellModel());
+                        wellOp->setDomainIndex(domainIndex_);
+                        flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
+                    }
+                    else {
+                        auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
+                        flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
+                    }
                 }
                 std::function<Vector()> weightCalculator = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex);
                 OPM_TIMEBLOCK(flexibleSolverCreate);
@@ -651,6 +665,10 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
         std::vector<int> overlapRows_;
         std::vector<int> interiorRows_;
 
+        bool is_nldd_local_solver_;
+
+        int domainIndex_ = -1;
+
         bool useWellConn_;
 
         std::vector<FlowLinearSolverParameters> parameters_;

opm/simulators/linalg/WellOperators.hpp

@@ -123,10 +123,46 @@ public:
         return wellMod_.numLocalWellsEnd();
     }
 
-private:
+protected:
     const WellModel& wellMod_;
 };
 
+
+template <class WellModel, class X, class Y>
+class DomainWellModelAsLinearOperator : public WellModelAsLinearOperator<WellModel, X, Y>
+{
+public:
+    using WBase = WellModelAsLinearOperator<WellModel, X, Y>;
+    using WBase::WBase;  // inherit all constructors from the base class
+    using field_type = typename WBase::field_type;
+    using PressureMatrix = typename WBase::PressureMatrix;
+
+    void setDomainIndex(int index) { domainIndex_ = index; }
+
+    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_);
+    }
+
+    void addWellPressureEquations(PressureMatrix& jacobian,
+                                  const X& weights,
+                                  const bool use_well_weights) const override
+    {
+        OPM_TIMEBLOCK(addWellPressureEquations);
+        this->wellMod_.addWellPressureEquationsDomain(jacobian, weights, use_well_weights, domainIndex_);
+    }
+
+private:
+    int domainIndex_ = -1;
+};
+
 /*!
    \brief Adapter to combine a matrix and another linear operator into
    a combined linear operator.

opm/simulators/wells/BlackoilWellModel.hpp

@@ -300,6 +300,8 @@ template<class Scalar> class WellContributions;
             // subtract B*inv(D)*C * x from A*x
             void apply(const BVector& x, BVector& Ax) const;
 
+            void applyDomain(const BVector& x, BVector& Ax, int domainIndex) const;
+
 #if COMPILE_BDA_BRIDGE
             // accumulate the contributions of all Wells in the WellContributions object
             void getWellContributions(WellContributions<Scalar>& x) const;
@@ -308,6 +310,8 @@ template<class Scalar> class WellContributions;
             // 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, int domainIndex) const;
+
             // Check if well equations is converged.
             ConvergenceReport getWellConvergence(const std::vector<Scalar>& B_avg, const bool checkWellGroupControls = false) const;
 
@@ -358,6 +362,9 @@ template<class Scalar> class WellContributions;
 
             void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const;
 
+            void addWellPressureEquationsDomain(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights, int domainIndex) const;
+
+
             void addWellPressureEquationsStruct(PressureMatrix& jacobian) const;
 
             void initGliftEclWellMap(GLiftEclWells &ecl_well_map);

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1749,9 +1749,36 @@ namespace Opm {
             x_local.resize(cells.size());
             Ax_local.resize(cells.size());
 
-            if (this->param_.nonlinear_solver_ == "nldd") {
-                // transfer global cells index to local subdomain cells index
+            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, int domainIndex) const
+    {
+        BVector x_local;
+        BVector Ax_local;
+
+        for (auto& well : well_container_) {
+            if (well_domain_.at(well->name()) == domainIndex) {
+                // Well equations B and C uses only the perforated cells, so need to apply on local vectors
+                auto cells = well->cells();
+                x_local.resize(cells.size());
+                Ax_local.resize(cells.size());
+                // transfer global cells index to local subdomain cells index
                 auto domain_cells = domains_cells_[well_domain_.at(well->name())];
 
                 // Assuming domain_cells is sorted
@@ -1765,18 +1792,18 @@ namespace Opm {
                         std::cerr << "Cell value " << cells[i] << " not found in domain_cells." << std::endl;
                     }
                 }
-            }
 
-            for (size_t i = 0; i < cells.size(); ++i) {
-                x_local[i] = x[cells[i]];
-                Ax_local[i] = Ax[cells[i]];
-            }
+                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);
+                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];
+                for (size_t i = 0; i < cells.size(); ++i) {
+                    // only need to update Ax
+                    Ax[cells[i]] = Ax_local[i];
+                }
             }
         }
     }
@@ -1840,6 +1867,27 @@ namespace Opm {
         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, 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>::
@@ -1867,6 +1915,28 @@ namespace Opm {
         }
     }
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    addWellPressureEquationsDomain(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights, int domainIndex) const
+    {
+        throw std::logic_error("CPRW is not yet implemented for NLDD subdomains");
+        // To fix this function, rdofs should be the size of the domain, and the nw should be the number of wells in the domain
+        // int nw =  this->numLocalWellsEnd(); // should number of wells in the domain
+        // int rdofs = local_num_cells_;  // should be the size of the domain
+        // for ( int i = 0; i < nw; i++ ){
+        //     int wdof = rdofs + i;
+        //     jacobian[wdof][wdof] = 1.0;// better scaling ?
+        // }
+
+        // for ( const auto& well : well_container_ ) {
+        //     if (well_domain_.at(well->name()) == domainIndex) {
+                    // weights should be the size of the domain
+        //         well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState());
+        //     }
+        // }
+    }
+
     template <typename TypeTag>
     void BlackoilWellModel<TypeTag>::
     addReservoirSourceTerms(GlobalEqVector& residual,