My Project
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 
39 namespace Opm {
40 
54 template <class TraitsT,
55  class GasOilMaterialLawT,
56  class OilWaterMaterialLawT,
57  class ParamsT = EclDefaultMaterialParams<TraitsT,
58  typename GasOilMaterialLawT::Params,
59  typename OilWaterMaterialLawT::Params> >
60 class EclDefaultMaterial : public TraitsT
61 {
62 public:
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  using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
141  values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
142  values[oilPhaseIdx] = 0;
143  values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
144 
145  Valgrind::CheckDefined(values[gasPhaseIdx]);
146  Valgrind::CheckDefined(values[oilPhaseIdx]);
147  Valgrind::CheckDefined(values[waterPhaseIdx]);
148  }
149 
150  /*
151  * Hysteresis parameters for oil-water
152  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
153  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
154  * \param params Parameters
155  */
156  static void oilWaterHysteresisParams(Scalar& pcSwMdc,
157  Scalar& krnSwMdc,
158  const Params& params)
159  {
160  pcSwMdc = params.oilWaterParams().pcSwMdc();
161  krnSwMdc = params.oilWaterParams().krnSwMdc();
162 
163  Valgrind::CheckDefined(pcSwMdc);
164  Valgrind::CheckDefined(krnSwMdc);
165  }
166 
167  /*
168  * Hysteresis parameters for oil-water
169  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
170  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
171  * \param params Parameters
172  */
173  static void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
174  const Scalar& krnSwMdc,
175  Params& params)
176  {
177  constexpr const double krwSw = 2.0; //Should not be used
178  params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc);
179  }
180 
181  /*
182  * Hysteresis parameters for gas-oil
183  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
184  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
185  * \param params Parameters
186  */
187  static void gasOilHysteresisParams(Scalar& pcSwMdc,
188  Scalar& krnSwMdc,
189  const Params& params)
190  {
191  const auto Swco = params.Swl();
192 
193  // Pretend oil saturation is 'Swco' larger than it really is in
194  // order to infer correct maximum Sg values in output layer.
195  pcSwMdc = std::min(params.gasOilParams().pcSwMdc() + Swco, Scalar{2.0});
196  krnSwMdc = std::min(params.gasOilParams().krnSwMdc() + Swco, Scalar{2.0});
197 
198  Valgrind::CheckDefined(pcSwMdc);
199  Valgrind::CheckDefined(krnSwMdc);
200  }
201 
202  /*
203  * Hysteresis parameters for gas-oil
204  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
205  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
206  * \param params Parameters
207  */
208  static void setGasOilHysteresisParams(const Scalar& pcSwMdc,
209  const Scalar& krnSwMdc,
210  Params& params)
211  {
212  // Maximum attainable oil saturation is 1-SWL
213  const auto Swco = params.Swl();
214  constexpr const double krwSw = 2.0; //Should not be used
215  params.gasOilParams().update(pcSwMdc - Swco, krwSw, krnSwMdc - Swco);
216  }
217 
227  template <class FluidState, class Evaluation = typename FluidState::Scalar>
228  static Evaluation pcgn(const Params& params,
229  const FluidState& fs)
230  {
231  // Maximum attainable oil saturation is 1-SWL.
232  const auto Sw = 1.0 - params.Swl() - decay<Evaluation>(fs.saturation(gasPhaseIdx));
233  return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
234  }
235 
245  template <class FluidState, class Evaluation = typename FluidState::Scalar>
246  static Evaluation pcnw(const Params& params,
247  const FluidState& fs)
248  {
249  const auto Sw = decay<Evaluation>(fs.saturation(waterPhaseIdx));
250  return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
251  }
252 
256  template <class ContainerT, class FluidState>
257  static void saturations(ContainerT& /*values*/,
258  const Params& /*params*/,
259  const FluidState& /*fluidState*/)
260  {
261  throw std::logic_error("Not implemented: saturations()");
262  }
263 
267  template <class FluidState, class Evaluation = typename FluidState::Scalar>
268  static Evaluation Sg(const Params& /*params*/,
269  const FluidState& /*fluidState*/)
270  {
271  throw std::logic_error("Not implemented: Sg()");
272  }
273 
277  template <class FluidState, class Evaluation = typename FluidState::Scalar>
278  static Evaluation Sn(const Params& /*params*/,
279  const FluidState& /*fluidState*/)
280  {
281  throw std::logic_error("Not implemented: Sn()");
282  }
283 
287  template <class FluidState, class Evaluation = typename FluidState::Scalar>
288  static Evaluation Sw(const Params& /*params*/,
289  const FluidState& /*fluidState*/)
290  {
291  throw std::logic_error("Not implemented: Sw()");
292  }
293 
309  template <class ContainerT, class FluidState>
310  static void relativePermeabilities(ContainerT& values,
311  const Params& params,
312  const FluidState& fluidState)
313  {
314  using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
315 
316  values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
317  values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
318  values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
319  }
320 
324  template <class FluidState, class Evaluation = typename FluidState::Scalar>
325  static Evaluation krg(const Params& params,
326  const FluidState& fluidState)
327  {
328  // Maximum attainable oil saturation is 1-SWL.
329  const Evaluation Sw = 1.0 - params.Swl() - decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
330  return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
331  }
332 
336  template <class FluidState, class Evaluation = typename FluidState::Scalar>
337  static Evaluation krw(const Params& params,
338  const FluidState& fluidState)
339  {
340  const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
341  return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
342  }
343 
347  template <class FluidState, class Evaluation = typename FluidState::Scalar>
348  static Evaluation krn(const Params& params,
349  const FluidState& fluidState)
350  {
351  const Scalar Swco = params.Swl();
352 
353  const Evaluation Sw =
354  max(Evaluation(Swco),
355  decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
356 
357  const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
358 
359  const Evaluation Sw_ow = Sg + Sw;
360  const Evaluation kro_ow = relpermOilInOilWaterSystem<Evaluation>(params, fluidState);
361  const Evaluation kro_go = relpermOilInOilGasSystem<Evaluation>(params, fluidState);
362 
363  // avoid the division by zero: chose a regularized kro which is used if Sw - Swco
364  // < epsilon/2 and interpolate between the oridinary and the regularized kro between
365  // epsilon and epsilon/2
366  constexpr const Scalar epsilon = 1e-5;
367  if (scalarValue(Sw_ow) - Swco < epsilon) {
368  const Evaluation kro2 = (kro_ow + kro_go)/2;
369  if (scalarValue(Sw_ow) - Swco > epsilon/2) {
370  const Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
371  const Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
372 
373  return kro2*alpha + kro1*(1 - alpha);
374  }
375 
376  return kro2;
377  }
378 
379  return (Sg*kro_go + (Sw - Swco)*kro_ow) / (Sw_ow - Swco);
380  }
381 
385  template <class Evaluation, class FluidState>
386  static Evaluation relpermOilInOilGasSystem(const Params& params,
387  const FluidState& fluidState)
388  {
389  const Evaluation Sw =
390  max(Evaluation{ params.Swl() },
391  decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
392 
393  const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
394  const Evaluation So_go = 1.0 - (Sg + Sw);
395 
396  return GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
397  }
398 
402  template <class Evaluation, class FluidState>
403  static Evaluation relpermOilInOilWaterSystem(const Params& params,
404  const FluidState& fluidState)
405  {
406  const Evaluation Sw =
407  max(Evaluation{ params.Swl() },
408  decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
409 
410  const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
411  const Evaluation Sw_ow = Sg + Sw;
412 
413  return OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
414  }
415 
423  template <class FluidState>
424  static void updateHysteresis(Params& params, const FluidState& fluidState)
425  {
426  const Scalar Swco = params.Swl();
427 
428  const Scalar Sw = clampSaturation(fluidState, waterPhaseIdx);
429  const Scalar So = clampSaturation(fluidState, oilPhaseIdx);
430  const Scalar Sg = clampSaturation(fluidState, gasPhaseIdx);
431 
432  if (params.inconsistentHysteresisUpdate()) {
433  // NOTE: the saturations which are passed to update the hysteresis curves are
434  // inconsistent with the ones used to calculate the relative permabilities. We do
435  // it like this anyway because (a) the saturation functions of opm-core do it
436  // this way (b) the simulations seem to converge better (which is not too much
437  // surprising actually, because the time step does not start on a kink in the
438  // solution) and (c) the Eclipse 100 simulator may do the same.
439  //
440  // Though be aware that from a physical perspective this is definitively
441  // incorrect!
442  params.oilWaterParams().update(/*pcSw=*/ Sw, //1.0 - So, (Effect is significant vs benchmark.)
443  /*krwSw=*/ 1.0 - So,
444  /*krnSw=*/ 1.0 - So);
445 
446  params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
447  /*krwSw=*/ 1.0 - Swco - Sg,
448  /*krnSw=*/ 1.0 - Swco - Sg);
449  }
450  else {
451  const Scalar Sw_ow = Sg + std::max(Swco, Sw);
452  const Scalar So_go = 1.0 - Sw_ow;
453 
454  params.oilWaterParams().update(/*pcSw=*/ Sw,
455  /*krwSw=*/ 1 - Sg,
456  /*krnSw=*/ Sw_ow);
457 
458  params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
459  /*krwSw=*/ So_go,
460  /*krnSw=*/ 1.0 - Swco - Sg);
461  }
462  }
463 
464  template <class FluidState>
465  static Scalar clampSaturation(const FluidState& fluidState, const int phaseIndex)
466  {
467  const auto sat = scalarValue(fluidState.saturation(phaseIndex));
468  return std::clamp(sat, Scalar{0.0}, Scalar{1.0});
469  }
470 };
471 } // namespace Opm
472 
473 #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:403
static Evaluation krg(const Params &params, const FluidState &fluidState)
The relative permeability of the gas phase.
Definition: EclDefaultMaterial.hpp:325
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:228
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:268
static Evaluation relpermOilInOilGasSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/gas system.
Definition: EclDefaultMaterial.hpp:386
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:257
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition: EclDefaultMaterial.hpp:278
static Evaluation krn(const Params &params, const FluidState &fluidState)
The relative permeability of the non-wetting (i.e., oil) phase.
Definition: EclDefaultMaterial.hpp:348
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:246
static constexpr bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: EclDefaultMaterial.hpp:107
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition: EclDefaultMaterial.hpp:310
static Evaluation krw(const Params &params, const FluidState &fluidState)
The relative permeability of the wetting phase.
Definition: EclDefaultMaterial.hpp:337
static Evaluation Sw(const Params &, const FluidState &)
The saturation of the wetting (i.e., water) phase.
Definition: EclDefaultMaterial.hpp:288
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
static void updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition: EclDefaultMaterial.hpp:424