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) {
93 retval = DiagonalPerm;
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)
184 if (kind == Invalid) {
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),
270 permeability_kind_(Invalid)
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>
593 void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignPorosity(
const Opm::Deck& deck,
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,
"PORO field must have the same size as the "
605 "logical cartesian size of the grid: "
606 << poro.size() <<
" != " << num_global_cells);
608 for (
int c = 0; c < int(porosity_.size()); ++c) {
609 porosity_[c] = poro[global_cell[c]];
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)
618 ntg_.assign(global_cell.size(), 1.0);
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);
630 for (
int c = 0; c < int(ntg_.size()); ++c) {
631 ntg_[c] = ntg[global_cell[c]];
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)
640 swcr_.assign(global_cell.size(), 0.0);
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);
653 for (
int c = 0; c < int(swcr_.size()); ++c) {
654 swcr_[c] = swcr[global_cell[c]];
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)
663 sowcr_.assign(global_cell.size(), 0.0);
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);
676 for (
int c = 0; c < int(sowcr_.size()); ++c) {
677 sowcr_[c] = sowcr[global_cell[c]];
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)
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);
692 permeability_.assign(dim * dim * global_cell.size(), 0.0);
694 std::vector<const std::vector<double>*> tensor;
697 const std::vector<double>
zero(num_global_cells, 0.0);
698 tensor.push_back(&
zero);
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);
718 if (tensor.size() > 1) {
719 const int nc = global_cell.size();
722 for (
int c = 0; c < nc; ++c, off += dim*dim) {
723 SharedPermTensor K(dim, dim, &permeability_[off]);
725 const int glob = global_cell[c];
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];
731 K(i,i) = std::max(K(i,i), perm_threshold);
734 permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
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)
746 const int nc = global_cell.size();
748 cell_to_rock_.assign(nc, 0);
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);
760 for (
int c = 0; c < nc; ++c) {
762 cell_to_rock_[c] = satnum[global_cell[c]] - 1;
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);
775 for (
int c = 0; c < nc; ++c) {
777 cell_to_rock_[c] = satnum[global_cell[c]] - 1;
785 template <
int dim,
class RPImpl,
class RockType>
786 void ReservoirPropertyCommon<dim, RPImpl, RockType>::readRocks(
const std::string& rock_list_file)
788 std::ifstream rl(rock_list_file.c_str());
790 OPM_THROW(std::runtime_error,
"Could not open file " << rock_list_file);
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) {
799 while (spec.empty()) {
800 std::getline(rl, spec);
802 rock_[i].read(dir, spec);
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