All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
PiecewiseLinearTwoPhaseMaterial.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_PIECEWISE_LINEAR_TWO_PHASE_MATERIAL_HPP
28 #define OPM_PIECEWISE_LINEAR_TWO_PHASE_MATERIAL_HPP
29 
31 
32 #include <opm/common/ErrorMacros.hpp>
33 #include <opm/common/Exceptions.hpp>
35 
36 #include <algorithm>
37 #include <cmath>
38 #include <cassert>
39 
40 namespace Opm {
51 template <class TraitsT, class ParamsT = PiecewiseLinearTwoPhaseMaterialParams<TraitsT> >
52 class PiecewiseLinearTwoPhaseMaterial : public TraitsT
53 {
54  typedef typename ParamsT::ValueVector ValueVector;
55 
56 public:
58  typedef TraitsT Traits;
59 
61  typedef ParamsT Params;
62 
64  typedef typename Traits::Scalar Scalar;
65 
67  static const int numPhases = Traits::numPhases;
68  static_assert(numPhases == 2,
69  "The piecewise linear two-phase capillary pressure law only"
70  "applies to the case of two fluid phases");
71 
74  static const bool implementsTwoPhaseApi = true;
75 
78  static const bool implementsTwoPhaseSatApi = true;
79 
82  static const bool isSaturationDependent = true;
83 
86  static const bool isPressureDependent = false;
87 
90  static const bool isTemperatureDependent = false;
91 
94  static const bool isCompositionDependent = false;
95 
99  template <class Container, class FluidState>
100  static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
101  {
102  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
103 
104  values[Traits::wettingPhaseIdx] = 0.0; // reference phase
105  values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
106  }
107 
112  template <class Container, class FluidState>
113  static void saturations(Container& /* values */, const Params& /* params */, const FluidState& /* fs */)
114  { OPM_THROW(std::logic_error, "Not implemented: saturations()"); }
115 
119  template <class Container, class FluidState>
120  static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
121  {
122  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
123 
124  values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
125  values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
126  }
127 
131  template <class FluidState, class Evaluation = typename FluidState::Scalar>
132  static Evaluation pcnw(const Params& params, const FluidState& fs)
133  {
134  const auto& Sw =
135  Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
136 
137  return twoPhaseSatPcnw(params, Sw);
138  }
139 
143  template <class Evaluation>
144  static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
145  { return eval_(params.SwPcwnSamples(), params.pcnwSamples(), Sw); }
146 
147  template <class Evaluation>
148  static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnw)
149  { return eval_(params.pcnwSamples(), params.SwPcwnSamples(), pcnw); }
150 
154  template <class FluidState, class Evaluation = typename FluidState::Scalar>
155  static Evaluation Sw(const Params& /* params */, const FluidState& /* fs */)
156  { OPM_THROW(std::logic_error, "Not implemented: Sw()"); }
157 
158  template <class Evaluation>
159  static Evaluation twoPhaseSatSw(const Params& /* params */, const Evaluation& /* pC */)
160  { OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); }
161 
166  template <class FluidState, class Evaluation = typename FluidState::Scalar>
167  static Evaluation Sn(const Params& params, const FluidState& fs)
168  { return 1 - Sw<FluidState, Scalar>(params, fs); }
169 
170  template <class Evaluation>
171  static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pC)
172  { return 1 - twoPhaseSatSw(params, pC); }
173 
178  template <class FluidState, class Evaluation = typename FluidState::Scalar>
179  static Evaluation krw(const Params& params, const FluidState& fs)
180  {
181  const auto& Sw =
182  Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
183 
184  return twoPhaseSatKrw(params, Sw);
185  }
186 
187  template <class Evaluation>
188  static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
189  { return eval_(params.SwKrwSamples(), params.krwSamples(), Sw); }
190 
191  template <class Evaluation>
192  static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
193  { return eval_(params.krwSamples(), params.SwKrwSamples(), krw); }
194 
199  template <class FluidState, class Evaluation = typename FluidState::Scalar>
200  static Evaluation krn(const Params& params, const FluidState& fs)
201  {
202  const auto& Sw =
203  Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
204 
205  return twoPhaseSatKrn(params, Sw);
206  }
207 
208  template <class Evaluation>
209  static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
210  { return eval_(params.SwKrnSamples(), params.krnSamples(), Sw); }
211 
212  template <class Evaluation>
213  static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
214  { return eval_(params.krnSamples(), params.SwKrnSamples(), krn); }
215 
216 private:
217  template <class Evaluation>
218  static Evaluation eval_(const ValueVector& xValues,
219  const ValueVector& yValues,
220  const Evaluation& x)
221  {
222  if (xValues.front() < xValues.back())
223  return evalAscending_(xValues, yValues, x);
224  return evalDescending_(xValues, yValues, x);
225  }
226 
227  template <class Evaluation>
228  static Evaluation evalAscending_(const ValueVector& xValues,
229  const ValueVector& yValues,
230  const Evaluation& x)
231  {
232  if (x <= xValues.front())
233  return yValues.front();
234  if (x >= xValues.back())
235  return yValues.back();
236 
237  size_t segIdx = findSegmentIndex_(xValues, Opm::scalarValue(x));
238 
239  Scalar x0 = xValues[segIdx];
240  Scalar x1 = xValues[segIdx + 1];
241 
242  Scalar y0 = yValues[segIdx];
243  Scalar y1 = yValues[segIdx + 1];
244 
245  Scalar m = (y1 - y0)/(x1 - x0);
246 
247  return y0 + (x - x0)*m;
248  }
249 
250  template <class Evaluation>
251  static Evaluation evalDescending_(const ValueVector& xValues,
252  const ValueVector& yValues,
253  const Evaluation& x)
254  {
255  if (x >= xValues.front())
256  return yValues.front();
257  if (x <= xValues.back())
258  return yValues.back();
259 
260  size_t segIdx = findSegmentIndexDescending_(xValues, Opm::scalarValue(x));
261 
262  Scalar x0 = xValues[segIdx];
263  Scalar x1 = xValues[segIdx + 1];
264 
265  Scalar y0 = yValues[segIdx];
266  Scalar y1 = yValues[segIdx + 1];
267 
268  Scalar m = (y1 - y0)/(x1 - x0);
269 
270  return y0 + (x - x0)*m;
271  }
272 
273  template <class Evaluation>
274  static Evaluation evalDeriv_(const ValueVector& xValues,
275  const ValueVector& yValues,
276  const Evaluation& x)
277  {
278  if (x <= xValues.front())
279  return 0.0;
280  if (x >= xValues.back())
281  return 0.0;
282 
283  size_t segIdx = findSegmentIndex_(xValues, Opm::scalarValue(x));
284 
285  Scalar x0 = xValues[segIdx];
286  Scalar x1 = xValues[segIdx + 1];
287 
288  Scalar y0 = yValues[segIdx];
289  Scalar y1 = yValues[segIdx + 1];
290 
291  return (y1 - y0)/(x1 - x0);
292  }
293 
294  static size_t findSegmentIndex_(const ValueVector& xValues, Scalar x)
295  {
296  assert(xValues.size() > 1); // we need at least two sampling points!
297  size_t n = xValues.size() - 1;
298  if (xValues.back() <= x)
299  return n - 1;
300  else if (x <= xValues.front())
301  return 0;
302 
303  // bisection
304  size_t lowIdx = 0, highIdx = n;
305  while (lowIdx + 1 < highIdx) {
306  size_t curIdx = (lowIdx + highIdx)/2;
307  if (xValues[curIdx] < x)
308  lowIdx = curIdx;
309  else
310  highIdx = curIdx;
311  }
312 
313  return lowIdx;
314  }
315 
316  static size_t findSegmentIndexDescending_(const ValueVector& xValues, Scalar x)
317  {
318  assert(xValues.size() > 1); // we need at least two sampling points!
319  size_t n = xValues.size() - 1;
320  if (x <= xValues.back())
321  return n;
322  else if (xValues.front() <= x)
323  return 0;
324 
325  // bisection
326  size_t lowIdx = 0, highIdx = n;
327  while (lowIdx + 1 < highIdx) {
328  size_t curIdx = (lowIdx + highIdx)/2;
329  if (xValues[curIdx] >= x)
330  lowIdx = curIdx;
331  else
332  highIdx = curIdx;
333  }
334 
335  return lowIdx;
336  }
337 };
338 } // namespace Opm
339 
340 #endif
static Evaluation Sn(const Params &params, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:167
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition...
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:94
static Evaluation krw(const Params &params, const FluidState &fs)
The relative permeability for the wetting phase of the porous medium.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:179
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:78
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:82
static void saturations(Container &, const Params &, const FluidState &)
The saturations of the fluid phases starting from their pressure differences.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:113
ParamsT Params
The type of the parameter objects for this law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:61
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:90
static Evaluation pcnw(const Params &params, const FluidState &fs)
The capillary pressure-saturation curve.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:132
static Evaluation krn(const Params &params, const FluidState &fs)
The relative permeability for the non-wetting phase of the porous medium.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:200
Implementation of a tabulated, piecewise linear capillary pressure law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:52
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation &Sw)
The saturation-capillary pressure curve.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:144
static Evaluation Sw(const Params &, const FluidState &)
The saturation-capillary pressure curve.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:155
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
The capillary pressure-saturation curve.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:100
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
The relative permeabilities.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:120
static const int numPhases
The number of fluid phases.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:67
Specification of the material parameters for a two-phase material law which uses a table and piecewis...
Traits::Scalar Scalar
The type of the scalar values for this law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:64
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure...
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:86
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:74
TraitsT Traits
The traits class for this material law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:58