Remove blackoilPhases and phaseUsage from the initialization code

Note 1: The initialization code now always consider 3 phases.
For 2-phase cases a trivial (0) state is returned.
Note 2: The initialization code does not compute a BlackoilStats,
but instead pass the initialization object with the initial state.

examples/compute_initial_state.cpp

@@ -28,6 +28,8 @@
 #include <opm/core/simulator/initStateEquil.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/core/props/phaseUsageFromDeck.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/utility/compressedToCartesian.hpp>
 
@@ -75,12 +77,45 @@ namespace
         std::copy(data.begin(), data.end(), std::ostream_iterator<double>(file, "\n"));
     }
 
+    /// Convert saturations from a vector of individual phase saturation vectors
+    /// to an interleaved format where all values for a given cell come before all
+    /// values for the next cell, all in a single vector.
+    template <class FluidSystem>
+    void convertSats(std::vector<double>& sat_interleaved, const std::vector< std::vector<double> >& sat, const Opm::PhaseUsage& pu)
+    {
+        assert(sat.size() == 3);
+        const auto nc = sat[0].size();
+        const auto np = sat_interleaved.size() / nc;
+        for (size_t c = 0; c < nc; ++c) {
+            if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                const int opos = pu.phase_pos[Opm::BlackoilPhases::Liquid];
+                const std::vector<double>& sat_p = sat[ FluidSystem::oilPhaseIdx];
+                sat_interleaved[np*c + opos] = sat_p[c];
+            }
+            if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                const int wpos = pu.phase_pos[Opm::BlackoilPhases::Aqua];
+                const std::vector<double>& sat_p = sat[ FluidSystem::waterPhaseIdx];
+                sat_interleaved[np*c + wpos] = sat_p[c];
+            }
+            if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                const int gpos = pu.phase_pos[Opm::BlackoilPhases::Vapour];
+                const std::vector<double>& sat_p = sat[ FluidSystem::gasPhaseIdx];
+                sat_interleaved[np*c + gpos] = sat_p[c];
+            }
+        }
+    }
+
 
 
 } // anon namespace
 
 
 
+
+
+
+
+
 // ----------------- Main program -----------------
 int
 main(int argc, char** argv)
@@ -119,7 +154,32 @@ try
     // Initialisation.
     //initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
     BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3);
-    initStateEquil(grid, materialLawManager, deck, eclipseState, grav, state);
+    typedef FluidSystems::BlackOil<double> FluidSystem;
+    FluidSystem::initFromDeck(deck, eclipseState);
+    PhaseUsage pu = phaseUsageFromDeck(deck);
+
+    typedef EQUIL::DeckDependent::InitialStateComputer<FluidSystem> ISC;
+
+    ISC isc(materialLawManager, eclipseState, grid, grav);
+
+    const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
+    const int oilpos = FluidSystem::oilPhaseIdx;
+    const int waterpos = FluidSystem::waterPhaseIdx;
+    const int ref_phase = oil ? oilpos : waterpos;
+
+    state.pressure() = isc.press()[ref_phase];
+    convertSats<FluidSystem>(state.saturation(), isc.saturation(), pu);
+
+    if (state.hasCellData(std::string("GASOILRATIO"))) {
+        std::vector<double>& rs = state.getCellData(std::string("GASOILRATIO"));
+        rs = isc.rs();
+    }
+    if (state.hasCellData(std::string("RV"))){
+        std::vector<double>& rv = state.getCellData(std::string("RV"));
+        rv = isc.rv();
+    }
+
+
 
     // Output.
     const std::string output_dir = param.getDefault<std::string>("output_dir", "output");
@@ -132,3 +192,5 @@ catch (const std::exception& e) {
     std::cerr << "Program threw an exception: " << e.what() << "\n";
     throw;
 }
+
+

opm/core/simulator/EquilibrationHelpers.hpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2017 IRIS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -20,7 +21,6 @@
 #ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
 #define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
 
-#include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/utility/linearInterpolation.hpp>
 #include <opm/core/utility/RegionMapping.hpp>
 #include <opm/core/utility/RootFinders.hpp>
@@ -478,13 +478,11 @@ namespace Opm
             EquilReg(const EquilRecord& rec,
                      std::shared_ptr<Miscibility::RsFunction> rs,
                      std::shared_ptr<Miscibility::RsFunction> rv,
-                     const int pvtIdx,
+                     const int pvtIdx)
-                     const PhaseUsage&  pu)
                 : rec_    (rec)
                 , rs_     (rs)
                 , rv_     (rv)
                 , pvtIdx_ (pvtIdx)
