added: StandardWellPrimaryVariables

this is a container class for the primary variables in
standard well

CMakeLists_files.cmake

@@ -100,6 +100,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/StandardWellEquations.cpp
   opm/simulators/wells/StandardWellEval.cpp
   opm/simulators/wells/StandardWellGeneric.cpp
+  opm/simulators/wells/StandardWellPrimaryVariables.cpp
   opm/simulators/wells/TargetCalculator.cpp
   opm/simulators/wells/VFPHelpers.cpp
   opm/simulators/wells/VFPProdProperties.cpp
@@ -386,6 +387,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/StandardWellEquations.hpp
   opm/simulators/wells/StandardWellEval.hpp
   opm/simulators/wells/StandardWellGeneric.hpp
+  opm/simulators/wells/StandardWellPrimaryVariables.hpp
   opm/simulators/wells/TargetCalculator.hpp
   opm/simulators/wells/VFPHelpers.hpp
   opm/simulators/wells/VFPInjProperties.hpp

opm/simulators/wells/StandardWellEval.cpp

@@ -51,27 +51,18 @@ StandardWellEval<FluidSystem,Indices,Scalar>::
 StandardWellEval(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
     : StandardWellGeneric<Scalar>(baseif)
     , baseif_(baseif)
+    , primary_variables_(baseif_)
     , F0_(numWellConservationEq)
     , linSys_(baseif_.parallelWellInfo())
 {
 }
 
-template<class FluidSystem, class Indices, class Scalar>
-void StandardWellEval<FluidSystem,Indices,Scalar>::
-initPrimaryVariablesEvaluation() const
-{
-    for (int eqIdx = 0; eqIdx < numWellEq_; ++eqIdx) {
-        primary_variables_evaluation_[eqIdx] =
-            EvalWell::createVariable(numWellEq_ + Indices::numEq, primary_variables_[eqIdx], Indices::numEq + eqIdx);
-    }
-}
-
 template<class FluidSystem, class Indices, class Scalar>
 typename StandardWellEval<FluidSystem,Indices,Scalar>::EvalWell
 StandardWellEval<FluidSystem,Indices,Scalar>::
 extendEval(const Eval& in) const
 {
-    EvalWell out(numWellEq_ + Indices::numEq, in.value());
+    EvalWell out(primary_variables_.numWellEq() + Indices::numEq, in.value());
     for(int eqIdx = 0; eqIdx < Indices::numEq;++eqIdx) {
         out.setDerivative(eqIdx, in.derivative(eqIdx));
     }
@@ -131,46 +122,46 @@ StandardWellEval<FluidSystem,Indices,Scalar>::
 wellVolumeFraction(const unsigned compIdx) const
 {
     if (FluidSystem::numActivePhases() == 1) {
-        return EvalWell(numWellEq_ + Indices::numEq, 1.0);
+        return EvalWell(primary_variables_.numWellEq() + Indices::numEq, 1.0);
     }
 
     if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
         if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
-            return primary_variables_evaluation_[WFrac];
+            return primary_variables_.evaluation_[WFrac];
         }
 
         if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
-            return primary_variables_evaluation_[GFrac];
+            return primary_variables_.evaluation_[GFrac];
         }
 
         if (Indices::enableSolvent && compIdx == (unsigned)Indices::contiSolventEqIdx) {
-            return primary_variables_evaluation_[SFrac];
+            return primary_variables_.evaluation_[SFrac];
         }
     }
     else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
         if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
-            return primary_variables_evaluation_[GFrac];
+            return primary_variables_.evaluation_[GFrac];
         }
     }
 
     // Oil or WATER fraction
-    EvalWell well_fraction(numWellEq_ + Indices::numEq, 1.0);
+    EvalWell well_fraction(primary_variables_.numWellEq() + Indices::numEq, 1.0);
     if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-            well_fraction -= primary_variables_evaluation_[WFrac];
+            well_fraction -= primary_variables_.evaluation_[WFrac];
         }
 
         if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-            well_fraction -= primary_variables_evaluation_[GFrac];
+            well_fraction -= primary_variables_.evaluation_[GFrac];
         }
 
         if (Indices::enableSolvent) {
-            well_fraction -= primary_variables_evaluation_[SFrac];
+            well_fraction -= primary_variables_.evaluation_[SFrac];
         }
     }
     else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))) {
 
-            well_fraction -= primary_variables_evaluation_[GFrac];
+            well_fraction -= primary_variables_.evaluation_[GFrac];
     }
 
     return well_fraction;
@@ -217,9 +208,9 @@ getQs(const int comp_idx) const
             //                         "Multi phase injectors are not supported, requested for well " + name());
             break;
         }
-        return inj_frac * primary_variables_evaluation_[WQTotal];
+        return inj_frac * primary_variables_.evaluation_[WQTotal];
     } else { // producers
-        return primary_variables_evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx);
+        return primary_variables_.evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx);
     }
 }
 
