diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp index b290f9be7..fd072798c 100644 --- a/opm/core/props/rock/RockCompressibility.cpp +++ b/opm/core/props/rock/RockCompressibility.cpp @@ -42,12 +42,16 @@ namespace Opm { if (deck.hasField("ROCKTAB")) { const table_t& rt = deck.getROCKTAB().rocktab_; - int n = rt[0][0].size(); + if (rt.size() != 1) { + THROW("Can only handle a single region in ROCKTAB."); + } + const int n = rt[0][0].size(); p_.resize(n); poromult_.resize(n); for (int i = 0; i < n; ++i) { p_[i] = rt[0][0][i]; poromult_[i] = rt[0][1][i]; + transmult_[i] = rt[0][2][i]; } } else if (deck.hasField("ROCK")) { const ROCK& r = deck.getROCK(); @@ -70,11 +74,39 @@ namespace Opm const double cpnorm = rock_comp_*(pressure - pref_); return (1.0 + cpnorm + 0.5*cpnorm*cpnorm); } else { - // return Opm::linearInterpolation(p_, poromult_, pressure); return Opm::linearInterpolation(p_, poromult_, pressure); } } + double RockCompressibility::poroMultDeriv(double pressure) const + { + if (p_.empty()) { + // Approximating poro multiplier with a quadratic curve, + // we must use its derivative. + return rock_comp_ + 2 * rock_comp_ * rock_comp_ * (pressure - pref_); + } else { + return Opm::linearInterpolationDerivative(p_, poromult_, pressure); + } + } + + double RockCompressibility::transMult(double pressure) const + { + if (p_.empty()) { + return 1.0; + } else { + return Opm::linearInterpolation(p_, transmult_, pressure); + } + } + + double RockCompressibility::transMultDeriv(double pressure) const + { + if (p_.empty()) { + return 0.0; + } else { + return Opm::linearInterpolationDerivative(p_, transmult_, pressure); + } + } + double RockCompressibility::rockComp(double pressure) const { if (p_.empty()) { diff --git a/opm/core/props/rock/RockCompressibility.hpp b/opm/core/props/rock/RockCompressibility.hpp index a5c9a983b..03567caed 100644 --- a/opm/core/props/rock/RockCompressibility.hpp +++ b/opm/core/props/rock/RockCompressibility.hpp @@ -47,12 +47,22 @@ namespace Opm /// Porosity multiplier. double poroMult(double pressure) const; + /// Derivative of porosity multiplier with respect to pressure. + double poroMultDeriv(double pressure) const; + + /// Transmissibility multiplier. + double transMult(double pressure) const; + + /// Derivative of transmissibility multiplier with respect to pressure. + double transMultDeriv(double pressure) const; + /// Rock compressibility = (d poro / d p)*(1 / poro). double rockComp(double pressure) const; private: std::vector p_; std::vector poromult_; + std::vector transmult_; double pref_; double rock_comp_; };