-                , pu_     (pu)
 
             {
             }
@@ -554,18 +552,12 @@ namespace Opm
             const int
             pvtIdx() const { return this->pvtIdx_; }
 
-            /**
-             * Retrieve active fluid phase summary.
-             */
-            const PhaseUsage&
-            phaseUsage() const { return this->pu_; }
 
         private:
             EquilRecord rec_;     /**< Equilibration data */
             std::shared_ptr<Miscibility::RsFunction> rs_;      /**< RS calculator */
             std::shared_ptr<Miscibility::RsFunction> rv_;      /**< RV calculator */
             const int pvtIdx_;
-            PhaseUsage  pu_;      /**< Active phase summary */
         };
 
 
@@ -588,7 +580,7 @@ namespace Opm
                 fluidState_.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
                 fluidState_.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
                 fluidState_.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
-                std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
+                std::fill(pc_, pc_ + FluidSystem::numPhases, 0.0);
 
             }
 
@@ -609,7 +601,7 @@ namespace Opm
             const int cell_;
             const double target_pc_;
             mutable SatOnlyFluidState fluidState_;
-            mutable double pc_[BlackoilPhases::MaxNumPhases];
+            mutable double pc_[FluidSystem::numPhases];
         };
 
         template <class FluidSystem, class MaterialLawManager>
@@ -722,7 +714,7 @@ namespace Opm
                 fluidState_.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
                 fluidState_.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
                 fluidState_.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
-                std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
+                std::fill(pc_, pc_ + FluidSystem::numPhases, 0.0);
             }
             double operator()(double s) const
             {
@@ -745,7 +737,7 @@ namespace Opm
             const int cell_;
             const double target_pc_;
             mutable SatOnlyFluidState fluidState_;
-            mutable double pc_[BlackoilPhases::MaxNumPhases];
+            mutable double pc_[FluidSystem::numPhases];
         };
 
 

opm/core/simulator/initStateEquil.hpp

@@ -24,9 +24,6 @@
 
 #include <opm/core/grid/GridHelpers.hpp>
 #include <opm/core/simulator/EquilibrationHelpers.hpp>
-#include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/props/BlackoilPhases.hpp>
-#include <opm/core/props/phaseUsageFromDeck.hpp>
 #include <opm/core/utility/RegionMapping.hpp>
 #include <opm/core/utility/extractPvtTableIndex.hpp>
 
@@ -40,6 +37,8 @@
 #include <opm/parser/eclipse/EclipseState/Tables/RsvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/data/SimulationDataContainer.hpp>
+
 
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 #include <opm/material/fluidstates/SimpleModularFluidState.hpp>
@@ -198,12 +197,12 @@ namespace Opm
             template<class Grid>
             inline
             std::vector<int>
