diff --git a/dumux/porousmediumflow/richards/implicit/volumevariables.hh b/dumux/porousmediumflow/richards/implicit/volumevariables.hh index c5d9feed9df6838466aa42c99951daa1ecff1308..bf993368f1dc85369b13887a77a33488c0db03a6 100644 --- a/dumux/porousmediumflow/richards/implicit/volumevariables.hh +++ b/dumux/porousmediumflow/richards/implicit/volumevariables.hh @@ -87,8 +87,9 @@ public: const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol); relativePermeabilityWetting_ = MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx)); - //! precompute the minimum capillary pressure (entry pressure) - minPc_ = MaterialLaw::pc(materialParams, /*Sw=*/ 1.0); + // precompute the minimum capillary pressure (entry pressure) + // needed to make sure we don't compute unphysical capillary pressures and thus saturations + minPc_ = MaterialLaw::endPointPc(materialParams); pn_ = problem.nonWettingReferencePressure(); porosity_ = problem.spatialParams().porosity(element, scv, elemSol); permeability_ = problem.spatialParams().permeability(element, scv, elemSol); @@ -113,8 +114,11 @@ public: fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]); // compute the capillary pressure to compute the saturation - // the material law takes care of cases where pc gets negative - const Scalar pc = problem.nonWettingReferencePressure() - fluidState.pressure(wPhaseIdx); + // make sure that we the capillary pressure is not smaller than the minimum pc + // this would possibly return unphysical values from regularized material laws + using std::max; + const Scalar pc = max(MaterialLaw::endPointPc(materialParams), + problem.nonWettingReferencePressure() - fluidState.pressure(wPhaseIdx)); const Scalar sw = MaterialLaw::sw(materialParams, pc); fluidState.setSaturation(wPhaseIdx, sw);