@@ -242,7 +233,7 @@ typename StandardWellEval<FluidSystem,Indices,Scalar>::EvalWell
 StandardWellEval<FluidSystem,Indices,Scalar>::
 wellSurfaceVolumeFraction(const int compIdx) const
 {
-    EvalWell sum_volume_fraction_scaled(numWellEq_ + Indices::numEq, 0.);
+    EvalWell sum_volume_fraction_scaled(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
     for (int idx = 0; idx < baseif_.numComponents(); ++idx) {
         sum_volume_fraction_scaled += wellVolumeFractionScaled(idx);
     }
@@ -278,13 +269,13 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
     if (baseif_.isInjector()) {
         switch (baseif_.wellEcl().injectorType()) {
         case InjectorType::WATER:
-            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]];
+            primary_variables_.value_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]];
             break;
         case InjectorType::GAS:
-            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]];
+            primary_variables_.value_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]];
             break;
         case InjectorType::OIL:
-            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]];
+            primary_variables_.value_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]];
             break;
         case InjectorType::MULTI:
             // Not supported.
@@ -293,19 +284,19 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
             break;
         }
     } else {
-            primary_variables_[WQTotal] = total_well_rate;
+            primary_variables_.value_[WQTotal] = total_well_rate;
     }
 
     if (std::abs(total_well_rate) > 0.) {
         if constexpr (has_wfrac_variable) {
-            primary_variables_[WFrac] = baseif_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
+            primary_variables_.value_[WFrac] = baseif_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
         }
         if constexpr (has_gfrac_variable) {
-            primary_variables_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (ws.surface_rates[pu.phase_pos[Gas]]
+            primary_variables_.value_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (ws.surface_rates[pu.phase_pos[Gas]]
                                                                                     - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate ;
         }
         if constexpr (Indices::enableSolvent) {
-            primary_variables_[SFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * ws.sum_solvent_rates() / total_well_rate ;
+            primary_variables_.value_[SFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * ws.sum_solvent_rates() / total_well_rate ;
         }
     } else { // total_well_rate == 0
         if (baseif_.isInjector()) {
@@ -314,9 +305,9 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
                 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                     auto phase = baseif_.wellEcl().getInjectionProperties().injectorType;
                     if (phase == InjectorType::WATER) {
-                        primary_variables_[WFrac] = 1.0;
+                        primary_variables_.value_[WFrac] = 1.0;
                     } else {
-                        primary_variables_[WFrac] = 0.0;
+                        primary_variables_.value_[WFrac] = 0.0;
                     }
                 }
             }
@@ -324,13 +315,13 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
                 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                     auto phase = baseif_.wellEcl().getInjectionProperties().injectorType;
                     if (phase == InjectorType::GAS) {
-                        primary_variables_[GFrac] = (1.0 - baseif_.rsRvInj());
+                        primary_variables_.value_[GFrac] = (1.0 - baseif_.rsRvInj());
                         if constexpr (Indices::enableSolvent) {
-                            primary_variables_[GFrac] = 1.0 - baseif_.rsRvInj() - baseif_.wsolvent();
+                            primary_variables_.value_[GFrac] = 1.0 - baseif_.rsRvInj() - baseif_.wsolvent();
-                            primary_variables_[SFrac] = baseif_.wsolvent();
+                            primary_variables_.value_[SFrac] = baseif_.wsolvent();
                         }
                     } else {
-                        primary_variables_[GFrac] = 0.0;
+                        primary_variables_.value_[GFrac] = 0.0;
                     }
                 }
             }
@@ -341,11 +332,11 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
         } else if (baseif_.isProducer()) { // producers
             // TODO: the following are not addressed for the solvent case yet
             if constexpr (has_wfrac_variable) {
-                primary_variables_[WFrac] = 1.0 / np;
+                primary_variables_.value_[WFrac] = 1.0 / np;
             }
 
             if constexpr (has_gfrac_variable) {
-                primary_variables_[GFrac] = 1.0 / np;
+                primary_variables_.value_[GFrac] = 1.0 / np;
             }
         } else {
             OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
@@ -354,7 +345,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
 
 
     // BHP
-    primary_variables_[Bhp] = ws.bhp;
+    primary_variables_.value_[Bhp] = ws.bhp;
 }
 
 template<class FluidSystem, class Indices, class Scalar>
@@ -369,10 +360,10 @@ updatePrimaryVariablesPolyMW(const BVectorWell& dwells) const
 
             const double relaxation_factor = 0.9;
             const double dx_wat_vel = dwells[0][wat_vel_index];
-            primary_variables_[wat_vel_index] -= relaxation_factor * dx_wat_vel;
+            primary_variables_.value_[wat_vel_index] -= relaxation_factor * dx_wat_vel;
 
             const double dx_pskin = dwells[0][pskin_index];
-            primary_variables_[pskin_index] -= relaxation_factor * dx_pskin;
+            primary_variables_.value_[pskin_index] -= relaxation_factor * dx_pskin;
         }
     }
 }
