Catch singular explosions in saturation_PHSU_pure

[msaitta-mpr]

Description of the Change

Sometimes, the Akasaka solver has issues. This is usually caught and then it is retried with a new omega. However, sometimes it goes bad because the J matrix is singular, and this is not caught because the error is not recalculated.

This PR recalculates the error term to prevent a bad result.

Benefits

This prevents saturation_PHSU_pure from providing wrong answers.

This closes #2245.

Possible Drawbacks

This slows calculations because it requires that the error term be recalculated an extra time for each update. However, it prevents bad results, which is usually worth it.

This does not solve the underlying issue of a singular matrix forming in the Akasaka solver. Therefore, we may want to look into that in the future.

Verification Process

I ran the following C++ code and verified that the results look smooth.

// Test issue #2245
#pragma once
#include "CoolProp.h"
#include "AbstractState.h"
#include <iostream>
#include <memory>
#include <vector>

std::vector <double> linspace(double start, double stop, int num_pts){
    std::vector<double> vals (num_pts);
    for (int ii=0; ii < num_pts; ii++){
        vals.at(ii) = start + (stop - start) * double(ii) / double(num_pts-1);
    }
    return vals;
}

int main()
{
    std::shared_ptr<CoolProp::AbstractState> AS = std::shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "Ammonia"));
    std::vector<double> p_vals = linspace(50000.0, 100000, 100);
    for (double p : p_vals) {
        AS->update(CoolProp::PQ_INPUTS, p, 1);
        std::cout << p << "  " << AS->T() << "  " << CoolProp::PropsSI("T", "P", p, "Q", 1.0, "Ammonia") << std::endl;
    }
    return 0;
}

Applicable Issues

Closes #2245

[danielkwalsh]

Maybe it would be a good idea to explain your intent of recomputing the error as a code comment? I don't fully understand the problem, but based on the message in your PR, it sounds like this extra check is superfluous provided that the Akasaka solver's Jacobian works as expected. When that issue is addressed in the future, I'm hoping that my proposed comment could be a reminder that the changes you made here could be omitted to cut out the extra error calculation you described. Just a thought.

[ibell]

Those additional checks should be very cheap, so if this fixes the iteration failures then I would say we should merge them.

[ibell]

But I also agree with @danielkwalsh that a few comments would be good

[msaitta-mpr]

I added a comment which I hope captures what you think should be communicated. If I am off, please let me know.

[danielkwalsh]

@msaitta-mpr yeah, I think this looks great. Thanks for adding this!

[danielkwalsh]

@ibell @msaitta-mpr Seems like there are two tests that are failing. Is that expected?

[msaitta-mpr]

These also fail in the master branch. I have been meaning to track down the cause and fix it, but I haven't had a chance.

[ibell]

Should I go ahead and merge, or wait?

[msaitta-mpr]

On a hunch, I reran my test case with much greater point density and found that there are still a few points that do not behave. So, while the PR has so far fixed a lot of these issues, there is still some underlying issue. I will try to look into it. I think it makes sense to leave the PR open for the moment. Hopefully these commits can serve as an interim patch for users who need a quick fix.

[ibell]

Hmm that's weird. I wonder what is going on. If this improves this, it might be worth it to merge it. But I also wonder whether there is something wrong with the ancillary functions it fails so hard. I don't recall similar problems with other EOS.

[danielkwalsh]

@ibell @msaitta-mpr Have we assessed whether this is an ammonia-specific issue, or does it exist for other fluids as well?

[msaitta-mpr]

Hmm that's weird. I wonder what is going on. If this improves this, it might be worth it to merge it. But I also wonder whether there is something wrong with the ancillary functions it fails so hard. I don't recall similar problems with other EOS.

Surprisingly, I do not think that it is the ancillary. Looking at a couple of the failing points, the ancillary is within 0.01% of the correct solution.

The issue seems to be in the update of SatL and SatV. A DT update is performed, but the resultant state has a p value that is not correct (i.e., what you get if you just give CoolProp these values directly).

[ibell]

I don't understand how the pressure can be not correct for input temperature and density. Can you figure out what is going on? Could a derivative be not correct somehow?

…

Message ID: ***@***.***>

[msaitta-mpr]

The "error" in the DT update seems to be due to the imposed phase on SatL and SatV. I added a final update without this phase specification. This now gives a completely smooth result.

[ibell]

This doesn't make any sense. The phases should always be set to liquid and vapor to avoid a nested phase equilibrium calculation. It doesn't make sense to not impose the phase. Can you see whether the obtained pressures are very different when doing a calc with and without the imposing of phase? I may need to dig in a bit here. It could be a problem with the EOS I fear.

[msaitta-mpr]

That is what I would have thought as well, but you can get some strange results even when doing DT updates. Here is some example Python:

import CoolProp
rho = 16425.375470309176
t = 337.43858705452624
state = CoolProp.AbstractState("HEOS", "AMMONIA")
state.update(CoolProp.DmolarT_INPUTS, rho, t)
print(state.p())
state.specify_phase(CoolProp.iphase_liquid)
state.update(CoolProp.DmolarT_INPUTS, rho, t)
print(state.p())

