SaturationSurrogate.cpp

Go to the documentation of this file.

// Implementation of CoolProp::sbtl::SaturationSurrogate.
//
// Lays out 11 cubic-spline curves over a log-uniform T knot mesh on
// [T_triple, 0.9999 * T_crit]:
//
//   * P_of_T_     (one curve — same for L/V by Maxwell equal-pressure)
//   * D_of_T_[2]  (liquid + vapor density)
//   * H_of_T_[2]
//   * S_of_T_[2]
//   * U_of_T_[2]
//   * T_of_logp_  (inverse: T as a function of log(p_sat))
//
// At construction we sweep the source backend via QT_INPUTS at each
// knot (Q=0 then Q=1) and record (P, D, H, S, U).  The 11 curves are
// then built independently from these column-extracted vectors.
//
// All splines are CoolProp::region::CubicSplineCurve — natural cubic,
// O(1) amortized eval via the indexed-search machinery — already used
// for boundary curves elsewhere in the codebase, so no new numerics
// here.  Per-property tight bounds (b_min, b_max) are computed at
// build time but currently unused (kept for parity with the existing
// CubicSplineCurve API; future range-check assertions can rely on
// them).
 
#include "CoolProp/sbtl/SaturationSurrogate.h"
 
#include <array>
#include <cmath>
#include <memory>
#include <stdexcept>
#include <utility>
#include <vector>
 
#include "CoolProp/AbstractState.h"
#include "CoolProp/region/CubicSplineCurve.h"
#include "CoolProp/DataStructures.h"
 
