added: MultisegmentWellAssemble

this handles assembly of the equation system for multisegment well.
start by moving assembleControlEq into the new class

CMakeLists_files.cmake

@@ -90,6 +90,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/GlobalWellInfo.cpp
   opm/simulators/wells/GroupState.cpp
   opm/simulators/wells/MSWellHelpers.cpp
+  opm/simulators/wells/MultisegmentWellAssemble.cpp
   opm/simulators/wells/MultisegmentWellEquations.cpp
   opm/simulators/wells/MultisegmentWellEval.cpp
   opm/simulators/wells/MultisegmentWellGeneric.cpp
@@ -372,6 +373,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/MSWellHelpers.hpp
   opm/simulators/wells/MultisegmentWell.hpp
   opm/simulators/wells/MultisegmentWell_impl.hpp
+  opm/simulators/wells/MultisegmentWellAssemble.hpp
   opm/simulators/wells/MultisegmentWellEquations.hpp
   opm/simulators/wells/MultisegmentWellEval.hpp
   opm/simulators/wells/MultisegmentWellGeneric.hpp

opm/simulators/wells/MultisegmentWellAssemble.cpp

@@ -0,0 +1,186 @@
+/*
+  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/MultisegmentWellAssemble.hpp>
+
+#include <opm/core/props/BlackoilPhases.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/MultisegmentWellEquations.hpp>
+#include <opm/simulators/wells/WellAssemble.hpp>
+#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
+#include <opm/simulators/wells/WellInterfaceIndices.hpp>
+
+namespace Opm {
+
+template<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellAssemble<FluidSystem,Indices,Scalar>::
+assembleControlEq(const WellState& well_state,
+                  const GroupState& group_state,
+                  const Schedule& schedule,
+                  const SummaryState& summaryState,
+                  const Well::InjectionControls& inj_controls,
+                  const Well::ProductionControls& prod_controls,
+                  const double rho,
+                  const EvalWell& wqTotal,
+                  const EvalWell& bhp,
+                  const int SPres,
+                  const std::function<EvalWell(const int)>& getQs,
+                  Equations& eqns,
+                  DeferredLogger& deferred_logger) const
+{
+    static constexpr int Gas = BlackoilPhases::Vapour;
+    static constexpr int Oil = BlackoilPhases::Liquid;
+    static constexpr int Water = BlackoilPhases::Aqua;
+
+    EvalWell control_eq(0.0);
+
+    const auto& well = well_.wellEcl();
+
+    auto getRates = [&]() {
+        std::vector<EvalWell> rates(3, 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 scaling factor to get injection rate,
+        const InjectorType injectorType = inj_controls.injector_type;
+        double scaling = 1.0;
+        const auto& pu = well_.phaseUsage();
+        switch (injectorType) {
+        case InjectorType::WATER:
+        {
+            scaling = well_.scalingFactor(pu.phase_pos[BlackoilPhases::Aqua]);
+            break;
+        }
+        case InjectorType::OIL:
+        {
+            scaling = well_.scalingFactor(pu.phase_pos[BlackoilPhases::Liquid]);
+            break;
+        }
+        case InjectorType::GAS:
+        {
+            scaling = well_.scalingFactor(pu.phase_pos[BlackoilPhases::Vapour]);
+            break;
+        }
+        default:
+            throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
+        }
+        const EvalWell injection_rate = wqTotal / scaling;
+        // 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.
+        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.
+        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
+    eqns.resWell_[0][SPres] = control_eq.value();
+    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+        eqns.duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
+    }
+}
+
+#define INSTANCE(...) \
+template class MultisegmentWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;
+
+// One phase
+INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
+INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
+INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>)
+
+// Two phase
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
+
+// Blackoil
+INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>)
+INSTANCE(BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>)
+INSTANCE(BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,true,0u,0u>)
+
+}

opm/simulators/wells/MultisegmentWellAssemble.hpp

@@ -0,0 +1,76 @@
+/*
+  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_MULTISEGMENTWELL_ASSEMBLE_HEADER_INCLUDED
+#define OPM_MULTISEGMENTWELL_ASSEMBLE_HEADER_INCLUDED
+
+#include <opm/input/eclipse/Schedule/Well/Well.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+
+#include <functional>
+
+namespace Opm
+{
+
+class DeferredLogger;
+class GroupState;
+template<class Scalar, int numWellEq, int numEq> class MultisegmentWellEquations;
+class Schedule;
+class SummaryState;
+template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
+class WellState;
+
+//! \brief Class handling assemble of the equation system for MultisegmentWell.
+template<class FluidSystem, class Indices, class Scalar>
+class MultisegmentWellAssemble
+{
+public:
+    static constexpr int numWellEq = Indices::numPhases+1;
+    using Equations = MultisegmentWellEquations<Scalar,numWellEq,Indices::numEq>;
+    using EvalWell = DenseAd::Evaluation<Scalar, numWellEq+Indices::numEq>;
+    //! \brief Constructor initializes reference to well.
+    MultisegmentWellAssemble(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well)
+        : well_(well)
+    {}
+
+    //! \brief Assemble control equation.
+    void assembleControlEq(const WellState& well_state,
+                           const GroupState& group_state,
+                           const Schedule& schedule,
+                           const SummaryState& summaryState,
+                           const Well::InjectionControls& inj_controls,
+                           const Well::ProductionControls& prod_controls,
+                           const double rho,
+                           const EvalWell& wqTotal,
+                           const EvalWell& bhp,
+                           const int SPres,
+                           const std::function<EvalWell(const int)>& getQs,
+                           Equations& eqns,
+                           DeferredLogger& deferred_logger) const;
+
+private:
+    const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well_; //!< Reference to well
+};
+
+}
+
+#endif // OPM_STANDARDWELL_ASSEMBLE_HEADER_INCLUDED

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -1107,120 +1107,6 @@ getSegmentSurfaceVolume(const EvalWell& temperature,
     return volume / vol_ratio;
 }
 
-template<typename FluidSystem, typename Indices, typename Scalar>
-void
-MultisegmentWellEval<FluidSystem,Indices,Scalar>::
-assembleControlEq(const WellState& well_state,
-                  const GroupState& group_state,
-                  const Schedule& schedule,
-                  const SummaryState& summaryState,
-                  const Well::InjectionControls& inj_controls,
-                  const Well::ProductionControls& prod_controls,
-                  const double rho,
-                  DeferredLogger& deferred_logger)
-{
-    static constexpr int Gas = BlackoilPhases::Vapour;
-    static constexpr int Oil = BlackoilPhases::Liquid;
-    static constexpr int Water = BlackoilPhases::Aqua;
-
-    EvalWell control_eq(0.0);
-
-    const auto& well = baseif_.wellEcl();
-
-    auto getRates = [&]() {
-        std::vector<EvalWell> rates(3, 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 scaling factor to get injection rate,
-        const InjectorType injectorType = inj_controls.injector_type;
-        double scaling = 1.0;
-        const auto& pu = baseif_.phaseUsage();
-        switch (injectorType) {
-        case InjectorType::WATER:
-        {
-            scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Aqua]);
-            break;
-        }
-        case InjectorType::OIL:
-        {
-            scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Liquid]);
-            break;
-        }
-        case InjectorType::GAS:
-        {
-            scaling = baseif_.scalingFactor(pu.phase_pos[BlackoilPhases::Vapour]);
-            break;
-        }
-        default:
-            throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
-        }
-        const EvalWell injection_rate = getWQTotal() / scaling;
-        // 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,
-                                                                     rho,
-                                                                     deferred_logger);
-        };
-        // Call generic implementation.
-        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,
-                                                                     rho,
-                                                                     deferred_logger);
-        };
-        // Call generic implementation.
-        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
-    linSys_.resWell_[0][SPres] = control_eq.value();
-    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
-        linSys_.duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
-    }
-}
-
 template<typename FluidSystem, typename Indices, typename Scalar>
 void
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::

opm/simulators/wells/MultisegmentWellEval.hpp

@@ -105,15 +105,6 @@ protected:
     void initMatrixAndVectors(const int num_cells);
     void initPrimaryVariablesEvaluation() const;
 
-    void assembleControlEq(const WellState& well_state,
-                           const GroupState& group_state,
-                           const Schedule& schedule,
-                           const SummaryState& summaryState,
-                           const Well::InjectionControls& inj_controls,
-                           const Well::ProductionControls& prod_controls,
-                           const double rho,
-                           DeferredLogger& deferred_logger);
-
     void assembleDefaultPressureEq(const int seg,
                                    WellState& well_state);
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -19,6 +19,7 @@
 */
 
 
+#include <opm/simulators/wells/MultisegmentWellAssemble.hpp>
 #include <opm/simulators/wells/WellBhpThpCalculator.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
@@ -1671,13 +1672,20 @@ namespace Opm
             if (seg == 0) { // top segment, pressure equation is the control equation
                 const auto& summaryState = ebosSimulator.vanguard().summaryState();
                 const Schedule& schedule = ebosSimulator.vanguard().schedule();
-                this->assembleControlEq(well_state,
+                std::function<EvalWell(const int)> gQ = [this](int a) { return this->getQs(a); };
+                MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(*this).
+                        assembleControlEq(well_state,
                                         group_state,
                                         schedule,
                                         summaryState,
                                         inj_controls,
                                         prod_controls,
                                         getRefDensity(),
+                                        this->getWQTotal(),
+                                        this->getBhp(),
+                                        SPres,
+                                        gQ,
+                                        this->linSys_,
                                         deferred_logger);
             } else {
                 const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();