alice_physics/

creep_longterm.rs

//! Long-Term Creep Prediction (Findley + Time-Temperature Superposition)
//!
//! Extends the Norton
//! power-law creep in `plastic.rs` with the **Findley** three-parameter
//! model that captures both the initial elastic strain and the transient
//! logarithmic-like creep phase common to polymers over months / years.
//!
//! # Findley model
//!
//! `ε(t) = ε_0 + m · tⁿ`
//!
//! - `ε_0` = instantaneous elastic strain (dimensionless)
// LIMITATION(COV-STRUCT-064): `n` is an **integer** here
//! - `m`, `n` = fitting constants; `n` is an **integer** here (`n_int`, 1 to 4,
//!   evaluated by repeated multiplication), unlike the fractional 0.15-0.30
//!   exponents of the literature
//!
//! The PLA preset is `ε_0 = 0.003`, `m = 8.3e-14 /h³`, `n = 3`: after 6 months
//! (4380 h) it gives `0.003 + 8.3e-14 · 4380³ ≈ 1.0 %` total strain, matching
//! the ALICE-Bamboo docs "PLA 棚 半年で反る" incident. The fit matches at 6 months
//! only (see [`FindleyParameters::pla_25c_moderate`]).
//!
//! # Time-Temperature Superposition (WLF)
//!
//! To predict lifetime at temperature `T` from data at reference `T_ref`,
//! shift the time axis by the Williams-Landel-Ferry factor `a_T`:
//!
//! `log a_T = −C_1·(T − T_ref) / (C_2 + T − T_ref)`
//!
//! `t_effective(T) = t / a_T`
//!
//! Universal WLF: `C_1 = 17.44`, `C_2 = 51.6` at `T_ref = T_g`.
//!
//! # References
//!
//! - Findley, Lai, Onaran, *Creep and Relaxation of Nonlinear Viscoelastic
//!   Materials*, Dover 1989. Ch 4-5.
//! - Williams, Landel, Ferry, "The temperature dependence of relaxation
//!   mechanisms in amorphous polymers", J. Am. Chem. Soc. 77(14), 1955.
//! - Bellehumeur (2004) — polymer creep parameters for FDM prints.
//!
//! # Integration status
//!
//! `FindleyParameters`, `FindleyParameters::pla_25c_moderate`, and
//! `predict_strain` are wired into `structural_solver.rs`. The WLF
//! subsystem (`WlfConstants`, `wlf_shift_factor`, `effective_time_at_temp`,
//! `CREEP_FROZEN_AT`) and the standalone `strain_at` are crate-internal and
//! reached through `predict_strain`; `CREEP_FROZEN_AT` is the saturation
//! sentinel `effective_time_at_temp` turns into "no effective time" (see its
//! doc for the threshold). `petg_25c_moderate` has no consumer (the solver
//! installs the PLA preset for PLA only and leaves creep unmodelled for
//! other materials; this PETG fit is uncalibrated) and carries an `ALLOW-UNWIRED` debt marker;
//! its oracle is a unit test in this module.

use crate::filament_db::MaterialProperties;
use crate::math::Fix128;
use crate::math_util::{exp_fix, EXP_OVERFLOW_SENTINEL};

/// `base^n` via repeated [`Fix128::checked_mul`], or `None` the moment any
/// intermediate product would overflow Fix128.
fn pow_checked(base: Fix128, n: u32) -> Option<Fix128> {
    let mut acc = Fix128::ONE;
    for _ in 0..n {
        acc = acc.checked_mul(base)?;
    }
    Some(acc)
}

// ============================================================================
// Findley model
// ============================================================================

/// Findley three-parameter creep model parameters.
///
/// `ε(t) = ε_0 + m · tⁿ` with `t` in hours.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct FindleyParameters {
    /// Instantaneous elastic strain (dimensionless).
    pub epsilon_0: Fix128,
    /// Transient coefficient `m` (dimensionless, calibrated per hour^n).
    pub m: Fix128,
    /// Time exponent `n` (integer for Fix128 pow: 1, 2, 3, or 4).
    /// n=1 → linear creep, n=2 → parabolic, n=3 → cubic. An integer stand-in for the literature's fractional 0.15-0.30; the presets use 3.
    pub n_int: u32,
}