-            equilnum(const Opm::Deck& deck,
+            equilnum(const Opm::EclipseState& eclipseState,
-                     const Opm::EclipseState& eclipseState,
                      const Grid&  G   )
             {
                 std::vector<int> eqlnum;
-                if (deck.hasKeyword("EQLNUM")) {
+
+                if (eclipseState.get3DProperties().hasDeckIntGridProperty("EQLNUM")) {
                     const int nc = UgGridHelpers::numCells(G);
                     eqlnum.resize(nc);
                     const std::vector<int>& e = 
@@ -223,45 +222,48 @@ namespace Opm
                 return eqlnum;
             }
 
-
+            template<class FluidSystem>
             class InitialStateComputer {
             public:
                 template<class MaterialLawManager, class Grid>
                 InitialStateComputer(std::shared_ptr<MaterialLawManager> materialLawManager,
-                                     const PhaseUsage& phaseUsage,
-                                     const Opm::Deck&             deck,
                                      const Opm::EclipseState& eclipseState,
                                      const Grid&                        G    ,
                                      const double grav = unit::gravity,
-                                     const std::vector<double>&  swat_init = {}
+                                     const bool applySwatInit = true
                                      )
-                    : pp_(phaseUsage.num_phases,
+                    : pp_(FluidSystem::numPhases,
                           std::vector<double>(UgGridHelpers::numCells(G))),
-                      sat_(phaseUsage.num_phases,
+                      sat_(FluidSystem::numPhases,
                           std::vector<double>(UgGridHelpers::numCells(G))),
                       rs_(UgGridHelpers::numCells(G)),
-                      rv_(UgGridHelpers::numCells(G)),
+                      rv_(UgGridHelpers::numCells(G))
-                      swat_init_(swat_init),
-                      phaseUsage_(phaseUsage)
-
                 {
 
-                    typedef FluidSystems::BlackOil<double> FluidSystem;
+                    //Check for presence of kw SWATINIT
-
+                    if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatInit) {
-                    // Initialize the fluid system
+                        const std::vector<double>& swat_init_ecl = eclipseState.
-                    FluidSystem::initFromDeck(deck, eclipseState);
+                                get3DProperties().getDoubleGridProperty("SWATINIT").getData();
+                        const int nc = UgGridHelpers::numCells(G);
+                        swat_init_.resize(nc);
+                        const int* gc = UgGridHelpers::globalCell(G);
+                        for (int c = 0; c < nc; ++c) {
+                            const int deck_pos = (gc == NULL) ? c : gc[c];
+                            swat_init_[c] = swat_init_ecl[deck_pos];
+                        }
+                    }
 
                     // Get the equilibration records.
                     const std::vector<EquilRecord> rec = getEquil(eclipseState);
                     const auto& tables = eclipseState.getTableManager();
                     // Create (inverse) region mapping.
-                    const RegionMapping<> eqlmap(equilnum(deck, eclipseState, G)); 
+                    const RegionMapping<> eqlmap(equilnum(eclipseState, G));
 
                     setRegionPvtIdx(G, eclipseState, eqlmap);
 
                     // Create Rs functions.
                     rs_func_.reserve(rec.size());
-                    if (deck.hasKeyword("DISGAS")) {
+                    if (FluidSystem::enableDissolvedGas()) {
                         const TableContainer& rsvdTables = tables.getRsvdTables();
                         for (size_t i = 0; i < rec.size(); ++i) {
                             if (eqlmap.cells(i).empty())
@@ -299,7 +301,7 @@ namespace Opm
                     }                    
 
                     rv_func_.reserve(rec.size());
-                    if (deck.hasKeyword("VAPOIL")) {
+                    if (FluidSystem::enableVaporizedOil()) {
                         const TableContainer& rvvdTables = tables.getRvvdTables();
                         for (size_t i = 0; i < rec.size(); ++i) {
                             if (eqlmap.cells(i).empty())
@@ -338,7 +340,7 @@ namespace Opm
                     }
                     
                     // Compute pressures, saturations, rs and rv factors.
-                    calcPressSatRsRv<FluidSystem>(eqlmap, rec, materialLawManager, G, grav);
+                    calcPressSatRsRv(eqlmap, rec, materialLawManager, G, grav);
 
                     // Modify oil pressure in no-oil regions so that the pressures of present phases can
                     // be recovered from the oil pressure and capillary relations.
@@ -366,7 +368,6 @@ namespace Opm
                 Vec rs_;
                 Vec rv_;
                 Vec swat_init_;
-                PhaseUsage phaseUsage_;
 
                 template<class Grid, class RMap>
                 void setRegionPvtIdx(const Grid& G, const Opm::EclipseState& eclipseState, const RMap& reg) {
@@ -381,7 +382,7 @@ namespace Opm
                     }
                 }
 
-                template <class FluidSystem, class RMap, class MaterialLawManager, class Grid>
+                template <class RMap, class MaterialLawManager, class Grid>
                 void
                 calcPressSatRsRv(const RMap&                       reg  ,
                                  const std::vector< EquilRecord >& rec  ,
@@ -398,23 +399,24 @@ namespace Opm
                             continue;
                         }
 
-                        const EqReg eqreg(rec[r],
+                        const EqReg eqreg(rec[r], rs_func_[r], rv_func_[r], regionPvtIdx_[r]);
-                                          rs_func_[r], rv_func_[r], regionPvtIdx_[r],
-                                          phaseUsage_);
                    
                         PVec pressures = phasePressures<FluidSystem>(G, eqreg, cells, grav);
                         const std::vector<double>& temp = temperature(G, eqreg, cells);
                         const PVec sat = phaseSaturations<FluidSystem>(G, eqreg, cells, materialLawManager, swat_init_, pressures);
 
-                        const int np = phaseUsage_.num_phases;
+                        const int np = FluidSystem::numPhases;
                         for (int p = 0; p < np; ++p) {
                             copyFromRegion(pressures[p], cells, pp_[p]);
                             copyFromRegion(sat[p], cells, sat_[p]);
                         }
-                        if (phaseUsage_.phase_used[BlackoilPhases::Liquid]
+
-                            && phaseUsage_.phase_used[BlackoilPhases::Vapour]) {
+                        const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
-                            const int oilpos = phaseUsage_.phase_pos[BlackoilPhases::Liquid];
+                        const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
-                            const int gaspos = phaseUsage_.phase_pos[BlackoilPhases::Vapour];
+
+                        if (oil && gas) {
+                            const int oilpos = FluidSystem::oilPhaseIdx;
+                            const int gaspos = FluidSystem::gasPhaseIdx;
                             const Vec rs_vals = computeRs(G, cells, pressures[oilpos], temp, *(rs_func_[r]), sat[gaspos]);
                             const Vec rv_vals = computeRs(G, cells, pressures[gaspos], temp, *(rv_func_[r]), sat[oilpos]);
                             copyFromRegion(rs_vals, cells, rs_);

opm/core/simulator/initStateEquil_impl.hpp

@@ -25,7 +25,7 @@
 #include <opm/core/grid.h>
 #include <opm/core/grid/GridHelpers.hpp>
 
-#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 
 #include <cassert>
 #include <cmath>
@@ -237,60 +237,6 @@ namespace Opm
             };
         } // namespace PhasePressODE
 
-        namespace PhaseUsed {
-            inline bool
-            water(const PhaseUsage& pu)
-            {
-                return bool(pu.phase_used[ Opm::BlackoilPhases::Aqua ]);
-            }
-
-            inline bool
-            oil(const PhaseUsage& pu)
-            {
-                return bool(pu.phase_used[ Opm::BlackoilPhases::Liquid ]);
-            }
-
-            inline bool
-            gas(const PhaseUsage& pu)
-            {
-                return bool(pu.phase_used[ Opm::BlackoilPhases::Vapour ]);
-            }
-        } // namespace PhaseUsed
-
-        namespace PhaseIndex {
-            inline int
-            water(const PhaseUsage& pu)
-            {
-                int i = -1;
-                if (PhaseUsed::water(pu)) {
-                    i = pu.phase_pos[ Opm::BlackoilPhases::Aqua ];
-                }
-
-                return i;
-            }
-
-            inline int
-            oil(const PhaseUsage& pu)
-            {
-                int i = -1;
-                if (PhaseUsed::oil(pu)) {
-                    i = pu.phase_pos[ Opm::BlackoilPhases::Liquid ];
-                }
-
-                return i;
-            }
-
-            inline int
-            gas(const PhaseUsage& pu)
-            {
-                int i = -1;
-                if (PhaseUsed::gas(pu)) {
-                    i = pu.phase_pos[ Opm::BlackoilPhases::Vapour ];
-                }
-
-                return i;
-            }
-        } // namespace PhaseIndex
 
         namespace PhasePressure {
             template <class Grid,
@@ -491,70 +437,66 @@ namespace Opm
                        const CellRange&                    cells,
                        std::vector< std::vector<double> >& press)
         {
-            const PhaseUsage& pu = reg.phaseUsage();
+            const bool water = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
+            const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
+            const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
+            const int oilpos = FluidSystem::oilPhaseIdx;
+            const int waterpos = FluidSystem::waterPhaseIdx;
+            const int gaspos = FluidSystem::gasPhaseIdx;
 
             if (reg.datum() > reg.zwoc()) { // Datum in water zone
                 double po_woc = -1;
                 double po_goc = -1;
 
-                if (PhaseUsed::water(pu)) {
+                if (water) {
-                    const int wix = PhaseIndex::water(pu);
                     PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
-                                         cells, press[ wix ]);
+                                         cells, press[ waterpos ]);
                 }
 
-                if (PhaseUsed::oil(pu)) {
+                if (oil) {
-                    const int oix = PhaseIndex::oil(pu);
                     PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
-                                       press[ oix ], po_woc, po_goc);
+                                       press[ oilpos ], po_woc, po_goc);
                 }
 
-                if (PhaseUsed::gas(pu)) {
+                if (gas) {
-                    const int gix = PhaseIndex::gas(pu);
                     PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
-                                       cells, press[ gix ]);
+                                       cells, press[ gaspos ]);
                 }
             } else if (reg.datum() < reg.zgoc()) { // Datum in gas zone
                 double po_woc = -1;
                 double po_goc = -1;
 
-                if (PhaseUsed::gas(pu)) {
+                if (gas) {
-                    const int gix = PhaseIndex::gas(pu);
                     PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
-                                       cells, press[ gix ]);
+                                       cells, press[ gaspos ]);
                 }
 
-                if (PhaseUsed::oil(pu)) {
+                if (oil) {
-                    const int oix = PhaseIndex::oil(pu);
                     PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
-                                       press[ oix ], po_woc, po_goc);
+                                       press[ oilpos ], po_woc, po_goc);
                 }
 
-                if (PhaseUsed::water(pu)) {
+                if (water) {
-                    const int wix = PhaseIndex::water(pu);
                     PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
-                                         cells, press[ wix ]);
+                                         cells, press[ waterpos ]);
                 }
             } else { // Datum in oil zone
                 double po_woc = -1;
                 double po_goc = -1;
 
-                if (PhaseUsed::oil(pu)) {
+                if (oil) {
-                    const int oix = PhaseIndex::oil(pu);
                     PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
-                                                    press[ oix ], po_woc, po_goc);
+                                                    press[ oilpos ], po_woc, po_goc);
                 }
 
-                if (PhaseUsed::water(pu)) {
+                if (water) {
-                    const int wix = PhaseIndex::water(pu);
                     PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
-                                         cells, press[ wix ]);
+                                         cells, press[ waterpos ]);
                 }
 
-                if (PhaseUsed::gas(pu)) {
+                if (gas) {
-                    const int gix = PhaseIndex::gas(pu);
                     PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
-                                       cells, press[ gix ]);
+                                       cells, press[ gaspos ]);
                 }
             }
         }
