36#ifndef OPENRS_ROCKANISOTROPICRELPERM_HEADER
37#define OPENRS_ROCKANISOTROPICRELPERM_HEADER
40#include <dune/common/fvector.hh>
41#include <opm/core/utility/NonuniformTableLinear.hpp>
60 OPM_THROW(std::runtime_error,
"RockAnisotropicRelperm cannot use J-scaling.");
66 OPM_THROW(std::runtime_error,
"RockAnisotropicRelperm cannot accept sigma and theta arguments.");
70 template <
template <
class>
class SP,
class OP>
71 void kr(
const int phase_index,
const double saturation, FullMatrix<double, SP, OP>& kr_value)
const
74 kr_value(0,0) = krxx_[phase_index](saturation);
75 kr_value(1,1) = kryy_[phase_index](saturation);
76 kr_value(2,2) = krzz_[phase_index](saturation);
79 template <
template <
class>
class SP,
class OP>
80 double capPress(
const FullMatrix<double, SP, OP>& ,
const double ,
const double saturation)
const
82 return cap_press_(saturation);
85 template <
template <
class>
class SP,
class OP>
86 double satFromCapPress(
const FullMatrix<double, SP, OP>& ,
const double ,
const double cp)
const
88 return cap_press_.inverse(cp);
91 void read(
const std::string& directory,
const std::string& specification)
94 std::istringstream specstream(specification);
95 std::string rockname[2];
96 specstream >> rockname[0] >> rockname[1];
97 for (
int i = 0; i < 2; ++i) {
98 std::string rockfilename = directory + rockname[i];
99 std::ifstream rock_stream(rockfilename.c_str());
101 OPM_THROW(std::runtime_error,
"Could not open file " << rockfilename);
103 readAnisoFormat(rock_stream, i);
108 void readAnisoFormat(std::istream& is,
int phase)
111 while (is.peek() ==
'#') {
112 is.ignore(std::numeric_limits<int>::max(),
'\n');
114 typedef Dune::FieldVector<double, 5> Data;
115 std::istream_iterator<Data> start(is);
116 std::istream_iterator<Data> end;
117 std::vector<Data> data(start, end);
119 OPM_THROW(std::runtime_error,
"Reading stopped but we're not at eof: something went wrong reading data.");
123 for (
int i = 2; i < 5; ++i) {
124 if (phase == 0 && data[0][i] != 0.0) {
127 }
else if (phase == 1 && data.back()[i] != 0.0) {
128 data.back()[i] = 0.0;
133 std::cout <<
"Warning: Relperm tables were modified to go to zero at the beginning/end." << std::endl;
135 std::vector<double> pcow, svals, krxx, kryy, krzz;
136 for (
int i = 0; i < int(data.size()); ++i) {
137 pcow.push_back(data[i][0]);
138 svals.push_back(data[i][1]);
139 krxx.push_back(data[i][2]);
140 kryy.push_back(data[i][3]);
141 krzz.push_back(data[i][4]);
143 krxx_[phase] = TabFunc(svals, krxx);
144 kryy_[phase] = TabFunc(svals, kryy);
145 krzz_[phase] = TabFunc(svals, krzz);
147 cap_press_ = TabFunc(svals, pcow);
150 assert(cap_press_ == TabFunc(svals, pcow));
154 typedef NonuniformTableLinear<double> TabFunc;
Definition: RockAnisotropicRelperm.hpp:54
double satFromCapPress(const FullMatrix< double, SP, OP > &, const double, const double cp) const
Definition: RockAnisotropicRelperm.hpp:86
void setSigmaAndTheta(const double, const double)
Definition: RockAnisotropicRelperm.hpp:64
void kr(const int phase_index, const double saturation, FullMatrix< double, SP, OP > &kr_value) const
Definition: RockAnisotropicRelperm.hpp:71
double capPress(const FullMatrix< double, SP, OP > &, const double, const double saturation) const
Definition: RockAnisotropicRelperm.hpp:80
void setUseJfunctionScaling(const bool use_j)
Definition: RockAnisotropicRelperm.hpp:57
void read(const std::string &directory, const std::string &specification)
Definition: RockAnisotropicRelperm.hpp:91
Definition: BlackoilFluid.hpp:32
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603