Make 2p gas oil compile and runs

-- use mapping from canonicalToActiveCompIdx from Ebos
-- add guards againts non-existing components

opm/autodiff/BlackoilModelEbos.hpp

@@ -635,7 +635,7 @@ namespace Opm {
 
                 PrimaryVariables& priVars = solution[ cell_idx ];
 
-                const double& dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)];
+                const double& dp = dx[cell_idx][Indices::pressureSwitchIdx];
                 double& p = priVars[Indices::pressureSwitchIdx];
                 const double& dp_rel_max = dpMaxRel();
                 const int sign_dp = dp > 0 ? 1: -1;
@@ -643,9 +643,9 @@ namespace Opm {
                 p = std::max(p, 0.0);
 
                 // Saturation updates.
-                const double dsw = active_[Water] ? dx[cell_idx][flowPhaseToEbosCompIdx(1)] : 0.0;
+                const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0;
                 const int xvar_ind = active_[Water] ? 2 : 1;
-                const double dxvar = active_[Gas] ? dx[cell_idx][flowPhaseToEbosCompIdx(xvar_ind)] : 0.0;
+                const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0;
 
                 double dso = 0.0;
                 double dsg = 0.0;
@@ -1533,51 +1533,36 @@ namespace Opm {
 
 
     public:
-        int ebosCompToFlowPhaseIdx( const int compIdx ) const
-        {
-            assert(compIdx < 3);
-            const int compToPhase[ 3 ] = { Oil, Water, Gas };
-            return compToPhase[ compIdx ];
-        }
-
-        int flowToEbosPvIdx( const int flowPv ) const
-        {
-            const int flowToEbos[] = {
-                Indices::pressureSwitchIdx,
-                Indices::waterSaturationIdx,
-                Indices::compositionSwitchIdx,
-                Indices::solventSaturationIdx
-            };
-
-            if (flowPv > 2 )
-                return flowPv;
-
-            return flowToEbos[ flowPv ];
-        }
-
         int flowPhaseToEbosCompIdx( const int phaseIdx ) const
         {
-            const int phaseToComp[] = {
-                FluidSystem::waterCompIdx,
-                FluidSystem::oilCompIdx,
-                FluidSystem::gasCompIdx
-            };
+            const auto& pu = phaseUsage_;
+            if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
+                return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+            if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
+                return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
+            if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
+                return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
 
-            if (phaseIdx > 2 )
-                return phaseIdx;
-
-            return phaseToComp[ phaseIdx ];
+            // for other phases return the index
+            return phaseIdx;
         }
 
-    private:
-
         int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
         {
+            const auto& pu = phaseUsage_;
+            if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
+                return FluidSystem::waterPhaseIdx;
+            if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
+                return FluidSystem::oilPhaseIdx;
+            if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
+                return FluidSystem::gasPhaseIdx;
+
             assert(phaseIdx < 3);
-            const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx};
-            return flowToEbos[ phaseIdx ];
+            // for other phases return the index
+            return phaseIdx;
         }
 
+    private:
 
         void updateRateConverter()
         {

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -900,7 +900,9 @@ protected:
             // set non-switching primary variables
             PrimaryVariables& cellPv = solution[ cellIdx ];
             // set water saturation
-            cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
+            if ( active[Water] ) {
+                cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
+            }
 
             if (has_solvent) {
                 cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
@@ -944,7 +946,9 @@ protected:
                             /*storeViscosity=*/false,
                             /*storeEnthalpy=*/false> SatOnlyFluidState;
                     SatOnlyFluidState fluidState;
-                    fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
+                    if ( active[Water] ) {
+                        fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
+                    }
                     fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
                     fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);
 

opm/autodiff/StandardWell.hpp

@@ -47,20 +47,21 @@ namespace Opm
         using typename Base::BlackoilIndices;
         using typename Base::PolymerModule;
 
+        using Base::numEq;
+
         // the positions of the primary variables for StandardWell
         // there are three primary variables, the second and the third ones are F_w and F_g
         // the first one can be total rate (G_t) or bhp, based on the control
-        enum WellVariablePositions {
-            XvarWell = 0,
-            WFrac = 1,
-            GFrac = 2,
-            SFrac = 3
-        };
+
+        static const bool gasoil = numEq == 2 && (BlackoilIndices::compositionSwitchIdx >= 0);
+        static const int XvarWell = 0;
+        static const int WFrac = gasoil? -1000: 1;
+        static const int GFrac = gasoil? 1: 2;
+        static const int SFrac = 3;
 
         using typename Base::Scalar;
         using typename Base::ConvergenceReport;
 
-        using Base::numEq;
 
         using Base::has_solvent;
         using Base::has_polymer;
@@ -296,6 +297,8 @@ namespace Opm
         void getMobility(const Simulator& ebosSimulator,
                          const int perf,
                          std::vector<EvalWell>& mob) const;
+
+        double scalingFactor(const int comp_idx) const;
     };
 
 }