@@ -621,7 +563,7 @@ namespace Opm
                     }
                 }
             }
-            const int np = reg.phaseUsage().num_phases;
+            const int np = FluidSystem::numPhases;  //reg.phaseUsage().num_phases;
 
             typedef std::vector<double> pval;
             std::vector<pval> press(np, pval(ncell, 0.0));
@@ -659,7 +601,7 @@ namespace Opm
                          const std::vector<double> swat_init,
                          std::vector< std::vector<double> >& phase_pressures)
         {
-            if (!reg.phaseUsage().phase_used[BlackoilPhases::Liquid]) {
+            if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                 OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases.");
             }
 
@@ -682,11 +624,11 @@ namespace Opm
             SatOnlyFluidState fluidState;
             typedef typename MaterialLawManager::MaterialLaw MaterialLaw;
 
-            const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
+            const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
-            const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
+            const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
-            const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+            const int oilpos = FluidSystem::oilPhaseIdx;
-            const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua];
+            const int waterpos = FluidSystem::waterPhaseIdx;
-            const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+            const int gaspos = FluidSystem::gasPhaseIdx;
             std::vector<double>::size_type local_index = 0;
             for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
                 const int cell = *ci;
@@ -769,7 +711,7 @@ namespace Opm
                 double threshold_sat = 1.0e-6;
 
                 double so = 1.0;