@@ -392,12 +383,12 @@ processFractions() const
         F[pu.phase_pos[Oil]] = 1.0;
 
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-            F[pu.phase_pos[Water]] = primary_variables_[WFrac];
+            F[pu.phase_pos[Water]] = primary_variables_.value_[WFrac];
             F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
         }
 
         if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-            F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
+            F[pu.phase_pos[Gas]] = primary_variables_.value_[GFrac];
             F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
         }
     }
@@ -405,7 +396,7 @@ processFractions() const
         F[pu.phase_pos[Water]] = 1.0;
 
         if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-            F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
+            F[pu.phase_pos[Gas]] = primary_variables_.value_[GFrac];
             F[pu.phase_pos[Water]] -= F[pu.phase_pos[Gas]];
         }
     }
@@ -415,7 +406,7 @@ processFractions() const
 
     [[maybe_unused]] double F_solvent;
     if constexpr (Indices::enableSolvent) {
-        F_solvent = primary_variables_[SFrac];
+        F_solvent = primary_variables_.value_[SFrac];
         F[pu.phase_pos[Oil]] -= F_solvent;
     }
 
@@ -465,14 +456,14 @@ processFractions() const
     }
 
     if constexpr (has_wfrac_variable) {
-        primary_variables_[WFrac] = F[pu.phase_pos[Water]];
+        primary_variables_.value_[WFrac] = F[pu.phase_pos[Water]];
     }
 
     if constexpr (has_gfrac_variable) {
-        primary_variables_[GFrac] = F[pu.phase_pos[Gas]];
+        primary_variables_.value_[GFrac] = F[pu.phase_pos[Gas]];
     }
     if constexpr (Indices::enableSolvent) {
-        primary_variables_[SFrac] = F_solvent;
+        primary_variables_.value_[SFrac] = F_solvent;
     }
 }
 
@@ -495,18 +486,18 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
 
         if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
             const int water_pos = pu.phase_pos[Water];
-            F[water_pos] = primary_variables_[WFrac];
+            F[water_pos] = primary_variables_.value_[WFrac];
             F[oil_pos] -= F[water_pos];
         }
 
         if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
             const int gas_pos = pu.phase_pos[Gas];
-            F[gas_pos] = primary_variables_[GFrac];
+            F[gas_pos] = primary_variables_.value_[GFrac];
             F[oil_pos] -= F[gas_pos];
         }
 
         if constexpr (Indices::enableSolvent) {
-            F_solvent = primary_variables_[SFrac];
+            F_solvent = primary_variables_.value_[SFrac];
             F[oil_pos] -= F_solvent;
         }
     }
@@ -516,7 +507,7 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
 
         if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
             const int gas_pos = pu.phase_pos[Gas];
-            F[gas_pos] = primary_variables_[GFrac];
+            F[gas_pos] = primary_variables_.value_[GFrac];
             F[water_pos] -= F[gas_pos];
         }
     }
@@ -546,12 +537,12 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
     }
 
     auto& ws = well_state.well(baseif_.indexOfWell());
-    ws.bhp = primary_variables_[Bhp];
+    ws.bhp = primary_variables_.value_[Bhp];
 
     // calculate the phase rates based on the primary variables
     // for producers, this is not a problem, while not sure for injectors here
     if (baseif_.isProducer()) {
-        const double g_total = primary_variables_[WQTotal];
+        const double g_total = primary_variables_.value_[WQTotal];
         for (int p = 0; p < baseif_.numPhases(); ++p) {
             ws.surface_rates[p] = g_total * F[p];
         }
@@ -561,13 +552,13 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
         }
         switch (baseif_.wellEcl().injectorType()) {
         case InjectorType::WATER:
-            ws.surface_rates[pu.phase_pos[Water]] = primary_variables_[WQTotal];
+            ws.surface_rates[pu.phase_pos[Water]] = primary_variables_.value_[WQTotal];
             break;
         case InjectorType::GAS:
-            ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_[WQTotal];
+            ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_.value_[WQTotal];
             break;
         case InjectorType::OIL:
-            ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_[WQTotal];
+            ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_.value_[WQTotal];
             break;
         case InjectorType::MULTI:
             // Not supported.
@@ -597,8 +588,8 @@ updateWellStateFromPrimaryVariablesPolyMW(WellState& well_state) const
         auto& perf_water_velocity = perf_data.water_velocity;
         auto& perf_skin_pressure = perf_data.skin_pressure;
         for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
-            perf_water_velocity[perf] = primary_variables_[Bhp + 1 + perf];
+            perf_water_velocity[perf] = primary_variables_.value_[Bhp + 1 + perf];
-            perf_skin_pressure[perf] = primary_variables_[Bhp + 1 + baseif_.numPerfs() + perf];
+            perf_skin_pressure[perf] = primary_variables_.value_[Bhp + 1 + baseif_.numPerfs() + perf];
         }
     }
 }
