engine: heat in a fluid step

A fluid domain takes `thermal` settings (the initial temperature, a
thermal condition per box face and per obstacle face, gravity with the
expansion coefficient for buoyancy) in model.fluid_domain.create and
.edit, and model.fluid_domain.thermal_boundary sets one boundary's
condition as model.fluid_domain.boundary does for the flow. A material
a domain with heat uses keeps needing its conductivity and specific
heat when it is replaced.

solve.fluid's summary adds the thermal relaxation time, the fluid's
temperature range and the heat into every wall (final and the mean
over the last quarter of the history); solve.probe gives the
temperature in the cells at points and, with `heat`, the walls' heat
histories. The form, the catalogue and docs/features.md follow.

Part of #157.

crates/engine/src/catalogue.rs

@@ -543,12 +543,13 @@ pub static CATALOGUE: &[Capability] = &[
         "The box of fluid cells, its fluid, and the obstacle bodies voxelised from their face-tagged surfaces",
         &["model.fluid_domain.edit", "model.fluid_domain.delete"],
     ),
-    ready(
+    ready_more(
         "model.fluid_domain.boundary",
         "Inlet & Outlet",
         Loads,
         "Flow",
-        "Velocity inlets, pressure outlets, walls, moving and slip walls, periodic faces",
+        "Velocity inlets, pressure outlets, walls, moving and slip walls, periodic faces; held temperatures, heat fluxes and adiabatic faces",
+        &["model.fluid_domain.thermal_boundary"],
     ),
     // ---- Interactions ----
     ready_more(
@@ -649,7 +650,13 @@ pub static CATALOGUE: &[Capability] = &[
         "Transient dynamics of an explicit step by central differences (contact and erosion coming in M10)",
     ),
     soon("solve.blast", "Blast & Impact", Solve, "Explicit", "M14", "EOS-driven fluids with fluid–structure contact (ALE or SPH)"),
-    ready("solve.fluid", "Lattice Boltzmann", Solve, "Fluids", "Lattice Boltzmann flow of a fluid step on the CPU, with wall forces"),
+    ready(
+        "solve.fluid",
+        "Lattice Boltzmann",
+        Solve,
+        "Fluids",
+        "Lattice Boltzmann flow of a fluid step on the CPU or a GPU, with wall forces, and the heat it carries",
+    ),
     soon("solve.fsi", "Fluid–Structure", Solve, "Fluids", "M13", "Coupled fluid and structure solve"),
     // ---- Results ----
     ready_more("results.deformed", "Deformed", Results, "Display", "The deformed shape of a solved step, to scale", &["results.clear"]),

crates/engine/src/cmd/define.rs

@@ -8,7 +8,7 @@
 
 use fieldcraft_model::{
     BoundaryCondition, Coupling, FaceCondition, FluidBoundary, FluidDomain, Load, Material, Model, PointMass, ReferencePoint, Section, Set, Spring,
-    Step,
+    Step, ThermalCondition,
 };
 use serde::de::DeserializeOwned;
 use serde_json::{Value, json};
@@ -288,7 +288,7 @@ definition_commands!(
     FluidDomain,
     "fluid_domain",
     "Fluid Domain",
-    r#"{"name": "Flow", "material": "Water" (a density and fluid behaviour), "min": [x, y, z], "max": [x, y, z], "spacing": 1.0 (cubic cells; every extent a whole number of spacings), "faces": {"x_min": {"type": "velocity", "velocity": [u, v, w], "profile": "uniform|parabolic"}, "x_max": {"type": "pressure", "pressure": 0}, "y_min": {"type": "wall"} | {"type": "moving_wall", "velocity": [u, v, w] (in the face's plane)} | {"type": "slip"} | {"type": "periodic"} (in opposite pairs), …} (default walls), "obstacles": ["Body1"] (solid bodies of the design or imported surfaces, by name), "surfaces": [{"face": {"body": "Body1", "face": "F2:top"}, "condition": {"type": "moving_wall", "velocity": [u, v, w]}}] (obstacle faces are walls unless listed), "initial_velocity": [u, v, w], "acceleration": [x, y, z] (gravity or a driving force per mass), "tolerance": 0.1 (obstacle triangulation; default a tenth of the spacing)}"#,
+    r#"{"name": "Flow", "material": "Water" (a density and fluid behaviour), "min": [x, y, z], "max": [x, y, z], "spacing": 1.0 (cubic cells; every extent a whole number of spacings), "faces": {"x_min": {"type": "velocity", "velocity": [u, v, w], "profile": "uniform|parabolic"}, "x_max": {"type": "pressure", "pressure": 0}, "y_min": {"type": "wall"} | {"type": "moving_wall", "velocity": [u, v, w] (in the face's plane)} | {"type": "slip"} | {"type": "periodic"} (in opposite pairs), …} (default walls), "obstacles": ["Body1"] (solid bodies of the design or imported surfaces, by name), "surfaces": [{"face": {"body": "Body1", "face": "F2:top"}, "condition": {"type": "moving_wall", "velocity": [u, v, w]}}] (obstacle faces are walls unless listed), "initial_velocity": [u, v, w], "acceleration": [x, y, z] (gravity or a driving force per mass), "tolerance": 0.1 (obstacle triangulation; default a tenth of the spacing), "thermal": {"initial_temperature": 20, "gravity": [0, -9810, 0], "expansion": 2.1e-4 (β per temperature; buoyancy needs both, Boussinesq), "faces": {"y_min": {"type": "temperature", "value": 80} | {"type": "flux", "value": 1000 (heat per area into the fluid)} | {"type": "adiabatic"}, …} (default adiabatic; a velocity inlet carries its held temperature in, a pressure outlet lets heat leave), "surfaces": [{"face": {"body": "Body1", "face": "F2:top"}, "condition": {"type": "temperature", "value": 80}}] (obstacle faces are adiabatic unless listed)} (heat in the fluid, carried by the flow and diffused with the material's conductivity and specific heat; left out, the flow is isothermal)}"#,
     super::forms::FLUID_DOMAIN,
     super::forms::FLUID_DOMAIN_DELETE,
     add_fluid_domain,
@@ -319,6 +319,24 @@ fn fluid_boundary(s: &mut Session, p: &Value) -> Result<Value> {
     Ok(json!({"domain": domain}))
 }
 
+pub const FLUID_THERMAL_BOUNDARY: CommandSpec = CommandSpec {
+    id: "model.fluid_domain.thermal_boundary",
+    label: "Temperature or Heat Flux",
+    params: r#"{"domain": "Flow", "at": "x_min" (a box face) | {"body": "Body1", "face": "F2:top"} (an obstacle face), "condition": {"type": "temperature", "value": 80} | {"type": "flux", "value": 1000 (heat per area into the fluid)} | {"type": "adiabatic"}} — sets one boundary's thermal condition; the domain needs `thermal` settings (model.fluid_domain.edit)"#,
+    form: Some(super::forms::FLUID_THERMAL_BOUNDARY),
+    kind: CommandKind::Model,
+    requires: &[Requirement::FluidDomain],
+    run: fluid_thermal_boundary,
+};
+
+fn fluid_thermal_boundary(s: &mut Session, p: &Value) -> Result<Value> {
+    let domain: &str = req(p, "domain")?;
+    let at: FluidBoundary = req(p, "at")?;
+    let condition: ThermalCondition = req(p, "condition")?;
+    s.model.set_fluid_thermal_boundary(domain, at, condition)?;
+    Ok(json!({"domain": domain}))
+}
+
 pub const STEP_CREATE: CommandSpec = CommandSpec {
     id: "model.step.create",
     label: "New Step",

crates/engine/src/cmd/forms.rs

@@ -608,10 +608,59 @@ const SURFACE_CONDITION: &[Variant] = &[
     ),
 ];
 
+const THERMAL_CONDITION: &[Variant] = &[
+    variant("adiabatic", "Adiabatic", &[]),
+    variant("temperature", "Temperature", &[num("value", "Temperature", D::Temperature).help("The temperature held at the boundary.").req()]),
+    variant("flux", "Heat flux", &[num("value", "Heat flux", D::HeatFlux).help("Heat per area into the fluid (negative out of it).").req()]),
+];
+
 const fn face(key: &'static str, label: &'static str) -> Field {
     variants(key, label, FACE_CONDITION).hint("wall")
 }
 
+const fn thermal_face(key: &'static str, label: &'static str) -> Field {
+    variants(key, label, THERMAL_CONDITION).hint("adiabatic")
+}
+
+const FLUID_THERMAL: Field = group(
+    "thermal",
+    "Heat",
+    &[
+        num("initial_temperature", "Initial temperature", D::Temperature)
+            .help("The temperature everywhere at the start, and the reference of the buoyancy.")
+            .def("0"),
+        xyz("gravity", "Gravity", D::Acceleration).help("Gravity for the buoyancy; zero for none. Needs the expansion coefficient."),
+        positive("expansion", "Expansion β", D::Expansion).help("Volumetric thermal expansion coefficient of the fluid, for the buoyancy."),
+        group(
+            "faces",
+            "Box faces",
+            &[
+                thermal_face("x_min", "X min").help("Thermal condition on the box face at the smallest X."),
+                thermal_face("x_max", "X max").help("Thermal condition on the box face at the largest X."),
+                thermal_face("y_min", "Y min").help("Thermal condition on the box face at the smallest Y."),
+                thermal_face("y_max", "Y max").help("Thermal condition on the box face at the largest Y."),
+                thermal_face("z_min", "Z min").help("Thermal condition on the box face at the smallest Z."),
+                thermal_face("z_max", "Z max").help("Thermal condition on the box face at the largest Z."),
+            ],
+        )
+        .help("The thermal condition on each face of the box; a velocity inlet carries its held temperature in."),
+        rows(
+            "surfaces",
+            "Obstacle faces",
+            &[
+                OBSTACLE_FACE,
+                variants("condition", "Condition", THERMAL_CONDITION)
+                    .help("The thermal condition on that face.")
+                    .req()
+                    .def(r#"{"type": "adiabatic"}"#),
+            ],
+        )
+        .help("Thermal conditions on single obstacle faces; every other obstacle face is adiabatic."),
+    ],
+)
+.help("Heat carried by the flow: needs the material's conductivity and specific heat. Off for an isothermal flow.")
+.hint("isothermal when off");
+
 pub const FLUID_DOMAIN: &[Field] = &[
     NAME,
     reference("material", "Fluid", C::Material).help("The fluid material: needs a density and a viscosity.").req(),
@@ -649,6 +698,7 @@ pub const FLUID_DOMAIN: &[Field] = &[
     positive("tolerance", "Obstacle tolerance", D::Length)
         .help("How close a cell must be to an obstacle surface to count as on it.")
         .hint("a tenth of the cell size"),
+    FLUID_THERMAL,
 ];
 
 pub const FLUID_BOUNDARY: &[Field] = &[
@@ -660,6 +710,15 @@ pub const FLUID_BOUNDARY: &[Field] = &[
     variants("condition", "Condition", FACE_CONDITION).help("The flow condition on that boundary.").req().def(r#"{"type": "velocity"}"#),
 ];
 
+pub const FLUID_THERMAL_BOUNDARY: &[Field] = &[
+    reference("domain", "Domain", C::FluidDomain).help("The fluid domain to change; it needs heat settings.").req(),
+    name_or("at", "Boundary", &BOX_FACE, "Obstacle face", OBSTACLE_FACE_FIELDS)
+        .help("A face of the domain's box, or a named face of an obstacle body.")
+        .req()
+        .def("\"x_min\""),
+    variants("condition", "Condition", THERMAL_CONDITION).help("The thermal condition on that boundary.").req().def(r#"{"type": "temperature"}"#),
+];
+
 const PROCEDURE: Field = variants(
     "procedure",
     "Procedure",
@@ -1972,6 +2031,7 @@ pub const SOLVE_PROBE: &[Field] = &[
     flag("elements", "Element results").help("Also read element results (stresses) around the nodes."),
     points("points", "Fluid points", D::Length, 10_000).help("Points in the fluid to read."),
     json("forces", "Fluid walls").help("Walls whose fluid forces to sum, as JSON.").hint("[\"x_min\", \"F2:top\"]"),
+    json("heat", "Fluid walls (heat)").help("Walls whose heat histories to read, as JSON.").hint("[\"y_min\", \"F2:top\"]"),
 ];
 
 // ---- Results ----

crates/engine/src/cmd/mod.rs

@@ -170,6 +170,7 @@ pub static ALL: &[CommandSpec] = &[
     define::FLUID_DOMAIN_EDIT,
     define::FLUID_DOMAIN_DELETE,
     define::FLUID_BOUNDARY,
+    define::FLUID_THERMAL_BOUNDARY,
     define::STEP_CREATE,
     define::STEP_EDIT,
     define::STEP_DELETE,

crates/engine/src/cmd/solve.rs

@@ -100,7 +100,7 @@ pub const HISTORY: CommandSpec = CommandSpec {
 pub const FLUID: CommandSpec = CommandSpec {
     id: "solve.fluid",
     label: "Solve Fluid Step",
-    params: r#"{"step": "Flow", "backend": "cpu", "max_cells": 16000000, "max_steps": 10000000, "max_cell_updates": 1e12, "max_frame_bytes": 2e9, "threads": 8} — runs a fluid step with the lattice Boltzmann solver (the only fluid step when `step` is left out): obstacles are triangulated from the open design or taken from imported surfaces; returns the time step, relaxation time, Mach number, mass and the final force on every wall; `backend` is `cpu` (the default: double precision, the same result on every machine) or `gpu` (the step on the GPU `FIELDCRAFT_GPU` picks, single precision, within the backend's stated tolerance of the CPU; an error when the build has none); the caps may only be lowered; `threads` as for solve.static"#,
+    params: r#"{"step": "Flow", "backend": "cpu", "max_cells": 16000000, "max_steps": 10000000, "max_cell_updates": 1e12, "max_frame_bytes": 2e9, "threads": 8} — runs a fluid step with the lattice Boltzmann solver (the only fluid step when `step` is left out): obstacles are triangulated from the open design or taken from imported surfaces; returns the time step, relaxation time, Mach number, mass, the final force on every wall and, with heat, the thermal relaxation time, the temperature range and the heat into every wall; `backend` is `cpu` (the default: double precision, the same result on every machine) or `gpu` (the step on the GPU `FIELDCRAFT_GPU` picks, single precision, within the backend's stated tolerance of the CPU; an error when the build has none); the caps may only be lowered; `threads` as for solve.static"#,
     form: Some(super::forms::SOLVE_FLUID),
     kind: CommandKind::Model,
     requires: &[Requirement::FluidDomain, Requirement::Step],
@@ -110,7 +110,7 @@ pub const FLUID: CommandSpec = CommandSpec {
 pub const PROBE: CommandSpec = CommandSpec {
     id: "solve.probe",
     label: "Probe Results",
-    params: r#"{"step": "Load", "nodes": [1, 2] (or "set": "fillet", a set's nodes), "mode": 1, "frame": 0} — the stored result at nodes (with `max_von_mises`, the largest among them, for static steps): displacement, rotation, reaction, reaction moment, stress and von Mises (static, with `elements` each element's results: solid, 2D and truss nodal stresses (and, in a nonlinear step, the equivalent plastic strain at the integration points of a plastic element), shell section forces (N, M, Q per length in the element's axes) and bottom, middle and top stresses, beam end forces (N, V_y, V_z, T, M_y, M_z in its axes)), the shape of mode `mode` with its rotations (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default) at `nodes` (or at the nodes nearest to points: "at": [[x, y, z], ...], lowest label on a tie, each entry then giving its `xyz` and `distance`) and, with `elements` (explicit), each element's stress, von Mises stress, strain, equivalent plastic strain, damage, temperature and whether it has eroded; a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity and pressure in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories"#,
+    params: r#"{"step": "Load", "nodes": [1, 2] (or "set": "fillet", a set's nodes), "mode": 1, "frame": 0} — the stored result at nodes (with `max_von_mises`, the largest among them, for static steps): displacement, rotation, reaction, reaction moment, stress and von Mises (static, with `elements` each element's results: solid, 2D and truss nodal stresses (and, in a nonlinear step, the equivalent plastic strain at the integration points of a plastic element), shell section forces (N, M, Q per length in the element's axes) and bottom, middle and top stresses, beam end forces (N, V_y, V_z, T, M_y, M_z in its axes)), the shape of mode `mode` with its rotations (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default) at `nodes` (or at the nodes nearest to points: "at": [[x, y, z], ...], lowest label on a tie, each entry then giving its `xyz` and `distance`) and, with `elements` (explicit), each element's stress, von Mises stress, strain, equivalent plastic strain, damage, temperature and whether it has eroded; a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity, pressure and (with heat) temperature in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories and {"step": "Flow", "heat": ["y_min"]} the heat into them over time (power, positive into the wall)"#,
     form: Some(super::forms::SOLVE_PROBE),
     kind: CommandKind::Model,
     requires: &[Requirement::Result],
@@ -520,6 +520,21 @@ fn fluid_summary(r: &FluidResult) -> Value {
         })
         .collect();
     let (m0, m1) = (r.frames.first().map_or(0.0, |f| f.mass), last.map_or(0.0, |f| f.mass));
+    let heat: Vec<Value> = r
+        .heat
+        .iter()
+        .filter(|h| h.power.iter().any(|v| *v != 0.0))
+        .map(|h| {
+            let recent = h.power.get(h.power.len().saturating_sub(tail)..).unwrap_or_default();
+            let mean = recent.iter().sum::<f64>() / recent.len().max(1) as f64;
+            json!({"name": h.name, "body": h.body, "final": h.power.last(), "mean": mean})
+        })
+        .collect();
+    let temperature = last.filter(|f| !f.temperature.is_empty()).map(|f| {
+        let fluid = f.temperature.iter().zip(&r.solid).filter(|(_, solid)| !**solid).map(|(t, _)| f64::from(*t));
+        let (lo, hi) = fluid.fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), t| (lo.min(t), hi.max(t)));
+        json!({"min": lo, "max": hi})
+    });
     json!({
         "step": r.step,
         "backend": r.stats.backend,
@@ -538,6 +553,9 @@ fn fluid_summary(r: &FluidResult) -> Value {
         "max_speed": max_speed,
         "mass": {"first": m0, "last": m1, "relative_change": if m0 != 0.0 { (m1 - m0) / m0 } else { 0.0 }},
         "forces": forces,
+        "tau_thermal": r.stats.tau_thermal,
+        "temperature": temperature,
+        "heat": heat,
         "warnings": r.warnings,
     })
 }
@@ -559,6 +577,25 @@ fn probe_fluid(r: &FluidResult, p: &Value) -> Result<Value> {
             .collect::<Result<_>>()?;
         return Ok(json!({"step": r.step, "times": r.history_times, "forces": out}));
     }
+    if let Some(names) = opt::<Vec<String>>(p, "heat")? {
+        if names.len() > MAX_PROBE_FORCES {
+            return Err(EngineError::Params(format!("`heat` lists at most {MAX_PROBE_FORCES} walls")));
+        }
+        if r.heat.is_empty() {
+            return Err(EngineError::Params(format!("step `{}` carries no heat: give its domain `thermal` settings", r.step)));
+        }
+        let out: Vec<Value> = names
+            .iter()
+            .map(|n| {
+                let h =
+                    r.heat.iter().find(|h| &h.name == n || h.body.as_ref().is_some_and(|b| format!("{b}/{}", h.name) == *n)).ok_or_else(|| {
+                        EngineError::Params(format!("no wall `{n}` (box faces x_min … z_max, or obstacle faces by name or body/name)"))
+                    })?;
+                Ok(json!({"name": h.name, "body": h.body, "power": h.power}))
+            })
+            .collect::<Result<_>>()?;
+        return Ok(json!({"step": r.step, "times": r.history_times, "heat": out}));
+    }
     let pts = points(p, "points", MAX_PROBE)?;
     let frame = match opt::<usize>(p, "frame")? {
         Some(k) if k < r.frames.len() => k,
@@ -569,8 +606,12 @@ fn probe_fluid(r: &FluidResult, p: &Value) -> Result<Value> {
     let out: Vec<Value> = pts
         .iter()
         .map(|&pt| match r.probe(frame, pt) {
-            Some((cell, velocity, pressure)) => {
+            Some((cell, velocity, pressure, temperature)) => {
-                json!({"point": pt, "cell": cell, "solid": r.solid.get(cell), "velocity": velocity, "pressure": pressure})
+                let mut v = json!({"point": pt, "cell": cell, "solid": r.solid.get(cell), "velocity": velocity, "pressure": pressure});
+                if let (Some(t), Some(obj)) = (temperature, v.as_object_mut()) {
+                    obj.insert("temperature".into(), json!(t));
+                }
+                v
             }
             None => json!({"point": pt, "cell": null}),
         })
@@ -845,6 +886,7 @@ pub(crate) mod tests {
     }
 
     use super::*;
+    use fieldcraft_model::{BoxFace, FluidBoundary, ThermalCondition};
 
     /// A cantilever 10 × 1 × 1 of ten hex8 elements with incompatible modes, clamped at x = 0,
     /// a total force of 1 down on the free end; a static and a frequency step.
@@ -1110,6 +1152,95 @@ pub(crate) mod tests {
         }
     }
 
+    /// The Couette gap with heat: the sliding top wall held at 1, the bottom at 0, the fluid's
+    /// diffusivity k / (ρ c) = 0.1.
+    fn heated_couette(s: &mut Session) {
+        couette(s);
+        heat_the_couette(s);
+    }
+
+    fn heat_the_couette(s: &mut Session) {
+        s.execute(
+            "model.material.edit",
+            &json!({"name": "Oil", "density": 1.0, "fluid": {"type": "newtonian", "viscosity": 0.1},
+            "thermal": {"conductivity": 0.1, "specific_heat": 1.0}}),
+        )
+        .unwrap();
+        s.execute(
+            "model.fluid_domain.edit",
+            &json!({"name": "Gap", "material": "Oil", "min": [0, 0, 0], "max": [0.125, 1, 0.125], "spacing": 0.125,
+            "faces": {"x_min": {"type": "periodic"}, "x_max": {"type": "periodic"}, "z_min": {"type": "periodic"}, "z_max": {"type": "periodic"},
+                "y_max": {"type": "moving_wall", "velocity": [0.1, 0, 0]}},
+            "thermal": {"initial_temperature": 0.5, "faces": {"y_min": {"type": "temperature", "value": 0}}}}),
+        )
+        .unwrap();
+        s.execute("model.fluid_domain.thermal_boundary", &json!({"domain": "Gap", "at": "y_max", "condition": {"type": "temperature", "value": 1}}))
+            .unwrap();
+    }
+
+    #[test]
+    fn a_fluid_step_carries_heat_between_its_walls() {
+        let mut s = Session::default();
+        heated_couette(&mut s);
+        let r = s.execute("solve.fluid", &json!({})).unwrap();
+        assert!(r["tau_thermal"].as_f64().unwrap() > 0.5, "{r}");
+        assert!((r["temperature"]["min"].as_f64().unwrap() - 0.0625).abs() < 1e-6, "{r}");
+        assert!((r["temperature"]["max"].as_f64().unwrap() - 0.9375).abs() < 1e-6, "{r}");
+        // Conduction across the gap: k ΔT / H over the wall's area, out of the hot wall and
+        // into the cold one.
+        let expected = 0.1 * 1.0 / 1.0 * 0.125 * 0.125;
+        let heat = r["heat"].as_array().unwrap();
+        let cold = heat.iter().find(|h| h["name"] == "y_min").unwrap();
+        let hot = heat.iter().find(|h| h["name"] == "y_max").unwrap();
+        assert!((cold["mean"].as_f64().unwrap() / expected - 1.0).abs() < 1e-6, "{cold}");
+        assert!((hot["mean"].as_f64().unwrap() / expected + 1.0).abs() < 1e-6, "{hot}");
+
+        let probe = s.execute("solve.probe", &json!({"step": "Shear", "points": [[0.06, 0.5625, 0.06]]})).unwrap();
+        assert!((probe["points"][0]["temperature"].as_f64().unwrap() - 0.5625).abs() < 1e-6, "{probe}");
+        let h = s.execute("solve.probe", &json!({"step": "Shear", "heat": ["y_min", "y_max"]})).unwrap();
+        assert_eq!(h["heat"].as_array().unwrap().len(), 2);
+        assert_eq!(h["times"].as_array().unwrap().len(), h["heat"][0]["power"].as_array().unwrap().len());
+        assert!(s.execute("solve.probe", &json!({"step": "Shear", "heat": ["nope"]})).is_err());
+
+        // The thermal boundary needs the domain's heat settings; the material must keep its
+        // heat properties while a domain carries heat.
+        let isothermal = json!({"name": "Gap", "material": "Oil", "min": [0, 0, 0], "max": [0.125, 1, 0.125], "spacing": 0.125});
+        s.execute("model.fluid_domain.edit", &isothermal).unwrap();
+        assert!(
+            s.execute("model.fluid_domain.thermal_boundary", &json!({"domain": "Gap", "at": "y_max", "condition": {"type": "adiabatic"}})).is_err()
+        );
+        let r = s.execute("solve.fluid", &json!({})).unwrap();
+        assert!(r["temperature"].is_null() && r["heat"].as_array().unwrap().is_empty(), "{r}");
+        assert!(s.execute("solve.probe", &json!({"step": "Shear", "heat": ["y_min"]})).is_err());
+        heat_the_couette(&mut s);
+        assert!(s.execute("model.material.edit", &json!({"name": "Oil", "density": 1.0, "fluid": {"type": "newtonian", "viscosity": 0.1}})).is_err());
+        // JSON has no NaN (`json!(f64::NAN)` is `null`): the model refuses one directly.
+        let nan = ThermalCondition::Temperature { value: f64::NAN };
+        assert!(s.model.set_fluid_thermal_boundary("Gap", FluidBoundary::Box(BoxFace::YMax), nan).is_err());
+        let hostile = [
+            json!({"domain": "Gap", "at": "y_max", "condition": {"type": "temperature", "value": "nan"}}),
+            json!({"domain": "Gap", "at": "top", "condition": {"type": "temperature", "value": 1}}),
+            json!({"domain": "Gap", "at": {"body": "Nope", "face": "F1"}, "condition": {"type": "temperature", "value": 1}}),
+            json!({"domain": "Nope", "at": "y_max", "condition": {"type": "temperature", "value": 1}}),
+            json!({"domain": "Gap", "at": "y_max", "condition": {"type": "oven"}}),
+            json!({"domain": "Gap", "at": "y_max"}),
+            json!(null),
+        ];
+        for p in &hostile {
+            let before = s.model.clone();
+            assert!(s.execute("model.fluid_domain.thermal_boundary", p).is_err(), "{p}");
+            assert_eq!(s.model, before, "{p}");
+        }
+        // A finite flux the lattice cannot hold, or one on a periodic face, is refused by the
+        // solve, not by the model.
+        for (at, value) in [("y_max", 1e308), ("x_min", 1.0)] {
+            let mut t = s.clone();
+            t.execute("model.fluid_domain.thermal_boundary", &json!({"domain": "Gap", "at": at, "condition": {"type": "flux", "value": value}}))
+                .unwrap();
+            assert!(t.execute("solve.fluid", &json!({})).is_err(), "{at}");
+        }
+    }
+
     #[test]
     fn a_fluid_domain_takes_its_obstacles_from_the_open_design() {
         use fieldcraft_geometry::sample::{Plate, plate_design};
@@ -1177,6 +1308,9 @@ pub(crate) mod tests {
             json!({"step": "Shear", "forces": ["nope"]}),
             json!({"step": "Shear", "forces": vec!["y_min"; 1000]}),
             json!({"step": "Shear", "forces": "y_min"}),
+            json!({"step": "Shear", "heat": ["y_min"]}),
+            json!({"step": "Shear", "heat": vec!["y_min"; 1000]}),
+            json!({"step": "Shear", "heat": 3}),
             json!({"step": "Shear", "nodes": [1]}),
             json!({"step": "Shear", "threads": 0}),
             json!({"step": "Shear", "threads": 100000}),

crates/engine/src/flow.rs

@@ -422,8 +422,8 @@ pub(crate) mod tests {
         let solid: Vec<bool> = (0..n).map(|c| matches!(grid.coords(c), [5 | 6, 1 | 2, 1 | 2])).collect();
         let velocity: Vec<[f32; 3]> = (0..n).map(|c| if solid[c] { [0.0; 3] } else { [grid.coords(c)[0] as f32 + 1.0, 0.0, 0.0] }).collect();
         let pressure: Vec<f32> = (0..n).map(|c| if solid[c] { 0.0 } else { 10.0 - grid.coords(c)[0] as f32 }).collect();
-        let at_rest = Frame { time: 0.0, step: 0, velocity: vec![[0.0; 3]; n], pressure: vec![0.0; n], mass: 1.0 };
+        let at_rest = Frame { time: 0.0, step: 0, velocity: vec![[0.0; 3]; n], pressure: vec![0.0; n], temperature: Vec::new(), mass: 1.0 };
-        let last = Frame { time: 1.0, step: 10, velocity, pressure, mass: 1.0 };
+        let last = Frame { time: 1.0, step: 10, velocity, pressure, temperature: Vec::new(), mass: 1.0 };
         let stats = FluidStats {
             cells: n,
             fluid_cells: n - 8,
@@ -438,6 +438,7 @@ pub(crate) mod tests {
             mach: 0.1,
             cell_updates: 0,
             backend: "cpu".into(),
+            tau_thermal: None,
         };
         FluidResult {
             step: "Flow".into(),
@@ -446,6 +447,7 @@ pub(crate) mod tests {
             frames: vec![at_rest, last],
             history_times: Vec::new(),
             forces: Vec::new(),
+            heat: Vec::new(),
             stats,
             warnings: Vec::new(),
         }

docs/architecture.md

@@ -27,7 +27,7 @@ crates in strictly lower layers.
 | L1 | `fieldcraft-render` | Display data without a GPU dependency: the visible faces, element edges and nodes of a model, the orbit camera, a CPU rasterizer for headless snapshots, geometry surfaces with their face boundaries, deformed shapes and banded contour plots of nodal results, slices, surface values and streamlines of flow fields on grids of cells, picking. | started |
 | L2 | `fieldcraft-solver` | Implicit solver: assembly from the resolved model, boundary conditions, sparse direct solution, linear static and modal analysis, stresses and reactions. | started |
 | L2 | `fieldcraft-explicit` | Explicit dynamics: central-difference time integration with lumped mass, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, rigid walls, large deformation with objective rates, plasticity, damage and erosion, energy account, histories and field frames. | started |
-| L2 | `fieldcraft-lbm` | Lattice Boltzmann flow solver on the CPU: D3Q19 two-relaxation-time scheme, bounce-back and slip walls, Zou–He inlets and outlets, obstacles voxelised from the face-tagged surface with exact predicates, momentum-exchange forces per face. Uses `mesh` for the predicates and the boundary surface. | started |
+| L2 | `fieldcraft-lbm` | Lattice Boltzmann flow solver on the CPU: D3Q19 two-relaxation-time scheme, bounce-back and slip walls, Zou–He inlets and outlets, obstacles voxelised from the face-tagged surface with exact predicates, momentum-exchange forces per face; a D3Q7 temperature lattice carried by the flow (held temperatures, heat fluxes, adiabatic faces, Boussinesq buoyancy, the heat into every face). Uses `mesh` for the predicates and the boundary surface. | started |
 | L3 | `fieldcraft-fsi` | Fluid–structure coupling between the L2 solvers (later). | planned |
 | L4 | `fieldcraft-doc` | The document: geometry, mesh, model, results, history. | planned |
 | L5 | `fieldcraft-io` | `.fieldcraft` files (docs/formats/fieldcraft.md): text model definition, binary mesh, geometry references with content hashes, the journal; atomic saves; autosave storage for crash recovery; read and write roots that confine an automation client's file access (`roots`). Bytes in, bytes out; the file system only in `fs`, `recovery` and `roots` (in the browser, files kept in memory: `memory`). Result files later. | started |
@@ -309,21 +309,26 @@ cantilever's free-vibration period against `solve.frequency`.
 
 ## Lattice Boltzmann solver
 
-`fieldcraft-lbm` solves a model's fluid steps (`solve.fluid`; `solve.probe` reads velocity and
+`fieldcraft-lbm` solves a model's fluid steps (`solve.fluid`; `solve.probe` reads velocity,
-pressure at points and wall force histories). Papers are cited in the module docs and in
+pressure and temperature at points and wall force and heat histories). Papers are cited in the
-`docs/solvers/literature.md` (Part I).
+module docs and in `docs/solvers/literature.md` (Part I).
 
-1. **Model.** A fluid material (density, Newtonian viscosity), a fluid domain (box, cubic cell
+1. **Model.** A fluid material (density, Newtonian viscosity; conductivity and specific heat
-   spacing, a condition on each box face: wall, moving wall, slip, velocity inlet with a uniform
+   for heat), a fluid domain (box, cubic cell spacing, a condition on each box face: wall,
-   or parabolic profile, pressure outlet, periodic pair; obstacle bodies and moving-wall
+   moving wall, slip, velocity inlet with a uniform or parabolic profile, pressure outlet,
-   conditions on their named faces; initial velocity; body acceleration) and a fluid step (end
+   periodic pair; obstacle bodies and moving-wall conditions on their named faces; initial
+   velocity; body acceleration; optional `thermal` settings: the initial temperature, a thermal
+   condition per box face and per named obstacle face — held temperature, heat flux or
+   adiabatic — and gravity with the expansion coefficient for buoyancy) and a fluid step (end
    time or steps, optional time step, frames, the TRT parameter Λ). Deck writers refuse fluid
    steps.
 2. **Setup.** Lattice units from the spacing, the time step and the density. Without a given
    time step the solver keeps the characteristic lattice velocity at 0.05 and τ ≤ 1. A lattice
    Mach number above 0.3, τ below 0.505 or a lattice acceleration above 0.01 is refused with the
    numbers and what to change; caps on cells, steps, cell updates and frame memory are checked
-   before anything is allocated.
+   before anything is allocated. With heat, the scalar relaxation time follows from the
+   diffusivity k / (ρ c) in lattice units (refused below 0.505, warned above 2) and a flux on an
+   inlet, outlet or periodic face is refused.
 3. **Grid.** The engine triangulates each obstacle body (or takes an imported surface) and
    prepares its `BoundarySurface`; cells whose centres have a positive winding number are solid,
    decided with exact `orient2d`/`orient3d` and one symbolic perturbation, so faces through cell
@@ -332,17 +337,27 @@ pressure at points and wall force histories). Papers are cited in the module doc
    direction with per-link rules for the cells next to walls and box faces; Zou–He on inlet and
    outlet layers. Wall links give the momentum exchanged with every box face and B-rep face
    (the face of the triangle nearest the link). A non-positive density or a lattice speed above
-   0.5 stops the run with an error.
+   0.5 stops the run with an error. With heat, a D3Q7 BGK scalar lattice (weights ¼ and ⅛,
-5. **Results.** Frames of velocity and gauge pressure per cell (single precision), force
+   c_s² = ¼) steps after the flow on the same links: it is advected by the cell velocity and
-   histories per wall and the fluid mass, in model units, stored per step like the implicit
+   diffuses with the scalar relaxation time; held temperatures are anti-bounce-back, adiabatic
-   results.
+   faces bounce-back, a flux is added to the bounced population; velocity inlets carry the
+   held temperature in at the cell's velocity, pressure outlets copy their neighbour. The heat
+   crossing every wall link is summed per face as the flow's force is. The temperature feeds
+   back through a per-cell Boussinesq force −β (T − T₀) g.
+5. **Results.** Frames of velocity, gauge pressure and (with heat) temperature per cell (single
+   precision), force and heat histories per wall and the fluid mass, in model units, stored
+   per step like the implicit results.
 
 Verification (`crates/lbm/tests/verification.rs`, analytic): plane Poiseuille exact to round-off
 with Λ = 3/16 and second order with Λ = 1/4; Couette exact with the wall shear μU/H; Taylor–Green
 decay second order; Poiseuille in a staircase pipe; Stokes drag of a simple cubic sphere array
 against Hasimoto (1959); a parabolic-inlet, pressure-outlet channel against Poiseuille; mass
-conserved to round-off; bit-for-bit reproducibility. `crates/lbm/tests/benchmark.rs` (ignored,
+conserved to round-off; bit-for-bit reproducibility. `crates/lbm/tests/thermal.rs` (heat):
-published reference values): Schäfer–Turek 2D-1.
+linear conduction between held plates to round-off; a flux wall against q H / k and the wall's
+heat balance; an inlet temperature carried through a plug flow; a periodic blob that keeps its
+total and returns after one lap; (ignored, minutes) the de Vahl Davis cavity at Ra 10³ and 10⁴
+against the published Nusselt numbers (both within 0.1 %). `crates/lbm/tests/benchmark.rs` (ignored, published
+reference values): Schäfer–Turek 2D-1.
 
 ## Commands
 

docs/features.md

@@ -108,7 +108,7 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands.
 | Boundary Conditions | Boundary Condition — Fixed, displacement or symmetry condition on a set or reference point | ready | `model.boundary_condition.create`, `model.boundary_condition.edit`, `model.boundary_condition.delete`, `model.boundary_condition.symmetry` |
 | Loads | Load — Force, moment, pressure, traction, gravity or body force | ready | `model.load.create`, `model.load.edit`, `model.load.delete` |
 | Flow | Fluid Domain — The box of fluid cells, its fluid, and the obstacle bodies voxelised from their face-tagged surfaces | ready | `model.fluid_domain.create`, `model.fluid_domain.edit`, `model.fluid_domain.delete` |
-| Flow | Inlet & Outlet — Velocity inlets, pressure outlets, walls, moving and slip walls, periodic faces | ready | `model.fluid_domain.boundary` |
+| Flow | Inlet & Outlet — Velocity inlets, pressure outlets, walls, moving and slip walls, periodic faces; held temperatures, heat fluxes and adiabatic faces | ready | `model.fluid_domain.boundary`, `model.fluid_domain.thermal_boundary` |
 
 ## Interactions
 
@@ -141,7 +141,7 @@ the control channel and AI agents over MCP (docs/mcp.md) run the same commands.
 | Implicit | Frequency — Natural frequencies and mode shapes of a frequency step | ready | `solve.frequency` |
 | Explicit | Explicit — Transient dynamics of an explicit step by central differences (contact and erosion coming in M10) | ready | `solve.explicit` |
 | Explicit | Blast & Impact — EOS-driven fluids with fluid–structure contact (ALE or SPH) | coming soon (M14) |  |
-| Fluids | Lattice Boltzmann — Lattice Boltzmann flow of a fluid step on the CPU, with wall forces | ready | `solve.fluid` |
+| Fluids | Lattice Boltzmann — Lattice Boltzmann flow of a fluid step on the CPU or a GPU, with wall forces, and the heat it carries | ready | `solve.fluid` |
 | Fluids | Fluid–Structure — Coupled fluid and structure solve | coming soon (M13) |  |
 
 ## Results