added: MultisegmentWellPrimaryVariables

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

CMakeLists_files.cmake

@@ -95,6 +95,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/MultisegmentWellEquations.cpp
   opm/simulators/wells/MultisegmentWellEval.cpp
   opm/simulators/wells/MultisegmentWellGeneric.cpp
+  opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
   opm/simulators/wells/ParallelWellInfo.cpp
   opm/simulators/wells/PerfData.cpp
   opm/simulators/wells/SegmentState.cpp
@@ -380,6 +381,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/MultisegmentWellEquations.hpp
   opm/simulators/wells/MultisegmentWellEval.hpp
   opm/simulators/wells/MultisegmentWellGeneric.hpp
+  opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
   opm/simulators/wells/ParallelWellInfo.hpp
   opm/simulators/wells/PerfData.hpp
   opm/simulators/wells/PerforationData.hpp

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -58,6 +58,7 @@ MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
     : MultisegmentWellGeneric<Scalar>(baseif)
     , baseif_(baseif)
     , linSys_(*this)
+    , primary_variables_(baseif)
     , upwinding_segments_(this->numberOfSegments(), 0)
     , segment_densities_(this->numberOfSegments(), 0.0)
     , segment_mass_rates_(this->numberOfSegments(), 0.0)
@@ -77,21 +78,6 @@ initMatrixAndVectors(const int num_cells)
 {
     linSys_.init(num_cells, baseif_.numPerfs(), baseif_.cells());
     primary_variables_.resize(this->numberOfSegments());
-    primary_variables_evaluation_.resize(this->numberOfSegments());
-}
-
-template<typename FluidSystem, typename Indices, typename Scalar>
-void
-MultisegmentWellEval<FluidSystem,Indices,Scalar>::
-initPrimaryVariablesEvaluation()
-{
-    for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
-        for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
-            primary_variables_evaluation_[seg][eq_idx] = 0.0;
-            primary_variables_evaluation_[seg][eq_idx].setValue(primary_variables_[seg][eq_idx]);
-            primary_variables_evaluation_[seg][eq_idx].setDerivative(eq_idx + Indices::numEq, 1.0);
-        }
-    }
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -204,13 +190,13 @@ processFractions(const int seg)
 
     if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
         const int water_pos = pu.phase_pos[Water];
-        fractions[water_pos] = primary_variables_[seg][WFrac];
+        fractions[water_pos] = primary_variables_.value_[seg][WFrac];
         fractions[oil_pos] -= fractions[water_pos];
     }
 
     if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
         const int gas_pos = pu.phase_pos[Gas];
-        fractions[gas_pos] = primary_variables_[seg][GFrac];
+        fractions[gas_pos] = primary_variables_.value_[seg][GFrac];
         fractions[oil_pos] -= fractions[gas_pos];
     }
 
@@ -247,11 +233,11 @@ processFractions(const int seg)
     }
 
     if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
-        primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]];
+        primary_variables_.value_[seg][WFrac] = fractions[pu.phase_pos[Water]];
     }
 
     if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
-        primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
+        primary_variables_.value_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
     }
 }
 
@@ -263,19 +249,19 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells,
                              const double dFLimit,
                              const double max_pressure_change)
 {
-    const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
+    const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_.value_;
 
     for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
         if (has_wfrac_variable) {
             const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
             const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit);
-            primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
+            primary_variables_.value_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
         }
 
         if (has_gfrac_variable) {
             const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
             const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit);
-            primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
+            primary_variables_.value_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
         }
 
         // handling the overshooting or undershooting of the fractions
@@ -285,19 +271,19 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells,
         {
             const int sign = dwells[seg][SPres] > 0.? 1 : -1;
             const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
-            primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5);
+            primary_variables_.value_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5);
         }
 
         // update the total rate // TODO: should we have a limitation of the total rate change?
         {
-            primary_variables_[seg][WQTotal] = old_primary_variables[seg][WQTotal] - relaxation_factor * dwells[seg][WQTotal];
+            primary_variables_.value_[seg][WQTotal] = old_primary_variables[seg][WQTotal] - relaxation_factor * dwells[seg][WQTotal];
 
             // make sure that no injector produce and no producer inject
             if (seg == 0) {
                 if (baseif_.isInjector()) {
-                    primary_variables_[seg][WQTotal] = std::max( primary_variables_[seg][WQTotal], 0.0);
+                    primary_variables_.value_[seg][WQTotal] = std::max( primary_variables_.value_[seg][WQTotal], 0.0);
                 } else {
-                    primary_variables_[seg][WQTotal] = std::min( primary_variables_[seg][WQTotal], 0.0);
+                    primary_variables_.value_[seg][WQTotal] = std::min( primary_variables_.value_[seg][WQTotal], 0.0);
                 }
             }
         }