@@ -620,7 +611,7 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells,
                              [[maybe_unused]] const double dFLimit,
                              const double dBHPLimit) const
 {
-    const std::vector<double> old_primary_variables = primary_variables_;
+    const std::vector<double> old_primary_variables = primary_variables_.value_;
 
     // for injectors, very typical one of the fractions will be one, and it is easy to get zero value
     // fractions. not sure what is the best way to handle it yet, so we just use 1.0 here
@@ -634,33 +625,33 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells,
         const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
         const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
         // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
-        primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
+        primary_variables_.value_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
     }
 
     if constexpr (has_gfrac_variable) {
         const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
         const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
-        primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
+        primary_variables_.value_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
     }
 
     if constexpr (Indices::enableSolvent) {
         const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
         const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
-        primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
+        primary_variables_.value_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
     }
 
     processFractions();
 
     // updating the total rates Q_t
     const double relaxation_factor_rate = this->relaxationFactorRate(old_primary_variables, dwells[0][WQTotal]);
-    primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
+    primary_variables_.value_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
 
     // updating the bottom hole pressure
     {
         const int sign1 = dwells[0][Bhp] > 0 ? 1: -1;
         const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit);
         // 1e5 to make sure bhp will not be below 1bar
-        primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5);
+        primary_variables_.value_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5);
     }
 }
 
