My Project
Loading...
Searching...
No Matches
EclDefaultMaterial.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_ECL_DEFAULT_MATERIAL_HPP
28#define OPM_ECL_DEFAULT_MATERIAL_HPP
29
31
34
35#include <algorithm>
36#include <stdexcept>
37#include <type_traits>
38
39namespace Opm {
40
54template <class TraitsT,
55 class GasOilMaterialLawT,
56 class OilWaterMaterialLawT,
57 class ParamsT = EclDefaultMaterialParams<TraitsT,
58 typename GasOilMaterialLawT::Params,
59 typename OilWaterMaterialLawT::Params> >
60class EclDefaultMaterial : public TraitsT
61{
62public:
63 using GasOilMaterialLaw = GasOilMaterialLawT;
64 using OilWaterMaterialLaw = OilWaterMaterialLawT;
65
66 // some safety checks
67 static_assert(TraitsT::numPhases == 3,
68 "The number of phases considered by this capillary pressure "
69 "law is always three!");
70 static_assert(GasOilMaterialLaw::numPhases == 2,
71 "The number of phases considered by the gas-oil capillary "
72 "pressure law must be two!");
73 static_assert(OilWaterMaterialLaw::numPhases == 2,
74 "The number of phases considered by the oil-water capillary "
75 "pressure law must be two!");
76 static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
77 typename OilWaterMaterialLaw::Scalar>::value,
78 "The two two-phase capillary pressure laws must use the same "
79 "type of floating point values.");
80
81 static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
82 "The gas-oil material law must implement the two-phase saturation "
83 "only API to for the default Ecl capillary pressure law!");
84 static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
85 "The oil-water material law must implement the two-phase saturation "
86 "only API to for the default Ecl capillary pressure law!");
87
88 using Traits = TraitsT;
89 using Params = ParamsT;
90 using Scalar = typename Traits::Scalar;
91
92 static constexpr int numPhases = 3;
93 static constexpr int waterPhaseIdx = Traits::wettingPhaseIdx;
94 static constexpr int oilPhaseIdx = Traits::nonWettingPhaseIdx;
95 static constexpr int gasPhaseIdx = Traits::gasPhaseIdx;
96
99 static constexpr bool implementsTwoPhaseApi = false;
100
103 static constexpr bool implementsTwoPhaseSatApi = false;
104
107 static constexpr bool isSaturationDependent = true;
108
111 static constexpr bool isPressureDependent = false;
112
115 static constexpr bool isTemperatureDependent = false;
116
119 static constexpr bool isCompositionDependent = false;
120
135 template <class ContainerT, class FluidState>
136 static void capillaryPressures(ContainerT& values,
137 const Params& params,
138 const FluidState& state)
139 {
140 OPM_TIMEFUNCTION_LOCAL();
141 using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
142 values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
143 values[oilPhaseIdx] = 0;
144 values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
145
146 Valgrind::CheckDefined(values[gasPhaseIdx]);
147 Valgrind::CheckDefined(values[oilPhaseIdx]);
148 Valgrind::CheckDefined(values[waterPhaseIdx]);
149 }
150
151 /*
152 * Hysteresis parameters for oil-water
153 * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
154 * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
155 * \param params Parameters
156 */
157 static void oilWaterHysteresisParams(Scalar& pcSwMdc,
158 Scalar& krnSwMdc,
159 const Params& params)
160 {
161 pcSwMdc = params.oilWaterParams().pcSwMdc();
162 krnSwMdc = params.oilWaterParams().krnSwMdc();
163
164 Valgrind::CheckDefined(pcSwMdc);
165 Valgrind::CheckDefined(krnSwMdc);
166 }
167
168 /*
169 * Hysteresis parameters for oil-water
170 * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
171 * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
172 * \param params Parameters
173 */
174 static void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
175 const Scalar& krnSwMdc,
176 Params& params)
177 {
178 constexpr const double krwSw = 2.0; //Should not be used
179 params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc);
180 }
181
182 /*
183 * Hysteresis parameters for gas-oil
184 * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
185 * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
186 * \param params Parameters
187 */
188 static void gasOilHysteresisParams(Scalar& pcSwMdc,
189 Scalar& krnSwMdc,
190 const Params& params)
191 {
192 const auto Swco = params.Swl();
193
194 // Pretend oil saturation is 'Swco' larger than it really is in
195 // order to infer correct maximum Sg values in output layer.
196 pcSwMdc = std::min(params.gasOilParams().pcSwMdc() + Swco, Scalar{2.0});
197 krnSwMdc = std::min(params.gasOilParams().krnSwMdc() + Swco, Scalar{2.0});
198
199 Valgrind::CheckDefined(pcSwMdc);
200 Valgrind::CheckDefined(krnSwMdc);
201 }
202
203
204 static Scalar trappedGasSaturation(const Params& params)
205 {
206
207 return params.gasOilParams().SnTrapped();
208 }
209 /*
210 * Hysteresis parameters for gas-oil
211 * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
212 * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
213 * \param params Parameters
214 */
215 static void setGasOilHysteresisParams(const Scalar& pcSwMdc,
216 const Scalar& krnSwMdc,
217 Params& params)
218 {
219 // Maximum attainable oil saturation is 1-SWL
220 const auto Swco = params.Swl();
221 constexpr const double krwSw = 2.0; //Should not be used
222 params.gasOilParams().update(pcSwMdc - Swco, krwSw, krnSwMdc - Swco);
223 }
224
234 template <class FluidState, class Evaluation = typename FluidState::Scalar>
235 static Evaluation pcgn(const Params& params,
236 const FluidState& fs)
237 {
238 OPM_TIMEFUNCTION_LOCAL();
239 // Maximum attainable oil saturation is 1-SWL.
240 const auto Sw = 1.0 - params.Swl() - decay<Evaluation>(fs.saturation(gasPhaseIdx));
241 return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
242 }
243
253 template <class FluidState, class Evaluation = typename FluidState::Scalar>
254 static Evaluation pcnw(const Params& params,
255 const FluidState& fs)
256 {
257 OPM_TIMEFUNCTION_LOCAL();
258 const auto Sw = decay<Evaluation>(fs.saturation(waterPhaseIdx));
259 return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
260 }
261
265 template <class ContainerT, class FluidState>
266 static void saturations(ContainerT& /*values*/,
267 const Params& /*params*/,
268 const FluidState& /*fluidState*/)
269 {
270 throw std::logic_error("Not implemented: saturations()");
271 }
272
276 template <class FluidState, class Evaluation = typename FluidState::Scalar>
277 static Evaluation Sg(const Params& /*params*/,
278 const FluidState& /*fluidState*/)
279 {
280 throw std::logic_error("Not implemented: Sg()");
281 }
282
286 template <class FluidState, class Evaluation = typename FluidState::Scalar>
287 static Evaluation Sn(const Params& /*params*/,
288 const FluidState& /*fluidState*/)
289 {
290 throw std::logic_error("Not implemented: Sn()");
291 }
292
296 template <class FluidState, class Evaluation = typename FluidState::Scalar>
297 static Evaluation Sw(const Params& /*params*/,
298 const FluidState& /*fluidState*/)
299 {
300 throw std::logic_error("Not implemented: Sw()");
301 }
302
318 template <class ContainerT, class FluidState>
319 static void relativePermeabilities(ContainerT& values,
320 const Params& params,
321 const FluidState& fluidState)
322 {
323 OPM_TIMEFUNCTION_LOCAL();
324 using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
325
326 values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
327 values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
328 values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
329 }
330
334 template <class FluidState, class Evaluation = typename FluidState::Scalar>
335 static Evaluation krg(const Params& params,
336 const FluidState& fluidState)
337 {
338 OPM_TIMEFUNCTION_LOCAL();
339 // Maximum attainable oil saturation is 1-SWL.
340 const Evaluation Sw = 1.0 - params.Swl() - decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
341 return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
342 }
343
347 template <class FluidState, class Evaluation = typename FluidState::Scalar>
348 static Evaluation krw(const Params& params,
349 const FluidState& fluidState)
350 {
351 OPM_TIMEFUNCTION_LOCAL();
352 const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
353 return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
354 }
355
359 template <class FluidState, class Evaluation = typename FluidState::Scalar>
360 static Evaluation krn(const Params& params,
361 const FluidState& fluidState)
362 {
363 OPM_TIMEFUNCTION_LOCAL();
364 const Scalar Swco = params.Swl();
365
366 const Evaluation Sw =
367 max(Evaluation(Swco),
368 decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
369
370 const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
371
372 const Evaluation Sw_ow = Sg + Sw;
373 const Evaluation kro_ow = relpermOilInOilWaterSystem<Evaluation>(params, fluidState);
374 const Evaluation kro_go = relpermOilInOilGasSystem<Evaluation>(params, fluidState);
375
376 // avoid the division by zero: chose a regularized kro which is used if Sw - Swco
377 // < epsilon/2 and interpolate between the oridinary and the regularized kro between
378 // epsilon and epsilon/2
379 constexpr const Scalar epsilon = 1e-5;
380 if (scalarValue(Sw_ow) - Swco < epsilon) {
381 const Evaluation kro2 = (kro_ow + kro_go)/2;
382 if (scalarValue(Sw_ow) - Swco > epsilon/2) {
383 const Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
384 const Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
385
386 return kro2*alpha + kro1*(1 - alpha);
387 }
388
389 return kro2;
390 }
391
392 return (Sg*kro_go + (Sw - Swco)*kro_ow) / (Sw_ow - Swco);
393 }
394
398 template <class Evaluation, class FluidState>
399 static Evaluation relpermOilInOilGasSystem(const Params& params,
400 const FluidState& fluidState)
401 {
402 OPM_TIMEFUNCTION_LOCAL();
403 const Evaluation Sw =
404 max(Evaluation{ params.Swl() },
405 decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
406
407 const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
408 const Evaluation So_go = 1.0 - (Sg + Sw);
409
410 return GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
411 }
412
416 template <class Evaluation, class FluidState>
417 static Evaluation relpermOilInOilWaterSystem(const Params& params,
418 const FluidState& fluidState)
419 {
420 OPM_TIMEFUNCTION_LOCAL();
421 const Evaluation Sw =
422 max(Evaluation{ params.Swl() },
423 decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
424
425 const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
426 const Evaluation Sw_ow = Sg + Sw;
427
428 return OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
429 }
430
438 template <class FluidState>
439 static bool updateHysteresis(Params& params, const FluidState& fluidState)
440 {
441 OPM_TIMEFUNCTION_LOCAL();
442 bool changed = false;
443 const Scalar Swco = params.Swl();
444
445 const Scalar Sw = clampSaturation(fluidState, waterPhaseIdx);
446 const Scalar So = clampSaturation(fluidState, oilPhaseIdx);
447 const Scalar Sg = clampSaturation(fluidState, gasPhaseIdx);
448
449 if (params.inconsistentHysteresisUpdate()) {
450 // NOTE: the saturations which are passed to update the hysteresis curves are
451 // inconsistent with the ones used to calculate the relative permabilities. We do
452 // it like this anyway because (a) the saturation functions of opm-core do it
453 // this way (b) the simulations seem to converge better (which is not too much
454 // surprising actually, because the time step does not start on a kink in the
455 // solution) and (c) the Eclipse 100 simulator may do the same.
456 //
457 // Though be aware that from a physical perspective this is definitively
458 // incorrect!
459 bool oilchanged = params.oilWaterParams().update(/*pcSw=*/ Sw, //1.0 - So, (Effect is significant vs benchmark.)
460 /*krwSw=*/ 1.0 - So,
461 /*krnSw=*/ 1.0 - So);
462
463 changed = changed || oilchanged;
464
465 bool gaschanged = params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
466 /*krwSw=*/ 1.0 - std::max(Swco, Sw) - Sg, //(Three phase behavior ...)
467 /*krnSw=*/ 1.0 - Swco - Sg);
468 changed = changed || gaschanged;
469 }
470 else {
471 const Scalar Sw_ow = Sg + std::max(Swco, Sw);
472 const Scalar So_go = 1.0 - Sw_ow;
473 bool oilchanged = params.oilWaterParams().update(/*pcSw=*/ Sw,
474 /*krwSw=*/ 1 - Sg,
475 /*krnSw=*/ Sw_ow);
476 changed = changed || oilchanged;
477 bool gaschanged = params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
478 /*krwSw=*/ So_go,
479 /*krnSw=*/ 1.0 - Swco - Sg);
480
481 changed = changed || gaschanged;
482 }
483 return changed;
484 }
485
486 template <class FluidState>
487 static Scalar clampSaturation(const FluidState& fluidState, const int phaseIndex)
488 {
489 OPM_TIMEFUNCTION_LOCAL();
490 const auto sat = scalarValue(fluidState.saturation(phaseIndex));
491 return std::clamp(sat, Scalar{0.0}, Scalar{1.0});
492 }
493};
494} // namespace Opm
495
496#endif
Default implementation for the parameters required by the default three-phase capillary pressure mode...
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Some templates to wrap the valgrind client request macros.
Implements the default three phase capillary pressure law used by the ECLipse simulator.
Definition EclDefaultMaterial.hpp:61
static Evaluation relpermOilInOilWaterSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/water system.
Definition EclDefaultMaterial.hpp:417
static Evaluation krg(const Params &params, const FluidState &fluidState)
The relative permeability of the gas phase.
Definition EclDefaultMaterial.hpp:335
static Evaluation pcgn(const Params &params, const FluidState &fs)
Capillary pressure between the gas and the non-wetting liquid (i.e., oil) phase.
Definition EclDefaultMaterial.hpp:235
static constexpr bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition EclDefaultMaterial.hpp:119
static constexpr bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition EclDefaultMaterial.hpp:115
static Evaluation Sg(const Params &, const FluidState &)
The saturation of the gas phase.
Definition EclDefaultMaterial.hpp:277
static Evaluation relpermOilInOilGasSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/gas system.
Definition EclDefaultMaterial.hpp:399
static constexpr bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition EclDefaultMaterial.hpp:99
static void capillaryPressures(ContainerT &values, const Params &params, const FluidState &state)
Implements the default three phase capillary pressure law used by the ECLipse simulator.
Definition EclDefaultMaterial.hpp:136
static void saturations(ContainerT &, const Params &, const FluidState &)
The inverse of the capillary pressure.
Definition EclDefaultMaterial.hpp:266
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition EclDefaultMaterial.hpp:287
static Evaluation krn(const Params &params, const FluidState &fluidState)
The relative permeability of the non-wetting (i.e., oil) phase.
Definition EclDefaultMaterial.hpp:360
static Evaluation pcnw(const Params &params, const FluidState &fs)
Capillary pressure between the non-wetting liquid (i.e., oil) and the wetting liquid (i....
Definition EclDefaultMaterial.hpp:254
static constexpr bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition EclDefaultMaterial.hpp:107
static bool updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition EclDefaultMaterial.hpp:439
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition EclDefaultMaterial.hpp:319
static Evaluation krw(const Params &params, const FluidState &fluidState)
The relative permeability of the wetting phase.
Definition EclDefaultMaterial.hpp:348
static Evaluation Sw(const Params &, const FluidState &)
The saturation of the wetting (i.e., water) phase.
Definition EclDefaultMaterial.hpp:297
static constexpr bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition EclDefaultMaterial.hpp:111
static constexpr bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition EclDefaultMaterial.hpp:103
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30