Entropyk/_bmad-output/implementation-artifacts/4-6-smart-initialization-heuristic.md

# Story 4.6: Smart Initialization Heuristic

Status: done

<!-- Note: Validation is optional. Run validate-create-story for quality check before dev-story. -->

## 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

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

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

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

- [x] Expose `SmartInitializer` in `lib.rs` (AC: #1)
  - [x] Add `pub mod initializer;` to `crates/solver/src/lib.rs`
  - [x] Re-export `SmartInitializer`, `InitializerConfig`, `AntoineCoefficients` from `lib.rs`

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

- [x] Integration test: cold start convergence (AC: #8)
  - [x] Build test system and use `SmartInitializer` to generate initial state
  - [x] Verify `FallbackSolver::with_initial_state` delegates to both sub-solvers
  - [x] 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:**

```rust
// 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()`:**

```rust
// 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:**

```rust
// 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:**

```rust
#[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):**
```rust
// 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)