From 893538c9f30eca8864d242dade58dbbd149fe762 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 13 Aug 2015 21:18:07 +0200 Subject: [PATCH 1/2] EclDefaultMaterial: use inconsistent saturations for the hysteresis update with this, I got slightly better performance than the opm-core master version if flow is tasked on simulating the full Norne deck. Be aware that from the physical POV, this is wrong. --- .../EclDefaultMaterial.hpp | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index ad5d475ab..b727527e5 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -327,18 +327,21 @@ public: { typedef MathToolbox FsToolbox; - Scalar Swco = params.Swl(); - - Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx)); + Scalar So = FsToolbox::value(fluidState.saturation(oilPhaseIdx)); Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx)); - Sw = std::min(1.0, std::max(Swco, Sw)); Sg = std::min(1.0, std::max(0.0, Sg)); - Scalar Sw_ow = Sg + Sw; - Scalar So_go = 1 - Sw_ow + Swco; - - params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow); - params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg); + // FIXME: the saturations which are passed to update the hysteresis curves are + // inconsistent with the ones used to calculate the relative permabilities. We do + // it like this anyway because (a) the saturation functions of opm-core do it + // this way (b) the simulations seem to converge better (which is not too much + // surprising actually, because the time step does not start on a kink in the + // solution) and (c) the Eclipse 100 simulator may do the same. + // + // Though be aware that from a physical perspective this is definitively + // incorrect! + params.oilWaterParams().update(/*pcSw=*/1 - So, /*krwSw=*/1 - So, /*krn_Sw=*/1 - So); + params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg); } }; } // namespace Opm From a8aaa168e2c62e9b246a5678b0b74a690271c633 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 14 Aug 2015 09:33:28 +0200 Subject: [PATCH 2/2] EclDefaultMaterial: allow the parameters to decide if consistent or inconsistent hysteresis are used --- .../EclDefaultMaterial.hpp | 43 ++++++++++++------- .../EclDefaultMaterialParams.hpp | 12 ++++++ 2 files changed, 39 insertions(+), 16 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index b727527e5..021297ebe 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -181,10 +181,7 @@ public: typedef MathToolbox FsToolbox; const auto& Sw = FsToolbox::template toLhs(fs.saturation(waterPhaseIdx)); - Valgrind::CheckDefined(Sw); - const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw); - Valgrind::CheckDefined(result); - return result; + return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw); } /*! @@ -327,21 +324,35 @@ public: { typedef MathToolbox FsToolbox; + Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx)); Scalar So = FsToolbox::value(fluidState.saturation(oilPhaseIdx)); Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx)); - Sg = std::min(1.0, std::max(0.0, Sg)); - // FIXME: the saturations which are passed to update the hysteresis curves are - // inconsistent with the ones used to calculate the relative permabilities. We do - // it like this anyway because (a) the saturation functions of opm-core do it - // this way (b) the simulations seem to converge better (which is not too much - // surprising actually, because the time step does not start on a kink in the - // solution) and (c) the Eclipse 100 simulator may do the same. - // - // Though be aware that from a physical perspective this is definitively - // incorrect! - params.oilWaterParams().update(/*pcSw=*/1 - So, /*krwSw=*/1 - So, /*krn_Sw=*/1 - So); - params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg); + if (params.inconsistentHysteresisUpdate()) { + Sg = std::min(1.0, std::max(0.0, Sg)); + // NOTE: the saturations which are passed to update the hysteresis curves are + // inconsistent with the ones used to calculate the relative permabilities. We do + // it like this anyway because (a) the saturation functions of opm-core do it + // this way (b) the simulations seem to converge better (which is not too much + // surprising actually, because the time step does not start on a kink in the + // solution) and (c) the Eclipse 100 simulator may do the same. + // + // Though be aware that from a physical perspective this is definitively + // incorrect! + params.oilWaterParams().update(/*pcSw=*/1 - So, /*krwSw=*/1 - So, /*krn_Sw=*/1 - So); + params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg); + } + else { + Scalar Swco = params.Swl(); + Sw = std::min(1.0, std::max(0.0, Sw)); + Sg = std::min(1.0, std::max(0.0, Sg)); + + Scalar Sw_ow = Sg + std::max(Swco, Sw); + Scalar So_go = 1 + Swco - Sw_ow; + + params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow); + params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg); + } } }; } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp index 8e331ce7f..dc43f9e71 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp @@ -124,6 +124,18 @@ public: Scalar Swl() const { assertFinalized_(); return Swl_; } + /*! + * \brief Specify whether inconsistent saturations should be used to update the + * hysteresis parameters. + * + * Returning 'true' is wrong from a physical point of view because the saturations + * which are used to update the hysteresis parameters are calculated differently than + * the ones used to calculate the relperms and capillary pressures. Since Eclipse + * E100 probably uses inconsistent saturations, we return true here anyway. + */ + bool inconsistentHysteresisUpdate() const + { return true; } + private: #ifndef NDEBUG void assertFinalized_() const