23#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
24#define OPM_INITSTATEEQUIL_HEADER_INCLUDED
26#include <opm/core/grid/GridHelpers.hpp>
28#include <opm/core/utility/RegionMapping.hpp>
29#include <opm/core/utility/extractPvtTableIndex.hpp>
31#include <opm/parser/eclipse/Units/Units.hpp>
32#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
33#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
34#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp>
35#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
36#include <opm/parser/eclipse/EclipseState/Tables/TableContainer.hpp>
37#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
38#include <opm/parser/eclipse/EclipseState/Tables/RsvdTable.hpp>
39#include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
40#include <opm/common/OpmLog/OpmLog.hpp>
41#include <opm/common/data/SimulationDataContainer.hpp>
43#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
44#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
45#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
57struct UnstructuredGrid;
110 template <
class Gr
id,
class Region,
class CellRange,
class Flu
idSystem>
111 std::vector< std::vector<double> >
114 const CellRange& cells,
115 const double grav = unit::gravity);
153 template <
class Flu
idSystem,
class Gr
id,
class Region,
class CellRange,
class MaterialLawManager>
154 std::vector< std::vector<double> >
157 const CellRange& cells,
158 MaterialLawManager& materialLawManager,
159 const std::vector<double> swat_init,
160 std::vector< std::vector<double> >& phase_pressures);
182 template <
class Gr
id,
class CellRangeType>
183 std::vector<double>
computeRs(
const Grid& grid,
184 const CellRangeType& cells,
185 const std::vector<double> oil_pressure,
187 const Miscibility::RsFunction& rs_func,
188 const std::vector<double> gas_saturation);
190 namespace DeckDependent {
192 std::vector<EquilRecord>
193 getEquil(
const Opm::EclipseState& state)
195 const auto& init = state.getInitConfig();
197 if( !init.hasEquil() ) {
198 OPM_THROW(std::domain_error,
"Deck does not provide equilibration data.");
201 const auto& equil = init.getEquil();
202 return { equil.begin(), equil.end() };
208 equilnum(
const Opm::EclipseState& eclipseState,
211 std::vector<int> eqlnum;
212 if (eclipseState.get3DProperties().hasDeckIntGridProperty(
"EQLNUM")) {
213 const int nc = UgGridHelpers::numCells(G);
215 const std::vector<int>& e =
216 eclipseState.get3DProperties().getIntGridProperty(
"EQLNUM").getData();
217 const int* gc = UgGridHelpers::globalCell(G);
218 for (
int cell = 0; cell < nc; ++cell) {
219 const int deck_pos = (gc == NULL) ? cell : gc[cell];
220 eqlnum[cell] = e[deck_pos] - 1;
226 eqlnum.assign(UgGridHelpers::numCells(G), 0);
232 template<
class Flu
idSystem>
233 class InitialStateComputer {
235 template<
class MaterialLawManager,
class Gr
id>
236 InitialStateComputer(MaterialLawManager& materialLawManager,
237 const Opm::EclipseState& eclipseState,
239 const double grav = unit::gravity,
240 const bool applySwatInit =
true
242 : pp_(FluidSystem::numPhases,
243 std::vector<double>(UgGridHelpers::numCells(G))),
244 sat_(FluidSystem::numPhases,
245 std::vector<double>(UgGridHelpers::numCells(G))),
246 rs_(UgGridHelpers::numCells(G)),
247 rv_(UgGridHelpers::numCells(G))
250 if (eclipseState.get3DProperties().hasDeckDoubleGridProperty(
"SWATINIT") && applySwatInit) {
251 const std::vector<double>& swat_init_ecl = eclipseState.
252 get3DProperties().getDoubleGridProperty(
"SWATINIT").getData();
253 const int nc = UgGridHelpers::numCells(G);
254 swat_init_.resize(nc);
255 const int* gc = UgGridHelpers::globalCell(G);
256 for (
int c = 0; c < nc; ++c) {
257 const int deck_pos = (gc == NULL) ? c : gc[c];
258 swat_init_[c] = swat_init_ecl[deck_pos];
262 const std::vector<EquilRecord> rec = getEquil(eclipseState);
263 const auto& tables = eclipseState.getTableManager();
265 const RegionMapping<> eqlmap(equilnum(eclipseState, G));
266 const int invalidRegion = -1;
267 regionPvtIdx_.resize(rec.size(), invalidRegion);
268 setRegionPvtIdx(G, eclipseState, eqlmap);
271 rs_func_.reserve(rec.size());
272 if (FluidSystem::enableDissolvedGas()) {
273 const TableContainer& rsvdTables = tables.getRsvdTables();
274 for (
size_t i = 0; i < rec.size(); ++i) {
275 if (eqlmap.cells(i).empty())
277 rs_func_.push_back(std::shared_ptr<Miscibility::RsVD<FluidSystem>>());
280 const int pvtIdx = regionPvtIdx_[i];
281 if (!rec[i].liveOilInitConstantRs()) {
282 if (rsvdTables.size() <= 0 ) {
283 OPM_THROW(std::runtime_error,
"Cannot initialise: RSVD table not available.");
285 const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
286 std::vector<double> depthColumn = rsvdTable.getColumn(
"DEPTH").vectorCopy();
287 std::vector<double> rsColumn = rsvdTable.getColumn(
"RS").vectorCopy();
288 rs_func_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
289 depthColumn , rsColumn));
291 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
292 OPM_THROW(std::runtime_error,
293 "Cannot initialise: when no explicit RSVD table is given, \n"
294 "datum depth must be at the gas-oil-contact. "
295 "In EQUIL region " << (i + 1) <<
" (counting from 1), this does not hold.");
297 const double p_contact = rec[i].datumDepthPressure();
298 const double T_contact = 273.15 + 20;
299 rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx, p_contact, T_contact));
303 for (
size_t i = 0; i < rec.size(); ++i) {
304 rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
308 rv_func_.reserve(rec.size());
309 if (FluidSystem::enableVaporizedOil()) {
310 const TableContainer& rvvdTables = tables.getRvvdTables();
311 for (
size_t i = 0; i < rec.size(); ++i) {
312 if (eqlmap.cells(i).empty())
314 rv_func_.push_back(std::shared_ptr<Miscibility::RvVD<FluidSystem>>());
317 const int pvtIdx = regionPvtIdx_[i];
318 if (!rec[i].wetGasInitConstantRv()) {
319 if (rvvdTables.size() <= 0) {
320 OPM_THROW(std::runtime_error,
"Cannot initialise: RVVD table not available.");
323 const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
324 std::vector<double> depthColumn = rvvdTable.getColumn(
"DEPTH").vectorCopy();
325 std::vector<double> rvColumn = rvvdTable.getColumn(
"RV").vectorCopy();
326 rv_func_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
327 depthColumn , rvColumn));
330 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
331 OPM_THROW(std::runtime_error,
332 "Cannot initialise: when no explicit RVVD table is given, \n"
333 "datum depth must be at the gas-oil-contact. "
334 "In EQUIL region " << (i + 1) <<
" (counting from 1), this does not hold.");
336 const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
337 const double T_contact = 273.15 + 20;
338 rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx ,p_contact, T_contact));
342 for (
size_t i = 0; i < rec.size(); ++i) {
343 rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
348 calcPressSatRsRv(eqlmap, rec, materialLawManager, G, grav);
354 typedef std::vector<double> Vec;
355 typedef std::vector<Vec> PVec;
357 const PVec& press()
const {
return pp_; }
358 const PVec& saturation()
const {
return sat_; }
359 const Vec& rs()
const {
return rs_; }
360 const Vec& rv()
const {
return rv_; }
363 typedef EquilReg EqReg;
364 std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
365 std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
366 std::vector<int> regionPvtIdx_;
373 template<
class Gr
id,
class RMap>
374 void setRegionPvtIdx(
const Grid& G,
const Opm::EclipseState& eclipseState,
const RMap& reg) {
376 std::vector<int> cellPvtRegionIdx;
377 extractPvtTableIndex(cellPvtRegionIdx, eclipseState, UgGridHelpers::numCells(G), UgGridHelpers::globalCell(G));
378 for (
const auto& r : reg.activeRegions()) {
379 const auto& cells = reg.cells(r);
380 const int cell = *(cells.begin());
381 regionPvtIdx_[r] = cellPvtRegionIdx[cell];
385 template <
class RMap,
class MaterialLawManager,
class Gr
id>
387 calcPressSatRsRv(
const RMap& reg ,
388 const std::vector< EquilRecord >& rec ,
389 MaterialLawManager& materialLawManager,
393 for (
const auto& r : reg.activeRegions()) {
394 const auto& cells = reg.cells(r);
397 OpmLog::warning(
"Equilibration region " + std::to_string(r + 1)
398 +
" has no active cells");
402 const EqReg eqreg(rec[r], rs_func_[r], rv_func_[r], regionPvtIdx_[r]);
404 PVec pressures = phasePressures<FluidSystem>(G, eqreg, cells, grav);
405 const std::vector<double>& temp =
temperature(G, eqreg, cells);
406 const PVec sat = phaseSaturations<FluidSystem>(G, eqreg, cells, materialLawManager, swat_init_, pressures);
408 const int np = FluidSystem::numPhases;
409 for (
int p = 0; p < np; ++p) {
410 copyFromRegion(pressures[p], cells, pp_[p]);
411 copyFromRegion(sat[p], cells, sat_[p]);
413 const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
414 const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
416 const int oilpos = FluidSystem::oilPhaseIdx;
417 const int gaspos = FluidSystem::gasPhaseIdx;
418 const Vec rs_vals =
computeRs(G, cells, pressures[oilpos], temp, *(rs_func_[r]), sat[gaspos]);
419 const Vec rv_vals =
computeRs(G, cells, pressures[gaspos], temp, *(rv_func_[r]), sat[oilpos]);
420 copyFromRegion(rs_vals, cells, rs_);
421 copyFromRegion(rv_vals, cells, rv_);
426 template <
class CellRangeType>
427 void copyFromRegion(
const Vec& source,
428 const CellRangeType& cells,
431 auto s = source.begin();
432 auto c = cells.begin();
433 const auto e = cells.end();
434 for (; c != e; ++c, ++s) {
435 destination[*c] = *s;
void gas(const Grid &G, const Region ®, const std::array< double, 2 > &span, const double grav, double &po_goc, const CellRange &cells, std::vector< double > &press)
Definition: initStateEquil_impl.hpp:378
void oil(const Grid &G, const Region ®, const std::array< double, 2 > &span, const double grav, const CellRange &cells, std::vector< double > &press, double &po_woc, double &po_goc)
Definition: initStateEquil_impl.hpp:319
std::vector< double > computeRs(const Grid &grid, const CellRangeType &cells, const std::vector< double > oil_pressure, const std::vector< double > &temperature, const Miscibility::RsFunction &rs_func, const std::vector< double > gas_saturation)
Definition: initStateEquil_impl.hpp:769
std::vector< std::vector< double > > phaseSaturations(const Grid &G, const Region ®, const CellRange &cells, MaterialLawManager &materialLawManager, const std::vector< double > swat_init, std::vector< std::vector< double > > &phase_pressures)
Definition: initStateEquil_impl.hpp:593
std::vector< std::vector< double > > phasePressures(const Grid &G, const Region ®, const CellRange &cells, const double grav)
Definition: initStateEquil_impl.hpp:510
std::vector< double > temperature(const Grid &, const Region &, const CellRange &cells)
Definition: initStateEquil_impl.hpp:583
Definition: AnisotropicEikonal.hpp:44