impl FindleyParameters {
    /// PLA at 25 °C under moderate stress (Bellehumeur 2004 fit):
    /// `ε₀ = 0.3 %`, and `m` chosen so the creep term reaches 0.7 %
    /// (total ≈ 1 %) at 6 months (4380 h): `m = 0.007 / 4380³ ≈ 8.3e-14`.
    ///
    /// **Warning**: integer `n = 3` is a much steeper time law than the
    /// literature's 0.25 — the fit matches at 6 months only, under-predicts
    /// before and over-predicts after (≈ 6 % total at 1 year). Before 1.2.0
    /// `m = 1e-12` gave 8.7 % at 6 months while the docstring claimed 1 %.
    #[must_use]
    pub fn pla_25c_moderate() -> Self {
        Self {
            epsilon_0: Fix128::from_ratio(3, 1000),
            m: Fix128::from_ratio(83, 1_000_000_000_000_000_i64),
            n_int: 3,
        }
    }

    /// PETG at 25 °C: `ε₀ = 0.2 %`, `m = 1e-13 /h³`, `n = 3`
    /// (`ε(4380 h) = 0.2 % + 1e-13 · 4380³ ≈ 1.04 %`).
    ///
    /// The elastic term is lower than PLA's (higher `T_g`), but `m` was not
    /// recalibrated when the PLA preset dropped to `8.3e-14` in 1.2.0, so
    /// this preset creeps *more* than `pla_25c_moderate` after
    /// `t = (0.001 / 1.7e-14)^(1/3) ≈ 3887 h`. Treat it as a placeholder
    /// fit until a PETG datasheet calibration replaces `m`.
    // ALLOW-DEAD: pub(crate) with no crate caller, same debt as the ALLOW-UNWIRED marker below
    // ALLOW-UNWIRED: wiring debt structural-pub-crate-residue (structural_solver::new installs the PLA preset for PLA only; this PETG fit is uncalibrated), oracle src/creep_longterm.rs tests::petg_preset_findley_matches_closed_form_at_six_months
    #[allow(dead_code)]
    #[must_use]
    pub(crate) fn petg_25c_moderate() -> Self {
        Self {
            epsilon_0: Fix128::from_ratio(2, 1000),
            m: Fix128::from_ratio(1, 10_000_000_000_000_i64),
            n_int: 3,
        }
    }

    /// Predicted total strain at time `t_hours` (dimensionless, crate-internal helper called by `predict_strain`).
    #[must_use]
    pub(crate) fn strain_at(&self, t_hours: Fix128) -> Fix128 {
        if t_hours <= Fix128::ZERO {
            return self.epsilon_0;
        }
        // 1) t^n by repeated multiplication, then times `m` — exact and
        // bit-identical to the pre-1.x behavior whenever it does not
        // overflow (the oracle tolerance is as tight as 5e-8 in that
        // range, so this path must stay untouched).
        if let Some(t_n) = pow_checked(t_hours, self.n_int) {
            if let Some(creep) = self.m.checked_mul(t_n) {
                return self.epsilon_0 + creep;
            }
        }
        // 2) `t_hours^n` overflows Fix128 well before `m * t_hours^n` (the
        // actual creep contribution) does, because the calibrated `m` is
        // many orders of magnitude below 1 (~8.3e-14 for the n=3 PLA
        // preset) — inside the documented WLF domain T_g..T_g+100 the
        // effective time can already be large enough to trip this.
        // Scaling `t_hours` by `m^(1/n)` first folds that smallness in
        // *before* the exponentiation: `(t_hours * m^(1/n))^n == m *
        // t_hours^n` exactly in real arithmetic, so the intermediate does
        // not leave the representable range as readily as the bare power
        // does.
        if self.n_int > 0 && self.m > Fix128::ZERO {
            let root = self
                .m
                .powf_pos(Fix128::ONE / Fix128::from_int(i64::from(self.n_int)));
            if let Some(scaled_t) = t_hours.checked_mul(root) {
                if let Some(creep) = pow_checked(scaled_t, self.n_int) {
                    return self.epsilon_0 + creep;
                }
            }
        }
        // 3) The effective time is astronomically large even after
        // rescaling (T far enough above T_g that the WLF shift factor is
        // minuscule) — the true creep strain genuinely exceeds Fix128's
        // representable range. Saturate instead of wrapping: the creep
        // term is always non-negative for t_hours > 0, so clamping at the
        // largest representable magnitude keeps the result monotone
        // instead of flipping sign. Returning the bare maximum (rather
        // than adding `epsilon_0` on top of it) avoids wrapping a second
        // time at the very ceiling it is meant to saturate at, and is
        // still >= every value branch 1/2 can produce (those are only
        // reached when `epsilon_0 + creep` already fits under that
        // ceiling).
        Fix128::from_raw(i64::MAX, u64::MAX)
    }
}

