ReservoirPropertyCommon_impl.hpp
1 //===========================================================================
2 //
3 // File: ReservoirPropertyCommon_impl.hpp
4 //
5 // Created: Mon Oct 26 08:29:09 2009
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // B�rd Skaflestad <bard.skaflestad@sintef.no>
9 //
10 // $Date$
11 //
12 // $Revision$
13 //
14 //===========================================================================
15 
16 /*
17  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18  Copyright 2009, 2010 Statoil ASA.
19 
20  This file is part of The Open Reservoir Simulator Project (OpenRS).
21 
22  OpenRS is free software: you can redistribute it and/or modify
23  it under the terms of the GNU General Public License as published by
24  the Free Software Foundation, either version 3 of the License, or
25  (at your option) any later version.
26 
27  OpenRS is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  GNU General Public License for more details.
31 
32  You should have received a copy of the GNU General Public License
33  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34 */
35 
36 #ifndef OPENRS_RESERVOIRPROPERTYCOMMON_IMPL_HEADER
37 #define OPENRS_RESERVOIRPROPERTYCOMMON_IMPL_HEADER
38 
39 #include <opm/output/eclipse/EclipseGridInspector.hpp>
40 #include <opm/parser/eclipse/Deck/Deck.hpp>
41 #include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
42 
43 #include <fstream>
44 #include <array>
45 
46 namespace Opm
47 {
48 
49  namespace {
50 
70  PermeabilityKind classifyPermeability(const Opm::Deck& deck)
71  {
72  const bool xx = deck.hasKeyword("PERMX" );
73  const bool xy = deck.hasKeyword("PERMXY");
74  const bool xz = deck.hasKeyword("PERMXZ");
75 
76  const bool yx = deck.hasKeyword("PERMYX");
77  const bool yy = deck.hasKeyword("PERMY" );
78  const bool yz = deck.hasKeyword("PERMYZ");
79 
80  const bool zx = deck.hasKeyword("PERMZX");
81  const bool zy = deck.hasKeyword("PERMZY");
82  const bool zz = deck.hasKeyword("PERMZ" );
83 
84  int num_cross_comp = xy + xz + yx + yz + zx + zy;
85  int num_comp = xx + yy + zz + num_cross_comp;
86  PermeabilityKind retval = None;
87  if (num_cross_comp > 0) {
88  retval = TensorPerm;
89  } else {
90  if (num_comp == 1) {
91  retval = ScalarPerm;
92  } else if (num_comp >= 2) {
93  retval = DiagonalPerm;
94  }
95  }
96 
97  bool ok = true;
98  if (num_comp > 0) {
99  // At least one tensor component specified on input.
100  // Verify that any remaining components are OK from a
101  // structural point of view. In particular, there
102  // must not be any cross-components (e.g., k_{xy})
103  // unless the corresponding diagonal component (e.g.,
104  // k_{xx}) is present as well...
105  //
106  ok = xx || !(xy || xz || yx || zx) ;
107  ok = ok && (yy || !(yx || yz || xy || zy));
108  ok = ok && (zz || !(zx || zy || xz || yz));
109  }
110  if (!ok) {
111  retval = Invalid;
112  }
113 
114  return retval;
115  }
116 
117 
138  void setScalarPermIfNeeded(std::array<int,9>& kmap,
139  int i, int j, int k)
140  {
141  if (kmap[j] == 0) { kmap[j] = kmap[i]; }
142  if (kmap[k] == 0) { kmap[k] = kmap[i]; }
143  }
144 
145 
178  inline
179  PermeabilityKind fillTensor(const Opm::Deck& deck,
180  std::vector<const std::vector<double>*>& tensor,
181  std::array<int,9>& kmap)
182  {
183  PermeabilityKind kind = classifyPermeability(deck);
184  if (kind == Invalid) {
185  OPM_THROW(std::runtime_error, "Invalid set of permeability fields given.");
186  }
187  assert (tensor.size() == 1);
188  for (int i = 0; i < 9; ++i) { kmap[i] = 0; }
189 
190  enum { xx, xy, xz, // 0, 1, 2
191  yx, yy, yz, // 3, 4, 5
192  zx, zy, zz }; // 6, 7, 8
193 
194  // -----------------------------------------------------------
195  // 1st row: [kxx, kxy, kxz]
196  if (deck.hasKeyword("PERMX" )) {
197  kmap[xx] = tensor.size();
198  tensor.push_back(&deck.getKeyword("PERMX").getSIDoubleData());
199 
200  setScalarPermIfNeeded(kmap, xx, yy, zz);
201  }
202  if (deck.hasKeyword("PERMXY")) {
203  kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry.
204  tensor.push_back(&deck.getKeyword("PERMXY").getSIDoubleData());
205  }
206  if (deck.hasKeyword("PERMXZ")) {
207  kmap[xz] = kmap[zx] = tensor.size(); // Enforce symmetry.
208  tensor.push_back(&deck.getKeyword("PERMXZ").getSIDoubleData());
209  }
210 
211  // -----------------------------------------------------------
212  // 2nd row: [kyx, kyy, kyz]
213  if (deck.hasKeyword("PERMYX")) {
214  kmap[yx] = kmap[xy] = tensor.size(); // Enforce symmetry.
215  tensor.push_back(&deck.getKeyword("PERMYX").getSIDoubleData());
216  }
217  if (deck.hasKeyword("PERMY" )) {
218  kmap[yy] = tensor.size();
219  tensor.push_back(&deck.getKeyword("PERMY").getSIDoubleData());
220 
221  setScalarPermIfNeeded(kmap, yy, zz, xx);
222  }
223  if (deck.hasKeyword("PERMYZ")) {
224  kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry.
225  tensor.push_back(&deck.getKeyword("PERMYZ").getSIDoubleData());
226  }
227 
228  // -----------------------------------------------------------
229  // 3rd row: [kzx, kzy, kzz]
230  if (deck.hasKeyword("PERMZX")) {
231  kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry.
232  tensor.push_back(&deck.getKeyword("PERMZX").getSIDoubleData());
233  }
234  if (deck.hasKeyword("PERMZY")) {
235  kmap[zy] = kmap[yz] = tensor.size(); // Enforce symmetry.
236  tensor.push_back(&deck.getKeyword("PERMZY").getSIDoubleData());
237  }
238  if (deck.hasKeyword("PERMZ" )) {
239  kmap[zz] = tensor.size();
240  tensor.push_back(&deck.getKeyword("PERMZ").getSIDoubleData());
241 
242  setScalarPermIfNeeded(kmap, zz, xx, yy);
243  }
244  return kind;
245  }
246 
247  } // anonymous namespace
248 
249 
250 
251 
252  // ----- Methods of ReservoirPropertyCommon -----
253 
254 
255 
256 
257  template <int dim, class RPImpl, class RockType>
259 #if 1
260  : density1_ (1013.9*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
261  density2_ ( 834.7*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
262  viscosity1_( 1.0*Opm::prefix::centi*Opm::unit::Poise),
263  viscosity2_( 3.0*Opm::prefix::centi*Opm::unit::Poise),
264 #else
265  : density1_ (1000.0*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
266  density2_ (1000.0*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
267  viscosity1_( 1000.0*Opm::prefix::centi*Opm::unit::Poise),
268  viscosity2_( 1000.0*Opm::prefix::centi*Opm::unit::Poise),
269 #endif
270  permeability_kind_(Invalid)
271  {
272  }
273 
274 
275  template <int dim, class RPImpl, class RockType>
277  const std::vector<int>& global_cell,
278  const double perm_threshold,
279  const std::string* rock_list_filename,
280  const bool use_jfunction_scaling,
281  const double sigma,
282  const double theta)
283  {
284  // This code is mostly copied from ReservoirPropertyCommon::init(...).
285  static_assert(dim == 3, "");
286 
287  permfield_valid_.assign(global_cell.size(),
288  std::vector<unsigned char>::value_type(0));
289 
290  assignPorosity (deck, global_cell);
291  assignNTG (deck, global_cell);
292  assignSWCR (deck, global_cell);
293  assignSOWCR (deck, global_cell);
294  assignPermeability(deck, global_cell, perm_threshold);
295  assignRockTable (deck, global_cell);
296 
297  if (rock_list_filename) {
298  readRocks(*rock_list_filename);
299  }
300 
301  // Added section. This is a hack, because not all rock classes
302  // may care about J-scaling. They still have to implement
303  // setUseJfunctionScaling() and setSigmaAndTheta(), though the
304  // latter may throw if called for a rock where it does not make sense.
305  int num_rocks = rock_.size();
306  for (int i = 0; i < num_rocks; ++i) {
307  rock_[i].setUseJfunctionScaling(use_jfunction_scaling);
308  if (use_jfunction_scaling) {
309  rock_[i].setSigmaAndTheta(sigma, theta);
310  }
311  }
312  // End of added section.
313 
314  asImpl().computeCflFactors();
315  }
316 
317 
318 
319  template <int dim, class RPImpl, class RockType>
321  const double uniform_poro,
322  const double uniform_perm)
323  {
324  permfield_valid_.assign(num_cells, std::vector<unsigned char>::value_type(1));
325  porosity_.assign(num_cells, uniform_poro);
326  permeability_.assign(dim*dim*num_cells, 0.0);
327  for (int i = 0; i < num_cells; ++i) {
328  SharedPermTensor K = permeabilityModifiable(i);
329  for (int dd = 0; dd < dim; ++dd) {
330  K(dd, dd) = uniform_perm;
331  }
332  }
333  cell_to_rock_.assign(num_cells, 0);
334  asImpl().computeCflFactors();
335  }
336 
337  template <int dim, class RPImpl, class RockType>
339  {
340  viscosity1_ = v1;
341  viscosity2_ = v2;
342  }
343 
344  template <int dim, class RPImpl, class RockType>
346  {
347  density1_ = d1;
348  density2_ = d2;
349  }
350 
351  template <int dim, class RPImpl, class RockType>
353  {
354  return viscosity1_;
355  }
356 
357 
358  template <int dim, class RPImpl, class RockType>
360  {
361  return viscosity2_;
362  }
363 
364 
365  template <int dim, class RPImpl, class RockType>
367  {
368  return density1_;
369  }
370 
371 
372  template <int dim, class RPImpl, class RockType>
374  {
375  return density2_;
376  }
377 
378 
379  template <int dim, class RPImpl, class RockType>
381  {
382  return porosity_[cell_index];
383  }
384 
385  template <int dim, class RPImpl, class RockType>
387  {
388  return ntg_[cell_index];
389  }
390 
391  template <int dim, class RPImpl, class RockType>
393  {
394  return swcr_[cell_index];
395  }
396 
397  template <int dim, class RPImpl, class RockType>
399  {
400  return sowcr_[cell_index];
401  }
402 
403 
404  template <int dim, class RPImpl, class RockType>
407  {
408  assert (permfield_valid_[cell_index]);
409 
410  const PermTensor K(dim, dim, &permeability_[dim*dim*cell_index]);
411  return K;
412  }
413 
414 
415  template <int dim, class RPImpl, class RockType>
418  {
419  // Typically only used for assigning synthetic perm values.
420  SharedPermTensor K(dim, dim, &permeability_[dim*dim*cell_index]);
421 
422  // Trust caller!
423  permfield_valid_[cell_index] = std::vector<unsigned char>::value_type(1);
424 
425  return K;
426  }
427 
428 
429  template <int dim, class RPImpl, class RockType>
430  template<class Vector>
431  void ReservoirPropertyCommon<dim, RPImpl, RockType>::phaseDensities(int /*cell_index*/, Vector& density) const
432  {
433  assert (density.size() >= NumberOfPhases);
434  density[0] = densityFirstPhase();
435  density[1] = densitySecondPhase();
436  }
437 
438 
439  template <int dim, class RPImpl, class RockType>
441  {
442  return density1_ - density2_;
443  }
444 
445 
446  template <int dim, class RPImpl, class RockType>
448  {
449  return cfl_factor_;
450  }
451 
452 
453  template <int dim, class RPImpl, class RockType>
455  {
456  return cfl_factor_gravity_;
457  }
458 
459 
460  template <int dim, class RPImpl, class RockType>
462  {
463  return cfl_factor_capillary_;
464  }
465 
466  template <int dim, class RPImpl, class RockType>
467  double ReservoirPropertyCommon<dim, RPImpl, RockType>::capillaryPressure(int cell_index, double saturation) const
468  {
469  if (rock_.size() > 0) {
470  int r = cell_to_rock_[cell_index];
471  return rock_[r].capPress(permeability(cell_index), porosity(cell_index), saturation);
472  } else {
473  // HACK ALERT!
474  // Use zero capillary pressure if no known rock table exists.
475  return 1e5*(1-saturation);
476  }
477  }
478 
479  template <int dim, class RPImpl, class RockType>
480  double ReservoirPropertyCommon<dim, RPImpl, RockType>::capillaryPressureDeriv(int cell_index, double saturation) const
481  {
482  if (rock_.size() > 0) {
483  int r = cell_to_rock_[cell_index];
484  double dpc = rock_[r].capPressDeriv(permeability(cell_index), porosity(cell_index), saturation);
485  return dpc;
486  } else {
487  // HACK ALERT!
488  // Use zero capillary pressure if no known rock table exists.
489  return -1.0e5;
490  }
491  }
492  template <int dim, class RPImpl, class RockType>
494  {
495  if (rock_.size() > 0) {
496  int r = cell_to_rock_[cell_index];
497  return rock_[r].s_min();
498  } else {
499  // HACK ALERT!
500  // Use zero as minimum saturation if no known rock table exists.
501  return 0;
502  }
503  }
504 
505  template <int dim, class RPImpl, class RockType>
507  {
508  if (rock_.size() > 0) {
509  int r = cell_to_rock_[cell_index];
510  return rock_[r].s_max();
511  } else {
512  // HACK ALERT!
513  // Use 1 as maximum saturation if no known rock table exists.
514  return 1;
515  }
516  }
517 
518  template <int dim, class RPImpl, class RockType>
520  {
521  if (rock_.size() > 0) {
522  int r = cell_to_rock_[cell_index];
523  return rock_[r].satFromCapPress(permeability(cell_index), porosity(cell_index), cap_press);
524  } else {
525  // HACK ALERT!
526  // Use a zero saturation if no known rock table exists.
527  return 0.0;
528  }
529  }
530 
531 
532  template <int dim, class RPImpl, class RockType>
534  {
535  int num_cells = porosity_.size();
536  // Write porosity.
537  {
538  std::string filename = grid_prefix + "-poro.dat";
539  std::ofstream file(filename.c_str());
540  if (!file) {
541  OPM_THROW(std::runtime_error, "Could not open file " << filename);
542  }
543  file << num_cells << '\n';
544  std::copy(porosity_.begin(), porosity_.end(), std::ostream_iterator<double>(file, "\n"));
545  }
546  // Write permeability.
547  {
548  std::string filename = grid_prefix + "-perm.dat";
549  std::ofstream file(filename.c_str());
550  if (!file) {
551  OPM_THROW(std::runtime_error, "Could not open file " << filename);
552  }
553  file << num_cells << '\n';
554  switch (permeability_kind_) {
555  case TensorPerm:
556  std::copy(permeability_.begin(), permeability_.end(), std::ostream_iterator<double>(file, "\n"));
557  break;
558  case DiagonalPerm:
559  for (int c = 0; c < num_cells; ++c) {
560  int index = c*dim*dim;
561  for (int dd = 0; dd < dim; ++dd) {
562  file << permeability_[index + (dim + 1)*dd] << ' ';
563  }
564  file << '\n';
565  }
566  break;
567  case ScalarPerm:
568  case None: // Treated like a scalar permeability.
569  for (int c = 0; c < num_cells; ++c) {
570  file << permeability_[c*dim*dim] << '\n';
571  }
572  break;
573  default:
574  OPM_THROW(std::runtime_error, "Cannot write invalid permeability.");
575  }
576  }
577  }
578 
579 
580  // ------ Private methods ------
581 
582 
583  template <int dim, class RPImpl, class RockType>
585  {
586  return static_cast<RPImpl&>(*this);
587  }
588 
589 
590 
591 
592  template <int dim, class RPImpl, class RockType>
593  void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignPorosity(const Opm::Deck& deck,
594  const std::vector<int>& global_cell)
595  {
596  porosity_.assign(global_cell.size(), 1.0);
597 
598  if (deck.hasKeyword("PORO")) {
599  Opm::EclipseGridInspector insp(deck);
600  std::array<int, 3> dims = insp.gridSize();
601  int num_global_cells = dims[0]*dims[1]*dims[2];
602  const std::vector<double>& poro = deck.getKeyword("PORO").getSIDoubleData();
603  if (int(poro.size()) != num_global_cells) {
604  OPM_THROW(std::runtime_error, "PORO field must have the same size as the "
605  "logical cartesian size of the grid: "
606  << poro.size() << " != " << num_global_cells);
607  }
608  for (int c = 0; c < int(porosity_.size()); ++c) {
609  porosity_[c] = poro[global_cell[c]];
610  }
611  }
612  }
613 
614  template <int dim, class RPImpl, class RockType>
615  void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignNTG(const Opm::Deck& deck,
616  const std::vector<int>& global_cell)
617  {
618  ntg_.assign(global_cell.size(), 1.0);
619 
620  if (deck.hasKeyword("NTG")) {
621  Opm::EclipseGridInspector insp(deck);
622  std::array<int, 3> dims = insp.gridSize();
623  int num_global_cells = dims[0]*dims[1]*dims[2];
624  const std::vector<double>& ntg = deck.getKeyword("NTG").getSIDoubleData();
625  if (int(ntg.size()) != num_global_cells) {
626  OPM_THROW(std::runtime_error, "NTG field must have the same size as the "
627  "logical cartesian size of the grid: "
628  << ntg.size() << " != " << num_global_cells);
629  }
630  for (int c = 0; c < int(ntg_.size()); ++c) {
631  ntg_[c] = ntg[global_cell[c]];
632  }
633  }
634  }
635 
636  template <int dim, class RPImpl, class RockType>
637  void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignSWCR(const Opm::Deck& deck,
638  const std::vector<int>& global_cell)
639  {
640  swcr_.assign(global_cell.size(), 0.0);
641 
642  if (deck.hasKeyword("SWCR")) {
643  Opm::EclipseGridInspector insp(deck);
644  std::array<int, 3> dims = insp.gridSize();
645  int num_global_cells = dims[0]*dims[1]*dims[2];
646  const std::vector<double>& swcr =
647  deck.getKeyword("SWCR").getSIDoubleData();
648  if (int(swcr.size()) != num_global_cells) {
649  OPM_THROW(std::runtime_error, "SWCR field must have the same size as the "
650  "logical cartesian size of the grid: "
651  << swcr.size() << " != " << num_global_cells);
652  }
653  for (int c = 0; c < int(swcr_.size()); ++c) {
654  swcr_[c] = swcr[global_cell[c]];
655  }
656  }
657  }
658 
659  template <int dim, class RPImpl, class RockType>
660  void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignSOWCR(const Opm::Deck& deck,
661  const std::vector<int>& global_cell)
662  {
663  sowcr_.assign(global_cell.size(), 0.0);
664 
665  if (deck.hasKeyword("SOWCR")) {
666  Opm::EclipseGridInspector insp(deck);
667  std::array<int, 3> dims = insp.gridSize();
668  int num_global_cells = dims[0]*dims[1]*dims[2];
669  const std::vector<double>& sowcr =
670  deck.getKeyword("SOWCR").getSIDoubleData();
671  if (int(sowcr.size()) != num_global_cells) {
672  OPM_THROW(std::runtime_error, "SOWCR field must have the same size as the "
673  "logical cartesian size of the grid: "
674  << sowcr.size() << " != " << num_global_cells);
675  }
676  for (int c = 0; c < int(sowcr_.size()); ++c) {
677  sowcr_[c] = sowcr[global_cell[c]];
678  }
679  }
680  }
681 
682  template <int dim, class RPImpl, class RockType>
683  void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignPermeability(const Opm::Deck& deck,
684  const std::vector<int>& global_cell,
685  double perm_threshold)
686  {
687  Opm::EclipseGridInspector insp(deck);
688  std::array<int, 3> dims = insp.gridSize();
689  int num_global_cells = dims[0]*dims[1]*dims[2];
690  assert (num_global_cells > 0);
691 
692  permeability_.assign(dim * dim * global_cell.size(), 0.0);
693 
694  std::vector<const std::vector<double>*> tensor;
695  tensor.reserve(10);
696 
697  const std::vector<double> zero(num_global_cells, 0.0);
698  tensor.push_back(&zero);
699 
700  static_assert(dim == 3, "");
701  std::array<int,9> kmap;
702  permeability_kind_ = fillTensor(deck, tensor, kmap);
703  for (int i = 1; i < int(tensor.size()); ++i) {
704  if (int(tensor[i]->size()) != num_global_cells) {
705  OPM_THROW(std::runtime_error, "All permeability fields must have the same size as the "
706  "logical cartesian size of the grid: "
707  << (tensor[i]->size()) << " != " << num_global_cells);
708  }
709  }
710 
711  // Assign permeability values only if such values are
712  // given in the input deck represented by 'deck'. In
713  // other words: Don't set any (arbitrary) default values.
714  // It is infinitely better to experience a reproducible
715  // crash than subtle errors resulting from a (poorly
716  // chosen) default value...
717  //
718  if (tensor.size() > 1) {
719  const int nc = global_cell.size();
720  int off = 0;
721 
722  for (int c = 0; c < nc; ++c, off += dim*dim) {
723  SharedPermTensor K(dim, dim, &permeability_[off]);
724  int kix = 0;
725  const int glob = global_cell[c];
726 
727  for (int i = 0; i < dim; ++i) {
728  for (int j = 0; j < dim; ++j, ++kix) {
729  K(i,j) = (*tensor[kmap[kix]])[glob];
730  }
731  K(i,i) = std::max(K(i,i), perm_threshold);
732  }
733 
734  permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
735  }
736  }
737  }
738 
739 
740 
741 
742  template <int dim, class RPImpl, class RockType>
743  void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignRockTable(const Opm::Deck& deck,
744  const std::vector<int>& global_cell)
745  {
746  const int nc = global_cell.size();
747 
748  cell_to_rock_.assign(nc, 0);
749 
750  if (deck.hasKeyword("SATNUM")) {
751  Opm::EclipseGridInspector insp(deck);
752  std::array<int, 3> dims = insp.gridSize();
753  int num_global_cells = dims[0]*dims[1]*dims[2];
754  const std::vector<int>& satnum = deck.getKeyword("SATNUM").getIntData();
755  if (int(satnum.size()) != num_global_cells) {
756  OPM_THROW(std::runtime_error, "SATNUM field must have the same size as the "
757  "logical cartesian size of the grid: "
758  << satnum.size() << " != " << num_global_cells);
759  }
760  for (int c = 0; c < nc; ++c) {
761  // Note: SATNUM is FORTRANish, ranging from 1 to n, therefore we subtract one.
762  cell_to_rock_[c] = satnum[global_cell[c]] - 1;
763  }
764  }
765  else if (deck.hasKeyword("ROCKTYPE")) {
766  Opm::EclipseGridInspector insp(deck);
767  std::array<int, 3> dims = insp.gridSize();
768  int num_global_cells = dims[0]*dims[1]*dims[2];
769  const std::vector<int>& satnum = deck.getKeyword("ROCKTYPE").getIntData();
770  if (int(satnum.size()) != num_global_cells) {
771  OPM_THROW(std::runtime_error, "ROCKTYPE field must have the same size as the "
772  "logical cartesian size of the grid: "
773  << satnum.size() << " != " << num_global_cells);
774  }
775  for (int c = 0; c < nc; ++c) {
776  // Note: ROCKTYPE is FORTRANish, ranging from 1 to n, therefore we subtract one.
777  cell_to_rock_[c] = satnum[global_cell[c]] - 1;
778  }
779  }
780  }
781 
782 
783 
784 
785  template <int dim, class RPImpl, class RockType>
786  void ReservoirPropertyCommon<dim, RPImpl, RockType>::readRocks(const std::string& rock_list_file)
787  {
788  std::ifstream rl(rock_list_file.c_str());
789  if (!rl) {
790  OPM_THROW(std::runtime_error, "Could not open file " << rock_list_file);
791  }
792  int num_rocks = -1;
793  rl >> num_rocks;
794  assert(num_rocks >= 1);
795  rock_.resize(num_rocks);
796  std::string dir(rock_list_file.begin(), rock_list_file.begin() + rock_list_file.find_last_of('/') + 1);
797  for (int i = 0; i < num_rocks; ++i) {
798  std::string spec;
799  while (spec.empty()) {
800  std::getline(rl, spec);
801  }
802  rock_[i].read(dir, spec);
803  }
804  }
805 
806 
807 
808 
809 } // namespace Opm
810 
811 
812 #endif // OPENRS_RESERVOIRPROPERTYCOMMON_IMPL_HEADER
double viscositySecondPhase() const
Viscosity of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:359
ImmutableCMatrix PermTensor
Tensor type for read-only access to permeability.
Definition: ReservoirPropertyCommon.hpp:62
void init(const Opm::Deck &, const std::vector< int > &global_cell, const double perm_threshold=0.0, const std::string *rock_list_filename=0, const bool use_jfunction_scaling=true, const double sigma=1.0, const double theta=0.0)
Initialize from a grdecl file.
Definition: ReservoirPropertyCommon_impl.hpp:276
PermeabilityKind
Enum for the kind of permeability field originally retrieved.
Definition: ReservoirPropertyCommon.hpp:50
double capillaryPressureDeriv(int c, double s) const
Derivative of Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:480
double s_min(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:493
double densitySecondPhase() const
Density of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:373
ReservoirPropertyCommon()
Default constructor.
Definition: ReservoirPropertyCommon_impl.hpp:258
SharedPermTensor permeabilityModifiable(int cell_index)
Read- and write-access to permeability.
Definition: ReservoirPropertyCommon_impl.hpp:417
double swcr(int cell_index) const
Read-access to swcr.
Definition: ReservoirPropertyCommon_impl.hpp:392
double densityDifference() const
Difference of densities.
Definition: ReservoirPropertyCommon_impl.hpp:440
double saturationFromCapillaryPressure(int cell_index, double cap_press) const
Inverse of the capillary pressure function.
Definition: ReservoirPropertyCommon_impl.hpp:519
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
double cflFactor() const
A factor useful in cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:447
void setViscosities(double v1, double v2)
Set viscosities of both faces.
Definition: ReservoirPropertyCommon_impl.hpp:338
A property class for incompressible two-phase flow.
Definition: ReservoirPropertyCommon.hpp:58
void phaseDensities(int, Vector &density) const
Densities for both phases.
Definition: ReservoirPropertyCommon_impl.hpp:431
double cflFactorCapillary() const
A factor useful in gravity cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:461
double s_max(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:506
double cflFactorGravity() const
A factor useful in gravity cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:454
double sowcr(int cell_index) const
Read-access to sowcr.
Definition: ReservoirPropertyCommon_impl.hpp:398
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
double porosity(int cell_index) const
Read-access to porosity.
Definition: ReservoirPropertyCommon_impl.hpp:380
double ntg(int cell_index) const
Read-access to ntg.
Definition: ReservoirPropertyCommon_impl.hpp:386
void writeSintefLegacyFormat(const std::string &grid_prefix) const
Write permeability and porosity in the Sintef legacy format.
Definition: ReservoirPropertyCommon_impl.hpp:533
PermTensor permeability(int cell_index) const
Read-access to permeability.
Definition: ReservoirPropertyCommon_impl.hpp:406
void setDensities(double d1, double d2)
Set densitities of both faces.
Definition: ReservoirPropertyCommon_impl.hpp:345
double viscosityFirstPhase() const
Viscosity of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:352
double capillaryPressure(int cell_index, double saturation) const
Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:467
SharedCMatrix SharedPermTensor
Tensor type for read and write access to permeability.
Definition: ReservoirPropertyCommon.hpp:66
double densityFirstPhase() const
Density of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:366