diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 98ebadbb3..fd8c35a13 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -22,7 +22,6 @@ #include #include -#include #include #include #include @@ -60,13 +59,6 @@ namespace Opm * \param[in] deck Simulation deck, used to obtain EQUIL and related data. * \param[in] gravity Acceleration of gravity, assumed to be in Z direction. */ - void initStateEquil(const UnstructuredGrid& grid, - const BlackoilPropertiesInterface& props, - const EclipseGridParser& deck, - const double gravity, - BlackoilState& state); - - void initStateEquil(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, const Opm::DeckConstPtr newParserDeck, @@ -188,54 +180,6 @@ namespace Opm const std::vector gas_saturation); namespace DeckDependent { - - inline - std::vector - getEquil(const EclipseGridParser& deck) - { - if (deck.hasField("EQUIL")) { - const EQUIL& eql = deck.getEQUIL(); - - typedef std::vector::size_type sz_t; - const sz_t nrec = eql.equil.size(); - - std::vector ret; - ret.reserve(nrec); - for (sz_t r = 0; r < nrec; ++r) { - const EquilLine& rec = eql.equil[r]; - - EquilRecord record = - { - { rec.datum_depth_ , - rec.datum_depth_pressure_ } - , - { rec.water_oil_contact_depth_ , - rec.oil_water_cap_pressure_ } - , - { rec.gas_oil_contact_depth_ , - rec.gas_oil_cap_pressure_ } - , - rec.live_oil_table_index_ - , - rec.wet_gas_table_index_ - , - rec.N_ - }; - if (record.N != 0) { - OPM_THROW(std::domain_error, - "kw EQUIL, item 9: Only N=0 supported."); - } - ret.push_back(record); - } - - return ret; - } - else { - OPM_THROW(std::domain_error, - "Deck does not provide equilibration data."); - } - } - inline std::vector getEquil(const Opm::DeckConstPtr newParserDeck) @@ -282,29 +226,6 @@ namespace Opm } } - - inline - std::vector - equilnum(const EclipseGridParser& deck, - const UnstructuredGrid& G ) - { - std::vector eqlnum; - if (deck.hasField("EQLNUM")) { - const std::vector& e = deck.getIntegerValue("EQLNUM"); - eqlnum.reserve(e.size()); - std::transform(e.begin(), e.end(), std::back_inserter(eqlnum), - std::bind2nd(std::minus(), 1)); - } - else { - // No explicit equilibration region. - // All cells in region zero. - eqlnum.assign(G.number_of_cells, 0); - } - - return eqlnum; - } - - inline std::vector equilnum(const Opm::DeckConstPtr newParserDeck, @@ -331,180 +252,7 @@ namespace Opm } - template - class InitialStateComputer; - - template <> - class InitialStateComputer { - public: - InitialStateComputer(const BlackoilPropertiesInterface& props, - const EclipseGridParser& deck , - const UnstructuredGrid& G , - const double grav = unit::gravity) - : pp_(props.numPhases(), - std::vector(G.number_of_cells)), - sat_(props.numPhases(), - std::vector(G.number_of_cells)), - rs_(G.number_of_cells), - rv_(G.number_of_cells) - { - // Get the equilibration records. - const std::vector rec = getEquil(deck); - - // Create (inverse) region mapping. - const RegionMapping<> eqlmap(equilnum(deck, G)); - - // Create Rs functions. - rs_func_.reserve(rec.size()); - if (deck.hasField("DISGAS")) { - for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i).begin()); - if (rec[i].live_oil_table_index > 0) { - if (deck.hasField("RSVD")) { - // TODO When this kw is actually parsed, also check for proper number of available tables - // For now, just use dummy ... - std::vector depth; depth.push_back(0.0); depth.push_back(100.0); - std::vector rs; rs.push_back(0.0); rs.push_back(100.0); - rs_func_.push_back(std::make_shared(props, cell, depth, rs)); - } else { - OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available."); - } - } else { - if (rec[i].goc.depth != rec[i].main.depth) { - OPM_THROW(std::runtime_error, - "Cannot initialise: when no explicit RSVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); - } - const double p_contact = rec[i].main.press; - rs_func_.push_back(std::make_shared(props, cell, p_contact)); - } - } - } else { - for (size_t i = 0; i < rec.size(); ++i) { - rs_func_.push_back(std::make_shared()); - } - } - - rv_func_.reserve(rec.size()); - if (deck.hasField("VAPOIL")) { - for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i).begin()); - if (rec[i].wet_gas_table_index > 0) { - if (deck.hasField("RVVD")) { - // TODO When this kw is actually parsed, also check for proper number of available tables - // For now, just use dummy ... - std::vector depth; depth.push_back(0.0); depth.push_back(100.0); - std::vector rv; rv.push_back(0.0); rv.push_back(0.0001); - rv_func_.push_back(std::make_shared(props, cell, depth, rv)); - } else { - OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available."); - } - } else { - if (rec[i].goc.depth != rec[i].main.depth) { - OPM_THROW(std::runtime_error, - "Cannot initialise: when no explicit RVVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); - } - const double p_contact = rec[i].main.press + rec[i].goc.press; - rv_func_.push_back(std::make_shared(props, cell, p_contact)); - } - } - } else { - for (size_t i = 0; i < rec.size(); ++i) { - rv_func_.push_back(std::make_shared()); - } - } - - // Compute pressures, saturations, rs and rv factors. - calcPressSatRsRv(eqlmap, rec, props, 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. - } - - typedef std::vector Vec; - typedef std::vector PVec; // One per phase. - - const PVec& press() const { return pp_; } - const PVec& saturation() const { return sat_; } - const Vec& rs() const { return rs_; } - const Vec& rv() const { return rv_; } - - private: - typedef DensityCalculator RhoCalc; - typedef EquilReg EqReg; - - std::vector< std::shared_ptr > rs_func_; - std::vector< std::shared_ptr > rv_func_; - - PVec pp_; - PVec sat_; - Vec rs_; - Vec rv_; - - template - void - calcPressSatRsRv(const RMap& reg , - const std::vector< EquilRecord >& rec , - const Opm::BlackoilPropertiesInterface& props, - const UnstructuredGrid& G , - const double grav) - { - typedef Miscibility::NoMixing NoMix; - - for (typename RMap::RegionId - r = 0, nr = reg.numRegions(); - r < nr; ++r) - { - const typename RMap::CellRange cells = reg.cells(r); - - const int repcell = *cells.begin(); - const RhoCalc calc(props, repcell); - const EqReg eqreg(rec[r], calc, - rs_func_[r], rv_func_[r], - props.phaseUsage()); - - PVec press = phasePressures(G, eqreg, cells, grav); - - const PVec sat = phaseSaturations(eqreg, cells, props, press); - - const int np = props.numPhases(); - for (int p = 0; p < np; ++p) { - copyFromRegion(press[p], cells, pp_[p]); - copyFromRegion(sat[p], cells, sat_[p]); - } - if (props.phaseUsage().phase_used[BlackoilPhases::Liquid] - && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) { - const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid]; - const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour]; - const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]); - const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]); - copyFromRegion(rs, cells, rs_); - copyFromRegion(rv, cells, rv_); - } - } - } - - template - void copyFromRegion(const Vec& source, - const CellRangeType& cells, - Vec& destination) - { - auto s = source.begin(); - auto c = cells.begin(); - const auto e = cells.end(); - for (; c != e; ++c, ++s) { - destination[*c] = *s; - } - } - - }; - - - template <> - class InitialStateComputer { + class InitialStateComputer { public: InitialStateComputer(const BlackoilPropertiesInterface& props, const Opm::DeckConstPtr newParserDeck, diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index d227d54d1..cd08a98ec 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -735,33 +735,13 @@ namespace Opm * \param[in] deck Simulation deck, used to obtain EQUIL and related data. * \param[in] gravity Acceleration of gravity, assumed to be in Z direction. */ - void initStateEquil(const UnstructuredGrid& grid, - const BlackoilPropertiesInterface& props, - const EclipseGridParser& deck, - const double gravity, - BlackoilState& state) - { - typedef Equil::DeckDependent::InitialStateComputer ISC; - ISC isc(props, deck, grid, gravity); - const auto pu = props.phaseUsage(); - 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() = convertSats(isc.saturation()); - state.gasoilratio() = isc.rs(); - state.rv() = isc.rv(); - // TODO: state.surfacevol() must be computed from s, rs, rv. - } - - void initStateEquil(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, const Opm::DeckConstPtr newParserDeck, const double gravity, BlackoilState& state) { - typedef Equil::DeckDependent::InitialStateComputer ISC; + typedef Equil::DeckDependent::InitialStateComputer ISC; ISC isc(props, newParserDeck, grid, gravity); const auto pu = props.phaseUsage(); const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]