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