My Project
Loading...
Searching...
No Matches
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#include <opm/common/TimingMacros.hpp>
33
34#include <stdexcept>
35#include <type_traits>
36
37
38#include <opm/common/utility/gpuDecorators.hpp>
39
40namespace Opm {
51template <class TraitsT, class ParamsT = PiecewiseLinearTwoPhaseMaterialParams<TraitsT> >
52class PiecewiseLinearTwoPhaseMaterial : public TraitsT
53{
54 using ValueVector = typename ParamsT::ValueVector;
55
56public:
58 using Traits = TraitsT;
59
61 using Params = ParamsT;
62
64 using Scalar = typename Traits::Scalar;
65
67 static constexpr 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 constexpr bool implementsTwoPhaseApi = true;
75
78 static constexpr bool implementsTwoPhaseSatApi = true;
79
82 static constexpr bool isSaturationDependent = true;
83
86 static constexpr bool isPressureDependent = false;
87
90 static constexpr bool isTemperatureDependent = false;
91
94 static constexpr bool isCompositionDependent = false;
95
99 template <class Container, class FluidState>
100 OPM_HOST_DEVICE static void capillaryPressures(Container& values, const Params& params, const FluidState& fs)
101 {
102 OPM_TIMEFUNCTION_LOCAL();
103 using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
104
105 values[Traits::wettingPhaseIdx] = 0.0; // reference phase
106 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
107 }
108
113 template <class Container, class FluidState>
114 OPM_HOST_DEVICE static void saturations(Container& /* values */, const Params& /* params */, const FluidState& /* fs */)
115 { throw std::logic_error("Not implemented: saturations()"); }
116
120 template <class Container, class FluidState>
121 OPM_HOST_DEVICE static void relativePermeabilities(Container& values, const Params& params, const FluidState& fs)
122 {
123 OPM_TIMEFUNCTION_LOCAL();
124 using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
125
126 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
127 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
128 }
129
133 template <class FluidState, class Evaluation = typename FluidState::Scalar>
134 OPM_HOST_DEVICE static Evaluation pcnw(const Params& params, const FluidState& fs)
135 {
136 OPM_TIMEFUNCTION_LOCAL();
137 const auto& Sw =
138 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
139
140 return twoPhaseSatPcnw(params, Sw);
141 }
142
146 template <class Evaluation>
147 OPM_HOST_DEVICE static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
148 {
149 OPM_TIMEFUNCTION_LOCAL();
150 return eval_(params.SwPcwnSamples(), params.pcwnSamples(), Sw);
151 }
152
153 template <class Evaluation>
154 static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnw)
155 { return eval_(params.pcwnSamples(), params.SwPcwnSamples(), pcnw); }
156
160 template <class FluidState, class Evaluation = typename FluidState::Scalar>
161 OPM_HOST_DEVICE static Evaluation Sw(const Params& /* params */, const FluidState& /* fs */)
162 { throw std::logic_error("Not implemented: Sw()"); }
163
164 template <class Evaluation>
165 OPM_HOST_DEVICE static Evaluation twoPhaseSatSw(const Params& /* params */, const Evaluation& /* pC */)
166 { throw std::logic_error("Not implemented: twoPhaseSatSw()"); }
167
172 template <class FluidState, class Evaluation = typename FluidState::Scalar>
173 OPM_HOST_DEVICE static Evaluation Sn(const Params& params, const FluidState& fs)
174 { return 1 - Sw<FluidState, Scalar>(params, fs); }
175
176 template <class Evaluation>
177 OPM_HOST_DEVICE static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pC)
178 { return 1 - twoPhaseSatSw(params, pC); }
179
184 template <class FluidState, class Evaluation = typename FluidState::Scalar>
185 OPM_HOST_DEVICE static Evaluation krw(const Params& params, const FluidState& fs)
186 {
187 OPM_TIMEFUNCTION_LOCAL();
188 const auto& Sw =
189 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
190
191 return twoPhaseSatKrw(params, Sw);
192 }
193
194 template <class Evaluation>
195 OPM_HOST_DEVICE static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
196 {
197 OPM_TIMEFUNCTION_LOCAL();
198 return eval_(params.SwKrwSamples(), params.krwSamples(), Sw);
199 }
200
201 template <class Evaluation>
202 OPM_HOST_DEVICE static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
203 {
204 OPM_TIMEFUNCTION_LOCAL();
205 return eval_(params.krwSamples(), params.SwKrwSamples(), krw);
206 }
207
212 template <class FluidState, class Evaluation = typename FluidState::Scalar>
213 OPM_HOST_DEVICE static Evaluation krn(const Params& params, const FluidState& fs)
214 {
215 OPM_TIMEFUNCTION_LOCAL();
216 const auto& Sw =
217 decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
218
219 return twoPhaseSatKrn(params, Sw);
220 }
221
222 template <class Evaluation>
223 OPM_HOST_DEVICE static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
224 {
225 OPM_TIMEFUNCTION_LOCAL();
226 return eval_(params.SwKrnSamples(), params.krnSamples(), Sw);
227 }
228
229 template <class Evaluation>
230 OPM_HOST_DEVICE static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
231 { return eval_(params.krnSamples(), params.SwKrnSamples(), krn); }
232
233 template <class Evaluation>
234 OPM_HOST_DEVICE static size_t findSegmentIndex(const ValueVector& xValues, const Evaluation& x){
235 return findSegmentIndex_(xValues, scalarValue(x));
236 }
237
238 template <class Evaluation>
239 OPM_HOST_DEVICE static size_t findSegmentIndexDescending(const ValueVector& xValues, const Evaluation& x){
240 return findSegmentIndexDescending_(xValues, scalarValue(x));
241 }
242
243 template <class Evaluation>
244 OPM_HOST_DEVICE static Evaluation eval(const ValueVector& xValues, const ValueVector& yValues, const Evaluation& x, unsigned segIdx){
245 Scalar x0 = xValues[segIdx];
246 Scalar x1 = xValues[segIdx + 1];
247
248 Scalar y0 = yValues[segIdx];
249 Scalar y1 = yValues[segIdx + 1];
250
251 Scalar m = (y1 - y0)/(x1 - x0);
252
253 return y0 + (x - x0)*m;
254 }
255
256private:
257 template <class Evaluation>
258 OPM_HOST_DEVICE static Evaluation eval_(const ValueVector& xValues,
259 const ValueVector& yValues,
260 const Evaluation& x)
261 {
262 OPM_TIMEFUNCTION_LOCAL();
263 if (xValues.front() < xValues.back())
264 return evalAscending_(xValues, yValues, x);
265 return evalDescending_(xValues, yValues, x);
266 }
267
268 template <class Evaluation>
269 OPM_HOST_DEVICE static Evaluation evalAscending_(const ValueVector& xValues,
270 const ValueVector& yValues,
271 const Evaluation& x)
272 {
273 OPM_TIMEFUNCTION_LOCAL();
274 if (x <= xValues.front())
275 return yValues.front();
276 if (x >= xValues.back())
277 return yValues.back();
278
279 size_t segIdx = findSegmentIndex_(xValues, scalarValue(x));
280
281 return eval(xValues, yValues, x, segIdx);
282 }
283
284 template <class Evaluation>
285 OPM_HOST_DEVICE static Evaluation evalDescending_(const ValueVector& xValues,
286 const ValueVector& yValues,
287 const Evaluation& x)
288 {
289 OPM_TIMEFUNCTION_LOCAL();
290 if (x >= xValues.front())
291 return yValues.front();
292 if (x <= xValues.back())
293 return yValues.back();
294
295 size_t segIdx = findSegmentIndexDescending_(xValues, scalarValue(x));
296
297 return eval(xValues, yValues, x, segIdx);
298 }
299
300 template <class Evaluation>
301 OPM_HOST_DEVICE static Evaluation evalDeriv_(const ValueVector& xValues,
302 const ValueVector& yValues,
303 const Evaluation& x)
304 {
305 OPM_TIMEFUNCTION_LOCAL();
306 if (x <= xValues.front())
307 return 0.0;
308 if (x >= xValues.back())
309 return 0.0;
310
311 size_t segIdx = findSegmentIndex_(xValues, scalarValue(x));
312
313 Scalar x0 = xValues[segIdx];
314 Scalar x1 = xValues[segIdx + 1];
315
316 Scalar y0 = yValues[segIdx];
317 Scalar y1 = yValues[segIdx + 1];
318
319 return (y1 - y0)/(x1 - x0);
320 }
321 template<class ScalarT>
322 OPM_HOST_DEVICE static size_t findSegmentIndex_(const ValueVector& xValues, const ScalarT& x)
323 {
324 OPM_TIMEFUNCTION_LOCAL();
325 assert(xValues.size() > 1); // we need at least two sampling points!
326 size_t n = xValues.size() - 1;
327 if (xValues.back() <= x)
328 return n - 1;
329 else if (x <= xValues.front())
330 return 0;
331
332 // bisection
333 size_t lowIdx = 0, highIdx = n;
334 while (lowIdx + 1 < highIdx) {
335 size_t curIdx = (lowIdx + highIdx)/2;
336 if (xValues[curIdx] < x)
337 lowIdx = curIdx;
338 else
339 highIdx = curIdx;
340 }
341
342 return lowIdx;
343 }
344
345 OPM_HOST_DEVICE static size_t findSegmentIndexDescending_(const ValueVector& xValues, Scalar x)
346 {
347 OPM_TIMEFUNCTION_LOCAL();
348 assert(xValues.size() > 1); // we need at least two sampling points!
349 size_t n = xValues.size() - 1;
350 if (x <= xValues.back())
351 return n;
352 else if (xValues.front() <= x)
353 return 0;
354
355 // bisection
356 size_t lowIdx = 0, highIdx = n;
357 while (lowIdx + 1 < highIdx) {
358 size_t curIdx = (lowIdx + highIdx)/2;
359 if (xValues[curIdx] >= x)
360 lowIdx = curIdx;
361 else
362 highIdx = curIdx;
363 }
364
365 return lowIdx;
366 }
367};
368
369} // namespace Opm
370
371#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:53
TraitsT Traits
The traits class for this material law.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:58
ParamsT Params
The type of the parameter objects for this law.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:61
static OPM_HOST_DEVICE void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
The capillary pressure-saturation curve.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:100
static OPM_HOST_DEVICE Evaluation Sn(const Params &params, const FluidState &fs)
Calculate the non-wetting phase saturations depending on the phase pressures.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:173
static constexpr int numPhases
The number of fluid phases.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:67
static OPM_HOST_DEVICE void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
The relative permeabilities.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:121
static constexpr bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition PiecewiseLinearTwoPhaseMaterial.hpp:78
static OPM_HOST_DEVICE Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation &Sw)
The saturation-capillary pressure curve.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:147
typename Traits::Scalar Scalar
The type of the scalar values for this law.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:64
static OPM_HOST_DEVICE Evaluation krw(const Params &params, const FluidState &fs)
The relative permeability for the wetting phase of the porous medium.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:185
static constexpr bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:94
static constexpr bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:82
static OPM_HOST_DEVICE void saturations(Container &, const Params &, const FluidState &)
The saturations of the fluid phases starting from their pressure differences.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:114
static OPM_HOST_DEVICE Evaluation krn(const Params &params, const FluidState &fs)
The relative permeability for the non-wetting phase of the porous medium.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:213
static constexpr bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:90
static constexpr bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:86
static OPM_HOST_DEVICE Evaluation Sw(const Params &, const FluidState &)
The saturation-capillary pressure curve.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:161
static constexpr bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:74
static OPM_HOST_DEVICE Evaluation pcnw(const Params &params, const FluidState &fs)
The capillary pressure-saturation curve.
Definition PiecewiseLinearTwoPhaseMaterial.hpp:134
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30