Issue in numerical backend for _T() under certain multiphase conditions

[jranalli]

Ok, I'm moving this discussion here because it's going to benefit from some conversation that shows code.

I've found a few cases where we are unable to calculate multiphase p&d properties, because of bracketing errors on the numerical inversion routines.

Here are two cases that fail:

sub = pm.get('mp.H2O')
p = 2497
d = 525
T = sub.T(p=p, d=d)  # expect 1026.xxx, get bracket error

p = 2497
d = 1090
T = sub.T(p=p, d=d)  # expect 291.xxx, get bracket error

The issue occurs because of code within the function mp1._T(), which uses _hybrid1() to compute an inversion seeking the temperature of interest using _p(T, d). For supercritical pressures, on line 2115 it calculates the limits for that inversion as the upper and lower limits of the temperature for the substance.

Because we're supercritical, there are cases where the target density will be under the dome if we really evaluate at the lowest possible temperature, and according to the documentation, _p() doesn't produce valid results under the dome. My observation for the first test case is that the reported value from _p(T=Tmin, d=525) is ridiculous and on the order of 10^20. So that's why we get a bracket error.

You may have expertise that can provide a different/better fix, but I've pushed one viable solution to my working branch. Basically, I'm computing new lower temperature limits for the _hybrid1() iteration that are based on temperature at which our density intersects the dome (if it does). Computing that value of ds requires its own inversion though and is therefore expensive. So I made it conditional on the return, such that it only tries the fix after we've already gotten a bracket error. See this commit.

I haven't exhaustively tested this out yet, but it does solve the problem for all of the individual conditions I've tested for water. It needs to be run on arrays and other substances before we can be sure.

Here's the test code I ran against for now pending further discussion:

import pyromat as pm
sub = pm.get('mp.H2O')
print(sub.T(p=2497, d=525))   # errors with the old method, works with the new, expect around 1026
print(sub.T(p=2497, d=1090))  # density doesn't intersect the dome, should work for old or new method, expect around 291
print(sub.T(p=2497, d=1100))  # really out of range, sub.d(p=2497, T=Tmin is around 1097), expect errors

[chmarti1]

This is a new problem created by an edit I made in the modifications for v2.2.0. I used to use _tditer() instead of _p(). _tditer() correctly resolves values under the dome, so these weird conditions don't crop up. I'll either transition back to _tditer() or I'll implement a fix to make _p() avoid the issue. _p() is preferred for numerical efficiency. _tditer() is apparently preferred for numerical stability.

[jranalli]

Well, if that's the case, maybe it's worth doing the try/except approach. _hybrid1 definitely handles the vast majority of cases and you can benefit from its efficiency. I only found these via a test code that pointwise crawls the space. So simplest approach for best of both worlds would be something like:

try:
  self._hybrid1(args)
except pm.utility.PMParamError as e:
    self._tditer(args)

Related by not identical issue, these types of cases crop up in the _d() calculation as well. I was able to resolve those as well by changing temperature limits. A few test cases for _d() are:

sub = pm.get('mp.CH4')
print(sub.d(p=0.11707696, T=624.375)) # T supercritical, p low
print(sub.d(p=9999.0, T=624.375)) # T&p supercritical

What I found here is that at high temperatures and various pressures, you get insane pressures at the high density boundary and it screws up the initiation of the iteration. My fix is uglier this time, but it might at least give you an idea of what kinds of changes I had to do to the iteration limits to make it work out.

Commit here. One part of the fix is to use the critical density as a limit when temperature is above the dome and pressure is below. The other part is to come up with a limit for when both p & T are above the dome. It works to compute the density at Tc & pmax, which would be an expected minimum value, but it results in an ugly recursion and I don't recommend using my exact method there.

[jranalli]

Since I'm documenting other numerics issues here, let me add these as well.

One involves calling things based on entropy.

sub = pm.get('mp.CH4')
ref = sub.d(T=223, s=-5.529)   # should be legit, 
# value expected is around 460 based on the forward calculation with p=2475 & T=223

The underlying codes do an iteration of _s based on density. But if you plot entropy as a function of density at that temperature you can see that it's got two possible density values that yield the entropy of interest and so your iteration code breaks. Better limits could maybe work here.

Something similar can happen with pressure and temperature for some cases. Here's the test condition

sub = pm.get('mp.N2')
sub.T(p=1.6e4, d=1355)  # errors out, should come out to be around 68

Leads to iteration of _p() based on T. Plotting pressure versus temperature shows that at very low temperatures there are two possible values of temperature that correspond to that pressure.

Just to be 100% clear, I'm not trying to nitpick here. I wrote some code to try to span the parameter space for all substances and I'm just diving into the test cases that errored out and trying to provide them so we can figure out what's going on with the numerics. No pressure as these aren't active cases I'm trying to use or anything.

[chmarti1]

Version 2.2.1 corrects the issue. For the time being, I have reverted to _tditer() instead of _p() for the error function. The alternative requires an algorithm to calculate the saturation temperature given the density, so we can establish a minimum temperature that keeps the iteration outside of the dome when the pressure is super-critical. That's really the right way to solve the problem efficiently. Writing that algorithm isn't horribly difficult, but implementing it requires some significant testing. This can remain a goal for the future.