opm/autodiff/StandardWell_impl.hpp

@@ -318,8 +318,7 @@ namespace Opm
             }
         }
 
-        std::vector<double> g = {1, 1, 0.01, 0.01};
-        return (wellVolumeFraction(compIdx) / g[compIdx]);
+        return (wellVolumeFraction(compIdx) / scalingFactor(compIdx) );
     }
 
 
@@ -331,11 +330,12 @@ namespace Opm
     StandardWell<TypeTag>::
     wellVolumeFraction(const int compIdx) const
     {
-        if (compIdx == Water) {
+        const auto pu = phaseUsage();
+        if (active()[Water] && compIdx == pu.phase_pos[Water]) {
             return primary_variables_evaluation_[WFrac];
         }
 
-        if (compIdx == Gas) {
+        if (active()[Gas] && compIdx == pu.phase_pos[Gas]) {
             return primary_variables_evaluation_[GFrac];
         }
 
@@ -390,7 +390,7 @@ namespace Opm
         EvalWell out = 0.0;
         out.setValue(in.value());
         for(int eqIdx = 0; eqIdx < numEq;++eqIdx) {
-            out.setDerivative(eqIdx, in.derivative(flowToEbosPvIdx(eqIdx)));
+            out.setDerivative(eqIdx, in.derivative(eqIdx));
         }
         return out;
     }
@@ -585,8 +585,8 @@ namespace Opm
                 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
                     if (!only_wells) {
                         // also need to consider the efficiency factor when manipulating the jacobians.
-                        ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
-                        duneB_[0][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx);
+                        ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s_effective.derivative(pvIdx);
+                        duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
                     }
                 }
 
@@ -613,7 +613,7 @@ namespace Opm
                 }
                 if (!only_wells) {
                     for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
-                        ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_poly.derivative(pvIdx);
+                        ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx);
                     }
                     ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value();
                 }
@@ -787,6 +787,7 @@ namespace Opm
         const int np = number_of_phases_;
         const double dBHPLimit = param.dbhp_max_rel_;
         const double dFLimit = param.dwell_fraction_max_;
+        const auto pu = phaseUsage();
 
         const std::vector<double> xvar_well_old = primary_variables_;
 
@@ -811,68 +812,68 @@ namespace Opm
         }
 
         assert(active()[ Oil ]);
-        F[Oil] = 1.0;
+        F[pu.phase_pos[Oil]] = 1.0;
 
         if (active()[ Water ]) {
-            F[Water] = primary_variables_[WFrac];
-            F[Oil] -= F[Water];
+            F[pu.phase_pos[Water]] = primary_variables_[WFrac];
+            F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
         }
 
         if (active()[ Gas ]) {
-            F[Gas] = primary_variables_[GFrac];
-            F[Oil] -= F[Gas];
+            F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
+            F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
         }
 
         double F_solvent = 0.0;
         if (has_solvent) {
             F_solvent = primary_variables_[SFrac];
-            F[Oil] -= F_solvent;
+            F[pu.phase_pos[Oil]] -= F_solvent;
         }
 
         if (active()[ Water ]) {
             if (F[Water] < 0.0) {
                 if (active()[ Gas ]) {
-                        F[Gas] /= (1.0 - F[Water]);
+                        F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]);
                 }
                 if (has_solvent) {
-                    F_solvent /= (1.0 - F[Water]);
+                    F_solvent /= (1.0 - F[pu.phase_pos[Water]]);
                 }
