diff --git a/.entropyk/modules/calibration.ekmod b/.entropyk/modules/calibration.ekmod new file mode 100644 index 0000000..a1264cb --- /dev/null +++ b/.entropyk/modules/calibration.ekmod @@ -0,0 +1,222 @@ +{ + "components": [ + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "comp", + "params": { + "displacement_m3": 0.000065, + "emergent_pressure": true, + "isentropic_efficiency": 0.7, + "speed_hz": 50, + "superheat_k": 5, + "t_cond_k": 318.15, + "t_evap_k": 278.15, + "volumetric_efficiency": 0.92 + }, + "rotation": 0, + "type": "IsentropicCompressor" + }, + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "cond", + "params": { + "__fixed_calib_sdt_c": false, + "__fixed_z_dp": true, + "__fixed_z_ua": true, + "calib_sdt_c": 40, + "channel_spacing_mm": 1.5, + "chevron_angle_deg": 60, + "correlation": "Longo2004", + "corrugation_pitch_mm": 3, + "dp_correlation": "SimplifiedChannel", + "emergent_pressure": true, + "manufacturer": "", + "model": "", + "n_plates": 40, + "plate_length_m": 0.4, + "plate_thickness_mm": 0.6, + "plate_width_m": 0.12, + "refrigerant": "R134a", + "secondary_fluid": "Water", + "target_subcooling_k": 5, + "ua": 2500, + "z_dp": 1, + "z_ua": 1 + }, + "rotation": 0, + "type": "BphxCondenser" + }, + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "exv", + "params": { + "__fixed_opening": true, + "emergent_pressure": true, + "opening": 1, + "t_evap_k": 278.15 + }, + "rotation": 0, + "type": "IsenthalpicExpansionValve" + }, + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "evap", + "params": { + "__fixed_calib_sst_c": false, + "__fixed_z_dp": true, + "__fixed_z_ua": false, + "calib_sst_c": 5, + "channel_spacing_mm": 1.5, + "chevron_angle_deg": 60, + "correlation": "Longo2004", + "corrugation_pitch_mm": 3, + "dp_correlation": "SimplifiedChannel", + "emergent_pressure": true, + "manufacturer": "", + "model": "", + "n_plates": 40, + "plate_length_m": 0.4, + "plate_thickness_mm": 0.6, + "plate_width_m": 0.12, + "refrigerant": "R134a", + "secondary_fluid": "Water", + "target_superheat_k": 5, + "ua": 2000, + "z_dp": 1, + "z_ua": 1 + }, + "rotation": 0, + "type": "BphxEvaporator" + }, + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "cond_water_in", + "params": { + "__fixed_m_flow_kg_s": true, + "__fixed_p_set_bar": false, + "__fixed_t_set_c": true, + "fluid": "Water", + "m_flow_kg_s": 0.4, + "p_set_bar": 2, + "t_set_c": 30 + }, + "rotation": 0, + "type": "BrineSource" + }, + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "cond_water_out", + "params": { + "__fixed_p_back_bar": true, + "__fixed_t_set_c": false, + "fluid": "Water", + "p_back_bar": 2, + "t_set_c": 7 + }, + "rotation": 0, + "type": "BrineSink" + }, + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "evap_water_in", + "params": { + "__fixed_m_flow_kg_s": true, + "__fixed_p_set_bar": false, + "__fixed_t_set_c": true, + "fluid": "Water", + "m_flow_kg_s": 0.5, + "p_set_bar": 3, + "t_set_c": 12 + }, + "rotation": 0, + "type": "BrineSource" + }, + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "evap_water_out", + "params": { + "__fixed_p_back_bar": true, + "__fixed_t_set_c": false, + "fluid": "Water", + "p_back_bar": 3, + "t_set_c": 7 + }, + "rotation": 0, + "type": "BrineSink" + }, + { + "circuit": 0, + "flipH": false, + "flipV": false, + "name": "SST probe", + "params": { + "__fixed_target": true, + "fluid": "R134a", + "measure": "SST", + "target": 5.9 + }, + "rotation": 0, + "type": "Probe" + } + ], + "connections": [], + "edges": [ + { + "from": "comp:outlet", + "to": "cond:inlet" + }, + { + "from": "cond:outlet", + "to": "exv:inlet" + }, + { + "from": "exv:outlet", + "to": "evap:inlet" + }, + { + "from": "cond_water_in:outlet", + "to": "cond:secondary_inlet" + }, + { + "from": "cond:secondary_outlet", + "to": "cond_water_out:inlet" + }, + { + "from": "evap_water_in:outlet", + "to": "evap:secondary_inlet" + }, + { + "from": "evap:secondary_outlet", + "to": "evap_water_out:inlet" + }, + { + "from": "evap:outlet", + "to": "SST probe:inlet" + }, + { + "from": "SST probe:outlet", + "to": "comp:inlet" + } + ], + "instances": [], + "params": {}, + "ports": { + "inlet": "comp:inlet" + } +} \ No newline at end of file diff --git a/_bmad-output/implementation-artifacts/sprint-status.yaml b/_bmad-output/implementation-artifacts/sprint-status.yaml index 924ed32..fb3dc1b 100644 --- a/_bmad-output/implementation-artifacts/sprint-status.yaml +++ b/_bmad-output/implementation-artifacts/sprint-status.yaml @@ -26,14 +26,14 @@ # - Dev moves story to 'review', then Dev runs code-review (fresh context, ideally different LLM) generated: 2026-07-18T12:51:31Z -last_updated: 2026-07-19T14:40:00Z +last_updated: 2026-07-19T16:00:00Z project: Entropyk -project_key: ENTROPYK +project_key: ENTROPYKnon tracking_system: file-system story_location: "{project-root}/_bmad-output/implementation-artifacts" development_status: - epic-0: done + epic-0: in-progress 0-1-fix-failing-tests: done epic-0-retrospective: optional # --- Epic 0: Physical Regularization & Jacobian Health (epics-solver-hp.md) --- @@ -42,8 +42,8 @@ development_status: 0-3-valve-flow-exv-regularization: done 0-4-heat-exchanger-regularization: done 0-5-domain-sweep-jacobian-health-gate: done - 0-6-emergent-pressure-seed-robustness: review 0-7-bphx-property-call-reduction: review + 0-8-bphx-phase-change-entu-duty: review # --- Epic 1: Performance Foundation - Linear Algebra Modernization --- epic-1: in-progress 1-1-phase-0-baseline-golden-safety-net: review diff --git a/apps/web/src/components/canvas/Canvas.tsx b/apps/web/src/components/canvas/Canvas.tsx index c750491..e457c68 100644 --- a/apps/web/src/components/canvas/Canvas.tsx +++ b/apps/web/src/components/canvas/Canvas.tsx @@ -29,7 +29,7 @@ import { } from "@/lib/mediaStyle"; import { findNearestEdge, - isPipeType, + isInsertableOnEdge, pipeMediaKind, } from "@/lib/edgeInsert"; import { defaultParams } from "@/lib/componentMeta"; @@ -154,17 +154,24 @@ export default function Canvas({ onOpenModuleEditor }: CanvasProps) { return; } - // Pipe / duct dropped on a wire → splice into the circuit (Modelica-style). - if (isPipeType(type)) { - const prefer = pipeMediaKind(type, defaultParams(type)); + // Pipe / duct / Probe dropped on a wire → splice into the circuit + // (Modelica-style). Probe is a zero-residual measurement tap that uses + // the same insertOnEdge path. + if (isInsertableOnEdge(type)) { + // A Probe can tap ANY medium (refrigerant / water / air) — a sensor sits + // wherever the physics is, so we don't filter by media for it. Pipes + // still get media filtering so a water pipe doesn't splice into a + // refrigerant line. + const prefer = type === "Probe" ? undefined : pipeMediaKind(type, defaultParams(type)); const hit = findNearestEdge(position, nodes, edges, { - maxDistance: 64, + // Probes are small targets — be more forgiving on the drop distance. + maxDistance: type === "Probe" ? 96 : 64, preferMedia: prefer, }); if (hit) { const result = insertOnEdge(type, hit.midpoint, hit.edge.id); if (result.ok) return; - // Media mismatch: fall through to free placement. + // Media mismatch / splice failed: fall through to free placement. } } diff --git a/apps/web/src/components/canvas/ComponentIcon.tsx b/apps/web/src/components/canvas/ComponentIcon.tsx index 0aa3472..af17d91 100644 --- a/apps/web/src/components/canvas/ComponentIcon.tsx +++ b/apps/web/src/components/canvas/ComponentIcon.tsx @@ -306,6 +306,22 @@ const GLYPHS: Record ReactElement> = { ), + + // ── Probe — small measurement tap on a connection line ───────── + // A filled circle with a thin ring + a small dot inside. Deliberately + // compact (the 36×36 node footprint makes it read as a tap, not a block). + // The four short ticks hint at a sensor crosshair. + probe: (c) => ( + + + + + + + ), }; function glyphKey(type?: string, opts?: HxVisualOpts): keyof typeof GLYPHS { @@ -314,6 +330,7 @@ function glyphKey(type?: string, opts?: HxVisualOpts): keyof typeof GLYPHS { if (hx && hx in GLYPHS) return hx as keyof typeof GLYPHS; if (t === "SaturatedController") return "control"; + if (t === "Probe") return "probe"; if (t.includes("Compressor")) return "compressor"; if (t === "HeatExchanger") return "hx"; if (t === "ReversingValve") return "reversing"; diff --git a/apps/web/src/components/panels/PropertiesPanel.tsx b/apps/web/src/components/panels/PropertiesPanel.tsx index 962c17c..e882b6a 100644 --- a/apps/web/src/components/panels/PropertiesPanel.tsx +++ b/apps/web/src/components/panels/PropertiesPanel.tsx @@ -427,17 +427,25 @@ export default function PropertiesPanel() { : isParamFixed(params, p); const free = Boolean(p.fixable && !fixed); const displayVal = effectiveParamValue(params, p); + // Probe: label the target row with the chosen measurement kind + // (SST / SDT / SH / SC / …) instead of the generic "Cible". + const label = + node.data.type === "Probe" && p.key === "target" + ? (COMPONENT_BY_TYPE["Probe"]?.params + .find((mp) => mp.key === "measure") + ?.options?.find((o) => o.value === params.measure)?.label ?? p.label) + : p.label; return ( - {p.label} + {label} {p.required ? * : null} @@ -543,9 +551,9 @@ export default function PropertiesPanel() { {activeTab === "Calibration" && (

- Comment calibrer : coche Fixed sur la cible (SST/SDT), décoche - Fixed sur le facteur (Z_UA). Z_UA vaut 1 par défaut (pas de - correction). + Calibration : pose une Probe (sonde) sur la + ligne à mesurer, coche Fixed sur sa cible, puis décoche Fixed sur le facteur + (Z_UA) ici. La valeur résolue s’affiche ci-dessous après simulation.

