CO2Store triggers taking square root of a negative value for high temperatures (>304 K)

[blattms]

This is the case for all activity models in ACTCO2S. In my case I experienced NANS in the well iterations. No warning or exceptions before this. Higher temperatures are clear due to the RTEMPV table.

If value 3 (non-iterated Duan-Sun according to Spycher and Pruess) is used we will end up in Brine_CO2.hpp#L129. For temperatures higher than 304K we will try to compute the power 1.5 on a negative value in CO2.hpp#L148, which has at least undefined behavior, see documentation of std::pow

The floating point exception occurs when we try

(gdb) 
#7  0x0000564814fdab70 in Opm::MathToolbox<double>::pow (base=-0.057369892561418778, exp=1.5) at 
  opm-common/opm/material/common/MathToolbox.hpp:271
271	    { return std::pow(base, exp); }

because the temperature is higher than the critical on in (1-T/critT).

The same holds for 2, because code only starts to iterate if temperature is over 372.15K, see Brine_CO2.hpp#L115. Unfortunately, we still end up in CO2.hpp#148 later.

Finally, value 1 does not work, because Brine_CO2.hpp#L579 actively prevents the fixed point iteration.

Is this intentional? Maybe @daavid00 knows?

We should probably do proper error/exception handling for this case. Is there any reason to not use the fixed point iterations for values 1 and 2?

@gdfldm There should probably be a prominent warning and a known issue in the manual.

[totto82]

This is suppose to not happen as there is if(highTemp) checks but if you get NANS in your solution it points to some issues. Could you share the model for faster debugging? I will look into this on Monday.

[daavid00]

Hi Markus, I did not implement that code myself, but as Tor said we will look at this issue next week.

[blattms]

Thanks a lot for the quick answers.

I can probably not share the model. As outlined above (see links) we do set hightemp wrong or at least consistently. I will prepare a draft PR with some obvious fixes (for value 2).

I do not know whether the nans are due to this. It was just uncovered when asking the os to send signals for floating point errors and catching this.

[totto82]

I do not know whether the nans are due to this.

That was good to know. I was really puzzled by this. Because the vaporization pressure is only used to check which parameters to use and nans here should not effect the results at all. An obvious fix is to only compute the vaporization pressure when you need it (i.e. when t_celsius < 31). That will remove the issue the way I read the code and also give a small improvement in performance for most cases as computing the vaporization pressure is none trivial. If you are not very eager to do it your self I can provide a fix on Monday.

[gdfldm]

There should probably be a prominent warning and a known issue in the manual.

I have created a PR opm-reference-manual#524

[blattms]

#4514 adds an explicit throw for wrong usage of the std::pow function when calculating the vapor pressure.

@totto82 that is good to know. I am a step further. I now see nan weights in cprw. Others have seen that too and tried a fix in OPM/opm-simulators#6053. I will do some further digging....

[blattms]

Just to clarify:
The actual problem is somewhere else. I got to that once I worked around this one and the division by zero in cprw (or switching to other weights there).
The actual nan that is hurting us arises during well assembly. The well is not stopped and has nonzero rate, but the rate at perforations that have a depth beyond a certain threshold are zero.

The above was when using ACTCO2S being 2

[blattms]

Well, it seems like we actually get to this place for ACTCO2S being 3 also when calculating the equilibrium. Namely, when determining the mole fractions used to compute the saturated gas dissolution factor

[totto82]

Can you check if this fix your problem? Shouldn't change any results, but it avoids unnecessary calls to CO2::vaporPressure() for T > TCritical.

#4515

[blattms]

Looks good to me. I'll try that.