diff --git a/crates/lbm-gpu/Cargo.toml b/crates/lbm-gpu/Cargo.toml index 8cdcc6c..e98bce2 100644 --- a/crates/lbm-gpu/Cargo.toml +++ b/crates/lbm-gpu/Cargo.toml @@ -10,11 +10,11 @@ workspace = true [dependencies] fieldcraft-lbm = { workspace = true } +fieldcraft-model = { workspace = true } fieldcraft-gpu = { workspace = true } wgpu = { workspace = true } pollster = { workspace = true } log = { workspace = true } [dev-dependencies] -fieldcraft-model = { workspace = true } fieldcraft-mesh = { workspace = true } diff --git a/crates/lbm-gpu/src/lbm.wgsl b/crates/lbm-gpu/src/lbm.wgsl index 2bda46c..d13a6a8 100644 --- a/crates/lbm-gpu/src/lbm.wgsl +++ b/crates/lbm-gpu/src/lbm.wgsl @@ -10,8 +10,18 @@ // the weight's): ρ = 1 + Σ h, momentum Σ h c (Σ w c = 0), and every rule below is written for h. // The CPU scheme (double precision, fixed sums) is the reference; see the crate's tests for the // tolerances. +// +// With heat, a D3Q7 scalar lattice (the rest direction and the six axes of D3Q19, weights ¼ and +// ⅛) rides on the flow in the buffers t and t2, stored as deviations from w_q T₀ (T₀ the initial +// temperature) so that f32 keeps the digits of the differences: θ = T − T₀ = Σ t. It collides +// (BGK) with the velocity the momentum collision left in `moments`, streams with the links' rules +// read through the group's thermal condition (anti-bounce-back for a held temperature, bounce- +// back for an adiabatic face, a flux added to the bounced population), takes its inlet +// temperatures and outlet neighbours, and gives each group the heat its links exchanged (in +// contrib.w and forces.w). The temperature feeds back as the buoyant acceleration −β g θ. const Q: u32 = 19u; +const QT: u32 = 7u; const WG: u32 = 256u; const KIND_PULL: u32 = 0u; const KIND_BOUNCE: u32 = 1u; @@ -33,6 +43,10 @@ struct Params { s_minus: f32, // The body acceleration (xyz) and whether there is one (w). acc: vec4, + // The buoyancy −β g in lattice units (xyz) and whether there is any (w). + buoy: vec4, + // The scalar lattice: its relaxation rate (x), T₀ (y), whether there is heat (z). + therm: vec4, }; @group(0) @binding(0) var params: Params; @@ -57,6 +71,14 @@ struct Params { @group(0) @binding(14) var forces: array>; // Per workgroup, the mass of its cells. @group(0) @binding(15) var partials: array; +// The scalar populations, A–B like f and g (t[q n + cell]; see the notes above). +@group(0) @binding(16) var t: array; +@group(0) @binding(17) var t2: array; +// Per force group: its thermal condition (x: 0 adiabatic, 1 temperature, 2 flux) and the value +// (y: the held temperature − T₀, or the flux as a temperature per step). +@group(0) @binding(18) var group_t: array>; +// Per open cell: whether it holds a temperature (x) and that temperature − T₀ (y). +@group(0) @binding(19) var open_t: array>; var red: array, 256>; // The largest squared speed of a workgroup's cells (collide): one atomic per workgroup, not @@ -99,6 +121,41 @@ fn w_of(q: u32) -> f32 { return 1.0 / 36.0; } +fn wt_of(q: u32) -> f32 { + if q == 0u { return 0.25; } + return 0.125; +} + +fn has_heat() -> bool { + return params.therm.z > 0.5; +} + +// The temperature of a cell above T₀: the sum of its scalar populations. +fn theta_of(cell: u32) -> f32 { + var s = 0.0; + for (var q = 0u; q < QT; q = q + 1u) { + s = s + t[q * params.n + cell]; + } + return s; +} + +// The body acceleration of a cell: the constant one plus the buoyant one. +fn acc_of(cell: u32) -> vec3 { + var a = params.acc.xyz; + if params.buoy.w > 0.5 { a = a + params.buoy.xyz * theta_of(cell); } + return a; +} + +fn forced() -> bool { + return params.acc.w > 0.5 || params.buoy.w > 0.5; +} + +// The scalar equilibrium w_q T (1 + 4 c_q·u), shifted by w_q T₀ like the populations. +fn eq_t(q: u32, theta: f32, u: vec3) -> f32 { + let cu = 4.0 * dot(vec3(c_of(q)), u); + return wt_of(q) * (theta * (1.0 + cu) + params.therm.y * cu); +} + fn cell_of(lid: u32, wg: vec3) -> u32 { return (wg.y * params.groups_x + wg.x) * WG + lid; } @@ -140,13 +197,14 @@ fn collide_cell(cell: u32, lid: u32) { atomicStore(&stats[1], 1u); return; } - let a = params.acc.xyz; + let a = acc_of(cell); let force = rho * a; let u = (j + 0.5 * force) / rho; let uu = dot(u, u); speed[lid] = uu; + if has_heat() { moments[cell] = vec4(u, rho); } let uf = dot(u, force); - let forced = params.acc.w > 0.5; + let forced = forced(); let wp = params.omega_plus; let wm = params.omega_minus; let sp = params.s_plus; @@ -242,7 +300,7 @@ fn open_boundaries(@builtin(local_invocation_index) lid: u32, @builtin(workgroup let cell = ou.x; let axis = ou.y; let sign = ot.w; - let a_half = 0.5 * params.acc.xyz; + let a_half = 0.5 * acc_of(cell); // The sums of f over the tangential (c_n = 0) and outgoing (c_n < 0) directions: the sums of // h plus the weights' (2/3 and 1/6 for D3Q19). var tangential = 2.0 / 3.0; @@ -307,6 +365,7 @@ fn link_forces(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) // cancels between a wall's links in exact arithmetic but is 10⁴ times the shear they carry) // is added on the host in double precision, so the sum keeps the deviations' digits. var v = vec3(0.0); + var heat = 0.0; if kind == KIND_BOUNCE { let out = f[opp(q) * n + cell]; let back = out + link_f[l]; @@ -315,7 +374,121 @@ fn link_forces(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) let out = f[q2 * n + lu.y]; v = vec3(c_of(q2)) * out - vec3(c_of(q)) * out; } - contrib[l] = vec4(v, 0.0); + // The heat the link gave its wall: what left the fluid minus what came back. + if has_heat() && q < QT && (kind == KIND_BOUNCE || kind == KIND_REFLECT) { + heat = t[opp(q) * n + cell] - thermal_rule(lu, q, kind, cell); + } + contrib[l] = vec4(v, heat); +} + +// ---- The scalar lattice ---- + +// BGK collision of the scalar lattice, in place in t, with the velocity of the momentum +// collision. +@compute @workgroup_size(256) +fn collide_thermal(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: vec3) { + let cell = cell_of(lid, wg); + if cell >= params.n { return; } + if !is_fluid(cell) { return; } + let n = params.n; + var tq: array; + var theta = 0.0; + for (var q = 0u; q < QT; q = q + 1u) { + let v = t[q * n + cell]; + tq[q] = v; + theta = theta + v; + } + let u = moments[cell].xyz; + let omega = params.therm.x; + for (var q = 0u; q < QT; q = q + 1u) { + t[q * n + cell] = tq[q] - omega * (tq[q] - eq_t(q, theta, u)); + } +} + +// Pull streaming of t into t2, first pass (as `stream`). +@compute @workgroup_size(256) +fn stream_thermal(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: vec3) { + let cell = cell_of(lid, wg); + if cell >= params.n { return; } + if !is_fluid(cell) { return; } + let n = params.n; + let nx = i32(params.nx); + let ny = i32(params.ny); + let nz = i32(params.nz); + let i = i32(cell % params.nx); + let jj = i32((cell / params.nx) % params.ny); + let k = i32(cell / (params.nx * params.ny)); + for (var q = 0u; q < QT; q = q + 1u) { + let c = c_of(q); + let si = (i - c.x + nx) % nx; + let sj = (jj - c.y + ny) % ny; + let sk = (k - c.z + nz) % nz; + let src = u32((sk * ny + sj) * nx + si); + t2[q * n + cell] = t[q * n + src]; + } +} + +// The scalar population a link's rule gives, from t: a pull as the flow's; at a wall or a slip +// face the group's thermal condition on the bounced population; an open placeholder. +fn thermal_rule(lu: vec4, q: u32, kind: u32, cell: u32) -> f32 { + let n = params.n; + if kind == KIND_PULL { return t[q * n + lu.y]; } + let back = t[opp(q) * n + cell]; + if kind == KIND_OPEN { return back; } + let cond = group_t[lu.z]; + let k = u32(cond.x); + if k == 1u { return -back + 2.0 * wt_of(q) * cond.y; } + if k == 2u { return back + cond.y; } + return back; +} + +// Second pass, one thread per link of the first seven directions. +@compute @workgroup_size(256) +fn stream_links_thermal(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: vec3) { + let l = cell_of(lid, wg); + if l >= params.links { return; } + let lu = link_u[l]; + let q = lu.x & 0xffu; + if q >= QT { return; } + let kind = (lu.x >> 8u) & 0xffu; + let cell = lu.w; + t2[q * params.n + cell] = thermal_rule(lu, q, kind, cell); +} + +// Inlets hold their temperature (the equilibrium at the cell's velocity); outlets take the +// streamed populations of the neighbour one cell inward. +@compute @workgroup_size(256) +fn open_thermal(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: vec3) { + let o = cell_of(lid, wg); + if o >= params.open { return; } + let n = params.n; + let ou = open_u[o]; + let cell = ou.x; + let held = open_t[o]; + if held.x > 0.5 { + let u = moments[cell].xyz; + for (var q = 0u; q < QT; q = q + 1u) { + t2[q * n + cell] = eq_t(q, held.y, u); + } + } else { + var step = 1; + if ou.y == 1u { step = i32(params.nx); } + if ou.y == 2u { step = i32(params.nx * params.ny); } + let src = u32(i32(cell) + i32(open_f[o].w) * step); + for (var q = 0u; q < QT; q = q + 1u) { + t2[q * n + cell] = t2[q * n + src]; + } + } +} + +// Flags a temperature that stopped being a number. +@compute @workgroup_size(256) +fn thermal_check(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: vec3) { + let cell = cell_of(lid, wg); + if cell >= params.n { return; } + if !is_fluid(cell) { return; } + let theta = theta_of(cell); + if theta != theta || abs(theta) > 1e30 { atomicStore(&stats[2], 1u); } } // One workgroup per group: the sum of its links' contributions. @@ -356,7 +529,7 @@ fn cell_moments(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id j = j + vec3(c_of(q)) * v; } let rho = 1.0 + drho; - let u = (j + 0.5 * rho * params.acc.xyz) / rho; + let u = (j + 0.5 * rho * acc_of(cell)) / rho; moments[cell] = vec4(u, rho); } diff --git a/crates/lbm-gpu/src/lib.rs b/crates/lbm-gpu/src/lib.rs index 7db16b5..9fecc00 100644 --- a/crates/lbm-gpu/src/lib.rs +++ b/crates/lbm-gpu/src/lib.rs @@ -1,7 +1,7 @@ //! The lattice Boltzmann solver on the GPU: [`GpuFlow`] steps a flow that //! [`fieldcraft_lbm::Flow`] prepared (grid, links, open cells, initial state) with the compute //! shaders of `lbm.wgsl`, which are the CPU scheme's collision, streaming rules, open boundaries -//! and momentum exchange written for the device; `fieldcraft_lbm::run` drives it through the +//! and momentum exchange, and with heat its temperature lattice, written for the device; `fieldcraft_lbm::run` drives it through the //! [`Backend`] trait in chunks of steps, so populations stay on the device and only frames, //! force histories and the stability check come back. //! diff --git a/crates/lbm-gpu/src/native.rs b/crates/lbm-gpu/src/native.rs index c94b3be..75f2432 100644 --- a/crates/lbm-gpu/src/native.rs +++ b/crates/lbm-gpu/src/native.rs @@ -7,6 +7,7 @@ use fieldcraft_lbm::backend::{Backend, Fields}; use fieldcraft_lbm::lattice::{C, Q, W}; use fieldcraft_lbm::setup::MAX_LATTICE_SPEED; use fieldcraft_lbm::{Error, Flow, Group, Result, Rule, Target}; +use fieldcraft_model::ThermalCondition; use wgpu::util::DeviceExt; const SHADER: &str = include_str!("lbm.wgsl"); @@ -16,6 +17,9 @@ const KIND_PULL: u32 = 0; const KIND_BOUNCE: u32 = 1; const KIND_REFLECT: u32 = 2; const KIND_OPEN: u32 = 3; +/// The scalar lattice's directions (the first seven of D3Q19) and weights. +const QT: usize = 7; +const WT: [f64; QT] = [0.25, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]; fn gpu_error(e: impl std::fmt::Display) -> Error { Error::Invalid(format!("GPU: {e}")) @@ -91,14 +95,18 @@ pub struct GpuFlow { forces: wgpu::Buffer, partials: wgpu::Buffer, partial_count: usize, + /// The scalar populations, A–B like the flow's, with heat; and T₀ they are shifted by. + thermal: Option<[wgpu::Buffer; 2]>, + initial_temperature: f64, fluid: Vec, fluid_cells: usize, groups: Vec, empty_obstacles: Vec, steps_done: u64, - /// The forces of the last step, once read back (zero before the first step). - forces_cache: Option>, - /// Each group's share of the exchange from the populations' weights, which the kernel leaves + /// The momentum (xyz) and heat (w) each group received in the last step, once read back + /// (zero before the first step). + exchange_cache: Option>, + /// Each group's share of the momentum from the populations' weights, which the kernel leaves /// out (`link_forces` sums the deviations): added to what it reads back. weight_forces: Vec<[f64; 3]>, } @@ -113,9 +121,6 @@ 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))); }; @@ -157,6 +162,20 @@ impl GpuFlow { } } } + // The scalar populations as deviations from w_q T₀ (see lbm.wgsl). + let thermal = flow.thermal(); + let t0 = thermal.map_or(0.0, |t| t.initial()); + let mut tpop = filled(if thermal.is_some() { QT * n } else { 0 }, 0.0f32, "temperature populations")?; + if thermal.is_some() { + for q in 0..QT { + let w = WT.get(q).copied().unwrap_or(0.0) * t0; + for cell in 0..n { + if let (Some(slot), true) = (tpop.get_mut(q * n + cell), flow.is_fluid(cell)) { + *slot = (flow.thermal_population(q, cell) - w) as f32; + } + } + } + } let mut fluid = filled(n, false, "cell flags")?; let mut fluid_u32 = filled(n, 0u32, "cell flags")?; for cell in 0..n { @@ -213,6 +232,16 @@ impl GpuFlow { } } let by_group: Vec = by_group.iter().map(|(_, l)| *l).collect(); + // Each group's thermal condition (held temperatures shifted by T₀). + let mut group_t = reserved(4 * group_count, "thermal conditions")?; + for g in 0..group_count { + let (kind, value) = match thermal.and_then(|t| t.groups().get(g)).copied().unwrap_or_default() { + ThermalCondition::Adiabatic {} => (0.0, 0.0), + ThermalCondition::Temperature { value } => (1.0, (value - t0) as f32), + ThermalCondition::Flux { value } => (2.0, value as f32), + }; + group_t.extend_from_slice(&[kind, value, 0.0, 0.0]); + } // The open cells. let open = flow.open_cells(); @@ -226,6 +255,13 @@ impl GpuFlow { open_u.extend_from_slice(&[oc.cell, oc.axis as u32, kind, 0]); open_f.extend_from_slice(&[target[0], target[1], target[2], oc.inward as f32]); } + let mut open_t = reserved(4 * open.len(), "open cell temperatures")?; + for k in 0..open.len() { + match thermal.and_then(|t| t.open().get(k)).copied().flatten() { + Some(value) => open_t.extend_from_slice(&[1.0, (value - t0) as f32, 0.0, 0.0]), + None => open_t.extend_from_slice(&[0.0; 4]), + } + } let (omega_plus, omega_minus) = flow.omegas(); let acc = flow.acceleration(); @@ -261,6 +297,19 @@ impl GpuFlow { ] { params.extend_from_slice(&v.to_le_bytes()); } + let buoyancy = thermal.map_or([0.0; 3], |t| t.buoyancy()); + for v in [ + buoyancy[0] as f32, + buoyancy[1] as f32, + buoyancy[2] as f32, + if buoyancy != [0.0; 3] { 1.0 } else { 0.0 }, + thermal.map_or(0.0, |t| t.omega()) as f32, + t0 as f32, + if thermal.is_some() { 1.0 } else { 0.0 }, + 0.0, + ] { + params.extend_from_slice(&v.to_le_bytes()); + } // Buffers, with the device's out-of-memory errors caught. let scopes = [wgpu::ErrorFilter::OutOfMemory, wgpu::ErrorFilter::Validation].map(|f| device.push_error_scope(f)); @@ -297,6 +346,14 @@ impl GpuFlow { let group_start = storage("lbm group start", &bytes_u32(&group_start)); let forces = empty("lbm forces", 16 * group_count as u64); let partials = empty("lbm partials", 4 * partial_count as u64); + let therm_bytes = elem * (QT * n) as u64; + let therm = if thermal.is_some() { + [storage("lbm t", &bytes_f32(&tpop)), empty("lbm t2", therm_bytes)] + } else { + [empty("lbm t", 16), empty("lbm t2", 16)] + }; + let group_t = storage("lbm thermal conditions", &bytes_f32(&group_t)); + let open_t = storage("lbm open temperatures", &bytes_f32(&open_t)); let module = device .create_shader_module(wgpu::ShaderModuleDescriptor { label: Some("fieldcraft lbm"), source: wgpu::ShaderSource::Wgsl(SHADER.into()) }); @@ -326,6 +383,10 @@ impl GpuFlow { entry(13, read), entry(14, write), entry(15, write), + entry(16, write), + entry(17, write), + entry(18, read), + entry(19, read), ], }); let bind = [0usize, 1].map(|k| { @@ -348,6 +409,10 @@ impl GpuFlow { wgpu::BindGroupEntry { binding: 13, resource: group_start.as_entire_binding() }, wgpu::BindGroupEntry { binding: 14, resource: forces.as_entire_binding() }, wgpu::BindGroupEntry { binding: 15, resource: partials.as_entire_binding() }, + wgpu::BindGroupEntry { binding: 16, resource: therm[k].as_entire_binding() }, + wgpu::BindGroupEntry { binding: 17, resource: therm[1 - k].as_entire_binding() }, + wgpu::BindGroupEntry { binding: 18, resource: group_t.as_entire_binding() }, + wgpu::BindGroupEntry { binding: 19, resource: open_t.as_entire_binding() }, ], }) }); @@ -357,7 +422,21 @@ impl GpuFlow { immediate_size: 0, }); let mut pipelines = HashMap::new(); - for name in ["collide", "stream", "stream_links", "open_boundaries", "link_forces", "group_forces", "cell_moments", "mass"] { + for name in [ + "collide", + "stream", + "stream_links", + "open_boundaries", + "link_forces", + "group_forces", + "cell_moments", + "mass", + "collide_thermal", + "stream_thermal", + "stream_links_thermal", + "open_thermal", + "thermal_check", + ] { let p = device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor { label: Some(name), layout: Some(&pl), @@ -393,12 +472,14 @@ impl GpuFlow { forces, partials, partial_count, + thermal: thermal.map(|_| therm), + initial_temperature: t0, fluid_cells: flow.fluid_cells(), fluid, groups: flow.groups().to_vec(), empty_obstacles: flow.empty_obstacles().to_vec(), steps_done: 0, - forces_cache: Some(vec![[0.0; 3]; group_count]), + exchange_cache: Some(vec![[0.0; 4]; group_count]), weight_forces, }) } @@ -452,6 +533,22 @@ impl GpuFlow { Ok(bytes.as_chunks::<4>().0.iter().take(count).map(|b| f32::from_le_bytes(*b)).collect()) } + /// The momentum (xyz) and heat (w) each group received in the last step. + fn exchange(&mut self) -> Result<&[[f64; 4]]> { + if self.exchange_cache.is_none() { + let v = self.read_f32(&self.forces, 4 * self.group_count)?; + self.exchange_cache = Some( + v.as_chunks::<4>() + .0 + .iter() + .zip(&self.weight_forces) + .map(|(g, w)| [f64::from(g[0]) + w[0], f64::from(g[1]) + w[1], f64::from(g[2]) + w[2], f64::from(g[3])]) + .collect(), + ); + } + Ok(self.exchange_cache.as_deref().unwrap_or_default()) + } + /// The moments of the current state into the moments buffer. fn compute_moments(&self) -> Result<()> { let mut enc = self.device.create_command_encoder(&wgpu::CommandEncoderDescriptor { label: Some("lbm moments") }); @@ -480,17 +577,31 @@ impl GpuFlow { if self.open_count > 0 { self.dispatch(&mut pass, self.current, "open_boundaries", self.open_groups)?; } + if self.thermal.is_some() { + self.dispatch(&mut pass, self.current, "collide_thermal", self.cells)?; + self.dispatch(&mut pass, self.current, "stream_thermal", self.cells)?; + if self.link_count > 0 { + self.dispatch(&mut pass, self.current, "stream_links_thermal", self.links_groups)?; + } + if self.open_count > 0 { + self.dispatch(&mut pass, self.current, "open_thermal", self.open_groups)?; + } + } self.current ^= 1; } - // The momentum exchange of the last step, from its post-collision populations. + // The momentum (and heat) exchange of the last step, from its post-collision + // populations. if self.link_count > 0 && self.group_count > 0 { self.dispatch(&mut pass, self.current ^ 1, "link_forces", self.links_groups)?; self.dispatch(&mut pass, self.current ^ 1, "group_forces", self.force_groups)?; } + if self.thermal.is_some() { + self.dispatch(&mut pass, self.current, "thermal_check", self.cells)?; + } } self.queue.submit([enc.finish()]); self.steps_done += steps; - self.forces_cache = None; + self.exchange_cache = None; let stats = self.read(&self.stats, 16)?; let word = |k: usize| stats.get(4 * k..4 * k + 4).and_then(|b| b.try_into().ok()).map(u32::from_le_bytes).unwrap_or(0); let max_u2 = f32::from_bits(word(0)); @@ -501,6 +612,12 @@ impl GpuFlow { max_u2.sqrt() ))); } + if word(2) != 0 { + return Err(Error::Unstable(format!( + "the heat went unstable by time step {}: a temperature stopped being a number; lower the time step, refine the spacing or check the heat fluxes", + self.steps_done + ))); + } Ok(()) } } @@ -531,12 +648,25 @@ 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(); - // No temperature lattice on the device yet (the CPU carries heat). - Ok(Fields { density, velocity, temperature: Vec::new() }) + // The temperature from the scalar populations, summed on the host: T₀ + Σ t (T₀ in + // solid cells, whose populations stay zero). + let temperature = match &self.thermal { + None => Vec::new(), + Some(therm) => { + let t = self.read_f32(therm.get(self.current).ok_or_else(|| Error::Invalid("GPU: no thermal buffer".into()))?, QT * self.n)?; + (0..self.n) + .map(|cell| self.initial_temperature + (0..QT).map(|q| t.get(q * self.n + cell).copied().map_or(0.0, f64::from)).sum::()) + .collect() + } + }; + Ok(Fields { density, velocity, temperature }) } fn heat(&mut self) -> Result> { - Ok(Vec::new()) + if self.thermal.is_none() { + return Ok(Vec::new()); + } + Ok(self.exchange()?.iter().map(|g| g[3]).collect()) } fn mass(&mut self) -> Result { @@ -553,19 +683,7 @@ impl Backend for GpuFlow { } fn forces(&mut self) -> Result> { - if let Some(f) = &self.forces_cache { - return Ok(f.clone()); - } - let v = self.read_f32(&self.forces, 4 * self.group_count)?; - let forces: Vec<[f64; 3]> = v - .as_chunks::<4>() - .0 - .iter() - .zip(&self.weight_forces) - .map(|(g, w)| [f64::from(g[0]) + w[0], f64::from(g[1]) + w[1], f64::from(g[2]) + w[2]]) - .collect(); - self.forces_cache = Some(forces.clone()); - Ok(forces) + Ok(self.exchange()?.iter().map(|g| [g[0], g[1], g[2]]).collect()) } fn groups(&self) -> &[Group] { diff --git a/crates/lbm-gpu/tests/verification.rs b/crates/lbm-gpu/tests/verification.rs index 5f70a6f..4b6c71f 100644 --- a/crates/lbm-gpu/tests/verification.rs +++ b/crates/lbm-gpu/tests/verification.rs @@ -3,9 +3,10 @@ //! deviations from their weights in single precision (measured on lavapipe at 2–5e-6, with //! a margin): velocity fields within 2e-5 of the largest speed, densities within 1e-6, wall //! forces within 5e-5 of their magnitude (a wall's shear within 5e-5 of itself), mass drift below -//! 1e-6. Skipped without a GPU adapter -//! unless `FIELDCRAFT_REQUIRE_GPU=1`, which makes a missing adapter a failure; `FIELDCRAFT_GPU` picks -//! the adapter (`hardware` refuses a software one). +//! 1e-6; with heat, temperatures within 2e-5 of the field's range and wall heats within 1e-4 of +//! the largest (a heat is a difference of populations, so it keeps fewer digits). Skipped without +//! a GPU adapter unless `FIELDCRAFT_REQUIRE_GPU=1`, which makes a missing adapter a failure; +//! `FIELDCRAFT_GPU` picks the adapter (`hardware` refuses a software one). #[path = "../../lbm/tests/common/mod.rs"] mod common; @@ -15,7 +16,9 @@ use fieldcraft_gpu::select::{Policy, gpu_required}; use fieldcraft_lbm::backend::Backend; use fieldcraft_lbm::{Flow, Setup}; use fieldcraft_lbm_gpu::GpuFlow; -use fieldcraft_model::{BoxFace, FaceCondition, Procedure, Profile}; +use fieldcraft_model::{ + BoxFace, Directional, FaceCondition, Fluid, FluidDomain, FluidThermal, Material, Model, Procedure, Profile, Step, Thermal, ThermalCondition, +}; const PERIODIC: FaceCondition = FaceCondition::Periodic {}; @@ -29,10 +32,12 @@ fn gpu(flow: &Flow) -> Option { Some(g) } Err(e) if gpu_required() => panic!("a GPU is required: {e}"), - Err(e) => { + // Only a missing device skips: a shader that does not compile is a failure anywhere. + Err(e) if e.to_string().contains("no GPU device") => { println!("skipping: {e}"); None } + Err(e) => panic!("{e}"), } } @@ -205,3 +210,139 @@ fn zou_he_channel_matches_the_cpu() { println!("Zou–He channel: CPU/GPU velocity {du:.2e}, density {drho:.2e}"); assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); } + +// ---- Heat ---- + +/// A fluid with heat: density ρ, kinematic viscosity ν, thermal diffusivity α (conductivity +/// ρ α, unit specific heat). +fn hot_model(density: f64, nu: f64, alpha: f64, domain: 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 +} + +/// The largest difference of the temperature fields over the CPU's temperature range, and of +/// the groups' heats over the largest CPU heat. +fn heat_difference(cpu: &mut Flow, gpu: &mut GpuFlow) -> (f64, f64) { + let (a, b) = (cpu.fields().unwrap(), gpu.fields().unwrap()); + assert_eq!(a.temperature.len(), b.temperature.len()); + let (lo, hi) = a.temperature.iter().fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), t| (lo.min(*t), hi.max(*t))); + let dt = a.temperature.iter().zip(&b.temperature).map(|(x, y)| (x - y).abs()).fold(0.0, f64::max) / (hi - lo).max(1e-30); + let (ha, hb) = (cpu.heat().unwrap(), gpu.heat().unwrap()); + assert_eq!(ha.len(), hb.len()); + let scale = ha.iter().map(|h| h.abs()).fold(0.0, f64::max).max(1e-30); + let dh = ha.iter().zip(&hb).map(|(x, y)| (x - y).abs()).fold(0.0, f64::max) / scale; + (dt, dh) +} + +/// Plates held at 20 and 21 with the fluid at rest: the steady temperature is linear between +/// them (anti-bounce-back puts the walls halfway between cells), the heat through both walls +/// is k ΔT A / H, and the GPU agrees with the CPU although it stores the small differences +/// about T₀ = 20 in single precision. +#[test] +fn conduction_between_plates_matches_the_analytic_profile_and_the_cpu() { + let (h, alpha, n) = (1.0, 0.05, 16usize); + 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 { initial_temperature: 20.0, ..FluidThermal::default() }; + t.faces.y_min = ThermalCondition::Temperature { value: 20.0 }; + t.faces.y_max = ThermalCondition::Temperature { value: 21.0 }; + d.thermal = Some(t); + let m = hot_model(1.0, 0.1, alpha, d, step(8.0 * h * h / alpha, dx * dx / (8.0 * alpha), None)); + let (setup, cpu) = prepare(&m, &[]); + let Some((mut cpu, mut gpu)) = both(&setup, cpu) else { return }; + let f = gpu.fields().unwrap(); + let grid = *cpu.grid(); + let mut err: f64 = 0.0; + for j in 0..n { + let expected = 20.0 + (j as f64 + 0.5) / n as f64; + err = err.max((f.temperature[grid.index([0, j, 0])] - expected).abs()); + } + let q = alpha * 1.0 / h * dx * dx; + let heat = gpu.heat().unwrap(); + let cold = cpu.groups().iter().position(|g| g.name == "y_min").unwrap(); + let hot = cpu.groups().iter().position(|g| g.name == "y_max").unwrap(); + // lattice heat (temperature × faces per step) to model units: ρ c h³ / dt + let scale = setup.thermal.as_ref().unwrap().heat_scale; + let (dt, dh) = heat_difference(&mut cpu, &mut gpu); + println!( + "conduction: profile error {err:.2e}; heat {:.6e} / {:.6e} vs {q:.6e}; GPU vs CPU: temperature {dt:.2e} heat {dh:.2e}", + heat[cold] * scale, + heat[hot] * scale + ); + assert!(err < 5e-5, "{err}"); + assert!((heat[cold] * scale / q - 1.0).abs() < 1e-4 && (heat[hot] * scale / q + 1.0).abs() < 1e-4, "{heat:?}"); + assert!(dt < 2e-5 && dh < 1e-4, "{dt} {dh}"); +} + +/// A heated channel: a velocity inlet carrying 1 in, a pressure outlet, a flux into the fluid +/// at the bottom wall, the top wall held at 0, periodic across. Every rule of the scalar +/// lattice (pulls, the three wall conditions, the inlet and the outlet) agrees with the CPU. +#[test] +fn a_heated_channel_matches_the_cpu() { + let (h, l, umax, nu, alpha) = (1.0, 2.0, 0.1, 0.05, 0.02); + 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: Profile::Parabolic }; + d.faces.x_max = FaceCondition::Pressure { pressure: 0.0 }; + d.faces.z_min = PERIODIC; + d.faces.z_max = PERIODIC; + let mut t = FluidThermal::default(); + t.faces.x_min = ThermalCondition::Temperature { value: 1.0 }; + t.faces.y_min = ThermalCondition::Flux { value: 0.05 }; + t.faces.y_max = ThermalCondition::Temperature { value: 0.0 }; + d.thermal = Some(t); + let m = hot_model(1.0, nu, alpha, d, step(2.0 * l / umax, 0.05 * dx / umax, None)); + let (setup, cpu) = prepare(&m, &[]); + let Some((mut cpu, mut gpu)) = both(&setup, cpu) else { return }; + let (du, drho) = field_difference(&mut cpu, &mut gpu); + let (dt, dh) = heat_difference(&mut cpu, &mut gpu); + println!("heated channel: GPU vs CPU: velocity {du:.2e} density {drho:.2e} temperature {dt:.2e} heat {dh:.2e}"); + assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); + assert!(dt < 2e-5 && dh < 1e-4, "{dt} {dh}"); +} + +/// A differentially heated cavity with buoyancy (the de Vahl Davis geometry at Ra 10³ on 16 +/// cells, over the first diffusive time): the temperature drives the flow on both backends +/// alike. +#[test] +fn natural_convection_matches_the_cpu() { + let (h, pr, alpha, n, u_max) = (1.0, 0.71, 1.0, 16usize, 3.65); + let nu = pr * alpha; + let g_beta = 1e3 * nu * alpha / (h * h * h); + let dx = h / n as f64; + 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); + 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(h * h / alpha, dt, None)); + let (setup, cpu) = prepare(&m, &[]); + let Some((mut cpu, mut gpu)) = both(&setup, cpu) else { return }; + let (du, drho) = field_difference(&mut cpu, &mut gpu); + let (dt, dh) = heat_difference(&mut cpu, &mut gpu); + // the fields are in lattice units: the time step puts the cavity's fastest flow at 0.05 + let speed = cpu.fields().unwrap().velocity.iter().map(|v| (v[0] * v[0] + v[1] * v[1]).sqrt()).fold(0.0, f64::max); + println!("cavity: fastest lattice speed {speed:.3e}; GPU vs CPU: velocity {du:.2e} density {drho:.2e} temperature {dt:.2e} heat {dh:.2e}"); + assert!(speed > 0.01, "the buoyancy drives a flow: {speed}"); + assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); + assert!(dt < 2e-5 && dh < 1e-4, "{dt} {dh}"); +} diff --git a/crates/lbm/src/flow.rs b/crates/lbm/src/flow.rs index dd84b2a..ed21ad1 100644 --- a/crates/lbm/src/flow.rs +++ b/crates/lbm/src/flow.rs @@ -143,6 +143,33 @@ pub struct ThermalState { heat: Vec, } +impl ThermalState { + /// The scalar relaxation rate 1 / τ_T. + pub fn omega(&self) -> f64 { + self.omega + } + + /// The initial temperature (and the reference of the buoyancy). + pub fn initial(&self) -> f64 { + self.initial + } + + /// `−β g` in lattice units; zero without buoyancy. + pub fn buoyancy(&self) -> [f64; 3] { + self.buoyancy + } + + /// Each force group's thermal condition (fluxes in lattice units). + pub fn groups(&self) -> &[ThermalCondition] { + &self.groups + } + + /// Per open cell, the temperature it holds, or `None` to take its neighbour's. + pub fn open(&self) -> &[Option] { + &self.open + } +} + /// The flow state. #[derive(Debug, Clone)] pub struct Flow { @@ -474,6 +501,16 @@ impl Flow { self.thermal.is_some() } + /// The temperature lattice and its settings, with heat. + pub fn thermal(&self) -> Option<&ThermalState> { + self.thermal.as_ref() + } + + /// Scalar population `q` (0..7) of a cell; zero without heat. + pub fn thermal_population(&self, q: usize, cell: usize) -> f64 { + self.thermal.as_ref().and_then(|t| t.g.get(q * self.stride + cell).copied()).unwrap_or(0.0) + } + /// 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) diff --git a/docs/architecture.md b/docs/architecture.md index 4b082dd..3d6703a 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -33,7 +33,7 @@ crates in strictly lower layers. | L5 | `fieldcraft-io` | `.fieldcraft` files (docs/formats/fieldcraft.md): text model definition, binary mesh, geometry references with content hashes, the journal; atomic saves; autosave storage for crash recovery; read and write roots that confine an automation client's file access (`roots`). Bytes in, bytes out; the file system only in `fs`, `recovery` and `roots` (in the browser, files kept in memory: `memory`). Result files later. | started | | L6 | `fieldcraft-engine` | Commands (`CommandSpec`), `Session::execute` with rollback and the panic guard, scripts; every file a command names checked against the session's roots (`Session::file_access`); headless pictures of the session (`view`: model, geometry, deformed shape and contours on the CPU rasterizer). | started | | L7 | `fieldcraft-gpu` | The viewport's GPU pipelines (lit faces, wide lines, MSAA, depth) without a window toolkit, the adapter policy (`FIELDCRAFT_GPU`: hardware GPU, software GPU, CPU rasterizer) and offscreen pictures (`fieldcraft-cli snapshot`). The window draws through the same pipelines. | started | -| L7 | `fieldcraft-lbm-gpu` | The lattice Boltzmann solver's GPU backend: the CPU scheme's step (TRT collision with Guo forcing, the per-link wall rules, Zou–He boundaries, momentum exchange) as wgpu compute shaders, driven through `lbm`'s `Backend` trait from the flow the CPU prepared; single precision, deterministic per device, checked against the CPU within its tests' tolerances; opt-in per solve (`solve.fluid` with `backend: "gpu"`). Opens its device with `gpu`'s adapter policy. | started | +| L7 | `fieldcraft-lbm-gpu` | The lattice Boltzmann solver's GPU backend: the CPU scheme's step (TRT collision with Guo forcing, the per-link wall rules, Zou–He boundaries, momentum exchange, and the temperature lattice with its wall conditions, buoyancy and heat exchange) as wgpu compute shaders, driven through `lbm`'s `Backend` trait from the flow the CPU prepared; single precision, deterministic per device, checked against the CPU within its tests' tolerances; opt-in per solve (`solve.fluid` with `backend: "gpu"`). Opens its device with `gpu`'s adapter policy. | started | | L7 | `fieldcraft-ui-egui`, `fieldcraft-mcp` | Desktop UI in the Craft family's layout (title bar with quick actions, workspace ribbon from the catalogue, browser, viewport with navigation bar and orientation gizmo, context panel, command search and command line, journal timeline; native file pickers as host services, background opening and solving; results display with legend and probe; offscreen capture of the whole window) and control channel requests (docs/control-protocol.md); MCP server over stdio (docs/mcp.md): tools over engine commands, a headless session or the running app through its control channel, read and write roots. | ui, mcp started | | testkit | `fieldcraft-testkit` | Test tooling, a dev-dependency only: the meshing test corpus (`tests/corpus/*/recipe.json`, analytic shapes) and its runner (`cargo xtask corpus`). | started | | apps | `fieldcraft`, `fieldcraft-cli`, `fieldcraft-web` | Desktop, command line, browser (docs/web.md). | started |