36#ifndef OPENRS_ROCKJFUNC_HEADER
37#define OPENRS_ROCKJFUNC_HEADER
39#include <dune/common/fvector.hh>
40#include <opm/core/utility/NonuniformTableLinear.hpp>
56 : use_jfunction_scaling_(true), sigma_cos_theta_(1.0)
62 use_jfunction_scaling_ = use_j;
67 sigma_cos_theta_ = sigma*std::cos(theta);
70 void krw(
const double saturation,
double& krw_value)
const
72 krw_value = krw_(saturation);
75 void kro(
const double saturation,
double& kro_value)
const
77 kro_value = kro_(saturation);
80 void dkrw(
const double saturation,
double& dkrw_value)
const
82 dkrw_value = krw_.derivative(saturation);
85 void dkro(
const double saturation,
double& dkro_value)
const
87 dkro_value = kro_.derivative(saturation);
99 template <
template <
class>
class SP,
class OP>
100 double capPress(
const FullMatrix<double, SP, OP>& perm,
const double poro,
const double saturation)
const
102 if (use_jfunction_scaling_) {
106 double sqrt_k_phi = std::sqrt(
trace(perm)/(perm.numRows()*poro));
107 return Jfunc_(saturation)*sigma_cos_theta_/sqrt_k_phi;
110 return Jfunc_(saturation);
114 template <
template <
class>
class SP,
class OP>
115 double capPressDeriv(
const FullMatrix<double, SP, OP>& perm,
const double poro,
const double saturation)
const
117 if (use_jfunction_scaling_) {
121 double sqrt_k_phi = std::sqrt(
trace(perm)/(perm.numRows()*poro));
122 return Jfunc_.derivative(saturation)*sigma_cos_theta_/sqrt_k_phi;
125 return Jfunc_.derivative(saturation);
129 template <
template <
class>
class SP,
class OP>
130 double satFromCapPress(
const FullMatrix<double, SP, OP>& perm,
const double poro,
const double cap_press)
const
133 if (use_jfunction_scaling_) {
137 double sqrt_k_phi = std::sqrt(
trace(perm)/(perm.numRows()*poro));
138 s = Jfunc_.inverse(cap_press*sqrt_k_phi/sigma_cos_theta_);
141 s = Jfunc_.inverse(cap_press);
143 s = std::min(s_max_, std::max(s_min_, s));
147 void read(
const std::string& directory,
const std::string& specification)
150 std::istringstream specstream(specification);
151 std::string rockname;
152 specstream >> rockname;
153 std::string rockfilename = directory + rockname;
154 std::ifstream rock_stream(rockfilename.c_str());
156 OPM_THROW(std::runtime_error,
"Could not open file " << rockfilename);
158 readStatoilFormat(rock_stream);
162 void readStatoilFormat(std::istream& is)
166 while ((c ==
'-' || c ==
'#') && is) {
168 std::string commentline;
169 std::getline(is, commentline);
175 std::string commentline;
176 std::getline(is, commentline);
184 OPM_THROW(std::runtime_error,
"Something went wrong while skipping optional comment header of rock file.");
186 typedef Dune::FieldVector<double, 4> Data;
187 std::istream_iterator<Data> start(is);
188 std::istream_iterator<Data> end;
189 std::vector<Data> data(start, end);
191 OPM_THROW(std::runtime_error,
"Reading stopped but we're not at eof: something went wrong reading data of rock file.");
193 std::vector<double> svals,
krw,
kro, Jfunc;
194 for (
int i = 0; i < int(data.size()); ++i) {
195 svals.push_back(data[i][0]);
196 krw.push_back(data[i][1]);
197 kro.push_back(data[i][2]);
198 Jfunc.push_back(data[i][3]);
202 std::cout <<
"Warning: krw table were modified to go to zero at the beginning." << std::endl;
204 if (
kro.back() != 0.0) {
206 std::cout <<
"Warning: kro table were modified to go to zero at the end." << std::endl;
208 s_min_ = svals.front();
209 s_max_ = svals.back();
210 krw_ = TabFunc(svals,
krw);
211 kro_ = TabFunc(svals,
kro);
212 Jfunc_ = TabFunc(svals, Jfunc);
213 std::vector<double> invJfunc(Jfunc);
214 std::reverse(invJfunc.begin(), invJfunc.end());
215 std::vector<double> invsvals(svals);
216 std::reverse(invsvals.begin(), invsvals.end());
217 invJfunc_ = TabFunc(invJfunc, invsvals);
220 typedef NonuniformTableLinear<double> TabFunc;
225 bool use_jfunction_scaling_;
226 double sigma_cos_theta_;
Definition: RockJfunc.hpp:53
double capPressDeriv(const FullMatrix< double, SP, OP > &perm, const double poro, const double saturation) const
Definition: RockJfunc.hpp:115
double satFromCapPress(const FullMatrix< double, SP, OP > &perm, const double poro, const double cap_press) const
Definition: RockJfunc.hpp:130
void kro(const double saturation, double &kro_value) const
Definition: RockJfunc.hpp:75
void setUseJfunctionScaling(const bool use_j)
Definition: RockJfunc.hpp:60
RockJfunc()
Definition: RockJfunc.hpp:55
void dkrw(const double saturation, double &dkrw_value) const
Definition: RockJfunc.hpp:80
double s_min() const
Definition: RockJfunc.hpp:90
double s_max() const
Definition: RockJfunc.hpp:94
void setSigmaAndTheta(const double sigma, const double theta)
Definition: RockJfunc.hpp:65
double capPress(const FullMatrix< double, SP, OP > &perm, const double poro, const double saturation) const
Definition: RockJfunc.hpp:100
void krw(const double saturation, double &krw_value) const
Definition: RockJfunc.hpp:70
void dkro(const double saturation, double &dkro_value) const
Definition: RockJfunc.hpp:85
void read(const std::string &directory, const std::string &specification)
Definition: RockJfunc.hpp:147
Definition: BlackoilFluid.hpp:32
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:639