My Project
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/input/eclipse/Deck/Deck.hpp>
41#include <opm/input/eclipse/Deck/DeckKeyword.hpp>
42
43#include <fstream>
44#include <array>
45
46namespace 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));
110 if (!ok) {
111 retval = Invalid;
112 }
114 return retval;
115 }
116
138 void setScalarPermIfNeeded(std::array<int,9>& kmap,
139 int i, int j, int k)
141 if (kmap[j] == 0) { kmap[j] = kmap[i]; }
142 if (kmap[k] == 0) { kmap[k] = kmap[i]; }
143 }
144
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["PERMX"].back().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["PERMXY"].back().getSIDoubleData());
205 }
206 if (deck.hasKeyword("PERMXZ")) {
207 kmap[xz] = kmap[zx] = tensor.size(); // Enforce symmetry.
208 tensor.push_back(&deck["PERMXZ"].back().getSIDoubleData());
209 }
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["PERMYX"].back().getSIDoubleData());
216 }
217 if (deck.hasKeyword("PERMY" )) {
218 kmap[yy] = tensor.size();
219 tensor.push_back(&deck["PERMY"].back().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["PERMYZ"].back().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["PERMZX"].back().getSIDoubleData());
233 }
234 if (deck.hasKeyword("PERMZY")) {
235 kmap[zy] = kmap[yz] = tensor.size(); // Enforce symmetry.
236 tensor.push_back(&deck["PERMZY"].back().getSIDoubleData());
237 }
238 if (deck.hasKeyword("PERMZ" )) {
239 kmap[zz] = tensor.size();
240 tensor.push_back(&deck["PERMZ"].back().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>
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["PORO"].back().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["NTG"].back().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["SWCR"].back().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["SOWCR"].back().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["SATNUM"].back().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["ROCKTYPE"].back().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
A property class for incompressible two-phase flow.
Definition: ReservoirPropertyCommon.hpp:59
double porosity(int cell_index) const
Read-access to porosity.
Definition: ReservoirPropertyCommon_impl.hpp:380
ReservoirPropertyCommon()
Default constructor.
Definition: ReservoirPropertyCommon_impl.hpp:258
double densitySecondPhase() const
Density of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:373
double cflFactorCapillary() const
A factor useful in gravity cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:461
void setViscosities(double v1, double v2)
Set viscosities of both faces.
Definition: ReservoirPropertyCommon_impl.hpp:338
double viscosityFirstPhase() const
Viscosity of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:352
double viscositySecondPhase() const
Viscosity of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:359
double saturationFromCapillaryPressure(int cell_index, double cap_press) const
Inverse of the capillary pressure function.
Definition: ReservoirPropertyCommon_impl.hpp:519
double densityDifference() const
Difference of densities.
Definition: ReservoirPropertyCommon_impl.hpp:440
PermTensor permeability(int cell_index) const
Read-access to permeability.
Definition: ReservoirPropertyCommon_impl.hpp:406
double densityFirstPhase() const
Density of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:366
double capillaryPressure(int cell_index, double saturation) const
Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:467
double capillaryPressureDeriv(int c, double s) const
Derivative of Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:480
ImmutableCMatrix PermTensor
Tensor type for read-only access to permeability.
Definition: ReservoirPropertyCommon.hpp:62
double swcr(int cell_index) const
Read-access to swcr.
Definition: ReservoirPropertyCommon_impl.hpp:392
double cflFactorGravity() const
A factor useful in gravity cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:454
void setDensities(double d1, double d2)
Set densitities of both faces.
Definition: ReservoirPropertyCommon_impl.hpp:345
double s_min(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:493
SharedCMatrix SharedPermTensor
Tensor type for read and write access to permeability.
Definition: ReservoirPropertyCommon.hpp:66
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
double cflFactor() const
A factor useful in cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:447
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
void phaseDensities(int, Vector &density) const
Densities for both phases.
Definition: ReservoirPropertyCommon_impl.hpp:431
double s_max(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:506
SharedPermTensor permeabilityModifiable(int cell_index)
Read- and write-access to permeability.
Definition: ReservoirPropertyCommon_impl.hpp:417
double sowcr(int cell_index) const
Read-access to sowcr.
Definition: ReservoirPropertyCommon_impl.hpp:398
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
PermeabilityKind
Enum for the kind of permeability field originally retrieved.
Definition: ReservoirPropertyCommon.hpp:50
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:602