-                F[Oil] /= (1.0 - F[Water]);
-                F[Water] = 0.0;
+                F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]);
+                F[pu.phase_pos[Water]] = 0.0;
             }
         }
 
         if (active()[ Gas ]) {
-            if (F[Gas] < 0.0) {
+            if (F[pu.phase_pos[Gas]] < 0.0) {
                 if (active()[ Water ]) {
-                    F[Water] /= (1.0 - F[Gas]);
+                    F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]);
                 }
                 if (has_solvent) {
-                    F_solvent /= (1.0 - F[Gas]);
+                    F_solvent /= (1.0 - F[pu.phase_pos[Gas]]);
                 }
-                F[Oil] /= (1.0 - F[Gas]);
-                F[Gas] = 0.0;
+                F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]);
+                F[pu.phase_pos[Gas]] = 0.0;
             }
         }
 
-        if (F[Oil] < 0.0) {
+        if (F[pu.phase_pos[Oil]] < 0.0) {
             if (active()[ Water ]) {
-                F[Water] /= (1.0 - F[Oil]);
+                F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]);
             }
             if (active()[ Gas ]) {
-                F[Gas] /= (1.0 - F[Oil]);
+                F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]);
             }
             if (has_solvent) {
-                F_solvent /= (1.0 - F[Oil]);
+                F_solvent /= (1.0 - F[pu.phase_pos[Oil]]);
             }
-            F[Oil] = 0.0;
+            F[pu.phase_pos[Oil]] = 0.0;
         }
 
         if (active()[ Water ]) {
-            primary_variables_[WFrac] = F[Water];
+            primary_variables_[WFrac] = F[pu.phase_pos[Water]];
         }
         if (active()[ Gas ]) {
-            primary_variables_[GFrac] = F[Gas];
+            primary_variables_[GFrac] = F[pu.phase_pos[Gas]];
         }
         if(has_solvent) {
             primary_variables_[SFrac] = F_solvent;
@@ -881,7 +882,7 @@ namespace Opm
         // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
         // More testing is needed to make sure this is correct for well groups and THP.
         if (has_solvent){
-            F[Gas] += F_solvent;
+            F[pu.phase_pos[Gas]] += F_solvent;
         }
 
             // The interpretation of the first well variable depends on the well control
@@ -892,7 +893,6 @@ namespace Opm
         const int current = well_state.currentControls()[index_of_well_];
         const double target_rate = well_controls_iget_target(wc, current);
 
-        std::vector<double> g = {1,1,0.01};
         if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
             const double* distr = well_controls_iget_distr(wc, current);
             for (int p = 0; p < np; ++p) {
@@ -904,7 +904,7 @@ namespace Opm
             }
         } else {
             for (int p = 0; p < np; ++p) {
-                F[p] /= g[p];
+                F[p] /= scalingFactor(p);
             }
         }
 
@@ -1225,7 +1225,7 @@ namespace Opm
         const int w = index_of_well_;
 
         //rs and rv are only used if both oil and gas is present
-        if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) {
+        if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_used[BlackoilPhases::Liquid]) {
             rsmax_perf.resize(nperf);
             rvmax_perf.resize(nperf);
         }
@@ -1874,12 +1874,17 @@ namespace Opm
         const int well_index = index_of_well_;
         const WellControls* wc = well_controls_;
         const double* distr = well_controls_get_current_distr(wc);
+        const auto pu = phaseUsage();
 
-        std::vector<double> g = {1.0, 1.0, 0.01};
+        std::vector<double> g(np);
         if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
             for (int phase = 0; phase < np; ++phase) {
                 g[phase] = distr[phase];
             }
