36#ifndef OPENRS_RESERVOIRPROPERTYCOMMON_IMPL_HEADER
37#define OPENRS_RESERVOIRPROPERTYCOMMON_IMPL_HEADER
39#include <opm/output/eclipse/EclipseGridInspector.hpp>
40#include <opm/input/eclipse/Deck/Deck.hpp>
41#include <opm/input/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]; }
180 std::vector<
const std::vector<double>*>& tensor,
181 std::array<int,9>& kmap)
185 OPM_THROW(std::runtime_error,
"Invalid set of permeability fields given.");
187 assert (tensor.size() == 1);
188 for (
int i = 0; i < 9; ++i) { kmap[i] = 0; }
196 if (
deck.hasKeyword(
"PERMX" )) {
197 kmap[xx] = tensor.size();
198 tensor.push_back(&
deck[
"PERMX"].back().getSIDoubleData());
200 setScalarPermIfNeeded(kmap, xx, yy, zz);
202 if (
deck.hasKeyword(
"PERMXY")) {
203 kmap[xy] = kmap[yx] = tensor.size();
204 tensor.push_back(&
deck[
"PERMXY"].back().getSIDoubleData());
206 if (
deck.hasKeyword(
"PERMXZ")) {
207 kmap[xz] = kmap[zx] = tensor.size();
208 tensor.push_back(&
deck[
"PERMXZ"].back().getSIDoubleData());
213 if (
deck.hasKeyword(
"PERMYX")) {
214 kmap[yx] = kmap[xy] = tensor.size();
215 tensor.push_back(&
deck[
"PERMYX"].back().getSIDoubleData());
217 if (
deck.hasKeyword(
"PERMY" )) {
218 kmap[yy] = tensor.size();
219 tensor.push_back(&
deck[
"PERMY"].back().getSIDoubleData());
221 setScalarPermIfNeeded(kmap, yy, zz, xx);
223 if (
deck.hasKeyword(
"PERMYZ")) {
224 kmap[yz] = kmap[zy] = tensor.size();
225 tensor.push_back(&
deck[
"PERMYZ"].back().getSIDoubleData());
230 if (
deck.hasKeyword(
"PERMZX")) {
231 kmap[zx] = kmap[xz] = tensor.size();
232 tensor.push_back(&
deck[
"PERMZX"].back().getSIDoubleData());
234 if (
deck.hasKeyword(
"PERMZY")) {
235 kmap[zy] = kmap[yz] = tensor.size();
236 tensor.push_back(&
deck[
"PERMZY"].back().getSIDoubleData());
238 if (
deck.hasKeyword(
"PERMZ" )) {
239 kmap[zz] = tensor.size();
240 tensor.push_back(&
deck[
"PERMZ"].back().getSIDoubleData());
242 setScalarPermIfNeeded(kmap, zz, xx, yy);
257 template <
int dim,
class RPImpl,
class RockType>
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),
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),
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,
285 static_assert(dim == 3,
"");
287 permfield_valid_.assign(global_cell.size(),
288 std::vector<unsigned char>::value_type(0));
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);
297 if (rock_list_filename) {
298 readRocks(*rock_list_filename);
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);
314 asImpl().computeCflFactors();
319 template <
int dim,
class RPImpl,
class RockType>
321 const double uniform_poro,
322 const double uniform_perm)
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) {
329 for (
int dd = 0; dd < dim; ++dd) {
330 K(dd, dd) = uniform_perm;
333 cell_to_rock_.assign(num_cells, 0);
334 asImpl().computeCflFactors();
337 template <
int dim,
class RPImpl,
class RockType>
344 template <
int dim,
class RPImpl,
class RockType>
351 template <
int dim,
class RPImpl,
class RockType>
358 template <
int dim,
class RPImpl,
class RockType>
365 template <
int dim,
class RPImpl,
class RockType>
372 template <
int dim,
class RPImpl,
class RockType>
379 template <
int dim,
class RPImpl,
class RockType>
382 return porosity_[cell_index];
385 template <
int dim,
class RPImpl,
class RockType>
388 return ntg_[cell_index];
391 template <
int dim,
class RPImpl,
class RockType>
394 return swcr_[cell_index];
397 template <
int dim,
class RPImpl,
class RockType>
400 return sowcr_[cell_index];
404 template <
int dim,
class RPImpl,
class RockType>
408 assert (permfield_valid_[cell_index]);
410 const PermTensor K(dim, dim, &permeability_[dim*dim*cell_index]);
415 template <
int dim,
class RPImpl,
class RockType>
423 permfield_valid_[cell_index] = std::vector<unsigned char>::value_type(1);
429 template <
int dim,
class RPImpl,
class RockType>
430 template<
class Vector>
433 assert (density.size() >= NumberOfPhases);
434 density[0] = densityFirstPhase();
435 density[1] = densitySecondPhase();
439 template <
int dim,
class RPImpl,
class RockType>
442 return density1_ - density2_;
446 template <
int dim,
class RPImpl,
class RockType>
453 template <
int dim,
class RPImpl,
class RockType>
456 return cfl_factor_gravity_;
460 template <
int dim,
class RPImpl,
class RockType>
463 return cfl_factor_capillary_;
466 template <
int dim,
class RPImpl,
class RockType>
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);
475 return 1e5*(1-saturation);
479 template <
int dim,
class RPImpl,
class RockType>
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);
492 template <
int dim,
class RPImpl,
class RockType>
495 if (rock_.size() > 0) {
496 int r = cell_to_rock_[cell_index];
497 return rock_[r].s_min();
505 template <
int dim,
class RPImpl,
class RockType>
508 if (rock_.size() > 0) {
509 int r = cell_to_rock_[cell_index];
510 return rock_[r].s_max();
518 template <
int dim,
class RPImpl,
class RockType>
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);
532 template <
int dim,
class RPImpl,
class RockType>
535 int num_cells = porosity_.size();
538 std::string filename = grid_prefix +
"-poro.dat";
539 std::ofstream file(filename.c_str());
541 OPM_THROW(std::runtime_error,
"Could not open file " + filename);
543 file << num_cells <<
'\n';
544 std::copy(porosity_.begin(), porosity_.end(), std::ostream_iterator<double>(file,
"\n"));
548 std::string filename = grid_prefix +
"-perm.dat";
549 std::ofstream file(filename.c_str());
551 OPM_THROW(std::runtime_error,
"Could not open file " + filename);
553 file << num_cells <<
'\n';
554 switch (permeability_kind_) {
556 std::copy(permeability_.begin(), permeability_.end(), std::ostream_iterator<double>(file,
"\n"));
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] <<
' ';
569 for (
int c = 0; c < num_cells; ++c) {
570 file << permeability_[c*dim*dim] <<
'\n';
574 OPM_THROW(std::runtime_error,
"Cannot write invalid permeability.");
583 template <
int dim,
class RPImpl,
class RockType>
586 return static_cast<RPImpl&
>(*this);
592 template <
int dim,
class RPImpl,
class RockType>
594 const std::vector<int>& global_cell)
596 porosity_.assign(global_cell.size(), 1.0);
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,
605 "PORO field must have the same size as the "
606 "logical cartesian size of the grid: " +
607 std::to_string(poro.size()) +
" != " +
608 std::to_string(num_global_cells));
610 for (
int c = 0; c < int(porosity_.size()); ++c) {
611 porosity_[c] = poro[global_cell[c]];
616 template <
int dim,
class RPImpl,
class RockType>
618 const std::vector<int>& global_cell)
620 ntg_.assign(global_cell.size(), 1.0);
622 if (
deck.hasKeyword(
"NTG")) {
623 Opm::EclipseGridInspector insp(
deck);
624 std::array<int, 3> dims = insp.gridSize();
625 int num_global_cells = dims[0]*dims[1]*dims[2];
626 const std::vector<double>& ntg =
deck[
"NTG"].back().getSIDoubleData();
627 if (
int(ntg.size()) != num_global_cells) {
628 OPM_THROW(std::runtime_error,
629 "NTG field must have the same size as the "
630 "logical cartesian size of the grid: " +
631 std::to_string(ntg.size()) +
" != " +
632 std::to_string(num_global_cells));
634 for (
int c = 0; c < int(ntg_.size()); ++c) {
635 ntg_[c] = ntg[global_cell[c]];
640 template <
int dim,
class RPImpl,
class RockType>
642 const std::vector<int>& global_cell)
644 swcr_.assign(global_cell.size(), 0.0);
646 if (
deck.hasKeyword(
"SWCR")) {
647 Opm::EclipseGridInspector insp(
deck);
648 std::array<int, 3> dims = insp.gridSize();
649 int num_global_cells = dims[0]*dims[1]*dims[2];
650 const std::vector<double>& swcr =
651 deck[
"SWCR"].back().getSIDoubleData();
652 if (
int(swcr.size()) != num_global_cells) {
653 OPM_THROW(std::runtime_error,
654 "SWCR field must have the same size as the "
655 "logical cartesian size of the grid: " +
656 std::to_string(swcr.size()) +
" != " +
657 std::to_string(num_global_cells));
659 for (
int c = 0; c < int(swcr_.size()); ++c) {
660 swcr_[c] = swcr[global_cell[c]];
665 template <
int dim,
class RPImpl,
class RockType>
667 const std::vector<int>& global_cell)
669 sowcr_.assign(global_cell.size(), 0.0);
671 if (
deck.hasKeyword(
"SOWCR")) {
672 Opm::EclipseGridInspector insp(
deck);
673 std::array<int, 3> dims = insp.gridSize();
674 int num_global_cells = dims[0]*dims[1]*dims[2];
675 const std::vector<double>& sowcr =
676 deck[
"SOWCR"].back().getSIDoubleData();
677 if (
int(sowcr.size()) != num_global_cells) {
678 OPM_THROW(std::runtime_error,
679 "SOWCR field must have the same size as the "
680 "logical cartesian size of the grid: " +
681 std::to_string(sowcr.size()) +
" != " +
682 std::to_string(num_global_cells));
684 for (
int c = 0; c < int(sowcr_.size()); ++c) {
685 sowcr_[c] = sowcr[global_cell[c]];
690 template <
int dim,
class RPImpl,
class RockType>
692 const std::vector<int>& global_cell,
693 double perm_threshold)
695 Opm::EclipseGridInspector insp(
deck);
696 std::array<int, 3> dims = insp.gridSize();
697 int num_global_cells = dims[0]*dims[1]*dims[2];
698 assert (num_global_cells > 0);
700 permeability_.assign(dim * dim * global_cell.size(), 0.0);
702 std::vector<const std::vector<double>*> tensor;
705 const std::vector<double>
zero(num_global_cells, 0.0);
706 tensor.push_back(&
zero);
708 static_assert(dim == 3,
"");
709 std::array<int,9> kmap;
710 permeability_kind_ = fillTensor(
deck, tensor, kmap);
711 for (
int i = 1; i < int(tensor.size()); ++i) {
712 if (
int(tensor[i]->size()) != num_global_cells) {
713 OPM_THROW(std::runtime_error,
714 "All permeability fields must have the same size as the "
715 "logical cartesian size of the grid: " +
716 std::to_string(tensor[i]->size()) +
" != " +
717 std::to_string(num_global_cells));
728 if (tensor.size() > 1) {
729 const int nc = global_cell.size();
732 for (
int c = 0; c < nc; ++c, off += dim*dim) {
735 const int glob = global_cell[c];
737 for (
int i = 0; i < dim; ++i) {
738 for (
int j = 0;
j < dim; ++
j, ++kix) {
739 K(i,
j) = (*tensor[kmap[kix]])[glob];
741 K(i,i) = std::max(K(i,i), perm_threshold);
744 permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
752 template <
int dim,
class RPImpl,
class RockType>
754 const std::vector<int>& global_cell)
756 const int nc = global_cell.size();
758 cell_to_rock_.assign(nc, 0);
760 if (
deck.hasKeyword(
"SATNUM")) {
761 Opm::EclipseGridInspector insp(
deck);
762 std::array<int, 3> dims = insp.gridSize();
763 int num_global_cells = dims[0]*dims[1]*dims[2];
764 const std::vector<int>&
satnum =
deck[
"SATNUM"].back().getIntData();
765 if (
int(
satnum.size()) != num_global_cells) {
766 OPM_THROW(std::runtime_error,
767 "SATNUM field must have the same size as the "
768 "logical cartesian size of the grid: " +
769 std::to_string(
satnum.size()) +
" != " +
770 std::to_string(num_global_cells));
772 for (
int c = 0; c < nc; ++c) {
774 cell_to_rock_[c] =
satnum[global_cell[c]] - 1;
777 else if (
deck.hasKeyword(
"ROCKTYPE")) {
778 Opm::EclipseGridInspector insp(
deck);
779 std::array<int, 3> dims = insp.gridSize();
780 int num_global_cells = dims[0]*dims[1]*dims[2];
781 const std::vector<int>&
satnum =
deck[
"ROCKTYPE"].back().getIntData();
782 if (
int(
satnum.size()) != num_global_cells) {
783 OPM_THROW(std::runtime_error,
784 "ROCKTYPE field must have the same size as the "
785 "logical cartesian size of the grid: " +
786 std::to_string(
satnum.size()) +
" != " +
787 std::to_string(num_global_cells));
789 for (
int c = 0; c < nc; ++c) {
791 cell_to_rock_[c] =
satnum[global_cell[c]] - 1;
799 template <
int dim,
class RPImpl,
class RockType>
802 std::ifstream rl(rock_list_file.c_str());
804 OPM_THROW(std::runtime_error,
"Could not open file " + rock_list_file);
808 assert(num_rocks >= 1);
809 rock_.resize(num_rocks);
810 std::string dir(rock_list_file.begin(), rock_list_file.begin() + rock_list_file.find_last_of(
'/') + 1);
811 for (
int i = 0; i < num_rocks; ++i) {
813 while (spec.empty()) {
814 std::getline(rl, spec);
816 rock_[i].read(dir, spec);
double porosity(int cell_index) const
Read-access to porosity.
Definition: ReservoirPropertyCommon_impl.hpp:380
ReservoirPropertyCommon()
Default constructor.
Definition: ReservoirPropertyCommon_impl.hpp:258
void assignSWCR(const Opm::Deck &deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:641
void assignPermeability(const Opm::Deck &deck, const std::vector< int > &global_cell, const double perm_threshold)
Definition: ReservoirPropertyCommon_impl.hpp:691
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
void assignRockTable(const Opm::Deck &deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:753
void assignPorosity(const Opm::Deck &deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:593
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
void assignNTG(const Opm::Deck &deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:617
void readRocks(const std::string &rock_list_file)
Definition: ReservoirPropertyCommon_impl.hpp:800
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
RPImpl & asImpl()
Definition: ReservoirPropertyCommon_impl.hpp:584
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. Use with caution.
Definition: ReservoirPropertyCommon_impl.hpp:417
void assignSOWCR(const Opm::Deck &deck, const std::vector< int > &global_cell)
Definition: ReservoirPropertyCommon_impl.hpp:666
double sowcr(int cell_index) const
Read-access to sowcr.
Definition: ReservoirPropertyCommon_impl.hpp:398
int j
Definition: elasticity_upscale_impl.hpp:658
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
std::vector< int > satnum
Definition: elasticity_upscale_impl.hpp:580
auto deck
Definition: elasticity_upscale_impl.hpp:590
Definition: ImplicitAssembly.hpp:43
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:602