This prints

2898759.978561249
56513.00000000532

whereas I would have expected the same value to be printed in either case.

Is the issue in the equation of state, or in the way that the pressure is found with some underlying phase assumption?

I will try to take a look into the intermediate derivative values.

[ibell]

Well, in your test case, the first set of inputs are two-phase, so that is why the results are different:

import CoolProp
rho = 16425.375470309176
t = 337.43858705452624
state = CoolProp.AbstractState("HEOS", "AMMONIA")
state.update(CoolProp.DmolarT_INPUTS, rho, t)
print(state.p(), state.Q())
state.specify_phase(CoolProp.iphase_liquid)
state.update(CoolProp.DmolarT_INPUTS, rho, t)
print(state.p(), state.Q())

yielding

2898759.978561249 0.04091462363285387
56513.00000000532 -1.0

and this is general true. If the inputs are two-phase you should expect different answers.

[danielkwalsh]

Hey, all. What's the status on this?

[msaitta-mpr]

@danielkwalsh Sorry for the slow progress on this, but I am a bit stuck. In short, the current version of the PR fixes the immediate issue, but it does so in an awkward way that might not be best for merging into master. If you need a short term fix, you can probably just use this branch as it seems to make all of my tests work.

To elaborate a bit further, there seem to have been two issues at play:

1. Sometimes, the saturation solver would diverge in the final iteration, causing a value that was incorrect to be passed. This was corrected in 41739c5 and c1dc498.
2. Sometimes, the saturation solver converges to a point where the provided residuals are all 0 (or very close to it). However, rather than finding points on the saturation curve, it has instead found points within the vapor dome. My solution to this was to force it to recompute the pressure, which seems to work. However, ibell believes that this might have an issue of a nested phase calculation. Therefore, there exists a way to determine that the found DT combination is incorrect, but I have not yet had the time to come up with a clever way of doing so.

[danielkwalsh]

@msaitta-mpr Perhaps I am missing something, but the saturation pressure applies to any point under the dome. Any isobar here will also be an isotherm, since the liquid and gaseous phases are in equilibrium. While I don't know the exact behavior of the solver being used here, nor do I understand the exact issue you are seeing, I could imagine that due to this degeneracy, asking for a state at a given pressure on the phase line could be multivalued due to this degeneracy, and consequently, the solver may give different results under different circumstances corresponding to different enthalpies. This is fine, since we are only interested in determining for a fluid at a given pressure, its corresponding vapor temperature, or vice versa. I'm happy to help try to debug this issue with you if you'd like, but unfortunately I don't understand the backend of the code well enough to do so on my own. Perhaps @ibell has some thoughts?

On the CoolProp low-level API page, I find the following:

Warning When specifying an imposed phase, it is absolutely critical that the input pair actually lie within the imposed phase region. If an incorrect phase is imposed for the given input pair, update() may throw unexpected errors or incorrect results may possibly be returned from the property functions. If the state point phase is not absolutely known, it is best to let CoolProp determine the phase at the cost of some computational efficiency.

[msaitta-mpr]

@danielkwalsh The low-level API page that you linked is key to this issue. We want to find the phase boundary. To do so, we perform repeated evaluations with two different states with imposed phases: one assumed liquid and the other gas. This works well, until they are evaluated within the two phase region. As the documentation states, they may return incorrect results in this region. The trick is that I do not have a good way to detect the phase to verify that they are not in the two phase region without recursing to another attempt at finding the phase boundary.

[danielkwalsh]

@msaitta-mpr @ibell It seems that a common approach used is to pre-calculate the phase boundary. See the "Vapor-Liquid Equilibrium" section of this link for more details.

[msaitta-mpr]

@ibell Looking at this again, I think that the issue actually is with the ancillary function. I went ahead and quickly remade the ancillary and things seem to be better now.

[danielkwalsh]

Hey, @msaitta-mpr, does this mean the issue is fixed?

[danielkwalsh]

Does the singular J matrix still crop up with the fixed ancillary? I'm wondering if this part of the code is still needed?

[msaitta-mpr]

I never actually checked whether it is still needed. I figure that it is a cheap enough check that it is probably worth leaving in.

[danielkwalsh]

Okay, great. I'm satisfied with that.

[danielkwalsh]

@msaitta-mpr @ibell, so what are the next steps to get this merged?

[ibell]

If everyone agrees we are ready I can merge it

[msaitta-mpr]

I believe that in its current state that it fixes the known issue and does not introduce any others. It raises the question of whether there are other ancillary issues hidden elsewhere, but I think that we should merge for now and reevaluate the ancillaries some other time.

[ibell]

@msaitta-mpr @danielkwalsh please feel free to introduce yourselves here: #2206

[danielkwalsh]

@ibell, @msaitta-mpr, thanks for merging! When do we get a new release as a new version?

[ibell]

No idea. It will depend on getting some other bugfixes/improvements together, enough that we can justify the time it takes to put together a release.