All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
32 #include <opm/common/Valgrind.hpp>
34 
35 #include <opm/common/Exceptions.hpp>
36 #include <opm/common/ErrorMacros.hpp>
37 
38 #include <algorithm>
39 
40 namespace Opm {
41 
55 template <class TraitsT,
56  class GasOilMaterialLawT,
57  class OilWaterMaterialLawT,
58  class ParamsT = EclDefaultMaterialParams<TraitsT,
59  typename GasOilMaterialLawT::Params,
60  typename OilWaterMaterialLawT::Params> >
61 class EclDefaultMaterial : public TraitsT
62 {
63 public:
64  typedef GasOilMaterialLawT GasOilMaterialLaw;
65  typedef OilWaterMaterialLawT OilWaterMaterialLaw;
66 
67  // some safety checks
68  static_assert(TraitsT::numPhases == 3,
69  "The number of phases considered by this capillary pressure "
70  "law is always three!");
71  static_assert(GasOilMaterialLaw::numPhases == 2,
72  "The number of phases considered by the gas-oil capillary "
73  "pressure law must be two!");
74  static_assert(OilWaterMaterialLaw::numPhases == 2,
75  "The number of phases considered by the oil-water capillary "
76  "pressure law must be two!");
77  static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
78  typename OilWaterMaterialLaw::Scalar>::value,
79  "The two two-phase capillary pressure laws must use the same "
80  "type of floating point values.");
81 
82  static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
83  "The gas-oil material law must implement the two-phase saturation "
84  "only API to for the default Ecl capillary pressure law!");
85  static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
86  "The oil-water material law must implement the two-phase saturation "
87  "only API to for the default Ecl capillary pressure law!");
88 
89  typedef TraitsT Traits;
90  typedef ParamsT Params;
91  typedef typename Traits::Scalar Scalar;
92 
93  static const int numPhases = 3;
94  static const int waterPhaseIdx = Traits::wettingPhaseIdx;
95  static const int oilPhaseIdx = Traits::nonWettingPhaseIdx;
96  static const int gasPhaseIdx = Traits::gasPhaseIdx;
97 
100  static const bool implementsTwoPhaseApi = false;
101 
104  static const bool implementsTwoPhaseSatApi = false;
105 
108  static const bool isSaturationDependent = true;
109 
112  static const bool isPressureDependent = false;
113 
116  static const bool isTemperatureDependent = false;
117 
120  static const bool isCompositionDependent = false;
121 
136  template <class ContainerT, class FluidState>
137  static void capillaryPressures(ContainerT& values,
138  const Params& params,
139  const FluidState& state)
140  {
141  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
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  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  pcSwMdc = params.gasOilParams().pcSwMdc();
193  krnSwMdc = params.gasOilParams().krnSwMdc();
194 
195  Valgrind::CheckDefined(pcSwMdc);
196  Valgrind::CheckDefined(krnSwMdc);
197  }
198 
199  /*
200  * Hysteresis parameters for gas-oil
201  * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
202  * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
203  * \param params Parameters
204  */
205  static void setGasOilHysteresisParams(const Scalar& pcSwMdc,
206  const Scalar& krnSwMdc,
207  Params& params)
208  {
209  const double krwSw = 2.0; //Should not be used
210  params.gasOilParams().update(pcSwMdc, krwSw, krnSwMdc);
211  }
212 
222  template <class FluidState, class Evaluation = typename FluidState::Scalar>
223  static Evaluation pcgn(const Params& params,
224  const FluidState& fs)
225  {
226  const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
227  return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
228  }
229 
239  template <class FluidState, class Evaluation = typename FluidState::Scalar>
240  static Evaluation pcnw(const Params& params,
241  const FluidState& fs)
242  {
243  const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
244  return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
245  }
246 
250  template <class ContainerT, class FluidState>
251  static void saturations(ContainerT& /*values*/,
252  const Params& /*params*/,
253  const FluidState& /*fluidState*/)
254  {
255  OPM_THROW(std::logic_error, "Not implemented: saturations()");
256  }
257 
261  template <class FluidState, class Evaluation = typename FluidState::Scalar>
262  static Evaluation Sg(const Params& /*params*/,
263  const FluidState& /*fluidState*/)
264  {
265  OPM_THROW(std::logic_error, "Not implemented: Sg()");
266  }
267 
271  template <class FluidState, class Evaluation = typename FluidState::Scalar>
272  static Evaluation Sn(const Params& /*params*/,
273  const FluidState& /*fluidState*/)
274  {
275  OPM_THROW(std::logic_error, "Not implemented: Sn()");
276  }
277 
281  template <class FluidState, class Evaluation = typename FluidState::Scalar>
282  static Evaluation Sw(const Params& /*params*/,
283  const FluidState& /*fluidState*/)
284  {
285  OPM_THROW(std::logic_error, "Not implemented: Sw()");
286  }
287 
303  template <class ContainerT, class FluidState>
304  static void relativePermeabilities(ContainerT& values,
305  const Params& params,
306  const FluidState& fluidState)
307  {
308  typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
309 
310  values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
311  values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
312  values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
313  }
314 
318  template <class FluidState, class Evaluation = typename FluidState::Scalar>
319  static Evaluation krg(const Params& params,
320  const FluidState& fluidState)
321  {
322  const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
323  return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
324  }
325 
329  template <class FluidState, class Evaluation = typename FluidState::Scalar>
330  static Evaluation krw(const Params& params,
331  const FluidState& fluidState)
332  {
333  const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
334  return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
335  }
336 
340  template <class FluidState, class Evaluation = typename FluidState::Scalar>
341  static Evaluation krn(const Params& params,
342  const FluidState& fluidState)
343  {
344  Scalar Swco = params.Swl();
345 
346  Evaluation Sw =
347  Opm::max(Evaluation(Swco),
348  Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
349  Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
350 
351  Evaluation Sw_ow = Sg + Sw;
352  Evaluation So_go = 1.0 - Sw_ow;
353  const Evaluation& kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
354  const Evaluation& kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
355 
356  // avoid the division by zero: chose a regularized kro which is used if Sw - Swco
357  // < epsilon/2 and interpolate between the oridinary and the regularized kro between
358  // epsilon and epsilon/2
359  const Scalar epsilon = 1e-5;
360  if (Opm::scalarValue(Sw_ow) - Swco < epsilon) {
361  Evaluation kro2 = (kro_ow + kro_go)/2;;
362  if (Opm::scalarValue(Sw_ow) - Swco > epsilon/2) {
363  Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
364  Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
365  return kro2*alpha + kro1*(1 - alpha);
366  }
367 
368  return kro2;
369  }
370  else
371  return (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
372  }
373 
381  template <class FluidState>
382  static void updateHysteresis(Params& params, const FluidState& fluidState)
383  {
384  Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
385  Scalar So = Opm::scalarValue(fluidState.saturation(oilPhaseIdx));
386  Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
387 
388  if (params.inconsistentHysteresisUpdate()) {
389  Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
390  // NOTE: the saturations which are passed to update the hysteresis curves are
391  // inconsistent with the ones used to calculate the relative permabilities. We do
392  // it like this anyway because (a) the saturation functions of opm-core do it
393  // this way (b) the simulations seem to converge better (which is not too much
394  // surprising actually, because the time step does not start on a kink in the
395  // solution) and (c) the Eclipse 100 simulator may do the same.
396  //
397  // Though be aware that from a physical perspective this is definitively
398  // incorrect!
399  params.oilWaterParams().update(/*pcSw=*/1 - So, /*krwSw=*/1 - So, /*krn_Sw=*/1 - So);
400  params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg);
401  }
402  else {
403  Scalar Swco = params.Swl();
404  Sw = std::min(Scalar(1.0), std::max(Scalar(0.0), Sw));
405  Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
406 
407  Scalar Sw_ow = Sg + std::max(Swco, Sw);
408  Scalar So_go = 1 + Sw_ow;
409 
410  params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
411  params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);
412  }
413  }
414 };
415 } // namespace Opm
416 
417 #endif
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
static const bool isCompositionDependent
Specify whether the quantities defined by this material law are dependent on the phase composition...
Definition: EclDefaultMaterial.hpp:120
static const bool implementsTwoPhaseApi
Specify whether this material law implements the two-phase convenience API.
Definition: EclDefaultMaterial.hpp:100
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:137
static Evaluation krg(const Params &params, const FluidState &fluidState)
The relative permeability of the gas phase.
Definition: EclDefaultMaterial.hpp:319
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:240
static Evaluation Sg(const Params &, const FluidState &)
The saturation of the gas phase.
Definition: EclDefaultMaterial.hpp:262
Default implementation for the parameters required by the default three-phase capillary pressure mode...
static void saturations(ContainerT &, const Params &, const FluidState &)
The inverse of the capillary pressure.
Definition: EclDefaultMaterial.hpp:251
static Evaluation krw(const Params &params, const FluidState &fluidState)
The relative permeability of the wetting phase.
Definition: EclDefaultMaterial.hpp:330
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:223
static const bool isTemperatureDependent
Specify whether the quantities defined by this material law are temperature dependent.
Definition: EclDefaultMaterial.hpp:116
static Evaluation krn(const Params &params, const FluidState &fluidState)
The relative permeability of the non-wetting (i.e., oil) phase.
Definition: EclDefaultMaterial.hpp:341
static const bool isPressureDependent
Specify whether the quantities defined by this material law are dependent on the absolute pressure...
Definition: EclDefaultMaterial.hpp:112
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition: EclDefaultMaterial.hpp:272
static const bool isSaturationDependent
Specify whether the quantities defined by this material law are saturation dependent.
Definition: EclDefaultMaterial.hpp:108
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition: EclDefaultMaterial.hpp:304
static Evaluation Sw(const Params &, const FluidState &)
The saturation of the wetting (i.e., water) phase.
Definition: EclDefaultMaterial.hpp:282
Implements the default three phase capillary pressure law used by the ECLipse simulator.
Definition: EclDefaultMaterial.hpp:61
static const bool implementsTwoPhaseSatApi
Specify whether this material law implements the two-phase convenience API which only depends on the ...
Definition: EclDefaultMaterial.hpp:104
static void updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition: EclDefaultMaterial.hpp:382