DAE‑safe bounding for CoolProp property calls in Julia

General Usage
, ,
1

I am building a large DAE-based vapor‑compression refrigeration loop in Julia using DifferentialEquations.jl and CoolProp for thermodynamic properties, and I am running into repeated CoolProp crashes during DAE initialization and Newton iterations. Although the physical solution is valid, the DAE solver temporarily explores nonphysical trial states (negative pressures, NaNs, or extreme enthalpies), which causes CoolProp property calls such as “(P,H)” or saturation queries “(P,Q)” to fail. I have been adding local “effective” projections using isfinite checks, pressure floors, and enthalpy bounds to keep property evaluations in a valid domain, but the challenge is ensuring these guards are applied consistently everywhere without overwriting the DAE state variables or breaking solver convergence. I’m looking for best practices or established patterns for making CoolProp calls DAE‑safe in Julia while preserving correct physics and solver robustness.

3 Likes
ModelingToolkit.jl: Handling Output/Input Dependency in a DAE Model
2

I’ve had great success adding a line search to the nonlinear solver, which backtracks (halves the proposed step) if it lands inside an infeasible domain. I’ve done this in my own implementations, but something similar should be possible using DifferentialEquations.jl.

For nominal integration, they have the option isoutofdomain and the other domain callbacks documented here: Frequently Asked Questions · DifferentialEquations.jl.

For initialization, I imagine you would have to change the nonlinear solver, possibly adding a line search. However, it is not transparent how to modify the DAE initialization solver in DifferentialEquations.jl.

1 Like
3

It has a multi-line search + trust region strategy by default. You can just pass any NonlinearSolve.jl algorithm into the struct for choosing a solver.

4

The TrustRegion NL solver does work when i reduce my DAE to steady state model. But when i tried to solve the same steady state model using IDA, the solver failed to initialize even when i fed in the solution generated from NL solver as the initial condition. I am thinking this might be related to the coolprop functions i called. Does anyone have any suggestions for formulating and solving large scale physical related DAE problems?

5

IDA isn’t very stable for this stuff. Biggest improvement is usually to move to mass matrix form. Other is to make sure the property models throw NaN when out of bounds.

6

Thank you for the reply. Mass matrix form is definitely a good suggestion and I am working on it. But the second suggestion is somewhat confusing to me. Could you please explain this more in detail?

7

One of the issues that can happen pretty often for this kind of thing is that for example f(x) can give an error if you’re outside the boundary. This then effects the solver because an error is unrecoverable. If you instead get f(x) => NaN, then you get something like du[i] = NaN which triggers a step rejection, i.e. a dt reduction, automatically, which can then cause the DAE solver to automatically recover. I don’t know if that’s what you’re hitting here because I don’t have code to run, but I’m guessing “running into repeated CoolProp crashes” means you’re seeing errors when x is outside the boundaries, so I’m mentioning that if you just check for whether the input is valid first and compute a NaN then that would cause the solving process to automatically have ways to adapt. You’d still need a decent enough initial guess so that the Newton does not start its initialization at a NaN, but most steps should be recoverable.

8

if you are using CoolProp.PropsSI, a non-erroring version would be the following:

using CoolProp_jll #compiled binary of CoolProp, CoolProp.jl is a wrapper of this binary
#modified from https://github.com/CoolProp/CoolProp.jl/blob/46a2496e77c673ea4782927eacb23364fa92e25d/src/CoolProp.jl#L80
function PropsSI_no_error(output::AbstractString, name1::AbstractString, value1::Real, name2::AbstractString, value2::Real, fluid::AbstractString)
    val = ccall( (:PropsSI, libcoolprop), Cdouble, (Cstring, Cstring, Cdouble, Cstring, Cdouble, Cstring), output, name1, value1, name2, value2, fluid)
    if val == Inf
        #error("CoolProp: ", get_global_param_string("errstring"))
        return NaN
    end
    return val
end