-                double pC[/*numPhases=*/BlackoilPhases::MaxNumPhases] = { 0.0, 0.0, 0.0 };
+                double pC[FluidSystem::numPhases] = { 0.0, 0.0, 0.0 };
 
                   if (water) {
                       double swu = scaledDrainageInfo.Swu;
@@ -850,94 +792,6 @@ namespace Opm
     } // namespace Equil
 
 
-    namespace Details
-    {
-        /// Convert saturations from a vector of individual phase saturation vectors
-        /// to an interleaved format where all values for a given cell come before all
-        /// values for the next cell, all in a single vector.
-        inline std::vector<double>
-        convertSats(const std::vector< std::vector<double> >& sat)
-        {
-            const auto np = sat.size();
-            const auto nc = sat[0].size();
-
-            std::vector<double> s(np * nc);
-
-            for (decltype(sat.size()) p = 0; p < np; ++p) {
-                const auto& sat_p = sat[p];
-                double*     sp    = & s[0*nc + p];
-
-                for (decltype(sat[0].size()) c = 0;
-                     c < nc; ++c, sp += np)
-                {
-                    *sp = sat_p[c];
-                }
-            }
-
-            return s;
-        }
-    } // namespace Details
-
-
-    /**
-     * Compute initial state by an equilibration procedure.
-     *
-     * The following state fields are modified:
-     *   pressure(),
-     *   saturation(),
-     *   surfacevol(),
-     *   gasoilratio(),
-     *   rv().
-     *
-     * \param[in] grid     Grid.
-     * \param[in] props    Property object, pvt and capillary properties are used.
-     * \param[in] deck     Simulation deck, used to obtain EQUIL and related data.
-     * \param[in] gravity  Acceleration of gravity, assumed to be in Z direction.
-     * \param[in] applySwatInit     Make it possible to not apply SWATINIT even if it
-     *                              is present in the deck
-     */
-    template<class MaterialLawManager, class Grid>
-    void initStateEquil(const Grid& grid,
-                        std::shared_ptr<MaterialLawManager> materialLawManager,
-                        const Opm::Deck& deck,
-                        const Opm::EclipseState& eclipseState,
-                        const double gravity,
-                        BlackoilState& state,
-                        bool applySwatinit = true)
-    {
-
-        typedef EQUIL::DeckDependent::InitialStateComputer ISC;
-
-        PhaseUsage pu = phaseUsageFromDeck(deck);
-
-        //Check for presence of kw SWATINIT
-        std::vector<double> swat_init = {};
-        if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatinit) {
-            const std::vector<double>& swat_init_ecl = eclipseState.
-                    get3DProperties().getDoubleGridProperty("SWATINIT").getData();
-            const int nc = UgGridHelpers::numCells(grid);
-            swat_init.resize(nc);
-            const int* gc = UgGridHelpers::globalCell(grid);
-            for (int c = 0; c < nc; ++c) {
-                const int deck_pos = (gc == NULL) ? c : gc[c];
-                swat_init[c] = swat_init_ecl[deck_pos];
-            }
-        }
-
-        ISC isc(materialLawManager, pu, deck, eclipseState, grid, gravity, swat_init);
-        const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
-            ? pu.phase_pos[BlackoilPhases::Liquid]
-            : pu.phase_pos[BlackoilPhases::Aqua];
-        state.pressure() = isc.press()[ref_phase];
-        state.saturation() = Details::convertSats(isc.saturation());
-        state.gasoilratio() = isc.rs();
-        state.rv() = isc.rv();
-
-        //initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
-    }
-
-
-
 } // namespace Opm
 
 #endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED

