lbm: a temperature lattice on the flow

A D3Q7 passive scalar rides on the D3Q19 flow: it is advected by the
cell velocity, diffuses with the fluid's thermal diffusivity (the
relaxation time follows from k / (rho c)), and feeds back into the
momentum through a per-cell Boussinesq force when the domain has an
expansion coefficient and gravity. Walls and obstacle faces hold a
temperature (anti-bounce-back), a heat flux (added to the bounced
population) or are adiabatic (bounce-back); velocity inlets carry the
held temperature in, pressure outlets copy their neighbour.

The heat crossing each face group is accumulated during streaming and
reported in physical units (W) in the result's heat histories, next to
the force histories; frames gain a temperature field and the stats
report the thermal relaxation time.

Closed-form checks: linear conduction between held plates, a flux wall
against q H / k with the wall's heat balance, an inlet temperature
carried through a plug flow, and a periodic blob that conserves its
total and returns after one lap. The de Vahl Davis cavity (Ra 1e3 and
1e4 against the published Nusselt numbers) is there as an ignored test
because it takes minutes.

The GPU backend declines domains with heat for now and asks for the
CPU; its kernels follow.

Part of #157.

crates/lbm-gpu/src/native.rs

@@ -113,6 +113,9 @@ impl GpuFlow {
     /// `flow` (its grid, rules and initial state: everything the step needs) on the device
     /// `policy` picks.
     pub fn new(flow: &Flow, policy: Policy) -> Result<GpuFlow> {
+        if flow.has_heat() {
+            return Err(Error::Invalid("heat in the fluid is not on the GPU yet: solve this step on the CPU (`backend: \"cpu\"`)".into()));
+        }
         let Some((adapter, tier)) = select::adapter(policy).map_err(gpu_error)? else {
             return Err(Error::Invalid(format!("{}=cpu: no GPU device; solve with `backend: \"cpu\"`", select::ENV)));
         };
@@ -528,7 +531,12 @@ impl Backend for GpuFlow {
         let m = self.read_f32(&self.moments, 4 * self.n)?;
         let density = m.as_chunks::<4>().0.iter().map(|v| f64::from(v[3])).collect();
         let velocity = m.as_chunks::<4>().0.iter().map(|v| [f64::from(v[0]), f64::from(v[1]), f64::from(v[2])]).collect();
-        Ok(Fields { density, velocity })
+        // No temperature lattice on the device yet (the CPU carries heat).
+        Ok(Fields { density, velocity, temperature: Vec::new() })
+    }
+
+    fn heat(&mut self) -> Result<Vec<f64>> {
+        Ok(Vec::new())
     }
 
     fn mass(&mut self) -> Result<f64> {

crates/lbm/src/backend.rs

@@ -13,11 +13,13 @@ use crate::flow::{Flow, Group};
 use crate::setup::Setup;
 use crate::{Error, Result};
 
-/// The density and velocity of every cell (lattice units; zero in solid cells).
+/// The density and velocity of every cell (lattice units; zero in solid cells), and the
+/// temperature (model units) when the flow carries heat (empty otherwise).
 #[derive(Debug, Clone, PartialEq)]
 pub struct Fields {
     pub density: Vec<f64>,
     pub velocity: Vec<[f64; 3]>,
+    pub temperature: Vec<f64>,
 }
 
 /// A stepper of one prepared flow.
@@ -35,6 +37,9 @@ pub trait Backend {
     fn mass(&mut self) -> Result<f64>;
     /// The momentum each force group received in the last step (lattice units: force per step).
     fn forces(&mut self) -> Result<Vec<[f64; 3]>>;
+    /// The heat each force group received in the last step (lattice units: temperature times
+    /// cell faces per step); empty without heat.
+    fn heat(&mut self) -> Result<Vec<f64>>;
     fn groups(&self) -> &[Group];
     fn fluid_cells(&self) -> usize;
     fn is_fluid(&self, cell: usize) -> bool;

crates/lbm/src/flow.rs

@@ -27,6 +27,20 @@
 //!   65:041203, 2002): a population `f*` leaving a fluid cell toward the wall and the one `f'`
 //!   coming back transfer `c_out f* − c_in f'` to the wall.
 //!
+//! # Heat
+//! With thermal settings a second lattice carries the temperature: D3Q7 (the rest direction and
+//! the six axis directions, the first seven of D3Q19), weights ¼ and ⅛, `c_s² = ¼`, equilibrium
+//! `w_q T (1 + 4 c_q·u)` with the velocity the momentum lattice just collided to, BGK collision
+//! with `τ_T = 4 α + ½` (Krüger et al. 2017, ch. 8: the passive-scalar double-distribution
+//! model). Its streaming follows the momentum lattice's links: a wall held at a temperature is
+//! anti-bounce-back (`g_q = −g*_q̄ + 2 w_q T_w`, which puts the temperature halfway between the
+//! cells), an adiabatic wall bounces back, a heat flux through a wall face adds the flux to the
+//! bounced population (one axis direction crosses each wall face, so the cell gains exactly that
+//! heat), an inlet holds its temperature (equilibrium at the cell's velocity) and an outlet takes
+//! its neighbour's. Buoyancy is the Boussinesq force per mass `−β (T − T₀) g`, added to the body
+//! force of every cell (Guo forcing as for the constant acceleration). The heat each wall face
+//! receives is summed over its links like the momentum exchange.
+//!
 //! # Layout and determinism
 //! Two arrays of populations, structure of arrays (`f[q n + cell]`), A–B pattern. Streaming
 //! pulls: every population of a cell comes from one source (a neighbour, a periodic image, the
@@ -36,7 +50,7 @@
 //! only its own populations; the wall forces are summed in link order. There are no parallel
 //! reductions, so the result is the same bit for bit on any number of threads and platform.
 
-use fieldcraft_model::{BoxFace, FaceCondition, Profile};
+use fieldcraft_model::{BoxFace, FaceCondition, Profile, ThermalCondition};
 
 use fieldcraft_parallel as parallel;
 
@@ -95,6 +109,40 @@ pub struct Group {
     pub body: Option<String>,
 }
 
+/// The scalar lattice: directions 0..7 of D3Q19 (rest and the axes).
+const QT: usize = 7;
+/// Its weights.
+const WT: [f64; QT] = [0.25, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125];
+
+/// The scalar equilibrium `w_q T (1 + c_q·u / c_s²)`, `c_s² = ¼`.
+#[inline(always)]
+fn equilibrium_t(q: usize, t: f64, u: [f64; 3]) -> f64 {
+    let cq = c(q);
+    WT.get(q).copied().unwrap_or(0.0) * t * (1.0 + 4.0 * (cq[0] * u[0] + cq[1] * u[1] + cq[2] * u[2]))
+}
+
+/// The temperature lattice and what it needs (see the module notes).
+#[derive(Debug, Clone)]
+pub struct ThermalState {
+    /// `g[q stride + cell]` and the streaming target, A–B like the momentum lattice.
+    g: Vec<f64>,
+    g2: Vec<f64>,
+    /// The temperature of every cell after the last step (the initial one in solid cells).
+    temperature: Vec<f64>,
+    omega: f64,
+    initial: f64,
+    /// `−β g` in lattice units; zero without buoyancy.
+    buoyancy: [f64; 3],
+    /// Each force group's thermal condition (fluxes in lattice units).
+    groups: Vec<ThermalCondition>,
+    /// Per open cell: the temperature it holds (an inlet), or `None` to take its neighbour's
+    /// (an outlet).
+    open: Vec<Option<f64>>,
+    /// The heat each group received in the last step (lattice units: temperature times cell
+    /// faces per step).
+    heat: Vec<f64>,
+}
+
 /// The flow state.
 #[derive(Debug, Clone)]
 pub struct Flow {
@@ -120,6 +168,13 @@ pub struct Flow {
     steps_done: u64,
     /// Obstacles that cover no cell centre (smaller than the spacing, or outside the box).
     empty_obstacles: Vec<String>,
+    /// The temperature lattice, with thermal settings.
+    thermal: Option<ThermalState>,
+    /// With buoyancy: every cell's body acceleration (the constant one plus the buoyant one),
+    /// refreshed from the temperature before each step.
+    cell_force: Vec<[f64; 3]>,
+    /// With heat: every cell's velocity from the last collision, for the scalar equilibrium.
+    cell_u: Vec<[f64; 3]>,
 }
 
 const NONE: u32 = u32::MAX;
@@ -328,6 +383,57 @@ impl Flow {
                 f[q * stride + cell] = if fluid[cell] { equilibrium(q, 1.0, setup.initial_velocity) } else { 0.0 };
             }
         }
+
+        // Heat: the scalar lattice at the initial temperature, each group's condition (the
+        // box faces, then every obstacle face), each open cell's temperature.
+        let (thermal, cell_force, cell_u) = match &setup.thermal {
+            None => (None, Vec::new(), Vec::new()),
+            Some(ts) => {
+                let mut group_conditions: Vec<ThermalCondition> = ts.faces.to_vec();
+                for faces in &setup.obstacles {
+                    group_conditions.extend(faces.thermal.iter().copied());
+                }
+                group_conditions.resize(groups.len(), ThermalCondition::Adiabatic {});
+                let open_t: Vec<Option<f64>> = open
+                    .iter()
+                    .map(|oc| {
+                        let face = BoxFace::of(oc.axis, oc.inward < 0.0);
+                        match (oc.target, ts.faces.get(face.index()).copied().unwrap_or_default()) {
+                            (_, ThermalCondition::Temperature { value }) => Some(value),
+                            (Target::Velocity(_), _) => Some(ts.initial_temperature),
+                            (Target::Density(_), _) => None,
+                        }
+                    })
+                    .collect();
+                let t_n = stride.checked_mul(QT).ok_or_else(|| Error::Limit("the grid is too large".into()))?;
+                let mut tg = filled(t_n, 0.0, "temperature populations")?;
+                let tg2 = filled(t_n, 0.0, "temperature populations")?;
+                for cell in 0..n {
+                    if fluid[cell] {
+                        for q in 0..QT {
+                            tg[q * stride + cell] = equilibrium_t(q, ts.initial_temperature, setup.initial_velocity);
+                        }
+                    }
+                }
+                let buoyant = ts.buoyancy != [0.0; 3];
+                let state = ThermalState {
+                    g: tg,
+                    g2: tg2,
+                    temperature: filled(n, ts.initial_temperature, "temperature")?,
+                    omega: 1.0 / ts.tau,
+                    initial: ts.initial_temperature,
+                    buoyancy: ts.buoyancy,
+                    groups: group_conditions,
+                    open: open_t,
+                    heat: vec![0.0; groups.len()],
+                };
+                (
+                    Some(state),
+                    if buoyant { filled(n, setup.acceleration, "cell forces")? } else { Vec::new() },
+                    filled(n, [0.0; 3], "cell velocities")?,
+                )
+            }
+        };
         Ok(Flow {
             grid,
             n,
@@ -346,9 +452,58 @@ impl Flow {
             fluid_cells,
             steps_done: 0,
             empty_obstacles,
+            thermal,
+            cell_force,
+            cell_u,
         })
     }
 
+    /// The body acceleration of a cell: the constant one, plus the buoyant one with heat.
+    #[inline(always)]
+    fn force_at(&self, cell: usize) -> [f64; 3] {
+        self.cell_force.get(cell).copied().unwrap_or(self.acceleration)
+    }
+
+    /// The temperature of a cell (model units), or `None` without heat.
+    pub fn temperature(&self, cell: usize) -> Option<f64> {
+        self.thermal.as_ref().and_then(|t| t.temperature.get(cell).copied())
+    }
+
+    /// Whether the flow carries heat.
+    pub fn has_heat(&self) -> bool {
+        self.thermal.is_some()
+    }
+
+    /// The heat each group received in the last step (lattice units), empty without heat.
+    pub fn heat(&self) -> &[f64] {
+        self.thermal.as_ref().map_or(&[], |t| &t.heat)
+    }
+
+    /// Set every fluid cell's temperature (and its populations to the equilibrium at its
+    /// velocity). Nothing without heat.
+    pub fn set_temperature(&mut self, temperature: impl Fn(usize) -> f64) {
+        let n = self.n;
+        let s = self.stride;
+        let velocities: Vec<[f64; 3]> = (0..n).map(|c| self.moments(c).map_or([0.0; 3], |m| m.1)).collect();
+        let fluid = &self.fluid;
+        let Some(t) = &mut self.thermal else { return };
+        for cell in 0..n {
+            if !fluid.get(cell).copied().unwrap_or(false) {
+                continue;
+            }
+            let value = temperature(cell);
+            if let Some(slot) = t.temperature.get_mut(cell) {
+                *slot = value;
+            }
+            let u = velocities.get(cell).copied().unwrap_or([0.0; 3]);
+            for q in 0..QT {
+                if let Some(g) = t.g.get_mut(q * s + cell) {
+                    *g = equilibrium_t(q, value, u);
+                }
+            }
+        }
+    }
+
     pub fn grid(&self) -> &Grid {
         &self.grid
     }
@@ -419,7 +574,8 @@ impl Flow {
                 continue;
             }
             let (rho, u) = state(cell);
-            let bare = [0, 1, 2].map(|a| u[a] - 0.5 * self.acceleration[a]);
+            let force = self.force_at(cell);
+            let bare = [0, 1, 2].map(|a| u[a] - 0.5 * force[a]);
             for q in 0..Q {
                 self.f[q * s + cell] = equilibrium(q, rho, bare);
             }
@@ -441,7 +597,8 @@ impl Flow {
             j[1] += fq * cq[1];
             j[2] += fq * cq[2];
         }
-        let u = [0, 1, 2].map(|a| (j[a] + 0.5 * rho * self.acceleration[a]) / rho);
+        let force = self.force_at(cell);
+        let u = [0, 1, 2].map(|a| (j[a] + 0.5 * rho * force[a]) / rho);
         Some((rho, u))
     }
 
@@ -461,6 +618,7 @@ impl Flow {
     /// One time step. Fails with [`Error::Unstable`] when a density stops being a positive
     /// number or a lattice speed exceeds [`MAX_LATTICE_SPEED`].
     pub fn step(&mut self) -> Result<()> {
+        self.refresh_buoyancy();
         let max_u2 = self.collide();
         if !(max_u2 <= MAX_LATTICE_SPEED * MAX_LATTICE_SPEED) {
             return Err(Error::Unstable(format!(
@@ -472,10 +630,177 @@ impl Flow {
         self.stream();
         self.open_boundaries();
         std::mem::swap(&mut self.f, &mut self.g);
+        if self.thermal.is_some() {
+            self.collide_thermal();
+            self.stream_thermal();
+            self.open_boundaries_thermal();
+            if let Some(t) = &mut self.thermal {
+                std::mem::swap(&mut t.g, &mut t.g2);
+            }
+            if !self.update_temperature() {
+                return Err(Error::Unstable(format!(
+                    "the heat went unstable at time step {}: a temperature stopped being a number; lower the time step, refine the spacing or check the heat fluxes",
+                    self.steps_done + 1
+                )));
+            }
+        }
         self.steps_done += 1;
         Ok(())
     }
 
+    /// With buoyancy: every cell's acceleration from its temperature.
+    fn refresh_buoyancy(&mut self) {
+        let Some(t) = &self.thermal else { return };
+        if self.cell_force.is_empty() {
+            return;
+        }
+        let (a, b, t0) = (self.acceleration, t.buoyancy, t.initial);
+        for (force, temp) in self.cell_force.iter_mut().zip(&t.temperature) {
+            let dt = temp - t0;
+            *force = [a[0] + b[0] * dt, a[1] + b[1] * dt, a[2] + b[2] * dt];
+        }
+    }
+
+    /// BGK collision of the scalar lattice, in place, with the velocity of the momentum
+    /// collision.
+    fn collide_thermal(&mut self) {
+        let n = self.n;
+        let s = self.stride;
+        let fluid = &self.fluid;
+        let cell_u = &self.cell_u;
+        let Some(t) = &mut self.thermal else { return };
+        let omega = t.omega;
+        let mut blocks = cell_blocks(&mut t.g, s, n);
+        parallel::for_each_mut(&mut blocks, |b, g| {
+            let len = g.first().map_or(0, |d| d.len());
+            for i in 0..len {
+                let cell = b * CELLS + i;
+                if !fluid.get(cell).copied().unwrap_or(false) {
+                    continue;
+                }
+                let mut temp = 0.0;
+                for q in 0..QT {
+                    temp += g.get(q).and_then(|d| d.get(i)).copied().unwrap_or(0.0);
+                }
+                let u = cell_u.get(cell).copied().unwrap_or([0.0; 3]);
+                for q in 0..QT {
+                    if let Some(x) = g.get_mut(q).and_then(|d| d.get_mut(i)) {
+                        *x -= omega * (*x - equilibrium_t(q, temp, u));
+                    }
+                }
+            }
+        });
+    }
+
+    /// Pull streaming of the scalar lattice into `g2` with the wall rules of the links (the
+    /// first seven directions), then the heat received per group.
+    fn stream_thermal(&mut self) {
+        let n = self.n;
+        let s = self.stride;
+        let (links, offset) = (&self.links, &self.offset);
+        let Some(t) = &mut self.thermal else { return };
+        let (g, groups) = (&t.g, &t.groups);
+        let mut blocks = cell_blocks(&mut t.g2, s, n);
+        parallel::for_each_mut(&mut blocks, |b, out| {
+            let start = b * CELLS;
+            for (q, dst) in out.iter_mut().enumerate().take(QT) {
+                let off = offset.get(q).copied().unwrap_or(0);
+                let src = g.get(q * s..q * s + n).unwrap_or_default();
+                for (i, x) in dst.iter_mut().enumerate() {
+                    if let Some(from) = (start + i).checked_add_signed(-off)
+                        && let Some(v) = src.get(from)
+                    {
+                        *x = *v;
+                    }
+                }
+            }
+            let first = links.partition_point(|l| (l.cell as usize) < start);
+            for link in links.get(first..).unwrap_or_default().iter().take_while(|l| (l.cell as usize) < start + CELLS) {
+                let (cell, q) = (link.cell as usize, link.q as usize);
+                if q >= QT {
+                    continue;
+                }
+                let v = thermal_rule(g, s, groups, cell, q, link.rule);
+                if let (Some(v), Some(x)) = (v, out.get_mut(q).and_then(|d| d.get_mut(cell - start))) {
+                    *x = v;
+                }
+            }
+        });
+        t.heat.fill(0.0);
+        for link in links {
+            let (cell, q) = (link.cell as usize, link.q as usize);
+            if q >= QT {
+                continue;
+            }
+            if let (Rule::Bounce { group, .. } | Rule::Reflect { group, .. }, Some(back)) =
+                (link.rule, thermal_rule(g, s, groups, cell, q, link.rule))
+            {
+                let out = g.get(opposite(q) * s + cell).copied().unwrap_or(0.0);
+                if let Some(h) = t.heat.get_mut(group as usize) {
+                    // what left the fluid toward the wall minus what came back
+                    *h += out - back;
+                }
+            }
+        }
+    }
+
+    /// Inlets hold their temperature; outlets take their neighbour's populations.
+    fn open_boundaries_thermal(&mut self) {
+        let s = self.stride;
+        let (open, offset, cell_u) = (&self.open, &self.offset, &self.cell_u);
+        let Some(t) = &mut self.thermal else { return };
+        for (oc, held) in open.iter().zip(&t.open) {
+            let cell = oc.cell as usize;
+            match held {
+                Some(value) => {
+                    let u = cell_u.get(cell).copied().unwrap_or([0.0; 3]);
+                    for q in 0..QT {
+                        if let Some(x) = t.g2.get_mut(q * s + cell) {
+                            *x = equilibrium_t(q, *value, u);
+                        }
+                    }
+                }
+                None => {
+                    // the neighbour one cell inward along the face's axis
+                    let axis_dir = 1 + 2 * oc.axis + usize::from(oc.inward < 0.0);
+                    let off = offset.get(axis_dir).copied().unwrap_or(0);
+                    if let Some(from) = cell.checked_add_signed(off) {
+                        for q in 0..QT {
+                            let v = t.g2.get(q * s + from).copied();
+                            if let (Some(v), Some(x)) = (v, t.g2.get_mut(q * s + cell)) {
+                                *x = v;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    /// Every fluid cell's temperature from its populations.
+    /// The temperature of every fluid cell from its populations; false when one is not finite.
+    fn update_temperature(&mut self) -> bool {
+        let n = self.n;
+        let s = self.stride;
+        let fluid = &self.fluid;
+        let Some(t) = &mut self.thermal else { return true };
+        let mut finite = true;
+        for cell in 0..n {
+            if !fluid.get(cell).copied().unwrap_or(false) {
+                continue;
+            }
+            let mut temp = 0.0;
+            for q in 0..QT {
+                temp += t.g.get(q * s + cell).copied().unwrap_or(0.0);
+            }
+            finite &= temp.is_finite();
+            if let Some(slot) = t.temperature.get_mut(cell) {
+                *slot = temp;
+            }
+        }
+        finite
+    }
+
     /// Two-relaxation-time collision with Guo forcing, in place, by blocks of [`CELLS`] cells
     /// in parallel. Returns the largest squared lattice speed, or NaN when a density is not
     /// positive.
@@ -483,19 +808,27 @@ impl Flow {
         let n = self.n;
         let (wp, wm) = (self.omega_plus, self.omega_minus);
         let (sp, sm) = (1.0 - 0.5 * wp, 1.0 - 0.5 * wm);
-        let a = self.acceleration;
+        let constant = self.acceleration;
-        let forced = a != [0.0; 3];
+        let cell_force = &self.cell_force;
+        let forced = constant != [0.0; 3] || !cell_force.is_empty();
         let fluid = &self.fluid;
         let mut blocks = cell_blocks(&mut self.f, self.stride, n);
         let mut maxima = vec![0.0f64; blocks.len()];
-        let mut work: Vec<(&mut Vec<&mut [f64]>, &mut f64)> = blocks.iter_mut().zip(maxima.iter_mut()).collect();
+        // with heat, each block's velocities go to its slice of cell_u (empty otherwise)
-        parallel::for_each_mut(&mut work, |b, (f, max_u2)| {
+        let mut u_blocks: Vec<&mut [[f64; 3]]> = self.cell_u.chunks_mut(CELLS).collect();
+        let mut work: Vec<(&mut Vec<&mut [f64]>, &mut f64, Option<&mut [[f64; 3]]>)> = Vec::with_capacity(blocks.len());
+        let mut u_iter = u_blocks.drain(..);
+        for (f, m) in blocks.iter_mut().zip(maxima.iter_mut()) {
+            work.push((f, m, u_iter.next()));
+        }
+        parallel::for_each_mut(&mut work, |b, (f, max_u2, u_out)| {
             let mut fq = [0.0; Q];
             let len = f.first().map_or(0, |d| d.len());
             for i in 0..len {
                 if !fluid.get(b * CELLS + i).copied().unwrap_or(false) {
                     continue;
                 }
+                let a = cell_force.get(b * CELLS + i).copied().unwrap_or(constant);
                 let (mut rho, mut j) = (0.0, [0.0; 3]);
                 for q in 0..Q {
                     let v = f.get(q).and_then(|d| d.get(i)).copied().unwrap_or(0.0);
@@ -514,6 +847,9 @@ impl Flow {
                 let u = [(j[0] + 0.5 * force[0]) / rho, (j[1] + 0.5 * force[1]) / rho, (j[2] + 0.5 * force[2]) / rho];
                 let uu = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
                 **max_u2 = max_u2.max(uu);
+                if let Some(slot) = u_out.as_mut().and_then(|u| u.get_mut(i)) {
+                    *slot = u;
+                }
                 let uf = u[0] * force[0] + u[1] * force[1] + u[2] * force[2];
                 let mut put = |q: usize, v: f64| {
                     if let Some(x) = f.get_mut(q).and_then(|d| d.get_mut(i)) {
@@ -624,9 +960,9 @@ impl Flow {
     /// Zou–He velocity and pressure boundaries on the streamed populations.
     fn open_boundaries(&mut self) {
         let s = self.stride;
-        let a_half = self.acceleration.map(|v| 0.5 * v);
         for oc in &self.open {
             let cell = oc.cell as usize;
+            let a_half = self.force_at(cell).map(|v| 0.5 * v);
             let axis = oc.axis;
             let sign = oc.inward;
             let ci = |q: usize| f64::from(C[q][axis]);
@@ -675,6 +1011,22 @@ impl Flow {
     }
 }
 
+/// What a scalar population of a cell next to a wall or open face becomes (`None`: streamed as
+/// usual): a periodic image pulls, an adiabatic wall (or a slip face) bounces back, a held
+/// temperature anti-bounces back, a flux adds its heat, an open face is set afterwards.
+fn thermal_rule(g: &[f64], s: usize, groups: &[ThermalCondition], cell: usize, q: usize, rule: Rule) -> Option<f64> {
+    let back = g.get(opposite(q) * s + cell).copied()?;
+    match rule {
+        Rule::Pull { from } => g.get(q * s + from as usize).copied(),
+        Rule::Bounce { group, .. } | Rule::Reflect { group, .. } => Some(match groups.get(group as usize).copied().unwrap_or_default() {
+            ThermalCondition::Adiabatic {} => back,
+            ThermalCondition::Temperature { value } => -back + 2.0 * WT.get(q).copied().unwrap_or(0.0) * value,
+            ThermalCondition::Flux { value } => back + value,
+        }),
+        Rule::Open => Some(back),
+    }
+}
+
 impl crate::backend::Backend for Flow {
     fn name(&self) -> String {
         "cpu".into()
@@ -701,7 +1053,12 @@ impl crate::backend::Backend for Flow {
                 *v = u;
             }
         }
-        Ok(crate::backend::Fields { density, velocity })
+        let temperature = self.thermal.as_ref().map(|t| t.temperature.clone()).unwrap_or_default();
+        Ok(crate::backend::Fields { density, velocity, temperature })
+    }
+
+    fn heat(&mut self) -> Result<Vec<f64>> {
+        Ok(Flow::heat(self).to_vec())
     }
 
     fn mass(&mut self) -> Result<f64> {

crates/lbm/src/lib.rs

@@ -56,9 +56,10 @@ pub mod solve;
 pub mod voxel;
 
 pub use backend::{Backend, Choice, Fields, register_backend};
-pub use flow::{Flow, Group, Link, OpenCell, Rule, Target};
+pub use flow::{Flow, Group, Link, OpenCell, Rule, Target, ThermalState};
+pub use setup::ThermalSetup;
 pub use setup::{Obstacle, Setup, fluid_step};
-pub use solve::{FluidResult, FluidStats, ForceHistory, Frame, run, solve_fluid, solve_fluid_on};
+pub use solve::{FluidResult, FluidStats, ForceHistory, Frame, HeatHistory, run, solve_fluid, solve_fluid_on};
 pub use voxel::Grid;
 
 #[derive(Debug, Clone, PartialEq, thiserror::Error)]

crates/lbm/src/setup.rs

@@ -26,7 +26,7 @@
 //! density stops the run with [`crate::Error::Unstable`].
 
 use fieldcraft_mesh::surface::BoundarySurface;
-use fieldcraft_model::{BoxFace, FaceCondition, FluidDomain, Model, Procedure, SurfaceCondition};
+use fieldcraft_model::{BoxFace, Directional, FaceCondition, FluidDomain, Model, Procedure, SurfaceCondition, ThermalCondition};
 
 use crate::voxel::Grid;
 use crate::{Error, Limits, Result};
@@ -66,6 +66,25 @@ pub struct ObstacleFaces {
     pub names: Vec<String>,
     /// Every face tag's wall velocity in lattice units (zero for walls at rest).
     pub velocities: Vec<[f64; 3]>,
+    /// Every face tag's thermal condition (fluxes in lattice units); adiabatic without heat.
+    pub thermal: Vec<ThermalCondition>,
+}
+
+/// Heat in the fluid, in lattice units: a D3Q7 scalar lattice for the temperature
+/// (`flow::ThermalState`).
+#[derive(Debug, Clone, PartialEq)]
+pub struct ThermalSetup {
+    /// The scalar lattice's relaxation time: its diffusivity is `(τ_T − ½) / 4`.
+    pub tau: f64,
+    /// The fluid's thermal diffusivity k / (ρ c), model units.
+    pub diffusivity: f64,
+    pub initial_temperature: f64,
+    /// `−β g` in lattice units: the buoyant acceleration per temperature above the initial one.
+    pub buoyancy: [f64; 3],
+    /// The box faces' conditions (fluxes in lattice units: temperature per step through a face).
+    pub faces: [ThermalCondition; 6],
+    /// Power per lattice heat flow (a temperature per step through one cell face): ρ c Δx³ / Δt.
+    pub heat_scale: f64,
 }
 
 /// A fluid step in lattice units.
@@ -93,6 +112,8 @@ pub struct Setup {
     pub initial_velocity: [f64; 3],
     /// The domain's obstacles, in its order.
     pub obstacles: Vec<ObstacleFaces>,
+    /// Heat in the fluid, when the domain has thermal settings.
+    pub thermal: Option<ThermalSetup>,
     pub warnings: Vec<String>,
 }
 
@@ -108,6 +129,21 @@ pub fn fluid_step<'m>(model: &'m Model, step: Option<&str>) -> Result<(&'m field
     Ok((s, domain))
 }
 
+/// A thermal condition in lattice units: a flux (power per area) becomes the temperature it adds
+/// to the cell behind the face per step.
+fn lattice_thermal(c: &ThermalCondition, flux_scale: f64) -> Result<ThermalCondition> {
+    Ok(match *c {
+        ThermalCondition::Flux { value } => {
+            let lattice = value * flux_scale;
+            if !lattice.is_finite() {
+                return Err(Error::Invalid(format!("a heat flux of {value} is too large for this lattice")));
+            }
+            ThermalCondition::Flux { value: lattice }
+        }
+        other => other,
+    })
+}
+
 impl Setup {
     /// Convert a fluid step, checking it against the stability limits and `limits`.
     pub fn new(model: &Model, step: Option<&str>, obstacles: &[Obstacle], limits: &Limits) -> Result<Setup> {
@@ -243,9 +279,68 @@ impl Setup {
             )));
         }
 
+        // Heat: the scalar lattice's relaxation time from the diffusivity, the conditions in
+        // lattice units (a flux is a temperature per step through a cell face).
+        let mut warnings = Vec::new();
+        let thermal = match &domain.thermal {
+            None => None,
+            Some(t) => {
+                let props = material.and_then(|m| m.thermal.as_ref());
+                let (Some(Directional::Isotropic(k)), Some(c)) = (props.and_then(|p| p.conductivity), props.and_then(|p| p.specific_heat)) else {
+                    return Err(Error::Invalid(format!(
+                        "material `{}` needs an isotropic conductivity and a specific heat for heat in the fluid",
+                        domain.material
+                    )));
+                };
+                let rho_c = density * c;
+                let diffusivity = k / rho_c;
+                if !(diffusivity > 0.0 && diffusivity.is_finite()) {
+                    return Err(Error::Invalid(format!("the thermal diffusivity k / (ρ c) of material `{}` is {diffusivity}", domain.material)));
+                }
+                let alpha = diffusivity * dt / (h * h);
+                let tau = 4.0 * alpha + 0.5;
+                if !(tau >= MIN_TAU) {
+                    return Err(Error::Invalid(format!(
+                        "the thermal relaxation time τ_T = {tau:.5} is below {MIN_TAU}: the grid cannot resolve diffusivity {diffusivity} at this time step (cell Péclet number {:.1}); refine the spacing or raise the time step",
+                        if alpha > 0.0 { lattice_velocity / alpha } else { f64::INFINITY }
+                    )));
+                }
+                if tau > 2.0 {
+                    warnings.push(format!("the thermal relaxation time τ_T = {tau:.3} is large: heat diffuses fast on this grid and time step; a shorter time step is more accurate"));
+                }
+                let flux_scale = dt / (rho_c * h);
+                let mut faces = [ThermalCondition::Adiabatic {}; 6];
+                for face in BoxFace::ALL {
+                    let c = t.faces.get(face);
+                    if matches!(c, ThermalCondition::Flux { .. })
+                        && matches!(
+                            domain.faces.get(face),
+                            FaceCondition::Velocity { .. } | FaceCondition::Pressure { .. } | FaceCondition::Periodic {}
+                        )
+                    {
+                        return Err(Error::Invalid(format!(
+                            "box face {}: a heat flux goes on a wall, not on an inlet, outlet or periodic face",
+                            face.name()
+                        )));
+                    }
+                    if let Some(slot) = faces.get_mut(face.index()) {
+                        *slot = lattice_thermal(c, flux_scale)?;
+                    }
+                }
+                let beta = t.expansion.unwrap_or(0.0);
+                Some(ThermalSetup {
+                    tau,
+                    diffusivity,
+                    initial_temperature: t.initial_temperature,
+                    buoyancy: t.gravity.map(|g| -beta * g * dt * dt / h),
+                    faces,
+                    heat_scale: rho_c * h * h * h / dt,
+                })
+            }
+        };
+
         // Obstacles, their face names and their faces' wall velocities.
         let mut domain_obstacles = Vec::with_capacity(domain.obstacles.len());
-        let mut warnings = Vec::new();
         for body in &domain.obstacles {
             let (index, o) = obstacles
                 .iter()
@@ -255,6 +350,25 @@ impl Setup {
             let tags = o.surface.tri_face().iter().copied().max().map_or(1, |m| m as usize + 1).max(o.faces.len());
             let names: Vec<String> = (0..tags).map(|t| o.faces.get(t).cloned().unwrap_or_else(|| format!("face{t}"))).collect();
             let mut velocities = vec![[0.0; 3]; tags];
+            let mut thermal_faces = vec![ThermalCondition::Adiabatic {}; tags];
+            if let (Some(t), Some(ts)) = (&domain.thermal, &thermal) {
+                // power per area → temperature per step through a cell face, as for the box faces
+                let flux_scale = h * h / ts.heat_scale;
+                for sb in t.surfaces.iter().filter(|sb| &sb.face.body == body) {
+                    let mut hit = false;
+                    for (tag, name) in names.iter().enumerate() {
+                        if sb.face.selects(&fieldcraft_model::GeometryRef::face(body.as_str(), name.as_str())?) {
+                            hit = true;
+                            if let Some(slot) = thermal_faces.get_mut(tag) {
+                                *slot = lattice_thermal(&sb.condition, flux_scale)?;
+                            }
+                        }
+                    }
+                    if !hit {
+                        return Err(Error::NotFound(format!("the {} is not a face of the obstacle's surface", sb.face)));
+                    }
+                }
+            }
             for sb in domain.surfaces.iter().filter(|sb| &sb.face.body == body) {
                 let mut hit = false;
                 for (tag, name) in names.iter().enumerate() {
@@ -269,7 +383,7 @@ impl Setup {
                     return Err(Error::NotFound(format!("the {} is not a face of the obstacle's surface", sb.face)));
                 }
             }
-            domain_obstacles.push(ObstacleFaces { index, names, velocities });
+            domain_obstacles.push(ObstacleFaces { index, names, velocities, thermal: thermal_faces });
         }
         if mach > 0.1 + 1e-12 {
             warnings.push(format!(
@@ -298,6 +412,7 @@ impl Setup {
             acceleration,
             initial_velocity: domain.initial_velocity.map(|c| c * dt / h),
             obstacles: domain_obstacles,
+            thermal,
             warnings,
             domain,
         })

crates/lbm/src/solve.rs

@@ -23,10 +23,21 @@ pub struct Frame {
     pub velocity: Vec<[f32; 3]>,
     /// Gauge pressure.
     pub pressure: Vec<f32>,
+    /// Temperature, when the flow carries heat (empty otherwise).
+    pub temperature: Vec<f32>,
     /// The total fluid mass (density times volume).
     pub mass: f64,
 }
 
+/// The heat the fluid gives one wall over time (power; positive into the wall).
+#[derive(Debug, Clone, PartialEq)]
+pub struct HeatHistory {
+    pub name: String,
+    pub body: Option<String>,
+    /// One value per [`FluidResult::history_times`].
+    pub power: Vec<f64>,
+}
+
 /// The force the fluid exerts on one wall over time.
 #[derive(Debug, Clone, PartialEq)]
 pub struct ForceHistory {
@@ -57,6 +68,8 @@ pub struct FluidStats {
     pub cell_updates: u64,
     /// What stepped the flow: `cpu`, or `gpu: <adapter>`.
     pub backend: String,
+    /// The scalar lattice's relaxation time, with heat.
+    pub tau_thermal: Option<f64>,
 }
 
 /// The result of a fluid step.
@@ -71,17 +84,19 @@ pub struct FluidResult {
     pub history_times: Vec<f64>,
     /// Every box wall and every obstacle face.
     pub forces: Vec<ForceHistory>,
+    /// The heat every wall receives, with heat (empty otherwise), the same walls and times.
+    pub heat: Vec<HeatHistory>,
     pub stats: FluidStats,
     pub warnings: Vec<String>,
 }
 
 impl FluidResult {
-    /// The cell containing a point and the frame's velocity and pressure there (`None` outside
+    /// The cell containing a point and the frame's velocity, pressure and temperature (`None`
-    /// the grid; zeros in a solid cell).
+    /// without heat) there (`None` outside the grid; zeros in a solid cell).
-    pub fn probe(&self, frame: usize, p: [f64; 3]) -> Option<(usize, [f32; 3], f32)> {
+    pub fn probe(&self, frame: usize, p: [f64; 3]) -> Option<(usize, [f32; 3], f32, Option<f32>)> {
         let cell = self.grid.cell_at(p)?;
         let f = self.frames.get(frame)?;
-        Some((cell, *f.velocity.get(cell)?, *f.pressure.get(cell)?))
+        Some((cell, *f.velocity.get(cell)?, *f.pressure.get(cell)?, f.temperature.get(cell).copied()))
     }
 }
 
@@ -106,7 +121,11 @@ fn frame(setup: &Setup, flow: &mut dyn Backend, step: u64) -> Result<Frame> {
     velocity.try_reserve_exact(n).map_err(|_| Error::Limit("not enough memory for an output frame".into()))?;
     pressure.try_reserve_exact(n).map_err(|_| Error::Limit("not enough memory for an output frame".into()))?;
     let vs = setup.velocity_scale();
-    let Fields { density, velocity: u } = flow.fields()?;
+    let Fields { density, velocity: u, temperature: t } = flow.fields()?;
+    let mut temperature = Vec::new();
+    if !t.is_empty() {
+        temperature.try_reserve_exact(n).map_err(|_| Error::Limit("not enough memory for an output frame".into()))?;
+    }
     for cell in 0..n {
         if flow.is_fluid(cell) {
             velocity.push(u.get(cell).copied().unwrap_or([0.0; 3]).map(|c| (c * vs) as f32));
@@ -115,9 +134,12 @@ fn frame(setup: &Setup, flow: &mut dyn Backend, step: u64) -> Result<Frame> {
             velocity.push([0.0; 3]);
             pressure.push(0.0);
         }
+        if !t.is_empty() {
+            temperature.push(t.get(cell).copied().unwrap_or(0.0) as f32);
+        }
     }
     let mass = flow.mass()? * setup.density * setup.grid.spacing.powi(3);
-    Ok(Frame { time: step as f64 * setup.time_step, step, velocity, pressure, mass })
+    Ok(Frame { time: step as f64 * setup.time_step, step, velocity, pressure, temperature, mass })
 }
 
 /// Run a prepared flow (fresh from [`Flow::new`], possibly with its own initial state, on any
@@ -127,7 +149,7 @@ pub fn run(setup: &Setup, flow: &mut dyn Backend, limits: &Limits) -> Result<Flu
     let n = setup.grid.len();
     let total = setup.steps;
     let frames = u64::from(setup.frames).min(total).max(1);
-    let frame_bytes = (frames + 1) as f64 * n as f64 * 16.0;
+    let frame_bytes = (frames + 1) as f64 * n as f64 * if setup.thermal.is_some() { 20.0 } else { 16.0 };
     if !(frame_bytes <= limits.max_frame_bytes as f64) {
         return Err(Error::Limit(format!(
             "{} frames of {n} cells need {:.0} MB, above the cap of {:.0} MB; write fewer frames or raise `max_frame_bytes`",
@@ -138,10 +160,12 @@ pub fn run(setup: &Setup, flow: &mut dyn Backend, limits: &Limits) -> Result<Flu
     }
     let every = total.div_ceil(MAX_HISTORY).max(1);
     let fs = setup.force_scale();
+    let hs = setup.thermal.as_ref().map_or(0.0, |t| t.heat_scale);
 
     let mut out_frames = vec![frame(setup, flow, 0)?];
     let mut history_times = Vec::new();
     let mut forces: Vec<Vec<[f64; 3]>> = vec![Vec::new(); flow.groups().len()];
+    let mut heats: Vec<Vec<f64>> = vec![Vec::new(); if setup.thermal.is_some() { flow.groups().len() } else { 0 }];
     let mut next_frame = 1u64;
     // Frame j (1..=frames) at step round(j · total / frames).
     let frame_step = |j: u64| (j * total + frames / 2) / frames;
@@ -165,6 +189,11 @@ pub fn run(setup: &Setup, flow: &mut dyn Backend, limits: &Limits) -> Result<Flu
             for (h, f) in forces.iter_mut().zip(flow.forces()?) {
                 h.push(f.map(|c| c * fs));
             }
+            if !heats.is_empty() {
+                for (h, q) in heats.iter_mut().zip(flow.heat()?) {
+                    h.push(q * hs);
+                }
+            }
         }
         while next_frame <= frames && k == frame_step(next_frame) {
             out_frames.push(frame(setup, flow, k)?);
@@ -198,6 +227,7 @@ pub fn run(setup: &Setup, flow: &mut dyn Backend, limits: &Limits) -> Result<Flu
         frames: out_frames,
         history_times,
         forces: flow.groups().iter().zip(forces).map(|(g, force)| ForceHistory { name: g.name.clone(), body: g.body.clone(), force }).collect(),
+        heat: flow.groups().iter().zip(heats).map(|(g, power)| HeatHistory { name: g.name.clone(), body: g.body.clone(), power }).collect(),
         stats: FluidStats {
             cells: n,
             fluid_cells: flow.fluid_cells(),
@@ -212,6 +242,7 @@ pub fn run(setup: &Setup, flow: &mut dyn Backend, limits: &Limits) -> Result<Flu
             mach: setup.mach,
             cell_updates: (flow.fluid_cells() as u64).saturating_mul(total),
             backend: flow.name(),
+            tau_thermal: setup.thermal.as_ref().map(|t| t.tau),
         },
         warnings,
     })

crates/lbm/tests/thermal.rs

@@ -0,0 +1,297 @@
+//! Heat in the fluid: the temperature lattice against closed-form solutions and a published
+//! benchmark (de Vahl Davis, heavy).
+
+mod common;
+
+use common::*;
+use fieldcraft_lbm::backend::Backend;
+use fieldcraft_lbm::{Limits, run};
+use fieldcraft_model::{BoxFace, Directional, FaceCondition, Fluid, FluidThermal, Material, Model, Procedure, Step, Thermal, ThermalCondition};
+
+const PERIODIC: FaceCondition = FaceCondition::Periodic {};
+
+/// A fluid with heat: density ρ, kinematic viscosity ν, thermal diffusivity α (conductivity
+/// ρ α, unit specific heat), so k / (ρ c) = α.
+fn hot_model(density: f64, nu: f64, alpha: f64, domain: fieldcraft_model::FluidDomain, procedure: Procedure) -> Model {
+    let mut m = Model::default();
+    m.add_material(
+        "Fluid",
+        Material {
+            density: Some(density),
+            fluid: Some(Fluid::Newtonian { viscosity: nu * density }),
+            thermal: Some(Thermal { conductivity: Some(Directional::Isotropic(density * alpha)), specific_heat: Some(1.0), ..Default::default() }),
+            ..Material::default()
+        },
+    )
+    .unwrap();
+    m.add_fluid_domain("Flow", domain).unwrap();
+    m.add_step(Step::new("Run", procedure), None).unwrap();
+    m
+}
+
+/// Plates at y = 0 and y = 1, periodic in x and z, no flow: held at 0 and 1, the steady
+/// temperature is linear, and anti-bounce-back puts the walls exactly halfway between cells.
+#[test]
+fn conduction_between_plates_is_linear() {
+    for n in [8usize, 16] {
+        let (h, alpha) = (1.0, 0.05);
+        let dx = h / n as f64;
+        let mut d = domain([0.0; 3], [dx, h, dx], dx);
+        d.faces.x_min = PERIODIC;
+        d.faces.x_max = PERIODIC;
+        d.faces.z_min = PERIODIC;
+        d.faces.z_max = PERIODIC;
+        let mut t = FluidThermal::default();
+        t.faces.y_min = ThermalCondition::Temperature { value: 0.0 };
+        t.faces.y_max = ThermalCondition::Temperature { value: 1.0 };
+        d.thermal = Some(t);
+        // τ_T = 1: Δt = Δx² / (8 α)
+        let dt = dx * dx / (8.0 * alpha);
+        let m = hot_model(1.0, 0.1, alpha, d, step(5.0 * h * h / alpha, dt, None));
+        let (setup, mut flow) = prepare(&m, &[]);
+        assert!((setup.thermal.as_ref().unwrap().tau - 1.0).abs() < 1e-12);
+        run_all(&setup, &mut flow);
+        let mut err: f64 = 0.0;
+        for j in 0..n {
+            let y = (j as f64 + 0.5) * dx;
+            let t = flow.temperature(j).unwrap();
+            err = err.max((t - y / h).abs());
+        }
+        println!("conduction between plates, {n} cells: max error {err:.2e}");
+        assert!(err < 1e-9, "{n}: {err}");
+    }
+}
+
+/// A flux q through the cold plate and a held temperature on the other: the steady drop over
+/// the gap is q H / k, and the heat the hot plate receives is the flux times its area.
+#[test]
+fn a_heat_flux_sets_the_gradient_and_the_wall_heat_balances() {
+    let (n, h, alpha, rho) = (16usize, 1.0, 0.05, 2.0);
+    let dx = h / n as f64;
+    let mut d = domain([0.0; 3], [dx, h, dx], dx);
+    d.faces.x_min = PERIODIC;
+    d.faces.x_max = PERIODIC;
+    d.faces.z_min = PERIODIC;
+    d.faces.z_max = PERIODIC;
+    let q = 0.3;
+    let mut t = FluidThermal::default();
+    t.faces.y_min = ThermalCondition::Flux { value: q };
+    t.faces.y_max = ThermalCondition::Temperature { value: 1.0 };
+    d.thermal = Some(t);
+    let dt = dx * dx / (8.0 * alpha);
+    let m = hot_model(rho, 0.1, alpha, d, step(8.0 * h * h / alpha, dt, None));
+    let setup = fieldcraft_lbm::Setup::new(&m, None, &[], &Limits::default()).unwrap();
+    let mut flow = fieldcraft_lbm::Flow::new(&setup, &[]).unwrap();
+    let r = run(&setup, &mut flow, &Limits::default()).unwrap();
+    let k = rho * alpha;
+    // the temperature halfway below the first cell (the wall) against the hot wall at 1
+    let t0 = f64::from(r.frames.last().unwrap().temperature[0]) + 0.5 * dx * q / k;
+    let drop = t0 - 1.0;
+    let expected = q * h / k;
+    println!("flux plate: drop {drop:.6} vs q H / k = {expected:.6}");
+    assert!((drop - expected).abs() < 1e-6 * expected, "{drop} vs {expected}");
+    let hot = r.heat.iter().find(|w| w.name == "y_max").unwrap();
+    let last = *hot.power.last().unwrap();
+    let area = dx * dx;
+    println!("heat into the hot wall: {last:.6e} vs q A = {:.6e}", q * area);
+    assert!((last - q * area).abs() < 1e-6 * q * area, "{last} vs {}", q * area);
+    let cold = r.heat.iter().find(|w| w.name == "y_min").unwrap();
+    assert!((cold.power.last().unwrap() + q * area).abs() < 1e-6 * q * area, "the flux face gives q A");
+}
+
+/// A channel with a velocity inlet at 1, slip (adiabatic) walls and a pressure outlet: the
+/// plug flow carries the inlet temperature through the whole channel. The flow starts at
+/// its steady state (an impulsive start in a frictionless channel rings for a long time).
+#[test]
+fn the_flow_carries_the_inlet_temperature() {
+    let (h, l, umax, nu, alpha) = (1.0, 4.0, 0.1, 0.05, 0.01);
+    let dx = h / 8.0;
+    let mut d = domain([0.0; 3], [l, h, dx], dx);
+    d.faces.x_min = FaceCondition::Velocity { velocity: [umax, 0.0, 0.0], profile: fieldcraft_model::Profile::Uniform };
+    d.faces.x_max = FaceCondition::Pressure { pressure: 0.0 };
+    d.faces.y_min = FaceCondition::Slip {};
+    d.faces.y_max = FaceCondition::Slip {};
+    d.faces.z_min = PERIODIC;
+    d.faces.z_max = PERIODIC;
+    let mut t = FluidThermal::default();
+    t.faces.x_min = ThermalCondition::Temperature { value: 1.0 };
+    d.thermal = Some(t);
+    let lattice_u = 0.05;
+    let m = hot_model(1.0, nu, alpha, d, step(4.0 * l / umax, lattice_u * dx / umax, None));
+    let (setup, mut flow) = prepare(&m, &[]);
+    flow.set_equilibrium(|_| (1.0, [lattice_u, 0.0, 0.0]));
+    run_all(&setup, &mut flow);
+    let grid = *flow.grid();
+    let (mut lo, mut hi) = (f64::INFINITY, f64::NEG_INFINITY);
+    for c in 0..grid.len() {
+        let t = flow.temperature(c).unwrap();
+        lo = lo.min(t);
+        hi = hi.max(t);
+    }
+    println!("inlet temperature carried through: {lo:.6}..{hi:.6}");
+    assert!(lo > 1.0 - 1e-3 && hi < 1.0 + 1e-3, "{lo}..{hi}");
+    let f = flow.fields().unwrap();
+    assert_eq!(f.temperature.len(), grid.len());
+}
+
+/// Natural convection in a square cavity (de Vahl Davis, Int. J. Numer. Meth. Fluids 3 (1983)
+/// 249–264): hot left wall at 1, cold right wall at 0, adiabatic top and bottom, Pr = 0.71.
+/// The mean Nusselt number on the hot wall at Ra = 10³ is 1.118 (10⁴: 2.243); `u_max` is the
+/// fastest flow in the cavity in units of α / H (3.65 and 19.6). Returns the Nusselt number
+/// and the steady power into the hot wall and into the cold wall (the mean over the last
+/// tenth of the run).
+fn cavity(ra: f64, n: usize, u_max: f64) -> (f64, f64, f64) {
+    let (h, pr) = (1.0, 0.71);
+    let dx = h / n as f64;
+    let alpha = 1.0;
+    let nu = pr * alpha;
+    // g β ΔT H³ / (ν α) = Ra, with ΔT = 1, H = 1
+    let g_beta = ra * nu * alpha / (h * h * h);
+    let mut d = domain([0.0; 3], [h, h, dx], dx);
+    d.faces.z_min = PERIODIC;
+    d.faces.z_max = PERIODIC;
+    let mut t = FluidThermal { initial_temperature: 0.5, gravity: [0.0, -g_beta, 0.0], expansion: Some(1.0), ..FluidThermal::default() };
+    t.faces.x_min = ThermalCondition::Temperature { value: 1.0 };
+    t.faces.x_max = ThermalCondition::Temperature { value: 0.0 };
+    d.thermal = Some(t);
+    // the fastest flow at lattice speed 0.05, and a lattice diffusivity of at most 1/4 (τ_T ≤ 1.5)
+    let dt = (0.05 * dx / (u_max * alpha / h)).min(dx * dx / (4.0 * alpha));
+    let m = hot_model(1.0, nu, alpha, d, step(4.0 * h * h / alpha, dt, None));
+    let setup = fieldcraft_lbm::Setup::new(&m, None, &[], &Limits::default()).unwrap();
+    let mut flow = fieldcraft_lbm::Flow::new(&setup, &[]).unwrap();
+    let r = run(&setup, &mut flow, &Limits::default()).unwrap();
+    // Nu = (heat into the fluid through the hot wall) H / (k ΔT A): the hot wall's history
+    // holds the power the wall receives (negative when it heats the fluid).
+    let steady = |name: &str| {
+        let w = r.heat.iter().find(|w| w.name == name).unwrap();
+        let m = w.power.len();
+        let tail = &w.power[m - m / 10..];
+        tail.iter().sum::<f64>() / tail.len() as f64
+    };
+    let (hot, cold) = (steady("x_min"), steady("x_max"));
+    let k = alpha;
+    let area = h * dx;
+    let nu_mean = -hot * h / (k * 1.0 * area);
+    println!("cavity Ra = {ra:.0e}, {n} cells: Nu = {nu_mean:.4}, hot wall {hot:.6e}, cold wall {cold:.6e}; {:?}", r.warnings);
+    (nu_mean, hot, cold)
+}
+
+/// At steady state the heat the hot wall gives the fluid is the heat the cold wall takes.
+fn heat_balances(hot: f64, cold: f64, tolerance: f64) {
+    assert!(hot < 0.0 && cold > 0.0, "the hot wall heats the fluid ({hot}) and the cold wall cools it ({cold})");
+    assert!((hot + cold).abs() < tolerance * cold, "heat in {hot:e} vs heat out {cold:e}");
+}
+
+/// The published cavity on a coarse grid, every run: Ra = 10³ on 16 cells is within 2 % of
+/// de Vahl Davis (the heavy test refines it), and the heat balances between the walls.
+#[test]
+fn natural_convection_cavity_is_close_on_a_coarse_grid() {
+    let (nu3, hot, cold) = cavity(1e3, 16, 3.65);
+    assert!((nu3 / 1.118 - 1.0).abs() < 0.02, "Ra 1e3: {nu3}");
+    heat_balances(hot, cold, 1e-3);
+}
+
+#[test]
+#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"]
+fn natural_convection_cavity_matches_de_vahl_davis() {
+    let (nu3, hot, cold) = cavity(1e3, 32, 3.65);
+    assert!((nu3 / 1.118 - 1.0).abs() < 0.01, "Ra 1e3: {nu3}");
+    heat_balances(hot, cold, 1e-3);
+    let (nu4, hot, cold) = cavity(1e4, 48, 19.6);
+    assert!((nu4 / 2.243 - 1.0).abs() < 0.01, "Ra 1e4: {nu4}");
+    heat_balances(hot, cold, 1e-3);
+}
+
+/// A periodic box with a uniform flow and a blob of heat: the total temperature is conserved
+/// to round-off and the blob moves with the flow.
+#[test]
+fn a_periodic_flow_conserves_and_carries_heat() {
+    let (l, u0, nu, alpha) = (4.0, 0.1, 0.05, 0.01);
+    let dx = l / 32.0;
+    let mut d = domain([0.0; 3], [l, dx, dx], dx);
+    for f in BoxFace::ALL {
+        *d.faces.get_mut(f) = PERIODIC;
+    }
+    d.initial_velocity = [u0, 0.0, 0.0];
+    d.thermal = Some(FluidThermal::default());
+    let m = hot_model(1.0, nu, alpha, d, step(l / u0, 0.05 * dx / u0, None));
+    let (setup, mut flow) = prepare(&m, &[]);
+    let grid = *flow.grid();
+    flow.set_temperature(|c| if grid.coords(c)[0] < 8 { 1.0 } else { 0.0 });
+    let total = |f: &fieldcraft_lbm::Flow| (0..grid.len()).map(|c| f.temperature(c).unwrap()).sum::<f64>();
+    let t0 = total(&flow);
+    run_all(&setup, &mut flow);
+    let t1 = total(&flow);
+    let row: Vec<String> = (0..32).map(|i| format!("{:.3}", flow.temperature(grid.index([i, 0, 0])).unwrap())).collect();
+    println!("periodic blob after one lap: total {t0:.6} -> {t1:.6}; {}", row.join(" "));
+    assert!((t1 - t0).abs() < 1e-10 * t0, "{t0} vs {t1}");
+}
+
+/// A sine of temperature in a periodic box drifting with a uniform flow: it keeps its shape,
+/// moves with the flow and decays as e^(−α k² t). At rest the decay converges to the
+/// diffusivity at second order. With the flow the D3Q7 lattice's linear equilibrium diffuses
+/// less along the flow by (τ_T − ½) u² in lattice units, 4 u² of the diffusivity (Chopard,
+/// Falcone and Latt, Eur. Phys. J. Special Topics 171 (2009) 245–249): the measured decay
+/// shows that deficit, which bounds the lattice speed a heat problem should run at.
+#[test]
+fn a_drifting_sine_decays_at_the_diffusivity() {
+    let (l, nu, alpha, amplitude) = (4.0, 0.05, 0.01, 0.25);
+    let k = 2.0 * std::f64::consts::PI / l;
+    let decay_time = 1.0 / (alpha * k * k);
+    let mut at_rest = Vec::new();
+    for (n, lattice_u) in [(32usize, 0.0), (64, 0.0), (64, 0.1), (64, 0.15)] {
+        let dx = l / n as f64;
+        // At rest τ_T = 1; with the flow at `lattice_u` the time step follows the speed.
+        let u0 = if lattice_u > 0.0 { 0.1 } else { 0.0 };
+        let dt = if lattice_u > 0.0 { lattice_u * dx / u0 } else { dx * dx / (8.0 * alpha) };
+        let mut d = domain([0.0; 3], [l, dx, dx], dx);
+        for f in BoxFace::ALL {
+            *d.faces.get_mut(f) = PERIODIC;
+        }
+        d.initial_velocity = [u0, 0.0, 0.0];
+        d.thermal = Some(FluidThermal { initial_temperature: 0.5, ..FluidThermal::default() });
+        let m = hot_model(1.0, nu, alpha, d, step(decay_time, dt, None));
+        let (setup, mut flow) = prepare(&m, &[]);
+        let tau = setup.thermal.as_ref().unwrap().tau;
+        let grid = *flow.grid();
+        let x = |c: usize| (grid.coords(c)[0] as f64 + 0.5) * dx;
+        flow.set_temperature(|c| 0.5 + amplitude * (k * x(c)).sin());
+        run_all(&setup, &mut flow);
+        let t = setup.steps as f64 * setup.time_step;
+        // The sine's amplitude and phase from the Fourier coefficient at k, then the shape.
+        let (mut a, mut b) = (0.0, 0.0);
+        for i in 0..n {
+            let c = grid.index([i, 0, 0]);
+            let temp = flow.temperature(c).unwrap();
+            a += temp * (k * x(c)).sin();
+            b += temp * (k * x(c)).cos();
+        }
+        let measured = (a * a + b * b).sqrt() * 2.0 / n as f64;
+        let drift = -b.atan2(a) / k;
+        let mut shape: f64 = 0.0;
+        for i in 0..n {
+            let c = grid.index([i, 0, 0]);
+            shape = shape.max((flow.temperature(c).unwrap() - 0.5 - measured * (k * (x(c) - drift)).sin()).abs());
+        }
+        let alpha_measured = -(measured / amplitude).ln() / (k * k * t);
+        let error = alpha_measured / alpha - 1.0;
+        let lag = (drift - u0 * t).rem_euclid(l);
+        let lag = lag.min(l - lag);
+        println!(
+            "drifting sine, {n} cells, lattice speed {lattice_u}: τ_T = {tau:.4}, α measured {alpha_measured:.5} vs {alpha} ({error:+.2e}), lag {lag:.2e}, shape error {shape:.2e}"
+        );
+        // A Fourier mode of a linear scheme on a periodic grid stays one, to round-off.
+        assert!(shape < 1e-9 * amplitude, "{n} cells at {lattice_u}: the sine lost its shape by {shape}");
+        assert!(lag < dx / 5.0, "{n} cells at {lattice_u}: drift {drift} vs {}", u0 * t);
+        if lattice_u > 0.0 {
+            let deficit = -4.0 * lattice_u * lattice_u;
+            assert!((error - deficit).abs() < 0.01, "{n} cells at {lattice_u}: α error {error} vs the deficit {deficit}");
+        } else {
+            assert!(error.abs() < 1e-3, "{n} cells at rest: α error {error}");
+            at_rest.push(error.abs());
+        }
+    }
+    let order = orders(&at_rest, 2.0);
+    println!("decay at rest converges at order {order:?}");
+    assert!(order.iter().all(|o| *o > 1.8), "{order:?}");
+}