namespace CoolProp {
namespace sbtl {
 
namespace {
 
// Knot mesh: Chebyshev-Lobatto distribution in [T_lo, T_hi] —
//
//   T_k = (T_lo + T_hi)/2 - (T_hi - T_lo)/2 * cos(pi * k / (N-1))
//
// concentrates knots at both endpoints, in particular near T_crit
// where the sat-density curve bends most sharply.  Empirically (Water
// at 64 knots): max relative error in saturated density drops from
// ~1.5e-6 (log-T spacing) to <1e-7, well inside the bead's 1e-6
// acceptance target.  Pressure / enthalpy / entropy / energy all
// already passed at 1e-6 with log-T spacing — the rho curve is the
// most curvature-concentrated one and is the budget binder here.
//
// (The bead suggested log-T spacing but the acceptance bar is on
// accuracy, not spacing law; this deviation is intentional.)
std::vector<double> chebyshev_lobatto_T_knots(double T_lo, double T_hi, std::size_t n) {
    // Callers are required to ensure n >= 2 (and build_from_source enforces
    // n >= kMinKnots == 4); the n==1 case would divide by zero below.
    std::vector<double> T(n);
    const double mid = 0.5 * (T_lo + T_hi);
    const double half = 0.5 * (T_hi - T_lo);
    const auto denom = static_cast<double>(n - 1);
    for (std::size_t k = 0; k < n; ++k) {
        const double theta = M_PI * static_cast<double>(k) / denom;
        T[k] = mid - half * std::cos(theta);
    }
    // Numerical safety: cos(0) == 1 and cos(pi) == -1 exactly in
    // IEEE-754, so endpoint values should already be T_lo / T_hi;
    // assign explicitly anyway to be robust against any future
    // change in the spacing formula.
    T.front() = T_lo;
    T.back() = T_hi;
    return T;
}
 
constexpr std::size_t kMinKnots = 4;
constexpr double kCriticalAvoidance = 0.9999;
 
}  // namespace
 
struct SaturationSurrogate::Impl
{
    std::unique_ptr<region::CubicSplineCurve> P_of_T;  // L/V identical (Maxwell)
    std::array<std::unique_ptr<region::CubicSplineCurve>, 2> D_of_T;
    std::array<std::unique_ptr<region::CubicSplineCurve>, 2> H_of_T;
    std::array<std::unique_ptr<region::CubicSplineCurve>, 2> S_of_T;
    std::array<std::unique_ptr<region::CubicSplineCurve>, 2> U_of_T;
    std::unique_ptr<region::CubicSplineCurve> T_of_logp;
    double T_lo = 0.0;
    double T_hi = 0.0;
    double p_lo = 0.0;
    double p_hi = 0.0;
    std::size_t n_knots = 0;
};
 
SaturationSurrogate::SaturationSurrogate() : impl_(std::make_unique<Impl>()) {}
SaturationSurrogate::~SaturationSurrogate() = default;
 
std::unique_ptr<SaturationSurrogate> SaturationSurrogate::build_from_source(::CoolProp::AbstractState& src, std::size_t n_knots) {
    if (n_knots < kMinKnots) {
        return nullptr;
    }
    double T_triple = 0.0;
    double T_crit = 0.0;
    try {
        T_triple = src.Ttriple();
        T_crit = src.T_critical();
    } catch (const std::exception&) {
        return nullptr;
    }
    // Triple point sometimes comes back below the source's actual
    // sat-curve floor (e.g. REFPROP reports a slightly different value
    // from its internal flash range) — bump up a hair so the first
    // QT_INPUTS probe doesn't reject.
    const double T_lo = T_triple * (1.0 + 1.0e-6);
    const double T_hi = T_crit * kCriticalAvoidance;
    if (!(T_lo > 0.0) || !(T_hi > T_lo)) {
        return nullptr;
    }
 
    const auto T_knots = chebyshev_lobatto_T_knots(T_lo, T_hi, n_knots);
 
    std::vector<double> p_col(n_knots);
    std::array<std::vector<double>, 2> D_col = {std::vector<double>(n_knots), std::vector<double>(n_knots)};
    std::array<std::vector<double>, 2> H_col = {std::vector<double>(n_knots), std::vector<double>(n_knots)};
    std::array<std::vector<double>, 2> S_col = {std::vector<double>(n_knots), std::vector<double>(n_knots)};
    std::array<std::vector<double>, 2> U_col = {std::vector<double>(n_knots), std::vector<double>(n_knots)};
 
    for (std::size_t k = 0; k < n_knots; ++k) {
        const double T = T_knots[k];
        double p_k = 0.0;
        for (int side : {0, 1}) {
            try {
                src.update(::CoolProp::QT_INPUTS, static_cast<double>(side), T);
            } catch (const std::exception&) {
                return nullptr;
            }
            if (side == 0) {
                p_k = src.p();
            }
            D_col[side][k] = src.rhomolar();
            H_col[side][k] = src.hmolar();
            S_col[side][k] = src.smolar();
            U_col[side][k] = src.umolar();
            // region::CubicSplineCurve::build doesn't reject non-finite
            // dependent values — a NaN/inf from the source would land
            // silently in the spline coefficients and corrupt every
            // eval downstream.  Bail on the first bad sample.
            if (!std::isfinite(D_col[side][k]) || !std::isfinite(H_col[side][k]) || !std::isfinite(S_col[side][k])
                || !std::isfinite(U_col[side][k])) {
                return nullptr;
            }
        }
        p_col[k] = p_k;
        // p_sat must be strictly increasing in T inside the dome; if
        // the source returned a non-monotone or non-finite sample, bail
        // — log(p) for the inverse curve will fail otherwise.
        if (!std::isfinite(p_k) || p_k <= 0.0) {
            return nullptr;
        }
        if (k > 0 && !(p_k > p_col[k - 1])) {
            return nullptr;
        }
    }
 
    auto surrogate = std::unique_ptr<SaturationSurrogate>(new SaturationSurrogate());
    auto& impl = *surrogate->impl_;
    try {
        impl.P_of_T = region::CubicSplineCurve::build(T_knots, p_col);
        for (int side : {0, 1}) {
            impl.D_of_T[side] = region::CubicSplineCurve::build(T_knots, D_col[side]);
            impl.H_of_T[side] = region::CubicSplineCurve::build(T_knots, H_col[side]);
            impl.S_of_T[side] = region::CubicSplineCurve::build(T_knots, S_col[side]);
            impl.U_of_T[side] = region::CubicSplineCurve::build(T_knots, U_col[side]);
        }
        std::vector<double> logp_col(n_knots);
        for (std::size_t k = 0; k < n_knots; ++k) {
            logp_col[k] = std::log(p_col[k]);
        }
        // T_knots is reused verbatim as the dependent column for the
        // inverse T(log p) spline — pass a copy because build() takes
        // ownership.
        impl.T_of_logp = region::CubicSplineCurve::build(logp_col, T_knots);
    } catch (const std::exception&) {
        return nullptr;
    }
    impl.T_lo = T_lo;
    impl.T_hi = T_hi;
    impl.p_lo = p_col.front();
    impl.p_hi = p_col.back();
    impl.n_knots = n_knots;
    return surrogate;
}
 
double SaturationSurrogate::get_T_from_p(double p) const {
    if (!(p > 0.0)) {
        throw std::out_of_range("SaturationSurrogate::get_T_from_p: p must be positive");
    }
    if (p < impl_->p_lo || p > impl_->p_hi) {
        throw std::out_of_range("SaturationSurrogate::get_T_from_p: p outside surrogate range");
    }
    return impl_->T_of_logp->eval(std::log(p));
}
 
double SaturationSurrogate::eval_sat(double T, char what, int side) const {
    if (side != 0 && side != 1) {
        throw std::invalid_argument("SaturationSurrogate::eval_sat: side must be 0 (liquid) or 1 (vapor)");
    }
    // NaN compares false against both bounds, so the range check below
    // silently lets non-finite T through into the spline evaluator.
    if (!std::isfinite(T)) {
        throw std::out_of_range("SaturationSurrogate::eval_sat: T must be finite");
    }
    if (T < impl_->T_lo || T > impl_->T_hi) {
        throw std::out_of_range("SaturationSurrogate::eval_sat: T outside surrogate range");
    }
    switch (what) {
        case 'P':
            return impl_->P_of_T->eval(T);
        case 'D':
            return impl_->D_of_T[side]->eval(T);
        case 'H':
            return impl_->H_of_T[side]->eval(T);
        case 'S':
            return impl_->S_of_T[side]->eval(T);
        case 'U':
            return impl_->U_of_T[side]->eval(T);
        default:
            throw std::invalid_argument("SaturationSurrogate::eval_sat: 'what' must be one of {P, D, H, S, U}");
    }
}
 
double SaturationSurrogate::T_min() const noexcept {
    return impl_->T_lo;
}
double SaturationSurrogate::T_max() const noexcept {
    return impl_->T_hi;
}
std::size_t SaturationSurrogate::n_knots() const noexcept {
    return impl_->n_knots;
}
 
}  // namespace sbtl
}  // namespace CoolProp