+        } else {
+            for (int phase = 0; phase < np; ++phase) {
+                g[phase] = scalingFactor(phase);
+            }
         }
 
         switch (well_controls_get_current_type(wc)) {
@@ -1909,13 +1914,13 @@ namespace Opm
         }
         if(std::abs(tot_well_rate) > 0) {
             if (active()[ Water ]) {
-                primary_variables_[WFrac] = g[Water] * well_state.wellRates()[np*well_index + Water] / tot_well_rate;
+                primary_variables_[WFrac] = g[pu.phase_pos[Water]] * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate;
             }
             if (active()[ Gas ]) {
-                primary_variables_[GFrac] = g[Gas] * (well_state.wellRates()[np*well_index + Gas] - well_state.solventWellRate(well_index)) / tot_well_rate ;
+                primary_variables_[GFrac] = g[ pu.phase_pos[Gas]] * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ;
             }
             if (has_solvent) {
-                primary_variables_[SFrac] = g[Gas] * well_state.solventWellRate(well_index) / tot_well_rate ;
+                primary_variables_[SFrac] = g[ pu.phase_pos[Gas]] * well_state.solventWellRate(well_index) / tot_well_rate ;
             }
         } else { // tot_well_rate == 0
             if (well_type_ == INJECTOR) {
@@ -1929,7 +1934,7 @@ namespace Opm
                 }
 
                 if (active()[Gas]) {
-                    if (distr[Gas] > 0.0) {
+                    if (distr[pu.phase_pos[Gas]] > 0.0) {
                         primary_variables_[GFrac] = 1.0 - wsolvent();
                         if (has_solvent) {
                             primary_variables_[SFrac] = wsolvent();
@@ -2053,4 +2058,25 @@ namespace Opm
          return thp;
     }
 
+    template<typename TypeTag>
+    double
+    StandardWell<TypeTag>::scalingFactor(const int phaseIdx) const
+    {
+        //std::vector<double> g = {1,1,0.01,0.01};
+        const auto& pu = phaseUsage();
+        if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
+            return 1.0;
+        if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
+            return 1.0;
+        if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
+            return 0.01;
+        if (has_solvent && phaseIdx == contiSolventEqIdx )
+            return 0.01;
+
+        // we should come this fare
+        assert(false);
+        return 1.0;
+    }
+
+
 }

opm/autodiff/StandardWellsDense_impl.hpp

@@ -302,9 +302,17 @@ namespace Opm {
     StandardWellsDense<TypeTag>::
     flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
     {
+        const auto& pu = phase_usage_;
+        if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
+            return FluidSystem::waterPhaseIdx;
+        if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
+            return FluidSystem::oilPhaseIdx;
+        if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
+            return FluidSystem::gasPhaseIdx;
+
         assert(phaseIdx < 3);
-        const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx };
-        return flowToEbos[ phaseIdx ];
+        // for other phases return the index
+        return phaseIdx;
     }
 
 

opm/autodiff/WellInterface_impl.hpp

@@ -224,47 +224,38 @@ namespace Opm
     WellInterface<TypeTag>::
     flowPhaseToEbosCompIdx( const int phaseIdx ) const
     {
-        const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx};
-        if (phaseIdx > 2 )
-            return phaseIdx;
-        return phaseToComp[ phaseIdx ];
+        const auto& pu = phaseUsage();
+        if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
+            return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
+        if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
+            return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
+        if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
+            return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+
+        // for other phases return the index
+        return phaseIdx;
     }
 
 
 
 
-
-    template<typename TypeTag>
-    int
-    WellInterface<TypeTag>::
-    flowToEbosPvIdx( const int flowPv ) const
-    {
-        const int flowToEbos[ 3 ] = {
-            BlackoilIndices::pressureSwitchIdx,
-            BlackoilIndices::waterSaturationIdx,
-            BlackoilIndices::compositionSwitchIdx
-        };
-
-        if (flowPv > 2 )
-            return flowPv;
-
-        return flowToEbos[ flowPv ];
-    }
-
-
-
-
-
     template<typename TypeTag>
     int
     WellInterface<TypeTag>::
     flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
     {
-        assert(phaseIdx < 3);
-        const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx };
-        return flowToEbos[ phaseIdx ];
-    }
+        const auto& pu = phaseUsage();
+        if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
+            return FluidSystem::waterPhaseIdx;
+        if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
+            return FluidSystem::oilPhaseIdx;
+        if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
+            return FluidSystem::gasPhaseIdx;
 
+        assert(phaseIdx < 3);
+        // for other phases return the index
+        return phaseIdx;
+    }