From 6011a42246fc17f5314a5b6126d5b17c76095dbc Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Fri, 18 Nov 2022 12:09:43 +0100
Subject: [PATCH] added: MultisegmentWellAssemble::assemblePerforationEq

extracted from MultisegmentWell::assembleWellEqWithoutIteration
---
 .../wells/MultisegmentWellAssemble.cpp        | 26 +++++++++++++++++++
 .../wells/MultisegmentWellAssemble.hpp        |  7 +++++
 .../wells/MultisegmentWell_impl.hpp           | 19 ++------------
 3 files changed, 35 insertions(+), 17 deletions(-)

diff --git a/opm/simulators/wells/MultisegmentWellAssemble.cpp b/opm/simulators/wells/MultisegmentWellAssemble.cpp
index f6a7c7254..7ff23370e 100644
--- a/opm/simulators/wells/MultisegmentWellAssemble.cpp
+++ b/opm/simulators/wells/MultisegmentWellAssemble.cpp
@@ -262,6 +262,32 @@ assembleInflowTerm(const int seg,
     // pressure derivative should be zero
 }
 
+template<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assemblePerforationEq(const int seg,
+                      const int cell_idx,
+                      const int comp_idx,
+                      const EvalWell& cq_s_effective,
+                      Equations& eqns) const
+{
+    // subtract sum of phase fluxes in the well equations.
+    eqns.resWell_[seg][comp_idx] += cq_s_effective.value();
+
+    // assemble the jacobians
+    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+        // also need to consider the efficiency factor when manipulating the jacobians.
+        eqns.duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // input in transformed matrix
+
+        // the index name for the D should be eq_idx / pv_idx
+        eqns.duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
+    }
+
+    for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
+        // also need to consider the efficiency factor when manipulating the jacobians.
+        eqns.duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
+    }
+}
+
 #define INSTANCE(...) \
 template class MultisegmentWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;
 
diff --git a/opm/simulators/wells/MultisegmentWellAssemble.hpp b/opm/simulators/wells/MultisegmentWellAssemble.hpp
index 43e768614..693aa0c0e 100644
--- a/opm/simulators/wells/MultisegmentWellAssemble.hpp
+++ b/opm/simulators/wells/MultisegmentWellAssemble.hpp
@@ -124,6 +124,13 @@ public:
                             const EvalWell& inlet_rate,
                             Equations& eqns) const;
 
+    //! \brief Assemble equation for a perforation.
+    void assemblePerforationEq(const int seg,
+                               const int cell_idx,
+                               const int comp_idx,
+                               const EvalWell& cq_s_effective,
+                               Equations& eqns) const;
+
 private:
     const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well_; //!< Reference to well
 };
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index 9e1e3f97a..dd08be69b 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -1626,23 +1626,8 @@ namespace Opm
 
                     this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
 
-                    // subtract sum of phase fluxes in the well equations.
-                    this->linSys_.resWell_[seg][comp_idx] += cq_s_effective.value();
-
-                    // assemble the jacobians
-                    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
-
-                        // also need to consider the efficiency factor when manipulating the jacobians.
-                        this->linSys_.duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // intput in transformed matrix
-
-                        // the index name for the D should be eq_idx / pv_idx
-                        this->linSys_.duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
-                    }
-
-                    for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
-                        // also need to consider the efficiency factor when manipulating the jacobians.
-                        this->linSys_.duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
-                    }
+                    MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
+                        assemblePerforationEq(seg, cell_idx, comp_idx, cq_s_effective, this->linSys_);
                 }
             }