sepehr/Entropyk
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Entropyk/_bmad-output/implementation-artifacts/4-6-smart-initialization-heuristic.md

20 KiB
Raw Blame History

Story 4.6: Smart Initialization Heuristic

Status: done

Story

As a R&D engineer (Marie), I want automatic initial guesses from source and sink temperatures, so that cold start convergence is fast and reliable.

Acceptance Criteria

  1. Antoine Equation Pressure Estimation (AC: #1)

    • Given source temperature T_source (evaporator side) and sink temperature T_sink (condenser side)
    • When calling SmartInitializer::estimate_pressures(fluid, T_source, T_sink)
    • Then evaporator pressure is estimated via Antoine equation at T_source - ΔT_approach
    • And condenser pressure is estimated via Antoine equation at T_sink + ΔT_approach
    • And both pressures are returned as Pressure NewType values

  2. Evaporator Pressure Below Critical (AC: #2)

    • Given any source temperature
    • When estimating evaporator pressure
    • Then P_evap < P_critical for the fluid
    • And if estimated pressure exceeds critical, it is clamped to 0.5 * P_critical

  3. Condenser Pressure from Sink + Approach ΔT (AC: #3)

    • Given sink temperature T_sink and configurable approach temperature difference ΔT_approach
    • When estimating condenser pressure
    • Then P_cond = P_sat(T_sink + ΔT_approach)
    • And default ΔT_approach = 5.0 K

  4. State Vector Population (AC: #4)

    • Given a finalized System and estimated pressures
    • When calling SmartInitializer::populate_state(system, P_evap, P_cond, h_guess)
    • Then the state vector is filled with alternating [P, h] per edge
    • And pressure assignment follows circuit topology (evaporator edges → P_evap, condenser edges → P_cond)
    • And enthalpy is initialized to a physically reasonable default per fluid

  5. Fluid-Agnostic Antoine Coefficients (AC: #5)

    • Given common refrigerants (R134a, R410A, R32, R744/CO2, R290/Propane)
    • When estimating saturation pressure
    • Then built-in Antoine coefficients are used (no external fluid backend required)
    • And results are within 5% of CoolProp saturation pressure for T ∈ [−40°C, +80°C]

  6. Graceful Fallback for Unknown Fluids (AC: #6)

    • Given an unrecognized fluid identifier
    • When calling estimate_pressures
    • Then a sensible default state is returned (e.g., P_evap = 5 bar, P_cond = 20 bar)
    • And a tracing::warn! is emitted with the unknown fluid name

  7. No Heap Allocation in Hot Path (AC: #7)

    • Given a pre-finalized System
    • When calling populate_state
    • Then no heap allocation occurs during state vector filling
    • And the state slice is written in-place

  8. Integration with Solver (AC: #8)

    • Given a FallbackSolver (or NewtonConfig / PicardConfig)
    • When the user calls solver.with_initial_state(state_vec)
    • Then the solver starts from the provided initial state instead of zeros
    • And the solver accepts Vec<f64> as the initial state

Tasks / Subtasks

  • Create crates/solver/src/initializer.rs module (AC: #1, #2, #3, #5, #6)

    • Define SmartInitializer struct with ΔT_approach: f64 field (default 5.0 K)
    • Implement InitializerConfig struct with dt_approach: f64 and fluid: FluidId
    • Implement AntoineCoefficients struct: A, B, C fields (log10 form, Pa units)
    • Implement antoine_pressure(T_kelvin: f64, coeffs: &AntoineCoefficients) -> f64 pure function
    • Add built-in coefficient table for: R134a, R410A, R32, R744 (CO2), R290 (Propane)
    • Implement SmartInitializer::estimate_pressures(fluid, T_source_K, T_sink_K) -> (Pressure, Pressure)
    • Clamp P_evap to 0.5 * P_critical if above critical (AC: #2)
    • Add tracing::warn! for unknown fluids with fallback defaults (AC: #6)

  • Implement SmartInitializer::populate_state (AC: #4, #7)

    • Accept system: &System, P_evap: Pressure, P_cond: Pressure, h_default: Enthalpy
    • Accept mutable state: &mut [f64] slice (no allocation)
    • Iterate over system.edge_indices() and fill state[2*i] = P, state[2*i+1] = h
    • Use P_evap for edges in circuit 0 (evaporator circuit), P_cond for circuit 1 (condenser circuit)
    • For single-circuit systems, alternate between P_evap and P_cond based on edge index parity

  • Add initial_state support to solver configs (AC: #8)

    • Add initial_state: Option<Vec<f64>> field to NewtonConfig
    • Add initial_state: Option<Vec<f64>> field to PicardConfig
    • Add with_initial_state(mut self, state: Vec<f64>) -> Self builder method to both
    • In NewtonConfig::solve(): if initial_state.is_some(), copy it into state before iteration
    • In PicardConfig::solve(): same pattern
    • Add with_initial_state to FallbackSolver (delegates to both newton/picard configs)

  • Expose SmartInitializer in lib.rs (AC: #1)

    • Add pub mod initializer; to crates/solver/src/lib.rs
    • Re-export SmartInitializer, InitializerConfig, AntoineCoefficients from lib.rs

  • Unit tests in initializer.rs (AC: #1–#7)

    • Test Antoine equation gives correct P_sat for R134a at 0°C (≈ 2.93 bar)
    • Test Antoine equation gives correct P_sat for R744 at 20°C (≈ 57.3 bar)
    • Test P_evap < P_critical for all built-in fluids at T_source = −40°C
    • Test P_cond = P_sat(T_sink + 5K) for default ΔT_approach
    • Test unknown fluid returns fallback (5 bar / 20 bar) with warn log
    • Test populate_state fills state vector correctly for 2-edge system
    • Test no allocation: verify populate_state works on pre-allocated slice

  • Integration test: cold start convergence (AC: #8)

    • Build test system and use SmartInitializer to generate initial state
    • Verify FallbackSolver::with_initial_state delegates to both sub-solvers
    • Assert convergence in ≤ iterations from zero-start for same system

Dev Notes

Epic Context

Epic 4: Intelligent Solver Engine — Solve any system with < 1s guarantee, Newton-Raphson ↔ Sequential Substitution fallback.

Story Dependencies:

  • Story 4.1 (Solver Trait Abstraction) — DONE: Solver trait, SolverError, ConvergedState defined
  • Story 4.2 (Newton-Raphson) — DONE: Full Newton with line search, timeout, divergence detection
  • Story 4.3 (Sequential Substitution) — DONE: Picard with relaxation, timeout, divergence detection
  • Story 4.4 (Intelligent Fallback) — DONE: FallbackSolver with Newton↔Picard switching
  • Story 4.5 (Time-Budgeted Solving) — IN-PROGRESS: TimeoutConfig, best-state tracking, ZOH
  • Story 4.7 (Convergence Criteria) — NEXT: Sparse Jacobian for multi-circuit

FRs covered: FR42 (Automatic Initialization Heuristic — Smart Guesser)

Architecture Context

Technical Stack:

  • entropyk_core::Temperature, entropyk_core::Pressure, entropyk_core::Enthalpy — NewType wrappers (MUST use, no bare f64 in public API)
  • entropyk_components::port::FluidId — fluid identifier (already in scope from system.rs)
  • tracing — structured logging (already in solver crate)
  • thiserror — error handling (already in solver crate)
  • No external fluid backend needed for Antoine estimation (pure math)

Code Structure:

  • NEW FILE: crates/solver/src/initializer.rs — SmartInitializer, InitializerConfig, AntoineCoefficients
  • MODIFY: crates/solver/src/solver.rs — add initial_state field to NewtonConfig, PicardConfig, FallbackSolver
  • MODIFY: crates/solver/src/lib.rs — add pub mod initializer; and re-exports

Relevant Architecture Decisions:

  • NewType pattern: Never use bare f64 for physical quantities in public API [Source: architecture.md]
  • No allocation in hot path: populate_state must write to pre-allocated &mut [f64] [Source: architecture.md NFR4]
  • Zero-panic policy: All operations return Result [Source: architecture.md]
  • tracing for logging: Never println! [Source: architecture.md]
  • #![deny(warnings)]: All new code must pass clippy [Source: architecture.md]

Developer Context

Existing Infrastructure to Leverage:

// crates/solver/src/solver.rs — EXISTING (Story 4.1–4.5)

pub struct NewtonConfig {
    pub max_iterations: usize,
    pub tolerance: f64,
    pub line_search: bool,
    pub timeout: Option<Duration>,
    pub use_numerical_jacobian: bool,
    pub line_search_armijo_c: f64,
    pub line_search_max_backtracks: usize,
    pub divergence_threshold: f64,
    pub timeout_config: TimeoutConfig,
    pub previous_state: Option<Vec<f64>>,
    // ADD: pub initial_state: Option<Vec<f64>>,
}

pub struct PicardConfig {
    pub max_iterations: usize,
    pub tolerance: f64,
    pub relaxation_factor: f64,
    pub timeout: Option<Duration>,
    pub divergence_threshold: f64,
    pub divergence_patience: usize,
    pub timeout_config: TimeoutConfig,
    pub previous_state: Option<Vec<f64>>,
    // ADD: pub initial_state: Option<Vec<f64>>,
}

pub struct FallbackSolver {
    pub config: FallbackConfig,
    pub newton_config: NewtonConfig,
    pub picard_config: PicardConfig,
}

// System state vector layout (from system.rs):
// [P_edge0, h_edge0, P_edge1, h_edge1, ...]
// Each edge contributes 2 entries: pressure (Pa) then enthalpy (J/kg)

How initial_state integrates into NewtonConfig::solve():

// EXISTING in NewtonConfig::solve():
let mut state: Vec<f64> = vec![0.0; n_state];

// CHANGE TO:
let mut state: Vec<f64> = if let Some(ref init) = self.initial_state {
    debug_assert_eq!(init.len(), n_state, "initial_state length mismatch");
    init.clone()
} else {
    vec![0.0; n_state]
};

Antoine Equation (log10 form):

log10(P_sat [Pa]) = A - B / (C + T [°C])

Built-in coefficients (valid range approximately −40°C to +80°C):

Fluid A B C P_critical (Pa)
R134a 8.9300 1766.0 243.0 4,059,280
R410A 9.1200 1885.0 243.0 4,901,200
R32 9.0500 1780.0 243.0 5,782,000
R744 9.8100 1347.8 273.0 7,377,300
R290 8.8700 1656.0 243.0 4,247,200

Note: These are approximate coefficients tuned for the −40°C to +80°C range. Accuracy within 5% of NIST is sufficient for initialization purposes — the solver will converge to the exact solution regardless.

SmartInitializer API Design:

// crates/solver/src/initializer.rs

use entropyk_core::{Pressure, Temperature, Enthalpy};
use entropyk_components::port::FluidId;
use crate::system::System;

pub struct AntoineCoefficients {
    pub a: f64,
    pub b: f64,
    pub c: f64,  // °C offset
    pub p_critical_pa: f64,
}

pub struct InitializerConfig {
    pub fluid: FluidId,
    /// Temperature approach difference for condenser (K). Default: 5.0
    pub dt_approach: f64,
}

impl Default for InitializerConfig {
    fn default() -> Self {
        Self {
            fluid: FluidId::new("R134a"),
            dt_approach: 5.0,
        }
    }
}

pub struct SmartInitializer {
    pub config: InitializerConfig,
}

impl SmartInitializer {
    pub fn new(config: InitializerConfig) -> Self { ... }

    /// Estimate (P_evap, P_cond) from source and sink temperatures.
    ///
    /// Returns Err if temperatures are physically impossible (e.g., T > T_critical).
    pub fn estimate_pressures(
        &self,
        t_source: Temperature,
        t_sink: Temperature,
    ) -> Result<(Pressure, Pressure), InitializerError> { ... }

    /// Fill a pre-allocated state vector with smart initial guesses.
    ///
    /// No heap allocation. `state` must have length == system.state_vector_len().
    pub fn populate_state(
        &self,
        system: &System,
        p_evap: Pressure,
        p_cond: Pressure,
        h_default: Enthalpy,
        state: &mut [f64],
    ) { ... }
}

InitializerError enum:

#[derive(Error, Debug, Clone, PartialEq)]
pub enum InitializerError {
    #[error("Temperature {temp_celsius:.1}°C exceeds critical temperature for {fluid}")]
    TemperatureAboveCritical { temp_celsius: f64, fluid: String },

    #[error("State slice length {actual} does not match system state vector length {expected}")]
    StateLengthMismatch { expected: usize, actual: usize },
}

populate_state algorithm:

For each edge i in system.edge_indices():
    circuit = system.edge_circuit(edge_i)
    p = if circuit.0 == 0 { p_evap } else { p_cond }
    state[2*i]     = p.to_pascals()
    state[2*i + 1] = h_default.to_joules_per_kg()

For single-circuit systems (all edges in circuit 0), use p_evap for all edges. The solver will quickly adjust pressures to the correct distribution.

Architecture Compliance

  • NewType pattern: Pressure::from_pascals(), Temperature::from_kelvin(), Enthalpy::from_joules_per_kg() — ALWAYS use these, never bare f64 in public API
  • No bare f64 in SmartInitializer public methods
  • tracing: tracing::warn! for unknown fluids, tracing::debug! for estimated pressures
  • Result<T, E>: estimate_pressures returns Result<(Pressure, Pressure), InitializerError>
  • approx: Use assert_relative_eq! in tests for floating-point comparisons (tolerance 5%)
  • Pre-allocation: populate_state takes &mut [f64], never allocates

Library/Framework Requirements

  • entropyk_core — Pressure, Temperature, Enthalpy NewTypes (already a dependency of solver crate)
  • entropyk_components — FluidId from port module (already a dependency)
  • thiserror — InitializerError derive (already in solver crate)
  • tracing — Structured logging (already in solver crate)
  • approx — Test assertions (already in dev-dependencies)

File Structure Requirements

New files:

  • crates/solver/src/initializer.rs — SmartInitializer, InitializerConfig, AntoineCoefficients, InitializerError

Modified files:

  • crates/solver/src/solver.rs — Add initial_state: Option<Vec<f64>> to NewtonConfig, PicardConfig; add with_initial_state() builder; update solve() to use it
  • crates/solver/src/lib.rs — Add pub mod initializer; and re-exports

Tests:

  • Unit tests in crates/solver/src/initializer.rs (inline #[cfg(test)] module)
  • Integration test in crates/solver/tests/ (cold start convergence test)

Testing Requirements

Unit Tests (initializer.rs):

  • test_antoine_r134a_at_0c: P_sat(0°C) ≈ 2.93 bar, assert within 5% of 293,000 Pa
  • test_antoine_r744_at_20c: P_sat(20°C) ≈ 57.3 bar, assert within 5%
  • test_p_evap_below_critical: For all built-in fluids at T_source = −40°C, P_evap < P_critical
  • test_p_cond_approach: P_cond = P_sat(T_sink + 5K) for default config
  • test_unknown_fluid_fallback: Returns (5 bar, 20 bar) with no panic
  • test_populate_state_2_edges: 2-edge system, state = [P_evap, h, P_cond, h]
  • test_populate_state_no_alloc: Verify function signature takes &mut [f64]

Integration Tests:

  • test_cold_start_convergence_r134a: 4-component cycle, smart init → FallbackSolver → converges in < 20 iterations

Performance:

  • populate_state must not allocate (verified by signature: &mut [f64])
  • Antoine computation: pure arithmetic, < 1 µs per call

Previous Story Intelligence (4.5)

Key Patterns Established:

  • Pre-allocated buffers: let mut buf: Vec<f64> = vec![0.0; n]; before loop, then buf.copy_from_slice(...) inside
  • Builder pattern: fn with_X(mut self, x: X) -> Self { self.field = x; self }
  • tracing::info! for solver events, tracing::debug! for per-iteration data, tracing::warn! for degraded behavior
  • ConvergedState::new(state, iterations, final_residual, status) constructor
  • residual_norm() helper: residuals.iter().map(|r| r * r).sum::<f64>().sqrt()

initial_state Integration Pattern (follow Story 4.5 previous_state pattern):

// In NewtonConfig::solve():
// BEFORE the iteration loop, replace:
//   let mut state: Vec<f64> = vec![0.0; n_state];
// WITH:
let mut state: Vec<f64> = self.initial_state
    .as_ref()
    .map(|s| {
        debug_assert_eq!(s.len(), n_state, "initial_state length mismatch");
        s.clone()
    })
    .unwrap_or_else(|| vec![0.0; n_state]);

Git Intelligence

Recent commits:

  • be70a7a — feat(core): implement physical types with NewType pattern (Pressure, Temperature, Enthalpy, MassFlow all available)
  • Stories 4.1–4.4 complete (Solver trait, Newton, Picard, FallbackSolver)
  • Story 4.5 in-progress (TimeoutConfig, best-state tracking, ZOH)

Project Structure Notes

  • Workspace root: /Users/sepehr/dev/Entropyk
  • Solver crate: crates/solver/src/ — add initializer.rs here
  • Core types: crates/core/src/ — Pressure, Temperature, Enthalpy already defined
  • Components: crates/components/src/port.rs — FluidId already defined
  • Existing modules in solver: coupling.rs, error.rs, graph.rs, jacobian.rs, lib.rs, solver.rs, system.rs

References

  • FR42: [Source: epics.md#FR42 — "System includes Automatic Initialization Heuristic (Smart Guesser) proposing coherent initial pressure values based on source/sink temperatures"]
  • Story 4.6 AC: [Source: epics.md#Story-4.6 — "pressures estimated via Antoine equation", "evaporator pressure < P_critical", "condenser pressure from T_sink + ΔT_approach"]
  • Architecture — NewType: [Source: architecture.md — "NewType pattern for all physical quantities (Pressure, Temperature, Enthalpy, MassFlow) - never bare f64 in public APIs"]
  • Architecture — No allocation: [Source: architecture.md NFR4 — "No dynamic allocation in solver loop (pre-calculated allocation only)"]
  • Architecture — tracing: [Source: architecture.md — "tracing for structured logging (never println!)"]
  • Architecture — thiserror: [Source: architecture.md — "thiserror for error handling with ThermoError enum"]
  • State vector layout: [Source: system.rs#state_layout — "[P_edge0, h_edge0, P_edge1, h_edge1, ...] — 2 per edge (pressure Pa, enthalpy J/kg)"]
  • System::edge_circuit(): [Source: system.rs#edge_circuit — returns CircuitId for an edge based on source node]
  • NFR1: [Source: prd.md — "Steady State convergence time < 1 second for standard cycle in Cold Start"]
  • NFR4: [Source: prd.md — "No dynamic allocation in solver loop (pre-calculated allocation only)"]

Dev Agent Record

Agent Model Used

Antigravity (Code Review + Fixes)

Debug Log References

  • Antoine A coefficient bug: coefficients were calibrated for log10(P [bar]) but equation expected log10(P [Pa]). Fixed by recalculating A from NIST reference P_sat values.
  • Existing tests (newton_raphson.rs, picard_sequential.rs, newton_convergence.rs) used struct literal syntax without ..Default::default() — broke when new fields added. Fixed.

Completion Notes List

  1. All 8 ACs implemented and verified by tests.
  2. Antoine coefficient bug found and fixed during code review (A values off by ~7 orders of magnitude in Pa output).
  3. 4 pre-existing compilation errors in older test files fixed (missing struct fields from Story 4.5/4.6 additions).
  4. New integration test file smart_initializer.rs added for AC #8 coverage.
  5. Full test suite: 140+ tests pass, 0 failures.

File List

  • crates/solver/src/initializer.rs (NEW)
  • crates/solver/src/solver.rs (MODIFIED — added initial_state to NewtonConfig, PicardConfig, FallbackSolver)
  • crates/solver/src/lib.rs (MODIFIED — added pub mod initializer; and re-exports)
  • crates/solver/tests/smart_initializer.rs (NEW — AC #8 integration tests)
  • crates/solver/tests/newton_raphson.rs (MODIFIED — fixed struct literal missing fields)
  • crates/solver/tests/newton_convergence.rs (MODIFIED — fixed struct literal missing fields)
  • crates/solver/tests/picard_sequential.rs (MODIFIED — fixed struct literal missing fields x2)