@@ -331,7 +317,7 @@ updatePrimaryVariables(const WellState& well_state)
         // calculate the total rate for each segment
         double total_seg_rate = 0.0;
         // the segment pressure
-        primary_variables_[seg][SPres] = segment_pressure[seg];
+        primary_variables_.value_[seg][SPres] = segment_pressure[seg];
         // TODO: under what kind of circustances, the following will be wrong?
         // the definition of g makes the gas phase is always the last phase
         for (int p = 0; p < baseif_.numPhases(); p++) {
@@ -345,15 +331,15 @@ updatePrimaryVariables(const WellState& well_state)
                 total_seg_rate = std::min(total_seg_rate, 0.);
             }
         }
-        primary_variables_[seg][WQTotal] = total_seg_rate;
+        primary_variables_.value_[seg][WQTotal] = total_seg_rate;
         if (std::abs(total_seg_rate) > 0.) {
             if (has_wfrac_variable) {
                 const int water_pos = pu.phase_pos[Water];
-                primary_variables_[seg][WFrac] = baseif_.scalingFactor(water_pos) * segment_rates[baseif_.numPhases() * seg + water_pos] / total_seg_rate;
+                primary_variables_.value_[seg][WFrac] = baseif_.scalingFactor(water_pos) * segment_rates[baseif_.numPhases() * seg + water_pos] / total_seg_rate;
             }
             if (has_gfrac_variable) {
                 const int gas_pos = pu.phase_pos[Gas];
-                primary_variables_[seg][GFrac] = baseif_.scalingFactor(gas_pos) * segment_rates[baseif_.numPhases() * seg + gas_pos] / total_seg_rate;
+                primary_variables_.value_[seg][GFrac] = baseif_.scalingFactor(gas_pos) * segment_rates[baseif_.numPhases() * seg + gas_pos] / total_seg_rate;
             }
         } else { // total_seg_rate == 0
             if (baseif_.isInjector()) {
@@ -362,27 +348,27 @@ updatePrimaryVariables(const WellState& well_state)
 
                 if (has_wfrac_variable) {
                     if (phase == InjectorType::WATER) {
-                        primary_variables_[seg][WFrac] = 1.0;
+                        primary_variables_.value_[seg][WFrac] = 1.0;
                     } else {
-                        primary_variables_[seg][WFrac] = 0.0;
+                        primary_variables_.value_[seg][WFrac] = 0.0;
                     }
                 }
 
                 if (has_gfrac_variable) {
                     if (phase == InjectorType::GAS) {
-                        primary_variables_[seg][GFrac] = 1.0;
+                        primary_variables_.value_[seg][GFrac] = 1.0;
                     } else {
-                        primary_variables_[seg][GFrac] = 0.0;
+                        primary_variables_.value_[seg][GFrac] = 0.0;
                     }
                 }
 
             } else if (baseif_.isProducer()) { // producers
                 if (has_wfrac_variable) {
-                    primary_variables_[seg][WFrac] = 1.0 / baseif_.numPhases();
+                    primary_variables_.value_[seg][WFrac] = 1.0 / baseif_.numPhases();
                 }
 
                 if (has_gfrac_variable) {
-                    primary_variables_[seg][GFrac] = 1.0 / baseif_.numPhases();
+                    primary_variables_.value_[seg][GFrac] = 1.0 / baseif_.numPhases();
                 }
             }
         }
@@ -396,21 +382,21 @@ volumeFraction(const int seg,
                const unsigned compIdx) const
 {
     if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
-        return primary_variables_evaluation_[seg][WFrac];
+        return primary_variables_.evaluation_[seg][WFrac];
     }
 
     if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
