Merge pull request #4276 from akva2/stdwell_assemble

Added: StandardWellAssemble

CMakeLists_files.cmake

@@ -96,6 +96,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/PerfData.cpp
   opm/simulators/wells/SegmentState.cpp
   opm/simulators/wells/SingleWellState.cpp
+  opm/simulators/wells/StandardWellAssemble.cpp
   opm/simulators/wells/StandardWellEquations.cpp
   opm/simulators/wells/StandardWellEval.cpp
   opm/simulators/wells/StandardWellGeneric.cpp
@@ -381,8 +382,10 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/SingleWellState.hpp
   opm/simulators/wells/StandardWell.hpp
   opm/simulators/wells/StandardWell_impl.hpp
+  opm/simulators/wells/StandardWellAssemble.hpp
   opm/simulators/wells/StandardWellEquations.hpp
   opm/simulators/wells/StandardWellEval.hpp
+  opm/simulators/wells/StandardWellGeneric.hpp
   opm/simulators/wells/TargetCalculator.hpp
   opm/simulators/wells/VFPHelpers.hpp
   opm/simulators/wells/VFPInjProperties.hpp

opm/simulators/wells/StandardWellAssemble.cpp

@@ -0,0 +1,337 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2016 - 2017 IRIS AS.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+#include <opm/simulators/wells/StandardWellAssemble.hpp>
+
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <opm/material/densead/DynamicEvaluation.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+
+#include <opm/models/blackoil/blackoilindices.hh>
+#include <opm/models/blackoil/blackoilonephaseindices.hh>
+#include <opm/models/blackoil/blackoiltwophaseindices.hh>
+
+#include <opm/simulators/wells/StandardWellEquations.hpp>
+#include <opm/simulators/wells/WellAssemble.hpp>
+#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
+#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
+
+namespace Opm {
+
+//! \brief Class administering assembler access to equation system.
+template<class Scalar, int numEq>
+class StandardWellEquationAccess {
+public:
+    //! \brief Constructor initializes reference to the equation system.
+    StandardWellEquationAccess(StandardWellEquations<Scalar,numEq>& eqns)
+        : eqns_(eqns)
+    {}
+
+    using BVectorWell = typename StandardWellEquations<Scalar,numEq>::BVectorWell;
+    using DiagMatWell = typename StandardWellEquations<Scalar,numEq>::DiagMatWell;
+    using OffDiatMatWell = typename StandardWellEquations<Scalar,numEq>::OffDiagMatWell;
+
+    //! \brief Returns a reference to residual vector.
+    BVectorWell& residual()
+    {
+        return eqns_.resWell_;
+    }
+
+    //! \brief Returns a reference to B matrix.
+    OffDiatMatWell& B()
+    {
+        return eqns_.duneB_;
+    }
+
+    //! \brief Returns a reference to C matrix.
+    OffDiatMatWell& C()
+    {
+        return eqns_.duneC_;
+    }
+
+    //! \brief Returns a reference to D matrix.
+    DiagMatWell& D()
+    {
+        return eqns_.duneD_;
+    }
+
+private:
+    StandardWellEquations<Scalar,numEq>& eqns_; //!< Reference to equation system
+};
+
+template<class FluidSystem, class Indices, class Scalar>
+template<class EvalWell>
+void
+StandardWellAssemble<FluidSystem,Indices,Scalar>::
+assembleControlEq(const WellState& well_state,
+                  const GroupState& group_state,
+                  const Schedule& schedule,
+                  const SummaryState& summaryState,
+                  const int numWellEq,
+                  const EvalWell& wqTotal,
+                  const EvalWell& bhp,
+                  const std::function<EvalWell(int)>& getQs,
+                  const double rho,
+                  const int Bhp,
+                  StandardWellEquations<Scalar,Indices::numEq>& eqns1,
+                  DeferredLogger& deferred_logger) const
+{
+    static constexpr int Water = BlackoilPhases::Aqua;
+    static constexpr int Oil = BlackoilPhases::Liquid;
+    static constexpr int Gas = BlackoilPhases::Vapour;
+    EvalWell control_eq(numWellEq + Indices::numEq, 0.0);
+
+    const auto& well = well_.wellEcl();
+
+    auto getRates = [&]() {
+        std::vector<EvalWell> rates(3, EvalWell(numWellEq + Indices::numEq, 0.0));
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+            rates[Water] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx));
+        }
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+            rates[Oil] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx));
+        }
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            rates[Gas] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx));
+        }
+        return rates;
+    };
+
+    if (well_.wellIsStopped()) {
+        control_eq = wqTotal;
+    } else if (well_.isInjector()) {
+        // Find injection rate.
+        const EvalWell injection_rate = wqTotal;
+                                                 // Setup function for evaluation of BHP from THP (used only if needed).
+        std::function<EvalWell()> bhp_from_thp = [&]() {
+            const auto rates = getRates();
+            return WellBhpThpCalculator(well_).calculateBhpFromThp(well_state,
+                                                                   rates,
+                                                                   well,
+                                                                   summaryState,
+                                                                   rho,
+                                                                   deferred_logger);
+        };
+
+             // Call generic implementation.
+        const auto& inj_controls = well.injectionControls(summaryState);
+        WellAssemble(well_).
+            assembleControlEqInj(well_state,
+                                 group_state,
+                                 schedule,
+                                 summaryState,
+                                 inj_controls,
+                                 bhp,
+                                 injection_rate,
+                                 bhp_from_thp,
+                                 control_eq,
+                                 deferred_logger);
+    } else {
+             // Find rates.
+        const auto rates = getRates();
+                                            // Setup function for evaluation of BHP from THP (used only if needed).
+        std::function<EvalWell()> bhp_from_thp = [&]() {
+            return WellBhpThpCalculator(well_).calculateBhpFromThp(well_state,
+                                                                   rates,
+                                                                   well,
+                                                                   summaryState,
+                                                                   rho,
+                                                                   deferred_logger);
+        };
+           // Call generic implementation.
+        const auto& prod_controls = well.productionControls(summaryState);
+        WellAssemble(well_).
+            assembleControlEqProd(well_state,
+                                  group_state,
+                                  schedule,
+                                  summaryState,
+                                  prod_controls,
+                                  bhp,
+                                  rates,
+                                  bhp_from_thp,
+                                  control_eq,
+                                  deferred_logger);
+    }
+
+    // using control_eq to update the matrix and residuals
+    // TODO: we should use a different index system for the well equations
+    StandardWellEquationAccess eqns(eqns1);
+    eqns.residual()[0][Bhp] = control_eq.value();
+    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+        eqns.D()[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
+    }
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+template<class EvalWell>
+void StandardWellAssemble<FluidSystem,Indices,Scalar>::
+assembleInjectivityEq(const EvalWell& eq_pskin,
+                      const EvalWell& eq_wat_vel,
+                      const int pskin_index,
+                      const int wat_vel_index,
+                      const int cell_idx,
+                      const int numWellEq,
+                      StandardWellEquations<Scalar,Indices::numEq>& eqns1) const
+{
+    StandardWellEquationAccess eqns(eqns1);
+    eqns.residual()[0][pskin_index] = eq_pskin.value();
+    eqns.residual()[0][wat_vel_index] = eq_wat_vel.value();
+    for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
+        eqns.D()[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
+        eqns.D()[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
+    }
+
+         // the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B
+    for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
+        eqns.B()[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
+    }
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+template<class EvalWell>
+void StandardWellAssemble<FluidSystem,Indices,Scalar>::
+assemblePerforationEq(const EvalWell& cq_s_effective,
+                      const int componentIdx,
+                      const int cell_idx,
+                      const int numWellEq,
+                      StandardWellEquations<Scalar,Indices::numEq>& eqns1) const
+{
+    StandardWellEquationAccess eqns(eqns1);
+
+    // subtract sum of phase fluxes in the well equations.
+    eqns.residual()[0][componentIdx] += cq_s_effective.value();
+
+       // assemble the jacobians
+    for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
+        // also need to consider the efficiency factor when manipulating the jacobians.
+        eqns.C()[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq); // intput in transformed matrix
+        eqns.D()[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
+    }
+
+    for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
+        eqns.B()[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
+    }
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+template<class EvalWell>
+void StandardWellAssemble<FluidSystem,Indices,Scalar>::
+assembleSourceEq(const EvalWell& resWell_loc,
+                 const int componentIdx,
+                 const int numWellEq,
+                 StandardWellEquations<Scalar,Indices::numEq>& eqns1) const
+{
+    StandardWellEquationAccess eqns(eqns1);
+    for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
+        eqns.D()[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
+    }
+    eqns.residual()[0][componentIdx] += resWell_loc.value();
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+template<class EvalWell>
+void StandardWellAssemble<FluidSystem,Indices,Scalar>::
+assembleZFracEq(const EvalWell& cq_s_zfrac_effective,
+                const int cell_idx,
+                const int numWellEq,
+                StandardWellEquations<Scalar,Indices::numEq>& eqns1) const
+{
+    StandardWellEquationAccess eqns(eqns1);
+    for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
+        eqns.C()[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
+    }
+}
+
+#define INSTANCE(Dim,...) \
+template class StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>; \
+template void \
+StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>:: \
+    assembleControlEq(const WellState&, \
+                      const GroupState&, \
+                      const Schedule&, \
+                      const SummaryState&, \
+                      const int, \
+                      const DenseAd::Evaluation<double,-1,Dim>&, \
+                      const DenseAd::Evaluation<double,-1,Dim>&, \
+                      const std::function<DenseAd::Evaluation<double,-1,Dim>(int)>&, \
+                      const double, \
+                      const int, \
+                      StandardWellEquations<double,__VA_ARGS__::numEq>&, \
+                      DeferredLogger&) const; \
+template void \
+StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>:: \
+    assembleInjectivityEq(const DenseAd::Evaluation<double,-1,Dim>&, \
+                          const DenseAd::Evaluation<double,-1,Dim>&, \
+                          const int, \
+                          const int, \
+                          const int, \
+                          const int, \
+                          StandardWellEquations<double,__VA_ARGS__::numEq>&) const; \
+template void \
+StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>:: \
+    assemblePerforationEq(const DenseAd::Evaluation<double,-1,Dim>&, \
+                          const int, \
+                          const int, \
+                          const int, \
+                          StandardWellEquations<double,__VA_ARGS__::numEq>&) const; \
+template void \
+StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>:: \
+    assembleZFracEq(const DenseAd::Evaluation<double,-1,Dim>&, \
+                    const int, \
+                    const int, \
+                    StandardWellEquations<double,__VA_ARGS__::numEq>&) const; \
+template void \
+StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>:: \
+  assembleSourceEq(const DenseAd::Evaluation<double,-1,Dim>&, \
+                   const int, \
+                   const int, \
+                   StandardWellEquations<double,__VA_ARGS__::numEq>&) const;
+
+// One phase
+INSTANCE(4u, BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
+INSTANCE(5u, BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
+INSTANCE(9u, BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>)
+
+     // Two phase
+INSTANCE(6u, BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>)
+INSTANCE(6u, BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
+INSTANCE(6u, BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>)
+INSTANCE(7u, BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>)
+INSTANCE(7u, BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,true,0u,2u,0u>)
+INSTANCE(7u, BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
+INSTANCE(7u, BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>)
+INSTANCE(7u, BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>)
+INSTANCE(8u, BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>)
+
+     // Blackoil
+INSTANCE(8u, BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>)
+INSTANCE(9u, BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>)
+INSTANCE(9u, BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>)
+INSTANCE(9u, BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>)
+INSTANCE(9u, BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>)
+INSTANCE(9u, BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>)
+INSTANCE(10u, BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>)
+INSTANCE(10u, BlackOilIndices<0u,0u,0u,1u,false,true,0u,0u>)
+INSTANCE(10u, BlackOilIndices<0u,0u,0u,1u,false,false,1u,0u>)
+
+}

opm/simulators/wells/StandardWellAssemble.hpp

@@ -0,0 +1,103 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2016 - 2017 IRIS AS.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef OPM_STANDARDWELL_ASSEMBLE_HEADER_INCLUDED
+#define OPM_STANDARDWELL_ASSEMBLE_HEADER_INCLUDED
+
+#include <functional>
+
+namespace Opm
+{
+
+class DeferredLogger;
+class GroupState;
+class Schedule;
+template<class Scalar, int numEq> class StandardWellEquations;
+class SummaryState;
+template<class FluidSystem> class WellInterfaceFluidSystem;
+class WellState;
+
+//! \brief Class handling assemble of the equation system for StandardWell.
+template<class FluidSystem, class Indices, class Scalar>
+class StandardWellAssemble
+{
+public:
+    //! \brief Constructor initializes reference to well.
+    StandardWellAssemble(const WellInterfaceFluidSystem<FluidSystem>& well)
+        : well_(well)
+    {}
+
+    //! \brief Assemble control equation.
+    template<class EvalWell>
+    void assembleControlEq(const WellState& well_state,
+                           const GroupState& group_state,
+                           const Schedule& schedule,
+                           const SummaryState& summaryState,
+                           const int numWellEq,
+                           const EvalWell& wqTotal,
+                           const EvalWell& bhp,
+                           const std::function<EvalWell(int)>& getQs,
+                           const double rho,
+                           const int Bhp,
+                           StandardWellEquations<Scalar,Indices::numEq>& eqns,
+                           DeferredLogger& deferred_logger) const;
+
+    //! \brief Assemble injectivity equation.
+    template<class EvalWell>
+    void assembleInjectivityEq(const EvalWell& eq_pskin,
+                               const EvalWell& eq_wat_vel,
+                               const int pskin_index,
+                               const int wat_vel_index,
+                               const int cell_idx,
+                               const int numWellEq,
+                               StandardWellEquations<Scalar,Indices::numEq>& eqns) const;
+
+    //! \brief Assemble equation for a perforation.
+    template<class EvalWell>
+    void assemblePerforationEq(const EvalWell& cq_s_effective,
+                               const int componentIdx,
+                               const int cell_idx,
+                               const int numWellEq,
+                               StandardWellEquations<Scalar,Indices::numEq>& eqns) const;
+
+    //! \brief Assemble equation for Z fraction.
+    template<class EvalWell>
+    void assembleZFracEq(const EvalWell& cq_s_zfrac_effective,
+                         const int cell_idx,
+                         const int numWellEq,
+                         StandardWellEquations<Scalar,Indices::numEq>& eqns) const;
+
+    //! \brief Assemble a source term.
+    template<class EvalWell>
+    void assembleSourceEq(const EvalWell& resWell_loc,
+                          const int componentIdx,
+                          const int numWellEq,
+                          StandardWellEquations<Scalar,Indices::numEq>& eqns) const;
+
+
+private:
+    const WellInterfaceFluidSystem<FluidSystem>& well_; //!< Reference to well
+};
+
+}
+
+#endif // OPM_STANDARDWELL_ASSEMBLE_HEADER_INCLUDED

opm/simulators/wells/StandardWellEquations.hpp

@@ -35,6 +35,7 @@ namespace Opm
 {
 
 class ParallelWellInfo;
+template<class Scalar, int numEq> class StandardWellEquationAccess;
 class WellContributions;
 class WellInterfaceGeneric;
 class WellState;
@@ -117,6 +118,15 @@ public:
     //! \brief Sum with off-process contribution.
     void sumDistributed(Parallel::Communication comm);
 
+    //! \brief Returns a const reference to the residual.
+    const BVectorWell& residual() const
+    {
+        return resWell_;
+    }
+
+private:
+    friend class StandardWellEquationAccess<Scalar,numEq>;
+
     // two off-diagonal matrices
     OffDiagMatWell duneB_;
     OffDiagMatWell duneC_;

opm/simulators/wells/StandardWellEval.cpp

@@ -32,7 +32,6 @@
 #include <opm/simulators/timestepping/ConvergenceReport.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
-#include <opm/simulators/wells/WellAssemble.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
 #include <opm/simulators/wells/WellConvergence.hpp>
 #include <opm/simulators/wells/WellInterfaceIndices.hpp>
@@ -358,100 +357,6 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
     primary_variables_[Bhp] = ws.bhp;
 }
 
-template<class FluidSystem, class Indices, class Scalar>
-void
-StandardWellEval<FluidSystem,Indices,Scalar>::
-assembleControlEq(const WellState& well_state,
-                  const GroupState& group_state,
-                  const Schedule& schedule,
-                  const SummaryState& summaryState,
-                  DeferredLogger& deferred_logger)
-{
-    static constexpr int Gas = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Gas;
-    static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
-    static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
-    EvalWell control_eq(numWellEq_ + Indices::numEq, 0.0);
-
-    const auto& well = baseif_.wellEcl();
-
-    auto getRates = [&]() {
-        std::vector<EvalWell> rates(3, EvalWell(numWellEq_ + Indices::numEq, 0.0));
-        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-            rates[Water] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx));
-        }
-        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
-            rates[Oil] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx));
-        }
-        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-            rates[Gas] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx));
-        }
-        return rates;
-    };
-
-    if (baseif_.wellIsStopped()) {
-        control_eq = getWQTotal();
-    } else if (baseif_.isInjector()) {
-        // Find injection rate.
-        const EvalWell injection_rate = getWQTotal();
-        // Setup function for evaluation of BHP from THP (used only if needed).
-        std::function<EvalWell()> bhp_from_thp = [&]() {
-            const auto rates = getRates();
-            return WellBhpThpCalculator(baseif_).calculateBhpFromThp(well_state,
-                                                                     rates,
-                                                                     well,
-                                                                     summaryState,
-                                                                     this->getRho(),
-                                                                     deferred_logger);
-        };
-        // Call generic implementation.
-        const auto& inj_controls = well.injectionControls(summaryState);
-        WellAssemble(baseif_).
-             assembleControlEqInj(well_state,
-                                  group_state,
-                                  schedule,
-                                  summaryState,
-                                  inj_controls,
-                                  getBhp(),
-                                  injection_rate,
-                                  bhp_from_thp,
-                                  control_eq,
-                                  deferred_logger);
-    } else {
-        // Find rates.
-        const auto rates = getRates();
-        // Setup function for evaluation of BHP from THP (used only if needed).
-        std::function<EvalWell()> bhp_from_thp = [&]() {
-             return WellBhpThpCalculator(baseif_).calculateBhpFromThp(well_state,
-                                                                      rates,
-                                                                      well,
-                                                                      summaryState,
-                                                                      this->getRho(),
-                                                                      deferred_logger);
-        };
-        // Call generic implementation.
-        const auto& prod_controls = well.productionControls(summaryState);
-        WellAssemble(baseif_).
-            assembleControlEqProd(well_state,
-                                  group_state,
-                                  schedule,
-                                  summaryState,
-                                  prod_controls,
-                                  getBhp(),
-                                  rates,
-                                  bhp_from_thp,
-                                  control_eq,
-                                  deferred_logger);
-    }
-
-    // using control_eq to update the matrix and residuals
-    // TODO: we should use a different index system for the well equations
-    this->linSys_.resWell_[0][Bhp] = control_eq.value();
-    for (int pv_idx = 0; pv_idx < numWellEq_; ++pv_idx) {
-        this->linSys_.duneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
-    }
-}
-
-
 template<class FluidSystem, class Indices, class Scalar>
 void
 StandardWellEval<FluidSystem,Indices,Scalar>::
