added: StandardWellAssemble

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

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,187 @@
+/*
+  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 {
+
+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>& eqns,
+                  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
+    eqns.resWell_[0][Bhp] = control_eq.value();
+    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
+        eqns.duneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + 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;
+
+// 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,71 @@
+/*
+  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;
+
+
+private:
+    const WellInterfaceFluidSystem<FluidSystem>& well_; //!< Reference to well
+};
+
+}
+
+#endif // OPM_STANDARDWELL_ASSEMBLE_HEADER_INCLUDED

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

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>
@@ -544,7 +545,18 @@ namespace Opm
 
         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.