-        return primary_variables_evaluation_[seg][GFrac];
+        return primary_variables_.evaluation_[seg][GFrac];
     }
 
     // Oil fraction
     EvalWell oil_fraction = 1.0;
     if (has_wfrac_variable) {
-        oil_fraction -= primary_variables_evaluation_[seg][WFrac];
+        oil_fraction -= primary_variables_.evaluation_[seg][WFrac];
     }
 
     if (has_gfrac_variable) {
-        oil_fraction -= primary_variables_evaluation_[seg][GFrac];
+        oil_fraction -= primary_variables_.evaluation_[seg][GFrac];
     }
     /* if (has_solvent) {
         oil_fraction -= primary_variables_evaluation_[seg][SFrac];
@@ -468,23 +454,23 @@ getSegmentRateUpwinding(const int seg,
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
                 && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx
                 && phase == InjectorType::WATER)
-            return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
+            return primary_variables_.evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
 
 
         if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
                 && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx
                 && phase == InjectorType::OIL)
-            return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
+            return primary_variables_.evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
 
         if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
                 && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx
                 && phase == InjectorType::GAS)
-            return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
+            return primary_variables_.evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
 
         return 0.0;
     }
 
-    const EvalWell segment_rate = primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg_upwind, comp_idx);
+    const EvalWell segment_rate = primary_variables_.evaluation_[seg][WQTotal] * volumeFractionScaled(seg_upwind, comp_idx);
 
     assert(segment_rate.derivative(SPres + Indices::numEq) == 0.);
 
@@ -682,7 +668,7 @@ typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 getSegmentPressure(const int seg) const
 {
-    return primary_variables_evaluation_[seg][SPres];
+    return primary_variables_.evaluation_[seg][SPres];
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -699,7 +685,7 @@ MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 getSegmentRate(const int seg,
                const int comp_idx) const
 {
-    return primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg, comp_idx);
+    return primary_variables_.evaluation_[seg][WQTotal] * volumeFractionScaled(seg, comp_idx);
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -715,7 +701,7 @@ typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 getSegmentWQTotal(const int seg) const
 {
-    return primary_variables_evaluation_[seg][WQTotal];
+    return primary_variables_.evaluation_[seg][WQTotal];
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -1223,13 +1209,13 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
 
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
             const int water_pos = pu.phase_pos[Water];
-            fractions[water_pos] = primary_variables_[seg][WFrac];
+            fractions[water_pos] = primary_variables_.value_[seg][WFrac];
             fractions[oil_pos] -= fractions[water_pos];
         }
 
         if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
             const int gas_pos = pu.phase_pos[Gas];
-            fractions[gas_pos] = primary_variables_[seg][GFrac];
+            fractions[gas_pos] = primary_variables_.value_[seg][GFrac];
             fractions[oil_pos] -= fractions[gas_pos];
         }
 
@@ -1246,7 +1232,7 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
         }
 
         // calculate the phase rates based on the primary variables
-        const double g_total = primary_variables_[seg][WQTotal];
+        const double g_total = primary_variables_.value_[seg][WQTotal];
         for (int p = 0; p < baseif_.numPhases(); ++p) {
             const double phase_rate = g_total * fractions[p];
             segment_rates[seg*baseif_.numPhases() + p] = phase_rate;
@@ -1256,7 +1242,7 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
         }
 
         // update the segment pressure
-        segment_pressure[seg] = primary_variables_[seg][SPres];
+        segment_pressure[seg] = primary_variables_.value_[seg][SPres];
 
         if (seg == 0) { // top segment
             ws.bhp = segment_pressure[seg];
@@ -1405,7 +1391,7 @@ assembleICDPressureEq(const int seg,
        (segment.segmentType() == Segment::SegmentType::VALVE) &&
        (segment.valve().status() == Opm::ICDStatus::SHUT) ) { // we use a zero rate equation to handle SHUT valve
         MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
-            assembleTrivialEq(seg, this->primary_variables_evaluation_[seg][WQTotal].value(), linSys_);
+            assembleTrivialEq(seg, this->primary_variables_.evaluation_[seg][WQTotal].value(), linSys_);
 
         auto& ws = well_state.well(baseif_.indexOfWell());
         ws.segments.pressure_drop_friction[seg] = 0.;
@@ -1629,14 +1615,14 @@ updateUpwindingSegments()
         // special treatment is needed for segment 0
         if (seg == 0) {
             // we are not supposed to have injecting producers and producing injectors
-            assert( ! (baseif_.isProducer() && primary_variables_evaluation_[seg][WQTotal] > 0.) );
+            assert( ! (baseif_.isProducer() && primary_variables_.evaluation_[seg][WQTotal] > 0.) );
-            assert( ! (baseif_.isInjector() && primary_variables_evaluation_[seg][WQTotal] < 0.) );
+            assert( ! (baseif_.isInjector() && primary_variables_.evaluation_[seg][WQTotal] < 0.) );
             upwinding_segments_[seg] = seg;
             continue;
         }
 
         // for other normal segments
-        if (primary_variables_evaluation_[seg][WQTotal] <= 0.) {
+        if (primary_variables_.evaluation_[seg][WQTotal] <= 0.) {
             upwinding_segments_[seg] = seg;
         } else {
             const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());

opm/simulators/wells/MultisegmentWellEval.hpp

@@ -24,6 +24,7 @@
 
 #include <opm/simulators/wells/MultisegmentWellEquations.hpp>
 #include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
+#include <opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp>
 
 #include <opm/material/densead/Evaluation.hpp>
 
@@ -103,7 +104,6 @@ protected:
     MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif);
 
     void initMatrixAndVectors(const int num_cells);
-    void initPrimaryVariablesEvaluation();
 
     void assembleDefaultPressureEq(const int seg,
                                    WellState& well_state);
@@ -217,12 +217,7 @@ protected:
 
     Equations linSys_; //!< The equation system
 
-    // the values for the primary varibles
+    MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar> primary_variables_;
-    // based on different solutioin strategies, the wells can have different primary variables
-    std::vector<std::array<double, numWellEq> > primary_variables_;
-
-    // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
-    std::vector<std::array<EvalWell, numWellEq> > primary_variables_evaluation_;
 
     // the upwinding segment for each segment based on the flow direction
     std::vector<int> upwinding_segments_;

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -0,0 +1,85 @@
+/*
+  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/MultisegmentWellPrimaryVariables.hpp>
+
+#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.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>
+
+namespace Opm {
+
+template<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
+resize(const int numSegments)
+{
+    value_.resize(numSegments);
+    evaluation_.resize(numSegments);
+}
+
+template<class FluidSystem, class Indices, class Scalar>
+void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
+init()
+{
+    for (size_t seg = 0; seg < value_.size(); ++seg) {
+        for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
+            evaluation_[seg][eq_idx] = 0.0;
+            evaluation_[seg][eq_idx].setValue(value_[seg][eq_idx]);
+            evaluation_[seg][eq_idx].setDerivative(eq_idx + Indices::numEq, 1.0);
+        }
+    }
+}
+
+#define INSTANCE(...) \
+template class MultisegmentWellPrimaryVariables<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,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,0u,false,false,1u,0u>)
+INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,2u,0u>)
+
+}

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -0,0 +1,95 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+
+  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_PRIMARY_VARIABLES_HEADER_INCLUDED
+#define OPM_MULTISEGMENTWELL_PRIMARY_VARIABLES_HEADER_INCLUDED
+
+#include <opm/material/densead/Evaluation.hpp>
+
+#include <array>
+#include <vector>
+
+namespace Opm
+{
+
+template<class FluidSystem, class Indices, class Scalar> class WellInterfaceIndices;
+
+template<class FluidSystem, class Indices, class Scalar>
+class MultisegmentWellPrimaryVariables
+{
+public:
+    // TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
+
+    // TODO: we need to have order for the primary variables and also the order for the well equations.
+    // sometimes, they are similar, while sometimes, they can have very different forms.
+
+    // 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      1      1                   -1000
+    // GFrac     2      1  -1000  -1000                   -1000
+    // Spres     3      2      2      2                       1
+
+    static constexpr bool has_water = (Indices::waterSwitchIdx >= 0);
+    static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0);
+    static constexpr bool has_oil = (Indices::numPhases - has_gas - has_water) > 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.
+    static constexpr bool has_wfrac_variable = has_water && Indices::numPhases > 1;
+    static constexpr bool has_gfrac_variable = has_gas && has_oil;
+
+    static constexpr int WQTotal = 0;
+    static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
+    static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
+    static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1;
+
+    //  the number of well equations  TODO: it should have a more general strategy for it
+    static constexpr int numWellEq = Indices::numPhases + 1;
+
+    using EvalWell = DenseAd::Evaluation<double, /*size=*/Indices::numEq + numWellEq>;
+
+    MultisegmentWellPrimaryVariables(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well)
+        : well_(well)
+    {}
+
+    //! \brief Resize values and evaluations.
+    void resize(const int numSegments);
+
+    //! \brief Initialize evaluations from values.
+    void init();
+
+    // the values for the primary varibles
+    // based on different solutioin strategies, the wells can have different primary variables
+    std::vector<std::array<double, numWellEq> > value_;
+
+    // the Evaluation for the well primary variables, which contain derivativles and are used in AD calculation
+    std::vector<std::array<EvalWell, numWellEq> > evaluation_;
+
+private:
+    const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well_; //!< Reference to well interface
+};
+
+}
+
+#endif // OPM_MULTISEGMENTWELL_PRIMARY_VARIABLES_HEADER_INCLUDED

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -130,7 +130,7 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     initPrimaryVariablesEvaluation()
     {
-        this->MSWEval::initPrimaryVariablesEvaluation();
+        this->primary_variables_.init();
     }
 
 
@@ -745,7 +745,7 @@ namespace Opm
                                                this->SPres,
                                                well_state);
     }
-    
+
 
     template<typename TypeTag>
     template<class Value>