tests/test_equil.cpp

@@ -23,13 +23,11 @@
 #include <opm/core/grid/cart_grid.h>
 #include <opm/core/grid/GridManager.hpp>
 #include <opm/core/utility/compressedToCartesian.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
 
-
-#include <opm/core/props/BlackoilPropertiesBasic.hpp>
-#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
-#include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
@@ -42,7 +40,6 @@
 
 #include <opm/core/pressure/msmfem/partition.h>
 
-#include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/parser/eclipse/Units/Units.hpp>
 
 #include <array>
@@ -55,11 +52,12 @@
 #include <vector>
 
 
+typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
 // Forward declaring the MaterialLawManager template.
 typedef Opm::ThreePhaseMaterialTraits<double,
-/*wettingPhaseIdx=*/Opm::BlackoilPhases::Aqua,
+/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
-/*nonWettingPhaseIdx=*/Opm::BlackoilPhases::Liquid,
+/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
-/*gasPhaseIdx=*/Opm::BlackoilPhases::Vapour> MaterialTraits;
+/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits;
 typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
 
 
@@ -71,8 +69,8 @@ typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
     BOOST_CHECK_CLOSE((value), (expected), (reltol)); \
 }
 
-typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
+
-Opm::PhaseUsage initDefaultFluidSystem() {
+void initDefaultFluidSystem() {
     std::vector<std::pair<double, double> > Bo = {
         { 101353, 1. },
         { 6.21542e+07, 1 }
@@ -139,17 +137,6 @@ Opm::PhaseUsage initDefaultFluidSystem() {
 
     FluidSystem::initEnd();
 
-    Opm::PhaseUsage pu;
-    pu.num_phases = 2;
-    // Might just as well assume water-oil.
-    pu.phase_used[Opm::BlackoilPhases::Aqua] = true;
-    pu.phase_used[Opm::BlackoilPhases::Liquid] = true;
-    pu.phase_used[Opm::BlackoilPhases::Vapour] = false;
-    pu.phase_pos[Opm::BlackoilPhases::Aqua] = 0;
-    pu.phase_pos[Opm::BlackoilPhases::Liquid] = 1;
-    pu.phase_pos[Opm::BlackoilPhases::Vapour] = 1; // Unused.
-    return pu;
-
 }
 
 BOOST_AUTO_TEST_SUITE ()
@@ -220,14 +207,13 @@ BOOST_AUTO_TEST_CASE (PhasePressure)
 
     auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 );
 
-    Opm::PhaseUsage pu = initDefaultFluidSystem();
+    initDefaultFluidSystem();
 
     Opm::EQUIL::EquilReg
             region(record,
                    std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
                    std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-                   0,
+                   0);
-                   pu);
 
     std::vector<int> cells(G->number_of_cells);
     std::iota(cells.begin(), cells.end(), 0);
@@ -251,8 +237,7 @@ BOOST_AUTO_TEST_CASE (CellSubset)
     std::shared_ptr<UnstructuredGrid>
         G(create_grid_cart3d(10, 1, 10), destroy_grid);
 
-    Opm::PhaseUsage pu = initDefaultFluidSystem();
+    initDefaultFluidSystem();
-
 
     Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
                                   mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
@@ -262,26 +247,22 @@ BOOST_AUTO_TEST_CASE (CellSubset)
         Opm::EQUIL::EquilReg(record[0],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-        0,
+        0)
-        pu)
         ,
         Opm::EQUIL::EquilReg(record[0],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-        0,
+        0)
-        pu)
         ,
         Opm::EQUIL::EquilReg(record[1],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-        0,
+        0)
-        pu)
         ,
         Opm::EQUIL::EquilReg(record[1],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-        0,
+        0)
-        pu)
     };
 
     const int cdim[] = { 2, 1, 2 };