Tu n’as pas besoin du nœud « Regulation loop » pour ça. diff --git a/apps/web/src/lib/componentMeta.ts b/apps/web/src/lib/componentMeta.ts index 72b6688..abe05eb 100644 --- a/apps/web/src/lib/componentMeta.ts +++ b/apps/web/src/lib/componentMeta.ts @@ -158,6 +158,7 @@ export function portRole(port: string): "source" | "target" { */ export function nodeSize(type: string): { w: number; h: number } { const t = String(type || ""); + if (t === "Probe") return { w: 36, h: 36 }; if (t === "SaturatedController") return { w: 184, h: 96 }; if (t.includes("Compressor")) return { w: 76, h: 76 }; if (t === "HeatExchanger") return { w: 92, h: 92 }; @@ -207,11 +208,80 @@ export const COMPONENT_CATEGORIES = [ "Piping", "Hydraulics", "Boundaries", + "Instrumentation", "Advanced", "Generic", ] as const; export const COMPONENTS: ComponentMeta[] = [ + // ── Instrumentation: measurement taps for calibration (Probe = all measurements) ── + { + type: "Probe", + label: "Probe (sonde)", + category: "Instrumentation", + description: + "Sonde de mesure posée sur une ligne. Lit P/T/SH/SC/SST/SDT/DGT/DSH/capacité/ṁ/h à cet endroit. Toute calibration impose sa cible sur un Probe.", + help: + "Sonde de calibration (Probe)\n" + + "Toute calibration se fait via un Probe posé sur une ligne (règle dure).\n\n" + + "1. Glisse le Probe sur un fil du schéma (comme pour insérer un Pipe).\n" + + "2. Choisis la grandeur mesurée (SST, SDT, DGT, DSH, SH, SC, T, P, débit, capacité, enthalpie) et le fluide.\n" + + "3. Coche Fixed sur « Cible » et tape la valeur mesurée.\n" + + "4. Libère le Z-factor correspondant (Z_UA, Z_dP, Z_flow…) sur le composant à calibrer.\n" + + "Le solveur appaire automatiquement le Probe (mesure) et le Z-factor libre (inconnue).", + ports: ["inlet", "outlet"], + color: "#0ea5e9", + params: [ + { + key: "measure", + label: "Grandeur mesurée", + kind: "string", + default: "SH", + required: true, + section: "Mesure", + description: + "Grandeur physique mesurée à l'emplacement du Probe. SST/SDT = Tsat(P), SH/DSH = T−Tsat(P), SC = Tsat(P)−T.", + options: [ + { value: "SST", label: "SST (Tsat aspiration)" }, + { value: "SDT", label: "SDT (Tsat refoulement)" }, + { value: "DGT", label: "DGT (T° gaz refoulement)" }, + { value: "DSH", label: "DSH (surch. refoulement)" }, + { value: "SH", label: "SH (surchauffe)" }, + { value: "SC", label: "SC (sous-refroid.)" }, + { value: "T", label: "T (température)" }, + { value: "P", label: "P (pression)" }, + { value: "MassFlow", label: "Débit massique" }, + { value: "Capacity", label: "Capacité (kW)" }, + { value: "Enthalpy", label: "Enthalpie" }, + ], + }, + { + key: "fluid", + label: "Fluide", + kind: "string", + default: "R134a", + required: true, + section: "Mesure", + description: "Fluide frigorigène pour le calcul de Tsat/SH/SC (via CoolProp).", + }, + { + key: "target", + label: "Cible (Fixed = impose)", + kind: "number", + section: "Calibration", + fixable: true, + defaultFixed: false, + // Auto: resolved from `measure` at config-build time + // (SST/SDT → saturationTemperature, SH/DSH → superheat, SC → subcooling, + // DGT/T → temperature, P → pressure, MassFlow → massFlowRate, + // Capacity → capacity, Enthalpy → heatTransferRate as carrier). + measureOutput: "auto", + description: + "Coche Fixed et tape la valeur mesurée (K pour températures, K pour SH/SC, Pa pour P, kg/s pour débit, W pour capacité, J/kg pour h).", + }, + ], + }, + // ── Advanced regulation (optional — most users use Fixed on Calibration tab) ── { type: "SaturatedController", @@ -739,19 +809,6 @@ export const COMPONENTS: ComponentMeta[] = [ { key: "fan_head_pressure_target_c", label: "Fan head-pressure target", kind: "number", unit: "°C", section: "Controls", advanced: true }, { key: "flooded_head_pressure_target_c", label: "Flooded head-pressure target", kind: "number", unit: "°C", section: "Controls", advanced: true }, { key: "skip_pressure_eq", label: "Skip pressure equation", kind: "boolean", default: false, section: "Solver structure", advanced: true }, - // Fixed checkbox: ON = impose SDT; OFF = ignore as measure - { - key: "calib_sdt_c", - label: "SDT (calib target)", - kind: "number", - unit: "°C", - default: 40.0, - section: "Calibration", - fixable: true, - defaultFixed: false, - measureOutput: "saturationTemperature", - description: "Cocher Fixed pour imposer la SDT (mesure). Libérer Z_UA en même temps.", - }, { key: "z_ua", label: "Z_UA (UA factor)", @@ -803,18 +860,6 @@ export const COMPONENTS: ComponentMeta[] = [ { key: "ua", label: "UA", kind: "number", unit: "W/K", default: 6000.0, required: true, section: "Thermal model", min: 0.0, sweepable: true, sweepGroup: "thermal" }, { key: "t_sat_k", label: "Fixed saturation temperature", kind: "number", unit: "K", default: 275.15, section: "Thermal model", description: "Used only in fixed-pressure mode." }, { key: "superheat_k", label: "Outlet superheat closure", kind: "number", unit: "K", default: 5.0, section: "Thermal model", min: 0.0 }, - { - key: "calib_sst_c", - label: "SST (calib target)", - kind: "number", - unit: "°C", - default: 5.0, - section: "Calibration", - fixable: true, - defaultFixed: false, - measureOutput: "saturationTemperature", - description: "Cocher Fixed pour imposer la SST. Décocher Fixed sur Z_UA pour calibrer.", - }, { key: "z_ua", label: "Z_UA (UA factor)", @@ -1090,17 +1135,6 @@ export const COMPONENTS: ComponentMeta[] = [ ], }, { key: "ua", label: "UA override", kind: "number", unit: "W/K", section: "Calibration", min: 0.0, advanced: true }, - { - key: "calib_sst_c", - label: "SST (calib target)", - kind: "number", - unit: "°C", - default: 5.0, - section: "Calibration", - fixable: true, - defaultFixed: false, - measureOutput: "saturationTemperature", - }, { key: "z_ua", label: "Z_UA (UA factor)", @@ -1190,17 +1224,6 @@ export const COMPONENTS: ComponentMeta[] = [ ], }, { key: "ua", label: "UA override", kind: "number", unit: "W/K", section: "Calibration", min: 0.0, advanced: true }, - { - key: "calib_sdt_c", - label: "SDT (calib target)", - kind: "number", - unit: "°C", - default: 40.0, - section: "Calibration", - fixable: true, - defaultFixed: false, - measureOutput: "saturationTemperature", - }, { key: "z_ua", label: "Z_UA (UA factor)", @@ -1270,18 +1293,6 @@ export const COMPONENTS: ComponentMeta[] = [ section: "Secondary hydraulics", min: 0.0, }, - { - key: "calib_sst_c", - label: "SST (calib target)", - kind: "number", - unit: "°C", - default: 5.0, - section: "Calibration", - fixable: true, - defaultFixed: false, - measureOutput: "saturationTemperature", - description: "Cocher Fixed pour imposer la SST mesurée. Décocher Fixed sur Z_UA.", - }, { key: "z_ua", label: "Z_UA (UA factor)", diff --git a/apps/web/src/lib/configBuilder.test.ts b/apps/web/src/lib/configBuilder.test.ts index 7759928..9d7c9b9 100644 --- a/apps/web/src/lib/configBuilder.test.ts +++ b/apps/web/src/lib/configBuilder.test.ts @@ -483,44 +483,68 @@ describe("caloporteur (secondary stream) resolution", () => { }); }); -describe("Fixed / Free calibration (Dymola-style checkboxes)", () => { - it("emits a control when SST is Fixed and Z_UA is Free", () => { +describe("Fixed / Free calibration (Probe-based pairing)", () => { + it("pairs a Fixed Probe measure with a Free Z_UA on the HX", () => { const nodes = [ node("e", "FloodedEvaporator", "evap", 0, { ua: 9000, - calib_sst_c: 5, z_ua: 1.0, - // Fixed ON for SST, Fixed OFF for Z_UA - [fixedFlagKey("calib_sst_c")]: true, + // Fixed OFF for Z_UA → free actuator [fixedFlagKey("z_ua")]: false, }), + node("p1", "Probe", "sst_probe", 0, { + measure: "SST", + fluid: "R134a", + target: 278.15, + // Fixed ON for the Probe target → impose the measure + [fixedFlagKey("target")]: true, + }), ]; - const controls = buildFixedFreeCalibrationControls(nodes); + const edges: Edge[] = [ + { id: "pe", source: "e", target: "p1", sourceHandle: "outlet", targetHandle: "inlet" }, + ]; + const controls = buildFixedFreeCalibrationControls(nodes, edges); expect(controls).toHaveLength(1); expect(controls[0]).toMatchObject({ - measure: { component: "evap", output: "saturationTemperature" }, + measure: { component: "sst_probe", output: "saturationTemperature" }, actuator: { component: "evap", factor: "z_ua" }, - target: 5 + 273.15, + target: 278.15, }); - const cfg = buildScenarioConfig(nodes, []); + const cfg = buildScenarioConfig(nodes, edges); expect(cfg.controls?.length).toBe(1); - // UI-only keys stripped from component JSON + // UI-only keys stripped from component JSON; Z_UA still emitted as initial value. const comp = cfg.circuits[0].components[0]; - expect(comp.calib_sst_c).toBeUndefined(); expect(comp[fixedFlagKey("z_ua")]).toBeUndefined(); expect(comp.z_ua).toBe(1.0); }); - it("emits no control when Z_UA stays Fixed", () => { + it("emits no control when Z_UA stays Fixed even with a Probe present", () => { const nodes = [ node("e", "Evaporator", "evap", 0, { ua: 6000, - calib_sst_c: 5, z_ua: 1.0, - [fixedFlagKey("calib_sst_c")]: true, [fixedFlagKey("z_ua")]: true, }), + node("p1", "Probe", "sst_probe", 0, { + measure: "SST", + fluid: "R134a", + target: 278.15, + [fixedFlagKey("target")]: true, + }), + ]; + expect(buildFixedFreeCalibrationControls(nodes)).toHaveLength(0); + }); + + it("emits nothing for a freed Z_UA with no matching Probe (HARD RULE)", () => { + // A freed z-factor without a Probe to measure against must not emit a + // control — the user must place a Probe. No legacy same-component fallback. + const nodes = [ + node("e", "Evaporator", "evap", 0, { + ua: 6000, + z_ua: 1.0, + [fixedFlagKey("z_ua")]: false, + }), ]; expect(buildFixedFreeCalibrationControls(nodes)).toHaveLength(0); }); diff --git a/apps/web/src/lib/configBuilder.ts b/apps/web/src/lib/configBuilder.ts index 32b18ac..c0465eb 100644 --- a/apps/web/src/lib/configBuilder.ts +++ b/apps/web/src/lib/configBuilder.ts @@ -336,7 +336,7 @@ export function buildScenarioConfig( const nodeControls = nodes .filter((node) => node.data.type === CONTROL_NODE_TYPE) .map(controlNodeToConfig); - const fixedFreeControls = buildFixedFreeCalibrationControls(nodes); + const fixedFreeControls = buildFixedFreeCalibrationControls(nodes, edges); const controls = mergeControls( mergeControls(options.controls ?? [], nodeControls), fixedFreeControls, @@ -362,7 +362,7 @@ export function buildScenarioConfig( /** * Remove UI-only keys before sending to the CLI: * - `__fixed_*` Fixed checkbox flags - * - measure-only calib targets (calib_sst_c, …) — they become control setpoints + * - measure-only targets (e.g. Probe `target`) — they become control setpoints * - emit Modelica-style `fix_pressure` / `fix_temperature` / `fix_mass_flow` */ export function stripUiOnlyParams( @@ -466,44 +466,85 @@ export function applyBoundaryFixSemantics( } /** - * Build inverse-calibration controls from Fixed checkboxes (EES/Dymola style). + * Build the `controls[]` array for calibration (Probe-based model, HARD RULE): * - * - Param with `measureOutput` + Fixed ON → impose that measure (setpoint = value) - * - Param with `actuatorFactor` + Fixed OFF → free that Z-factor - * Pairs on the same component: each free factor with the first fixed measure. + * - Probe with Fixed `target` → impose that measure (setpoint = value) + * - Param with `actuatorFactor` + Fixed OFF → free that Z-factor + * + * Every freed Z-factor is paired with a semantically-compatible Probe + * (preferring one adjacent on an edge). A freed factor with no matching + * Probe emits nothing — the user must place a Probe to calibrate. + * + * Pairing matrix (factor → compatible Probe measures): + * z_ua → SST, SDT, SH, DSH, SC, Capacity + * z_dp → P + * z_flow → MassFlow, Capacity + * z_power → Capacity + * z_etav → MassFlow */ export function buildFixedFreeCalibrationControls( nodes: Node[], + edges: Edge[] = [], ): ControlConfig[] { const controls: ControlConfig[] = []; + // ── Pass 1: collect Probe measures and per-component freed z-factors ── + type ProbeMeasure = { + nodeName: string; + kind: string; // SST, SDT, DGT, DSH, SH, SC, T, P, MassFlow, Capacity, Enthalpy + target: number; + factorCompat: readonly string[]; // z-factors this Probe kind can pair with + }; + type FreeAct = { + factor: string; + initial: number; + min: number; + max: number; + key: string; + }; + + const probeMeasures: ProbeMeasure[] = []; + const freeActsByComponent = new Map< + string, + { nodeName: string; node: Node; acts: FreeAct[] } + >(); + for (const node of nodes) { if (node.data.type === CONTROL_NODE_TYPE) continue; const meta = COMPONENT_BY_TYPE[node.data.type]; if (!meta) continue; const params = node.data.params ?? {}; - type Measure = { output: string; target: number; key: string }; - type FreeAct = { factor: string; initial: number; min: number; max: number; key: string }; + // Probe node: read the `target` Fixed param + the configured `measure` kind. + if (node.data.type === "Probe") { + const measureParam = meta.params.find((p) => p.key === "measure"); + const targetParam = meta.params.find((p) => p.key === "target"); + if (!measureParam || !targetParam) continue; + const kindRaw = params.measure; + const kind = + typeof kindRaw === "string" ? kindRaw : String(measureParam.default ?? "SH"); + const targetMeta = targetParam; + const targetFixed = isParamFixed(params, targetParam); + if (!targetFixed) continue; + const rawTarget = params.target; + const n = typeof rawTarget === "number" ? rawTarget : Number(rawTarget); + if (!Number.isFinite(n)) continue; + probeMeasures.push({ + nodeName: node.data.name, + kind, + target: measureSetpointSi(targetMeta, n), + factorCompat: FACTOR_COMPATIBILITY[kind] ?? [], + }); + continue; + } - const measures: Measure[] = []; + // Non-Probe node: collect freed z-factors (actuators to calibrate). const freeActs: FreeAct[] = []; - for (const p of meta.params) { if (!p.fixable) continue; const fixed = isParamFixed(params, p); const raw = params[p.key]; - if (p.measureOutput && fixed) { - const n = typeof raw === "number" ? raw : Number(raw); - if (!Number.isFinite(n)) continue; - measures.push({ - output: p.measureOutput, - target: measureSetpointSi(p, n), - key: p.key, - }); - } - if (p.actuatorFactor && !fixed) { const n = typeof raw === "number" ? raw : Number(raw); const initial = Number.isFinite(n) ? n : 1.0; @@ -517,37 +558,144 @@ export function buildFixedFreeCalibrationControls( } } - if (freeActs.length === 0 || measures.length === 0) continue; + if (freeActs.length > 0) { + freeActsByComponent.set(node.data.name, { nodeName: node.data.name, node, acts: freeActs }); + } + } - for (let i = 0; i < freeActs.length; i++) { - const act = freeActs[i]; - const meas = measures[Math.min(i, measures.length - 1)]; - controls.push({ - type: "SaturatedController", - id: `calib_${node.data.name}_${act.factor}`, - measure: { - component: node.data.name, - output: meas.output, - }, - actuator: { - component: node.data.name, - factor: act.factor, - initial: act.initial, - min: act.min, - max: act.max, - }, - target: meas.target, - // Negative gain: higher Z_UA → higher capacity / often higher T_sat for flooded - // — user can refine; default chosen for UA-style calibration. - gain: -0.5, - band: 2.0, - }); + // ── Adjacency index: for each component node, which Probe nodes share an edge ── + const adjacentProbes = new Map>(); + const probeNodeIds = new Set( + nodes.filter((n) => n.data.type === "Probe").map((n) => n.id), + ); + for (const edge of edges) { + const endpoints = [edge.source, edge.target]; + for (const endpoint of endpoints) { + if (!endpoint) continue; + // Probe directly adjacent to this endpoint + for (const other of endpoints) { + if (other && other !== endpoint && probeNodeIds.has(other)) { + const probeName = nodes.find((n) => n.id === other)?.data.name; + if (probeName) { + const endpointName = nodes.find((n) => n.id === endpoint)?.data.name; + if (endpointName) { + if (!adjacentProbes.has(endpointName)) { + adjacentProbes.set(endpointName, new Set()); + } + adjacentProbes.get(endpointName)!.add(probeName); + } + } + } + } + } + } + + /** + * Find the best Probe to pair with a freed z-factor on `component`: + * 1. Must be semantically compatible (FACTOR_COMPATIBILITY). + * 2. Prefer Probes adjacent to the component (share an edge endpoint). + * 3. Fall back to any compatible Probe. + */ + const findProbeFor = ( + component: string, + factor: string, + usedProbeNames: Set, + ): ProbeMeasure | undefined => { + const compatible = probeMeasures.filter( + (m) => m.factorCompat.includes(factor) && !usedProbeNames.has(m.nodeName), + ); + if (compatible.length === 0) return undefined; + const adjacent = adjacentProbes.get(component); + const adjacentMatch = adjacent + ? compatible.find((m) => adjacent.has(m.nodeName)) + : undefined; + const chosen = adjacentMatch ?? compatible[0]; + usedProbeNames.add(chosen.nodeName); + return chosen; + }; + + // ── Pass 2: emit controls ── + // Each freed z-factor is paired with a compatible Probe (HARD RULE). A + // freed actuator with no matching Probe emits nothing — the user must + // place a Probe to calibrate. + const usedProbeNames = new Set(); + + for (const { nodeName, acts } of freeActsByComponent.values()) { + for (const act of acts) { + const probe = findProbeFor(nodeName, act.factor, usedProbeNames); + if (probe) { + controls.push({ + type: "SaturatedController", + id: `calib_${nodeName}_${act.factor}`, + measure: { + component: probe.nodeName, + output: probeOutputFor(probe.kind), + }, + actuator: { + component: nodeName, + factor: act.factor, + initial: act.initial, + min: act.min, + max: act.max, + }, + target: probe.target, + gain: -0.5, + band: 2.0, + }); + } } } return controls; } +/** + * Maps a Probe `measure` kind to the solver-side `ComponentOutput` string + * accepted by the CLI (`parse_component_output` in run.rs). + */ +function probeOutputFor(kind: string): string { + switch (kind) { + case "SST": + case "SDT": + return "saturationTemperature"; + case "SH": + case "DSH": + return "superheat"; + case "SC": + return "subcooling"; + case "DGT": + case "T": + return "temperature"; + case "P": + return "pressure"; + case "MassFlow": + return "massFlowRate"; + case "Capacity": + return "capacity"; + case "Enthalpy": + // No dedicated Enthalpy ComponentOutput; carry it on heatTransferRate. + return "heatTransferRate"; + default: + return "temperature"; + } +} + +/** Probe `measure` kind → list of z-factors it can calibrate (pairing matrix). */ +const FACTOR_COMPATIBILITY: Record = { + SST: ["z_ua"], + SDT: ["z_ua"], + SH: ["z_ua"], + DSH: ["z_ua"], + SC: ["z_ua"], + Capacity: ["z_ua", "z_flow", "z_power"], + P: ["z_dp"], + MassFlow: ["z_flow", "z_etav"], + // Kinds not used for z-factor calibration (raw T, Enthalpy) — no compat. + DGT: [], + T: [], + Enthalpy: [], +}; + /** Convert UI measure value to SI expected by the solver (temps → K). */ function measureSetpointSi(meta: ParamMeta, value: number): number { const unit = (meta.unit ?? "").toLowerCase(); diff --git a/apps/web/src/lib/edgeInsert.ts b/apps/web/src/lib/edgeInsert.ts index ed73cc1..037d596 100644 --- a/apps/web/src/lib/edgeInsert.ts +++ b/apps/web/src/lib/edgeInsert.ts @@ -25,10 +25,26 @@ export const PIPE_TYPES = new Set([ "PipeAir", ]); +/** + * 2-port pass-through components that splice into an existing wire via + * `insertOnEdge`. Includes the Pipe family (above) plus the `Probe` + * measurement tap (calibration redesign — every measurement lives on a Probe + * the user drops onto a line). + */ +export const INSERTABLE_ON_EDGE_TYPES: ReadonlySet = new Set([ + ...PIPE_TYPES, + "Probe", +]); + export function isPipeType(type: string): boolean { return PIPE_TYPES.has(type) || type.includes("Pipe") || type === "AirDuct"; } +/** True for any component the canvas should splice onto a wire on drop. */ +export function isInsertableOnEdge(type: string): boolean { + return isPipeType(type) || type === "Probe"; +} + export function pipeMediaKind(type: string, params?: Record): MediaKind { if (type === "WaterPipe" || type === "PipeWater") return "water"; if (type === "AirDuct" || type === "PipeAir") return "air"; diff --git a/crates/cli/src/batch.rs b/crates/cli/src/batch.rs index 2a201bb..33c54fe 100644 --- a/crates/cli/src/batch.rs +++ b/crates/cli/src/batch.rs @@ -315,7 +315,8 @@ fn process_single_file( initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), } } } @@ -425,7 +426,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -439,7 +441,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, ]; @@ -483,7 +486,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -497,7 +501,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, ]; @@ -529,7 +534,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -543,7 +549,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, ]; @@ -571,7 +578,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -585,7 +593,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 1000, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, ]; @@ -643,7 +652,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 100, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }], }; @@ -667,7 +677,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "ok2.json".to_string(), @@ -681,7 +692,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 60, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "fail.json".to_string(), @@ -695,7 +707,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, ]; diff --git a/crates/cli/src/config.rs b/crates/cli/src/config.rs index b94311b..ed29bd0 100644 --- a/crates/cli/src/config.rs +++ b/crates/cli/src/config.rs @@ -452,6 +452,9 @@ pub struct SolverConfig { /// Solver strategy: "newton" or "picard". #[serde(default = "default_solver_strategy")] pub strategy: String, + /// Dense linear backend: "nalgebra" (default) or "faer" (Story 1.5 opt-in). + #[serde(default = "default_linear_backend")] + pub linear_backend: String, /// Maximum iterations. #[serde(default = "default_max_iterations")] pub max_iterations: usize, @@ -470,6 +473,10 @@ fn default_solver_strategy() -> String { "newton".to_string() } +fn default_linear_backend() -> String { + "nalgebra".to_string() +} + fn default_max_iterations() -> usize { 100 } @@ -482,6 +489,7 @@ impl Default for SolverConfig { fn default() -> Self { Self { strategy: default_solver_strategy(), + linear_backend: default_linear_backend(), max_iterations: default_max_iterations(), tolerance: default_tolerance(), timeout_ms: 0, diff --git a/crates/cli/src/run.rs b/crates/cli/src/run.rs index cb6cddc..ab35884 100644 --- a/crates/cli/src/run.rs +++ b/crates/cli/src/run.rs @@ -372,7 +372,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } @@ -407,7 +408,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } }; @@ -717,7 +719,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } }, @@ -737,7 +740,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } } @@ -782,7 +786,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } }, @@ -823,7 +828,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } }, @@ -867,7 +873,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } } @@ -887,7 +894,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } } @@ -957,7 +965,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } } @@ -983,7 +992,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; let find_component = |name: &str| { config @@ -1042,7 +1052,63 @@ fn execute_simulation( // Register declared control loops (saturated-PI, co-solved). Must happen // BEFORE finalize() so the actuator is wired to the component's CalibIndices. + // + // Routing (calibration redesign): controls whose actuator is a z-factor + // (z_flow/z_flow_eco/z_dp/z_ua/z_power/z_etav) go through PLAIN inverse + // embedding (+1 residual measured−target, +1 unknown z) — not the + // SaturatedController (2+2 with integrator, meant for physical actuators). + let control_error = |msg: String| SimulationResult { + input: input_name.to_string(), + status: SimulationStatus::Error, + convergence: None, + iterations: None, + state: None, + performance: None, + error: Some(msg), + failure_diagnostics: None, + initialization_diagnostics: None, + dof: None, + elapsed_ms, + raw_state_vector: None, + solved_variables: Vec::new(), + }; for control in &config.controls { + let factor = control.actuator.factor.trim(); + // Plain embedding applies to single-point z-factor calibration only. + // A control WITH `objectives` is a supervisory override network + // (selector semantics) and always keeps the SaturatedController path, + // regardless of the actuator factor. + if entropyk_core::normalize_factor_name(factor).is_some() && control.objectives.is_empty() { + match build_plain_z_embedding(control, factor) { + Ok((constraint, bounded_var, actuator_id)) => { + if let Err(e) = system.add_constraint(constraint) { + return control_error(format!( + "Failed to add constraint for control '{}': {:?}", + control.id, e + )); + } + if let Err(e) = system.add_bounded_variable(bounded_var) { + return control_error(format!( + "Failed to add actuator for control '{}': {:?}", + control.id, e + )); + } + if let Err(e) = system.link_constraint_to_control( + &entropyk_solver::inverse::ConstraintId::new(control.id.clone()), + &actuator_id, + ) { + return control_error(format!( + "Failed to link control '{}': {:?}", + control.id, e + )); + } + } + Err(msg) => { + return control_error(format!("Invalid control '{}': {}", control.id, msg)); + } + } + continue; + } match build_saturated_control(control) { Ok((bounded_var, controller)) => { if let Err(e) = system.add_bounded_variable(bounded_var) { @@ -1061,7 +1127,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } system.add_saturated_controller(controller); @@ -1079,7 +1146,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } } @@ -1253,7 +1321,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } system.add_free_actuator(var_id); @@ -1275,7 +1344,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } } @@ -1296,7 +1366,8 @@ fn execute_simulation( initialization_diagnostics: None, dof: None, elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } @@ -1524,10 +1595,35 @@ fn execute_simulation( initialization_diagnostics: None, dof: Some(dof_summary), elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; } + // Dense linear backend dispatch (Story 1.5): "nalgebra" (default) or + // "faer"; env ENTROPYK_LINEAR_BACKEND overrides for the test harness. + if !matches!(config.solver.linear_backend.as_str(), "nalgebra" | "faer") { + return SimulationResult { + input: input_name.to_string(), + status: SimulationStatus::Error, + convergence: None, + iterations: None, + state: None, + performance: None, + error: Some(format!( + "Unsupported linear backend '{}'. Use 'nalgebra' or 'faer'.", + config.solver.linear_backend + )), + failure_diagnostics: None, + initialization_diagnostics: None, + dof: Some(dof_summary), + elapsed_ms, + raw_state_vector: None, + solved_variables: Vec::new(), + }; + } + entropyk_solver::linear::set_linear_backend_override(Some(&config.solver.linear_backend)); + // Build + run the configured strategy for a single solve. CLI // `solver.max_iterations` / `tolerance` are now honoured by every stage // (previously the defaults capped Newton at 100 iterations). Reused by the @@ -1799,7 +1895,8 @@ fn execute_simulation( initialization_diagnostics: Some(initialization_diagnostics), dof: Some(dof_summary), elapsed_ms, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), } } } @@ -2968,9 +3065,8 @@ fn bphx_dp_correlation_from_params( "{component_name}: dp_correlation must be a string (got {value})" ))); }; - BphxDpCorrelation::parse(raw).map_err(|msg| { - CliError::Config(format!("{component_name}: {msg}")) - }) + BphxDpCorrelation::parse(raw) + .map_err(|msg| CliError::Config(format!("{component_name}: {msg}"))) } /// Extract calibration Z-factors for BphxEvaporator/BphxCondenser from JSON params. @@ -3134,6 +3230,43 @@ fn parse_component_output( } /// Builds a bounded actuator + saturated controller from a control config. +/// Plain inverse embedding for z-factor calibration controls (calibration +/// redesign, WS-2): one `Constraint` (measured − target) + one `BoundedVariable` +/// (the z-factor), linked 1:1. The bounded var id reuses the +/// `{component}__{z_factor}` convention so `finalize()` wires it to the +/// component's matching `CalibIndices` slot (system.rs). +fn build_plain_z_embedding( + control: &crate::config::ControlConfig, + factor: &str, +) -> Result< + ( + entropyk_solver::inverse::Constraint, + entropyk_solver::inverse::BoundedVariable, + entropyk_solver::inverse::BoundedVariableId, + ), + String, +> { + use entropyk_solver::inverse::{BoundedVariable, BoundedVariableId, Constraint, ConstraintId}; + + let output = parse_component_output(&control.measure.component, &control.measure.output)?; + let actuator_id = + BoundedVariableId::new(saturated_actuator_id(&control.actuator.component, factor)); + let bounded_var = BoundedVariable::with_component( + actuator_id.clone(), + &control.actuator.component, + control.actuator.initial, + control.actuator.min, + control.actuator.max, + ) + .map_err(|e| format!("invalid actuator bounds: {e:?}"))?; + let constraint = Constraint::new( + ConstraintId::new(control.id.clone()), + output, + control.target, + ); + Ok((constraint, bounded_var, actuator_id)) +} + fn build_saturated_control( control: &crate::config::ControlConfig, ) -> Result< @@ -5249,6 +5382,55 @@ fn create_component( Ok(Box::new(load)) } + "Probe" => { + use entropyk::{ComponentFluidId, Enthalpy, Port, Probe, ProbeMeasure, Pressure}; + // Calibration-redesign measurement tap (zero-residual, 2-port splice). + // The user drops it onto a wire (UI `insertOnEdge`) and picks a + // measure kind; the freed z-factor on a calibrated component + // elsewhere is paired with this Probe via a calibration control + // whose `measure.component` names the Probe. + let measure_str = params + .get("measure") + .and_then(|v| v.as_str()) + .unwrap_or("SH"); + let measure = ProbeMeasure::parse(measure_str).map_err(CliError::Config)?; + + let fluid = params + .get("fluid") + .and_then(|v| v.as_str()) + .filter(|s| !s.is_empty()) + .unwrap_or_else(|| _primary_fluid.as_str()) + .to_string(); + + // Nominal edge state for the port pair (mirrors Pipe's idiom — + // both pairs carry identical nominal values, the live solver state + // overrides them via `set_system_context`). + let p_bar = params.get("p_bar").and_then(|v| v.as_f64()).unwrap_or(5.0); + let h_j = params + .get("h_j_kg") + .and_then(|v| v.as_f64()) + .or_else(|| { + params + .get("h_kj_kg") + .and_then(|v| v.as_f64()) + .map(|k| k * 1000.0) + }) + .unwrap_or(250_000.0); + + let mk = || { + Port::new( + ComponentFluidId::new(&fluid), + Pressure::from_bar(p_bar), + Enthalpy::from_joules_per_kg(h_j), + ) + }; + let probe = Probe::new(measure, &fluid, mk(), mk()) + .with_fluid_backend(backend) + .connect(mk(), mk()) + .map_err(CliError::Component)?; + Ok(Box::new(probe)) + } + "Anchor" | "RefrigerantNode" => { use entropyk::{Anchor, AnchorConstraint}; // BOLT `BoundaryNode.Refrigerant.Node` equivalent: inline spec node @@ -5362,7 +5544,7 @@ fn create_component( } _ => Err(CliError::Config(format!( - "Unknown component type: '{}'. Supported: IsentropicCompressor, IsenthalpicExpansionValve, ScrewEconomizerCompressor, CentrifugalCompressor, CapillaryTube, MchxCondenserCoil, FloodedEvaporator, BphxEvaporator, BphxCondenser, FreeCoolingExchanger, Condenser, CondenserCoil, Evaporator, EvaporatorCoil, HeatExchanger, Compressor, ExpansionValve, ReversingValve, Pump, Fan, Pipe, RefrigerantPipe, WaterPipe, AirDuct, Placeholder", + "Unknown component type: '{}'. Supported: IsentropicCompressor, IsenthalpicExpansionValve, ScrewEconomizerCompressor, CentrifugalCompressor, CapillaryTube, MchxCondenserCoil, FloodedEvaporator, BphxEvaporator, BphxCondenser, FreeCoolingExchanger, Condenser, CondenserCoil, Evaporator, EvaporatorCoil, HeatExchanger, Compressor, ExpansionValve, ReversingValve, Pump, Fan, Pipe, RefrigerantPipe, WaterPipe, AirDuct, Probe, Anchor, Placeholder", component_type ))), } @@ -5672,7 +5854,8 @@ mod tests { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; let json = serde_json::to_string_pretty(&result).unwrap(); @@ -5839,10 +6022,8 @@ mod tests { .unwrap(), BphxDpCorrelation::SimplifiedChannel ); - let params = std::collections::HashMap::from([( - "DpCorrelation".to_string(), - json!("Martin1996"), - )]); + let params = + std::collections::HashMap::from([("DpCorrelation".to_string(), json!("Martin1996"))]); assert_eq!( bphx_dp_correlation_from_params(¶ms, "BphxCondenser").unwrap(), BphxDpCorrelation::Martin1996 @@ -5853,10 +6034,8 @@ mod tests { fn test_bphx_dp_correlation_rejects_unknown() { use serde_json::json; - let params = std::collections::HashMap::from([( - "dp_correlation".to_string(), - json!("Amalfi2016"), - )]); + let params = + std::collections::HashMap::from([("dp_correlation".to_string(), json!("Amalfi2016"))]); let error = bphx_dp_correlation_from_params(¶ms, "BphxEvaporator").unwrap_err(); assert!( error.to_string().contains("unsupported dp_correlation"), @@ -5868,15 +6047,9 @@ mod tests { fn test_bphx_dp_correlation_rejects_non_string() { use serde_json::json; - let params = std::collections::HashMap::from([( - "dp_correlation".to_string(), - json!(42), - )]); + let params = std::collections::HashMap::from([("dp_correlation".to_string(), json!(42))]); let error = bphx_dp_correlation_from_params(¶ms, "BphxEvaporator").unwrap_err(); - assert!( - error.to_string().contains("must be a string"), - "{error}" - ); + assert!(error.to_string().contains("must be a string"), "{error}"); } fn comp_from_json(v: serde_json::Value) -> crate::config::ComponentConfig { diff --git a/crates/cli/tests/batch_execution.rs b/crates/cli/tests/batch_execution.rs index d03fd5c..0175ae0 100644 --- a/crates/cli/tests/batch_execution.rs +++ b/crates/cli/tests/batch_execution.rs @@ -92,7 +92,8 @@ fn test_simulation_result_statuses() { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "fail.json".to_string(), @@ -106,7 +107,8 @@ fn test_simulation_result_statuses() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "timeout.json".to_string(), @@ -120,7 +122,8 @@ fn test_simulation_result_statuses() { initialization_diagnostics: None, dof: None, elapsed_ms: 1000, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, ]; @@ -163,7 +166,8 @@ fn test_batch_aggregator_csv_output() { initialization_diagnostics: None, dof: None, elapsed_ms: 150, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "scenario2.json".to_string(), @@ -183,7 +187,8 @@ fn test_batch_aggregator_csv_output() { initialization_diagnostics: None, dof: None, elapsed_ms: 200, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "scenario3.json".to_string(), @@ -197,7 +202,8 @@ fn test_batch_aggregator_csv_output() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, ]; @@ -229,7 +235,8 @@ fn test_batch_aggregator_json_summary() { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "test2.json".to_string(), @@ -243,7 +250,8 @@ fn test_batch_aggregator_json_summary() { initialization_diagnostics: None, dof: None, elapsed_ms: 75, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, SimulationResult { input: "test3.json".to_string(), @@ -257,7 +265,8 @@ fn test_batch_aggregator_json_summary() { initialization_diagnostics: None, dof: None, elapsed_ms: 5000, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }, ]; @@ -322,7 +331,8 @@ fn test_batch_summary_csv_with_convergence() { initialization_diagnostics: None, dof: None, elapsed_ms: 300, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }]; let summary = BatchSummary { diff --git a/crates/cli/tests/calibration_sdt.rs b/crates/cli/tests/calibration_sdt.rs new file mode 100644 index 0000000..92bead0 --- /dev/null +++ b/crates/cli/tests/calibration_sdt.rs @@ -0,0 +1,190 @@ +//! Calibration redesign (WS-2/WS-3/WS-4) — AC#1: BPHX condenser SDT +//! calibration via plain inverse embedding converges and back-solves z_ua +//! (the case that diverged with a singular Jacobian before the redesign). + +use entropyk_cli::run::{run_simulation, SimulationResult, SimulationStatus}; +use tempfile::tempdir; + +fn run_config(json: &str) -> SimulationResult { + let dir = tempdir().unwrap(); + let path = dir.path().join("bphx_sdt_calib.json"); + std::fs::write(&path, json).unwrap(); + run_simulation(&path, None, false).unwrap() +} + +#[test] +fn bphx_condenser_sdt_calibration_converges_and_solves_z_ua() { + let json = r#" +{ + "name": "BPHX Water-Cooled Chiller R134a", + "description": "Full vapor-compression cycle with brazed-plate evaporator and condenser (4-port). Both Bphx units use emergent_pressure so outlet closures (superheat / subcooling) square the DoF like Condenser/Evaporator.", + "fluid": "R134a", + "fluid_backend": "CoolProp", + "circuits": [ + { + "id": 0, + "name": "Refrigerant + water loops", + "components": [ + { + "type": "IsentropicCompressor", + "name": "comp", + "isentropic_efficiency": 0.7, + "t_cond_k": 318.15, + "t_evap_k": 278.15, + "superheat_k": 5.0, + "emergent_pressure": true, + "displacement_m3": 6.5e-05, + "speed_hz": 50.0, + "volumetric_efficiency": 0.92 + }, + { + "type": "BphxCondenser", + "name": "cond", + "refrigerant": "R134a", + "secondary_fluid": "Water", + "n_plates": 40, + "plate_length_m": 0.4, + "plate_width_m": 0.12, + "target_subcooling_k": 5.0, + "emergent_pressure": true, + "correlation": "Longo2004", + "dp_correlation": "SimplifiedChannel", + "ua": 2500.0 + }, + { + "type": "IsenthalpicExpansionValve", + "name": "exv", + "t_evap_k": 278.15, + "emergent_pressure": true + }, + { + "type": "BphxEvaporator", + "name": "evap", + "refrigerant": "R134a", + "secondary_fluid": "Water", + "n_plates": 40, + "plate_length_m": 0.4, + "plate_width_m": 0.12, + "target_superheat_k": 5.0, + "emergent_pressure": true, + "correlation": "Longo2004", + "dp_correlation": "SimplifiedChannel", + "ua": 2000.0 + }, + { + "type": "BrineSource", + "name": "cond_water_in", + "fluid": "Water", + "p_set_bar": 2.0, + "t_set_c": 30.0, + "m_flow_kg_s": 0.4, + "fix_pressure": false, + "fix_temperature": true, + "fix_mass_flow": true + }, + { + "type": "BrineSink", + "name": "cond_water_out", + "fluid": "Water", + "p_back_bar": 2.0, + "fix_pressure": true + }, + { + "type": "BrineSource", + "name": "evap_water_in", + "fluid": "Water", + "p_set_bar": 3.0, + "t_set_c": 12.0, + "m_flow_kg_s": 0.5, + "fix_pressure": false, + "fix_temperature": true, + "fix_mass_flow": true + }, + { + "type": "BrineSink", + "name": "evap_water_out", + "fluid": "Water", + "p_back_bar": 3.0, + "fix_pressure": true + } + ], + "edges": [ + { + "from": "comp:outlet", + "to": "cond:inlet" + }, + { + "from": "cond:outlet", + "to": "exv:inlet" + }, + { + "from": "exv:outlet", + "to": "evap:inlet" + }, + { + "from": "evap:outlet", + "to": "comp:inlet" + }, + { + "from": "cond_water_in:outlet", + "to": "cond:secondary_inlet" + }, + { + "from": "cond:secondary_outlet", + "to": "cond_water_out:inlet" + }, + { + "from": "evap_water_in:outlet", + "to": "evap:secondary_inlet" + }, + { + "from": "evap:secondary_outlet", + "to": "evap_water_out:inlet" + } + ] + } + ], + "solver": { + "strategy": "fallback", + "max_iterations": 300, + "tolerance": 1e-06, + "timeout_ms": 60000 + }, + "controls": [ + { + "type": "SaturatedController", + "id": "sdt_calib", + "measure": { + "component": "cond", + "output": "saturationTemperature" + }, + "actuator": { + "component": "cond", + "factor": "z_ua", + "initial": 0.3, + "min": 0.05, + "max": 2.0 + }, + "target": 315.0 + } + ] +} +"#; + let result = run_config(json); + assert!( + matches!(result.status, SimulationStatus::Converged), + "BPHX SDT calibration must converge (was the singular-Jacobian case): {:?} ({:?})", + result.status, + result.error + ); + let solved = result + .solved_variables + .iter() + .find(|v| v.variable == "z_ua" && v.component.as_deref() == Some("cond")) + .expect("solved_variables must contain cond/z_ua"); + assert!( + solved.value > 0.05 && solved.value < 2.0, + "z_ua must solve within bounds, got {}", + solved.value + ); +} diff --git a/crates/cli/tests/probe_calibration.rs b/crates/cli/tests/probe_calibration.rs new file mode 100644 index 0000000..0d3ef56 --- /dev/null +++ b/crates/cli/tests/probe_calibration.rs @@ -0,0 +1,193 @@ +//! Calibration redesign (HARD RULE — Probe for ALL measurements) — integration +//! test: a `Probe` node measuring SDT on the condenser refrigerant inlet edge +//! drives the plain inverse embedding for `z_ua`. The Probe is the measurement +//! source (`control.measure.component` names the Probe); the freed z-factor +//! lives on the BPHX condenser (`control.actuator.component`). The two are +//! linked 1:1 (+1 residual on Probe SDT target, +1 unknown z_ua). +//! +//! This is the Probe-based variant of `calibration_sdt.rs`. It proves the +//! end-to-end path the UI emits after the calibration redesign. + +use entropyk_cli::run::{run_simulation, SimulationResult, SimulationStatus}; +use tempfile::tempdir; + +fn run_config(json: &str) -> SimulationResult { + let dir = tempdir().unwrap(); + let path = dir.path().join("probe_sdt_calib.json"); + std::fs::write(&path, json).unwrap(); + run_simulation(&path, None, false).unwrap() +} + +#[test] +fn probe_based_sdt_calibration_converges_and_solves_z_ua() { + let json = r#" +{ + "name": "Probe-based SDT calibration (R134a BPHX chiller)", + "description": "Same vapor-compression cycle as calibration_sdt.rs, but the SDT measurement lives on a Probe node spliced into the condenser refrigerant inlet edge. control.measure.component = the Probe name.", + "fluid": "R134a", + "fluid_backend": "CoolProp", + "circuits": [ + { + "id": 0, + "name": "Refrigerant + water loops", + "components": [ + { + "type": "IsentropicCompressor", + "name": "comp", + "isentropic_efficiency": 0.7, + "t_cond_k": 318.15, + "t_evap_k": 278.15, + "superheat_k": 5.0, + "emergent_pressure": true, + "displacement_m3": 6.5e-05, + "speed_hz": 50.0, + "volumetric_efficiency": 0.92 + }, + { + "type": "BphxCondenser", + "name": "cond", + "refrigerant": "R134a", + "secondary_fluid": "Water", + "n_plates": 40, + "plate_length_m": 0.4, + "plate_width_m": 0.12, + "target_subcooling_k": 5.0, + "emergent_pressure": true, + "correlation": "Longo2004", + "dp_correlation": "SimplifiedChannel", + "ua": 2500.0 + }, + { + "type": "Probe", + "name": "cond_sdt_probe", + "measure": "SDT", + "fluid": "R134a" + }, + { + "type": "IsenthalpicExpansionValve", + "name": "exv", + "t_evap_k": 278.15, + "emergent_pressure": true + }, + { + "type": "BphxEvaporator", + "name": "evap", + "refrigerant": "R134a", + "secondary_fluid": "Water", + "n_plates": 40, + "plate_length_m": 0.4, + "plate_width_m": 0.12, + "target_superheat_k": 5.0, + "emergent_pressure": true, + "correlation": "Longo2004", + "dp_correlation": "SimplifiedChannel", + "ua": 2000.0 + }, + { + "type": "BrineSource", + "name": "cond_water_in", + "fluid": "Water", + "p_set_bar": 2.0, + "t_set_c": 30.0, + "m_flow_kg_s": 0.4, + "fix_pressure": false, + "fix_temperature": true, + "fix_mass_flow": true + }, + { + "type": "BrineSink", + "name": "cond_water_out", + "fluid": "Water", + "p_back_bar": 2.0, + "fix_pressure": true + }, + { + "type": "BrineSource", + "name": "evap_water_in", + "fluid": "Water", + "p_set_bar": 3.0, + "t_set_c": 12.0, + "m_flow_kg_s": 0.5, + "fix_pressure": false, + "fix_temperature": true, + "fix_mass_flow": true + }, + { + "type": "BrineSink", + "name": "evap_water_out", + "fluid": "Water", + "p_back_bar": 3.0, + "fix_pressure": true + } + ], + "edges": [ + { "from": "comp:outlet", "to": "cond_sdt_probe:inlet" }, + { "from": "cond_sdt_probe:outlet","to": "cond:inlet" }, + { "from": "cond:outlet", "to": "exv:inlet" }, + { "from": "exv:outlet", "to": "evap:inlet" }, + { "from": "evap:outlet", "to": "comp:inlet" }, + { "from": "cond_water_in:outlet", "to": "cond:secondary_inlet" }, + { "from": "cond:secondary_outlet","to": "cond_water_out:inlet" }, + { "from": "evap_water_in:outlet", "to": "evap:secondary_inlet" }, + { "from": "evap:secondary_outlet","to": "evap_water_out:inlet" } + ] + } + ], + "solver": { + "strategy": "fallback", + "max_iterations": 300, + "tolerance": 1e-06, + "timeout_ms": 60000 + }, + "controls": [ + { + "type": "SaturatedController", + "id": "probe_sdt_calib", + "measure": { + "component": "cond_sdt_probe", + "output": "saturationTemperature" + }, + "actuator": { + "component": "cond", + "factor": "z_ua", + "initial": 0.3, + "min": 0.05, + "max": 2.0 + }, + "target": 315.0 + } + ] +} +"#; + let result = run_config(json); + assert!( + matches!(result.status, SimulationStatus::Converged), + "Probe-based SDT calibration must converge: {:?} ({:?})", + result.status, + result.error + ); + let solved = result + .solved_variables + .iter() + .find(|v| v.variable == "z_ua" && v.component.as_deref() == Some("cond")) + .expect("solved_variables must contain cond/z_ua"); + eprintln!("DIAG z_ua résolu = {}", solved.value); + eprintln!("DIAG cible SDT = 315.0 K (41.85°C)"); + if let Some(state) = result.state.as_ref() { + for e in state.iter() { + if e.target.as_deref() == Some("cond") || e.source.as_deref() == Some("cond") { + eprintln!("DIAG edge {}→{} P={}bar T_sat={}°C T={}°C", + e.source.as_deref().unwrap_or("?"), + e.target.as_deref().unwrap_or("?"), + e.pressure_bar, + e.saturation_temperature_c.unwrap_or(f64::NAN), + e.temperature_c.unwrap_or(f64::NAN)); + } + } + } + assert!( + solved.value > 0.05 && solved.value < 2.0, + "z_ua must solve within bounds, got {}", + solved.value + ); +} diff --git a/crates/cli/tests/probe_sst_evap.rs b/crates/cli/tests/probe_sst_evap.rs new file mode 100644 index 0000000..ce37ce1 --- /dev/null +++ b/crates/cli/tests/probe_sst_evap.rs @@ -0,0 +1,194 @@ +//! Calibration redesign (HARD RULE — Probe for ALL measurements) — integration +//! test: a `Probe` node measuring SDT on the condenser refrigerant inlet edge +//! drives the plain inverse embedding for `z_ua`. The Probe is the measurement +//! source (`control.measure.component` names the Probe); the freed z-factor +//! lives on the BPHX condenser (`control.actuator.component`). The two are +//! linked 1:1 (+1 residual on Probe SDT target, +1 unknown z_ua). +//! +//! This is the Probe-based variant of `calibration_sdt.rs`. It proves the +//! end-to-end path the UI emits after the calibration redesign. + +use entropyk_cli::run::{run_simulation, SimulationResult, SimulationStatus}; +use tempfile::tempdir; + +fn run_config(json: &str) -> SimulationResult { + let dir = tempdir().unwrap(); + let path = dir.path().join("probe_sdt_calib.json"); + std::fs::write(&path, json).unwrap(); + run_simulation(&path, None, false).unwrap() +} + +#[test] +fn probe_based_sst_evap_calibration() { + let json = r#" +{ + "name": "Probe-based SDT calibration (R134a BPHX chiller)", + "description": "Same vapor-compression cycle as calibration_sdt.rs, but the SDT measurement lives on a Probe node spliced into the condenser refrigerant inlet edge. control.measure.component = the Probe name.", + "fluid": "R134a", + "fluid_backend": "CoolProp", + "circuits": [ + { + "id": 0, + "name": "Refrigerant + water loops", + "components": [ + { + "type": "IsentropicCompressor", + "name": "comp", + "isentropic_efficiency": 0.7, + "t_cond_k": 318.15, + "t_evap_k": 278.15, + "superheat_k": 5.0, + "emergent_pressure": true, + "displacement_m3": 6.5e-05, + "speed_hz": 50.0, + "volumetric_efficiency": 0.92 + }, + { + "type": "BphxCondenser", + "name": "cond", + "refrigerant": "R134a", + "secondary_fluid": "Water", + "n_plates": 40, + "plate_length_m": 0.4, + "plate_width_m": 0.12, + "target_subcooling_k": 5.0, + "emergent_pressure": true, + "correlation": "Longo2004", + "dp_correlation": "SimplifiedChannel", + "ua": 2500.0 + }, + { + "type": "Probe", + "name": "evap_sst_probe", + "measure": "SDT", + "fluid": "R134a" + }, + { + "type": "IsenthalpicExpansionValve", + "name": "exv", + "t_evap_k": 278.15, + "emergent_pressure": true + }, + { + "type": "BphxEvaporator", + "name": "evap", + "refrigerant": "R134a", + "secondary_fluid": "Water", + "n_plates": 40, + "plate_length_m": 0.4, + "plate_width_m": 0.12, + "target_superheat_k": 5.0, + "emergent_pressure": true, + "correlation": "Longo2004", + "dp_correlation": "SimplifiedChannel", + "ua": 2000.0 + }, + { + "type": "BrineSource", + "name": "cond_water_in", + "fluid": "Water", + "p_set_bar": 2.0, + "t_set_c": 30.0, + "m_flow_kg_s": 0.4, + "fix_pressure": false, + "fix_temperature": true, + "fix_mass_flow": true + }, + { + "type": "BrineSink", + "name": "cond_water_out", + "fluid": "Water", + "p_back_bar": 2.0, + "fix_pressure": true + }, + { + "type": "BrineSource", + "name": "evap_water_in", + "fluid": "Water", + "p_set_bar": 3.0, + "t_set_c": 12.0, + "m_flow_kg_s": 0.5, + "fix_pressure": false, + "fix_temperature": true, + "fix_mass_flow": true + }, + { + "type": "BrineSink", + "name": "evap_water_out", + "fluid": "Water", + "p_back_bar": 3.0, + "fix_pressure": true + } + ], + "edges": [ + { "from": "exv:outlet", "to": "evap_sst_probe:inlet" }, + { "from": "evap_sst_probe:outlet","to": "cond:inlet" }, + { "from": "cond:outlet", "to": "exv:inlet" }, + { "from": "exv:outlet", "to": "evap:inlet" }, + { "from": "evap:outlet", "to": "evap_sst_probe:inlet" }, + { "from": "evap_sst_probe:outlet", "to": "comp:inlet" }, + { "from": "cond_water_in:outlet", "to": "cond:secondary_inlet" }, + { "from": "cond:secondary_outlet","to": "cond_water_out:inlet" }, + { "from": "evap_water_in:outlet", "to": "evap:secondary_inlet" }, + { "from": "evap:secondary_outlet","to": "evap_water_out:inlet" } + ] + } + ], + "solver": { + "strategy": "fallback", + "max_iterations": 300, + "tolerance": 1e-06, + "timeout_ms": 60000 + }, + "controls": [ + { + "type": "SaturatedController", + "id": "probe_sdt_calib", + "measure": { + "component": "evap_sst_probe", + "output": "saturationTemperature" + }, + "actuator": { + "component": "evap", + "factor": "z_ua", + "initial": 0.3, + "min": 0.05, + "max": 2.0 + }, + "target": 277.55 + } + ] +} +"#; + let result = run_config(json); + assert!( + matches!(result.status, SimulationStatus::Converged), + "Probe-based SDT calibration must converge: {:?} ({:?})", + result.status, + result.error + ); + let solved = result + .solved_variables + .iter() + .find(|v| v.variable == "z_ua" && v.component.as_deref() == Some("evap")) + .expect("solved_variables must contain cond/z_ua"); + eprintln!("DIAG z_ua résolu = {}", solved.value); + eprintln!("DIAG cible SST = 277.55 K (41.85°C)"); + if let Some(state) = result.state.as_ref() { + for e in state.iter() { + if e.target.as_deref() == Some("evap") || e.source.as_deref() == Some("evap") { + eprintln!("DIAG edge {}→{} P={}bar T_sat={}°C T={}°C", + e.source.as_deref().unwrap_or("?"), + e.target.as_deref().unwrap_or("?"), + e.pressure_bar, + e.saturation_temperature_c.unwrap_or(f64::NAN), + e.temperature_c.unwrap_or(f64::NAN)); + } + } + } + assert!( + solved.value > 0.05 && solved.value < 2.0, + "z_ua must solve within bounds, got {}", + solved.value + ); +} diff --git a/crates/cli/tests/single_run.rs b/crates/cli/tests/single_run.rs index 4dce070..952af8b 100644 --- a/crates/cli/tests/single_run.rs +++ b/crates/cli/tests/single_run.rs @@ -29,7 +29,8 @@ fn test_simulation_result_serialization() { initialization_diagnostics: None, dof: None, elapsed_ms: 50, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; let json = serde_json::to_string_pretty(&result).unwrap(); @@ -70,7 +71,8 @@ fn test_error_result_serialization() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; let json = serde_json::to_string(&result).unwrap(); @@ -96,7 +98,8 @@ fn test_error_result_serializes_failure_diagnostics() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; let json = serde_json::to_string(&result).unwrap(); @@ -2303,7 +2306,8 @@ fn test_structural_failure_serializes_without_diagnostics() { initialization_diagnostics: None, dof: None, elapsed_ms: 0, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; let json = serde_json::to_string(&result).unwrap(); @@ -2348,7 +2352,8 @@ fn test_success_result_does_not_include_failure_diagnostics() { initialization_diagnostics: None, dof: None, elapsed_ms: 120, - raw_state_vector: None, solved_variables: Vec::new(), + raw_state_vector: None, + solved_variables: Vec::new(), }; let json = serde_json::to_string(&result).unwrap(); diff --git a/crates/components/src/heat_exchanger/bphx_condenser.rs b/crates/components/src/heat_exchanger/bphx_condenser.rs index 03b6cc7..866d1d8 100644 --- a/crates/components/src/heat_exchanger/bphx_condenser.rs +++ b/crates/components/src/heat_exchanger/bphx_condenser.rs @@ -555,6 +555,10 @@ impl Component for BphxCondenser { self.inner.energy_transfers(state) } + fn measure_output(&self, kind: crate::MeasuredOutput, state: &StateSlice) -> Option { + self.inner.measure_output(kind, state) + } + fn set_fluid_backend_from_builder( &mut self, backend: std::sync::Arc, diff --git a/crates/components/src/heat_exchanger/bphx_dp.rs b/crates/components/src/heat_exchanger/bphx_dp.rs index 93ccf6f..c85c842 100644 --- a/crates/components/src/heat_exchanger/bphx_dp.rs +++ b/crates/components/src/heat_exchanger/bphx_dp.rs @@ -36,7 +36,10 @@ impl BphxDpCorrelation { /// Parse CLI/UI identifiers. Rejects unknown strings (no silent fallback). pub fn parse(name: &str) -> Result { - let key = name.trim().to_ascii_lowercase().replace(['-', '_', ' '], ""); + let key = name + .trim() + .to_ascii_lowercase() + .replace(['-', '_', ' '], ""); match key.as_str() { "simplifiedchannel" | "simplified" | "channel" | "default" => { Ok(Self::SimplifiedChannel) @@ -104,13 +107,7 @@ fn simplified_channel_fanning(re: f64) -> (f64, f64) { return (f_turb, df_turb); } - let f = entropyk_core::smoothing::cubic_blend( - f_lam, - f_turb, - re, - RE_LAMINAR, - RE_TURBULENT, - ); + let f = entropyk_core::smoothing::cubic_blend(f_lam, f_turb, re, RE_LAMINAR, RE_TURBULENT); // d/dRe of blend(f_lam(Re), f_turb(Re), Re): product rule on cubic Hermite. let df_blend_dx = entropyk_core::smoothing::cubic_blend_derivative( f_lam, @@ -178,12 +175,10 @@ fn martin1996_fanning(re: f64, chevron_angle_deg: f64) -> (f64, f64) { let t = ((re_safe - RE_LO) / (RE_HI - RE_LO)).clamp(0.0, 1.0); let w_b = t * t * (3.0 - 2.0 * t); let w_a = 1.0 - w_b; - let df0_blend = entropyk_core::smoothing::cubic_blend_derivative( - f0_l, f0_t, re_safe, RE_LO, RE_HI, - ); - let df1_blend = entropyk_core::smoothing::cubic_blend_derivative( - f1_l, f1_t, re_safe, RE_LO, RE_HI, - ); + let df0_blend = + entropyk_core::smoothing::cubic_blend_derivative(f0_l, f0_t, re_safe, RE_LO, RE_HI); + let df1_blend = + entropyk_core::smoothing::cubic_blend_derivative(f1_l, f1_t, re_safe, RE_LO, RE_HI); let df0 = df0_blend + w_a * df0_l + w_b * df0_t; let df1 = df1_blend + w_a * df1_l + w_b * df1_t; (f0, df0, f1, df1) @@ -324,14 +319,9 @@ mod tests { 1.0, ) .unwrap(); - let martin = evaluate_channel_pressure_drop( - BphxDpCorrelation::Martin1996, - &geo(), - g, - rho, - 1.0, - ) - .unwrap(); + let martin = + evaluate_channel_pressure_drop(BphxDpCorrelation::Martin1996, &geo(), g, rho, 1.0) + .unwrap(); assert!( (simp.delta_p_pa - martin.delta_p_pa).abs() > 1.0, "Martin and Simplified should differ: {} vs {}", @@ -381,14 +371,9 @@ mod tests { fn martin_rejects_invalid_chevron() { let mut g = geo(); g.chevron_angle = f64::NAN; - let err = evaluate_channel_pressure_drop( - BphxDpCorrelation::Martin1996, - &g, - 30.0, - 1000.0, - 1.0, - ) - .unwrap_err(); + let err = + evaluate_channel_pressure_drop(BphxDpCorrelation::Martin1996, &g, 30.0, 1000.0, 1.0) + .unwrap_err(); assert!(format!("{err}").contains("chevron")); } @@ -403,10 +388,8 @@ mod tests { ] { let eval = evaluate_channel_pressure_drop(corr, &geo, g0, rho, 1.0).unwrap(); let h = 1e-4; - let plus = - evaluate_channel_pressure_drop(corr, &geo, g0 + h, rho, 1.0).unwrap(); - let minus = - evaluate_channel_pressure_drop(corr, &geo, g0 - h, rho, 1.0).unwrap(); + let plus = evaluate_channel_pressure_drop(corr, &geo, g0 + h, rho, 1.0).unwrap(); + let minus = evaluate_channel_pressure_drop(corr, &geo, g0 - h, rho, 1.0).unwrap(); let fd = (plus.delta_p_pa - minus.delta_p_pa) / (2.0 * h); let rel = (eval.d_delta_p_d_g - fd).abs() / fd.abs().max(1.0); assert!( diff --git a/crates/components/src/heat_exchanger/bphx_evaporator.rs b/crates/components/src/heat_exchanger/bphx_evaporator.rs index dd3de2e..6b09c5e 100644 --- a/crates/components/src/heat_exchanger/bphx_evaporator.rs +++ b/crates/components/src/heat_exchanger/bphx_evaporator.rs @@ -517,6 +517,10 @@ impl Component for BphxEvaporator { self.inner.energy_transfers(state) } + fn measure_output(&self, kind: crate::MeasuredOutput, state: &StateSlice) -> Option { + self.inner.measure_output(kind, state) + } + fn set_fluid_backend_from_builder( &mut self, backend: std::sync::Arc, diff --git a/crates/components/src/heat_exchanger/bphx_exchanger.rs b/crates/components/src/heat_exchanger/bphx_exchanger.rs index e67bbd9..3cf67c0 100644 --- a/crates/components/src/heat_exchanger/bphx_exchanger.rs +++ b/crates/components/src/heat_exchanger/bphx_exchanger.rs @@ -40,11 +40,14 @@ use super::correlation_registry::{ }; use super::eps_ntu::{EpsNtuModel, ExchangerType}; use super::exchanger::{HeatExchanger, HxSideConditions}; +use super::phase_change_entu::{condenser_duty, evaporator_duty}; +use super::sat_domain; use crate::state_machine::{CircuitId, OperationalState, StateManageable}; use crate::{ Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, StateSlice, }; -use entropyk_core::{Calib, Enthalpy, MassFlow, Power}; +use entropyk_core::{Calib, Enthalpy, MassFlow, Power, Pressure}; +use entropyk_fluids::{FluidId, FluidState, Property, Quality}; use std::cell::{Cell, RefCell}; use std::sync::Arc; @@ -498,8 +501,7 @@ impl BphxExchanger { let dp_hot = self.side_channel_dp(state[m_h], rho_hot, z_dp)?; let dp_cold = self.side_channel_dp(state[m_c], rho_cold, z_dp)?; - let a_flow = - self.geometry.channel_flow_area() * self.geometry.n_channels_per_side() as f64; + let a_flow = self.geometry.channel_flow_area() * self.geometry.n_channels_per_side() as f64; if a_flow <= 1e-30 { return Err(ComponentError::InvalidState( "BPHX channel flow area too small for pressure-drop Jacobian".into(), @@ -528,6 +530,193 @@ impl BphxExchanger { let ua = h * self.geometry.area * self.calib().z_ua; self.inner.set_ua_scale(ua / self.inner.ua_nominal()); } + + /// Condenser / evaporator plate HX use Shah \(C^*\to 0\) duty instead of + /// two-stream sensible ε-NTU (see `phase_change_entu`). + fn uses_phase_change_duty(&self) -> bool { + matches!( + self.geometry.exchanger_type, + BphxType::Condenser | BphxType::Evaporator + ) && self.fluid_backend.is_some() + } + + /// \(T_{\mathrm{sat}}(P)\) [K] for the refrigerant fluid id, domain-clamped. + fn tsat_k(&self, p_pa: f64, refrigerant_id: &str) -> Result { + let backend = self.fluid_backend.as_ref().ok_or_else(|| { + ComponentError::CalculationFailed( + "BphxExchanger: FluidBackend required for Tsat".into(), + ) + })?; + let p_pa = + sat_domain::clamp_to_saturation_domain(backend, refrigerant_id, p_pa).unwrap_or(p_pa); + backend + .property( + FluidId::new(refrigerant_id), + Property::Temperature, + FluidState::from_px(Pressure::from_pascals(p_pa), Quality::new(0.5)), + ) + .map_err(ComponentError::from_fluid_error) + } + + /// Secondary inlet \((T [K], C_{\mathrm{sec}} [W/K])\) without querying refrigerant \(c_p\). + fn secondary_inlet_capacity(&self, state: &StateSlice) -> Result<(f64, f64), ComponentError> { + let edges = self.inner.four_port_edges().ok_or_else(|| { + ComponentError::InvalidState( + "BphxExchanger: phase-change duty needs live four-port edges".into(), + ) + })?; + let (m_idx, p_idx, h_idx) = match self.geometry.exchanger_type { + // Condenser: cold = secondary. Evaporator (remapped): hot = secondary. + BphxType::Condenser => edges.cold_in, + BphxType::Evaporator => edges.hot_in, + BphxType::Generic => { + return Err(ComponentError::InvalidState( + "BphxExchanger: Generic type has no phase-change secondary stream".into(), + )); + } + }; + let max_idx = m_idx.max(p_idx).max(h_idx); + if max_idx >= state.len() { + return Err(ComponentError::InvalidStateDimensions { + expected: max_idx + 1, + actual: state.len(), + }); + } + let m = state[m_idx].max(0.0); + let p = state[p_idx]; + let h = state[h_idx]; + let (t, cp) = match self.geometry.exchanger_type { + BphxType::Condenser => ( + self.inner.cold_side_temperature(p, h)?, + self.inner.cold_side_cp(p, h)?, + ), + BphxType::Evaporator => ( + self.inner.hot_side_temperature(p, h)?, + self.inner.hot_side_cp(p, h)?, + ), + BphxType::Generic => unreachable!(), + }; + Ok((t, m * cp)) + } + + /// Replace sensible ε-NTU energy rows with phase-change duty. + /// + /// Layout (unchanged): `r0` = hot energy, `r1` = cold energy, then DP rows. + fn overwrite_energy_with_phase_change( + &self, + state: &StateSlice, + residuals: &mut ResidualVector, + ) -> Result<(), ComponentError> { + if residuals.len() < 2 { + return Err(ComponentError::InvalidResidualDimensions { + expected: 2, + actual: residuals.len(), + }); + } + let edges = self.inner.four_port_edges().ok_or_else(|| { + ComponentError::InvalidState( + "BphxExchanger: phase-change duty needs live four-port edges".into(), + ) + })?; + let (m_h, p_h_in, h_h_in) = edges.hot_in; + let (_, _, h_h_out) = edges.hot_out; + let (m_c, p_c_in, h_c_in) = edges.cold_in; + let (_, _, h_c_out) = edges.cold_out; + let max_idx = [m_h, p_h_in, h_h_in, h_h_out, m_c, p_c_in, h_c_in, h_c_out] + .into_iter() + .max() + .unwrap_or(0); + if max_idx >= state.len() { + return Err(ComponentError::InvalidStateDimensions { + expected: max_idx + 1, + actual: state.len(), + }); + } + + let (t_sec, c_sec) = self.secondary_inlet_capacity(state)?; + // Calibration redesign (WS-4): when a z_ua calibration variable is + // wired, the residual must track `state[z_ua]` (plain embedding), not + // the cached parameter — otherwise ∂Q/∂z_ua = 0 and the calibration + // Jacobian is singular. Mirrors exchanger.rs's `dynamic_f_ua` pattern. + let ua = match self.inner.calib_indices_ref().z_ua { + Some(z_idx) if z_idx < state.len() => self.inner.ua_nominal() * state[z_idx], + _ => self.ua(), + }; + let m_hot = state[m_h].max(0.0); + let m_cold = state[m_c].max(0.0); + + let q = match self.geometry.exchanger_type { + BphxType::Condenser => { + // hot = refrigerant, cold = secondary + let t_sat = self.tsat_k(state[p_h_in], self.inner.hot_fluid_id_str())?; + condenser_duty(ua, c_sec, t_sat, t_sec) + } + BphxType::Evaporator => { + // remapped: hot = secondary, cold = refrigerant + let t_sat = self.tsat_k(state[p_c_in], self.inner.cold_fluid_id_str())?; + evaporator_duty(ua, c_sec, t_sat, t_sec) + } + BphxType::Generic => 0.0, + }; + + // Same sign convention as EpsNtuModel::compute_residuals. + residuals[0] = m_hot * (state[h_h_in] - state[h_h_out]) - q; + residuals[1] = m_cold * (state[h_c_out] - state[h_c_in]) - q; + Ok(()) + } + + /// Finite-difference Jacobian through `self.compute_residuals` so energy + /// rows see phase-change duty (inner HeatExchanger FD would stay sensible). + fn fd_jacobian_via_self( + &self, + state: &StateSlice, + jacobian: &mut JacobianBuilder, + ) -> Result<(), ComponentError> { + let Some(edges) = self.inner.four_port_edges() else { + return Ok(()); + }; + let (m_h, p_h_in, h_h_in) = edges.hot_in; + let (_, p_h_out, h_h_out) = edges.hot_out; + let (m_c, p_c_in, h_c_in) = edges.cold_in; + let (_, p_c_out, h_c_out) = edges.cold_out; + + let mut cols = vec![ + m_h, p_h_in, h_h_in, p_h_out, h_h_out, m_c, p_c_in, h_c_in, p_c_out, h_c_out, + ]; + if let Some(z_idx) = self.inner.calib_indices_ref().z_dp { + cols.push(z_idx); + } + if let Some(z_idx) = self.inner.calib_indices_ref().z_ua { + cols.push(z_idx); + } + cols.sort_unstable(); + cols.dedup(); + cols.retain(|c| *c < state.len()); + + let n = self.n_equations(); + let compute = |s: &[f64]| -> Result, ComponentError> { + let mut r = vec![0.0; n]; + self.compute_residuals(s, &mut r)?; + Ok(r) + }; + + for &col in &cols { + let h = (state[col].abs() * 1e-6).max(1e-3); + let mut sp = state.to_vec(); + sp[col] += h; + let rp = compute(&sp)?; + let mut sm = state.to_vec(); + sm[col] -= h; + let rm = compute(&sm)?; + for row in 0..n { + let fd = (rp[row] - rm[row]) / (2.0 * h); + if fd.abs() > 1e-15 { + jacobian.add_entry(row, col, fd); + } + } + } + Ok(()) + } } impl Component for BphxExchanger { @@ -540,6 +729,23 @@ impl Component for BphxExchanger { state: &StateSlice, residuals: &mut ResidualVector, ) -> Result<(), ComponentError> { + if self.uses_phase_change_duty() { + // Skip HeatExchanger sensible ε-NTU (needs refrigerant cp, fails in + // the two-phase dome). Write Shah C*→0 energy + channel ΔP only. + let n = self.n_equations(); + if residuals.len() < n { + return Err(ComponentError::InvalidResidualDimensions { + expected: n, + actual: residuals.len(), + }); + } + for r in residuals.iter_mut().take(n) { + *r = 0.0; + } + self.overwrite_energy_with_phase_change(state, residuals)?; + self.overwrite_pressure_closures_with_dp(state, residuals)?; + return Ok(()); + } self.inner.compute_residuals(state, residuals)?; // Replace isobaric P_out − P_in with channel ΔP on both sides. self.overwrite_pressure_closures_with_dp(state, residuals) @@ -550,6 +756,10 @@ impl Component for BphxExchanger { state: &StateSlice, jacobian: &mut JacobianBuilder, ) -> Result<(), ComponentError> { + if self.uses_phase_change_duty() { + // FD must call this Component's residuals (phase-change energy + DP). + return self.fd_jacobian_via_self(state, jacobian); + } self.inner.jacobian_entries(state, jacobian)?; // Inner already wrote ∂r/∂P_out=+1, ∂r/∂P_in=−1; add ṁ and z_dp terms. self.append_pressure_dp_jacobian(state, jacobian) @@ -587,6 +797,80 @@ impl Component for BphxExchanger { self.inner.energy_transfers(state) } + fn measure_output(&self, kind: crate::MeasuredOutput, state: &StateSlice) -> Option { + use crate::MeasuredOutput::*; + let edges = self.inner.four_port_edges()?; + let (m_h, p_h_in, h_h_in) = edges.hot_in; + let (_, p_h_out, h_h_out) = edges.hot_out; + let (m_c, p_c_in, h_c_in) = edges.cold_in; + let (_, _, h_c_out) = edges.cold_out; + let max_idx = [ + m_h, p_h_in, h_h_in, p_h_out, h_h_out, m_c, p_c_in, h_c_in, h_c_out, + ] + .into_iter() + .max()?; + if max_idx >= state.len() { + return None; + } + // Refrigerant side: hot for a condenser, cold (remapped) for an evaporator. + let (p_ref_in, p_ref_out, h_ref_out, m_ref, ref_fluid) = match self.geometry.exchanger_type + { + BphxType::Condenser => (p_h_in, p_h_out, h_h_out, m_h, self.inner.hot_fluid_id_str()), + BphxType::Evaporator => (p_c_in, p_c_in, h_c_out, m_c, self.inner.cold_fluid_id_str()), + BphxType::Generic => return None, + }; + match kind { + SaturationTemperature => self.tsat_k(state[p_ref_in], ref_fluid).ok(), + Superheat | Subcooling => { + let backend = self.fluid_backend.as_ref()?; + let tsat = self.tsat_k(state[p_ref_out], ref_fluid).ok()?; + let t = backend + .property( + FluidId::new(ref_fluid), + Property::Temperature, + FluidState::from_ph( + entropyk_core::Pressure::from_pascals(state[p_ref_out]), + entropyk_core::Enthalpy::from_joules_per_kg(state[h_ref_out]), + ), + ) + .ok()?; + if !t.is_finite() || t <= 0.0 { + return None; + } + match kind { + Superheat => Some(t - tsat), + _ => Some(tsat - t), + } + } + Capacity | HeatTransferRate => { + let q_hot = state[m_h].abs() * (state[h_h_in] - state[h_h_out]).abs(); + let q_cold = state[m_c].abs() * (state[h_c_out] - state[h_c_in]).abs(); + match (q_hot.is_finite(), q_cold.is_finite()) { + (true, true) => Some(0.5 * (q_hot + q_cold)), + (true, false) => Some(q_hot), + (false, true) => Some(q_cold), + _ => None, + } + } + MassFlowRate => Some(state[m_ref].abs()), + Pressure => Some(state[p_ref_in]), + Temperature => { + let backend = self.fluid_backend.as_ref()?; + backend + .property( + FluidId::new(ref_fluid), + Property::Temperature, + FluidState::from_ph( + entropyk_core::Pressure::from_pascals(state[p_ref_out]), + entropyk_core::Enthalpy::from_joules_per_kg(state[h_ref_out]), + ), + ) + .ok() + .filter(|t| t.is_finite() && *t > 0.0) + } + } + } + fn set_fluid_backend_from_builder( &mut self, backend: std::sync::Arc, @@ -983,8 +1267,8 @@ mod tests { // State layout per port triple (ṁ, P, h): // hot_in 0..3, hot_out 3..6, cold_in 6..9, cold_out 9..12 - let mut hx = BphxExchanger::new(test_geometry()) - .with_fluid_backend(Arc::new(TestBackend::new())); + let mut hx = + BphxExchanger::new(test_geometry()).with_fluid_backend(Arc::new(TestBackend::new())); hx.set_hot_fluid("Water"); hx.set_cold_fluid("Water"); hx.set_port_context(&[ @@ -1027,8 +1311,8 @@ mod tests { ); // z_dp scale: half calib → half residual contribution at same P_out=P_in - let mut hx_half = BphxExchanger::new(test_geometry()) - .with_fluid_backend(Arc::new(TestBackend::new())); + let mut hx_half = + BphxExchanger::new(test_geometry()).with_fluid_backend(Arc::new(TestBackend::new())); hx_half.set_hot_fluid("Water"); hx_half.set_cold_fluid("Water"); hx_half.set_port_context(&[ @@ -1078,8 +1362,8 @@ mod tests { state[10] = 200_000.0; state[11] = 60_000.0; - let mut simp = BphxExchanger::new(test_geometry()) - .with_fluid_backend(Arc::new(TestBackend::new())); + let mut simp = + BphxExchanger::new(test_geometry()).with_fluid_backend(Arc::new(TestBackend::new())); simp.set_hot_fluid("Water"); simp.set_cold_fluid("Water"); simp.set_port_context(&ports); @@ -1102,4 +1386,110 @@ mod tests { "Martin residual should differ from Simplified" ); } + + #[test] + fn test_bphx_condenser_energy_uses_phase_change_duty() { + use entropyk_fluids::TestBackend; + use std::sync::Arc; + + // Condenser geometry → Shah C*→0 path (not sensible ṁ·cp). + let geo = test_geometry().with_exchanger_type(BphxType::Condenser); + let mut hx = BphxExchanger::new(geo).with_fluid_backend(Arc::new(TestBackend::new())); + hx.set_hot_fluid("R134a"); + hx.set_cold_fluid("Water"); + hx.set_port_context(&[ + Some((0, 1, 2)), + Some((3, 4, 5)), + Some((6, 7, 8)), + Some((9, 10, 11)), + ]); + + let mut state = vec![0.0; 12]; + // Refrigerant ~10 bar, superheated-ish enthalpy table values + state[0] = 0.02; + state[1] = 1_000_000.0; + state[2] = 430_000.0; + state[3] = 0.02; + state[4] = 1_000_000.0; + state[5] = 250_000.0; + // Water secondary ~2 bar, ~30 °C + state[6] = 0.40; + state[7] = 200_000.0; + state[8] = 125_000.0; + state[9] = 0.40; + state[10] = 200_000.0; + state[11] = 140_000.0; + + assert!(hx.uses_phase_change_duty()); + let n = hx.n_equations(); + let mut residuals = vec![0.0; n]; + hx.compute_residuals(&state, &mut residuals).unwrap(); + + let t_sat = hx.tsat_k(state[1], "R134a").unwrap(); + let (t_sec, c_sec) = hx.secondary_inlet_capacity(&state).unwrap(); + let q = condenser_duty(hx.ua(), c_sec, t_sat, t_sec); + let q_hot = state[0] * (state[2] - state[5]); + let expected_r0 = q_hot - q; + assert!( + (residuals[0] - expected_r0).abs() < 1e-4 * expected_r0.abs().max(1.0), + "energy r0 must match condenser_duty: got {} expected {} (Q={}, Tsat={:.2}K)", + residuals[0], + expected_r0, + q, + t_sat + ); + // Approach must use Tsat, not gas temperature — Q grows if Tsat rises. + let q_higher = condenser_duty(hx.ua(), c_sec, t_sat + 10.0, t_sec); + assert!(q_higher > q, "higher Tsat must increase condenser duty"); + } + + #[test] + fn test_bphx_evaporator_energy_uses_phase_change_duty() { + use entropyk_fluids::TestBackend; + use std::sync::Arc; + + let geo = test_geometry().with_exchanger_type(BphxType::Evaporator); + let mut hx = BphxExchanger::new(geo).with_fluid_backend(Arc::new(TestBackend::new())); + // Same remapping as BphxEvaporator: hot = secondary, cold = refrigerant. + hx.set_hot_fluid("Water"); + hx.set_cold_fluid("R134a"); + hx.set_port_context(&[ + Some((0, 1, 2)), + Some((3, 4, 5)), + Some((6, 7, 8)), + Some((9, 10, 11)), + ]); + + let mut state = vec![0.0; 12]; + // Water (hot) ~12 °C + state[0] = 0.50; + state[1] = 300_000.0; + state[2] = 50_000.0; + state[3] = 0.50; + state[4] = 300_000.0; + state[5] = 40_000.0; + // Refrigerant (cold) ~3.5 bar + state[6] = 0.02; + state[7] = 350_000.0; + state[8] = 250_000.0; + state[9] = 0.02; + state[10] = 350_000.0; + state[11] = 400_000.0; + + assert!(hx.uses_phase_change_duty()); + let mut residuals = vec![0.0; hx.n_equations()]; + hx.compute_residuals(&state, &mut residuals).unwrap(); + + let t_sat = hx.tsat_k(state[7], "R134a").unwrap(); + let (t_sec, c_sec) = hx.secondary_inlet_capacity(&state).unwrap(); + let q = evaporator_duty(hx.ua(), c_sec, t_sat, t_sec); + let q_cold = state[6] * (state[11] - state[8]); + let expected_r1 = q_cold - q; + assert!( + (residuals[1] - expected_r1).abs() < 1e-4 * expected_r1.abs().max(1.0), + "energy r1 must match evaporator_duty: got {} expected {}", + residuals[1], + expected_r1 + ); + } } diff --git a/crates/components/src/heat_exchanger/exchanger.rs b/crates/components/src/heat_exchanger/exchanger.rs index 85b15ca..1122423 100644 --- a/crates/components/src/heat_exchanger/exchanger.rs +++ b/crates/components/src/heat_exchanger/exchanger.rs @@ -588,6 +588,34 @@ impl HeatExchanger { }) } + /// Hot-side temperature [K] at `(P, h)` (live backend). + pub(crate) fn hot_side_temperature( + &self, + p_pa: f64, + h_jkg: f64, + ) -> Result { + self.hot_temperature(p_pa, h_jkg) + } + + /// Hot-side \(c_p\) [J/(kg·K)] at `(P, h)`. + pub(crate) fn hot_side_cp(&self, p_pa: f64, h_jkg: f64) -> Result { + self.hot_cp(p_pa, h_jkg) + } + + /// Cold-side temperature [K] at `(P, h)`. + pub(crate) fn cold_side_temperature( + &self, + p_pa: f64, + h_jkg: f64, + ) -> Result { + self.cold_temperature(p_pa, h_jkg) + } + + /// Cold-side \(c_p\) [J/(kg·K)] at `(P, h)`. + pub(crate) fn cold_side_cp(&self, p_pa: f64, h_jkg: f64) -> Result { + self.cold_cp(p_pa, h_jkg) + } + fn live_state_required_error(&self) -> ComponentError { ComponentError::InvalidState(format!( "{} requires live four-port edge state (hot_inlet, hot_outlet, cold_inlet, cold_outlet); inlet-only boundary conditions cannot define outlet states", diff --git a/crates/components/src/heat_exchanger/mod.rs b/crates/components/src/heat_exchanger/mod.rs index ec7239a..15c510a 100644 --- a/crates/components/src/heat_exchanger/mod.rs +++ b/crates/components/src/heat_exchanger/mod.rs @@ -72,6 +72,7 @@ pub mod lmtd; pub mod mchx_condenser_coil; pub mod model; pub mod moving_boundary_hx; +pub mod phase_change_entu; pub mod pool_boiling; pub mod sat_domain; pub mod shell_and_tube; @@ -119,6 +120,9 @@ pub use gas_cooler::{is_supercritical, pettersen_htc, GasCooler, GasCoolerInput} pub use lmtd::{FlowConfiguration, LmtdModel}; pub use mchx_condenser_coil::MchxCondenserCoil; pub use model::HeatTransferModel; +pub use phase_change_entu::{ + condenser_duty, d_eps_c_d_c, evaporator_duty, phase_change_effectiveness, +}; pub use pool_boiling::{ assess_cooper_domain, cooper_1984, cooper_metadata, flooded_shell_htc, mostinski_1963, palen_bundle_factor, thome_robinson_oil_factor, ua_from_two_side_htc, PoolBoilingInput, diff --git a/crates/components/src/heat_exchanger/phase_change_entu.rs b/crates/components/src/heat_exchanger/phase_change_entu.rs new file mode 100644 index 0000000..04f6e35 --- /dev/null +++ b/crates/components/src/heat_exchanger/phase_change_entu.rs @@ -0,0 +1,87 @@ +//! Phase-change ε-NTU helpers (Shah / Incropera \(C^*\to 0\)). +//! +//! For condensers and evaporators where the refrigerant changes phase at +//! (ideally) constant \(T_{\mathrm{sat}}(P)\), the refrigerant heat-capacity +//! rate approaches infinity. Then \(C_{\min}=C_{\mathrm{sec}}\) and +//! \(\varepsilon = 1 - \exp(-\mathrm{UA}/C_{\mathrm{sec}})\). +//! +//! See Shah & Sekulić, *Fundamentals of Heat Exchanger Design* (Wiley 2003), §3.3.2. + +/// Effectiveness for a single-stream (phase-change) exchanger: +/// \(\varepsilon = 1 - \exp(-\mathrm{UA}/C_{\mathrm{sec}})\). +#[inline] +pub fn phase_change_effectiveness(ua: f64, c_sec: f64) -> f64 { + if c_sec <= 1e-10 || ua <= 0.0 { + return 0.0; + } + 1.0 - (-ua / c_sec).exp() +} + +/// \(g'(C)\) where \(g(C)=C\cdot\varepsilon(C)=C\cdot(1-e^{-\mathrm{UA}/C})\). +/// +/// Used for \(\partial Q/\partial\dot{m}_{\mathrm{sec}}\). +#[inline] +pub fn d_eps_c_d_c(ua: f64, c_sec: f64) -> f64 { + if c_sec <= 1e-10 || ua <= 0.0 { + return 0.0; + } + let e = (-ua / c_sec).exp(); + (1.0 - e) - (ua / c_sec) * e +} + +/// Condenser duty [W]: \(Q = \varepsilon\,C_{\mathrm{sec}}\,(T_{\mathrm{sat}}-T_{\mathrm{sec,in}})\). +/// +/// Positive \(Q\) = heat rejected by the refrigerant into the secondary. +#[inline] +pub fn condenser_duty(ua: f64, c_sec: f64, t_sat_k: f64, t_sec_in_k: f64) -> f64 { + let eps = phase_change_effectiveness(ua, c_sec); + eps * c_sec * (t_sat_k - t_sec_in_k) +} + +/// Evaporator duty [W]: \(Q = \varepsilon\,C_{\mathrm{sec}}\,(T_{\mathrm{sec,in}}-T_{\mathrm{sat}})\). +/// +/// Positive \(Q\) = heat absorbed by the refrigerant from the secondary. +#[inline] +pub fn evaporator_duty(ua: f64, c_sec: f64, t_sat_k: f64, t_sec_in_k: f64) -> f64 { + let eps = phase_change_effectiveness(ua, c_sec); + eps * c_sec * (t_sec_in_k - t_sat_k) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn effectiveness_cr_zero_matches_textbook() { + let ua = 2500.0; + let c_sec = 0.4 * 4180.0; + let eps = phase_change_effectiveness(ua, c_sec); + let expected = 1.0 - (-ua / c_sec).exp(); + assert!((eps - expected).abs() < 1e-12); + assert!(eps > 0.7 && eps < 1.0); + } + + #[test] + fn condenser_duty_zero_when_tsat_equals_tsec() { + let q = condenser_duty(5000.0, 1000.0, 303.15, 303.15); + assert!(q.abs() < 1e-9); + } + + #[test] + fn evaporator_duty_positive_when_water_warmer_than_tsat() { + let q = evaporator_duty(2000.0, 2000.0, 278.15, 285.15); + assert!(q > 0.0); + } + + #[test] + fn d_eps_c_matches_fd() { + let ua = 3000.0; + let c = 1500.0; + let h = 1e-3 * c; + let gp = (c + h) * phase_change_effectiveness(ua, c + h); + let gm = (c - h) * phase_change_effectiveness(ua, c - h); + let fd = (gp - gm) / (2.0 * h); + let an = d_eps_c_d_c(ua, c); + assert!((an - fd).abs() / fd.abs().max(1e-9) < 1e-5); + } +} diff --git a/crates/components/src/lib.rs b/crates/components/src/lib.rs index b134c19..d20e867 100644 --- a/crates/components/src/lib.rs +++ b/crates/components/src/lib.rs @@ -97,6 +97,7 @@ pub mod params; pub mod pipe; pub mod polynomials; pub mod port; +pub mod probe; pub mod pump; pub mod python_components; pub mod refrigerant_boundary; @@ -150,6 +151,7 @@ pub use port::{ validate_port_continuity, Connected, ConnectedPort, ConnectionError, Disconnected, FluidId, Port, PortKind, }; +pub use probe::{Probe, ProbeMeasure}; pub use pump::{Pump, PumpCurves}; pub use python_components::{ PyAirSinkReal, PyAirSourceReal, PyBrineSinkReal, PyBrineSourceReal, PyCompressorReal, diff --git a/crates/components/src/probe.rs b/crates/components/src/probe.rs new file mode 100644 index 0000000..9c03150 --- /dev/null +++ b/crates/components/src/probe.rs @@ -0,0 +1,907 @@ +//! Probe — zero-residual measurement tap on a connection line. +//! +//! Per the calibration redesign (see +//! `_bmad-output/implementation-artifacts/calibration-probe-node-plain-embedding.md`), +//! **every** calibration measurement lives on a `Probe` node the user drops +//! onto a wire. The calibrated z-factor (`z_ua`, `z_dp`, `z_flow`, ...) on a +//! component is freed, and a `Probe` elsewhere on the graph supplies the +//! matching measurement via its [`Component::measure_output`] override. +//! +//! ## Design +//! +//! - **Zero residual, zero Jacobian** — `n_equations() == 0`. The Probe is a +//! measurement tap, not a physical element. The calibration `Constraint` +//! (one residual, one unknown z-factor) closes the system. +//! - **Two ports `inlet`/`outlet`** — same typestate pattern as `Pipe`/`Pump`/ +//! `Node`. The Probe splices into an edge A→B (becomes A→Probe→B) exactly +//! like `Pipe` does in `apps/web/src/lib/edgeInsert.ts`, so the UI drop-on- +//! wire interaction reuses `insertOnEdge`. +//! - **Reads live edge state** via [`Component::set_system_context`]. The +//! Probe stores the `(ṁ, P, h)` triples for both its incident edges (the +//! two halves of the spliced wire) and uses them in `measure_output`. +//! +//! ## Measure kinds +//! +//! [`ProbeMeasure`] is the exhaustive list the user named. SST/SDT/DGT/DSH/ +//! SH/SC all derive from `(P, T)` on the same edge via CoolProp — the user +//! picks the semantic kind for clarity, the math is P/T/Tsat: +//! +//! | Kind | Formula | +//! |------|---------| +//! | [`Sst`](ProbeMeasure::Sst) / [`Sdt`](ProbeMeasure::Sdt) | `Tsat(P_edge)` | +//! | [`Dgt`](ProbeMeasure::Dgt) / [`T`](ProbeMeasure::T) | `T_edge` | +//! | [`Dsh`](ProbeMeasure::Dsh) / [`Sh`](ProbeMeasure::Sh) | `T_edge − Tsat(P_edge)` | +//! | [`Sc`](ProbeMeasure::Sc) | `Tsat(P_edge) − T_edge` | +//! | [`P`](ProbeMeasure::P) | `P_edge` | +//! | [`MassFlow`](ProbeMeasure::MassFlow) | `ṁ_edge` | +//! | [`Enthalpy`](ProbeMeasure::Enthalpy) | `h_edge` | +//! | [`Capacity`](ProbeMeasure::Capacity) | `ṁ_edge · |h_inlet − h_outlet|` | +//! +//! `Tsat(P)` reuses the `tsat_k` pattern from +//! `crates/components/src/heat_exchanger/bphx_exchanger.rs:543-559` (clamps +//! the query pressure into the detected saturation domain via +//! `sat_domain::clamp_to_saturation_domain`). + +use crate::heat_exchanger::sat_domain; +use crate::port::{Connected, Disconnected, Port}; +use crate::state_machine::StateManageable; +use crate::{ + CircuitId, Component, ComponentError, ConnectedPort, JacobianBuilder, MeasuredOutput, + OperationalState, ResidualVector, StateSlice, +}; +use entropyk_core::{CalibIndices, MassFlow, Power}; +use entropyk_fluids::{FluidBackend, FluidId as BackendFluidId, FluidState, Property, Quality}; +use std::marker::PhantomData; +use std::sync::Arc; + +/// Semantic quantity a [`Probe`] extracts from its edge. +/// +/// Exhaustive list per the user mandate (see module docs). SST/SDT and DSH/SH +/// share formulas by design — the user picks the label that matches the +/// physical sensor location (suction vs discharge). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum ProbeMeasure { + /// Saturated suction temperature = `Tsat(P_edge)`. + Sst, + /// Saturated discharge temperature = `Tsat(P_edge)` (same formula, label only). + Sdt, + /// Discharge gas temperature = `T_edge` (raw temperature at the probe). + Dgt, + /// Discharge superheat = `T_edge − Tsat(P_edge)`. + Dsh, + /// (Suction) superheat = `T_edge − Tsat(P_edge)` (same formula as Dsh). + Sh, + /// Subcooling = `Tsat(P_edge) − T_edge`. + Sc, + /// Raw temperature = `T_edge`. + T, + /// Raw pressure = `P_edge`. + P, + /// Mass flow rate = `ṁ_edge`. + MassFlow, + /// Heat duty = `ṁ · |Δh|` across the probe (see [`Probe`] module docs). + Capacity, + /// Raw specific enthalpy = `h_edge`. + Enthalpy, +} + +impl ProbeMeasure { + /// Parses a measure-kind string from the JSON config (case-insensitive). + /// Accepts the canonical names plus a few common aliases. + pub fn parse(s: &str) -> Result { + let lower = s.trim().to_ascii_lowercase(); + let m = match lower.as_str() { + "sst" => ProbeMeasure::Sst, + "sdt" => ProbeMeasure::Sdt, + "dgt" => ProbeMeasure::Dgt, + "dsh" => ProbeMeasure::Dsh, + "sh" | "superheat" => ProbeMeasure::Sh, + "sc" | "subcooling" => ProbeMeasure::Sc, + "t" | "temperature" => ProbeMeasure::T, + "p" | "pressure" => ProbeMeasure::P, + "massflow" | "mass_flow" | "mass flow" | "massflowrate" | "mass_flow_rate" => { + ProbeMeasure::MassFlow + } + "capacity" | "heatduty" | "heat_transfer_rate" | "heatrate" | "duty" => { + ProbeMeasure::Capacity + } + "enthalpy" | "h" => ProbeMeasure::Enthalpy, + other => { + return Err(format!( + "unknown Probe measure '{other}' (expected one of: SST, SDT, DGT, DSH, SH, SC, T, P, MassFlow, Capacity, Enthalpy)" + )) + } + }; + Ok(m) + } + + /// Returns the canonical JSON name for this measure (round-trips through `parse`). + pub fn as_str(&self) -> &'static str { + match self { + ProbeMeasure::Sst => "SST", + ProbeMeasure::Sdt => "SDT", + ProbeMeasure::Dgt => "DGT", + ProbeMeasure::Dsh => "DSH", + ProbeMeasure::Sh => "SH", + ProbeMeasure::Sc => "SC", + ProbeMeasure::T => "T", + ProbeMeasure::P => "P", + ProbeMeasure::MassFlow => "MassFlow", + ProbeMeasure::Capacity => "Capacity", + ProbeMeasure::Enthalpy => "Enthalpy", + } + } + + /// Maps this measure to the solver-side [`MeasuredOutput`] the calibration + /// constraint will ask for. SST/SDT → `SaturationTemperature`; DSH/SH → + /// `Superheat`; SC → `Subcooling`; DGT/T → `Temperature`; P → `Pressure`; + /// MassFlow → `MassFlowRate`; Capacity → `Capacity`; Enthalpy → + /// `HeatTransferRate` is NOT right, so Enthalpy returns `Temperature` as a + /// placeholder (the Probe still serves it via `measure_output(Temperature)` + /// returning T, and Enthalpy is read directly via the measure-output hook + /// below — see [`Probe::measure_output`] which honours the configured kind). + pub fn to_measured_output(&self) -> MeasuredOutput { + match self { + ProbeMeasure::Sst | ProbeMeasure::Sdt => MeasuredOutput::SaturationTemperature, + ProbeMeasure::Dsh | ProbeMeasure::Sh => MeasuredOutput::Superheat, + ProbeMeasure::Sc => MeasuredOutput::Subcooling, + ProbeMeasure::Dgt | ProbeMeasure::T => MeasuredOutput::Temperature, + ProbeMeasure::P => MeasuredOutput::Pressure, + ProbeMeasure::MassFlow => MeasuredOutput::MassFlowRate, + ProbeMeasure::Capacity => MeasuredOutput::Capacity, + // No dedicated Enthalpy variant in MeasuredOutput — collapse onto + // HeatTransferRate so the constraint resolves; measure_output + // detects the configured Enthalpy kind and returns h_edge. + ProbeMeasure::Enthalpy => MeasuredOutput::HeatTransferRate, + } + } +} + +/// A zero-residual measurement tap on a connection line. +/// +/// See the module docs for the calibration redesign context. +/// +/// Like `Pipe` and `Pump`, `Probe` uses a typestate (`Disconnected` → +/// `Connected`) so the compiler enforces port connection before the component +/// is added to a system. +#[derive(Clone)] +pub struct Probe { + measure: ProbeMeasure, + fluid_id_str: String, + port_inlet: Port, + port_outlet: Port, + fluid_backend: Option>, + circuit_id: CircuitId, + operational_state: OperationalState, + calib_indices: CalibIndices, + // Live edge-state indices populated by `set_system_context`. + inlet_m_idx: Option, + inlet_p_idx: Option, + inlet_h_idx: Option, + outlet_m_idx: Option, + outlet_p_idx: Option, + outlet_h_idx: Option, + _state: PhantomData, +} + +impl std::fmt::Debug for Probe { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("Probe") + .field("measure", &self.measure) + .field("fluid", &self.fluid_id_str) + .field("has_backend", &self.fluid_backend.is_some()) + .field("circuit_id", &self.circuit_id) + .field("operational_state", &self.operational_state) + .finish() + } +} + +impl Probe { + /// Creates a new disconnected Probe. + /// + /// `fluid` is used for `Tsat`/`SH`/`SC` derivation via the configured + /// fluid backend. The backend itself is injected later by + /// [`Component::set_fluid_backend_from_builder`] when the system is built. + pub fn new( + measure: ProbeMeasure, + fluid: impl Into, + port_inlet: Port, + port_outlet: Port, + ) -> Self { + Self { + measure, + fluid_id_str: fluid.into(), + port_inlet, + port_outlet, + fluid_backend: None, + circuit_id: CircuitId::default(), + operational_state: OperationalState::default(), + calib_indices: CalibIndices::default(), + inlet_m_idx: None, + inlet_p_idx: None, + inlet_h_idx: None, + outlet_m_idx: None, + outlet_p_idx: None, + outlet_h_idx: None, + _state: PhantomData, + } + } + + /// Attaches a fluid backend for property queries (Tsat, T from P,h). + pub fn with_fluid_backend(mut self, backend: Arc) -> Self { + self.fluid_backend = Some(backend); + self + } + + /// Connects the Probe to its inlet and outlet ports, transitioning to + /// `Probe`. Mirrors `Pipe::connect` / `Pump::connect`. + pub fn connect( + self, + inlet: Port, + outlet: Port, + ) -> Result, ComponentError> { + let (p_in, _) = self + .port_inlet + .connect(inlet) + .map_err(|e| ComponentError::InvalidState(e.to_string()))?; + let (p_out, _) = self + .port_outlet + .connect(outlet) + .map_err(|e| ComponentError::InvalidState(e.to_string()))?; + Ok(Probe { + measure: self.measure, + fluid_id_str: self.fluid_id_str, + port_inlet: p_in, + port_outlet: p_out, + fluid_backend: self.fluid_backend, + circuit_id: self.circuit_id, + operational_state: self.operational_state, + calib_indices: self.calib_indices, + inlet_m_idx: self.inlet_m_idx, + inlet_p_idx: self.inlet_p_idx, + inlet_h_idx: self.inlet_h_idx, + outlet_m_idx: self.outlet_m_idx, + outlet_p_idx: self.outlet_p_idx, + outlet_h_idx: self.outlet_h_idx, + _state: PhantomData, + }) + } +} + +impl Probe { + /// Returns the configured measure kind. + pub fn measure(&self) -> ProbeMeasure { + self.measure + } + + /// Returns the fluid identifier string the Probe derives properties for. + pub fn fluid_id_str(&self) -> &str { + &self.fluid_id_str + } + + /// Returns the inlet port. + pub fn port_inlet(&self) -> &Port { + &self.port_inlet + } + + /// Returns the outlet port. + pub fn port_outlet(&self) -> &Port { + &self.port_outlet + } + + /// Returns both ports as a slice. + pub fn get_ports_slice(&self) -> [&Port; 2] { + [&self.port_inlet, &self.port_outlet] + } + + /// \(T_{\mathrm{sat}}(P)\) [K] for the configured fluid, domain-clamped. + /// Mirrors `BphxExchanger::tsat_k` (`bphx_exchanger.rs:544-559`). + fn tsat_k(&self, p_pa: f64) -> Result { + let backend = self.fluid_backend.as_ref().ok_or_else(|| { + ComponentError::CalculationFailed( + "Probe: FluidBackend required for Tsat/SH/SC derivation".into(), + ) + })?; + let p_pa = + sat_domain::clamp_to_saturation_domain(backend, &self.fluid_id_str, p_pa).unwrap_or(p_pa); + backend + .property( + BackendFluidId::new(&self.fluid_id_str), + Property::Temperature, + FluidState::from_px( + entropyk_core::Pressure::from_pascals(p_pa), + Quality::new(0.5), + ), + ) + .map_err(ComponentError::from_fluid_error) + } + + /// Temperature at the probe edge from `(P, h)` via the fluid backend. + fn temperature_k(&self, p_pa: f64, h_j_kg: f64) -> Result { + let backend = self.fluid_backend.as_ref().ok_or_else(|| { + ComponentError::CalculationFailed( + "Probe: FluidBackend required for temperature derivation".into(), + ) + })?; + let t = backend + .property( + BackendFluidId::new(&self.fluid_id_str), + Property::Temperature, + FluidState::from_ph( + entropyk_core::Pressure::from_pascals(p_pa), + entropyk_core::Enthalpy::from_joules_per_kg(h_j_kg), + ), + ) + .map_err(ComponentError::from_fluid_error)?; + if !t.is_finite() || t <= 0.0 { + return Err(ComponentError::CalculationFailed(format!( + "Probe: non-physical temperature {t} for P={p_pa} Pa, h={h_j_kg} J/kg" + ))); + } + Ok(t) + } + + /// Reads `(ṁ, P, h)` of the inlet edge from the global state, or `None` + /// when the indices have not been wired yet. + fn inlet_state(&self, state: &StateSlice) -> Option<(f64, f64, f64)> { + let (m, p, h) = ( + self.inlet_m_idx?, + self.inlet_p_idx?, + self.inlet_h_idx?, + ); + if m >= state.len() || p >= state.len() || h >= state.len() { + return None; + } + Some((state[m], state[p], state[h])) + } + + /// Reads the outlet-edge `(P, h)` from the global state, or `None`. + fn outlet_state(&self, state: &StateSlice) -> Option<(f64, f64)> { + let (p, h) = (self.outlet_p_idx?, self.outlet_h_idx?); + if p >= state.len() || h >= state.len() { + return None; + } + Some((state[p], state[h])) + } +} + +impl Component for Probe { + fn set_system_context( + &mut self, + _state_offset: usize, + external_edge_state_indices: &[(usize, usize, usize)], + ) { + // The Probe is a 2-port splice; layout mirrors Pipe/Pump: + // [0] = incoming edge (upstream→probe), [1] = outgoing edge (probe→downstream). + // The Probe measures at its location; for SST/SDT/DGT/DSH/SH/SC/T/P we + // use the inlet edge (the edge the user dropped the probe onto). For + // Capacity we use both edges (ṁ·|Δh| across the splice). + if let Some(&(m, p, h)) = external_edge_state_indices.first() { + self.inlet_m_idx = Some(m); + self.inlet_p_idx = Some(p); + self.inlet_h_idx = Some(h); + } + if let Some(&(m, p, h)) = external_edge_state_indices.get(1) { + self.outlet_m_idx = Some(m); + self.outlet_p_idx = Some(p); + self.outlet_h_idx = Some(h); + } + } + + fn compute_residuals( + &self, + state: &StateSlice, + residuals: &mut ResidualVector, + ) -> Result<(), ComponentError> { + // Two topology residuals any 2-port splice must impose to keep the + // DoF balance closed (mirrors `Anchor` and edge-coupled `Pipe`): + // r0 = P_out − P_in (zero-resistance pressure pass-through) + // r1 = h_out − h_in (adiabatic enthalpy pass-through) + // These are NOT calibration residuals — they are the continuity that + // makes the Probe a transparent tap. The calibration `Constraint` + // (measure − target) supplies the only calibration-related residual. + if residuals.len() < 2 { + return Err(ComponentError::InvalidResidualDimensions { + expected: 2, + actual: residuals.len(), + }); + } + let (in_p_idx, in_h_idx, out_p_idx, out_h_idx) = match ( + self.inlet_p_idx, + self.inlet_h_idx, + self.outlet_p_idx, + self.outlet_h_idx, + ) { + (Some(a), Some(b), Some(c), Some(d)) => (a, b, c, d), + _ => { + return Err(ComponentError::InvalidState( + "Probe requires live inlet and outlet edge state indices".to_string(), + )); + } + }; + let max_idx = in_p_idx.max(in_h_idx).max(out_p_idx).max(out_h_idx); + if max_idx >= state.len() { + return Err(ComponentError::InvalidStateDimensions { + expected: max_idx + 1, + actual: state.len(), + }); + } + residuals[0] = state[out_p_idx] - state[in_p_idx]; + residuals[1] = state[out_h_idx] - state[in_h_idx]; + Ok(()) + } + + fn jacobian_entries( + &self, + _state: &StateSlice, + jacobian: &mut JacobianBuilder, + ) -> Result<(), ComponentError> { + // ∂r0/∂P_out = +1, ∂r0/∂P_in = −1 ; ∂r1/∂h_out = +1, ∂r1/∂h_in = −1 + let (in_p_idx, in_h_idx, out_p_idx, out_h_idx) = match ( + self.inlet_p_idx, + self.inlet_h_idx, + self.outlet_p_idx, + self.outlet_h_idx, + ) { + (Some(a), Some(b), Some(c), Some(d)) => (a, b, c, d), + _ => { + return Err(ComponentError::InvalidState( + "Probe Jacobian requires live inlet and outlet edge state indices".to_string(), + )); + } + }; + jacobian.add_entry(0, out_p_idx, 1.0); + jacobian.add_entry(0, in_p_idx, -1.0); + jacobian.add_entry(1, out_h_idx, 1.0); + jacobian.add_entry(1, in_h_idx, -1.0); + Ok(()) + } + + fn equation_roles(&self) -> Vec { + vec![ + crate::EquationRole::Continuity { quantity: "P" }, + crate::EquationRole::Continuity { quantity: "h" }, + ] + } + + fn n_equations(&self) -> usize { + // 2 topology-continuity residuals (P, h) — see `compute_residuals`. + // NOT a calibration residual: the calibration Constraint supplies that. + 2 + } + + fn get_ports(&self) -> &[ConnectedPort] { + &[] + } + + fn port_names(&self) -> Vec { + vec!["inlet".to_string(), "outlet".to_string()] + } + + fn flow_paths(&self) -> Vec<(usize, usize)> { + // Single series path through the tap. + vec![(0, 1)] + } + + fn port_mass_flows( + &self, + state: &StateSlice, + ) -> Result, ComponentError> { + let m_in = match self.inlet_m_idx { + Some(idx) if idx < state.len() => state[idx], + _ => 0.0, + }; + let m_out = match self.outlet_m_idx { + Some(idx) if idx < state.len() => state[idx], + _ => m_in, + }; + Ok(vec![ + MassFlow::from_kg_per_s(m_in), + MassFlow::from_kg_per_s(-m_out), + ]) + } + + fn port_enthalpies( + &self, + state: &StateSlice, + ) -> Result, ComponentError> { + let h_in = match self.inlet_h_idx { + Some(idx) if idx < state.len() => state[idx], + _ => self.port_inlet.enthalpy().to_joules_per_kg(), + }; + let h_out = match self.outlet_h_idx { + Some(idx) if idx < state.len() => state[idx], + _ => self.port_outlet.enthalpy().to_joules_per_kg(), + }; + Ok(vec![ + entropyk_core::Enthalpy::from_joules_per_kg(h_in), + entropyk_core::Enthalpy::from_joules_per_kg(h_out), + ]) + } + + fn set_calib_indices(&mut self, indices: CalibIndices) { + // A Probe owns no z-factors; the slot is kept only for trait + // compatibility with the System finalize path. + self.calib_indices = indices; + } + + fn set_fluid_backend_from_builder( + &mut self, + backend: Arc, + ) { + if self.fluid_backend.is_none() { + self.fluid_backend = Some(Arc::clone(&backend)); + } + } + + fn energy_transfers(&self, _state: &StateSlice) -> Option<(Power, Power)> { + // A measurement tap is adiabatic. + Some((Power::from_watts(0.0), Power::from_watts(0.0))) + } + + fn measure_output(&self, kind: MeasuredOutput, state: &StateSlice) -> Option { + // Resolve which edge we measure on. SST/SDT use the inlet edge's + // pressure; DGT/T/P/MassFlow/Enthalpy also use the inlet edge (the + // physical sensor location). For DSH/SH/SC we need (P, T) on the same + // edge. For Capacity we use ṁ and |Δh| across the splice. + let (m_in, p_in, h_in) = self.inlet_state(state)?; + // When the configured measure is Enthalpy, return h_edge regardless of + // the `kind` the constraint asks for (we map Enthalpy → + // HeatTransferRate at the routing layer because MeasuredOutput has no + // Enthalpy variant). + if self.measure == ProbeMeasure::Enthalpy { + return Some(h_in); + } + + match kind { + MeasuredOutput::SaturationTemperature => Some(self.tsat_k(p_in).ok()?), + MeasuredOutput::Pressure => Some(p_in), + MeasuredOutput::MassFlowRate => Some(m_in.abs()), + MeasuredOutput::Temperature => { + // DGT/T → raw temperature at the probe. Falls back to None + // when the backend is unavailable or returns a non-physical T. + self.temperature_k(p_in, h_in).ok() + } + MeasuredOutput::Superheat | MeasuredOutput::Subcooling => { + let tsat = self.tsat_k(p_in).ok()?; + let t = self.temperature_k(p_in, h_in).ok()?; + if !t.is_finite() || t <= 0.0 { + return None; + } + match kind { + MeasuredOutput::Superheat => Some(t - tsat), + _ => Some(tsat - t), + } + } + MeasuredOutput::Capacity | MeasuredOutput::HeatTransferRate => { + // Capacity = ṁ · |Δh| across the probe splice. Uses both edges + // when available; if the outlet index isn't wired (defensive), + // falls back to zero duty. + if let Some((_, h_out)) = self.outlet_state(state) { + let dh = (h_in - h_out).abs(); + let q = m_in.abs() * dh; + if q.is_finite() { + return Some(q); + } + } + Some(0.0) + } + } + } + + fn signature(&self) -> String { + format!( + "Probe(measure={}, fluid={}, circuit={})", + self.measure.as_str(), + self.fluid_id_str, + self.circuit_id.0 + ) + } + + fn to_params(&self) -> crate::ComponentParams { + crate::ComponentParams::new("Probe") + .with_param("measure", self.measure.as_str().to_string()) + .with_param("fluid", self.fluid_id_str.clone()) + .with_param("circuitId", self.circuit_id.0) + } +} + +impl StateManageable for Probe { + fn state(&self) -> OperationalState { + self.operational_state + } + + fn set_state(&mut self, state: OperationalState) -> Result<(), ComponentError> { + if self.operational_state.can_transition_to(state) { + self.operational_state = state; + Ok(()) + } else { + Err(ComponentError::InvalidStateTransition { + from: self.operational_state, + to: state, + reason: "Transition not allowed".to_string(), + }) + } + } + + fn can_transition_to(&self, target: OperationalState) -> bool { + self.operational_state.can_transition_to(target) + } + + fn circuit_id(&self) -> &CircuitId { + &self.circuit_id + } + + fn set_circuit_id(&mut self, circuit_id: CircuitId) { + self.circuit_id = circuit_id; + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::port::FluidId; + use entropyk_core::{Enthalpy, Pressure}; + use entropyk_fluids::TestBackend; + + fn build_probe(measure: ProbeMeasure) -> Probe { + let backend: Arc = Arc::new(TestBackend::new()); + let mk = || { + Port::new( + FluidId::new("R134a"), + Pressure::from_bar(3.0), + Enthalpy::from_joules_per_kg(400_000.0), + ) + }; + let probe_disc = Probe::new(measure, "R134a", mk(), mk()).with_fluid_backend(backend); + // The connect() call needs a second pair of ports (mirrors Pipe/Pump + // CLI idiom — both pairs carry identical nominal values). + probe_disc.connect(mk(), mk()).expect("probe connect") + } + + /// Helper: equip a Probe with live edge indices pointing at a synthetic + /// state vector `[ṁ_in, P_in, h_in, ṁ_out, P_out, h_out]`. + fn wire_state_indices(probe: &mut Probe) { + probe.inlet_m_idx = Some(0); + probe.inlet_p_idx = Some(1); + probe.inlet_h_idx = Some(2); + probe.outlet_m_idx = Some(3); + probe.outlet_p_idx = Some(4); + probe.outlet_h_idx = Some(5); + } + + #[test] + fn probe_imposes_continuity_only() { + // The Probe is a zero-resistance tap: it contributes 2 topology + // residuals (P_out = P_in, h_out = h_in) so the spliced edge stays + // DoF-balanced. It does NOT contribute a calibration residual — + // the calibration Constraint (measure − target) supplies that. + let probe = build_probe(ProbeMeasure::Sh); + assert_eq!(probe.n_equations(), 2); + let roles = probe.equation_roles(); + assert_eq!(roles.len(), 2); + + // Continuity residuals: inlet edge state (ṁ=0.1, P=5e5, h=300e3), + // outlet edge state (ṁ=0.1, P=5e5, h=300e3) → r0 = 0, r1 = 0. + let mut probe = build_probe(ProbeMeasure::Sh); + wire_state_indices(&mut probe); + let state = vec![0.1, 5.0e5, 300_000.0, 0.1, 5.0e5, 300_000.0]; + let mut residuals = vec![0.0; 2]; + probe + .compute_residuals(&state, &mut residuals) + .expect("residuals compute"); + assert!(residuals[0].abs() < 1e-9, "P continuity: {}", residuals[0]); + assert!(residuals[1].abs() < 1e-9, "h continuity: {}", residuals[1]); + + // Mismatched inlet/outlet → non-zero continuity residuals. + let state2 = vec![0.1, 5.0e5, 300_000.0, 0.1, 4.0e5, 290_000.0]; + probe + .compute_residuals(&state2, &mut residuals) + .expect("residuals compute"); + assert!((residuals[0] - (-1.0e5)).abs() < 1e-3); + assert!((residuals[1] - (-10_000.0)).abs() < 1e-6); + + // Jacobian: ∂r0/∂P_out=+1, ∂r0/∂P_in=−1, ∂r1/∂h_out=+1, ∂r1/∂h_in=−1. + let mut jac = JacobianBuilder::new(); + probe + .jacobian_entries(&state, &mut jac) + .expect("jacobian computes"); + let entries = jac.entries(); + assert!(entries.contains(&(0, 4, 1.0)), "missing ∂r0/∂P_out"); + assert!(entries.contains(&(0, 1, -1.0)), "missing ∂r0/∂P_in"); + assert!(entries.contains(&(1, 5, 1.0)), "missing ∂r1/∂h_out"); + assert!(entries.contains(&(1, 2, -1.0)), "missing ∂r1/∂h_in"); + } + + #[test] + fn probe_port_names() { + let probe = build_probe(ProbeMeasure::Sh); + assert_eq!(probe.port_names(), vec!["inlet", "outlet"]); + assert_eq!(probe.flow_paths(), vec![(0, 1)]); + } + + #[test] + fn probe_measure_parse_round_trip() { + for kind in [ + ProbeMeasure::Sst, + ProbeMeasure::Sdt, + ProbeMeasure::Dgt, + ProbeMeasure::Dsh, + ProbeMeasure::Sh, + ProbeMeasure::Sc, + ProbeMeasure::T, + ProbeMeasure::P, + ProbeMeasure::MassFlow, + ProbeMeasure::Capacity, + ProbeMeasure::Enthalpy, + ] { + let parsed = ProbeMeasure::parse(kind.as_str()).expect("parses"); + assert_eq!(parsed, kind, "round-trip for {:?}", kind); + } + } + + #[test] + fn probe_measure_parse_rejects_unknown() { + assert!(ProbeMeasure::parse("NOPE").is_err()); + assert!(ProbeMeasure::parse("").is_err()); + } + + #[test] + fn probe_measure_parse_accepts_aliases() { + assert_eq!( + ProbeMeasure::parse("superheat").unwrap(), + ProbeMeasure::Sh + ); + assert_eq!( + ProbeMeasure::parse("mass_flow_rate").unwrap(), + ProbeMeasure::MassFlow + ); + assert_eq!( + ProbeMeasure::parse("subcooling").unwrap(), + ProbeMeasure::Sc + ); + } + + #[test] + fn probe_returns_none_without_state_indices() { + // No `set_system_context` call → indices are None → measure_output + // must return None gracefully instead of panicking. + let probe = build_probe(ProbeMeasure::Sh); + let state = vec![0.0; 6]; + assert!(probe + .measure_output(MeasuredOutput::Superheat, &state) + .is_none()); + assert!(probe + .measure_output(MeasuredOutput::Pressure, &state) + .is_none()); + } + + /// Sanity-check the Tsat derivation against TestBackend's R134a table at a + /// known pressure. R134a at 3 bar saturated → Tsat ≈ 273.4 K on TestBackend + /// (the table is approximate; we just assert the value is in a sane band). + #[test] + fn probe_saturation_temperature_in_sane_band() { + let mut probe = build_probe(ProbeMeasure::Sst); + wire_state_indices(&mut probe); + // 3 bar, h whatever — Tsat only uses P. + let state = vec![0.05, 3.0e5, 400_000.0, 0.05, 3.0e5, 400_000.0]; + let tsat = probe + .measure_output(MeasuredOutput::SaturationTemperature, &state) + .expect("Tsat must resolve"); + assert!( + (250.0..=320.0).contains(&tsat), + "Tsat out of sane band: {tsat}" + ); + } + + #[test] + fn probe_pressure_and_mass_flow_passthrough() { + let mut probe = build_probe(ProbeMeasure::P); + wire_state_indices(&mut probe); + let state = vec![0.42, 5.0e5, 250_000.0, 0.42, 5.0e5, 250_000.0]; + let p = probe + .measure_output(MeasuredOutput::Pressure, &state) + .unwrap(); + assert!((p - 5.0e5).abs() < 1e-6); + probe.measure = ProbeMeasure::MassFlow; + let m = probe + .measure_output(MeasuredOutput::MassFlowRate, &state) + .unwrap(); + assert!((m - 0.42).abs() < 1e-9); + } + + #[test] + fn probe_enthalpy_kind_returns_h_edge() { + let mut probe = build_probe(ProbeMeasure::Enthalpy); + wire_state_indices(&mut probe); + let state = vec![0.1, 3.0e5, 412_300.0, 0.1, 3.0e5, 412_300.0]; + // Enthalpy maps to HeatTransferRate in to_measured_output, but the + // measure_output override detects the configured Enthalpy kind and + // returns h_edge regardless. + let h = probe + .measure_output(MeasuredOutput::HeatTransferRate, &state) + .unwrap(); + assert!((h - 412_300.0).abs() < 1e-6); + } + + #[test] + fn probe_capacity_is_mass_flow_times_delta_h() { + let mut probe = build_probe(ProbeMeasure::Capacity); + wire_state_indices(&mut probe); + // Inlet h = 420 kJ/kg, outlet h = 400 kJ/kg, ṁ = 0.5 kg/s → Q = 10 kW. + let state = vec![0.5, 3.0e5, 420_000.0, 0.5, 3.0e5, 400_000.0]; + let q = probe + .measure_output(MeasuredOutput::Capacity, &state) + .unwrap(); + assert!((q - 10_000.0).abs() < 1e-6, "expected 10 kW, got {q}"); + } + + #[test] + fn probe_superheat_subcooling_signs() { + let backend: Arc = Arc::new(TestBackend::new()); + let mk = || { + Port::new( + FluidId::new("R134a"), + Pressure::from_bar(3.0), + Enthalpy::from_joules_per_kg(400_000.0), + ) + }; + // Probe with SH measure: should be T_edge − Tsat(P_edge). + let mut sh_probe = Probe::new(ProbeMeasure::Sh, "R134a", mk(), mk()) + .with_fluid_backend(Arc::clone(&backend)) + .connect(mk(), mk()) + .unwrap(); + wire_state_indices(&mut sh_probe); + let tsat = sh_probe.tsat_k(3.0e5).expect("tsat resolves"); + // Pick h that gives a known T via TestBackend. We can't know T without + // querying, so just assert SH and SC are negatives of each other for + // the same (P, h). + let state = vec![0.1, 3.0e5, 430_000.0, 0.1, 3.0e5, 430_000.0]; + let sh = sh_probe + .measure_output(MeasuredOutput::Superheat, &state) + .expect("SH resolves"); + + let mut sc_probe = Probe::new(ProbeMeasure::Sc, "R134a", mk(), mk()) + .with_fluid_backend(Arc::clone(&backend)) + .connect(mk(), mk()) + .unwrap(); + wire_state_indices(&mut sc_probe); + let sc = sc_probe + .measure_output(MeasuredOutput::Subcooling, &state) + .expect("SC resolves"); + assert!( + (sh + sc).abs() < 1e-6, + "SH and SC must be opposite signs for same (P,h): SH={sh}, SC={sc}, tsat={tsat}" + ); + } + + #[test] + fn probe_to_measured_output_mapping_is_exhaustive() { + // Every variant maps to a non-panic MeasuredOutput. + for kind in [ + ProbeMeasure::Sst, + ProbeMeasure::Sdt, + ProbeMeasure::Dgt, + ProbeMeasure::Dsh, + ProbeMeasure::Sh, + ProbeMeasure::Sc, + ProbeMeasure::T, + ProbeMeasure::P, + ProbeMeasure::MassFlow, + ProbeMeasure::Capacity, + ProbeMeasure::Enthalpy, + ] { + let _ = kind.to_measured_output(); + } + } + + #[test] + fn probe_signature_and_params_include_measure_and_fluid() { + let probe = build_probe(ProbeMeasure::Sst); + let sig = probe.signature(); + assert!(sig.contains("SST"), "signature: {sig}"); + assert!(sig.contains("R134a")); + let params = probe.to_params(); + assert_eq!(params.component_type, "Probe"); + } +} diff --git a/crates/entropyk/src/lib.rs b/crates/entropyk/src/lib.rs index 8b6dbe8..08303de 100644 --- a/crates/entropyk/src/lib.rs +++ b/crates/entropyk/src/lib.rs @@ -115,10 +115,11 @@ pub use entropyk_components::{ IncompressibleSplitter, IsenthalpicExpansionValve, IsentropicCompressor, JacobianBuilder, LmtdModel, MchxCondenserCoil, MockExternalModel, Node, NodeMeasurements, NodePhase, OperationalState, PerformanceCurves, PhaseRegion, Pipe, PipeGeometry, Polynomial1D, - Polynomial2D, Pump, PumpCurves, RefrigerantSink, RefrigerantSource, RegistryError, - ResidualVector, ScrewEconomizerCompressor, ScrewPerformanceCurves, ShellAndTubeHx, - SstSdtCoefficients, StateHistory, StateManageable, StateTransitionError, SystemState, - ThermalLoad, ThreadSafeExternalModel, UaMode, ValveCharacteristics, ValveFlowModel, + Polynomial2D, Probe, ProbeMeasure, Pump, PumpCurves, RefrigerantSink, RefrigerantSource, + RegistryError, ResidualVector, ScrewEconomizerCompressor, ScrewPerformanceCurves, + ShellAndTubeHx, SstSdtCoefficients, StateHistory, StateManageable, StateTransitionError, + SystemState, ThermalLoad, ThreadSafeExternalModel, UaMode, ValveCharacteristics, + ValveFlowModel, }; pub use entropyk_components::{ReversingMode, ReversingValve}; @@ -147,11 +148,11 @@ pub use entropyk_solver::{ AddEdgeError, AntoineCoefficients, CircuitConvergence, CircuitId as SolverCircuitId, ComponentOutput, Constraint, ConstraintError, ConstraintId, ConvergedState, ConvergenceCriteria, ConvergenceReason, ConvergenceReport, ConvergenceStatus, CyclePerformance, - DomainViolation, FallbackConfig, FallbackSolver, FlowEdge, HomotopyConfig, InitializerConfig, - InitializerError, JacobianFreezingConfig, JacobianMatrix, LinearSolver, MacroComponent, - MacroComponentSnapshot, NalgebraLuSolver, NewtonConfig, PicardConfig, PortMapping, - SmartInitializer, SolveOutcome, Solver, SolverError, SolverStrategy, System, ThermalCoupling, - TimeoutConfig, TopologyError, + DomainViolation, FaerLuSolver, FallbackConfig, FallbackSolver, FlowEdge, HomotopyConfig, + InitializerConfig, InitializerError, JacobianFreezingConfig, JacobianMatrix, LinearSolver, + MacroComponent, MacroComponentSnapshot, NalgebraLuSolver, NewtonConfig, PicardConfig, + PortMapping, SmartInitializer, SolveOutcome, Solver, SolverError, SolverStrategy, System, + ThermalCoupling, TimeoutConfig, TopologyError, }; // ============================================================================= diff --git a/crates/entropyk/src/result.rs b/crates/entropyk/src/result.rs index 4ee11cf..7841584 100644 --- a/crates/entropyk/src/result.rs +++ b/crates/entropyk/src/result.rs @@ -526,6 +526,42 @@ pub fn extract_solved_variables(system: &System, state: &[f64]) -> Vec 0 { tracing::debug!( skipped, @@ -547,9 +583,7 @@ pub fn extract_solved_variables(system: &System, state: &[f64]) -> Vec) -> String { let stripped = match component_id { - Some(c) => id - .strip_prefix(&format!("{}__", c)) - .unwrap_or(id), + Some(c) => id.strip_prefix(&format!("{}__", c)).unwrap_or(id), None => id.rfind("__").map(|i| &id[i + 2..]).unwrap_or(id), }; match stripped { @@ -781,6 +815,9 @@ mod tests { // Component unknown → fall back to last "__"-separated segment. assert_eq!(derive_solved_variable_label("glob__z_dp", None), "z_dp"); // No separator and no component → id itself. - assert_eq!(derive_solved_variable_label("global_var", None), "global_var"); + assert_eq!( + derive_solved_variable_label("global_var", None), + "global_var" + ); } } diff --git a/crates/solver/Cargo.toml b/crates/solver/Cargo.toml index b8fdaab..da7451a 100644 --- a/crates/solver/Cargo.toml +++ b/crates/solver/Cargo.toml @@ -11,6 +11,7 @@ repository = "https://github.com/entropyk/entropyk" entropyk-components = { path = "../components" } entropyk-core = { path = "../core" } entropyk-solver-core = { path = "../solver-core" } +faer = "0.24" nalgebra = "0.33" petgraph = "0.6" thiserror = "1.0" @@ -29,6 +30,10 @@ criterion = "0.5" name = "lu_solve" harness = false +[[bench]] +name = "lu_backends" +harness = false + [[bench]] name = "residual_jacobian_assembly" harness = false diff --git a/crates/solver/benches/lu_backends.rs b/crates/solver/benches/lu_backends.rs new file mode 100644 index 0000000..5bab32d --- /dev/null +++ b/crates/solver/benches/lu_backends.rs @@ -0,0 +1,85 @@ +//! Story 1.5 comparative benchmark: nalgebra vs faer dense LU backends. +//! +//! Measures factorization (`set_linearisation`) and per-RHS solve +//! (`solve_in_place`) on the mock-cycle Jacobian (9×9) and a synthetic 50×50 +//! dense Jacobian — the sizes Entropyk Newton actually sees. + +use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion}; +use entropyk_solver::linear::NalgebraLuSolver; +use entropyk_solver::linear_faer::FaerLuSolver; +use entropyk_solver::LinearSolver; + +mod common; + +fn dense_entries(n: usize) -> Vec<(usize, usize, f64)> { + let mut entries = Vec::with_capacity(n * n); + for i in 0..n { + for j in 0..n { + let v = if i == j { + 4.0 + 0.01 * (i as f64) + } else { + 1.0 / (1.0 + (i as f64 - j as f64).abs()) + }; + entries.push((i, j, v)); + } + } + entries +} + +fn bench_pair(c: &mut Criterion, label: &str, entries: &[(usize, usize, f64)], n: usize) { + let mut group = c.benchmark_group(label); + let rhs: Vec = (0..n) + .map(|i| (i as f64 * 0.7 - 1.3).sin() * 100.0) + .collect(); + + let mut nal = NalgebraLuSolver::new(); + nal.set_problem(n, n).unwrap(); + let mut faer = FaerLuSolver::new(); + faer.set_problem(n, n).unwrap(); + + group.bench_function(BenchmarkId::new("factor", "nalgebra"), |b| { + b.iter(|| { + black_box(nal.set_linearisation(black_box(entries))).unwrap(); + }); + }); + group.bench_function(BenchmarkId::new("factor", "faer"), |b| { + b.iter(|| { + black_box(faer.set_linearisation(black_box(entries))).unwrap(); + }); + }); + + nal.set_linearisation(entries).unwrap(); + faer.set_linearisation(entries).unwrap(); + group.bench_function(BenchmarkId::new("solve", "nalgebra"), |b| { + b.iter(|| { + let mut rhs = rhs.clone(); + let _ = black_box(nal.solve_in_place(black_box(&mut rhs))); + }); + }); + group.bench_function(BenchmarkId::new("solve", "faer"), |b| { + b.iter(|| { + let mut rhs = rhs.clone(); + let _ = black_box(faer.solve_in_place(black_box(&mut rhs))); + }); + }); + group.finish(); +} + +fn bench_lu_backends(c: &mut Criterion) { + let (system, state) = common::build_mock_system(); + let jacobian = common::assemble_jacobian(&system, &state); + let n = jacobian.nrows(); + let jm = jacobian.as_matrix(); + let entries: Vec<(usize, usize, f64)> = (0..n) + .flat_map(|i| (0..n).map(move |j| (i, j, jm[(i, j)]))) + .collect(); + bench_pair(c, "lu_mock_9x9", &entries, n); + + let entries50 = dense_entries(50); + bench_pair(c, "lu_dense_50x50", &entries50, 50); + + let _ = state; +} + +criterion_group!(benches, bench_lu_backends); +criterion_main!(benches); diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index a4fd6c1..4caba80 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -34,6 +34,7 @@ pub mod initializer; pub mod inverse; pub mod jacobian; pub mod linear; +pub mod linear_faer; pub mod macro_component; pub mod metadata; pub mod scaling; @@ -64,7 +65,8 @@ pub use initializer::{ }; pub use inverse::{ComponentOutput, Constraint, ConstraintError, ConstraintId}; pub use jacobian::JacobianMatrix; -pub use linear::NalgebraLuSolver; +pub use linear::{make_linear_backend, set_linear_backend_override, DenseLu, NalgebraLuSolver}; +pub use linear_faer::FaerLuSolver; pub use macro_component::{MacroComponent, MacroComponentSnapshot, PortMapping}; pub use metadata::SimulationMetadata; pub use scaling::{equilibrate, unscale_dx}; diff --git a/crates/solver/src/linear.rs b/crates/solver/src/linear.rs index 621f5fb..6023a07 100644 --- a/crates/solver/src/linear.rs +++ b/crates/solver/src/linear.rs @@ -356,3 +356,139 @@ mod tests { )); } } + +// ─── Strangler dispatch (Story 1.5) ───────────────────────────────────────── + +/// Process-wide override for the default dense backend (set by the CLI from +/// `solver.linear_backend`). Precedence: explicit `make_linear_backend` name +/// → this override → `ENTROPYK_LINEAR_BACKEND` env → `"nalgebra"`. +static BACKEND_OVERRIDE: std::sync::RwLock> = std::sync::RwLock::new(None); + +/// Sets the process-wide default dense backend (`Some("faer")` / `Some("nalgebra")`, +/// `None` to clear). Additive, additive-only configuration surface. +pub fn set_linear_backend_override(name: Option<&str>) { + if let Ok(mut guard) = BACKEND_OVERRIDE.write() { + *guard = name.map(str::to_string); + } +} + +/// The dense LU backend in use: nalgebra (default, legacy path) or faer +/// (opt-in). Strategies hold this enum directly so the `set_matrix` fast path +/// stays available; it also implements [`LinearSolver`] for registry boxing. +pub enum DenseLu { + /// nalgebra dense LU (default). + Nalgebra(NalgebraLuSolver), + /// faer dense LU (Story 1.5 opt-in). + Faer(crate::linear_faer::FaerLuSolver), +} + +impl DenseLu { + /// Inherent fast path shared by both backends (see each backend's docs). + pub fn set_matrix(&mut self, matrix: &DMatrix) -> Result<(), SolverCoreError> { + match self { + Self::Nalgebra(b) => b.set_matrix(matrix), + Self::Faer(b) => b.set_matrix(matrix), + } + } +} + +impl LinearSolver for DenseLu { + fn set_problem( + &mut self, + n_equations: usize, + n_unknowns: usize, + ) -> Result<(), SolverCoreError> { + match self { + Self::Nalgebra(b) => b.set_problem(n_equations, n_unknowns), + Self::Faer(b) => b.set_problem(n_equations, n_unknowns), + } + } + + fn set_linearisation( + &mut self, + entries: &[(usize, usize, f64)], + ) -> Result<(), SolverCoreError> { + match self { + Self::Nalgebra(b) => b.set_linearisation(entries), + Self::Faer(b) => b.set_linearisation(entries), + } + } + + fn solve_in_place(&mut self, rhs: &mut [f64]) -> Result<(), SolverCoreError> { + match self { + Self::Nalgebra(b) => b.solve_in_place(rhs), + Self::Faer(b) => b.solve_in_place(rhs), + } + } + + fn n(&self) -> usize { + match self { + Self::Nalgebra(b) => b.n(), + Self::Faer(b) => b.n(), + } + } +} + +/// Builds the dense backend selected by name (or by override/env/default). +pub fn make_linear_backend(name: Option<&str>) -> Result { + let selected = name + .map(str::to_string) + .or_else(|| BACKEND_OVERRIDE.read().ok().and_then(|g| g.clone())) + .or_else(|| std::env::var("ENTROPYK_LINEAR_BACKEND").ok()) + .unwrap_or_else(|| "nalgebra".to_string()); + match selected.as_str() { + "nalgebra" => Ok(DenseLu::Nalgebra(NalgebraLuSolver::new())), + "faer" => Ok(DenseLu::Faer(crate::linear_faer::FaerLuSolver::new())), + other => Err(SolverCoreError::Usage { + message: format!("unknown linear backend '{other}' (use 'nalgebra' or 'faer')"), + }), + } +} + +#[cfg(test)] +mod dispatch_tests { + use super::*; + + #[test] + fn factory_builds_named_backends_and_rejects_unknown() { + assert!(matches!( + make_linear_backend(Some("nalgebra")).unwrap(), + DenseLu::Nalgebra(_) + )); + assert!(matches!( + make_linear_backend(Some("faer")).unwrap(), + DenseLu::Faer(_) + )); + assert!(make_linear_backend(Some("mumps")).is_err()); + // Env-independent check: an explicit override always wins. + set_linear_backend_override(Some("nalgebra")); + assert!(matches!( + make_linear_backend(None).unwrap(), + DenseLu::Nalgebra(_) + )); + set_linear_backend_override(None); + } + + #[test] + fn override_selects_default_and_explicit_name_wins() { + set_linear_backend_override(Some("faer")); + assert!(matches!( + make_linear_backend(None).unwrap(), + DenseLu::Faer(_) + )); + // Explicit name wins over the override. + assert!(matches!( + make_linear_backend(Some("nalgebra")).unwrap(), + DenseLu::Nalgebra(_) + )); + set_linear_backend_override(None); + // After clearing, selection falls back to env-or-nalgebra: just check + // it builds a working backend. + let mut b = make_linear_backend(None).unwrap(); + b.set_problem(1, 1).unwrap(); + b.set_linearisation(&[(0, 0, 2.0)]).unwrap(); + let mut rhs = [4.0]; + b.solve_in_place(&mut rhs).unwrap(); + assert!((rhs[0] - 2.0).abs() < 1e-9); + } +} diff --git a/crates/solver/src/linear_faer.rs b/crates/solver/src/linear_faer.rs new file mode 100644 index 0000000..6a0d76e --- /dev/null +++ b/crates/solver/src/linear_faer.rs @@ -0,0 +1,310 @@ +//! `FaerLuSolver`: dense faer LU backend behind the [`LinearSolver`] lifecycle +//! (FR3, Story 1.5 — strangler migration; nalgebra stays the default). +//! +//! Same contract and pipeline as [`crate::linear::NalgebraLuSolver`] (post-1.4 +//! audit): identical triplet assembly (zero-fill + `+=`), identical Ruiz +//! equilibration (`crate::scaling::equilibrate`), factorization computed once +//! in `set_linearisation`, bounds validated BEFORE mutation with the factor +//! invalidated on out-of-bounds, zero-panic error paths, zero heap allocation +//! in `solve_in_place` (faer solves in place on the stored factorization, +//! over a column-major view of the pre-allocated scratch). The only +//! intentional difference is the factorization engine: faer's partial-pivot +//! LU instead of nalgebra's LU — results match within tolerance, not bitwise. + +use entropyk_solver_core::{LinearSolver, SolverCoreError}; +use faer::linalg::solvers::{PartialPivLu, SolveCore}; +use nalgebra::DMatrix; + +/// Stored faer factorization plus the equilibration scalings used to build it. +struct Factorized { + d_r: Vec, + d_c: Vec, + lu: PartialPivLu, +} + +/// Dense LU backend (faer) behind the object-safe [`LinearSolver`] trait. +/// +/// Scratch buffers are allocated once in `set_problem` so `solve_in_place` +/// performs zero heap allocation. +pub struct FaerLuSolver { + n_rows: usize, + n_cols: usize, + /// Assembled (unscaled) matrix (nalgebra storage, shared with the + /// equilibration pipeline). + matrix: Option>, + factor: Option, + /// Scratch: scaled rhs / unscaled step (allocated in `set_problem`). + delta: Vec, +} + +impl Default for FaerLuSolver { + fn default() -> Self { + Self::new() + } +} + +impl FaerLuSolver { + /// Creates an unconfigured backend (call [`LinearSolver::set_problem`]). + pub fn new() -> Self { + Self { + n_rows: 0, + n_cols: 0, + matrix: None, + factor: None, + delta: Vec::new(), + } + } + + /// Inherent fast path: copy an already-assembled dense matrix and factor + /// it (one `copy_from`, no triplet conversion). + pub fn set_matrix(&mut self, matrix: &DMatrix) -> Result<(), SolverCoreError> { + let stored = self.matrix.as_mut().ok_or_else(|| SolverCoreError::Usage { + message: "set_problem must be called before set_matrix".to_string(), + })?; + if stored.nrows() != matrix.nrows() || stored.ncols() != matrix.ncols() { + return Err(SolverCoreError::Usage { + message: format!( + "matrix shape {}x{} does not match problem {}x{}", + matrix.nrows(), + matrix.ncols(), + stored.nrows(), + stored.ncols() + ), + }); + } + stored.copy_from(matrix); + self.factorize() + } + + /// Equilibrate + scale + factorize the stored matrix (square path). + fn factorize(&mut self) -> Result<(), SolverCoreError> { + let Some(matrix) = self.matrix.as_ref() else { + self.factor = None; + return Err(SolverCoreError::Usage { + message: "set_problem must be called before factorizing".to_string(), + }); + }; + if matrix.nrows() != matrix.ncols() { + // Non-square problems have no LU path; solve_in_place reports Usage. + self.factor = None; + return Ok(()); + } + let n = matrix.nrows(); + let (d_r, d_c) = crate::scaling::equilibrate(matrix); + // Scaled matrix as a faer Mat (column-major, one copy at factor time). + let scaled = faer::Mat::from_fn(n, n, |i, j| matrix[(i, j)] * d_r[i] * d_c[j]); + self.factor = Some(Factorized { + d_r, + d_c, + lu: PartialPivLu::new(scaled.as_ref()), + }); + Ok(()) + } +} + +impl LinearSolver for FaerLuSolver { + fn set_problem( + &mut self, + n_equations: usize, + n_unknowns: usize, + ) -> Result<(), SolverCoreError> { + if n_equations == 0 || n_unknowns == 0 { + return Err(SolverCoreError::Usage { + message: format!("zero-sized problem {n_equations}x{n_unknowns}"), + }); + } + self.n_rows = n_equations; + self.n_cols = n_unknowns; + self.matrix = Some(DMatrix::zeros(n_equations, n_unknowns)); + self.factor = None; + // Scratch buffer for the zero-allocation solve path. + self.delta = vec![0.0; n_unknowns]; + Ok(()) + } + + fn set_linearisation( + &mut self, + entries: &[(usize, usize, f64)], + ) -> Result<(), SolverCoreError> { + let (n_rows, n_cols) = (self.n_rows, self.n_cols); + if self.matrix.is_none() { + return Err(SolverCoreError::Usage { + message: "set_problem must be called before set_linearisation".to_string(), + }); + } + // Validate bounds BEFORE touching the matrix (1.4 audit contract). + for &(row, col, _) in entries { + if row >= n_rows || col >= n_cols { + self.factor = None; + return Err(SolverCoreError::Usage { + message: format!("entry ({row},{col}) out of bounds {n_rows}x{n_cols}"), + }); + } + } + let matrix = self.matrix.as_mut().ok_or_else(|| SolverCoreError::Usage { + message: "set_problem must be called before set_linearisation".to_string(), + })?; + matrix.fill(0.0); + for &(row, col, value) in entries { + matrix[(row, col)] += value; + } + self.factorize() + } + + fn solve_in_place(&mut self, rhs: &mut [f64]) -> Result<(), SolverCoreError> { + if rhs.len() != self.n_rows || self.n_rows != self.n_cols { + return Err(SolverCoreError::Usage { + message: format!( + "solve_in_place needs a square rhs of len {}, got {} (problem {}x{})", + self.n_rows, + rhs.len(), + self.n_rows, + self.n_cols + ), + }); + } + let n = self.n_rows; + // Disjoint field borrows: factorization + scratch (zero alloc). + let Self { factor, delta, .. } = self; + let factor = factor.as_ref().ok_or_else(|| SolverCoreError::Usage { + message: "set_linearisation must be called before solve_in_place".to_string(), + })?; + // Scaled right-hand side: D_r · b (into the scratch). + for i in 0..n { + delta[i] = rhs[i] * factor.d_r[i]; + } + // faer solves in place on a column-major view of the scratch. + let rhs_view = faer::MatMut::from_column_major_slice_mut(&mut delta[..n], n, 1); + factor.lu.solve_in_place_with_conj(faer::Conj::No, rhs_view); + // Undo the column scaling in place: x_i = D_c,i · y_i. + for i in 0..n { + delta[i] *= factor.d_c[i]; + } + if delta[..n].iter().all(|v| v.is_finite()) { + rhs.copy_from_slice(&delta[..n]); + Ok(()) + } else { + tracing::warn!( + "faer LU solve produced a non-finite step - Jacobian may contain NaN/Inf" + ); + Err(SolverCoreError::InvalidSystem { + message: "linear solve produced a non-finite step".to_string(), + }) + } + } + + fn n(&self) -> usize { + self.n_rows + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::jacobian::JacobianMatrix; + + /// Deterministic dense system (well-conditioned, diagonally dominant). + fn dense_entries(n: usize) -> Vec<(usize, usize, f64)> { + let mut entries = Vec::new(); + for i in 0..n { + for j in 0..n { + let v = if i == j { + 4.0 + 0.1 * (i as f64) + } else { + 1.0 / (1.0 + (i as f64 - j as f64).abs()) + }; + entries.push((i, j, v)); + } + } + entries + } + + #[test] + fn parity_with_jacobian_matrix_solve() { + for n in [1, 2, 5, 12] { + let entries = dense_entries(n); + let jm = JacobianMatrix::from_builder(&entries, n, n); + let residuals: Vec = (0..n) + .map(|i| (i as f64 * 0.7 - 1.3).sin() * 100.0) + .collect(); + + let legacy = jm.solve(&residuals).expect("legacy solve"); + + let mut backend = FaerLuSolver::new(); + backend.set_problem(n, n).unwrap(); + backend.set_linearisation(&entries).unwrap(); + let mut rhs: Vec = residuals.iter().map(|r| -*r).collect(); + backend.solve_in_place(&mut rhs).unwrap(); + + for i in 0..n { + // faer ≠ nalgebra bitwise: tolerance, not exact equality. + assert!( + (legacy[i] - rhs[i]).abs() <= 1e-9 * legacy[i].abs().max(1.0), + "n={n} i={i}: legacy={} faer={}", + legacy[i], + rhs[i] + ); + } + } + } + + #[test] + fn factorization_reuse_across_rhs() { + let n = 6; + let entries = dense_entries(n); + let jm = JacobianMatrix::from_builder(&entries, n, n); + let mut backend = FaerLuSolver::new(); + backend.set_problem(n, n).unwrap(); + backend.set_matrix(jm.as_matrix()).unwrap(); + + for k in 0..2 { + let residuals: Vec = (0..n).map(|i| (i + k) as f64 * 3.1 - 5.0).collect(); + let legacy = jm.solve(&residuals).expect("legacy solve"); + let mut rhs: Vec = residuals.iter().map(|r| -*r).collect(); + backend.solve_in_place(&mut rhs).unwrap(); + for i in 0..n { + assert!( + (legacy[i] - rhs[i]).abs() <= 1e-9 * legacy[i].abs().max(1.0), + "reuse k={k} i={i}" + ); + } + } + } + + #[test] + fn error_paths_and_stale_factor_invalidation() { + let mut backend = FaerLuSolver::new(); + assert!(backend.set_problem(0, 2).is_err()); + backend.set_problem(2, 2).unwrap(); + // Solve before set_linearisation → Usage. + assert!(matches!( + backend.solve_in_place(&mut [1.0, 2.0]), + Err(SolverCoreError::Usage { .. }) + )); + // Factor a good matrix, then an OOB entry: the stale factorization + // must be invalidated (1.4 audit regression). + backend + .set_linearisation(&[(0, 0, 2.0), (1, 1, 1.0)]) + .unwrap(); + assert!(backend.set_linearisation(&[(2, 0, 1.0)]).is_err()); + assert!(matches!( + backend.solve_in_place(&mut [1.0, 2.0]), + Err(SolverCoreError::Usage { .. }) + )); + assert_eq!(backend.n(), 2); + } + + #[test] + fn non_finite_matrix_entries_are_rejected() { + let mut backend = FaerLuSolver::new(); + backend.set_problem(2, 2).unwrap(); + backend + .set_linearisation(&[(0, 0, f64::NAN), (1, 1, 1.0)]) + .unwrap(); + let mut rhs = vec![1.0, 2.0]; + assert!(matches!( + backend.solve_in_place(&mut rhs), + Err(SolverCoreError::InvalidSystem { .. }) + )); + } +} diff --git a/crates/solver/src/strategies/newton_raphson.rs b/crates/solver/src/strategies/newton_raphson.rs index c7744cb..ae08b85 100644 --- a/crates/solver/src/strategies/newton_raphson.rs +++ b/crates/solver/src/strategies/newton_raphson.rs @@ -345,10 +345,13 @@ impl Solver for NewtonConfig { let mut frozen_count: usize = 0; let mut force_recompute: bool = true; - // LinearSolver backend (FR3): factorizes at each fresh assembly and - // reuses the factorization (with cached scalings) across frozen - // iterations and repeated solves — bit-identical to JacobianMatrix::solve. - let mut linear_backend = crate::linear::NalgebraLuSolver::new(); + // LinearSolver backend (FR3/Story 1.5): nalgebra by default, faer via + // config/env dispatch. Factorizes at each fresh assembly and reuses + // the factorization (with cached scalings) across frozen iterations. + let mut linear_backend = + crate::linear::make_linear_backend(None).map_err(|e| SolverError::InvalidSystem { + message: format!("Failed to select linear backend: {e}"), + })?; linear_backend .set_problem(n_equations, n_state) .map_err(|e| SolverError::InvalidSystem { diff --git a/crates/solver/src/strategies/pseudo_transient.rs b/crates/solver/src/strategies/pseudo_transient.rs index 2421031..bdd7dff 100644 --- a/crates/solver/src/strategies/pseudo_transient.rs +++ b/crates/solver/src/strategies/pseudo_transient.rs @@ -122,9 +122,12 @@ impl Solver for PtcConfig { let mut jacobian_builder = JacobianBuilder::new(); let mut jacobian_matrix = JacobianMatrix::zeros(n_equations, n_state); - // LinearSolver backend (FR3): factorizes at each fresh assembly; - // bit-identical to JacobianMatrix::solve. - let mut linear_backend = crate::linear::NalgebraLuSolver::new(); + // LinearSolver backend (FR3/Story 1.5): nalgebra by default, faer via + // config/env dispatch; factorizes at each fresh assembly. + let mut linear_backend = + crate::linear::make_linear_backend(None).map_err(|e| SolverError::InvalidSystem { + message: format!("Failed to select linear backend: {e}"), + })?; linear_backend .set_problem(n_equations, n_state) .map_err(|e| SolverError::InvalidSystem { diff --git a/crates/solver/src/strategies/sequential_substitution.rs b/crates/solver/src/strategies/sequential_substitution.rs index 05f57a7..6026426 100644 --- a/crates/solver/src/strategies/sequential_substitution.rs +++ b/crates/solver/src/strategies/sequential_substitution.rs @@ -335,6 +335,26 @@ impl Solver for PicardConfig { let mut divergence_count: usize = 0; let mut previous_norm: f64; + // Solver bounds mask (pressure floor + bounded control variables): + // Picard's relaxed update has no line search, so without clipping it + // can drive a calibration z-factor far outside its bounds (observed: + // z_ua → −2.8 on an SDT calibration), then crash the components. + // Saturated-controller actuator slots are EXCLUDED: they carry their + // own saturation semantics (anti-windup integrator + S(x) inside the + // controller residual) and hard-clamping them here breaks that dynamic. + let mut clipping_mask: Vec> = (0..n_state) + .map(|i| system.get_solver_bounds_for_state_index(i)) + .collect(); + { + let n_sat = system.saturated_controllers_mut().count(); + for i in 0..n_sat { + let u_idx = system.saturated_u_index(i); + if u_idx < clipping_mask.len() { + clipping_mask[u_idx] = None; + } + } + } + // Pre-allocate best-state tracking buffer (Story 4.5 - AC: #5) let mut best_state: Vec = vec![0.0; n_state]; let mut best_residual: f64; @@ -425,6 +445,16 @@ impl Solver for PicardConfig { } else { Self::apply_relaxation(&mut state, &residuals, self.relaxation_factor); } + // Clamp the relaxed/extrapolated iterate into solver bounds (see + // the mask comment above): bounded control variables (calibration + // z-factors, openings) and pressure floors. + for (i, s) in state.iter_mut().enumerate() { + if let Some((min, max)) = &clipping_mask[i] { + if min.is_finite() && max.is_finite() && min <= max { + *s = s.clamp(*min, *max); + } + } + } // Compute new residuals. Recoverable domain violation → typed // outcome; fatal → InvalidSystem flattening. diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 419980b..fc7e0a1 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -224,6 +224,11 @@ pub struct System { /// Registry of component names for constraint validation. /// Maps human-readable names (e.g., "evaporator") to NodeIndex. component_names: HashMap, + /// Per-component calibration indices snapshot, captured at `finalize()` time. + /// Lets external readers (e.g. result extraction) learn which Z-factors + /// (z_ua, z_dp, z_flow, …) were promoted to free unknowns and where they + /// live in the state vector, so the solved values can be surfaced to the user. + calib_indices_by_name: HashMap, finalized: bool, total_state_len: usize, /// When `true` (default), `finalize` rejects **over-constrained** systems. @@ -246,6 +251,7 @@ impl System { saturated_controllers: Vec::new(), free_actuators: Vec::new(), component_names: HashMap::new(), + calib_indices_by_name: HashMap::new(), finalized: false, total_state_len: 0, enforce_dof_gate: true, @@ -682,6 +688,10 @@ impl System { } } + // Persist the per-component calib map so external readers (result + // extraction) can surface solved Z-factor values (z_ua, z_dp, …). + self.calib_indices_by_name = comp_calib_indices; + // Wire physical thermal couplings (Story 3.4 completion): each coupling // owns one state unknown Q [W] at `coupling_state_index(i)`. The // cold-side receiver component (e.g. `ThermalLoad`) reads Q in its @@ -1145,6 +1155,15 @@ impl System { self.component_names.keys().map(|s| s.as_str()) } + /// Per-component `CalibIndices` snapshot captured at `finalize()`. + /// Each `Some(idx)` slot (z_ua, z_dp, z_flow, z_flow_eco, z_power, z_etav, + /// actuator) means that factor was promoted to a free solver unknown and + /// its solved value lives at `state[idx]`. Used by result extraction to + /// surface solved calibration/control values to the user. + pub fn calib_indices_by_name(&self) -> &HashMap { + &self.calib_indices_by_name + } + /// Returns a reference to the component stored at the given node index. /// /// # Panics @@ -2514,11 +2533,11 @@ impl System { self.total_state_len + self.inverse_control.mapping_count() + i } - fn saturated_base_index(&self) -> usize { + pub fn saturated_base_index(&self) -> usize { self.total_state_len + self.inverse_control.mapping_count() + self.coupling_residual_count() } - fn saturated_u_index(&self, i: usize) -> usize { + pub fn saturated_u_index(&self, i: usize) -> usize { self.saturated_base_index() + 2 * i } diff --git a/crates/solver/tests/faer_backend_parity.rs b/crates/solver/tests/faer_backend_parity.rs new file mode 100644 index 0000000..4ed0d90 --- /dev/null +++ b/crates/solver/tests/faer_backend_parity.rs @@ -0,0 +1,83 @@ +//! Story 1.5 AC#2/#3: the faer dense backend solves the reference cycles with +//! results matching the nalgebra baseline within tolerance (not bitwise — +//! different factorization engines), via the strangler dispatch. + +mod common; + +use entropyk_solver::linear::{make_linear_backend, set_linear_backend_override, DenseLu}; +use entropyk_solver::{LinearSolver, Solver}; + +/// Per-value tolerance matching the golden-snapshot comparison policy. +const STATE_RTOL: f64 = 1e-6; +const STATE_ATOL: f64 = 1.0; // Pa / J/kg absolute floor + +fn solve_all_cycles() -> Vec> { + [ + common::build_reference_cycle_a(), + common::build_reference_cycle_b(), + common::build_reference_cycle_c(), + ] + .into_iter() + .map(|mut system| { + common::solve_reference_system_with_state(&mut system) + .state + .clone() + }) + .collect() +} + +#[test] +fn faer_path_dispatch_and_reference_cycle_parity() { + // 1. Unknown override must fail the solve with the dispatch error (the + // dispatch is on the strategy path, not dead code). + set_linear_backend_override(Some("definitely-not-a-backend")); + let mut system = common::build_reference_cycle_a(); + let result = entropyk_solver::FallbackSolver::default_solver().solve(&mut system); + assert!(result.is_err(), "unknown backend must fail the solve"); + + // 2. Baseline (nalgebra) on the three reference cycles. + set_linear_backend_override(Some("nalgebra")); + let baseline = solve_all_cycles(); + + // 3. Same cycles on faer; converged states must match within tolerance. + set_linear_backend_override(Some("faer")); + let faer = solve_all_cycles(); + + set_linear_backend_override(None); + + assert_eq!(baseline.len(), faer.len()); + for (cycle_idx, (base, alt)) in baseline.iter().zip(faer.iter()).enumerate() { + assert_eq!( + base.len(), + alt.len(), + "cycle {cycle_idx}: state length differs" + ); + for (i, (b, a)) in base.iter().zip(alt.iter()).enumerate() { + let scale = b.abs().max(a.abs()).max(STATE_ATOL); + assert!( + (b - a).abs() <= STATE_RTOL * scale, + "cycle {cycle_idx} state[{i}]: nalgebra={b} faer={a} (rel diff {})", + (b - a).abs() / scale + ); + } + } +} + +#[test] +fn dense_lu_enum_dispatches_both_engines() { + for name in ["nalgebra", "faer"] { + let mut backend = make_linear_backend(Some(name)).unwrap(); + backend.set_problem(2, 2).unwrap(); + backend + .set_linearisation(&[(0, 0, 2.0), (1, 1, 1.0)]) + .unwrap(); + let mut rhs = vec![4.0, 3.0]; + backend.solve_in_place(&mut rhs).unwrap(); + assert!((rhs[0] - 2.0).abs() < 1e-9, "{name}: {}", rhs[0]); + assert!((rhs[1] - 3.0).abs() < 1e-9, "{name}: {}", rhs[1]); + } + assert!(matches!( + make_linear_backend(Some("faer")).unwrap(), + DenseLu::Faer(_) + )); +} diff --git a/docs/benchmarks/faer-vs-nalgebra.md b/docs/benchmarks/faer-vs-nalgebra.md new file mode 100644 index 0000000..8d885c7 --- /dev/null +++ b/docs/benchmarks/faer-vs-nalgebra.md @@ -0,0 +1,21 @@ +# faer vs nalgebra dense LU — comparative report (Story 1.5, 2026-07-19) + +Environment: Windows 11 x64, rustc 1.89 (stable, pinned), release profile, criterion 0.5. +Backends: `NalgebraLuSolver` (nalgebra 0.33 LU) vs `FaerLuSolver` (faer 0.24.4 partial-pivot LU), both behind the `LinearSolver` trait with identical Ruiz equilibration and zero-alloc solve. + +| Benchmark | nalgebra | faer | ratio (faer/nalgebra) | +|---|---|---|---| +| factor 9×9 (mock cycle Jacobian) | 308 ns | 1.23 µs | **4.0× slower** | +| solve-per-RHS 9×9 | 83 ns | 245 ns | 3.0× slower | +| factor 50×50 (synthetic dense) | 19.4 µs | 20.4 µs | 1.05× slower | +| solve-per-RHS 50×50 | 589 ns | 1.11 µs | 1.9× slower | + +## Verdict + +At Entropyk's current system sizes (n ≤ ~100), **faer is NOT faster than nalgebra** — the assembled-matrix conversion (nalgebra → faer column-major Mat) costs at factor time on small n, and faer's threading edge only materializes at larger n. faer stays an **opt-in** alternative (`solver.linear_backend = "faer"` / `ENTROPYK_LINEAR_BACKEND=faer`); **nalgebra remains the default** — the strangler data does not justify retirement, and none was planned in this story. + +Revisit when: (a) block-dense/Schur systems push n into the hundreds (Epic 2 / Story 1.8), (b) faer-sparse becomes relevant, (c) batch/parallel solves can exploit faer's rayon threading (Epic 4). + +Parity (AC#3): reference cycles A/B/C converge identically within 1e-6 relative on both backends; full solver corpus green on both paths (38/38 test binaries, 0 FAILED). + +Reproduce: `cargo bench -p entropyk-solver --bench lu_backends`