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>
71 const bool xx = deck->hasKeyword(
"PERMX" );
72 const bool xy = deck->hasKeyword(
"PERMXY");
73 const bool xz = deck->hasKeyword(
"PERMXZ");
75 const bool yx = deck->hasKeyword(
"PERMYX");
76 const bool yy = deck->hasKeyword(
"PERMY" );
77 const bool yz = deck->hasKeyword(
"PERMYZ");
79 const bool zx = deck->hasKeyword(
"PERMZX");
80 const bool zy = deck->hasKeyword(
"PERMZY");
81 const bool zz = deck->hasKeyword(
"PERMZ" );
83 int num_cross_comp = xy + xz + yx + yz + zx + zy;
84 int num_comp = xx + yy + zz + num_cross_comp;
86 if (num_cross_comp > 0) {
91 }
else if (num_comp >= 2) {
105 ok = xx || !(xy || xz || yx || zx) ;
106 ok = ok && (yy || !(yx || yz || xy || zy));
107 ok = ok && (zz || !(zx || zy || xz || yz));
137 void setScalarPermIfNeeded(std::array<int,9>& kmap,
140 if (kmap[j] == 0) { kmap[j] = kmap[i]; }
141 if (kmap[k] == 0) { kmap[k] = kmap[i]; }
178 std::vector<
const std::vector<double>*>& tensor,
179 std::array<int,9>& kmap)
183 OPM_THROW(std::runtime_error,
"Invalid set of permeability fields given.");
185 assert (tensor.size() == 1);
186 for (
int i = 0; i < 9; ++i) { kmap[i] = 0; }
194 if (deck->hasKeyword(
"PERMX" )) {
195 kmap[xx] = tensor.size();
196 tensor.push_back(&deck->getKeyword(
"PERMX")->getSIDoubleData());
198 setScalarPermIfNeeded(kmap, xx, yy, zz);
200 if (deck->hasKeyword(
"PERMXY")) {
201 kmap[xy] = kmap[yx] = tensor.size();
202 tensor.push_back(&deck->getKeyword(
"PERMXY")->getSIDoubleData());
204 if (deck->hasKeyword(
"PERMXZ")) {
205 kmap[xz] = kmap[zx] = tensor.size();
206 tensor.push_back(&deck->getKeyword(
"PERMXZ")->getSIDoubleData());
211 if (deck->hasKeyword(
"PERMYX")) {
212 kmap[yx] = kmap[xy] = tensor.size();
213 tensor.push_back(&deck->getKeyword(
"PERMYX")->getSIDoubleData());
215 if (deck->hasKeyword(
"PERMY" )) {
216 kmap[yy] = tensor.size();
217 tensor.push_back(&deck->getKeyword(
"PERMY")->getSIDoubleData());
219 setScalarPermIfNeeded(kmap, yy, zz, xx);
221 if (deck->hasKeyword(
"PERMYZ")) {
222 kmap[yz] = kmap[zy] = tensor.size();
223 tensor.push_back(&deck->getKeyword(
"PERMYZ")->getSIDoubleData());
228 if (deck->hasKeyword(
"PERMZX")) {
229 kmap[zx] = kmap[xz] = tensor.size();
230 tensor.push_back(&deck->getKeyword(
"PERMZX")->getSIDoubleData());
232 if (deck->hasKeyword(
"PERMZY")) {
233 kmap[zy] = kmap[yz] = tensor.size();
234 tensor.push_back(&deck->getKeyword(
"PERMZY")->getSIDoubleData());
236 if (deck->hasKeyword(
"PERMZ" )) {
237 kmap[zz] = tensor.size();
238 tensor.push_back(&deck->getKeyword(
"PERMZ")->getSIDoubleData());
240 setScalarPermIfNeeded(kmap, zz, xx, yy);
255 template <
int dim,
class RPImpl,
class RockType>
258 : density1_ (1013.9*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
259 density2_ ( 834.7*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
260 viscosity1_( 1.0*Opm::prefix::centi*Opm::unit::Poise),
261 viscosity2_( 3.0*Opm::prefix::centi*Opm::unit::Poise),
263 : density1_ (1000.0*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
264 density2_ (1000.0*Opm::unit::kilogram/Opm::unit::cubic(Opm::unit::meter)),
265 viscosity1_( 1000.0*Opm::prefix::centi*Opm::unit::Poise),
266 viscosity2_( 1000.0*Opm::prefix::centi*Opm::unit::Poise),
273 template <
int dim,
class RPImpl,
class RockType>
275 const std::vector<int>& global_cell,
276 const double perm_threshold,
277 const std::string* rock_list_filename,
278 const bool use_jfunction_scaling,
283 static_assert(dim == 3,
"");
285 permfield_valid_.assign(global_cell.size(),
286 std::vector<unsigned char>::value_type(0));
288 assignPorosity (deck, global_cell);
289 assignNTG (deck, global_cell);
290 assignSWCR (deck, global_cell);
291 assignSOWCR (deck, global_cell);
292 assignPermeability(deck, global_cell, perm_threshold);
293 assignRockTable (deck, global_cell);
295 if (rock_list_filename) {
296 readRocks(*rock_list_filename);
303 int num_rocks = rock_.size();
304 for (
int i = 0; i < num_rocks; ++i) {
305 rock_[i].setUseJfunctionScaling(use_jfunction_scaling);
306 if (use_jfunction_scaling) {
307 rock_[i].setSigmaAndTheta(sigma, theta);
312 asImpl().computeCflFactors();
317 template <
int dim,
class RPImpl,
class RockType>
319 const double uniform_poro,
320 const double uniform_perm)
322 permfield_valid_.assign(num_cells, std::vector<unsigned char>::value_type(1));
323 porosity_.assign(num_cells, uniform_poro);
324 permeability_.assign(dim*dim*num_cells, 0.0);
325 for (
int i = 0; i < num_cells; ++i) {
327 for (
int dd = 0; dd < dim; ++dd) {
328 K(dd, dd) = uniform_perm;
331 cell_to_rock_.assign(num_cells, 0);
332 asImpl().computeCflFactors();
335 template <
int dim,
class RPImpl,
class RockType>
342 template <
int dim,
class RPImpl,
class RockType>
349 template <
int dim,
class RPImpl,
class RockType>
356 template <
int dim,
class RPImpl,
class RockType>
363 template <
int dim,
class RPImpl,
class RockType>
370 template <
int dim,
class RPImpl,
class RockType>
377 template <
int dim,
class RPImpl,
class RockType>
380 return porosity_[cell_index];
383 template <
int dim,
class RPImpl,
class RockType>
386 return ntg_[cell_index];
389 template <
int dim,
class RPImpl,
class RockType>
392 return swcr_[cell_index];
395 template <
int dim,
class RPImpl,
class RockType>
398 return sowcr_[cell_index];
402 template <
int dim,
class RPImpl,
class RockType>
406 assert (permfield_valid_[cell_index]);
408 const PermTensor K(dim, dim, &permeability_[dim*dim*cell_index]);
413 template <
int dim,
class RPImpl,
class RockType>
421 permfield_valid_[cell_index] = std::vector<unsigned char>::value_type(1);
427 template <
int dim,
class RPImpl,
class RockType>
428 template<
class Vector>
431 assert (density.size() >= NumberOfPhases);
432 density[0] = densityFirstPhase();
433 density[1] = densitySecondPhase();
437 template <
int dim,
class RPImpl,
class RockType>
440 return density1_ - density2_;
444 template <
int dim,
class RPImpl,
class RockType>
451 template <
int dim,
class RPImpl,
class RockType>
454 return cfl_factor_gravity_;
458 template <
int dim,
class RPImpl,
class RockType>
461 return cfl_factor_capillary_;
464 template <
int dim,
class RPImpl,
class RockType>
467 if (rock_.size() > 0) {
468 int r = cell_to_rock_[cell_index];
469 return rock_[r].capPress(permeability(cell_index), porosity(cell_index), saturation);
473 return 1e5*(1-saturation);
477 template <
int dim,
class RPImpl,
class RockType>
480 if (rock_.size() > 0) {
481 int r = cell_to_rock_[cell_index];
482 double dpc = rock_[r].capPressDeriv(permeability(cell_index), porosity(cell_index), saturation);
490 template <
int dim,
class RPImpl,
class RockType>
493 if (rock_.size() > 0) {
494 int r = cell_to_rock_[cell_index];
495 return rock_[r].s_min();
503 template <
int dim,
class RPImpl,
class RockType>
506 if (rock_.size() > 0) {
507 int r = cell_to_rock_[cell_index];
508 return rock_[r].s_max();
516 template <
int dim,
class RPImpl,
class RockType>
519 if (rock_.size() > 0) {
520 int r = cell_to_rock_[cell_index];
521 return rock_[r].satFromCapPress(permeability(cell_index), porosity(cell_index), cap_press);
530 template <
int dim,
class RPImpl,
class RockType>
533 int num_cells = porosity_.size();
536 std::string filename = grid_prefix +
"-poro.dat";
537 std::ofstream file(filename.c_str());
539 OPM_THROW(std::runtime_error,
"Could not open file " << filename);
541 file << num_cells <<
'\n';
542 std::copy(porosity_.begin(), porosity_.end(), std::ostream_iterator<double>(file,
"\n"));
546 std::string filename = grid_prefix +
"-perm.dat";
547 std::ofstream file(filename.c_str());
549 OPM_THROW(std::runtime_error,
"Could not open file " << filename);
551 file << num_cells <<
'\n';
552 switch (permeability_kind_) {
554 std::copy(permeability_.begin(), permeability_.end(), std::ostream_iterator<double>(file,
"\n"));
557 for (
int c = 0; c < num_cells; ++c) {
558 int index = c*dim*dim;
559 for (
int dd = 0; dd < dim; ++dd) {
560 file << permeability_[index + (dim + 1)*dd] <<
' ';
567 for (
int c = 0; c < num_cells; ++c) {
568 file << permeability_[c*dim*dim] <<
'\n';
572 OPM_THROW(std::runtime_error,
"Cannot write invalid permeability.");
581 template <
int dim,
class RPImpl,
class RockType>
584 return static_cast<RPImpl&
>(*this);
590 template <
int dim,
class RPImpl,
class RockType>
592 const std::vector<int>& global_cell)
594 porosity_.assign(global_cell.size(), 1.0);
596 if (deck->hasKeyword(
"PORO")) {
597 Opm::EclipseGridInspector insp(deck);
598 std::array<int, 3> dims = insp.gridSize();
599 int num_global_cells = dims[0]*dims[1]*dims[2];
600 const std::vector<double>& poro = deck->getKeyword(
"PORO")->getSIDoubleData();
601 if (
int(poro.size()) != num_global_cells) {
602 OPM_THROW(std::runtime_error,
"PORO field must have the same size as the "
603 "logical cartesian size of the grid: "
604 << poro.size() <<
" != " << num_global_cells);
606 for (
int c = 0; c < int(porosity_.size()); ++c) {
607 porosity_[c] = poro[global_cell[c]];
612 template <
int dim,
class RPImpl,
class RockType>
614 const std::vector<int>& global_cell)
616 ntg_.assign(global_cell.size(), 1.0);
618 if (deck->hasKeyword(
"NTG")) {
619 Opm::EclipseGridInspector insp(deck);
620 std::array<int, 3> dims = insp.gridSize();
621 int num_global_cells = dims[0]*dims[1]*dims[2];
622 const std::vector<double>& ntg = deck->getKeyword(
"NTG")->getSIDoubleData();
623 if (
int(ntg.size()) != num_global_cells) {
624 OPM_THROW(std::runtime_error,
"NTG field must have the same size as the "
625 "logical cartesian size of the grid: "
626 << ntg.size() <<
" != " << num_global_cells);
628 for (
int c = 0; c < int(ntg_.size()); ++c) {
629 ntg_[c] = ntg[global_cell[c]];
634 template <
int dim,
class RPImpl,
class RockType>
636 const std::vector<int>& global_cell)
638 swcr_.assign(global_cell.size(), 0.0);
640 if (deck->hasKeyword(
"SWCR")) {
641 Opm::EclipseGridInspector insp(deck);
642 std::array<int, 3> dims = insp.gridSize();
643 int num_global_cells = dims[0]*dims[1]*dims[2];
644 const std::vector<double>& swcr =
645 deck->getKeyword(
"SWCR")->getSIDoubleData();
646 if (
int(swcr.size()) != num_global_cells) {
647 OPM_THROW(std::runtime_error,
"SWCR field must have the same size as the "
648 "logical cartesian size of the grid: "
649 << swcr.size() <<
" != " << num_global_cells);
651 for (
int c = 0; c < int(swcr_.size()); ++c) {
652 swcr_[c] = swcr[global_cell[c]];
657 template <
int dim,
class RPImpl,
class RockType>
659 const std::vector<int>& global_cell)
661 sowcr_.assign(global_cell.size(), 0.0);
663 if (deck->hasKeyword(
"SOWCR")) {
664 Opm::EclipseGridInspector insp(deck);
665 std::array<int, 3> dims = insp.gridSize();
666 int num_global_cells = dims[0]*dims[1]*dims[2];
667 const std::vector<double>& sowcr =
668 deck->getKeyword(
"SOWCR")->getSIDoubleData();
669 if (
int(sowcr.size()) != num_global_cells) {
670 OPM_THROW(std::runtime_error,
"SOWCR field must have the same size as the "
671 "logical cartesian size of the grid: "
672 << sowcr.size() <<
" != " << num_global_cells);
674 for (
int c = 0; c < int(sowcr_.size()); ++c) {
675 sowcr_[c] = sowcr[global_cell[c]];
680 template <
int dim,
class RPImpl,
class RockType>
682 const std::vector<int>& global_cell,
683 double perm_threshold)
685 Opm::EclipseGridInspector insp(deck);
686 std::array<int, 3> dims = insp.gridSize();
687 int num_global_cells = dims[0]*dims[1]*dims[2];
688 assert (num_global_cells > 0);
690 permeability_.assign(dim * dim * global_cell.size(), 0.0);
692 std::vector<const std::vector<double>*> tensor;
695 const std::vector<double>
zero(num_global_cells, 0.0);
696 tensor.push_back(&zero);
698 static_assert(dim == 3,
"");
699 std::array<int,9> kmap;
700 permeability_kind_ = fillTensor(deck, tensor, kmap);
701 for (
int i = 1; i < int(tensor.size()); ++i) {
702 if (
int(tensor[i]->size()) != num_global_cells) {
703 OPM_THROW(std::runtime_error,
"All permeability fields must have the same size as the "
704 "logical cartesian size of the grid: "
705 << (tensor[i]->size()) <<
" != " << num_global_cells);
716 if (tensor.size() > 1) {
717 const int nc = global_cell.size();
720 for (
int c = 0; c < nc; ++c, off += dim*dim) {
723 const int glob = global_cell[c];
725 for (
int i = 0; i < dim; ++i) {
726 for (
int j = 0; j < dim; ++j, ++kix) {
727 K(i,j) = (*tensor[kmap[kix]])[glob];
729 K(i,i) = std::max(K(i,i), perm_threshold);
732 permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
740 template <
int dim,
class RPImpl,
class RockType>
742 const std::vector<int>& global_cell)
744 const int nc = global_cell.size();
746 cell_to_rock_.assign(nc, 0);
748 if (deck->hasKeyword(
"SATNUM")) {
749 Opm::EclipseGridInspector insp(deck);
750 std::array<int, 3> dims = insp.gridSize();
751 int num_global_cells = dims[0]*dims[1]*dims[2];
752 const std::vector<int>& satnum = deck->getKeyword(
"SATNUM")->getIntData();
753 if (
int(satnum.size()) != num_global_cells) {
754 OPM_THROW(std::runtime_error,
"SATNUM field must have the same size as the "
755 "logical cartesian size of the grid: "
756 << satnum.size() <<
" != " << num_global_cells);
758 for (
int c = 0; c < nc; ++c) {
760 cell_to_rock_[c] = satnum[global_cell[c]] - 1;
763 else if (deck->hasKeyword(
"ROCKTYPE")) {
764 Opm::EclipseGridInspector insp(deck);
765 std::array<int, 3> dims = insp.gridSize();
766 int num_global_cells = dims[0]*dims[1]*dims[2];
767 const std::vector<int>& satnum = deck->getKeyword(
"ROCKTYPE")->getIntData();
768 if (
int(satnum.size()) != num_global_cells) {
769 OPM_THROW(std::runtime_error,
"ROCKTYPE field must have the same size as the "
770 "logical cartesian size of the grid: "
771 << satnum.size() <<
" != " << num_global_cells);
773 for (
int c = 0; c < nc; ++c) {
775 cell_to_rock_[c] = satnum[global_cell[c]] - 1;
783 template <
int dim,
class RPImpl,
class RockType>
786 std::ifstream rl(rock_list_file.c_str());
788 OPM_THROW(std::runtime_error,
"Could not open file " << rock_list_file);
792 assert(num_rocks >= 1);
793 rock_.resize(num_rocks);
794 std::string dir(rock_list_file.begin(), rock_list_file.begin() + rock_list_file.find_last_of(
'/') + 1);
795 for (
int i = 0; i < num_rocks; ++i) {
797 while (spec.empty()) {
798 std::getline(rl, spec);
800 rock_[i].read(dir, spec);
810 #endif // OPENRS_RESERVOIRPROPERTYCOMMON_IMPL_HEADER
double ntg(int cell_index) const
Read-access to ntg.
Definition: ReservoirPropertyCommon_impl.hpp:384
void setDensities(double d1, double d2)
Set densitities of both faces.
Definition: ReservoirPropertyCommon_impl.hpp:343
double cflFactorGravity() const
A factor useful in gravity cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:452
double viscositySecondPhase() const
Viscosity of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:357
PermTensor permeability(int cell_index) const
Read-access to permeability.
Definition: ReservoirPropertyCommon_impl.hpp:404
void assignRockTable(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:741
double densitySecondPhase() const
Density of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:371
Definition: BlackoilFluid.hpp:31
double porosity(int cell_index) const
Read-access to porosity.
Definition: ReservoirPropertyCommon_impl.hpp:378
double swcr(int cell_index) const
Read-access to swcr.
Definition: ReservoirPropertyCommon_impl.hpp:390
double cflFactorCapillary() const
A factor useful in gravity cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:459
void readRocks(const std::string &rock_list_file)
Definition: ReservoirPropertyCommon_impl.hpp:784
double viscosityFirstPhase() const
Viscosity of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:350
SharedCMatrix SharedPermTensor
Tensor type for read and write access to permeability.
Definition: ReservoirPropertyCommon.hpp:66
double s_min(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:491
PermeabilityKind
Enum for the kind of permeability field originally retrieved.
Definition: ReservoirPropertyCommon.hpp:50
double capillaryPressure(int cell_index, double saturation) const
Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:465
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:274
double densityDifference() const
Difference of densities.
Definition: ReservoirPropertyCommon_impl.hpp:438
double densityFirstPhase() const
Density of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:364
void assignPermeability(Opm::DeckConstPtr deck, const std::vector< int > &global_cell, const double perm_threshold)
Definition: ReservoirPropertyCommon_impl.hpp:681
Definition: ReservoirPropertyCommon.hpp:50
double s_max(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:504
double saturationFromCapillaryPressure(int cell_index, double cap_press) const
Inverse of the capillary pressure function.
Definition: ReservoirPropertyCommon_impl.hpp:517
Definition: ReservoirPropertyCommon.hpp:50
double cflFactor() const
A factor useful in cfl computations.
Definition: ReservoirPropertyCommon_impl.hpp:445
RPImpl & asImpl()
Definition: ReservoirPropertyCommon_impl.hpp:582
Definition: ReservoirPropertyCommon.hpp:50
void assignPorosity(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:591
void phaseDensities(int, Vector &density) const
Densities for both phases.
Definition: ReservoirPropertyCommon_impl.hpp:429
void assignSOWCR(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:658
ImmutableCMatrix PermTensor
Tensor type for read-only access to permeability.
Definition: ReservoirPropertyCommon.hpp:62
ReservoirPropertyCommon()
Default constructor.
Definition: ReservoirPropertyCommon_impl.hpp:256
Definition: ReservoirPropertyCommon.hpp:50
void assignNTG(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:613
void assignSWCR(Opm::DeckConstPtr deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:635
void setViscosities(double v1, double v2)
Set viscosities of both faces.
Definition: ReservoirPropertyCommon_impl.hpp:336
Definition: ReservoirPropertyCommon.hpp:50
double sowcr(int cell_index) const
Read-access to sowcr.
Definition: ReservoirPropertyCommon_impl.hpp:396
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
SharedPermTensor permeabilityModifiable(int cell_index)
Read- and write-access to permeability. Use with caution.
Definition: ReservoirPropertyCommon_impl.hpp:415
void writeSintefLegacyFormat(const std::string &grid_prefix) const
Write permeability and porosity in the Sintef legacy format.
Definition: ReservoirPropertyCommon_impl.hpp:531
double capillaryPressureDeriv(int c, double s) const
Derivative of Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:478