Prevent crashes near critical density due to saturation calc

[msaitta-mpr]

Description of the Change

When updating an HEOS abstract state from density and energy, crashes (ValueError thrown) can occur when the density is very near the critical density. This PR fixes this behavior and prevents the exception.

HSU_D_flash is called when updating with D and U (as well as other combinations). As one of the first steps in the flash routine, the saturation points are determined at the given density (saturation_D_pure). However, when near the critical density, instabilities in the solver routines cause an exception. This PR utilizes a few strategies to prevent this:

• When inverting the ancillary curve, if the density is outside of the ancillary curve's domain, extrapolate in the secant solver.
• When using the 2D Newton solver within saturation_D_pure, prevent the calculation of a negative temperature or density by dynamically changing omega (underrelaxation factor). I.e., check what the new values would be and if they were going to be negative, set the omega value to instead cause a bisection.
• Wrap the calls to saturation_D_pure in a try block and repeat the solution with a smaller starting omega value (and more allowed iterations) a few times in an attempt to get even better stability. Cases that were already solving fine will require no further computation.

The "hole" in the results is illustrated in the attached image (colormap of temperature).

Benefits

Currently, CoolProp will fail for many fluids when updating from density/energy when the density is near the critical density. Updates from density/energy are a fairly common method in thermal hydraulic code development and crossing the isochor of the critical density is not uncommon. Prior to this PR, CoolProp would throw an exception, usually leading to a crash of the user's program. This allows the solver to continue without impacting calculations at other parts of the domain.

Possible Drawbacks

I foresee two potential drawbacks:

• Because of the underrelaxation factor, certain flash calculation could now take significantly (1000x or more) longer. However, this should only occur for cases which previous would have thrown an exception. In most cases, a slower solution is better than a faster failure.
• Extrapolating the ancillary curves may cause poor guesses for further solution, but a poor guess is better than no guess.

Verification Process

I evaluated several thousand density/energy combinations near the critical point for both hydrogen and helium (see c++ code below). In almost all cases, rather than throwing an exception, a good value was created. A few points fail in the Brent solver. These failures are due to a two-phase state where the temperature is within 0.01 of the critical temperature because HSU_D_flash_twophase applies a 0.01 bracket limit to the temperature. I have not removed this for fear of causing bracketing issues in other parts of the code.

Furthermore, when linked to another code, the code produces good, expected results - both near and removed from the critical point. This indicates that there is likely no regression failure.

Applicable Issues

Closes #2154

Test code

/*
This program demonstrates how the topic-2154 PR fixes failures of HSU_D
(update from DU) when density is near the critical density for certain fluids.
*/
#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;
}

void test_range(std::shared_ptr<CoolProp::AbstractState> state)
{
    double density = state->rhomass_critical();
    double tcrit = state->T_critical();
    state->update(CoolProp::DmassT_INPUTS, density, tcrit);
    double energy = state->umass();
    state->update(CoolProp::DmassUmass_INPUTS, density, energy);
    std::vector<double> density_vals = linspace(density*0.75, density*1.5, 5000);
    std::vector<double> energy_vals = linspace(energy*0.5, energy*10., 200);
    for (double ee: energy_vals){
        for (double d : density_vals){
            try{
                state->update(CoolProp::DmassUmass_INPUTS, d, ee);
            }
            catch(CoolProp::ValueError e){
                std::cout << "BAD rho=" << d << " e=" << ee << " error=" << e.what() << std::endl;
            }
            catch(CoolProp::CoolPropBaseError e){
                std::cout << "BAD rho=" << d << " e=" << ee << " error=" << e.what() << std::endl;
            }
        }
    }
}

int main()
{
    std::cout << " === Test Helium ===" << std::endl;
    std::shared_ptr<CoolProp::AbstractState> helium = std::shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "Helium"));
    test_range(helium);

    std::cout << " === Test Hydrogen ===" << std::endl;
    std::shared_ptr<CoolProp::AbstractState> hydrogen = std::shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "Hydrogen"));
    test_range(hydrogen);

    std::cout << " == Test was successful if you saw no 'BAD' output" << "\n";
    std::cout << " == Brent failures are due to 2-phase temperature within 0.01K of critical (HSU_D_flash_twophase)" << std::endl;
    return EXIT_SUCCESS;
}

Sample Plot

[CLAassistant]

All committers have signed the CLA.

[ibell]

This is amazing!! I have known of this problem for a long time now, but didn't have the bandwidth to do anything about it. Can you also double-check that you have fixed DQ inputs? I added some editorial comments on your PR.

[ibell]

Overall, looks good, I'd like to keep the function signature the same as much as possible.

[ibell]

I think the max_arguments should be a field in the options class rather than a new argument. That's why I made that options class in the first place.

[ibell]

Empty comment, please remove

[ibell]

Same comment about the max_iter argument

[msaitta-mpr]

I was unaware that DQ inputs were a problem, so I had not looked into it. At your prompting, I took a look and saw that DQ inputs fail in the Brent solver (not bracketing) in many places - see the attached image. This PR does not fix that issue (same results with and without the modification), and I suspect that they are unrelated.

I can see that this could be an issue and I can try to look into it when I get some more time. Is there an existing issue number that I should track against?

[ibell]

What does your picture look like after your fix? If complete, I'll merge. Overall your fixes seem reasonable from my looking at it.

[msaitta-mpr]

My fix does not resolve this, as the failure for Dmass_Q inputs is unrelated. I will provide some discussion below, but suggest that you merge this change and that the DQ issue be solved in a separate PR.

DmassQ fails for densities greater than the critical density. This is illustrated in the plot below where the red points indicate failures to converge. All of these failures indicate that the initial guess of the Brent solver does not bracket. Note that in this plot the consecutive dots form lines of constant quality.

Recall that the Brent solver that is used here iterates on temperature and seeks the error in the quality to be reduced to 0. The initial guesses for the Brent solver are the triple point temperature (lower bound) and the critical temperature (upper bound). Initially, these seem like good guesses because a two-phase fluid should obviously be within this temperature range. However, at these temperatures, we do not bracket the desired quality. Instead, we find that both temperatures will give a quality that is too low! The triple point temperature evaluation is reasonable (typically giving a quality in the range [0, 1e-3]), but the critical point is not. The critical point quality is reported as a negative number. We can plot the change in quality with respect to temperature to investigate further. We see that there is a sudden reversal and that the estimated quality goes negative once the density/temperature combination goes outside of the two-phase region.

My initial guess of how to correct this is to add another solver loop to provide an upper limit of the temperature guess to be the saturation temperature for the given density, but I have not yet implemented this.

Red line indicates the true quality that we are seeking.

[ibell]

Ok, let's merge this and then if you want to dig into the DQ that would be great. The DQ routines in REFPROP also fail sometimes.