fix well model for gasoil thermal

opm/simulators/wells/MultisegmentWell.hpp

@@ -48,6 +48,7 @@ namespace Opm
 
         /// the number of reservior equations
         using Base::numEq;
+        using Base::numPhases;
 
         using Base::has_solvent;
         using Base::has_polymer;
@@ -69,7 +70,7 @@ namespace Opm
         static constexpr int SPres = has_gas + has_water + 1;
 
         //  the number of well equations  TODO: it should have a more general strategy for it
-        static const int numWellEq = getPropValue<TypeTag, Properties::EnablePolymer>() ? numEq : numEq + 1;
+        static const int numWellEq = numPhases + 1;
 
         using typename Base::Scalar;
 

opm/simulators/wells/StandardWell.hpp

@@ -76,6 +76,7 @@ namespace Opm
         using GasLiftHandler = Opm::GasLiftRuntime<TypeTag>;
 
         using Base::numEq;
+        using Base::numPhases;
 
         using Base::has_solvent;
         using Base::has_zFraction;
@@ -88,16 +89,10 @@ namespace Opm
         using FoamModule = Opm::BlackOilFoamModule<TypeTag>;
         using BrineModule = Opm::BlackOilBrineModule<TypeTag>;
 
-        // polymer concentration and temperature are already known by the well, so
+        static const int numSolventEq = Indices::numSolvents;
-        // polymer and energy conservation do not need to be considered explicitly
-        static const int numPolymerEq = Indices::numPolymers;
-        static const int numEnergyEq = Indices::numEnergy;
-        static const int numFoamEq = Indices::numFoam;
-        static const int numBrineEq = Indices::numBrine;
-        static const int numExtbos = Indices::numExtbos;
 
         // number of the conservation equations
-        static const int numWellConservationEq = numEq - numPolymerEq - numEnergyEq - numFoamEq - numBrineEq - numExtbos;
+        static const int numWellConservationEq = numPhases + numSolventEq;
         // number of the well control equations
         static const int numWellControlEq = 1;
         // number of the well equations that will always be used
@@ -113,7 +108,7 @@ namespace Opm
         // 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.
-        static const bool gasoil = numEq == 2 && (Indices::compositionSwitchIdx >= 0);
+        static const bool gasoil = numPhases == 2 && (Indices::compositionSwitchIdx >= 0);
         static const int WQTotal = 0;
         static const int WFrac = gasoil? -1000: 1;
         static const int GFrac = gasoil? 1: 2;
@@ -331,6 +326,7 @@ namespace Opm
         // protected functions from the Base class
         using Base::getAllowCrossFlow;
         using Base::flowPhaseToEbosCompIdx;
+        using Base::flowPhaseToEbosPhaseIdx;
         using Base::ebosCompIdxToFlowCompIdx;
         using Base::wsalt;
         using Base::wsolvent;

opm/simulators/wells/StandardWell_impl.hpp

@@ -4170,7 +4170,8 @@ namespace Opm
             // Note: E100's notion of PI value phase mobility includes
             // the reciprocal FVF.
             const auto connMob =
-                mobility[ flowPhaseToEbosCompIdx(p) ].value() * fs.invB(p).value();
+                mobility[ flowPhaseToEbosCompIdx(p) ].value()
+                    * fs.invB(flowPhaseToEbosPhaseIdx(p)).value();
 
             connPI[p] = connPICalc(connMob);
         }
@@ -4227,6 +4228,6 @@ namespace Opm
 
         const auto zero   = EvalWell { this->numWellEq_ + this->numEq, 0.0 };
         const auto mt     = std::accumulate(mobility.begin(), mobility.end(), zero);
-        connII[phase_pos] = connIICalc(mt.value() * fs.invB(phase_pos).value());
+        connII[phase_pos] = connIICalc(mt.value() * fs.invB(flowPhaseToEbosPhaseIdx(phase_pos)).value());
     }
 } // namespace Opm

opm/simulators/wells/WellInterface.hpp

@@ -85,6 +85,7 @@ namespace Opm
         using RateVector = GetPropType<TypeTag, Properties::RateVector>;
 
         static const int numEq = Indices::numEq;
+        static const int numPhases = Indices::numPhases;
         typedef double Scalar;
 
         typedef Dune::FieldVector<Scalar, numEq    > VectorBlockType;
@@ -434,6 +435,8 @@ namespace Opm
 
         int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
 
+        int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
+
         int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const;
 
         double wsolvent() const;

opm/simulators/wells/WellInterface_impl.hpp

@@ -297,6 +297,23 @@ namespace Opm
         return phaseIdx;
     }
 
+    template<typename TypeTag>
+    int
+    WellInterface<TypeTag>::
+    flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
+    {
+        const auto& pu = phaseUsage();
+        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && pu.phase_pos[Water] == phaseIdx)
+            return FluidSystem::waterPhaseIdx;
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && pu.phase_pos[Oil] == phaseIdx)
+            return FluidSystem::oilPhaseIdx;
+        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx)
+            return FluidSystem::gasPhaseIdx;
+
+        // for other phases return the index
+        return phaseIdx;
+    }
+
     template<typename TypeTag>
     int
     WellInterface<TypeTag>::