@@ -676,8 +667,8 @@ getWellConvergence(const WellState& well_state,
                    std::vector<double>& res,
                    DeferredLogger& deferred_logger) const
 {
-    res.resize(numWellEq_);
+    res.resize(this->primary_variables_.numWellEq());
-    for (int eq_idx = 0; eq_idx < numWellEq_; ++eq_idx) {
+    for (int eq_idx = 0; eq_idx < this->primary_variables_.numWellEq(); ++eq_idx) {
         // magnitude of the residual matters
         res[eq_idx] = std::abs(this->linSys_.residual()[0][eq_idx]);
     }
@@ -931,22 +922,21 @@ init(std::vector<double>& perf_depth,
     }
 
     // counting/updating primary variable numbers
+    int numWellEq = primary_variables_.numWellEq();
     if (has_polymermw) {
         if (baseif_.isInjector()) {
             // adding a primary variable for water perforation rate per connection
-            numWellEq_ += baseif_.numPerfs();
+            numWellEq += baseif_.numPerfs();
             // adding a primary variable for skin pressure per connection
-            numWellEq_ += baseif_.numPerfs();
+            numWellEq += baseif_.numPerfs();
         }
     }
 
     // with the updated numWellEq_, we can initialize the primary variables and matrices now
-    primary_variables_.resize(numWellEq_, 0.0);
+    primary_variables_.resize(numWellEq);
-    primary_variables_evaluation_.resize(numWellEq_, EvalWell{numWellEq_ + Indices::numEq, 0.0});
 
     // setup sparsity pattern for the matrices
-    this->linSys_.init(num_cells, this->numWellEq_,
+    this->linSys_.init(num_cells, numWellEq, baseif_.numPerfs(), baseif_.cells());
-                       baseif_.numPerfs(), baseif_.cells());
 }
 
 #define INSTANCE(...) \

opm/simulators/wells/StandardWellEval.hpp

@@ -25,6 +25,7 @@
 
 #include <opm/simulators/wells/StandardWellEquations.hpp>
 #include <opm/simulators/wells/StandardWellGeneric.hpp>
+#include <opm/simulators/wells/StandardWellPrimaryVariables.hpp>
 
 #include <opm/material/densead/DynamicEvaluation.hpp>
 
@@ -47,57 +48,19 @@ template<class FluidSystem, class Indices, class Scalar>
 class StandardWellEval : public StandardWellGeneric<Scalar>
 {
 protected:
-    // number of the conservation equations
+    using PrimaryVariables = StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>;
-    static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
+    static constexpr int Bhp = PrimaryVariables::Bhp;
-    // number of the well control equations
+    static constexpr int WQTotal= PrimaryVariables::WQTotal;
-    static constexpr int numWellControlEq = 1;
+    static constexpr int numWellConservationEq = PrimaryVariables::numWellConservationEq;
-    // number of the well equations that will always be used
+
-    // based on the solution strategy, there might be other well equations be introduced
+    static constexpr bool has_wfrac_variable = PrimaryVariables::has_wfrac_variable;
+    static constexpr bool has_gfrac_variable = PrimaryVariables::has_gfrac_variable;
+    static constexpr int WFrac = PrimaryVariables::WFrac;
+    static constexpr int GFrac = PrimaryVariables::GFrac;
+    static constexpr int SFrac = PrimaryVariables::SFrac;
 
 public:
-    static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
+    using EvalWell = typename PrimaryVariables::EvalWell;
-
-protected:
-    // the index for Bhp in primary variables and also the index of well control equation
-    // they both will be the last one in their respective system.
-    // TODO: we should have indices for the well equations and well primary variables separately
-    static constexpr int Bhp = numStaticWellEq - numWellControlEq;
-
-    // the positions of the primary variables for StandardWell
-    // the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g,
-    // which represent the fraction of Water and Gas based on the weighted total rate, the last one is BHP.
-    // correspondingly, we have four well equations for blackoil model, the first three are mass
-    // converstation equations, and the last one is the well control equation.
-    // primary variables related to other components, will be before the Bhp and after F_g.
-    // well control equation is always the last well equation.
-    // TODO: in the current implementation, we use the well rate as the first primary variables for injectors,
-    // instead of G_t.
-
-    // Table showing the primary variable indices, depending on what phases are present:
-    //
-    //         WOG     OG     WG     WO    W/O/G (single phase)
-    // WQTotal   0      0      0      0                       0
-    // WFrac     1  -1000     -1000   1                   -1000
-    // GFrac     2      1      1  -1000                   -1000
-    // Spres     3      2      2      2                       1
-
-    static const int WQTotal = 0; 
-    // In the implementation, one should use has_wfrac_variable
-    // rather than has_water to check if you should do something
-    // with the variable at the WFrac location, similar for GFrac.
-    // (following implementation MultisegmentWellEval.hpp)
-    static const bool waterEnabled = Indices::waterEnabled;
-    static const bool gasEnabled = Indices::gasEnabled;
-    static const bool oilEnabled = Indices::oilEnabled;
-
-    static constexpr bool has_wfrac_variable = Indices::waterEnabled && Indices::oilEnabled;
-    static constexpr bool has_gfrac_variable = Indices::gasEnabled && Indices::numPhases > 1;
-    static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
-    static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
-    static constexpr int SFrac = !Indices::enableSolvent ? -1000 : 3;
-
-public:
-    using EvalWell = DenseAd::DynamicEvaluation<Scalar, numStaticWellEq + Indices::numEq + 1>;
     using Eval = DenseAd::Evaluation<Scalar, Indices::numEq>;
     using BVectorWell = typename StandardWellGeneric<Scalar>::BVectorWell;
 
@@ -110,16 +73,14 @@ protected:
 
     const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif_;
 
-    void initPrimaryVariablesEvaluation() const;
-
     const EvalWell& getBhp() const
     {
-        return primary_variables_evaluation_[Bhp];
+        return primary_variables_.evaluation_[Bhp];
     }
 
     const EvalWell& getWQTotal() const
     {
-        return primary_variables_evaluation_[WQTotal];
+        return primary_variables_.evaluation_[WQTotal];
     }
 
     EvalWell extendEval(const Eval& in) const;
@@ -177,16 +138,7 @@ protected:
 
     void updateWellStateFromPrimaryVariablesPolyMW(WellState& well_state) const;
 
-    // total number of the well equations and primary variables
+    mutable PrimaryVariables primary_variables_; //!< Primary variables for well
-    // there might be extra equations be used, numWellEq will be updated during the initialization
-    int numWellEq_ = numStaticWellEq;
-
-    // the values for the primary varibles
-    // based on different solutioin strategies, the wells can have different primary variables
-    mutable std::vector<double> primary_variables_;
-
-    // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
-    mutable std::vector<EvalWell> primary_variables_evaluation_;
 
     // the saturations in the well bore under surface conditions at the beginning of the time step
     std::vector<double> F0_;

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -0,0 +1,88 @@
+/*
+  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/StandardWellPrimaryVariables.hpp>
+
+#include <opm/material/densead/Evaluation.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/WellInterfaceIndices.hpp>
+
+namespace Opm {
+
+template<class FluidSystem, class Indices, class Scalar>
+void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
+init()
+{
+    for (int eqIdx = 0; eqIdx < numWellEq_; ++eqIdx) {
+        evaluation_[eqIdx] =
+            EvalWell::createVariable(numWellEq_ + Indices::numEq,
+                                     value_[eqIdx],
+                                     Indices::numEq + eqIdx);
+
+    }
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
+resize(const int numWellEq)
+{
+    value_.resize(numWellEq, 0.0);
+    evaluation_.resize(numWellEq, EvalWell{numWellEq + Indices::numEq, 0.0});
+    numWellEq_ = numWellEq;
+}
+
+#define INSTANCE(...) \
+template class StandardWellPrimaryVariables<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,1u,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,1u,0u,false,true,0u,2u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>)
+INSTANCE(BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>)
+
+// Blackoil
+INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,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,0u,true,false,0u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,true,0u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,false,1u,0u>)
+
+}

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -0,0 +1,116 @@
+/*
+  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_PRIMARY_VARIABLES_HEADER_INCLUDED
+#define OPM_STANDARDWELL_PRIMARY_VARIABLES_HEADER_INCLUDED
+
+#include <opm/material/densead/DynamicEvaluation.hpp>
+
+#include <vector>
+
+namespace Opm
+{
+
+template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
+
+//! \brief Class holding primary variables for StandardWell.
+template<class FluidSystem, class Indices, class Scalar>
+class StandardWellPrimaryVariables {
+protected:
+    // the positions of the primary variables for StandardWell
+    // the first one is the weighted total rate (WQ_t), the second and the third ones are F_w and F_g,
+    // which represent the fraction of Water and Gas based on the weighted total rate, the last one is BHP.
+    // correspondingly, we have four well equations for blackoil model, the first three are mass
+    // converstation equations, and the last one is the well control equation.
+    // primary variables related to other components, will be before the Bhp and after F_g.
+    // well control equation is always the last well equation.
+    // TODO: in the current implementation, we use the well rate as the first primary variables for injectors,
+    // instead of G_t.
+
+    // Table showing the primary variable indices, depending on what phases are present:
+    //
+    //         WOG     OG     WG     WO    W/O/G (single phase)
+    // WQTotal   0      0      0      0                       0
+    // WFrac     1  -1000     -1000   1                   -1000
+    // GFrac     2      1      1  -1000                   -1000
+    // Spres     3      2      2      2                       1
+
+    //! \brief Number of the well control equations.
+    static constexpr int numWellControlEq = 1;
+
+public:
+    //! \brief Number of the conservation equations.
+    static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
+
+    //! \brief Number of the well equations that will always be used.
+    //! \details Based on the solution strategy, there might be other well equations be introduced.
+    static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
+
+    static const int WQTotal = 0; //!< The index for the weights total rate
+
+    //! \brief The index for Bhp in primary variables and the index of well control equation.
+    //! \details They both will be the last one in their respective system.
+    //! \todo: We should have indices for the well equations and well primary variables separately.
+    static constexpr int Bhp = numStaticWellEq - numWellControlEq;
+
+    static constexpr bool has_wfrac_variable = Indices::waterEnabled && Indices::oilEnabled;
+    static constexpr bool has_gfrac_variable = Indices::gasEnabled && Indices::numPhases > 1;
+    static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
+    static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
+    static constexpr int SFrac = !Indices::enableSolvent ? -1000 : 3;
+
+    //! \brief Evaluation for the well equations.
+    using EvalWell = DenseAd::DynamicEvaluation<Scalar, numStaticWellEq + Indices::numEq + 1>;
+
+    //! \brief Constructor initializes reference to well interface.
+    StandardWellPrimaryVariables(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well)
+        : well_(well)
+    {}
+
+    //! \brief The values for the primary variables.
+    //! \details Based on different solution strategies, the wells can have different primary variables.
+    std::vector<double> value_;
+
+    //! \brief The Evaluation for the well primary variables.
+    //! \details Contain derivatives and are used in AD calculation
+    std::vector<EvalWell> evaluation_;
+
+    //! \brief Initialize evaluations from values.
+    void init();
+
+    //! \brief Resize values and evaluations.
+    void resize(const int numWellEq);
+
+    //! \brief Returns number of well equations.
+    int numWellEq() const { return numWellEq_; }
+
+private:
+    const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well_; //!< Reference to well interface
+
+    //! \brief Total number of the well equations and primary variables.
+    //! \details There might be extra equations be used, numWellEq will be updated during the initialization
+    int numWellEq_ = numStaticWellEq;
+};
+
+}
+
+#endif // OPM_STANDARDWELL_PRIMARY_VARIABLES_HEADER_INCLUDED

opm/simulators/wells/StandardWell_impl.hpp

@@ -77,7 +77,7 @@ namespace Opm
     void StandardWell<TypeTag>::
     initPrimaryVariablesEvaluation() const
     {
-        this->StdWellEval::initPrimaryVariablesEvaluation();
+        this->primary_variables_.init();
     }
 
 
@@ -105,7 +105,7 @@ namespace Opm
         const EvalWell rv = this->extendEval(fs.Rv());
         const EvalWell rvw = this->extendEval(fs.Rvw());
 
-        std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->numWellEq_ + Indices::numEq, 0.0});
+        std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
         for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
                 continue;
@@ -125,16 +125,16 @@ namespace Opm
             }
         }
 
-        EvalWell skin_pressure = EvalWell{this->numWellEq_ + Indices::numEq, 0.0};
+        EvalWell skin_pressure = EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
         if (has_polymermw) {
             if (this->isInjector()) {
                 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
-                skin_pressure = this->primary_variables_evaluation_[pskin_index];
+                skin_pressure = this->primary_variables_.evaluation_[pskin_index];
             }
         }
 
         // surface volume fraction of fluids within wellbore
-        std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->numWellEq_ + Indices::numEq});
+        std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->primary_variables_.numWellEq() + Indices::numEq});
         for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
             cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx);
         }
@@ -199,7 +199,7 @@ namespace Opm
         if (has_polymermw) {
             if (this->isInjector()) {
                 const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
-                skin_pressure = getValue(this->primary_variables_evaluation_[pskin_index]);
+                skin_pressure = getValue(this->primary_variables_.evaluation_[pskin_index]);
             }
         }
 
@@ -465,9 +465,9 @@ namespace Opm
         auto& perf_rates = perf_data.phase_rates;
         for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
             // Calculate perforation quantities.
-            std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
+            std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
-            EvalWell water_flux_s{this->numWellEq_ + Indices::numEq, 0.0};
+            EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
-            EvalWell cq_s_zfrac_effective{this->numWellEq_ + Indices::numEq, 0.0};
+            EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
             calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
 
             // Equation assembly for this perforation.
@@ -487,7 +487,7 @@ namespace Opm
                     assemblePerforationEq(cq_s_effective,
                                           componentIdx,
                                           cell_idx,
-                                          this->numWellEq_,
+                                          this->primary_variables_.numWellEq(),
                                           this->linSys_);
 
                 // Store the perforation phase flux for later usage.
@@ -503,7 +503,7 @@ namespace Opm
                 StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
                     assembleZFracEq(cq_s_zfrac_effective,
                                     cell_idx,
-                                    this->numWellEq_,
+                                    this->primary_variables_.numWellEq(),
                                     this->linSys_);
             }
         }
@@ -526,7 +526,7 @@ namespace Opm
         for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
             // TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
             // since all the rates are under surface condition
-            EvalWell resWell_loc(this->numWellEq_ + Indices::numEq, 0.0);
+            EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
             if (FluidSystem::numActivePhases() > 1) {
                 assert(dt > 0);
                 resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt;
@@ -535,7 +535,7 @@ namespace Opm
             StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
                 assembleSourceEq(resWell_loc,
                                  componentIdx,
-                                 this->numWellEq_,
+                                 this->primary_variables_.numWellEq(),
                                  this->linSys_);
         }
 
@@ -545,7 +545,7 @@ namespace Opm
         StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
             assembleControlEq(well_state, group_state,
                               schedule, summaryState,
-                              this->numWellEq_,
+                              this->primary_variables_.numWellEq(),
                               this->getWQTotal(),
                               this->getBhp(),
                               gQ,
@@ -582,7 +582,7 @@ namespace Opm
         const EvalWell& bhp = this->getBhp();
         const int cell_idx = this->well_cells_[perf];
         const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
-        std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
+        std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
         getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
 
         double perf_dis_gas_rate = 0.;
@@ -627,7 +627,7 @@ namespace Opm
                 }
 
                 // convert to reservoir conditions
-                EvalWell cq_r_thermal(this->numWellEq_ + Indices::numEq, 0.);
+                EvalWell cq_r_thermal(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
                 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
                 const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
                 if ( !both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx ) {
@@ -917,7 +917,7 @@ namespace Opm
             // for the cases related to polymer molecular weight, we assume fully mixing
             // as a result, the polymer and water share the same viscosity
             if constexpr (!Base::has_polymermw) {
-                std::vector<EvalWell> mob_eval(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
+                std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
                 updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
                 for (size_t i = 0; i < mob.size(); ++i) {
                     mob[i] = getValue(mob_eval[i]);
@@ -960,7 +960,7 @@ namespace Opm
 
         updateExtraPrimaryVariables(dwells);
 
-        for (double v : this->primary_variables_) {
+        for (double v : this->primary_variables_.value_) {
             if(!isfinite(v))
                 OPM_DEFLOG_THROW(NumericalIssue, "Infinite primary variable after newton update well: " << this->name(),  deferred_logger);
         }
@@ -1535,7 +1535,7 @@ namespace Opm
                 return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
             };
 
-            std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
+            std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
             getMobilityEval(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
 
             const auto& fs = fluidState(subsetPerfID);
@@ -1684,7 +1684,7 @@ namespace Opm
         // We assemble the well equations, then we check the convergence,
         // which is why we do not put the assembleWellEq here.
         BVectorWell dx_well(1);
-        dx_well[0].resize(this->numWellEq_);
+        dx_well[0].resize(this->primary_variables_.numWellEq());
         this->linSys_.solve( dx_well);
 
         updateWellState(dx_well, well_state, deferred_logger);
@@ -1751,7 +1751,7 @@ namespace Opm
         if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         BVectorWell xw(1);
-        xw[0].resize(this->numWellEq_);
+        xw[0].resize(this->primary_variables_.numWellEq());
 
         this->linSys_.recoverSolutionWell(x, xw);
         updateWellState(xw, well_state, deferred_logger);
@@ -2040,12 +2040,12 @@ namespace Opm
                 const auto& water_velocity = perf_data.water_velocity;
                 const auto& skin_pressure = perf_data.skin_pressure;
                 for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
-                    this->primary_variables_[Bhp + 1 + perf] = water_velocity[perf];
+                    this->primary_variables_.value_[Bhp + 1 + perf] = water_velocity[perf];
-                    this->primary_variables_[Bhp + 1 + this->number_of_perforations_ + perf] = skin_pressure[perf];
+                    this->primary_variables_.value_[Bhp + 1 + this->number_of_perforations_ + perf] = skin_pressure[perf];
                 }
             }
         }
-        for (double v : this->primary_variables_) {
+        for (double v : this->primary_variables_.value_) {
             if(!isfinite(v))
                 OPM_DEFLOG_THROW(NumericalIssue, "Infinite primary variable after update from wellState well: " << this->name(),  deferred_logger);
         }
@@ -2097,7 +2097,7 @@ namespace Opm
             const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
             const EvalWell& bhp = this->getBhp();
 
-            std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
+            std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
             double perf_dis_gas_rate = 0.;
             double perf_vap_oil_rate = 0.;
             double perf_vap_wat_rate = 0.;
@@ -2171,9 +2171,9 @@ namespace Opm
                 OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger);
             }
             const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
-            const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
+            const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
             // the skin pressure when injecting water, which also means the polymer concentration is zero
-            EvalWell pskin_water(this->numWellEq_ + Indices::numEq, 0.0);
+            EvalWell pskin_water(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
             pskin_water = water_table_func.eval(throughput_eval, water_velocity);
             return pskin_water;
         } else {
@@ -2206,9 +2206,9 @@ namespace Opm
             }
             const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
             const double reference_concentration = skprpolytable.refConcentration;
-            const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
+            const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
             // the skin pressure when injecting water, which also means the polymer concentration is zero
-            EvalWell pskin_poly(this->numWellEq_ + Indices::numEq, 0.0);
+            EvalWell pskin_poly(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
             pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
             if (poly_inj_conc == reference_concentration) {
                 return sign * pskin_poly;
@@ -2237,8 +2237,8 @@ namespace Opm
         if constexpr (Base::has_polymermw) {
             const int table_id = this->well_ecl_.getPolymerProperties().m_plymwinjtable;
             const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
-            const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
+            const EvalWell throughput_eval(this->primary_variables_.numWellEq() + Indices::numEq, throughput);
-            EvalWell molecular_weight(this->numWellEq_ + Indices::numEq, 0.);
+            EvalWell molecular_weight(this->primary_variables_.numWellEq() + Indices::numEq, 0.);
             if (this->wpolymer() == 0.) { // not injecting polymer
                 return molecular_weight;
             }
@@ -2264,7 +2264,7 @@ namespace Opm
                 auto& ws = well_state.well(this->index_of_well_);
                 auto& perf_water_throughput = ws.perf_data.water_throughput;
                 for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
-                    const double perf_water_vel = this->primary_variables_[Bhp + 1 + perf];
+                    const double perf_water_vel = this->primary_variables_.value_[Bhp + 1 + perf];
                     // we do not consider the formation damage due to water flowing from reservoir into wellbore
                     if (perf_water_vel > 0.) {
                         perf_water_throughput[perf] += perf_water_vel * dt;
@@ -2295,7 +2295,7 @@ namespace Opm
 
         // water rate is update to use the form from water velocity, since water velocity is
         // a primary variable now
-        cq_s[water_comp_idx] = area * this->primary_variables_evaluation_[wat_vel_index] * b_w;
+        cq_s[water_comp_idx] = area * this->primary_variables_.evaluation_[wat_vel_index] * b_w;
     }
 
 
@@ -2320,7 +2320,7 @@ namespace Opm
         const int wat_vel_index = Bhp + 1 + perf;
 
         // equation for the water velocity
-        const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
+        const EvalWell eq_wat_vel = this->primary_variables_.evaluation_[wat_vel_index] - water_velocity;
 
         const auto& ws = well_state.well(this->index_of_well_);
         const auto& perf_data = ws.perf_data;
@@ -2328,12 +2328,12 @@ namespace Opm
         const double throughput = perf_water_throughput[perf];
         const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
 
-        EvalWell poly_conc(this->numWellEq_ + Indices::numEq, 0.0);
+        EvalWell poly_conc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
         poly_conc.setValue(this->wpolymer());
 
         // equation for the skin pressure
-        const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
+        const EvalWell eq_pskin = this->primary_variables_.evaluation_[pskin_index]
-                                  - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
+                                  - pskin(throughput, this->primary_variables_.evaluation_[wat_vel_index], poly_conc, deferred_logger);
 
         StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
                 assembleInjectivityEq(eq_pskin,
@@ -2341,7 +2341,7 @@ namespace Opm
                                       pskin_index,
                                       wat_vel_index,
                                       cell_idx,
-                                      this->numWellEq_,
+                                      this->primary_variables_.numWellEq(),
                                       this->linSys_);
     }
 
@@ -2382,7 +2382,7 @@ namespace Opm
         EvalWell cq_s_polymw = cq_s_poly;
         if (this->isInjector()) {
             const int wat_vel_index = Bhp + 1 + perf;
-            const EvalWell water_velocity = this->primary_variables_evaluation_[wat_vel_index];
+            const EvalWell water_velocity = this->primary_variables_.evaluation_[wat_vel_index];
             if (water_velocity > 0.) { // injecting
                 const auto& ws = well_state.well(this->index_of_well_);
                 const auto& perf_water_throughput = ws.perf_data.water_throughput;
@@ -2674,7 +2674,7 @@ namespace Opm
                              deferred_logger);
         }
 
-        const auto zero   = EvalWell { this->numWellEq_ + Indices::numEq, 0.0 };
+        const auto zero   = EvalWell{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
         const auto mt     = std::accumulate(mobility.begin(), mobility.end(), zero);
         connII[phase_pos] = connIICalc(mt.value() * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
     }