00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef OPM_PIECEWISE_LINEAR_TWO_PHASE_MATERIAL_HPP
00028 #define OPM_PIECEWISE_LINEAR_TWO_PHASE_MATERIAL_HPP
00029
00030 #include "PiecewiseLinearTwoPhaseMaterialParams.hpp"
00031
00032 #include <opm/common/ErrorMacros.hpp>
00033 #include <opm/common/Exceptions.hpp>
00034 #include <opm/material/common/MathToolbox.hpp>
00035
00036 #include <algorithm>
00037 #include <cmath>
00038 #include <cassert>
00039
00040 namespace Opm {
00051 template <class TraitsT, class ParamsT = PiecewiseLinearTwoPhaseMaterialParams<TraitsT> >
00052 class PiecewiseLinearTwoPhaseMaterial : public TraitsT
00053 {
00054 typedef typename ParamsT::ValueVector ValueVector;
00055
00056 public:
00058 typedef TraitsT Traits;
00059
00061 typedef ParamsT Params;
00062
00064 typedef typename Traits::Scalar Scalar;
00065
00067 static const int numPhases = Traits::numPhases;
00068 static_assert(numPhases == 2,
00069 "The piecewise linear two-phase capillary pressure law only"
00070 "applies to the case of two fluid phases");
00071
00074 static const bool implementsTwoPhaseApi = true;
00075
00078 static const bool implementsTwoPhaseSatApi = true;
00079
00082 static const bool isSaturationDependent = true;
00083
00086 static const bool isPressureDependent = false;
00087
00090 static const bool isTemperatureDependent = false;
00091
00094 static const bool isCompositionDependent = false;
00095
00099 template <class Container, class FluidState>
00100 static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
00101 {
00102 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00103
00104 values[Traits::wettingPhaseIdx] = 0.0;
00105 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
00106 }
00107
00112 template <class Container, class FluidState>
00113 static void saturations(Container& , const Params& , const FluidState& )
00114 { OPM_THROW(std::logic_error, "Not implemented: saturations()"); }
00115
00119 template <class Container, class FluidState>
00120 static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
00121 {
00122 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00123
00124 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
00125 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
00126 }
00127
00131 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00132 static Evaluation pcnw(const Params& params, const FluidState& fs)
00133 {
00134 const auto& Sw =
00135 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00136
00137 return twoPhaseSatPcnw(params, Sw);
00138 }
00139
00143 template <class Evaluation>
00144 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
00145 { return eval_(params.SwPcwnSamples(), params.pcnwSamples(), Sw); }
00146
00147 template <class Evaluation>
00148 static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnw)
00149 { return eval_(params.pcnwSamples(), params.SwPcwnSamples(), pcnw); }
00150
00154 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00155 static Evaluation Sw(const Params& , const FluidState& )
00156 { OPM_THROW(std::logic_error, "Not implemented: Sw()"); }
00157
00158 template <class Evaluation>
00159 static Evaluation twoPhaseSatSw(const Params& , const Evaluation& )
00160 { OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); }
00161
00166 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00167 static Evaluation Sn(const Params& params, const FluidState& fs)
00168 { return 1 - Sw<FluidState, Scalar>(params, fs); }
00169
00170 template <class Evaluation>
00171 static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pC)
00172 { return 1 - twoPhaseSatSw(params, pC); }
00173
00178 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00179 static Evaluation krw(const Params& params, const FluidState& fs)
00180 {
00181 const auto& Sw =
00182 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00183
00184 return twoPhaseSatKrw(params, Sw);
00185 }
00186
00187 template <class Evaluation>
00188 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
00189 { return eval_(params.SwKrwSamples(), params.krwSamples(), Sw); }
00190
00191 template <class Evaluation>
00192 static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
00193 { return eval_(params.krwSamples(), params.SwKrwSamples(), krw); }
00194
00199 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00200 static Evaluation krn(const Params& params, const FluidState& fs)
00201 {
00202 const auto& Sw =
00203 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00204
00205 return twoPhaseSatKrn(params, Sw);
00206 }
00207
00208 template <class Evaluation>
00209 static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
00210 { return eval_(params.SwKrnSamples(), params.krnSamples(), Sw); }
00211
00212 template <class Evaluation>
00213 static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
00214 { return eval_(params.krnSamples(), params.SwKrnSamples(), krn); }
00215
00216 private:
00217 template <class Evaluation>
00218 static Evaluation eval_(const ValueVector& xValues,
00219 const ValueVector& yValues,
00220 const Evaluation& x)
00221 {
00222 if (xValues.front() < xValues.back())
00223 return evalAscending_(xValues, yValues, x);
00224 return evalDescending_(xValues, yValues, x);
00225 }
00226
00227 template <class Evaluation>
00228 static Evaluation evalAscending_(const ValueVector& xValues,
00229 const ValueVector& yValues,
00230 const Evaluation& x)
00231 {
00232 if (x <= xValues.front())
00233 return yValues.front();
00234 if (x >= xValues.back())
00235 return yValues.back();
00236
00237 size_t segIdx = findSegmentIndex_(xValues, Opm::scalarValue(x));
00238
00239 Scalar x0 = xValues[segIdx];
00240 Scalar x1 = xValues[segIdx + 1];
00241
00242 Scalar y0 = yValues[segIdx];
00243 Scalar y1 = yValues[segIdx + 1];
00244
00245 Scalar m = (y1 - y0)/(x1 - x0);
00246
00247 return y0 + (x - x0)*m;
00248 }
00249
00250 template <class Evaluation>
00251 static Evaluation evalDescending_(const ValueVector& xValues,
00252 const ValueVector& yValues,
00253 const Evaluation& x)
00254 {
00255 if (x >= xValues.front())
00256 return yValues.front();
00257 if (x <= xValues.back())
00258 return yValues.back();
00259
00260 size_t segIdx = findSegmentIndexDescending_(xValues, Opm::scalarValue(x));
00261
00262 Scalar x0 = xValues[segIdx];
00263 Scalar x1 = xValues[segIdx + 1];
00264
00265 Scalar y0 = yValues[segIdx];
00266 Scalar y1 = yValues[segIdx + 1];
00267
00268 Scalar m = (y1 - y0)/(x1 - x0);
00269
00270 return y0 + (x - x0)*m;
00271 }
00272
00273 template <class Evaluation>
00274 static Evaluation evalDeriv_(const ValueVector& xValues,
00275 const ValueVector& yValues,
00276 const Evaluation& x)
00277 {
00278 if (x <= xValues.front())
00279 return 0.0;
00280 if (x >= xValues.back())
00281 return 0.0;
00282
00283 size_t segIdx = findSegmentIndex_(xValues, Opm::scalarValue(x));
00284
00285 Scalar x0 = xValues[segIdx];
00286 Scalar x1 = xValues[segIdx + 1];
00287
00288 Scalar y0 = yValues[segIdx];
00289 Scalar y1 = yValues[segIdx + 1];
00290
00291 return (y1 - y0)/(x1 - x0);
00292 }
00293
00294 static size_t findSegmentIndex_(const ValueVector& xValues, Scalar x)
00295 {
00296 assert(xValues.size() > 1);
00297 size_t n = xValues.size() - 1;
00298 if (xValues.back() <= x)
00299 return n - 1;
00300 else if (x <= xValues.front())
00301 return 0;
00302
00303
00304 size_t lowIdx = 0, highIdx = n;
00305 while (lowIdx + 1 < highIdx) {
00306 size_t curIdx = (lowIdx + highIdx)/2;
00307 if (xValues[curIdx] < x)
00308 lowIdx = curIdx;
00309 else
00310 highIdx = curIdx;
00311 }
00312
00313 return lowIdx;
00314 }
00315
00316 static size_t findSegmentIndexDescending_(const ValueVector& xValues, Scalar x)
00317 {
00318 assert(xValues.size() > 1);
00319 size_t n = xValues.size() - 1;
00320 if (x <= xValues.back())
00321 return n;
00322 else if (xValues.front() <= x)
00323 return 0;
00324
00325
00326 size_t lowIdx = 0, highIdx = n;
00327 while (lowIdx + 1 < highIdx) {
00328 size_t curIdx = (lowIdx + highIdx)/2;
00329 if (xValues[curIdx] >= x)
00330 lowIdx = curIdx;
00331 else
00332 highIdx = curIdx;
00333 }
00334
00335 return lowIdx;
00336 }
00337 };
00338 }
00339
00340 #endif