From dc997185d92553818ce0d964cbf1c9489aa72165 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 20 Nov 2017 08:47:41 +0100 Subject: [PATCH] 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 | 64 +++++- opm/core/simulator/EquilibrationHelpers.hpp | 20 +- opm/core/simulator/initStateEquil.hpp | 70 +++---- opm/core/simulator/initStateEquil_impl.hpp | 212 +++----------------- tests/test_equil.cpp | 156 +++++++------- 5 files changed, 205 insertions(+), 317 deletions(-) diff --git a/examples/compute_initial_state.cpp b/examples/compute_initial_state.cpp index 67c00a363..05f22e590 100644 --- a/examples/compute_initial_state.cpp +++ b/examples/compute_initial_state.cpp @@ -28,6 +28,8 @@ #include #include #include +#include +#include #include #include @@ -75,12 +77,45 @@ namespace std::copy(data.begin(), data.end(), std::ostream_iterator(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 + void convertSats(std::vector& sat_interleaved, const std::vector< std::vector >& 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& 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& 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& 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 FluidSystem; + FluidSystem::initFromDeck(deck, eclipseState); + PhaseUsage pu = phaseUsageFromDeck(deck); + + typedef EQUIL::DeckDependent::InitialStateComputer 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(state.saturation(), isc.saturation(), pu); + + if (state.hasCellData(std::string("GASOILRATIO"))) { + std::vector& rs = state.getCellData(std::string("GASOILRATIO")); + rs = isc.rs(); + } + if (state.hasCellData(std::string("RV"))){ + std::vector& rv = state.getCellData(std::string("RV")); + rv = isc.rv(); + } + + // Output. const std::string output_dir = param.getDefault("output_dir", "output"); @@ -132,3 +192,5 @@ catch (const std::exception& e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; throw; } + + diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index 1c2643265..eff050102 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/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 #include #include #include @@ -478,13 +478,11 @@ namespace Opm EquilReg(const EquilRecord& rec, std::shared_ptr rs, std::shared_ptr rv, - const int pvtIdx, - const PhaseUsage& pu) + const int pvtIdx) : 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 rs_; /**< RS calculator */ std::shared_ptr 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 @@ -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]; }; diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 4d8956078..d3f1110d4 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -24,9 +24,6 @@ #include #include -#include -#include -#include #include #include @@ -40,6 +37,8 @@ #include #include #include +#include + #include #include @@ -198,12 +197,12 @@ namespace Opm template inline std::vector - equilnum(const Opm::Deck& deck, - const Opm::EclipseState& eclipseState, + equilnum(const Opm::EclipseState& eclipseState, const Grid& G ) { std::vector eqlnum; - if (deck.hasKeyword("EQLNUM")) { + + if (eclipseState.get3DProperties().hasDeckIntGridProperty("EQLNUM")) { const int nc = UgGridHelpers::numCells(G); eqlnum.resize(nc); const std::vector& e = @@ -223,45 +222,48 @@ namespace Opm return eqlnum; } - + template class InitialStateComputer { public: template InitialStateComputer(std::shared_ptr materialLawManager, - const PhaseUsage& phaseUsage, - const Opm::Deck& deck, const Opm::EclipseState& eclipseState, const Grid& G , const double grav = unit::gravity, - const std::vector& swat_init = {} + const bool applySwatInit = true ) - : pp_(phaseUsage.num_phases, + : pp_(FluidSystem::numPhases, std::vector(UgGridHelpers::numCells(G))), - sat_(phaseUsage.num_phases, + sat_(FluidSystem::numPhases, std::vector(UgGridHelpers::numCells(G))), rs_(UgGridHelpers::numCells(G)), - rv_(UgGridHelpers::numCells(G)), - swat_init_(swat_init), - phaseUsage_(phaseUsage) - + rv_(UgGridHelpers::numCells(G)) { - typedef FluidSystems::BlackOil FluidSystem; - - // Initialize the fluid system - FluidSystem::initFromDeck(deck, eclipseState); + //Check for presence of kw SWATINIT + if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatInit) { + const std::vector& swat_init_ecl = 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 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(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 void setRegionPvtIdx(const Grid& G, const Opm::EclipseState& eclipseState, const RMap& reg) { @@ -381,7 +382,7 @@ namespace Opm } } - template + template void calcPressSatRsRv(const RMap& reg , const std::vector< EquilRecord >& rec , @@ -398,23 +399,24 @@ namespace Opm continue; } - const EqReg eqreg(rec[r], - rs_func_[r], rv_func_[r], regionPvtIdx_[r], - phaseUsage_); + const EqReg eqreg(rec[r], rs_func_[r], rv_func_[r], regionPvtIdx_[r]); PVec pressures = phasePressures(G, eqreg, cells, grav); const std::vector& temp = temperature(G, eqreg, cells); const PVec sat = phaseSaturations(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 int oilpos = phaseUsage_.phase_pos[BlackoilPhases::Liquid]; - const int gaspos = phaseUsage_.phase_pos[BlackoilPhases::Vapour]; + + const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + + 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_); diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index 42314faf0..a2435c47e 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -25,7 +25,7 @@ #include #include -#include +#include #include #include @@ -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 >& 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)) { - const int wix = PhaseIndex::water(pu); + if (water) { PhasePressure::water(G, reg, span, grav, po_woc, - cells, press[ wix ]); + cells, press[ waterpos ]); } - if (PhaseUsed::oil(pu)) { - const int oix = PhaseIndex::oil(pu); + if (oil) { PhasePressure::oil(G, reg, span, grav, cells, - press[ oix ], po_woc, po_goc); + press[ oilpos ], po_woc, po_goc); } - if (PhaseUsed::gas(pu)) { - const int gix = PhaseIndex::gas(pu); + if (gas) { PhasePressure::gas(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)) { - const int gix = PhaseIndex::gas(pu); + if (gas) { PhasePressure::gas(G, reg, span, grav, po_goc, - cells, press[ gix ]); + cells, press[ gaspos ]); } - if (PhaseUsed::oil(pu)) { - const int oix = PhaseIndex::oil(pu); + if (oil) { PhasePressure::oil(G, reg, span, grav, cells, - press[ oix ], po_woc, po_goc); + press[ oilpos ], po_woc, po_goc); } - if (PhaseUsed::water(pu)) { - const int wix = PhaseIndex::water(pu); + if (water) { PhasePressure::water(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)) { - const int oix = PhaseIndex::oil(pu); + if (oil) { PhasePressure::oil(G, reg, span, grav, cells, - press[ oix ], po_woc, po_goc); + press[ oilpos ], po_woc, po_goc); } - if (PhaseUsed::water(pu)) { - const int wix = PhaseIndex::water(pu); + if (water) { PhasePressure::water(G, reg, span, grav, po_woc, - cells, press[ wix ]); + cells, press[ waterpos ]); } - if (PhaseUsed::gas(pu)) { - const int gix = PhaseIndex::gas(pu); + if (gas) { PhasePressure::gas(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 pval; std::vector press(np, pval(ncell, 0.0)); @@ -659,7 +601,7 @@ namespace Opm const std::vector swat_init, std::vector< std::vector >& 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 gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour]; - const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid]; - const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua]; - const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour]; + const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); + const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + const int oilpos = FluidSystem::oilPhaseIdx; + const int waterpos = FluidSystem::waterPhaseIdx; + const int gaspos = FluidSystem::gasPhaseIdx; std::vector::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 - convertSats(const std::vector< std::vector >& sat) - { - const auto np = sat.size(); - const auto nc = sat[0].size(); - - std::vector 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 - void initStateEquil(const Grid& grid, - std::shared_ptr 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 swat_init = {}; - if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatinit) { - const std::vector& 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 diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index c5a6069a3..0ce421404 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -23,13 +23,11 @@ #include #include #include +#include - -#include -#include -#include #include #include +#include #include @@ -42,7 +40,6 @@ #include -#include #include #include @@ -55,11 +52,12 @@ #include +typedef Opm::FluidSystems::BlackOil FluidSystem; // Forward declaring the MaterialLawManager template. typedef Opm::ThreePhaseMaterialTraits MaterialTraits; +/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx, +/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx, +/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits; typedef Opm::EclMaterialLawManager MaterialLawManager; @@ -71,8 +69,8 @@ typedef Opm::EclMaterialLawManager MaterialLawManager; BOOST_CHECK_CLOSE((value), (expected), (reltol)); \ } -typedef Opm::FluidSystems::BlackOil FluidSystem; -Opm::PhaseUsage initDefaultFluidSystem() { + +void initDefaultFluidSystem() { std::vector > 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(), std::make_shared(), - 0, - pu); + 0); std::vector cells(G->number_of_cells); std::iota(cells.begin(), cells.end(), 0); @@ -251,8 +237,7 @@ BOOST_AUTO_TEST_CASE (CellSubset) std::shared_ptr 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(), std::make_shared(), - 0, - pu) + 0) , Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), - 0, - pu) + 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), - 0, - pu) + 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), - 0, - pu) + 0) }; 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(), std::make_shared(), - 0, - pu) + 0) , Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), - 0, - pu) + 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), - 0, - pu) + 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), - 0, - pu) + 0) }; std::vector eqlnum(G->number_of_cells); @@ -436,9 +413,12 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) auto materialLawManager = std::make_shared(); materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx); - Opm::PhaseUsage pu = phaseUsageFromDeck(deck); + typedef Opm::FluidSystems::BlackOil 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 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->initFromDeck(deck, eclipseState, compressedToCartesianIdx); - Opm::PhaseUsage pu = phaseUsageFromDeck(deck); + typedef Opm::FluidSystems::BlackOil FluidSystem; + + // Initialize the fluid system + FluidSystem::initFromDeck(deck, eclipseState); + + Opm::EQUIL::DeckDependent::InitialStateComputer 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->initFromDeck(deck, eclipseState, compressedToCartesianIdx); - Opm::PhaseUsage pu = phaseUsageFromDeck(deck); + typedef Opm::FluidSystems::BlackOil 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 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->initFromDeck(deck, eclipseState, compressedToCartesianIdx); - Opm::PhaseUsage pu = phaseUsageFromDeck(deck); + typedef Opm::FluidSystems::BlackOil 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 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->initFromDeck(deck, eclipseState, compressedToCartesianIdx); - Opm::PhaseUsage pu = phaseUsageFromDeck(deck); + typedef Opm::FluidSystems::BlackOil 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 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->initFromDeck(deck, eclipseState, compressedToCartesianIdx); - Opm::PhaseUsage pu = phaseUsageFromDeck(deck); + typedef Opm::FluidSystems::BlackOil 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 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(); materialLawManagerScaled->initFromDeck(deck, eclipseState, compressedToCartesianIdx); - Opm::BlackoilState state( Opm::UgGridHelpers::numCells( grid ) , Opm::UgGridHelpers::numFaces( grid ) , 3); - // reference saturations const std::vector 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 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 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 pc_original = state.saturation(); const int numCells = Opm::UgGridHelpers::numCells(grid); + std::vector pc_original(numCells * FluidSystem::numPhases); for (int c = 0; c < numCells; ++c) { std::vector pc = {0,0,0}; - double sw = sats[3*c + 0]; - double so = sats[3*c + 1]; - double sg = sats[3*c + 2]; + double sw = s[0][c]; + double so = s[1][c]; + 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; - initStateEquil(grid, materialLawManagerScaled, deck, eclipseState, 9.81, state_scaled, true); - + Opm::EQUIL::DeckDependent::InitialStateComputer compScaled(materialLawManagerScaled, eclipseState, grid, 9.81, true); // don't apply swatinit - Opm::BlackoilState state_unscaled = state; - initStateEquil(grid, materialLawManager, deck, eclipseState, 9.81, state_unscaled, false); + Opm::EQUIL::DeckDependent::InitialStateComputer compUnscaled(materialLawManager, eclipseState, grid, 9.81, false); // compute pc - std::vector pc_scaled= state.saturation(); + std::vector pc_scaled(numCells * FluidSystem::numPhases); for (int c = 0; c < numCells; ++c) { std::vector pc = {0,0,0}; - double sw = state_scaled.saturation().data()[3*c + 0]; - double so = state_scaled.saturation().data()[3*c + 1]; - double sg = state_scaled.saturation().data()[3*c + 2]; + double sw = compScaled.saturation().data()[0][c]; + double so = compScaled.saturation().data()[1][c]; + 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 pc_unscaled= state.saturation(); + std::vector pc_unscaled(numCells * FluidSystem::numPhases); for (int c = 0; c < numCells; ++c) { std::vector pc = {0,0,0}; - double sw = state_unscaled.saturation().data()[3*c + 0]; - double so = state_unscaled.saturation().data()[3*c + 1]; - double sg = state_unscaled.saturation().data()[3*c + 2]; + double sw = compUnscaled.saturation().data()[0][c]; + double so = compUnscaled.saturation().data()[1][c]; + 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(state_scaled.saturation()[3*i + phase], swatinit[phase][i], reltol); + CHECK(compUnscaled.saturation()[phase][i], s[phase][i], reltol); + CHECK(compScaled.saturation()[phase][i], swatinit[phase][i], reltol); } } }