// ============================================================================
// Time-Temperature Superposition
// ============================================================================

/// WLF constants (crate-internal).
#[derive(Clone, Copy, Debug)]
pub(crate) struct WlfConstants {
    /// C_1 (dimensionless).
    pub(crate) c1: Fix128,
    /// C_2 (°C).
    pub(crate) c2: Fix128,
}

impl WlfConstants {
    /// Universal WLF (`C_1 = 17.44`, `C_2 = 51.6`), applied at `T_ref = T_g` (crate-internal).
    #[must_use]
    pub(crate) fn universal() -> Self {
        Self {
            c1: Fix128::from_ratio(1744, 100),
            c2: Fix128::from_ratio(516, 10),
        }
    }
}

/// Sentinel re-export of [`crate::math_util::EXP_OVERFLOW_SENTINEL`], used to
/// signal "creep frozen" (temperature well below reference, crate-internal).
pub(crate) const CREEP_FROZEN_AT: Fix128 = EXP_OVERFLOW_SENTINEL;

/// WLF shift factor `a_T` (dimensionless multiplier). For `T > T_ref` returns
/// a value < 1 (creep accelerates); for `T < T_ref` returns a value > 1.
///
/// Uses the shared [`crate::math_util::exp_fix`] deterministic exponential.
/// Very cold conditions saturate at `CREEP_FROZEN_AT` (crate-internal helper).
#[must_use]
pub(crate) fn wlf_shift_factor(temp_c: Fix128, t_ref_c: Fix128, wlf: &WlfConstants) -> Fix128 {
    let dt = temp_c - t_ref_c;
    let denom = wlf.c2 + dt;
    if denom.is_zero() {
        return Fix128::ONE;
    }
    // log10 a_T = -C_1 · dT / (C_2 + dT)
    let log_at = Fix128::ZERO - wlf.c1 * dt / denom;
    let ln10 = Fix128::from_ratio(230_258_509, 100_000_000);
    let y = log_at * ln10;
    exp_fix(y)
}

/// Effective time (hours) at temperature `t` relative to reference `t_ref`,
/// using the WLF shift factor.
///
/// When the shift factor saturates at [`CREEP_FROZEN_AT`] (the `exp_fix`
/// overflow sentinel, reached once
/// `ln10 · C₁ · (T_ref − T) / (C₂ − (T_ref − T)) ≥ 20`, i.e.
/// `T ≤ T_ref − 20·C₂ / (20 + C₁·ln10) ≈ T_ref − 17.16 °C` for the universal
/// constants) creep is frozen and the effective time is exactly zero rather
/// than `t / sentinel`.
#[must_use]
pub(crate) fn effective_time_at_temp(
    t_hours: Fix128,
    temp_c: Fix128,
    t_ref_c: Fix128,
    wlf: &WlfConstants,
) -> Fix128 {
    let dt = temp_c - t_ref_c;
    if dt > Fix128::ZERO {
        // above the reference a_T < 1 can fall below the 2^-64 resolution
        // (a_T = 4e-18 at T_ref + 13200 for the universal constants): divide by
        // it as t * exp(-ln a_T) instead, -ln a_T = ln10 C1 dT / (C2 + dT) being
        // below ln10 C1 = 40.2; `exp_fix` covers arguments below 20, so a
        // larger one is taken as the cube of exp(x / 3)
        let ln10 = Fix128::from_ratio(230_258_509, 100_000_000);
        let neg_ln_at = ln10 * wlf.c1 * dt / (wlf.c2 + dt);
        let speedup = if neg_ln_at < Fix128::from_int(20) {
            exp_fix(neg_ln_at)
        } else {
            let third = exp_fix(neg_ln_at / Fix128::from_int(3));
            third * third * third
        };
        return t_hours
            .checked_mul(speedup)
            .unwrap_or(Fix128::from_raw(i64::MAX, u64::MAX));
    }
    let a_t = wlf_shift_factor(temp_c, t_ref_c, wlf);
    if a_t.is_zero() {
        return Fix128::ZERO;
    }
    if a_t >= CREEP_FROZEN_AT {
        return Fix128::ZERO;
    }
    t_hours / a_t
}

// ============================================================================
// Aggregate helper
// ============================================================================

