initStateEquil.hpp
Go to the documentation of this file.
1 /*
2  Copyright 2014 SINTEF ICT, Applied Mathematics.
3  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2015 NTNU
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 #ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
23 #define OPM_INITSTATEEQUIL_HEADER_INCLUDED
24 
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>
41 
42 #include <array>
43 #include <cassert>
44 #include <utility>
45 #include <vector>
46 
52 struct UnstructuredGrid;
53 
54 namespace Opm
55 {
56 
74  template<class Grid>
75  void initStateEquil(const Grid& grid,
76  const BlackoilPropertiesInterface& props,
77  const Opm::Deck& deck,
78  const Opm::EclipseState& eclipseState,
79  const double gravity,
80  BlackoilState& state,
81  bool applySwatInit = true);
82 
83 
91  namespace EQUIL {
92 
131  template <class Grid, class Region, class CellRange>
132  std::vector< std::vector<double> >
133  phasePressures(const Grid& G,
134  const Region& reg,
135  const CellRange& cells,
136  const double grav = unit::gravity);
137 
138 
139 
164  template <class Grid, class Region, class CellRange>
165  std::vector< std::vector<double> >
166  phaseSaturations(const Grid& grid,
167  const Region& reg,
168  const CellRange& cells,
169  BlackoilPropertiesFromDeck& props,
170  const std::vector<double> swat_init,
171  std::vector< std::vector<double> >& phase_pressures);
172 
173 
174 
193  template <class Grid, 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);
200 
201  namespace DeckDependent {
202  inline
203  std::vector<EquilRecord>
204  getEquil(const Opm::EclipseState& state)
205  {
206  const auto& init = state.getInitConfig();
207 
208  if( !init.hasEquil() ) {
209  OPM_THROW(std::domain_error, "Deck does not provide equilibration data.");
210  }
211 
212  const auto& equil = init.getEquil();
213  return { equil.begin(), equil.end() };
214  }
215 
216  template<class Grid>
217  inline
218  std::vector<int>
219  equilnum(const Opm::Deck& deck,
220  const Opm::EclipseState& eclipseState,
221  const Grid& G )
222  {
223  std::vector<int> eqlnum;
224  if (deck.hasKeyword("EQLNUM")) {
225  const int nc = UgGridHelpers::numCells(G);
226  eqlnum.resize(nc);
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;
233  }
234  }
235  else {
236  // No explicit equilibration region.
237  // All cells in region zero.
238  eqlnum.assign(UgGridHelpers::numCells(G), 0);
239  }
240 
241  return eqlnum;
242  }
243 
244 
245  class InitialStateComputer {
246  public:
247  template<class Grid>
248  InitialStateComputer(BlackoilPropertiesInterface& props,
249  const Opm::Deck& deck,
250  const Opm::EclipseState& eclipseState,
251  const Grid& G ,
252  const double grav = unit::gravity,
253  const std::vector<double>& swat_init = {}
254  )
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)
262 
263  {
264  // Get the equilibration records.
265  const std::vector<EquilRecord> rec = getEquil(eclipseState);
266  const auto& tables = eclipseState.getTableManager();
267  // Create (inverse) region mapping.
268  const RegionMapping<> eqlmap(equilnum(deck, eclipseState, G));
269 
270  // Create Rs functions.
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())
276  {
277  rs_func_.push_back(std::shared_ptr<Miscibility::RsVD>());
278  continue;
279  }
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.");
284  }
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,
289  cell,
290  depthColumn , rsColumn));
291  } else {
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.");
297  }
298  const double p_contact = rec[i].datumDepthPressure();
299  const double T_contact = 273.15 + 20; // standard temperature for now
300  rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact, T_contact));
301  }
302  }
303  } else {
304  for (size_t i = 0; i < rec.size(); ++i) {
305  rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
306  }
307  }
308 
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())
314  {
315  rv_func_.push_back(std::shared_ptr<Miscibility::RvVD>());
316  continue;
317  }
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.");
322  }
323 
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();
327 
328  rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props,
329  cell,
330  depthColumn , rvColumn));
331 
332  } else {
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.");
338  }
339  const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
340  const double T_contact = 273.15 + 20; // standard temperature for now
341  rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact, T_contact));
342  }
343  }
344  } else {
345  for (size_t i = 0; i < rec.size(); ++i) {
346  rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
347  }
348  }
349 
350  // Compute pressures, saturations, rs and rv factors.
351  calcPressSatRsRv(eqlmap, rec, props, G, grav);
352 
353  // Modify oil pressure in no-oil regions so that the pressures of present phases can
354  // be recovered from the oil pressure and capillary relations.
355  }
356 
357  typedef std::vector<double> Vec;
358  typedef std::vector<Vec> PVec; // One per phase.
359 
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_; }
364 
365  private:
366  typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
367  typedef EquilReg<RhoCalc> EqReg;
368 
369  std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
370  std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
371 
372  PVec pp_;
373  PVec sat_;
374  Vec rs_;
375  Vec rv_;
376  Vec swat_init_;
377 
378  template <class RMap, class Grid>
379  void
380  calcPressSatRsRv(const RMap& reg ,
381  const std::vector< EquilRecord >& rec ,
383  const Grid& G ,
384  const double grav)
385  {
386  for (const auto& r : reg.activeRegions()) {
387  const auto& cells = reg.cells(r);
388  if (cells.empty())
389  {
390  OpmLog::warning("Equilibration region " + std::to_string(r + 1)
391  + " has no active cells");
392  continue;
393  }
394  const int repcell = *cells.begin();
395 
396  const RhoCalc calc(props, repcell);
397  const EqReg eqreg(rec[r], calc,
398  rs_func_[r], rv_func_[r],
399  props.phaseUsage());
400 
401  PVec pressures = phasePressures(G, eqreg, cells, grav);
402  const std::vector<double>& temp = temperature(G, eqreg, cells);
403 
404  const PVec sat = phaseSaturations(G, eqreg, cells, props, swat_init_, pressures);
405 
406  const int np = props.numPhases();
407  for (int p = 0; p < np; ++p) {
408  copyFromRegion(pressures[p], cells, pp_[p]);
409  copyFromRegion(sat[p], cells, sat_[p]);
410  }
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_);
419  }
420  }
421  }
422 
423  template <class CellRangeType>
424  void copyFromRegion(const Vec& source,
425  const CellRangeType& cells,
426  Vec& destination)
427  {
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;
433  }
434  }
435 
436  };
437  } // namespace DeckDependent
438  } // namespace EQUIL
439 } // namespace Opm
440 
441 #include <opm/core/simulator/initStateEquil_impl.hpp>
442 
443 #endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED
void initStateEquil(const Grid &grid, const BlackoilPropertiesInterface &props, const Opm::Deck &deck, const Opm::EclipseState &eclipseState, const double gravity, BlackoilState &state, bool applySwatInit=true)
Compute initial state by an equilibration procedure.
Definition: AnisotropicEikonal.cpp:446
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