@@ -348,33 +329,29 @@ BOOST_AUTO_TEST_CASE (RegMapping)
     Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
                                   mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
 
-    Opm::PhaseUsage pu = initDefaultFluidSystem();
+    initDefaultFluidSystem();
 
     Opm::EQUIL::EquilReg region[] =
     {
         Opm::EQUIL::EquilReg(record[0],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-        0,
+        0)
-        pu)
         ,
         Opm::EQUIL::EquilReg(record[0],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-        0,
+        0)
-        pu)
         ,
         Opm::EQUIL::EquilReg(record[1],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-        0,
+        0)
-        pu)
         ,
         Opm::EQUIL::EquilReg(record[1],
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
         std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
-        0,
+        0)
-        pu)
         };
 
     std::vector<int> eqlnum(G->number_of_cells);
@@ -436,9 +413,12 @@ BOOST_AUTO_TEST_CASE (DeckAllDead)
     auto materialLawManager = std::make_shared<MaterialLawManager>();
     materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
 
-    Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
+    typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, *grid, 10.0);
+    // Initialize the fluid system
+    FluidSystem::initFromDeck(deck, eclipseState);
+
+    Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, *grid, 10.0);
     const auto& pressures = comp.press();
     BOOST_REQUIRE(pressures.size() == 3);
     BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
@@ -538,9 +518,13 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
     auto materialLawManager = std::make_shared<MaterialLawManager>();
     materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
 
-    Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
+    typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
+
+    // Initialize the fluid system
+    FluidSystem::initFromDeck(deck, eclipseState);
+
+    Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 10.0);
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 10.0);
     const auto& pressures = comp.press();
     BOOST_REQUIRE(pressures.size() == 3);
     BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -585,9 +569,13 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
     auto materialLawManager = std::make_shared<MaterialLawManager>();
     materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
 
-    Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
+    typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
+    // Initialize the fluid system
+    FluidSystem::initFromDeck(deck, eclipseState);
+
+
+    Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE(pressures.size() == 3);
     BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -654,9 +642,11 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
     auto materialLawManager = std::make_shared<MaterialLawManager>();
     materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
 
-    Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
+    typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
+    // Initialize the fluid system
+    FluidSystem::initFromDeck(deck, eclipseState);
+    Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE(pressures.size() == 3);
     BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -740,9 +730,12 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
     auto materialLawManager = std::make_shared<MaterialLawManager>();
     materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
 
-    Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
+    typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
+    // Initialize the fluid system
+    FluidSystem::initFromDeck(deck, eclipseState);
+
+    Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE(pressures.size() == 3);
     BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -828,9 +821,12 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
     auto materialLawManager = std::make_shared<MaterialLawManager>();
     materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
 
-    Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
+    typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
+    // Initialize the fluid system
+    FluidSystem::initFromDeck(deck, eclipseState);
+
+    Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE(pressures.size() == 3);
     BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -940,8 +936,6 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
     auto materialLawManagerScaled = std::make_shared<MaterialLawManager>();
     materialLawManagerScaled->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
 