/// Predict long-term strain under constant stress and temperature.
///
/// Combines `FindleyParameters` with a WLF time–temperature shift **above
/// `T_g`** (the domain the universal constants are fitted for, Williams,
/// Landel & Ferry 1955: `T_g ≤ T ≤ T_g + 100 °C`), referenced to `T_g`. At
/// or below `T_g` the parameters are used as calibrated: the crate presets
/// (`pla_25c_moderate`, …) are room-temperature fits, and the glassy-regime
/// temperature dependence follows Arrhenius with an activation energy the
/// parameter set does not carry (a 2.0 field). Before 1.2.0 the WLF shift
/// was applied *below* `T_g` too, where the universal constants give
/// `a_T ≈ 10^37` at `T_g − 35 °C` — the effective time collapsed to zero and
/// room-temperature PLA reported no creep at all
/// (`tests/engineering_oracles_solid.rs`).
#[must_use]
pub fn predict_strain(
    parameters: &FindleyParameters,
    material: &MaterialProperties,
    t_hours: Fix128,
    operating_temp_c: Fix128,
) -> Fix128 {
    let t_g = material.glass_transition_c;
    if t_g.is_zero() || operating_temp_c <= t_g {
        return parameters.strain_at(t_hours);
    }
    let wlf = WlfConstants::universal();
    let t_eff = effective_time_at_temp(t_hours, operating_temp_c, t_g, &wlf);
    parameters.strain_at(t_eff)
}

// ============================================================================
// Tests
// ============================================================================

#[cfg(test)]
mod tests {
    use super::*;

    fn approx_eq(a: Fix128, b: Fix128, tol: Fix128) -> bool {
        let d = if a > b { a - b } else { b - a };
        d <= tol
    }

    #[test]
    fn strain_at_zero_time_is_epsilon_0() {
        let p = FindleyParameters::pla_25c_moderate();
        assert_eq!(p.strain_at(Fix128::ZERO), p.epsilon_0);
    }

    #[test]
    fn strain_monotonic_increasing() {
        let p = FindleyParameters::pla_25c_moderate();
        let e1 = p.strain_at(Fix128::from_int(100));
        let e2 = p.strain_at(Fix128::from_int(1000));
        assert!(e2 > e1);
    }

    #[test]
    fn petg_lower_creep_than_pla() {
        let pla = FindleyParameters::pla_25c_moderate();
        let petg = FindleyParameters::petg_25c_moderate();
        let t = Fix128::from_int(1000);
        assert!(petg.strain_at(t) < pla.strain_at(t));
    }

    #[test]
    fn wlf_shift_at_ref_temp_is_unity() {
        let wlf = WlfConstants::universal();
        let a = wlf_shift_factor(Fix128::from_int(60), Fix128::from_int(60), &wlf);
        assert!(approx_eq(a, Fix128::ONE, Fix128::from_ratio(1, 1000)));
    }

    #[test]
    fn wlf_shift_above_ref_less_than_unity() {
        let wlf = WlfConstants::universal();
        // T > Tref → creep accelerates → a_T < 1 → effective time > actual
        let a = wlf_shift_factor(Fix128::from_int(70), Fix128::from_int(60), &wlf);
        assert!(a < Fix128::ONE);
    }

    #[test]
    fn wlf_shift_below_ref_greater_than_unity() {
        let wlf = WlfConstants::universal();
        // T < Tref → creep slows → a_T > 1
        let a = wlf_shift_factor(Fix128::from_int(50), Fix128::from_int(60), &wlf);
        assert!(a > Fix128::ONE);
    }

    #[test]
    fn effective_time_zero_a_t_returns_zero() {
        // Manually construct pathological WLF whose numerator forces log→−∞
        // Skip this edge case here — normal usage never triggers.
        let wlf = WlfConstants::universal();
        let t_eff = effective_time_at_temp(
            Fix128::from_int(100),
            Fix128::from_int(60),
            Fix128::from_int(60),
            &wlf,
        );
        assert!(approx_eq(
            t_eff,
            Fix128::from_int(100),
            Fix128::from_ratio(1, 100)
        ));
    }

    #[test]
    fn predict_strain_scales_with_time_near_tg() {
        // At 55°C (close to PLA Tg=60), the WLF shift factor is close to
        // unity and creep is significant. At 25°C (35°C below Tg) creep is
        // essentially frozen — the model correctly returns near-zero
        // additional strain.
        let params = FindleyParameters::pla_25c_moderate();
        let material = MaterialProperties::pla();
        let e_short = predict_strain(
            &params,
            &material,
            Fix128::from_int(1000),
            Fix128::from_int(55),
        );
        let e_long = predict_strain(
            &params,
            &material,
            Fix128::from_int(10_000),
            Fix128::from_int(55),
        );
        assert!(e_long > e_short);
    }

