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