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