initState_impl.hpp
Go to the documentation of this file.
1 /*
2  Copyright 2012 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_INITSTATE_IMPL_HEADER_INCLUDED
21 #define OPM_INITSTATE_IMPL_HEADER_INCLUDED
22 
23 
24 #include <opm/core/utility/parameters/ParameterGroup.hpp>
25 #include <opm/core/grid.h>
26 #include <opm/core/grid/GridHelpers.hpp>
27 #include <opm/common/ErrorMacros.hpp>
28 #include <opm/core/utility/MonotCubicInterpolator.hpp>
29 #include <opm/core/utility/Units.hpp>
30 #include <opm/core/props/IncompPropertiesInterface.hpp>
31 #include <opm/core/props/BlackoilPropertiesInterface.hpp>
32 #include <opm/core/props/phaseUsageFromDeck.hpp>
33 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
34 
35 #include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
36 
37 #include <iostream>
38 #include <cmath>
39 
40 namespace Opm
41 {
42 
43  namespace
44  {
45 #ifdef __clang__
46 #pragma clang diagnostic push
47 #pragma clang diagnostic ignored "-Wunneeded-internal-declaration"
48 #endif /* __clang__ */
49  // Find the cells that are below and above a depth.
50  // TODO: add 'anitialiasing', obtaining a more precise split
51  // by f. ex. subdividing cells cut by the split depths.
52  template<class T>
53  void cellsBelowAbove(int number_of_cells,
54  T begin_cell_centroids,
55  int dimension,
56  const double depth,
57  std::vector<int>& below,
58  std::vector<int>& above)
59  {
60  below.reserve(number_of_cells);
61  above.reserve(number_of_cells);
62  for (int c = 0; c < number_of_cells; ++c) {
63  const double z = UgGridHelpers
64  ::getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c,
65  dimension),
66  dimension-1);
67  if (z > depth) {
68  below.push_back(c);
69  } else {
70  above.push_back(c);
71  }
72  }
73  }
74 #ifdef __clang__
75 #pragma clang diagnostic pop
76 #endif /* __clang__ */
77 
78  enum WaterInit { WaterBelow, WaterAbove };
79 
80  // Initialize saturations so that there is water below woc,
81  // and oil above.
82  // If invert is true, water is instead above, oil below.
83  template <class T, class Props, class State>
84  void initWaterOilContact(int number_of_cells,
85  T begin_cell_centroids,
86  int dimensions,
87  const Props& props,
88  const double woc,
89  const WaterInit waterinit,
90  State& state)
91  {
92  // Find out which cells should have water and which should have oil.
93  std::vector<int> water;
94  std::vector<int> oil;
95  // Assuming that water should go below the woc, but warning if oil is heavier.
96  // if (props.density()[0] < props.density()[1]) {
97  // std::cout << "*** warning: water density is less than oil density, "
98  // "but initialising water below woc." << std::endl;
99  // }
100  switch (waterinit) {
101  case WaterBelow:
102  cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions,
103  woc, water, oil);
104  break;
105  case WaterAbove:
106  cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions,
107  woc, oil, water);
108  }
109  // Set saturations.
110  state.setFirstSat(oil, props, State::MinSat);
111  state.setFirstSat(water, props, State::MaxSat);
112  }
113 
114 
115  // Initialize hydrostatic pressures depending only on gravity,
116  // (constant) phase densities and a water-oil contact depth.
117  // The pressure difference between points is equal to
118  // delta_p = delta_z * gravity * rho
119  // where rho is the (assumed constant) oil density above the
120  // woc, water density below woc.
121  // Note that even if there is (immobile) water in the oil zone,
122  // it does not contribute to the pressure difference.
123  // Note that by manipulating the densities parameter,
124  // it is possible to initialise with water on top instead of
125  // on the bottom etc.
126  template <class T, class State>
127  void initHydrostaticPressure(int number_of_cells,
128  T begin_cell_centroids,
129  int dimensions,
130  const double* densities,
131  const double woc,
132  const double gravity,
133  const double datum_z,
134  const double datum_p,
135  State& state)
136  {
137  std::vector<double>& p = state.pressure();
138  // Compute pressure at woc
139  const double rho_datum = datum_z > woc ? densities[0] : densities[1];
140  const double woc_p = datum_p + (woc - datum_z)*gravity*rho_datum;
141  for (int c = 0; c < number_of_cells; ++c) {
142  // Compute pressure as delta from woc pressure.
143  const double z = UgGridHelpers::
144  getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
145  dimensions-1);
146  const double rho = z > woc ? densities[0] : densities[1];
147  p[c] = woc_p + (z - woc)*gravity*rho;
148  }
149  }
150 
151 
152  // Facade to initHydrostaticPressure() taking a property object,
153  // for similarity to initHydrostaticPressure() for compressible fluids.
154  template <class T, class State>
155  void initHydrostaticPressure(int number_of_cells,
156  T begin_cell_centroids,
157  int dimensions,
158  const IncompPropertiesInterface& props,
159  const double woc,
160  const double gravity,
161  const double datum_z,
162  const double datum_p,
163  State& state)
164  {
165  const double* densities = props.density();
166  initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
167  densities, woc, gravity, datum_z, datum_p, state);
168  }
169 
170 
171 
172  // Helper functor for initHydrostaticPressure() for compressible fluids.
173  struct Density
174  {
175  const BlackoilPropertiesInterface& props_;
176  Density(const BlackoilPropertiesInterface& props) : props_(props) {}
177  double operator()(const double pressure, const double temperature, const int phase)
178  {
179  assert(props_.numPhases() == 2);
180  const double surfvol[2][2] = { { 1.0, 0.0 },
181  { 0.0, 1.0 } };
182  // We do not handle multi-region PVT/EQUIL at this point.
183  const int cell = 0;
184  double A[4] = { 0.0 };
185  props_.matrix(1, &pressure, &temperature, surfvol[phase], &cell, A, 0);
186  double rho[2] = { 0.0 };
187  props_.density(1, A, &cell, rho);
188  return rho[phase];
189  }
190  };
191 
192  // Initialize hydrostatic pressures depending only on gravity,
193  // phase densities that may vary with pressure and a water-oil
194  // contact depth. The pressure ODE is given as
195  // \grad p = \rho g \grad z
196  // where rho is the oil density above the woc, water density
197  // below woc. Note that even if there is (immobile) water in
198  // the oil zone, it does not contribute to the pressure there.
199  template <class T, class State>
200  void initHydrostaticPressure(int number_of_cells,
201  T begin_cell_centroids,
202  int dimensions,
203  const BlackoilPropertiesInterface& props,
204  const double woc,
205  const double gravity,
206  const double datum_z,
207  const double datum_p,
208  State& state)
209  {
210  assert(props.numPhases() == 2);
211 
212  // Obtain max and min z for which we will need to compute p.
213  double zlim[2] = { 1e100, -1e100 };
214  for (int c = 0; c < number_of_cells; ++c) {
215  const double z = UgGridHelpers::
216  getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
217  dimensions-1);;
218  zlim[0] = std::min(zlim[0], z);
219  zlim[1] = std::max(zlim[1], z);
220  }
221 
222  // We'll use a minimum stepsize of 1/100 of the z range.
223  const double hmin = (zlim[1] - zlim[0])/100.0;
224 
225  // Store p(z) results in an associative array.
226  std::map<double, double> press_by_z;
227  press_by_z[datum_z] = datum_p;
228 
229  // Set up density evaluator.
230  Density rho(props);
231 
232  // Solve the ODE from datum_z to woc.
233  int phase = (datum_z > woc) ? 0 : 1;
234  int num_steps = int(std::ceil(std::fabs(woc - datum_z)/hmin));
235  double pval = datum_p;
236  double Tval = 273.15 + 20; // standard temperature
237  double zval = datum_z;
238  double h = (woc - datum_z)/double(num_steps);
239  for (int i = 0; i < num_steps; ++i) {
240  zval += h;
241  const double dp = rho(pval, Tval, phase)*gravity;
242  pval += h*dp;
243  press_by_z[zval] = pval;
244  }
245  double woc_p = pval;
246 
247  // Solve the ODE from datum_z to the end of the interval.
248  double z_end = (datum_z > woc) ? zlim[1] : zlim[0];
249  num_steps = int(std::ceil(std::fabs(z_end - datum_z)/hmin));
250  pval = datum_p;
251  zval = datum_z;
252  h = (z_end - datum_z)/double(num_steps);
253  for (int i = 0; i < num_steps; ++i) {
254  zval += h;
255  const double dp = rho(pval, Tval, phase)*gravity;
256  pval += h*dp;
257  press_by_z[zval] = pval;
258  }
259 
260  // Solve the ODE from woc to the other end of the interval.
261  // Switching phase and z_end.
262  phase = (phase + 1) % 2;
263  z_end = (datum_z > woc) ? zlim[0] : zlim[1];
264  pval = woc_p;
265  zval = woc;
266  h = (z_end - datum_z)/double(num_steps);
267  for (int i = 0; i < num_steps; ++i) {
268  zval += h;
269  const double dp = rho(pval, Tval, phase)*gravity;
270  pval += h*dp;
271  press_by_z[zval] = pval;
272  }
273 
274  // Create monotone spline for interpolating solution.
275  std::vector<double> zv, pv;
276  zv.reserve(press_by_z.size());
277  pv.reserve(press_by_z.size());
278  for (std::map<double, double>::const_iterator it = press_by_z.begin();
279  it != press_by_z.end(); ++it) {
280  zv.push_back(it->first);
281  pv.push_back(it->second);
282  }
283  MonotCubicInterpolator press(zv, pv);
284 
285  // Evaluate pressure at each cell centroid.
286  std::vector<double>& p = state.pressure();
287  for (int c = 0; c < number_of_cells; ++c) {
288  const double z = UgGridHelpers::
289  getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
290  dimensions-1);
291  p[c] = press(z);
292  }
293  }
294 
295  // Initialize face pressures to distance-weighted average of adjacent cell pressures.
296  template <class State, class FaceCells, class FCI, class CCI>
297  void initFacePressure(int dimensions,
298  int number_of_faces,
299  FaceCells face_cells,
300  FCI begin_face_centroids,
301  CCI begin_cell_centroids,
302  State& state)
303  {
304  const std::vector<double>& cp = state.pressure();
305  std::vector<double>& fp = state.facepressure();
306  for (int f = 0; f < number_of_faces; ++f) {
307  double dist[2] = { 0.0, 0.0 };
308  double press[2] = { 0.0, 0.0 };
309  int bdy_idx = -1;
310  for (int j = 0; j < 2; ++j) {
311  const int c = face_cells(f, j);
312  if (c >= 0) {
313  dist[j] = 0.0;
314  for (int dd = 0; dd < dimensions; ++dd) {
315  double diff = UgGridHelpers::
316  getCoordinate(UgGridHelpers::increment(begin_face_centroids, f,
317  dimensions),
318  dd)
319  - UgGridHelpers::
320  getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c,
321  dimensions),
322  dd);
323  dist[j] += diff*diff;
324  }
325  dist[j] = std::sqrt(dist[j]);
326  press[j] = cp[c];
327  } else {
328  bdy_idx = j;
329  }
330  }
331  if (bdy_idx == -1) {
332  fp[f] = press[0]*(dist[1]/(dist[0] + dist[1])) + press[1]*(dist[0]/(dist[0] + dist[1]));
333  } else {
334  fp[f] = press[(bdy_idx + 1) % 2];
335  }
336  }
337  }
338 
339 
340  } // anonymous namespace
341 
342 
344  template <class State>
345  void initStateBasic(const UnstructuredGrid& grid,
346  const IncompPropertiesInterface& props,
347  const parameter::ParameterGroup& param,
348  const double gravity,
349  State& state)
350  {
351  initStateBasic(grid.number_of_cells, grid.global_cell, grid.cartdims,
352  grid.number_of_faces, UgGridHelpers::faceCells(grid),
353  grid.face_centroids, grid.cell_centroids,
354  grid.dimensions, props, param,
355  gravity, state);
356  }
357 
358  template <class FaceCells, class CCI, class FCI, class State>
359  void initStateBasic(int number_of_cells,
360  const int* global_cell,
361  const int* cartdims,
362  int number_of_faces,
363  FaceCells face_cells,
364  FCI begin_face_centroids,
365  CCI begin_cell_centroids,
366  int dimensions,
367  const IncompPropertiesInterface& props,
368  const parameter::ParameterGroup& param,
369  const double gravity,
370  State& state)
371 {
372  const int num_phases = props.numPhases();
373  if (num_phases != 2) {
374  OPM_THROW(std::runtime_error, "initStateTwophaseBasic(): currently handling only two-phase scenarios.");
375  }
376  state.init(number_of_cells, number_of_faces, num_phases);
377  const int num_cells = props.numCells();
378  // By default: initialise water saturation to minimum everywhere.
379  std::vector<int> all_cells(num_cells);
380  for (int i = 0; i < num_cells; ++i) {
381  all_cells[i] = i;
382  }
383  state.setFirstSat(all_cells, props, State::MinSat);
384  const bool convection_testcase = param.getDefault("convection_testcase", false);
385  const bool segregation_testcase = param.getDefault("segregation_testcase", false);
386  if (convection_testcase) {
387  // Initialise water saturation to max in the 'left' part.
388  std::vector<int> left_cells;
389  left_cells.reserve(num_cells/2);
390  for (int cell = 0; cell < num_cells; ++cell) {
391  const int gc = global_cell == 0 ? cell : global_cell[cell];
392  bool left = (gc % cartdims[0]) < cartdims[0]/2;
393  if (left) {
394  left_cells.push_back(cell);
395  }
396  }
397  state.setFirstSat(left_cells, props, State::MaxSat);
398  const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
399  std::fill(state.pressure().begin(), state.pressure().end(), init_p);
400  } else if (segregation_testcase) {
401  // Warn against error-prone usage.
402  if (gravity == 0.0) {
403  std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl;
404  }
405  if (cartdims[2] <= 1) {
406  std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl;
407  }
408  // Initialise water saturation to max *above* water-oil contact.
409  const double woc = param.get<double>("water_oil_contact");
410  initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
411  props, woc, WaterAbove, state);
412  // Initialise pressure to hydrostatic state.
413  const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
414  double dens[2] = { props.density()[1], props.density()[0] };
415  initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
416  dens, woc, gravity, woc, ref_p, state);
417  } else if (param.has("water_oil_contact")) {
418  // Warn against error-prone usage.
419  if (gravity == 0.0) {
420  std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
421  }
422  if (cartdims[2] <= 1) {
423  std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
424  }
425  // Initialise water saturation to max below water-oil contact.
426  const double woc = param.get<double>("water_oil_contact");
427  initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
428  props, woc, WaterBelow, state);
429  // Initialise pressure to hydrostatic state.
430  const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
431  initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
432  props.density(), woc, gravity, woc, ref_p, state);
433  } else if (param.has("init_saturation")) {
434  // Initialise water saturation to init_saturation parameter.
435  const double init_saturation = param.get<double>("init_saturation");
436  for (int cell = 0; cell < num_cells; ++cell) {
437  state.saturation()[2*cell] = init_saturation;
438  state.saturation()[2*cell + 1] = 1.0 - init_saturation;
439  }
440  // Initialise pressure to hydrostatic state.
441  const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
442  const double rho = props.density()[0]*init_saturation + props.density()[1]*(1.0 - init_saturation);
443  const double dens[2] = { rho, rho };
444  const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1);
445  initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
446  dens, ref_z, gravity, ref_z, ref_p, state);
447  } else {
448  // Use default: water saturation is minimum everywhere.
449  // Initialise pressure to hydrostatic state.
450  const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
451  const double rho = props.density()[1];
452  const double dens[2] = { rho, rho };
453  const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1);
454  initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
455  dens, ref_z, gravity, ref_z, ref_p, state);
456  }
457 
458  // Finally, init face pressures.
459  initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
460  begin_cell_centroids, state);
461  }
462 
463 
465  template <class State>
466  void initStateBasic(const UnstructuredGrid& grid,
467  const BlackoilPropertiesInterface& props,
468  const parameter::ParameterGroup& param,
469  const double gravity,
470  State& state)
471  {
472  initStateBasic(grid.number_of_cells, grid.global_cell, grid.cartdims,
473  grid.number_of_faces, UgGridHelpers::faceCells(grid),
474  grid.face_centroids, grid.cell_centroids, grid.dimensions,
475  props, param, gravity, state);
476  }
477 
478  template <class FaceCells, class FCI, class CCI, class State>
479  void initStateBasic(int number_of_cells,
480  const int* global_cell,
481  const int* cartdims,
482  int number_of_faces,
483  FaceCells face_cells,
484  FCI begin_face_centroids,
485  CCI begin_cell_centroids,
486  int dimensions,
487  const BlackoilPropertiesInterface& props,
488  const parameter::ParameterGroup& param,
489  const double gravity,
490  State& state)
491  {
492  // TODO: Refactor to exploit similarity with IncompProp* case.
493  const int num_phases = props.numPhases();
494  if (num_phases != 2) {
495  OPM_THROW(std::runtime_error, "initStateTwophaseBasic(): currently handling only two-phase scenarios.");
496  }
497  state.init(number_of_cells, number_of_faces, num_phases);
498  const int num_cells = props.numCells();
499  // By default: initialise water saturation to minimum everywhere.
500  std::vector<int> all_cells(num_cells);
501  for (int i = 0; i < num_cells; ++i) {
502  all_cells[i] = i;
503  }
504  state.setFirstSat(all_cells, props, State::MinSat);
505  const bool convection_testcase = param.getDefault("convection_testcase", false);
506  if (convection_testcase) {
507  // Initialise water saturation to max in the 'left' part.
508  std::vector<int> left_cells;
509  left_cells.reserve(num_cells/2);
510  for (int cell = 0; cell < num_cells; ++cell) {
511  const int gc = global_cell == 0 ? cell : global_cell[cell];
512  bool left = (gc % cartdims[0]) < cartdims[0]/2;
513  if (left) {
514  left_cells.push_back(cell);
515  }
516  }
517  state.setFirstSat(left_cells, props, State::MaxSat);
518  const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
519  std::fill(state.pressure().begin(), state.pressure().end(), init_p);
520  } else if (param.has("water_oil_contact")) {
521  // Warn against error-prone usage.
522  if (gravity == 0.0) {
523  std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
524  }
525  if (cartdims[2] <= 1) {
526  std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
527  }
528  // Initialise water saturation to max below water-oil contact.
529  const double woc = param.get<double>("water_oil_contact");
530  initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
531  props, woc, WaterBelow, state);
532  // Initialise pressure to hydrostatic state.
533  const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
534  initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
535  props, woc, gravity, woc, ref_p, state);
536  } else {
537  // Use default: water saturation is minimum everywhere.
538  // Initialise pressure to hydrostatic state.
539  const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
540  const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1);
541  const double woc = -1e100;
542  initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
543  props, woc, gravity, ref_z, ref_p, state);
544  }
545 
546  // Finally, init face pressures.
547  initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
548  begin_cell_centroids, state);
549  }
550 
551 
553  template <class Props, class State>
554  void initStateFromDeck(const UnstructuredGrid& grid,
555  const Props& props,
556  Opm::DeckConstPtr deck,
557  const double gravity,
558  State& state)
559  {
560  initStateFromDeck(grid.number_of_cells, grid.global_cell,
561  grid.number_of_faces, UgGridHelpers::faceCells(grid),
562  grid.face_centroids, grid.cell_centroids, grid.dimensions,
563  props, deck, gravity, state);
564  }
566  template <class FaceCells, class FCI, class CCI, class Props, class State>
567  void initStateFromDeck(int number_of_cells,
568  const int* global_cell,
569  int number_of_faces,
570  FaceCells face_cells,
571  FCI begin_face_centroids,
572  CCI begin_cell_centroids,
573  int dimensions,
574  const Props& props,
575  Opm::DeckConstPtr deck,
576  const double gravity,
577  State& state)
578  {
579  const int num_phases = props.numPhases();
580  const PhaseUsage pu = phaseUsageFromDeck(deck);
581  if (num_phases != pu.num_phases) {
582  OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, "
583  "found " << pu.num_phases << " phases in deck.");
584  }
585  state.init(number_of_cells, number_of_faces, num_phases);
586  if (deck->hasKeyword("EQUIL") && deck->hasKeyword("PRESSURE")) {
587  OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
588  "condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
589  }
590 
591  if (deck->hasKeyword("EQUIL")) {
592  if (num_phases != 2) {
593  OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
594  }
595  if (pu.phase_used[BlackoilPhases::Vapour]) {
596  OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
597  }
598  // Set saturations depending on oil-water contact.
599  EquilWrapper equil(deck->getKeyword("EQUIL"));
600  if (equil.numRegions() != 1) {
601  OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
602  }
603  const double woc = equil.waterOilContactDepth(0);
604  initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
605  props, woc, WaterBelow, state);
606  // Set pressure depending on densities and depths.
607  const double datum_z = equil.datumDepth(0);
608  const double datum_p = equil.datumDepthPressure(0);
609  initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
610  props, woc, gravity, datum_z, datum_p, state);
611  } else if (deck->hasKeyword("PRESSURE")) {
612  // Set saturations from SWAT/SGAS, pressure from PRESSURE.
613  std::vector<double>& s = state.saturation();
614  std::vector<double>& p = state.pressure();
615  const std::vector<double>& p_deck = deck->getKeyword("PRESSURE")->getSIDoubleData();
616  const int num_cells = number_of_cells;
617  if (num_phases == 2) {
618  if (!pu.phase_used[BlackoilPhases::Aqua]) {
619  // oil-gas: we require SGAS
620  if (!deck->hasKeyword("SGAS")) {
621  OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS keyword in 2-phase init");
622  }
623  const std::vector<double>& sg_deck = deck->getKeyword("SGAS")->getSIDoubleData();
624  const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
625  const int opos = pu.phase_pos[BlackoilPhases::Liquid];
626  for (int c = 0; c < num_cells; ++c) {
627  int c_deck = (global_cell == NULL) ? c : global_cell[c];
628  s[2*c + gpos] = sg_deck[c_deck];
629  s[2*c + opos] = 1.0 - sg_deck[c_deck];
630  p[c] = p_deck[c_deck];
631  }
632  } else {
633  // water-oil or water-gas: we require SWAT
634  if (!deck->hasKeyword("SWAT")) {
635  OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init");
636  }
637  const std::vector<double>& sw_deck = deck->getKeyword("SWAT")->getSIDoubleData();
638  const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
639  const int nwpos = (wpos + 1) % 2;
640  for (int c = 0; c < num_cells; ++c) {
641  int c_deck = (global_cell == NULL) ? c : global_cell[c];
642  s[2*c + wpos] = sw_deck[c_deck];
643  s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
644  p[c] = p_deck[c_deck];
645  }
646  }
647  } else if (num_phases == 3) {
648  const bool has_swat_sgas = deck->hasKeyword("SWAT") && deck->hasKeyword("SGAS");
649  if (!has_swat_sgas) {
650  OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
651  }
652  const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
653  const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
654  const int opos = pu.phase_pos[BlackoilPhases::Liquid];
655  const std::vector<double>& sw_deck = deck->getKeyword("SWAT")->getSIDoubleData();
656  const std::vector<double>& sg_deck = deck->getKeyword("SGAS")->getSIDoubleData();
657  for (int c = 0; c < num_cells; ++c) {
658  int c_deck = (global_cell == NULL) ? c : global_cell[c];
659  s[3*c + wpos] = sw_deck[c_deck];
660  s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
661  s[3*c + gpos] = sg_deck[c_deck];
662  p[c] = p_deck[c_deck];
663  }
664  } else {
665  OPM_THROW(std::runtime_error, "initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
666  }
667  } else {
668  OPM_THROW(std::runtime_error, "initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS.");
669  }
670 
671  // Finally, init face pressures.
672  initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
673  begin_cell_centroids, state);
674  }
675 
679  template <class Props, class State>
680  void initBlackoilSurfvol(const UnstructuredGrid& grid,
681  const Props& props,
682  State& state)
683  {
684  initBlackoilSurfvol(grid.number_of_cells, props, state);
685  }
686 
687  template <class Props, class State>
688  void initBlackoilSurfvol(int number_of_cells,
689  const Props& props,
690  State& state)
691  {
692  state.surfacevol() = state.saturation();
693  const int np = props.numPhases();
694  const int nc = number_of_cells;
695  std::vector<double> allA(nc*np*np);
696  std::vector<int> allcells(nc);
697  for (int c = 0; c < nc; ++c) {
698  allcells[c] = c;
699  }
700  // Assuming that using the saturation as z argument here does not change
701  // the outcome. This is not guaranteed unless we have only a single phase
702  // per cell.
703  props.matrix(nc, &state.pressure()[0], &state.temperature()[0], &state.surfacevol()[0], &allcells[0], &allA[0], 0);
704  for (int c = 0; c < nc; ++c) {
705  // Using z = As
706  double* z = &state.surfacevol()[c*np];
707  const double* s = &state.saturation()[c*np];
708  const double* A = &allA[c*np*np];
709 
710  for (int row = 0; row < np; ++row) { z[row] = 0.0; }
711 
712  for (int col = 0; col < np; ++col) {
713  for (int row = 0; row < np; ++row) {
714  // Recall: A has column-major ordering.
715  z[row] += A[row + np*col]*s[col];
716  }
717  }
718  }
719  }
723  template <class Props, class State>
724  void initBlackoilSurfvolUsingRSorRV(const UnstructuredGrid& grid,
725  const Props& props,
726  State& state)
727  {
728  initBlackoilSurfvolUsingRSorRV(grid.number_of_cells, props, state);
729  }
730 
731  template <class Props, class State>
732  void initBlackoilSurfvolUsingRSorRV(int number_of_cells,
733  const Props& props,
734  State& state)
735  {
736 
737  if (props.numPhases() != 3) {
738  OPM_THROW(std::runtime_error, "initBlackoilSurfvol() is only supported in three-phase simulations.");
739  }
740  const std::vector<double>& rs = state.gasoilratio();
741  const std::vector<double>& rv = state.rv();
742 
743  //make input for computation of the A matrix
744  state.surfacevol() = state.saturation();
745  const PhaseUsage pu = props.phaseUsage();
746 
747  const int np = props.numPhases();
748 
749  std::vector<double> allA_a(number_of_cells*np*np);
750  std::vector<double> allA_l(number_of_cells*np*np);
751  std::vector<double> allA_v(number_of_cells*np*np);
752 
753  std::vector<int> allcells(number_of_cells);
754  std::vector<double> z_init(number_of_cells*np);
755  std::fill(z_init.begin(),z_init.end(),0.0);
756 
757  for (int c = 0; c < number_of_cells; ++c) {
758  allcells[c] = c;
759  }
760 
761  std::vector<double> capPressures(number_of_cells*np);
762  props.capPress(number_of_cells,&state.saturation()[0],&allcells[0],&capPressures[0],NULL);
763 
764  std::vector<double> Pw(number_of_cells);
765  std::vector<double> Pg(number_of_cells);
766 
767  for (int c = 0; c < number_of_cells; ++c){
768  Pw[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Aqua];
769  Pg[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Vapour];
770  }
771 
772 
773  double z_tmp;
774 
775  // Water phase
776  if(pu.phase_used[BlackoilPhases::Aqua])
777  for (int c = 0; c < number_of_cells ; ++c){
778  for (int p = 0; p < np ; ++p){
779  if (p == BlackoilPhases::Aqua)
780  z_tmp = 1;
781  else
782  z_tmp = 0;
783 
784  z_init[c*np + p] = z_tmp;
785  }
786  }
787  props.matrix(number_of_cells, &state.pressure()[0], &state.temperature()[0], &z_init[0], &allcells[0], &allA_a[0], 0);
788 
789  // Liquid phase
790  if(pu.phase_used[BlackoilPhases::Liquid]){
791  for (int c = 0; c < number_of_cells ; ++c){
792  for (int p = 0; p < np ; ++p){
793  if(p == BlackoilPhases::Vapour){
794  if(state.saturation()[np*c + p] > 0)
795  z_tmp = 1e10;
796  else
797  z_tmp = rs[c];
798  }
799  else if(p == BlackoilPhases::Liquid)
800  z_tmp = 1;
801  else
802  z_tmp = 0;
803 
804  z_init[c*np + p] = z_tmp;
805 
806  }
807  }
808  }
809  props.matrix(number_of_cells, &state.pressure()[0], &state.temperature()[0], &z_init[0], &allcells[0], &allA_l[0], 0);
810 
811  if(pu.phase_used[BlackoilPhases::Vapour]){
812  for (int c = 0; c < number_of_cells ; ++c){
813  for (int p = 0; p < np ; ++p){
814  if(p == BlackoilPhases::Liquid){
815  if(state.saturation()[np*c + p] > 0)
816  z_tmp = 1e10;
817  else
818  z_tmp = rv[c];
819  }
820  else if(p == BlackoilPhases::Vapour)
821  z_tmp = 1;
822  else
823  z_tmp = 0;
824 
825  z_init[c*np + p] = z_tmp;
826 
827  }
828  }
829  }
830  props.matrix(number_of_cells, &state.pressure()[0], &state.temperature()[0], &z_init[0], &allcells[0], &allA_v[0], 0);
831 
832  for (int c = 0; c < number_of_cells; ++c) {
833  // Using z = As
834  double* z = &state.surfacevol()[c*np];
835  const double* s = &state.saturation()[c*np];
836  const double* A_a = &allA_a[c*np*np];
837  const double* A_l = &allA_l[c*np*np];
838  const double* A_v = &allA_v[c*np*np];
839 
840  for (int row = 0; row < np; ++row) { z[row] = 0.0; }
841 
842  for (int col = 0; col < np; ++col) {
843  z[0] += A_a[0 + np*col]*s[col];
844  z[1] += A_l[1 + np*col]*s[col];
845  z[2] += A_v[2 + np*col]*s[col];
846 
847  }
848  double ztmp = z[2];
849  z[2] += z[1]*rs[c];
850  z[1] += ztmp*rv[c];
851 
852  }
853  }
854 
856  template <class Props, class State>
857  void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
858  const Props& props,
859  Opm::DeckConstPtr deck,
860  const double gravity,
861  State& state)
862  {
863  initBlackoilStateFromDeck(grid.number_of_cells, grid.global_cell,
864  grid.number_of_faces, UgGridHelpers::faceCells(grid),
865  grid.face_centroids, grid.cell_centroids,grid.dimensions,
866  props, deck, gravity, state);
867  }
868 
870  template <class FaceCells, class FCI, class CCI, class Props, class State>
871  void initBlackoilStateFromDeck(int number_of_cells,
872  const int* global_cell,
873  int number_of_faces,
874  FaceCells face_cells,
875  FCI begin_face_centroids,
876  CCI begin_cell_centroids,
877  int dimensions,
878  const Props& props,
879  Opm::DeckConstPtr deck,
880  const double gravity,
881  State& state)
882  {
883  initStateFromDeck(number_of_cells, global_cell, number_of_faces,
884  face_cells, begin_face_centroids, begin_cell_centroids,
885  dimensions, props, deck, gravity, state);
886  if (deck->hasKeyword("RS")) {
887  const std::vector<double>& rs_deck = deck->getKeyword("RS")->getSIDoubleData();
888  const int num_cells = number_of_cells;
889  for (int c = 0; c < num_cells; ++c) {
890  int c_deck = (global_cell == NULL) ? c : global_cell[c];
891  state.gasoilratio()[c] = rs_deck[c_deck];
892  }
893  initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
894  computeSaturation(props,state);
895  } else if (deck->hasKeyword("RV")){
896  const std::vector<double>& rv_deck = deck->getKeyword("RV")->getSIDoubleData();
897  const int num_cells = number_of_cells;
898  for (int c = 0; c < num_cells; ++c) {
899  int c_deck = (global_cell == NULL) ? c : global_cell[c];
900  state.rv()[c] = rv_deck[c_deck];
901  }
902  initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
903  computeSaturation(props,state);
904  }
905  else {
906  state.gasoilratio() = std::vector<double>(number_of_cells, 0.0);
907  state.rv() = std::vector<double>(number_of_cells, 0.0);
908  initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
909  computeSaturation(props,state);
910  }
911  }
912 
913 } // namespace Opm
914 
915 
916 #endif // OPM_INITSTATE_IMPL_HEADER_INCLUDED
Opm::IncompPropertiesInterface
Definition: IncompPropertiesInterface.hpp:35
Opm::IncompPropertiesInterface::numPhases
virtual int numPhases() const =0
UnstructuredGrid
Definition: grid.h:98
Opm::UgGridHelpers::getCoordinate
double getCoordinate(T *cc, int i)
Get the i-th corrdinate of a centroid.
Definition: GridHelpers.hpp:337
Opm::UgGridHelpers::faceCells
FaceCellTraits< UnstructuredGrid >::Type faceCells(const UnstructuredGrid &grid)
Get the face to cell mapping of a grid.
Opm
Definition: AnisotropicEikonal.hpp:43
ParameterGroup.hpp
Opm::BlackoilPhases::Vapour
Definition: BlackoilPhases.hpp:32
Opm::parameter::ParameterGroup
Definition: ParameterGroup.hpp:109
MonotCubicInterpolator.hpp
Opm::PhaseUsage::phase_used
int phase_used[MaxNumPhases]
Definition: BlackoilPhases.hpp:39
Opm::IncompPropertiesInterface::density
virtual const double * density() const =0
Opm::parameter::ParameterGroup::getDefault
T getDefault(const std::string &name, const T &default_value) const
This method is used to read a parameter from the parameter group.
Definition: ParameterGroup_impl.hpp:226
Opm::BlackoilPropertiesInterface::numPhases
virtual int numPhases() const =0
Opm::Details::PhaseUsed::water
bool water(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:252
Opm::parameter::ParameterGroup::get
T get(const std::string &name) const
This method is used to read a parameter from the parameter group.
Definition: ParameterGroup_impl.hpp:173
Opm::unit::barsa
const double barsa
Definition: Units.hpp:136
Opm::BlackoilPropertiesInterface
Definition: BlackoilPropertiesInterface.hpp:37
Opm::PhaseUsage::phase_pos
int phase_pos[MaxNumPhases]
Definition: BlackoilPhases.hpp:40
UnstructuredGrid::number_of_faces
int number_of_faces
Definition: grid.h:111
Opm::PhaseUsage::num_phases
int num_phases
Definition: BlackoilPhases.hpp:38
UnstructuredGrid::number_of_cells
int number_of_cells
Definition: grid.h:109
UnstructuredGrid::face_centroids
double * face_centroids
Definition: grid.h:168
Opm::UgGridHelpers::dimensions
int dimensions(const UnstructuredGrid &grid)
Get the dimensions of a grid.
IncompPropertiesInterface.hpp
Opm::initBlackoilSurfvol
void initBlackoilSurfvol(const UnstructuredGrid &grid, const Props &props, State &state)
Definition: initState_impl.hpp:680
UnstructuredGrid::dimensions
int dimensions
Definition: grid.h:106
Opm::parameter::ParameterGroup::has
bool has(const std::string &name) const
This method checks if there is something with name name in the parameter gropup.
Opm::unit::gravity
const double gravity
Definition: Units.hpp:120
UnstructuredGrid::cell_centroids
double * cell_centroids
Definition: grid.h:192
Opm::computeSaturation
void computeSaturation(const BlackoilPropertiesInterface &props, BlackoilState &state)
Computes saturation from surface volume densities.
Opm::UgGridHelpers::increment
T * increment(T *cc, int i, int dim)
Increment an iterator over an array that reresents a dense row-major matrix with dims columns...
Definition: GridHelpers.hpp:318
GridHelpers.hpp
Opm::BlackoilPropertiesInterface::numCells
virtual int numCells() const =0
Opm::BlackoilPhases::Liquid
Definition: BlackoilPhases.hpp:32
UnstructuredGrid::cartdims
int cartdims[3]
Definition: grid.h:227
grid
const UnstructuredGrid & grid
Definition: ColumnExtract.hpp:31
UnstructuredGrid::global_cell
int * global_cell
Definition: grid.h:214
Opm::Details::PhaseUsed::oil
bool oil(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:258
Opm::initStateFromDeck
void initStateFromDeck(const UnstructuredGrid &grid, const Props &props, Opm::DeckConstPtr deck, const double gravity, State &state)
Initialize a state from input deck.
Definition: initState_impl.hpp:554
Opm::initBlackoilSurfvolUsingRSorRV
void initBlackoilSurfvolUsingRSorRV(const UnstructuredGrid &grid, const Props &props, State &state)
Definition: initState_impl.hpp:724
Opm::IncompPropertiesInterface::numCells
virtual int numCells() const =0
Opm::phaseUsageFromDeck
PhaseUsage phaseUsageFromDeck(Opm::EclipseStateConstPtr eclipseState)
Definition: phaseUsageFromDeck.hpp:35
Opm::initStateBasic
void initStateBasic(const UnstructuredGrid &grid, const IncompPropertiesInterface &props, const parameter::ParameterGroup &param, const double gravity, State &state)
Initialize a twophase state from parameters.
Definition: initState_impl.hpp:345
BlackoilPropertiesInterface.hpp
Opm::BlackoilPhases::Aqua
Definition: BlackoilPhases.hpp:32
Opm::Equil::temperature
std::vector< double > temperature(const Grid &, const Region &, const CellRange &cells)
Definition: initStateEquil_impl.hpp:667
Opm::PhaseUsage
Definition: BlackoilPhases.hpp:36
Opm::initBlackoilStateFromDeck
void initBlackoilStateFromDeck(const UnstructuredGrid &grid, const Props &props, Opm::DeckConstPtr deck, const double gravity, State &state)
Initialize a blackoil state from input deck.
Definition: initState_impl.hpp:857
Units.hpp