My Project
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 
33 
34 #include <algorithm>
35 #include <cmath>
36 #include <cassert>
37 
38 namespace Opm {
49 template <class TraitsT, class ParamsT = PiecewiseLinearTwoPhaseMaterialParams<TraitsT> >
50 class PiecewiseLinearTwoPhaseMaterial : public TraitsT
51 {
52  typedef typename ParamsT::ValueVector ValueVector;
53 
54 public:
56  typedef TraitsT Traits;
57 
59  typedef ParamsT Params;
60 
62  typedef typename Traits::Scalar Scalar;
63 
65  static const int numPhases = Traits::numPhases;
66  static_assert(numPhases == 2,
67  "The piecewise linear two-phase capillary pressure law only"
68  "applies to the case of two fluid phases");
69 
72  static const bool implementsTwoPhaseApi = true;
73 
76  static const bool implementsTwoPhaseSatApi = true;
77 
80  static const bool isSaturationDependent = true;
81 
84  static const bool isPressureDependent = false;
85 
88  static const bool isTemperatureDependent = false;
89 
92  static const bool isCompositionDependent = false;
93 
97  template <class Container, class FluidState>
98  static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
99  {
100  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
101 
102  values[Traits::wettingPhaseIdx] = 0.0; // reference phase
103  values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
104  }
105 
110  template <class Container, class FluidState>
111  static void saturations(Container& /* values */, const Params& /* params */, const FluidState& /* fs */)
112  { throw std::logic_error("Not implemented: saturations()"); }
113 
117  template <class Container, class FluidState>
118  static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
119  {
120  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
121 
122  values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
123  values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
124  }
125 
129  template <class FluidState, class Evaluation = typename FluidState::Scalar>
130  static Evaluation pcnw(const Params& params, const FluidState& fs)
131  {
132  const auto& Sw =
133  decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
134 
135  return twoPhaseSatPcnw(params, Sw);
136  }
137 
141  template <class Evaluation>
142  static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
143  { return eval_(params.SwPcwnSamples(), params.pcnwSamples(), Sw); }
144 
145  template <class Evaluation>
146  static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnw)
147  { return eval_(params.pcnwSamples(), params.SwPcwnSamples(), pcnw); }
148 
152  template <class FluidState, class Evaluation = typename FluidState::Scalar>
153  static Evaluation Sw(const Params& /* params */, const FluidState& /* fs */)
154  { throw std::logic_error("Not implemented: Sw()"); }
155 
156  template <class Evaluation>
157  static Evaluation twoPhaseSatSw(const Params& /* params */, const Evaluation& /* pC */)
158  { throw std::logic_error("Not implemented: twoPhaseSatSw()"); }
159 
164  template <class FluidState, class Evaluation = typename FluidState::Scalar>
165  static Evaluation Sn(const Params& params, const FluidState& fs)
166  { return 1 - Sw<FluidState, Scalar>(params, fs); }
167 
168  template <class Evaluation>
169  static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pC)
170  { return 1 - twoPhaseSatSw(params, pC); }
171 
176  template <class FluidState, class Evaluation = typename FluidState::Scalar>
177  static Evaluation krw(const Params& params, const FluidState& fs)
178  {
179  const auto& Sw =
180  decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
181 
182  return twoPhaseSatKrw(params, Sw);
183  }
184 
185  template <class Evaluation>
186  static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
187  { return eval_(params.SwKrwSamples(), params.krwSamples(), Sw); }
188 
189  template <class Evaluation>
190  static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
191  { return eval_(params.krwSamples(), params.SwKrwSamples(), krw); }
192 
197  template <class FluidState, class Evaluation = typename FluidState::Scalar>
198  static Evaluation krn(const Params& params, const FluidState& fs)
199  {
200  const auto& Sw =
201  decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
202 
203  return twoPhaseSatKrn(params, Sw);
204  }
205 
206  template <class Evaluation>
207  static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
208  { return eval_(params.SwKrnSamples(), params.krnSamples(), Sw); }
209 
210  template <class Evaluation>
211  static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
212  { return eval_(params.krnSamples(), params.SwKrnSamples(), krn); }
213 
214 private:
215  template <class Evaluation>
216  static Evaluation eval_(const ValueVector& xValues,
217  const ValueVector& yValues,
218  const Evaluation& x)
219  {
220  if (xValues.front() < xValues.back())
221  return evalAscending_(xValues, yValues, x);
222  return evalDescending_(xValues, yValues, x);
223  }
224 
225  template <class Evaluation>
226  static Evaluation evalAscending_(const ValueVector& xValues,
227  const ValueVector& yValues,
228  const Evaluation& x)
229  {
230  if (x <= xValues.front())
231  return yValues.front();
232  if (x >= xValues.back())
233  return yValues.back();
234 
235  size_t segIdx = findSegmentIndex_(xValues, scalarValue(x));
236 
237  Scalar x0 = xValues[segIdx];
238  Scalar x1 = xValues[segIdx + 1];
239 
240  Scalar y0 = yValues[segIdx];
241  Scalar y1 = yValues[segIdx + 1];
242 
243  Scalar m = (y1 - y0)/(x1 - x0);
244 
245  return y0 + (x - x0)*m;
246  }
247 
248  template <class Evaluation>
249  static Evaluation evalDescending_(const ValueVector& xValues,
250  const ValueVector& yValues,
251  const Evaluation& x)
252  {
253  if (x >= xValues.front())
254  return yValues.front();
255  if (x <= xValues.back())
256  return yValues.back();
257 
258  size_t segIdx = findSegmentIndexDescending_(xValues, scalarValue(x));
259 
260  Scalar x0 = xValues[segIdx];
261  Scalar x1 = xValues[segIdx + 1];
262 
263  Scalar y0 = yValues[segIdx];
264  Scalar y1 = yValues[segIdx + 1];
265 
266  Scalar m = (y1 - y0)/(x1 - x0);
267 
268  return y0 + (x - x0)*m;
269  }
270 
271  template <class Evaluation>
272  static Evaluation evalDeriv_(const ValueVector& xValues,
273  const ValueVector& yValues,
274  const Evaluation& x)
275  {
276  if (x <= xValues.front())
277  return 0.0;
278  if (x >= xValues.back())
279  return 0.0;
280 
281  size_t segIdx = findSegmentIndex_(xValues, scalarValue(x));
282 
283  Scalar x0 = xValues[segIdx];
284  Scalar x1 = xValues[segIdx + 1];
285 
286  Scalar y0 = yValues[segIdx];
287  Scalar y1 = yValues[segIdx + 1];
288 
289  return (y1 - y0)/(x1 - x0);
290  }
291 
292  static size_t findSegmentIndex_(const ValueVector& xValues, Scalar x)
293  {
294  assert(xValues.size() > 1); // we need at least two sampling points!
295  size_t n = xValues.size() - 1;
296  if (xValues.back() <= x)
297  return n - 1;
298  else if (x <= xValues.front())
299  return 0;
300 
301  // bisection
302  size_t lowIdx = 0, highIdx = n;
303  while (lowIdx + 1 < highIdx) {
304  size_t curIdx = (lowIdx + highIdx)/2;
305  if (xValues[curIdx] < x)
306  lowIdx = curIdx;
307  else
308  highIdx = curIdx;
309  }
310 
311  return lowIdx;
312  }
313 
314  static size_t findSegmentIndexDescending_(const ValueVector& xValues, Scalar x)
315  {
316  assert(xValues.size() > 1); // we need at least two sampling points!
317  size_t n = xValues.size() - 1;
318  if (x <= xValues.back())
319  return n;
320  else if (xValues.front() <= x)
321  return 0;
322 
323  // bisection
324  size_t lowIdx = 0, highIdx = n;
325  while (lowIdx + 1 < highIdx) {
326  size_t curIdx = (lowIdx + highIdx)/2;
327  if (xValues[curIdx] >= x)
328  lowIdx = curIdx;
329  else
330  highIdx = curIdx;
331  }
332 
333  return lowIdx;
334  }
335 };
336 } // namespace Opm
337 
338 #endif
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Specification of the material parameters for a two-phase material law which uses a table and piecewis...
Implementation of a tabulated, piecewise linear capillary pressure law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:51
static Evaluation krn(const Params &params, const FluidState &fs)
The relative permeability for the non-wetting phase of the porous medium.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:198
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation &Sw)
The saturation-capillary pressure curve.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:142
TraitsT Traits
The traits class for this material law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:56
static Evaluation krw(const Params &params, const FluidState &fs)
The relative permeability for the wetting phase of the porous medium.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:177
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:72
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
The relative permeabilities.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:118
static Evaluation Sn(const Params &params, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:165
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:88
static void saturations(Container &, const Params &, const FluidState &)
The saturations of the fluid phases starting from their pressure differences.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:111
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:92
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:80
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:76
static Evaluation pcnw(const Params &params, const FluidState &fs)
The capillary pressure-saturation curve.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:130
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:84
static const int numPhases
The number of fluid phases.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:65
Traits::Scalar Scalar
The type of the scalar values for this law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:62
ParamsT Params
The type of the parameter objects for this law.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:59
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
The capillary pressure-saturation curve.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:98
static Evaluation Sw(const Params &, const FluidState &)
The saturation-capillary pressure curve.
Definition: PiecewiseLinearTwoPhaseMaterial.hpp:153