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 
37 namespace Opm {
38 
52 template <class TraitsT,
53  class GasOilMaterialLawT,
54  class OilWaterMaterialLawT,
55  class ParamsT = EclDefaultMaterialParams<TraitsT,
56  typename GasOilMaterialLawT::Params,
57  typename OilWaterMaterialLawT::Params> >
58 class EclDefaultMaterial : public TraitsT
59 {
60 public:
61  typedef GasOilMaterialLawT GasOilMaterialLaw;
62  typedef OilWaterMaterialLawT OilWaterMaterialLaw;
63 
64  // some safety checks
65  static_assert(TraitsT::numPhases == 3,
66  "The number of phases considered by this capillary pressure "
67  "law is always three!");
68  static_assert(GasOilMaterialLaw::numPhases == 2,
69  "The number of phases considered by the gas-oil capillary "
70  "pressure law must be two!");
71  static_assert(OilWaterMaterialLaw::numPhases == 2,
72  "The number of phases considered by the oil-water capillary "
73  "pressure law must be two!");
74  static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
75  typename OilWaterMaterialLaw::Scalar>::value,
76  "The two two-phase capillary pressure laws must use the same "
77  "type of floating point values.");
78 
79  static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
80  "The gas-oil material law must implement the two-phase saturation "
81  "only API to for the default Ecl capillary pressure law!");
82  static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
83  "The oil-water material law must implement the two-phase saturation "
84  "only API to for the default Ecl capillary pressure law!");
85 
86  typedef TraitsT Traits;
87  typedef ParamsT Params;
88  typedef typename Traits::Scalar Scalar;
89 
90  static const int numPhases = 3;
91  static const int waterPhaseIdx = Traits::wettingPhaseIdx;
92  static const int oilPhaseIdx = Traits::nonWettingPhaseIdx;
93  static const int gasPhaseIdx = Traits::gasPhaseIdx;
94 
97  static const bool implementsTwoPhaseApi = false;
98 
101  static const bool implementsTwoPhaseSatApi = false;
102 
105  static const bool isSaturationDependent = true;
106 
109  static const bool isPressureDependent = false;
110 
113  static const bool isTemperatureDependent = false;
114 
117  static const bool isCompositionDependent = false;
118 
133  template <class ContainerT, class FluidState>
134  static void capillaryPressures(ContainerT& values,
135  const Params& params,
136  const FluidState& state)
137  {
138  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
139  values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
140  values[oilPhaseIdx] = 0;
141  values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
142 
143  Valgrind::CheckDefined(values[gasPhaseIdx]);
144  Valgrind::CheckDefined(values[oilPhaseIdx]);
145  Valgrind::CheckDefined(values[waterPhaseIdx]);
146  }
147 
148  /*
149  * Hysteresis parameters for oil-water
150  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
151  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
152  * \param params Parameters
153  */
154  static void oilWaterHysteresisParams(Scalar& pcSwMdc,
155  Scalar& krnSwMdc,
156  const Params& params)
157  {
158  pcSwMdc = params.oilWaterParams().pcSwMdc();
159  krnSwMdc = params.oilWaterParams().krnSwMdc();
160 
161  Valgrind::CheckDefined(pcSwMdc);
162  Valgrind::CheckDefined(krnSwMdc);
163  }
164 
165  /*
166  * Hysteresis parameters for oil-water
167  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
168  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
169  * \param params Parameters
170  */
171  static void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
172  const Scalar& krnSwMdc,
173  Params& params)
174  {
175  const double krwSw = 2.0; //Should not be used
176  params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc);
177  }
178 
179  /*
180  * Hysteresis parameters for gas-oil
181  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
182  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
183  * \param params Parameters
184  */
185  static void gasOilHysteresisParams(Scalar& pcSwMdc,
186  Scalar& krnSwMdc,
187  const Params& params)
188  {
189  const auto Swco = params.Swl();
190 
191  // Pretend oil saturation is 'Swco' larger than it really is in
192  // order to infer correct maximum Sg values in output layer.
193  pcSwMdc = std::min(params.gasOilParams().pcSwMdc() + Swco, Scalar{2.0});
194  krnSwMdc = std::min(params.gasOilParams().krnSwMdc() + Swco, Scalar{2.0});
195 
196  Valgrind::CheckDefined(pcSwMdc);
197  Valgrind::CheckDefined(krnSwMdc);
198  }
199 
200  /*
201  * Hysteresis parameters for gas-oil
202  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
203  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
204  * \param params Parameters
205  */
206  static void setGasOilHysteresisParams(const Scalar& pcSwMdc,
207  const Scalar& krnSwMdc,
208  Params& params)
209  {
210  // Maximum attainable oil saturation is 1-SWL
211  const auto Swco = params.Swl();
212  const double krwSw = 2.0; //Should not be used
213  params.gasOilParams().update(pcSwMdc - Swco, krwSw, krnSwMdc - Swco);
214  }
215 
225  template <class FluidState, class Evaluation = typename FluidState::Scalar>
226  static Evaluation pcgn(const Params& params,
227  const FluidState& fs)
228  {
229  // Maximum attainable oil saturation is 1-SWL.
230  const auto Sw = 1.0 - params.Swl() - decay<Evaluation>(fs.saturation(gasPhaseIdx));
231  return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
232  }
233 
243  template <class FluidState, class Evaluation = typename FluidState::Scalar>
244  static Evaluation pcnw(const Params& params,
245  const FluidState& fs)
246  {
247  const auto Sw = decay<Evaluation>(fs.saturation(waterPhaseIdx));
248  return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
249  }
250 
254  template <class ContainerT, class FluidState>
255  static void saturations(ContainerT& /*values*/,
256  const Params& /*params*/,
257  const FluidState& /*fluidState*/)
258  {
259  throw std::logic_error("Not implemented: saturations()");
260  }
261 
265  template <class FluidState, class Evaluation = typename FluidState::Scalar>
266  static Evaluation Sg(const Params& /*params*/,
267  const FluidState& /*fluidState*/)
268  {
269  throw std::logic_error("Not implemented: Sg()");
270  }
271 
275  template <class FluidState, class Evaluation = typename FluidState::Scalar>
276  static Evaluation Sn(const Params& /*params*/,
277  const FluidState& /*fluidState*/)
278  {
279  throw std::logic_error("Not implemented: Sn()");
280  }
281 
285  template <class FluidState, class Evaluation = typename FluidState::Scalar>
286  static Evaluation Sw(const Params& /*params*/,
287  const FluidState& /*fluidState*/)
288  {
289  throw std::logic_error("Not implemented: Sw()");
290  }
291 
307  template <class ContainerT, class FluidState>
308  static void relativePermeabilities(ContainerT& values,
309  const Params& params,
310  const FluidState& fluidState)
311  {
312  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
313 
314  values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
315  values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
316  values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
317  }
318 
322  template <class FluidState, class Evaluation = typename FluidState::Scalar>
323  static Evaluation krg(const Params& params,
324  const FluidState& fluidState)
325  {
326  // Maximum attainable oil saturation is 1-SWL.
327  const Evaluation Sw = 1.0 - params.Swl() - decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
328  return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
329  }
330 
334  template <class FluidState, class Evaluation = typename FluidState::Scalar>
335  static Evaluation krw(const Params& params,
336  const FluidState& fluidState)
337  {
338  const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
339  return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
340  }
341 
345  template <class FluidState, class Evaluation = typename FluidState::Scalar>
346  static Evaluation krn(const Params& params,
347  const FluidState& fluidState)
348  {
349  const Scalar Swco = params.Swl();
350 
351  const Evaluation Sw =
352  max(Evaluation(Swco),
353  decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
354 
355  const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
356 
357  const Evaluation Sw_ow = Sg + Sw;
358  const Evaluation kro_ow = relpermOilInOilWaterSystem<Evaluation>(params, fluidState);
359  const Evaluation kro_go = relpermOilInOilGasSystem<Evaluation>(params, fluidState);
360 
361  // avoid the division by zero: chose a regularized kro which is used if Sw - Swco
362  // < epsilon/2 and interpolate between the oridinary and the regularized kro between
363  // epsilon and epsilon/2
364  const Scalar epsilon = 1e-5;
365  if (scalarValue(Sw_ow) - Swco < epsilon) {
366  const Evaluation kro2 = (kro_ow + kro_go)/2;
367  if (scalarValue(Sw_ow) - Swco > epsilon/2) {
368  const Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
369  const Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
370 
371  return kro2*alpha + kro1*(1 - alpha);
372  }
373 
374  return kro2;
375  }
376 
377  return (Sg*kro_go + (Sw - Swco)*kro_ow) / (Sw_ow - Swco);
378  }
379 
383  template <class Evaluation, class FluidState>
384  static Evaluation relpermOilInOilGasSystem(const Params& params,
385  const FluidState& fluidState)
386  {
387  const Evaluation Sw =
388  max(Evaluation{ params.Swl() },
389  decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
390 
391  const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
392  const Evaluation So_go = 1.0 - (Sg + Sw);
393 
394  return GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
395  }
396 
400  template <class Evaluation, class FluidState>
401  static Evaluation relpermOilInOilWaterSystem(const Params& params,
402  const FluidState& fluidState)
403  {
404  const Evaluation Sw =
405  max(Evaluation{ params.Swl() },
406  decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
407 
408  const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
409  const Evaluation Sw_ow = Sg + Sw;
410 
411  return OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
412  }
413 
421  template <class FluidState>
422  static void updateHysteresis(Params& params, const FluidState& fluidState)
423  {
424  const Scalar Swco = params.Swl();
425 
426  const Scalar Sw = clampSaturation(fluidState, waterPhaseIdx);
427  const Scalar So = clampSaturation(fluidState, oilPhaseIdx);
428  const Scalar Sg = clampSaturation(fluidState, gasPhaseIdx);
429 
430  if (params.inconsistentHysteresisUpdate()) {
431  // NOTE: the saturations which are passed to update the hysteresis curves are
432  // inconsistent with the ones used to calculate the relative permabilities. We do
433  // it like this anyway because (a) the saturation functions of opm-core do it
434  // this way (b) the simulations seem to converge better (which is not too much
435  // surprising actually, because the time step does not start on a kink in the
436  // solution) and (c) the Eclipse 100 simulator may do the same.
437  //
438  // Though be aware that from a physical perspective this is definitively
439  // incorrect!
440  params.oilWaterParams().update(/*pcSw=*/ 1.0 - So,
441  /*krwSw=*/ 1.0 - So,
442  /*krnSw=*/ 1.0 - So);
443 
444  params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
445  /*krwSw=*/ 1.0 - Swco - Sg,
446  /*krnSw=*/ 1.0 - Swco - Sg);
447  }
448  else {
449  const Scalar Sw_ow = Sg + std::max(Swco, Sw);
450  const Scalar So_go = 1.0 - Sw_ow;
451 
452  params.oilWaterParams().update(/*pcSw=*/ Sw,
453  /*krwSw=*/ 1 - Sg,
454  /*krnSw=*/ Sw_ow);
455 
456  params.gasOilParams().update(/*pcSw=*/ 1.0 - Swco - Sg,
457  /*krwSw=*/ So_go,
458  /*krnSw=*/ 1.0 - Swco - Sg);
459  }
460  }
461 
462  template <class FluidState>
463  static Scalar clampSaturation(const FluidState& fluidState, const int phaseIndex)
464  {
465  const auto sat = scalarValue(fluidState.saturation(phaseIndex));
466  return std::clamp(sat, Scalar{0.0}, Scalar{1.0});
467  }
468 };
469 } // namespace Opm
470 
471 #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:59
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: EclDefaultMaterial.hpp:101
static Evaluation relpermOilInOilWaterSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/water system.
Definition: EclDefaultMaterial.hpp:401
static Evaluation krg(const Params &params, const FluidState &fluidState)
The relative permeability of the gas phase.
Definition: EclDefaultMaterial.hpp:323
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:226
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure.
Definition: EclDefaultMaterial.hpp:109
static Evaluation Sg(const Params &, const FluidState &)
The saturation of the gas phase.
Definition: EclDefaultMaterial.hpp:266
static Evaluation relpermOilInOilGasSystem(const Params &params, const FluidState &fluidState)
The relative permeability of oil in oil/gas system.
Definition: EclDefaultMaterial.hpp:384
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:134
static void saturations(ContainerT &, const Params &, const FluidState &)
The inverse of the capillary pressure.
Definition: EclDefaultMaterial.hpp:255
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition.
Definition: EclDefaultMaterial.hpp:117
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition: EclDefaultMaterial.hpp:276
static Evaluation krn(const Params &params, const FluidState &fluidState)
The relative permeability of the non-wetting (i.e., oil) phase.
Definition: EclDefaultMaterial.hpp:346
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:244
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: EclDefaultMaterial.hpp:97
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: EclDefaultMaterial.hpp:105
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition: EclDefaultMaterial.hpp:308
static Evaluation krw(const Params &params, const FluidState &fluidState)
The relative permeability of the wetting phase.
Definition: EclDefaultMaterial.hpp:335
static Evaluation Sw(const Params &, const FluidState &)
The saturation of the wetting (i.e., water) phase.
Definition: EclDefaultMaterial.hpp:286
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: EclDefaultMaterial.hpp:113
static void updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition: EclDefaultMaterial.hpp:422