36#ifndef OPENRS_RESERVOIRPROPERTYCOMMON_IMPL_HEADER
37#define OPENRS_RESERVOIRPROPERTYCOMMON_IMPL_HEADER
39#include <opm/core/io/eclipse/EclipseGridInspector.hpp>
40#include <opm/parser/eclipse/Deck/Deck.hpp>
41#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
72 const bool xx = deck->hasKeyword(
"PERMX" );
73 const bool xy = deck->hasKeyword(
"PERMXY");
74 const bool xz = deck->hasKeyword(
"PERMXZ");
76 const bool yx = deck->hasKeyword(
"PERMYX");
77 const bool yy = deck->hasKeyword(
"PERMY" );
78 const bool yz = deck->hasKeyword(
"PERMYZ");
80 const bool zx = deck->hasKeyword(
"PERMZX");
81 const bool zy = deck->hasKeyword(
"PERMZY");
82 const bool zz = deck->hasKeyword(
"PERMZ" );
84 int num_cross_comp = xy + xz + yx + yz + zx + zy;
85 int num_comp = xx + yy + zz + num_cross_comp;
87 if (num_cross_comp > 0) {
92 }
else if (num_comp >= 2) {
106 ok = xx || !(xy || xz || yx || zx) ;
107 ok = ok && (yy || !(yx || yz || xy || zy));
108 ok = ok && (zz || !(zx || zy || xz || yz));
138 void setScalarPermIfNeeded(std::array<int,9>& kmap,
141 if (kmap[j] == 0) { kmap[j] = kmap[i]; }
142 if (kmap[k] == 0) { kmap[k] = kmap[i]; }
179 std::vector<
const std::vector<double>*>& tensor,
180 std::array<int,9>& kmap)
184 OPM_THROW(std::runtime_error,
"Invalid set of permeability fields given.");
186 assert (tensor.size() == 1);
187 for (
int i = 0; i < 9; ++i) { kmap[i] = 0; }
195 if (deck->hasKeyword(
"PERMX" )) {
196 kmap[xx] = tensor.size();
197 tensor.push_back(&deck->getKeyword(
"PERMX").getSIDoubleData());
199 setScalarPermIfNeeded(kmap, xx, yy, zz);
201 if (deck->hasKeyword(
"PERMXY")) {
202 kmap[xy] = kmap[yx] = tensor.size();
203 tensor.push_back(&deck->getKeyword(
"PERMXY").getSIDoubleData());
205 if (deck->hasKeyword(
"PERMXZ")) {
206 kmap[xz] = kmap[zx] = tensor.size();
207 tensor.push_back(&deck->getKeyword(
"PERMXZ").getSIDoubleData());
212 if (deck->hasKeyword(
"PERMYX")) {
213 kmap[yx] = kmap[xy] = tensor.size();
214 tensor.push_back(&deck->getKeyword(
"PERMYX").getSIDoubleData());
216 if (deck->hasKeyword(
"PERMY" )) {
217 kmap[yy] = tensor.size();
218 tensor.push_back(&deck->getKeyword(
"PERMY").getSIDoubleData());
220 setScalarPermIfNeeded(kmap, yy, zz, xx);
222 if (deck->hasKeyword(
"PERMYZ")) {
223 kmap[yz] = kmap[zy] = tensor.size();
224 tensor.push_back(&deck->getKeyword(
"PERMYZ").getSIDoubleData());
229 if (deck->hasKeyword(
"PERMZX")) {
230 kmap[zx] = kmap[xz] = tensor.size();
231 tensor.push_back(&deck->getKeyword(
"PERMZX").getSIDoubleData());
233 if (deck->hasKeyword(
"PERMZY")) {
234 kmap[zy] = kmap[yz] = tensor.size();
235 tensor.push_back(&deck->getKeyword(
"PERMZY").getSIDoubleData());
237 if (deck->hasKeyword(
"PERMZ" )) {
238 kmap[zz] = tensor.size();
239 tensor.push_back(&deck->getKeyword(
"PERMZ").getSIDoubleData());
241 setScalarPermIfNeeded(kmap, zz, xx, yy);
256 template <
int dim,
class RPImpl,
class RockType>
259 : density1_ (1013.9*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
260 density2_ ( 834.7*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
261 viscosity1_( 1.0*Opm::prefix::centi*Opm::unit::Poise),
262 viscosity2_( 3.0*Opm::prefix::centi*Opm::unit::Poise),
264 : density1_ (1000.0*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
265 density2_ (1000.0*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
266 viscosity1_( 1000.0*Opm::prefix::centi*Opm::unit::Poise),
267 viscosity2_( 1000.0*Opm::prefix::centi*Opm::unit::Poise),
274 template <
int dim,
class RPImpl,
class RockType>
276 const std::vector<int>& global_cell,
277 const double perm_threshold,
278 const std::string* rock_list_filename,
279 const bool use_jfunction_scaling,
284 static_assert(dim == 3,
"");
286 permfield_valid_.assign(global_cell.size(),
287 std::vector<unsigned char>::value_type(0));
289 assignPorosity (deck, global_cell);
290 assignNTG (deck, global_cell);
291 assignSWCR (deck, global_cell);
292 assignSOWCR (deck, global_cell);
293 assignPermeability(deck, global_cell, perm_threshold);
294 assignRockTable (deck, global_cell);
296 if (rock_list_filename) {
297 readRocks(*rock_list_filename);
304 int num_rocks = rock_.size();
305 for (
int i = 0; i < num_rocks; ++i) {
306 rock_[i].setUseJfunctionScaling(use_jfunction_scaling);
307 if (use_jfunction_scaling) {
308 rock_[i].setSigmaAndTheta(sigma, theta);
313 asImpl().computeCflFactors();
318 template <
int dim,
class RPImpl,
class RockType>
320 const double uniform_poro,
321 const double uniform_perm)
323 permfield_valid_.assign(num_cells, std::vector<unsigned char>::value_type(1));
324 porosity_.assign(num_cells, uniform_poro);
325 permeability_.assign(dim*dim*num_cells, 0.0);
326 for (
int i = 0; i < num_cells; ++i) {
328 for (
int dd = 0; dd < dim; ++dd) {
329 K(dd, dd) = uniform_perm;
332 cell_to_rock_.assign(num_cells, 0);
333 asImpl().computeCflFactors();
336 template <
int dim,
class RPImpl,
class RockType>
343 template <
int dim,
class RPImpl,
class RockType>
350 template <
int dim,
class RPImpl,
class RockType>
357 template <
int dim,
class RPImpl,
class RockType>
364 template <
int dim,
class RPImpl,
class RockType>
371 template <
int dim,
class RPImpl,
class RockType>
378 template <
int dim,
class RPImpl,
class RockType>
381 return porosity_[cell_index];
384 template <
int dim,
class RPImpl,
class RockType>
387 return ntg_[cell_index];
390 template <
int dim,
class RPImpl,
class RockType>
393 return swcr_[cell_index];
396 template <
int dim,
class RPImpl,
class RockType>
399 return sowcr_[cell_index];
403 template <
int dim,
class RPImpl,
class RockType>
407 assert (permfield_valid_[cell_index]);
409 const PermTensor K(dim, dim, &permeability_[dim*dim*cell_index]);
414 template <
int dim,
class RPImpl,
class RockType>
422 permfield_valid_[cell_index] = std::vector<unsigned char>::value_type(1);
428 template <
int dim,
class RPImpl,
class RockType>
429 template<
class Vector>
432 assert (density.size() >= NumberOfPhases);
433 density[0] = densityFirstPhase();
434 density[1] = densitySecondPhase();
438 template <
int dim,
class RPImpl,
class RockType>
441 return density1_ - density2_;
445 template <
int dim,
class RPImpl,
class RockType>
452 template <
int dim,
class RPImpl,
class RockType>
455 return cfl_factor_gravity_;
459 template <
int dim,
class RPImpl,
class RockType>
462 return cfl_factor_capillary_;
465 template <
int dim,
class RPImpl,
class RockType>
468 if (rock_.size() > 0) {
469 int r = cell_to_rock_[cell_index];
470 return rock_[r].capPress(permeability(cell_index), porosity(cell_index), saturation);
474 return 1e5*(1-saturation);
478 template <
int dim,
class RPImpl,
class RockType>
481 if (rock_.size() > 0) {
482 int r = cell_to_rock_[cell_index];
483 double dpc = rock_[r].capPressDeriv(permeability(cell_index), porosity(cell_index), saturation);
491 template <
int dim,
class RPImpl,
class RockType>
494 if (rock_.size() > 0) {
495 int r = cell_to_rock_[cell_index];
496 return rock_[r].s_min();
504 template <
int dim,
class RPImpl,
class RockType>
507 if (rock_.size() > 0) {
508 int r = cell_to_rock_[cell_index];
509 return rock_[r].s_max();
517 template <
int dim,
class RPImpl,
class RockType>
520 if (rock_.size() > 0) {
521 int r = cell_to_rock_[cell_index];
522 return rock_[r].satFromCapPress(permeability(cell_index), porosity(cell_index), cap_press);
531 template <
int dim,
class RPImpl,
class RockType>
534 int num_cells = porosity_.size();
537 std::string filename = grid_prefix +
"-poro.dat";
538 std::ofstream file(filename.c_str());
540 OPM_THROW(std::runtime_error,
"Could not open file " << filename);
542 file << num_cells <<
'\n';
543 std::copy(porosity_.begin(), porosity_.end(), std::ostream_iterator<double>(file,
"\n"));
547 std::string filename = grid_prefix +
"-perm.dat";
548 std::ofstream file(filename.c_str());
550 OPM_THROW(std::runtime_error,
"Could not open file " << filename);
552 file << num_cells <<
'\n';
553 switch (permeability_kind_) {
555 std::copy(permeability_.begin(), permeability_.end(), std::ostream_iterator<double>(file,
"\n"));
558 for (
int c = 0; c < num_cells; ++c) {
559 int index = c*dim*dim;
560 for (
int dd = 0; dd < dim; ++dd) {
561 file << permeability_[index + (dim + 1)*dd] <<
' ';
568 for (
int c = 0; c < num_cells; ++c) {
569 file << permeability_[c*dim*dim] <<
'\n';
573 OPM_THROW(std::runtime_error,
"Cannot write invalid permeability.");
582 template <
int dim,
class RPImpl,
class RockType>
585 return static_cast<RPImpl&
>(*this);
591 template <
int dim,
class RPImpl,
class RockType>
593 const std::vector<int>& global_cell)
595 porosity_.assign(global_cell.size(), 1.0);
597 if (deck->hasKeyword(
"PORO")) {
598 Opm::EclipseGridInspector insp(deck);
599 std::array<int, 3> dims = insp.gridSize();
600 int num_global_cells = dims[0]*dims[1]*dims[2];
601 const std::vector<double>& poro = deck->getKeyword(
"PORO").getSIDoubleData();
602 if (
int(poro.size()) != num_global_cells) {
603 OPM_THROW(std::runtime_error,
"PORO field must have the same size as the "
604 "logical cartesian size of the grid: "
605 << poro.size() <<
" != " << num_global_cells);
607 for (
int c = 0; c < int(porosity_.size()); ++c) {
608 porosity_[c] = poro[global_cell[c]];
613 template <
int dim,
class RPImpl,
class RockType>
615 const std::vector<int>& global_cell)
617 ntg_.assign(global_cell.size(), 1.0);
619 if (deck->hasKeyword(
"NTG")) {
620 Opm::EclipseGridInspector insp(deck);
621 std::array<int, 3> dims = insp.gridSize();
622 int num_global_cells = dims[0]*dims[1]*dims[2];
623 const std::vector<double>& ntg = deck->getKeyword(
"NTG").getSIDoubleData();
624 if (
int(ntg.size()) != num_global_cells) {
625 OPM_THROW(std::runtime_error,
"NTG field must have the same size as the "
626 "logical cartesian size of the grid: "
627 << ntg.size() <<
" != " << num_global_cells);
629 for (
int c = 0; c < int(ntg_.size()); ++c) {
630 ntg_[c] = ntg[global_cell[c]];
635 template <
int dim,
class RPImpl,
class RockType>
637 const std::vector<int>& global_cell)
639 swcr_.assign(global_cell.size(), 0.0);
641 if (deck->hasKeyword(
"SWCR")) {
642 Opm::EclipseGridInspector insp(deck);
643 std::array<int, 3> dims = insp.gridSize();
644 int num_global_cells = dims[0]*dims[1]*dims[2];
645 const std::vector<double>& swcr =
646 deck->getKeyword(
"SWCR").getSIDoubleData();
647 if (
int(swcr.size()) != num_global_cells) {
648 OPM_THROW(std::runtime_error,
"SWCR field must have the same size as the "
649 "logical cartesian size of the grid: "
650 << swcr.size() <<
" != " << num_global_cells);
652 for (
int c = 0; c < int(swcr_.size()); ++c) {
653 swcr_[c] = swcr[global_cell[c]];
658 template <
int dim,
class RPImpl,
class RockType>
660 const std::vector<int>& global_cell)
662 sowcr_.assign(global_cell.size(), 0.0);
664 if (deck->hasKeyword(
"SOWCR")) {
665 Opm::EclipseGridInspector insp(deck);
666 std::array<int, 3> dims = insp.gridSize();
667 int num_global_cells = dims[0]*dims[1]*dims[2];
668 const std::vector<double>& sowcr =
669 deck->getKeyword(
"SOWCR").getSIDoubleData();
670 if (
int(sowcr.size()) != num_global_cells) {
671 OPM_THROW(std::runtime_error,
"SOWCR field must have the same size as the "
672 "logical cartesian size of the grid: "
673 << sowcr.size() <<
" != " << num_global_cells);
675 for (
int c = 0; c < int(sowcr_.size()); ++c) {
676 sowcr_[c] = sowcr[global_cell[c]];
681 template <
int dim,
class RPImpl,
class RockType>
683 const std::vector<int>& global_cell,
684 double perm_threshold)
686 Opm::EclipseGridInspector insp(deck);
687 std::array<int, 3> dims = insp.gridSize();
688 int num_global_cells = dims[0]*dims[1]*dims[2];
689 assert (num_global_cells > 0);
691 permeability_.assign(dim * dim * global_cell.size(), 0.0);
693 std::vector<const std::vector<double>*> tensor;
696 const std::vector<double>
zero(num_global_cells, 0.0);
697 tensor.push_back(&
zero);
699 static_assert(dim == 3,
"");
700 std::array<int,9> kmap;
701 permeability_kind_ = fillTensor(deck, tensor, kmap);
702 for (
int i = 1; i < int(tensor.size()); ++i) {
703 if (
int(tensor[i]->size()) != num_global_cells) {
704 OPM_THROW(std::runtime_error,
"All permeability fields must have the same size as the "
705 "logical cartesian size of the grid: "
706 << (tensor[i]->size()) <<
" != " << num_global_cells);
717 if (tensor.size() > 1) {
718 const int nc = global_cell.size();
721 for (
int c = 0; c < nc; ++c, off += dim*dim) {
724 const int glob = global_cell[c];
726 for (
int i = 0; i < dim; ++i) {
727 for (
int j = 0; j < dim; ++j, ++kix) {
728 K(i,j) = (*tensor[kmap[kix]])[glob];
730 K(i,i) = std::max(K(i,i), perm_threshold);
733 permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
741 template <
int dim,
class RPImpl,
class RockType>
743 const std::vector<int>& global_cell)
745 const int nc = global_cell.size();
747 cell_to_rock_.assign(nc, 0);
749 if (deck->hasKeyword(
"SATNUM")) {
750 Opm::EclipseGridInspector insp(deck);
751 std::array<int, 3> dims = insp.gridSize();
752 int num_global_cells = dims[0]*dims[1]*dims[2];
753 const std::vector<int>& satnum = deck->getKeyword(
"SATNUM").getIntData();
754 if (
int(satnum.size()) != num_global_cells) {
755 OPM_THROW(std::runtime_error,
"SATNUM field must have the same size as the "
756 "logical cartesian size of the grid: "
757 << satnum.size() <<
" != " << num_global_cells);
759 for (
int c = 0; c < nc; ++c) {
761 cell_to_rock_[c] = satnum[global_cell[c]] - 1;
764 else if (deck->hasKeyword(
"ROCKTYPE")) {
765 Opm::EclipseGridInspector insp(deck);
766 std::array<int, 3> dims = insp.gridSize();
767 int num_global_cells = dims[0]*dims[1]*dims[2];
768 const std::vector<int>& satnum = deck->getKeyword(
"ROCKTYPE").getIntData();
769 if (
int(satnum.size()) != num_global_cells) {
770 OPM_THROW(std::runtime_error,
"ROCKTYPE field must have the same size as the "
771 "logical cartesian size of the grid: "
772 << satnum.size() <<
" != " << num_global_cells);
774 for (
int c = 0; c < nc; ++c) {
776 cell_to_rock_[c] = satnum[global_cell[c]] - 1;
784 template <
int dim,
class RPImpl,
class RockType>
787 std::ifstream rl(rock_list_file.c_str());
789 OPM_THROW(std::runtime_error,
"Could not open file " << rock_list_file);
793 assert(num_rocks >= 1);
794 rock_.resize(num_rocks);
795 std::string dir(rock_list_file.begin(), rock_list_file.begin() + rock_list_file.find_last_of(
'/') + 1);
796 for (
int i = 0; i < num_rocks; ++i) {
798 while (spec.empty()) {
799 std::getline(rl, spec);
801 rock_[i].read(dir, spec);
double porosity(int cell_index) const
Read-access to porosity.
Definition: ReservoirPropertyCommon_impl.hpp:379
ReservoirPropertyCommon()
Default constructor.
Definition: ReservoirPropertyCommon_impl.hpp:257
double densitySecondPhase() const
Density of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:372
double cflFactorCapillary() const
A factor useful in gravity cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:460
void assignPorosity(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:592
void setViscosities(double v1, double v2)
Set viscosities of both faces.
Definition: ReservoirPropertyCommon_impl.hpp:337
double viscosityFirstPhase() const
Viscosity of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:351
double viscositySecondPhase() const
Viscosity of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:358
double saturationFromCapillaryPressure(int cell_index, double cap_press) const
Inverse of the capillary pressure function.
Definition: ReservoirPropertyCommon_impl.hpp:518
void assignNTG(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:614
double densityDifference() const
Difference of densities.
Definition: ReservoirPropertyCommon_impl.hpp:439
void assignSWCR(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:636
PermTensor permeability(int cell_index) const
Read-access to permeability.
Definition: ReservoirPropertyCommon_impl.hpp:405
double densityFirstPhase() const
Density of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:365
double capillaryPressure(int cell_index, double saturation) const
Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:466
double capillaryPressureDeriv(int c, double s) const
Derivative of Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:479
ImmutableCMatrix PermTensor
Tensor type for read-only access to permeability.
Definition: ReservoirPropertyCommon.hpp:62
void readRocks(const std::string &rock_list_file)
Definition: ReservoirPropertyCommon_impl.hpp:785
void init(Opm::DeckConstPtr 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:275
double swcr(int cell_index) const
Read-access to swcr.
Definition: ReservoirPropertyCommon_impl.hpp:391
double cflFactorGravity() const
A factor useful in gravity cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:453
void setDensities(double d1, double d2)
Set densitities of both faces.
Definition: ReservoirPropertyCommon_impl.hpp:344
void assignPermeability(Opm::DeckConstPtr deck, const std::vector< int > &global_cell, const double perm_threshold)
Definition: ReservoirPropertyCommon_impl.hpp:682
double s_min(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:492
void assignRockTable(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:742
void assignSOWCR(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:659
SharedCMatrix SharedPermTensor
Tensor type for read and write access to permeability.
Definition: ReservoirPropertyCommon.hpp:66
RPImpl & asImpl()
Definition: ReservoirPropertyCommon_impl.hpp:583
double ntg(int cell_index) const
Read-access to ntg.
Definition: ReservoirPropertyCommon_impl.hpp:385
void writeSintefLegacyFormat(const std::string &grid_prefix) const
Write permeability and porosity in the Sintef legacy format.
Definition: ReservoirPropertyCommon_impl.hpp:532
double cflFactor() const
A factor useful in cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:446
void phaseDensities(int, Vector &density) const
Densities for both phases.
Definition: ReservoirPropertyCommon_impl.hpp:430
double s_max(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:505
SharedPermTensor permeabilityModifiable(int cell_index)
Read- and write-access to permeability. Use with caution.
Definition: ReservoirPropertyCommon_impl.hpp:416
double sowcr(int cell_index) const
Read-access to sowcr.
Definition: ReservoirPropertyCommon_impl.hpp:397
Definition: BlackoilFluid.hpp:32
PermeabilityKind
Enum for the kind of permeability field originally retrieved.
Definition: ReservoirPropertyCommon.hpp:50
@ Invalid
Definition: ReservoirPropertyCommon.hpp:50
@ ScalarPerm
Definition: ReservoirPropertyCommon.hpp:50
@ None
Definition: ReservoirPropertyCommon.hpp:50
@ DiagonalPerm
Definition: ReservoirPropertyCommon.hpp:50
@ TensorPerm
Definition: ReservoirPropertyCommon.hpp:50
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603