-    Opm::BlackoilState state( Opm::UgGridHelpers::numCells( grid ) , Opm::UgGridHelpers::numFaces( grid ) , 3);
-
     // reference saturations
     const std::vector<double> s[3]{
         { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.42528761746004229, 0.77462669821009045, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
@@ -955,18 +949,6 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
         { 0, 0, 0, 0.014813991154779993, 0.78525420807446045, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0 },
         { 0.8, 0.8, 0.8, 0.78518600884522005, 0.014745791925539575, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
     };
-    std::vector<double> sats = state.saturation();
-    for (int phase = 0; phase < 3; ++phase) {
-        for (size_t i = 0; i < 20; ++i) {
-            sats[3*i + phase] = s[phase][i];
-        }
-    }
-    std::vector<double> sats_swatinit = state.saturation();
-    for (int phase = 0; phase < 3; ++phase) {
-        for (size_t i = 0; i < 20; ++i) {
-            sats_swatinit[3*i + phase] = swatinit[phase][i];
-        }
-    }
 
     // Adjust oil pressure according to gas saturation and cap pressure
     typedef Opm::SimpleModularFluidState<double,
@@ -989,13 +971,13 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
     FluidSystem::initFromDeck(deck, eclipseState);
 
     // reference pcs
-    std::vector<double> pc_original = state.saturation();
     const int numCells = Opm::UgGridHelpers::numCells(grid);
+    std::vector<double> pc_original(numCells * FluidSystem::numPhases);
     for (int c = 0; c < numCells; ++c) {
         std::vector<double> pc = {0,0,0};
-        double sw = sats[3*c + 0];
+        double sw = s[0][c];
-        double so = sats[3*c + 1];
+        double so = s[1][c];
-        double sg = sats[3*c + 2];
+        double sg = s[2][c];
         fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
         fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
         fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg);
@@ -1024,22 +1006,18 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
     pc_scaled_truth[3*11 + 0] =  5364.1;
 
     // compute the initial state
-
     // apply swatinit
-    Opm::BlackoilState state_scaled = state;
+    Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> compScaled(materialLawManagerScaled, eclipseState, grid, 9.81, true);
-    initStateEquil(grid, materialLawManagerScaled, deck, eclipseState, 9.81, state_scaled, true);
-
     // don't apply swatinit
-    Opm::BlackoilState state_unscaled = state;
+    Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> compUnscaled(materialLawManager, eclipseState, grid, 9.81, false);
-    initStateEquil(grid, materialLawManager, deck, eclipseState, 9.81, state_unscaled, false);
 
     // compute pc
-    std::vector<double> pc_scaled= state.saturation();
+    std::vector<double> pc_scaled(numCells * FluidSystem::numPhases);
     for (int c = 0; c < numCells; ++c) {
         std::vector<double> pc = {0,0,0};
-        double sw = state_scaled.saturation().data()[3*c + 0];
+        double sw = compScaled.saturation().data()[0][c];
-        double so = state_scaled.saturation().data()[3*c + 1];
+        double so = compScaled.saturation().data()[1][c];
-        double sg = state_scaled.saturation().data()[3*c + 2];
+        double sg = compScaled.saturation().data()[2][c];
 
         fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
         fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
@@ -1050,12 +1028,12 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
         pc_scaled[3*c + 1] = 0.0;
         pc_scaled[3*c + 2] = pc[FluidSystem::oilPhaseIdx] + pc[FluidSystem::gasPhaseIdx];
     }
-    std::vector<double> pc_unscaled= state.saturation();
+    std::vector<double> pc_unscaled(numCells * FluidSystem::numPhases);
     for (int c = 0; c < numCells; ++c) {
         std::vector<double> pc = {0,0,0};
-        double sw = state_unscaled.saturation().data()[3*c + 0];
+        double sw = compUnscaled.saturation().data()[0][c];
-        double so = state_unscaled.saturation().data()[3*c + 1];
+        double so = compUnscaled.saturation().data()[1][c];
-        double sg = state_unscaled.saturation().data()[3*c + 2];
+        double sg = compUnscaled.saturation().data()[2][c];
 
         fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
         fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
@@ -1079,8 +1057,8 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
 
     for (int phase = 0; phase < 3; ++phase) {
         for (size_t i = 0; i < 20; ++i) {
-            CHECK(state_unscaled.saturation()[3*i + phase], s[phase][i], reltol);
+            CHECK(compUnscaled.saturation()[phase][i], s[phase][i], reltol);
-            CHECK(state_scaled.saturation()[3*i + phase], swatinit[phase][i], reltol);
+            CHECK(compScaled.saturation()[phase][i], swatinit[phase][i], reltol);
         }
     }
 }