@@ -774,7 +679,7 @@ getWellConvergence(const WellState& well_state,
     res.resize(numWellEq_);
     for (int eq_idx = 0; eq_idx < numWellEq_; ++eq_idx) {
         // magnitude of the residual matters
-        res[eq_idx] = std::abs(this->linSys_.resWell_[0][eq_idx]);
+        res[eq_idx] = std::abs(this->linSys_.residual()[0][eq_idx]);
     }
 
     std::vector<double> well_flux_residual(baseif_.numComponents());
@@ -811,7 +716,7 @@ getWellConvergence(const WellState& well_state,
     WellConvergence(baseif_).
         checkConvergenceControlEq(well_state,
                                   {1.e3, 1.e4, 1.e-4, 1.e-6, maxResidualAllowed},
-                                  std::abs(this->linSys_.resWell_[0][Bhp]),
+                                  std::abs(this->linSys_.residual()[0][Bhp]),
                                   report,
                                   deferred_logger);
 

opm/simulators/wells/StandardWellEval.hpp

@@ -133,12 +133,6 @@ protected:
     static double relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
                                                     const BVectorWell& dwells);
 
-    void assembleControlEq(const WellState& well_state,
-                           const GroupState& group_state,
-                           const Schedule& schedule,
-                           const SummaryState& summaryState,
-                           DeferredLogger& deferred_logger);
-
     // computing the accumulation term for later use in well mass equations
     void computeAccumWell();
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -20,6 +20,7 @@
 */
 
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
+#include <opm/simulators/wells/StandardWellAssemble.hpp>
 #include <opm/simulators/wells/VFPHelpers.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
 #include <opm/simulators/wells/WellConvergence.hpp>
@@ -482,19 +483,12 @@ namespace Opm
 
                 connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
 
-                // subtract sum of phase fluxes in the well equations.
-                this->linSys_.resWell_[0][componentIdx] += cq_s_effective.value();
-
-                // assemble the jacobians
-                for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-                    // also need to consider the efficiency factor when manipulating the jacobians.
-                    this->linSys_.duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq); // intput in transformed matrix
-                    this->linSys_.duneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
-                }
-
-                for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
-                    this->linSys_.duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
-                }
+                StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
+                    assemblePerforationEq(cq_s_effective,
+                                          componentIdx,
+                                          cell_idx,
+                                          this->numWellEq_,
+                                          this->linSys_);
 
                 // Store the perforation phase flux for later usage.
                 if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
