diff --git a/crates/lbm-gpu/src/native.rs b/crates/lbm-gpu/src/native.rs index 92bf3d8..c94b3be 100644 --- a/crates/lbm-gpu/src/native.rs +++ b/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 { + 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> { + Ok(Vec::new()) } fn mass(&mut self) -> Result { diff --git a/crates/lbm/src/backend.rs b/crates/lbm/src/backend.rs index c9dafc9..e80d930 100644 --- a/crates/lbm/src/backend.rs +++ b/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, pub velocity: Vec<[f64; 3]>, + pub temperature: Vec, } /// A stepper of one prepared flow. @@ -35,6 +37,9 @@ pub trait Backend { fn mass(&mut self) -> Result; /// The momentum each force group received in the last step (lattice units: force per step). fn forces(&mut self) -> Result>; + /// 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>; fn groups(&self) -> &[Group]; fn fluid_cells(&self) -> usize; fn is_fluid(&self, cell: usize) -> bool; diff --git a/crates/lbm/src/flow.rs b/crates/lbm/src/flow.rs index 83061b4..dd84b2a 100644 --- a/crates/lbm/src/flow.rs +++ b/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, } +/// 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, + g2: Vec, + /// The temperature of every cell after the last step (the initial one in solid cells). + temperature: Vec, + 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, + /// Per open cell: the temperature it holds (an inlet), or `None` to take its neighbour's + /// (an outlet). + open: Vec>, + /// The heat each group received in the last step (lattice units: temperature times cell + /// faces per step). + heat: Vec, +} + /// 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, + /// The temperature lattice, with thermal settings. + thermal: Option, + /// 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 = 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> = 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 { + 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 forced = a != [0.0; 3]; + let constant = self.acceleration; + 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(); - parallel::for_each_mut(&mut work, |b, (f, max_u2)| { + // with heat, each block's velocities go to its slice of cell_u (empty otherwise) + 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 { + 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> { + Ok(Flow::heat(self).to_vec()) } fn mass(&mut self) -> Result { diff --git a/crates/lbm/src/lib.rs b/crates/lbm/src/lib.rs index 7721b26..6444c52 100644 --- a/crates/lbm/src/lib.rs +++ b/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)] diff --git a/crates/lbm/src/setup.rs b/crates/lbm/src/setup.rs index 9e9cd6e..6a316a7 100644 --- a/crates/lbm/src/setup.rs +++ b/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, /// 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, +} + +/// 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, + /// Heat in the fluid, when the domain has thermal settings. + pub thermal: Option, pub warnings: Vec, } @@ -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 { + 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 { @@ -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 = (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, }) diff --git a/crates/lbm/src/solve.rs b/crates/lbm/src/solve.rs index f47dbc2..0ff1d8c 100644 --- a/crates/lbm/src/solve.rs +++ b/crates/lbm/src/solve.rs @@ -23,10 +23,21 @@ pub struct Frame { pub velocity: Vec<[f32; 3]>, /// Gauge pressure. pub pressure: Vec, + /// Temperature, when the flow carries heat (empty otherwise). + pub temperature: Vec, /// 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, + /// One value per [`FluidResult::history_times`]. + pub power: Vec, +} + /// 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: `. pub backend: String, + /// The scalar lattice's relaxation time, with heat. + pub tau_thermal: Option, } /// The result of a fluid step. @@ -71,17 +84,19 @@ pub struct FluidResult { pub history_times: Vec, /// Every box wall and every obstacle face. pub forces: Vec, + /// The heat every wall receives, with heat (empty otherwise), the same walls and times. + pub heat: Vec, pub stats: FluidStats, pub warnings: Vec, } impl FluidResult { - /// The cell containing a point and the frame's velocity and pressure there (`None` outside - /// the grid; zeros in a solid cell). - pub fn probe(&self, frame: usize, p: [f64; 3]) -> Option<(usize, [f32; 3], f32)> { + /// The cell containing a point and the frame's velocity, pressure and temperature (`None` + /// 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, Option)> { 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 { 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 { 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 Result> = vec![Vec::new(); flow.groups().len()]; + let mut heats: Vec> = 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 Result Result 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::() / 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::(); + let t0 = total(&flow); + run_all(&setup, &mut flow); + let t1 = total(&flow); + let row: Vec = (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:?}"); +}