22 #ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
23 #define OPM_INITSTATEEQUIL_HEADER_INCLUDED
25 #include <opm/core/grid/GridHelpers.hpp>
26 #include <opm/core/simulator/EquilibrationHelpers.hpp>
27 #include <opm/core/simulator/BlackoilState.hpp>
28 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
29 #include <opm/core/props/BlackoilPhases.hpp>
30 #include <opm/core/utility/RegionMapping.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>
52 struct UnstructuredGrid;
75 void initStateEquil(
const Grid& grid,
76 const BlackoilPropertiesInterface& props,
77 const Opm::Deck& deck,
78 const Opm::EclipseState& eclipseState,
81 bool applySwatInit =
true);
131 template <
class Gr
id,
class Region,
class CellRange>
132 std::vector< std::vector<double> >
133 phasePressures(
const Grid& G,
135 const CellRange& cells,
136 const double grav = unit::gravity);
164 template <
class Gr
id,
class Region,
class CellRange>
165 std::vector< std::vector<double> >
166 phaseSaturations(
const Grid& grid,
168 const CellRange& cells,
169 BlackoilPropertiesFromDeck& props,
170 const std::vector<double> swat_init,
171 std::vector< std::vector<double> >& phase_pressures);
193 template <
class Gr
id,
class CellRangeType>
194 std::vector<double>
computeRs(
const Grid& grid,
195 const CellRangeType& cells,
196 const std::vector<double> oil_pressure,
197 const std::vector<double>& temperature,
198 const Miscibility::RsFunction& rs_func,
199 const std::vector<double> gas_saturation);
201 namespace DeckDependent {
203 std::vector<EquilRecord>
204 getEquil(
const Opm::EclipseState& state)
206 const auto& init = state.getInitConfig();
208 if( !init.hasEquil() ) {
209 OPM_THROW(std::domain_error,
"Deck does not provide equilibration data.");
212 const auto& equil = init.getEquil();
213 return { equil.begin(), equil.end() };
219 equilnum(
const Opm::Deck& deck,
220 const Opm::EclipseState& eclipseState,
223 std::vector<int> eqlnum;
224 if (deck.hasKeyword(
"EQLNUM")) {
225 const int nc = UgGridHelpers::numCells(G);
227 const std::vector<int>& e =
228 eclipseState.get3DProperties().getIntGridProperty(
"EQLNUM").getData();
229 const int* gc = UgGridHelpers::globalCell(G);
230 for (
int cell = 0; cell < nc; ++cell) {
231 const int deck_pos = (gc == NULL) ? cell : gc[cell];
232 eqlnum[cell] = e[deck_pos] - 1;
238 eqlnum.assign(UgGridHelpers::numCells(G), 0);
245 class InitialStateComputer {
248 InitialStateComputer(BlackoilPropertiesInterface& props,
249 const Opm::Deck& deck,
250 const Opm::EclipseState& eclipseState,
252 const double grav = unit::gravity,
253 const std::vector<double>& swat_init = {}
255 : pp_(props.numPhases(),
256 std::vector<double>(UgGridHelpers::numCells(G))),
257 sat_(props.numPhases(),
258 std::vector<double>(UgGridHelpers::numCells(G))),
259 rs_(UgGridHelpers::numCells(G)),
260 rv_(UgGridHelpers::numCells(G)),
261 swat_init_(swat_init)
265 const std::vector<EquilRecord> rec = getEquil(eclipseState);
266 const auto& tables = eclipseState.getTableManager();
268 const RegionMapping<> eqlmap(equilnum(deck, eclipseState, G));
271 rs_func_.reserve(rec.size());
272 if (deck.hasKeyword(
"DISGAS")) {
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>());
280 const int cell = *(eqlmap.cells(i).begin());
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>(props,
290 depthColumn , rsColumn));
292 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
293 OPM_THROW(std::runtime_error,
294 "Cannot initialise: when no explicit RSVD table is given, \n"
295 "datum depth must be at the gas-oil-contact. "
296 "In EQUIL region " << (i + 1) <<
" (counting from 1), this does not hold.");
298 const double p_contact = rec[i].datumDepthPressure();
299 const double T_contact = 273.15 + 20;
300 rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact, T_contact));
304 for (
size_t i = 0; i < rec.size(); ++i) {
305 rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
309 rv_func_.reserve(rec.size());
310 if (deck.hasKeyword(
"VAPOIL")) {
311 const TableContainer& rvvdTables = tables.getRvvdTables();
312 for (
size_t i = 0; i < rec.size(); ++i) {
313 if (eqlmap.cells(i).empty())
315 rv_func_.push_back(std::shared_ptr<Miscibility::RvVD>());
318 const int cell = *(eqlmap.cells(i).begin());
319 if (!rec[i].wetGasInitConstantRv()) {
320 if (rvvdTables.size() <= 0) {
321 OPM_THROW(std::runtime_error,
"Cannot initialise: RVVD table not available.");
324 const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
325 std::vector<double> depthColumn = rvvdTable.getColumn(
"DEPTH").vectorCopy();
326 std::vector<double> rvColumn = rvvdTable.getColumn(
"RV").vectorCopy();
328 rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props,
330 depthColumn , rvColumn));
333 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
334 OPM_THROW(std::runtime_error,
335 "Cannot initialise: when no explicit RVVD table is given, \n"
336 "datum depth must be at the gas-oil-contact. "
337 "In EQUIL region " << (i + 1) <<
" (counting from 1), this does not hold.");
339 const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
340 const double T_contact = 273.15 + 20;
341 rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact, T_contact));
345 for (
size_t i = 0; i < rec.size(); ++i) {
346 rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
351 calcPressSatRsRv(eqlmap, rec, props, G, grav);
357 typedef std::vector<double> Vec;
358 typedef std::vector<Vec> PVec;
360 const PVec& press()
const {
return pp_; }
361 const PVec& saturation()
const {
return sat_; }
362 const Vec& rs()
const {
return rs_; }
363 const Vec& rv()
const {
return rv_; }
366 typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
367 typedef EquilReg<RhoCalc> EqReg;
369 std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
370 std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
378 template <
class RMap,
class Gr
id>
380 calcPressSatRsRv(
const RMap& reg ,
381 const std::vector< EquilRecord >& rec ,
386 for (
const auto& r : reg.activeRegions()) {
387 const auto& cells = reg.cells(r);
390 OpmLog::warning(
"Equilibration region " + std::to_string(r + 1)
391 +
" has no active cells");
394 const int repcell = *cells.begin();
396 const RhoCalc calc(props, repcell);
397 const EqReg eqreg(rec[r], calc,
398 rs_func_[r], rv_func_[r],
401 PVec pressures = phasePressures(G, eqreg, cells, grav);
402 const std::vector<double>& temp = temperature(G, eqreg, cells);
404 const PVec sat = phaseSaturations(G, eqreg, cells, props, swat_init_, pressures);
407 for (
int p = 0; p < np; ++p) {
408 copyFromRegion(pressures[p], cells, pp_[p]);
409 copyFromRegion(sat[p], cells, sat_[p]);
411 if (props.
phaseUsage().phase_used[BlackoilPhases::Liquid]
412 && props.
phaseUsage().phase_used[BlackoilPhases::Vapour]) {
413 const int oilpos = props.
phaseUsage().phase_pos[BlackoilPhases::Liquid];
414 const int gaspos = props.
phaseUsage().phase_pos[BlackoilPhases::Vapour];
415 const Vec rs_vals =
computeRs(G, cells, pressures[oilpos], temp, *(rs_func_[r]), sat[gaspos]);
416 const Vec rv_vals =
computeRs(G, cells, pressures[gaspos], temp, *(rv_func_[r]), sat[oilpos]);
417 copyFromRegion(rs_vals, cells, rs_);
418 copyFromRegion(rv_vals, cells, rv_);
423 template <
class CellRangeType>
424 void copyFromRegion(
const Vec& source,
425 const CellRangeType& cells,
428 auto s = source.begin();
429 auto c = cells.begin();
430 const auto e = cells.end();
431 for (; c != e; ++c, ++s) {
432 destination[*c] = *s;
441 #include <opm/core/simulator/initStateEquil_impl.hpp>
443 #endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED
Abstract base class for blackoil fluid and reservoir properties.
Definition: BlackoilPropertiesInterface.hpp:37
virtual int numPhases() const =0
virtual PhaseUsage phaseUsage() const =0
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)
Compute initial Rs values.
Definition: initStateEquil_impl.hpp:822