@@ -506,9 +500,11 @@ namespace Opm
             }
 
             if constexpr (has_zFraction) {
-                for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-                    this->linSys_.duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
-                }
+                StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
+                    assembleZFracEq(cq_s_zfrac_effective,
+                                    cell_idx,
+                                    this->numWellEq_,
+                                    this->linSys_);
             }
         }
         // Update the connection
@@ -536,15 +532,27 @@ namespace Opm
                 resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt;
             }
             resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
-            for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-                this->linSys_.duneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
-            }
-            this->linSys_.resWell_[0][componentIdx] += resWell_loc.value();
+            StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
+                assembleSourceEq(resWell_loc,
+                                 componentIdx,
+                                 this->numWellEq_,
+                                 this->linSys_);
         }
 
         const auto& summaryState = ebosSimulator.vanguard().summaryState();
         const Schedule& schedule = ebosSimulator.vanguard().schedule();
-        this->assembleControlEq(well_state, group_state, schedule, summaryState, deferred_logger);
+        std::function<EvalWell(int)> gQ = [this](int a) { return this->getQs(a); };
+        StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
+            assembleControlEq(well_state, group_state,
+                              schedule, summaryState,
+                              this->numWellEq_,
+                              this->getWQTotal(),
+                              this->getBhp(),
+                              gQ,
+                              this->getRho(),
+                              Bhp,
+                              this->linSys_,
+                              deferred_logger);
 
 
         // do the local inversion of D.
@@ -2313,7 +2321,6 @@ namespace Opm
 
         // equation for the water velocity
         const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
-        this->linSys_.resWell_[0][wat_vel_index] = eq_wat_vel.value();
 
         const auto& ws = well_state.well(this->index_of_well_);
         const auto& perf_data = ws.perf_data;
@@ -2328,16 +2335,14 @@ namespace Opm
         const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
                                   - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
 
-        this->linSys_.resWell_[0][pskin_index] = eq_pskin.value();
-        for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-            this->linSys_.duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
-            this->linSys_.duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
-        }
-
-        // the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B
-        for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
-            this->linSys_.duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
-        }
+        StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
+                assembleInjectivityEq(eq_pskin,
+                                      eq_wat_vel,
+                                      pskin_index,
+                                      wat_vel_index,
+                                      cell_idx,
+                                      this->numWellEq_,
+                                      this->linSys_);
     }