avoid negative rs/rv max

[totto82]

Negative rs/rv max leads to failure in std::clamp inferDissolvedVaporisedRatio in debug mode.

[totto82]

jenkins build this please

[akva2]

This fixes https://ci.opm-project.org/job/opm-simulators-static-analysis/394/testReport/junit/(root)/debug_iterator/compareECLFiles_flow_9_3D_GINJ_GAS_MAX_EXPORT_MSW/

[totto82]

@akva2 merge? Or should we dig deeper to find the reason for rssat < 0?

[akva2]

I'd be worried we are burying some issue but I trust you more than my knee jerk so if you say it is okay..

[GitPaean]

Is this during well potential calculation with 1 atm bhp limit?

And also, do we know which OPM code trigger the error message?

/usr/include/c++/12/bits/stl_algo.h:3623: constexpr const _Tp& std::clamp(const _Tp&, const _Tp&, const _Tp&) [with _Tp = double]: Assertion '!(__hi < __lo)' failed.

[GitPaean]

Thanks for @akva2 pointing out.

On my computer, it only triggers if we use the old 16-day time step limits and it runs through with default setting of the master branch.

As far as I can see, there is something seriously wrong with this well during the occasion when it crashes.

previously, the bhp of the well INJ2 is 235 bar, but at the moment it crashes, the segment 0 has a pressure of about 1.4 bar.

the segment pressures for the whole well look the following,

It does not add up there the pressure difference between the segments (but it is not converged either, so it might be fine theoretically with huge residual). And it is almost a gas well through the segments.

I am just curious while there is no objection to the proposed solution. Maybe some debug output will be nice.

[totto82]

jenkins build this please

[totto82]

jenkins build this please

[bska]

While I agree that a bit of resilience against unexpected values is nice, this PR will only hide a deeper issue. If we get negative Rs or Rv values then we should probably abort the simulation.

[atgeirr]

If we get negative Rs or Rv values then we should probably abort the simulation.

I would think it quite possible to get this when extrapolating to low pressures. Therefore I do not think it is necessary a bug, other than failing to guard against it, which is what this PR tries to address.

[bska]

If we get negative Rs or Rv values then we should probably abort the simulation.

I would think it quite possible to get this when extrapolating to low pressures.

PR OPM/opm-common#3779 was supposed to fix that particular problem. If it still exists then we definitely have a bug.

Therefore I do not think it is necessary a bug, other than failing to guard against it, which is what this PR tries to address.

No, I'd still like to understand what's happening here. This PR just hides a problem.

[atgeirr]

PR OPM/opm-common#3779 was supposed to fix that particular problem. If it still exists then we definitely have a bug.

Good point, I agree with your assessment then.

[GitPaean]

I am testing the master branch, it runs through the 9_3D_GINJ_GAS_MAX_EXPORT_MSW.DATA without problem.

It is still running the proved gcc12 data file, it was reported problem at report step 24 (with 4 mpi processes), now it is at the report step 550.

I am trying to pick this topic top, any suggestion how to reveal the problem?

[akva2]

you need to build with stdlib assertions.

[GitPaean]

We have a negative rs here and the simulation crashes with 9_3D_GINJ_GAS_MAX_EXPORT_MSW

template <typename Scalar, typename Rates>
std::pair<Scalar, Scalar>
dissolvedVaporisedRatio(const int    io,
                        const int    ig,
                        const Scalar rs,
                        const Scalar rv,
                        const Rates& surface_rates)
{
    if ((io < 0) || (ig < 0)) {
        return { rs, rv };
    }
    auto eps = std::copysign(1.0e-15, surface_rates[io]);
    const auto Rs = surface_rates[ig] / (surface_rates[io] + eps);

    eps = std::copysign(1.0e-15, surface_rates[ig]);
    const auto Rv = surface_rates[io] / (surface_rates[ig] + eps);

    return {
        std::clamp(static_cast<Scalar>(Rs), Scalar{0.0}, rs),
        std::clamp(static_cast<Scalar>(Rv), Scalar{0.0}, rv)
    };
}

you need to build with stdlib assertions.

Thanks for the information. What happens when building without stdlib assertions, it just runs with kind of undefined behavoir?

[akva2]

yup

[GitPaean]

And it happens with 1 bar bhp constraint,

and the pressure drop between the segments for this well looks weird, while in theory is possible but most likely problematic.

But the problem is that we get negative rsMax with 1 bar pressure.

[GitPaean]

The weird pressure drop is from an incomplete initialization of the well states after changing to 1 bar bhp.

[GitPaean]

And the negative value for rsMax simply from the extrapolation at the lower pressure end. The PVTO table has a lower pressure end at 50 bar.

[GitPaean]

PVTG looks like patched til 1 bar, probably with https://github.com/OPM/opm-common/pull/3779. Maybe we should do it for the keyword PVTO?

[totto82]

PVTG looks like patched til 1 bar, probably with https://github.com/OPM/opm-common/pull/3779. Maybe we should do it for the keyword PVTO?

That sounds like a good idea. I will keep the PR open (if somebody needs a quick fix) but change it to draft so nobody merges it.