    #[test]
    fn predict_strain_above_tg_much_higher() {
        let params = FindleyParameters::pla_25c_moderate();
        let material = MaterialProperties::pla();
        let e_room = predict_strain(
            &params,
            &material,
            Fix128::from_int(1000),
            Fix128::from_int(25),
        );
        let e_hot = predict_strain(
            &params,
            &material,
            Fix128::from_int(1000),
            Fix128::from_int(65),
        );
        // At 65°C (above Tg 60), creep should be worse
        assert!(e_hot > e_room);
    }

    /// Oracle: Findley `ε(t) = ε₀ + m·tⁿ` with the PETG preset's constants
    /// (`ε₀ = 0.002`, `m = 1e-13`, `n = 3`) at 6 months, evaluated by hand:
    /// `4380³ = 84 027 672 000`, `m·t³ = 0.0084027672`, total
    /// `0.0104027672`. Routed through `predict_strain` at 25 °C, below the
    /// PETG `T_g = 80 °C`, so no WLF shift applies. Tolerance: `m = 1e-13`
    /// is `1 844 674.4` Fix128 ulp (`2⁻⁶⁴`), stored truncated, so the creep
    /// term carries a relative error of `0.4 / 1.8e6 = 2.2e-7`
    /// (`1.9e-9` absolute); `1e-8` leaves no room for a wrong constant.
    #[test]
    fn petg_preset_findley_matches_closed_form_at_six_months() {
        let got = predict_strain(
            &FindleyParameters::petg_25c_moderate(),
            &MaterialProperties::petg(),
            Fix128::from_int(4380),
            Fix128::from_int(25),
        )
        .to_f64();
        let want = 0.002 + 1e-13 * 84_027_672_000.0;
        assert!((got - want).abs() < 1e-8, "got {got}, want {want}");
    }

    /// Degenerate time for the PETG preset: `strain_at` documents
    /// `t ≤ 0 → ε₀`, so zero and negative hours return exactly `0.002`.
    #[test]
    fn petg_preset_degenerate_time_returns_epsilon_0() {
        let p = FindleyParameters::petg_25c_moderate();
        let m = MaterialProperties::petg();
        let t25 = Fix128::from_int(25);
        assert_eq!(
            predict_strain(&p, &m, Fix128::ZERO, t25),
            Fix128::from_ratio(2, 1000)
        );
        assert_eq!(
            predict_strain(&p, &m, Fix128::from_int(-1), t25),
            Fix128::from_ratio(2, 1000)
        );
    }

    /// Oracle: the frozen threshold derived by hand from the WLF form and the
    /// `exp_fix` saturation point. `a_T = exp(ln10 · C₁ · |ΔT| / (C₂ − |ΔT|))`
    /// saturates once the exponent reaches 20, i.e.
    /// `|ΔT| ≥ 20·C₂ / (20 + C₁·ln10) = 1032 / 60.158 = 17.155 °C`.
    /// At `ΔT = −17.5 °C` the exponent is `2.302585·17.44·17.5/34.1 = 20.61`
    /// (frozen, effective time exactly 0); at `ΔT = −17 °C` it is
    /// `2.302585·17.44·17/34.6 = 19.729` (not frozen,
    /// `t_eff = t · exp(−19.729)`).
    #[test]
    fn effective_time_is_exactly_zero_once_the_wlf_shift_saturates() {
        let wlf = WlfConstants::universal();
        let t = Fix128::from_int(4380);
        let t_ref = Fix128::from_int(60);
        let frozen = effective_time_at_temp(t, t_ref - Fix128::from_ratio(35, 2), t_ref, &wlf);
        assert_eq!(frozen, Fix128::ZERO, "ΔT = −17.5 °C must be frozen");
        let not_frozen = effective_time_at_temp(t, t_ref - Fix128::from_int(17), t_ref, &wlf);
        let exponent = core::f64::consts::LN_10 * 17.44 * 17.0 / (51.6 - 17.0);
        let want = 4380.0 * crate::det_math::exp64(-exponent);
        let got = not_frozen.to_f64();
        assert!(
            got > 0.0 && ((got - want) / want).abs() < 1e-6,
            "ΔT = −17 °C: got {got:e}, want {want:e}"
        );
    }

    #[test]
    fn universal_wlf_constants() {
        let w = WlfConstants::universal();
        assert!(approx_eq(
            w.c1,
            Fix128::from_ratio(1744, 100),
            Fix128::from_ratio(1, 100)
        ));
        assert!(approx_eq(
            w.c2,
            Fix128::from_ratio(516, 10),
            Fix128::from_ratio(1, 100)
        ));
    }
}