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,
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>
617 void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignNTG(
const Opm::Deck& deck,
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>
641 void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignSWCR(
const Opm::Deck& deck,
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>
666 void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignSOWCR(
const Opm::Deck& deck,
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>
691 void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignPermeability(
const Opm::Deck& deck,
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) {
733 SharedPermTensor K(dim, dim, &permeability_[off]);
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>
753 void ReservoirPropertyCommon<dim, RPImpl, RockType>::assignRockTable(
const Opm::Deck& deck,
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>
800 void ReservoirPropertyCommon<dim, RPImpl, RockType>::readRocks(
const std::string& rock_list_file)
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);