From 60fd0998f55ee08afba6067a8724269cf6fc82a0 Mon Sep 17 00:00:00 2001 From: Christian Svensson Date: Sat, 17 Oct 2026 16:30:06 +0200 Subject: [PATCH] lbm-gpu: stream in place with the AA pattern The collide and stream passes each moved the whole population field through two buffers (304 bytes per cell update). The AA pattern (Bailey et al. 2009; Wittmann et al. 2013) keeps one buffer: even steps collide each cell in place with its directions swapped, odd steps read a cell's incoming populations at its neighbours and write the outgoing ones where the next step reads them, so each thread writes the addresses it read and streaming happens inside the collision (152 bytes). A cell's coordinates are computed once and neighbours wrap with a compare and an add. The per-link wall rules gather their sources before a scatter writes them, so no link overwrites another's post-collision value, and the momentum exchange is taken in the gather of the step whose forces are read back. Options::streaming picks the layout; Auto, what a solve uses, streams in place, which measured faster on every device tried: 1.5x on an RTX 4090, 1.85x on the GB10 and 1.45x on an M4 Max at periodic 128^3. Options::profile times each kernel with timestamp queries for the throughput test, which now runs both layouts. The verification tests run in both layouts, with two new ones comparing fields, mass, forces and heat with the CPU after odd and even step counts across submission boundaries. Part of #562. --- crates/lbm-gpu/src/lbm.wgsl | 154 ++++++-- crates/lbm-gpu/src/lib.rs | 2 +- crates/lbm-gpu/src/native.rs | 235 +++++++++--- crates/lbm-gpu/tests/throughput.rs | 43 ++- crates/lbm-gpu/tests/verification.rs | 517 ++++++++++++++++----------- docs/architecture.md | 2 +- docs/lbm-gpu-performance.md | 49 +++ 7 files changed, 698 insertions(+), 304 deletions(-) create mode 100644 docs/lbm-gpu-performance.md diff --git a/crates/lbm-gpu/src/lbm.wgsl b/crates/lbm-gpu/src/lbm.wgsl index 2349dbf..1af944c 100644 --- a/crates/lbm-gpu/src/lbm.wgsl +++ b/crates/lbm-gpu/src/lbm.wgsl @@ -2,8 +2,15 @@ // time collision with Guo forcing, pull streaming with the per-link rules the CPU prepared // (bounce-back with Ladd's moving-wall term, specular reflection, open placeholders), the Zou–He // open boundaries, and the momentum exchange per force group. Populations are stored direction- -// major (f[q n + cell]) in two buffers, A–B: a step collides `f` in place, streams `f` into `g` -// and applies the open boundaries to `g`; the next step binds the buffers the other way round. +// major. The default AA pattern uses one population buffer (Bailey et al., ICPP 2009; +// Wittmann et al., Comput. Math. Appl. 65 (2013), 924–935). Even steps read a cell and +// write its opposite directions; odd steps read opposite directions at the upstream +// neighbours and write each outgoing direction at the downstream neighbour. Each thread +// writes exactly the addresses it read, so distinct cells never race. Physical coordinates +// wrap; the existing per-link rules still supply walls, reflection and open placeholders. +// A gather pass saves those boundary values before a scatter pass overrides their targets: +// overwriting in the gather itself could destroy another link's post-collision source. +// The optional A–B diagnostic path retains the original two-buffer collide/pull scheme. // // Single precision, with every population stored as its deviation from the weight, h_q = f_q − w_q // (Krüger et al. 2017, §13.3.4: the deviations are small, so f32 keeps their digits instead of @@ -20,6 +27,8 @@ // 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 θ. +override AA_MODE: bool = true; + const Q: u32 = 19u; const QT: u32 = 7u; const WG: u32 = 256u; @@ -47,6 +56,8 @@ struct Params { buoy: vec4, // The scalar lattice: its relaxation rate (x), T₀ (y), whether there is heat (z). therm: vec4, + // Current AA orientation (x: 0 normal, 1 opposite directions at upstream neighbours). + phase: vec4, }; @group(0) @binding(0) var params: Params; @@ -71,7 +82,7 @@ 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). +// The scalar populations remain A–B (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 @@ -164,6 +175,60 @@ fn is_fluid(cell: u32) -> bool { return fluid[cell] != 0u; } +// A cell's grid coordinates (integer divisions once per cell, not once per direction). +fn coords(cell: u32) -> vec3 { + let plane = params.nx * params.ny; + return vec3(i32(cell % params.nx), i32((cell % plane) / params.nx), i32(cell / plane)); +} + +// The cell at `p`, wrapped periodically. A lattice step moves at most one cell per axis, so a +// compare and an add per axis replace the remainder. +fn wrapped(p: vec3) -> u32 { + let size = vec3(i32(params.nx), i32(params.ny), i32(params.nz)); + var w = select(p, p + size, p < vec3(0)); + w = select(w, w - size, w >= size); + return u32((w.z * size.y + w.y) * size.x + w.x); +} + +// Where population q of the cell at `xyz` is stored in the given AA orientation: in place when +// even, reversed at the upstream neighbour when odd. +fn population_index(q: u32, cell: u32, xyz: vec3, phase: u32) -> u32 { + if phase == 0u { return q * params.n + cell; } + return opp(q) * params.n + wrapped(xyz - c_of(q)); +} + +fn population(q: u32, cell: u32, xyz: vec3) -> f32 { + if AA_MODE { return f[population_index(q, cell, xyz, params.phase.x)]; } + return f[q * params.n + cell]; +} + +// Collision's output at its source cell, before the boundary rules. +fn post_index(q: u32, cell: u32, xyz: vec3) -> u32 { + if !AA_MODE { return q * params.n + cell; } + if params.phase.x == 0u { return opp(q) * params.n + cell; } + return q * params.n + wrapped(xyz + c_of(q)); +} + +fn post(q: u32, cell: u32) -> f32 { + return f[post_index(q, cell, aa_coords(cell))]; +} + +fn streamed(q: u32, cell: u32, xyz: vec3) -> f32 { + if AA_MODE { return f[population_index(q, cell, xyz, params.phase.x ^ 1u)]; } + return g[q * params.n + cell]; +} + +fn store_streamed(q: u32, cell: u32, xyz: vec3, v: f32) { + if AA_MODE { f[population_index(q, cell, xyz, params.phase.x ^ 1u)] = v; } + else { g[q * params.n + cell] = v; } +} + +// The coordinates the AA addressing needs; none in the two-buffer layout. +fn aa_coords(cell: u32) -> vec3 { + if AA_MODE { return coords(cell); } + return vec3(0); +} + // Two-relaxation-time collision with Guo forcing, in place in f (Ginzburg et al. 2008; Guo et // al. 2002; Krüger et al. 2017 ch. 10), and the chunk's largest speed. @compute @workgroup_size(256) @@ -183,11 +248,14 @@ fn collide(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: fn collide_cell(cell: u32, lid: u32) { let n = params.n; + // Only the odd AA orientation addresses neighbours. + var xyz = vec3(0); + if AA_MODE && params.phase.x == 1u { xyz = coords(cell); } var fq: array; var drho = 0.0; var j = vec3(0.0); for (var q = 0u; q < Q; q = q + 1u) { - let v = f[q * n + cell]; + let v = population(q, cell, xyz); fq[q] = v; drho = drho + v; j = j + vec3(c_of(q)) * v; @@ -213,7 +281,7 @@ fn collide_cell(cell: u32, lid: u32) { let eq0 = w_of(0u) * (drho - 1.5 * rho * uu); var out0 = fq[0] - wp * (fq[0] - eq0); if forced { out0 = out0 + sp * w_of(0u) * (-3.0 * uf); } - f[cell] = out0; + f[post_index(0u, cell, xyz)] = out0; for (var q = 1u; q < Q; q = q + 2u) { let o = q + 1u; let cv = vec3(c_of(q)); @@ -232,8 +300,8 @@ fn collide_cell(cell: u32, lid: u32) { s_plus = sp * w * (9.0 * cu * cf - 3.0 * uf); s_minus = sm * w * 3.0 * cf; } - f[q * n + cell] = fq[q] - relax_plus - relax_minus + s_plus + s_minus; - f[o * n + cell] = fq[o] - relax_plus + relax_minus + s_plus - s_minus; + f[post_index(q, cell, xyz)] = fq[q] - relax_plus - relax_minus + s_plus + s_minus; + f[post_index(o, cell, xyz)] = fq[o] - relax_plus + relax_minus + s_plus - s_minus; } } @@ -266,7 +334,17 @@ fn stream(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: // write the same population. @compute @workgroup_size(256) fn stream_links(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: vec3) { - let l = cell_of(lid, wg); + gather_link(cell_of(lid, wg), false); +} + +// The same on the step whose wall forces are read back (a chunk's last): the AA layout +// overwrites the post-collision sources in the scatter, so the momentum exchange is taken here. +@compute @workgroup_size(256) +fn stream_links_exchange(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wg: vec3) { + gather_link(cell_of(lid, wg), true); +} + +fn gather_link(l: u32, exchange: bool) { if l >= params.links { return; } let n = params.n; let lu = link_u[l]; @@ -276,16 +354,30 @@ fn stream_links(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id let cell = lu.w; var v = 0.0; if kind == KIND_PULL { - v = f[q * n + lu.y]; + v = post(q, lu.y); } else if kind == KIND_BOUNCE { - v = f[opp(q) * n + cell] + link_f[l]; + v = post(opp(q), cell) + link_f[l]; } else if kind == KIND_REFLECT { // The mirrored direction has the same weight, so h reflects like f. - v = f[q2 * n + lu.y]; + v = post(q2, lu.y); } else { - v = f[opp(q) * n + cell]; + v = post(opp(q), cell); } - g[q * n + cell] = v; + if AA_MODE { + g[l] = v; + if exchange { contrib[l] = vec4(momentum_exchange(l, lu, q, kind, q2), 0.0); } + } else { + g[q * n + cell] = v; + } +} + +// The gather above must complete before any AA boundary population is overwritten. +@compute @workgroup_size(256) +fn scatter_links(@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]; + store_streamed(lu.x & 0xffu, lu.w, aa_coords(lu.w), g[l]); } // Zou–He velocity and pressure boundaries on the streamed populations in g (Zou & He 1997; @@ -301,6 +393,7 @@ fn open_boundaries(@builtin(local_invocation_index) lid: u32, @builtin(workgroup let axis = ou.y; let sign = ot.w; let a_half = 0.5 * acc_of(cell); + let xyz = aa_coords(cell); // The sums of h over the tangential (c_n = 0) and outgoing (c_n < 0) directions. The sums of // f add the weights' (2/3 and 1/6 for D3Q19, so tangential + 2 outgoing gains exactly 1): the // weight is added after the deviations are summed, not before, which in single precision would @@ -308,7 +401,7 @@ fn open_boundaries(@builtin(local_invocation_index) lid: u32, @builtin(workgroup var tangential = 0.0; var outgoing = 0.0; for (var q = 0u; q < Q; q = q + 1u) { - let v = g[q * n + cell]; + let v = streamed(q, cell, xyz); let cn = f32(c_of(q)[axis]) * sign; if cn == 0.0 { tangential = tangential + v; @@ -333,24 +426,37 @@ fn open_boundaries(@builtin(local_invocation_index) lid: u32, @builtin(workgroup let c = c_of(q); if f32(c[axis]) * sign > 0.0 { let cu = dot(vec3(c), bare); - g[q * n + cell] = g[opp(q) * n + cell] + 6.0 * w_of(q) * rho * cu; + store_streamed(q, cell, xyz, streamed(opp(q), cell, xyz) + 6.0 * w_of(q) * rho * cu); } } for (var t = 0u; t < 3u; t = t + 1u) { if t == axis { continue; } var m = -rho * bare[t]; for (var q = 0u; q < Q; q = q + 1u) { - m = m + g[q * n + cell] * f32(c_of(q)[t]); + m = m + streamed(q, cell, xyz) * f32(c_of(q)[t]); } for (var q = 1u; q < Q; q = q + 1u) { let c = c_of(q); if f32(c[axis]) * sign > 0.0 && c[t] != 0 { - g[q * n + cell] = g[q * n + cell] - 0.5 * m * f32(c[t]); + store_streamed(q, cell, xyz, streamed(q, cell, xyz) - 0.5 * m * f32(c[t])); } } } } +fn momentum_exchange(l: u32, lu: vec4, q: u32, kind: u32, q2: u32) -> vec3 { + var v = vec3(0.0); + if kind == KIND_BOUNCE { + let out = post(opp(q), lu.w); + let back = out + link_f[l]; + v = -vec3(c_of(q)) * (out + back); + } else if kind == KIND_REFLECT { + let out = post(q2, lu.y); + v = vec3(c_of(q2)) * out - vec3(c_of(q)) * out; + } + return v; +} + // The momentum each link transferred to its group in the last step, from the post-collision // populations in f (Ladd 1994; Mei et al. 2002): c_out f* − c_in f'. @compute @workgroup_size(256) @@ -366,16 +472,9 @@ fn link_forces(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) // The deviations only: the weights' share of the exchange (2 w_q c_q per bounced link, which // 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 v = contrib[l].xyz; + if !AA_MODE { v = momentum_exchange(l, lu, q, kind, q2); } var heat = 0.0; - if kind == KIND_BOUNCE { - let out = f[opp(q) * n + cell]; - let back = out + link_f[l]; - v = -vec3(c_of(q)) * (out + back); - } else if kind == KIND_REFLECT { - let out = f[q2 * n + lu.y]; - v = vec3(c_of(q2)) * out - vec3(c_of(q)) * out; - } // 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); @@ -525,8 +624,9 @@ fn cell_moments(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id let n = params.n; var drho = 0.0; var j = vec3(0.0); + let xyz = aa_coords(cell); for (var q = 0u; q < Q; q = q + 1u) { - let v = f[q * n + cell]; + let v = population(q, cell, xyz); drho = drho + v; j = j + vec3(c_of(q)) * v; } diff --git a/crates/lbm-gpu/src/lib.rs b/crates/lbm-gpu/src/lib.rs index 9fecc00..c3e2407 100644 --- a/crates/lbm-gpu/src/lib.rs +++ b/crates/lbm-gpu/src/lib.rs @@ -18,7 +18,7 @@ #[cfg(not(target_arch = "wasm32"))] mod native; #[cfg(not(target_arch = "wasm32"))] -pub use native::{GpuFlow, register}; +pub use native::{GpuFlow, KernelTime, Options, Streaming, register}; /// In the browser there is no device backend yet (WebGPU later): registering is a no-op and /// `solve.fluid` with `backend: "gpu"` reports that none is available. diff --git a/crates/lbm-gpu/src/native.rs b/crates/lbm-gpu/src/native.rs index 75f2432..f997891 100644 --- a/crates/lbm-gpu/src/native.rs +++ b/crates/lbm-gpu/src/native.rs @@ -31,6 +31,9 @@ const MAX_CHUNK: u64 = 64; /// How long one wait for the device may take before the solve gives up on it. const DEVICE_TIMEOUT: std::time::Duration = std::time::Duration::from_secs(120); +/// A bounded chunk's dispatches, shared by the normal and timestamped paths. +type Dispatch = (usize, &'static str, (u32, u32)); + /// A vector of `len` copies of `value`, or [`Error::Limit`] when memory runs out (the host /// arrays are sized by the grid: 16 M cells is over a gigabyte of populations). fn filled(len: usize, value: T, what: &str) -> Result> { @@ -72,13 +75,48 @@ pub fn register() { fieldcraft_lbm::register_backend(Box::new(|flow, _setup| GpuFlow::from_env(flow).map(|g| Box::new(g) as Box))); } +/// How the momentum populations stream (see `lbm.wgsl`). +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum Streaming { + /// In place: measured faster on every device tried (RTX 4090, GB10, M4 Max; #563). + #[default] + Auto, + /// One population buffer, collided and streamed in place (the AA pattern). + InPlace, + /// Two population buffers: collide in place, then pull into the other (A–B). + TwoBuffers, +} + +/// How a [`GpuFlow`] runs; the defaults are what a solve uses. +#[derive(Clone, Copy, Debug, Default)] +pub struct Options { + pub streaming: Streaming, + /// Time every kernel with timestamp queries ([`GpuFlow::kernel_times`]). Each dispatch then + /// gets its own pass and every chunk waits for its timestamps: for measurements, not solves. + pub profile: bool, +} + +/// GPU time per kernel, summed over the dispatches since the flow was created (with +/// [`Options::profile`]). +#[derive(Clone, Debug, Default, PartialEq)] +pub struct KernelTime { + pub dispatches: u64, + pub seconds: f64, +} + /// A prepared flow on the device. pub struct GpuFlow { device: wgpu::Device, queue: wgpu::Queue, + /// With [`Options::profile`]: the GPU time of each kernel so far. + kernel_times: Option>>, + /// How long creating the compute pipelines (compiling the shader) took. + pipeline_seconds: f64, + aa: bool, adapter: String, pipelines: HashMap<&'static str, wgpu::ComputePipeline>, - /// Bind groups with (f, g) = (buffer 0, buffer 1) and the other way round. + /// AB swaps populations; AA keeps one population buffer and changes its orientation. + /// Both modes swap the separate scalar lattice buffers. bind: [wgpu::BindGroup; 2], /// Which bind group's `f` holds the state. current: usize, @@ -95,7 +133,7 @@ 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. + /// The scalar populations remain A–B with heat; and T₀ they are shifted by. thermal: Option<[wgpu::Buffer; 2]>, initial_temperature: f64, fluid: Vec, @@ -121,10 +159,20 @@ 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::with_options(flow, policy, Options::default()) + } + + /// [`GpuFlow::new`] with `options`. + pub fn with_options(flow: &Flow, policy: Policy, options: Options) -> Result { 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))); }; let info = adapter.get_info(); + let profile = options.profile; + let features = if profile { wgpu::Features::TIMESTAMP_QUERY } else { wgpu::Features::empty() }; + if !adapter.features().contains(features) { + return Err(gpu_error("profiling needs timestamp queries, which this adapter does not offer")); + } let have = adapter.limits(); let mut limits = wgpu::Limits::default().using_resolution(have.clone()); limits.max_buffer_size = have.max_buffer_size; @@ -132,13 +180,17 @@ impl GpuFlow { limits.max_storage_buffers_per_shader_stage = have.max_storage_buffers_per_shader_stage.max(limits.max_storage_buffers_per_shader_stage); let (device, queue) = pollster::block_on(adapter.request_device(&wgpu::DeviceDescriptor { label: Some("fieldcraft lbm"), - required_features: wgpu::Features::empty(), + required_features: features, required_limits: limits, ..Default::default() })) .map_err(gpu_error)?; let adapter_name = format!("{}{}", info.name, if tier == select::Tier::Software { " (software)" } else { "" }); + let aa = match options.streaming { + Streaming::Auto | Streaming::InPlace => true, + Streaming::TwoBuffers => false, + }; let n = flow.grid().len(); let elem = 4u64; let pop_bytes = elem * (Q * n) as u64; @@ -328,13 +380,18 @@ impl GpuFlow { mapped_at_creation: false, }) }; - let params = device.create_buffer_init(&wgpu::util::BufferInitDescriptor { - label: Some("lbm params"), - contents: ¶ms, - usage: wgpu::BufferUsages::UNIFORM, + let params = [0u32, 1].map(|phase| { + let mut bytes = params.clone(); + bytes.extend_from_slice(&bytes_u32(&[phase, 0, 0, 0])); + device.create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("lbm params"), + contents: &bytes, + usage: wgpu::BufferUsages::UNIFORM, + }) }); let fluid_buf = storage("lbm fluid", &bytes_u32(&fluid_u32)); - let pops = [storage("lbm f", &bytes_f32(&f)), empty("lbm g", pop_bytes)]; + // AA reuses binding g for boundary scratch only, never another full population field. + let pops = [storage("lbm f", &bytes_f32(&f)), empty("lbm g / AA boundary scratch", if aa { 4 * links.len() as u64 } else { pop_bytes })]; let link_u = storage("lbm links", &bytes_u32(&link_u)); let link_f = storage("lbm link du", &bytes_f32(&link_f)); let open_u = storage("lbm open cells", &bytes_u32(&open_u)); @@ -394,10 +451,10 @@ impl GpuFlow { label: Some("fieldcraft lbm"), layout: &layout, entries: &[ - wgpu::BindGroupEntry { binding: 0, resource: params.as_entire_binding() }, + wgpu::BindGroupEntry { binding: 0, resource: params[k].as_entire_binding() }, wgpu::BindGroupEntry { binding: 1, resource: fluid_buf.as_entire_binding() }, - wgpu::BindGroupEntry { binding: 2, resource: pops[k].as_entire_binding() }, - wgpu::BindGroupEntry { binding: 3, resource: pops[1 - k].as_entire_binding() }, + wgpu::BindGroupEntry { binding: 2, resource: pops[if aa { 0 } else { k }].as_entire_binding() }, + wgpu::BindGroupEntry { binding: 3, resource: pops[if aa { 1 } else { 1 - k }].as_entire_binding() }, wgpu::BindGroupEntry { binding: 5, resource: link_u.as_entire_binding() }, wgpu::BindGroupEntry { binding: 6, resource: link_f.as_entire_binding() }, wgpu::BindGroupEntry { binding: 7, resource: open_u.as_entire_binding() }, @@ -422,10 +479,13 @@ impl GpuFlow { immediate_size: 0, }); let mut pipelines = HashMap::new(); + let pipeline_start = std::time::Instant::now(); for name in [ "collide", "stream", "stream_links", + "stream_links_exchange", + "scatter_links", "open_boundaries", "link_forces", "group_forces", @@ -442,7 +502,7 @@ impl GpuFlow { layout: Some(&pl), module: &module, entry_point: Some(name), - compilation_options: Default::default(), + compilation_options: wgpu::PipelineCompilationOptions { constants: &[("AA_MODE", if aa { 1.0 } else { 0.0 })], ..Default::default() }, cache: None, }); pipelines.insert(name, p); @@ -452,9 +512,13 @@ impl GpuFlow { return Err(gpu_error(e)); } } + let pipeline_seconds = pipeline_start.elapsed().as_secs_f64(); Ok(GpuFlow { device, queue, + kernel_times: profile.then(Default::default), + pipeline_seconds, + aa, adapter: adapter_name, pipelines, bind, @@ -489,6 +553,21 @@ impl GpuFlow { &self.adapter } + /// Whether the populations stream in place (see [`Streaming`]). + pub fn streams_in_place(&self) -> bool { + self.aa + } + + /// How long creating the compute pipelines took, in seconds. + pub fn pipeline_seconds(&self) -> f64 { + self.pipeline_seconds + } + + /// With [`Options::profile`], the GPU time of each kernel so far, by kernel name. + pub fn kernel_times(&self) -> Option> { + self.kernel_times.as_ref().map(|t| t.borrow().clone()) + } + fn pipeline(&self, name: &str) -> Result<&wgpu::ComputePipeline> { self.pipelines.get(name).ok_or_else(|| Error::Invalid(format!("GPU: no kernel {name}"))) } @@ -500,6 +579,57 @@ impl GpuFlow { Ok(()) } + /// Timestamped passes are diagnostic: their extra pass boundaries and readback are not + /// used for throughput measurements. Normally every chunk remains one compute pass. + fn submit(&self, label: &str, commands: &[Dispatch]) -> Result<()> { + let mut enc = self.device.create_command_encoder(&wgpu::CommandEncoderDescriptor { label: Some(label) }); + let Some(times) = &self.kernel_times else { + { + let mut pass = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { label: Some(label), timestamp_writes: None }); + for &(bind, name, groups) in commands { + self.dispatch(&mut pass, bind, name, groups)?; + } + } + self.queue.submit([enc.finish()]); + return Ok(()); + }; + let count = u32::try_from(commands.len() * 2).map_err(|_| gpu_error("too many timestamp queries"))?; + if count == 0 || count > wgpu::QUERY_SET_MAX_QUERIES { + return Err(gpu_error("timestamp query count is outside the device limit")); + } + let query = self.device.create_query_set(&wgpu::QuerySetDescriptor { label: Some(label), ty: wgpu::QueryType::Timestamp, count }); + let result = self.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("lbm timestamps"), + size: u64::from(count) * 8, + usage: wgpu::BufferUsages::QUERY_RESOLVE | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, + }); + for (i, &(bind, name, groups)) in commands.iter().enumerate() { + let mut pass = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some(name), + timestamp_writes: Some(wgpu::ComputePassTimestampWrites { + query_set: &query, + beginning_of_pass_write_index: Some(2 * i as u32), + end_of_pass_write_index: Some(2 * i as u32 + 1), + }), + }); + self.dispatch(&mut pass, bind, name, groups)?; + } + enc.resolve_query_set(&query, 0..count, &result, 0); + self.queue.submit([enc.finish()]); + let bytes = self.read(&result, u64::from(count) * 8)?; + let period = f64::from(self.queue.get_timestamp_period()); + let mut totals = times.borrow_mut(); + for ((_, name, _), pair) in commands.iter().zip(bytes.as_chunks::<16>().0) { + let [begin, end] = pair.as_chunks::<8>().0 else { return Err(gpu_error("invalid timestamp pair")) }; + let ticks = u64::from_le_bytes(*end).checked_sub(u64::from_le_bytes(*begin)).ok_or_else(|| gpu_error("invalid GPU timestamps"))?; + let total = totals.entry(name).or_default(); + total.dispatches += 1; + total.seconds += ticks as f64 * period * 1e-9; + } + Ok(()) + } + /// Copy `size` bytes of `src` back to the host. fn read(&self, src: &wgpu::Buffer, size: u64) -> Result> { let size = size.max(16); @@ -551,13 +681,7 @@ impl GpuFlow { /// 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") }); - { - let mut pass = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { label: Some("lbm moments"), timestamp_writes: None }); - self.dispatch(&mut pass, self.current, "cell_moments", self.cells)?; - } - self.queue.submit([enc.finish()]); - Ok(()) + self.submit("lbm moments", &[(self.current, "cell_moments", self.cells)]) } } @@ -565,41 +689,45 @@ impl GpuFlow { /// One submission of `steps` (≤ MAX_CHUNK) time steps, then the stability check. fn advance_chunk(&mut self, steps: u64) -> Result<()> { self.queue.write_buffer(&self.stats, 0, &bytes_u32(&[0, 0, 0, 0])); - let mut enc = self.device.create_command_encoder(&wgpu::CommandEncoderDescriptor { label: Some("lbm steps") }); - { - let mut pass = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { label: Some("lbm steps"), timestamp_writes: None }); - for _ in 0..steps { - self.dispatch(&mut pass, self.current, "collide", self.cells)?; - self.dispatch(&mut pass, self.current, "stream", self.cells)?; - if self.link_count > 0 { - self.dispatch(&mut pass, self.current, "stream_links", self.links_groups)?; - } - 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; + let mut commands = Vec::with_capacity(8 * MAX_CHUNK as usize + 3); + for step in 0..steps { + commands.push((self.current, "collide", self.cells)); + if !self.aa { + commands.push((self.current, "stream", self.cells)); } - // 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.link_count > 0 { + // The last step's momentum exchange is the one read back. + let exchange = self.aa && step + 1 == steps; + commands.push((self.current, if exchange { "stream_links_exchange" } else { "stream_links" }, self.links_groups)); + if self.aa { + commands.push((self.current, "scatter_links", self.links_groups)); + } + } + if self.open_count > 0 { + commands.push((self.current, "open_boundaries", self.open_groups)); } if self.thermal.is_some() { - self.dispatch(&mut pass, self.current, "thermal_check", self.cells)?; + commands.push((self.current, "collide_thermal", self.cells)); + commands.push((self.current, "stream_thermal", self.cells)); + if self.link_count > 0 { + commands.push((self.current, "stream_links_thermal", self.links_groups)); + } + if self.open_count > 0 { + commands.push((self.current, "open_thermal", self.open_groups)); + } } + self.current ^= 1; } - self.queue.submit([enc.finish()]); + // The momentum (and heat) exchange of the last step, from its post-collision + // populations. + if self.link_count > 0 && self.group_count > 0 { + commands.push((self.current ^ 1, "link_forces", self.links_groups)); + commands.push((self.current ^ 1, "group_forces", self.force_groups)); + } + if self.thermal.is_some() { + commands.push((self.current, "thermal_check", self.cells)); + } + self.submit("lbm steps", &commands)?; self.steps_done += steps; self.exchange_cache = None; let stats = self.read(&self.stats, 16)?; @@ -624,7 +752,7 @@ impl GpuFlow { impl Backend for GpuFlow { fn name(&self) -> String { - format!("gpu: {}", self.adapter) + format!("gpu: {}, {} streaming", self.adapter, if self.aa { "in-place" } else { "two-buffer" }) } fn advance(&mut self, steps: u64) -> Result<()> { @@ -671,12 +799,7 @@ impl Backend for GpuFlow { fn mass(&mut self) -> Result { self.compute_moments()?; - let mut enc = self.device.create_command_encoder(&wgpu::CommandEncoderDescriptor { label: Some("lbm mass") }); - { - let mut pass = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { label: Some("lbm mass"), timestamp_writes: None }); - self.dispatch(&mut pass, self.current, "mass", self.cells)?; - } - self.queue.submit([enc.finish()]); + self.submit("lbm mass", &[(self.current, "mass", self.cells)])?; let partials = self.read_f32(&self.partials, self.partial_count)?; // Summed in a fixed order on the host, so the result is deterministic on this device. Ok(partials.iter().map(|v| f64::from(*v)).sum()) diff --git a/crates/lbm-gpu/tests/throughput.rs b/crates/lbm-gpu/tests/throughput.rs index 35d9940..751b7c9 100644 --- a/crates/lbm-gpu/tests/throughput.rs +++ b/crates/lbm-gpu/tests/throughput.rs @@ -1,6 +1,8 @@ //! Throughput of the GPU backend (ignored; run with //! `cargo test --release -p fieldcraft-lbm-gpu --test throughput -- --ignored --nocapture`), -//! on the cases of `fieldcraft-lbm`'s throughput test. +//! on the cases of `fieldcraft-lbm`'s throughput test, streaming in place and with two buffers +//! (what `Streaming::Auto` picks on this device is marked). `FIELDCRAFT_LBM_GPU_PROFILE=1` +//! adds each kernel's GPU time (timestamp queries; the throughput is then not representative). #[path = "../../lbm/tests/common/mod.rs"] mod common; @@ -8,7 +10,7 @@ mod common; use common::*; use fieldcraft_gpu::select::Policy; use fieldcraft_lbm::backend::Backend; -use fieldcraft_lbm_gpu::GpuFlow; +use fieldcraft_lbm_gpu::{GpuFlow, Options, Streaming}; use fieldcraft_model::{BoxFace, FaceCondition}; #[test] @@ -41,16 +43,33 @@ fn cell_updates_per_second() { }, ); let (_, flow) = prepare(&m, &obstacles); - let Ok(mut gpu) = GpuFlow::new(&flow, Policy::from_env().expect("FIELDCRAFT_GPU names a policy")) else { - println!("skipping: no GPU"); - return; + let policy = Policy::from_env().expect("FIELDCRAFT_GPU names a policy"); + let profile = std::env::var("FIELDCRAFT_LBM_GPU_PROFILE").is_ok_and(|v| v == "1"); + let Ok(auto) = GpuFlow::new(&flow, policy) else { + println!("skipping {name}, {l}³: no GPU or not enough memory"); + continue; }; - // One warm-up chunk, then the timed run: the chunk's steps are one submission. - gpu.advance(1).unwrap(); - let start = std::time::Instant::now(); - gpu.advance(steps).unwrap(); - let _ = gpu.mass().unwrap(); - let s = start.elapsed().as_secs_f64(); - println!("{name}, {l}³ on {}: {:.1} M cell updates/s", gpu.adapter_name(), (gpu.fluid_cells() as f64 * steps as f64) / s / 1e6); + let picked = auto.streams_in_place(); + drop(auto); + for (streaming, in_place) in [(Streaming::InPlace, true), (Streaming::TwoBuffers, false)] { + let mut gpu = GpuFlow::with_options(&flow, policy, Options { streaming, profile }).unwrap(); + // One warm-up chunk, then the timed run: the chunk's steps are one submission. + gpu.advance(1).unwrap(); + let start = std::time::Instant::now(); + gpu.advance(steps).unwrap(); + let _ = gpu.mass().unwrap(); + let s = start.elapsed().as_secs_f64(); + println!( + "{name}, {l}³ on {}, {}{}: {:.1} M cell updates/s (pipelines {:.2} s)", + gpu.adapter_name(), + if in_place { "in place" } else { "two buffers" }, + if in_place == picked { " (auto)" } else { "" }, + (gpu.fluid_cells() as f64 * steps as f64) / s / 1e6, + gpu.pipeline_seconds() + ); + for (kernel, t) in gpu.kernel_times().unwrap_or_default() { + println!(" {kernel}: {} dispatches, {:.3} ms each", t.dispatches, 1e3 * t.seconds / t.dispatches.max(1) as f64); + } + } } } diff --git a/crates/lbm-gpu/tests/verification.rs b/crates/lbm-gpu/tests/verification.rs index 4b6c71f..916a963 100644 --- a/crates/lbm-gpu/tests/verification.rs +++ b/crates/lbm-gpu/tests/verification.rs @@ -15,7 +15,7 @@ use common::*; use fieldcraft_gpu::select::{Policy, gpu_required}; use fieldcraft_lbm::backend::Backend; use fieldcraft_lbm::{Flow, Setup}; -use fieldcraft_lbm_gpu::GpuFlow; +use fieldcraft_lbm_gpu::{GpuFlow, Options, Streaming}; use fieldcraft_model::{ BoxFace, Directional, FaceCondition, Fluid, FluidDomain, FluidThermal, Material, Model, Procedure, Profile, Step, Thermal, ThermalCondition, }; @@ -24,11 +24,11 @@ const PERIODIC: FaceCondition = FaceCondition::Periodic {}; /// The flow on the GPU `FIELDCRAFT_GPU` picks, or `None` to skip (a failure when a GPU is /// required). -fn gpu(flow: &Flow) -> Option { +fn gpu(flow: &Flow, streaming: Streaming) -> Option { let policy = Policy::from_env().expect("FIELDCRAFT_GPU names a policy"); - match GpuFlow::new(flow, policy) { + match GpuFlow::with_options(flow, policy, Options { streaming, profile: false }) { Ok(g) => { - println!("GPU: {}", g.adapter_name()); + println!("GPU: {} ({streaming:?})", g.adapter_name()); Some(g) } Err(e) if gpu_required() => panic!("a GPU is required: {e}"), @@ -42,8 +42,8 @@ fn gpu(flow: &Flow) -> Option { } /// Both backends through every step; the CPU's fields and the GPU's. -fn both(setup: &Setup, mut cpu: Flow) -> Option<(Flow, GpuFlow)> { - let mut gpu = gpu(&cpu)?; +fn both(setup: &Setup, mut cpu: Flow, streaming: Streaming) -> Option<(Flow, GpuFlow)> { + let mut gpu = gpu(&cpu, streaming)?; cpu.advance(setup.steps).unwrap(); gpu.advance(setup.steps).unwrap(); Some((cpu, gpu)) @@ -66,149 +66,165 @@ fn field_difference(cpu: &mut Flow, gpu: &mut GpuFlow) -> (f64, f64) { #[test] fn plane_poiseuille_matches_the_analytic_profile_and_the_cpu() { - let (h, nu, umax, n) = (1.0, 1.0 / 6.0, 0.1, 16usize); - let dx = h / n as f64; - let a = 8.0 * nu * umax / (h * h); - 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; - d.acceleration = [a, 0.0, 0.0]; - let m = model(1.0, nu, d, step(4.0 * h * h / nu, dx * dx, Some(3.0 / 16.0))); - let (setup, flow) = prepare(&m, &[]); - let Some((mut cpu, mut gpu)) = both(&setup, flow) else { return }; - let f = gpu.fields().unwrap(); - let mut err: f64 = 0.0; - for j in 0..n { - let y = (j as f64 + 0.5) * dx; - let exact = a / (2.0 * nu) * y * (h - y); - err = err.max((f.velocity[j][0] * setup.velocity_scale() - exact).abs() / umax); + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + let (h, nu, umax, n) = (1.0, 1.0 / 6.0, 0.1, 16usize); + let dx = h / n as f64; + let a = 8.0 * nu * umax / (h * h); + 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; + d.acceleration = [a, 0.0, 0.0]; + let m = model(1.0, nu, d, step(4.0 * h * h / nu, dx * dx, Some(3.0 / 16.0))); + let (setup, flow) = prepare(&m, &[]); + let Some((mut cpu, mut gpu)) = both(&setup, flow, streaming) else { continue }; + let f = gpu.fields().unwrap(); + let mut err: f64 = 0.0; + for j in 0..n { + let y = (j as f64 + 0.5) * dx; + let exact = a / (2.0 * nu) * y * (h - y); + err = err.max((f.velocity[j][0] * setup.velocity_scale() - exact).abs() / umax); + } + let (du, drho) = field_difference(&mut cpu, &mut gpu); + println!("Poiseuille on the GPU: max error {err:.2e} of the peak; CPU/GPU velocity {du:.2e}, density {drho:.2e}"); + assert!(err < 2e-5, "{err}"); + assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); } - let (du, drho) = field_difference(&mut cpu, &mut gpu); - println!("Poiseuille on the GPU: max error {err:.2e} of the peak; CPU/GPU velocity {du:.2e}, density {drho:.2e}"); - assert!(err < 2e-5, "{err}"); - assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); } #[test] fn couette_wall_forces_match_the_cpu() { - let (h, nu, uw, n) = (1.0, 0.1, 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; - d.faces.y_max = FaceCondition::MovingWall { velocity: [uw, 0.0, 0.0] }; - let m = model(1.0, nu, d, step(3.0 * h * h / nu, 0.5 * dx * dx / (3.0 * nu), None)); - let (setup, flow) = prepare(&m, &[]); - let Some((mut cpu, mut gpu)) = both(&setup, flow) else { return }; - let (fc, fg) = (Backend::forces(&mut cpu).unwrap(), gpu.forces().unwrap()); - let bottom = BoxFace::ALL.iter().position(|f| f.name() == "y_min").unwrap(); - let top = BoxFace::ALL.iter().position(|f| f.name() == "y_max").unwrap(); - let tau = nu * uw / h / setup.force_scale() * dx * dx; - println!("Couette: wall force CPU {:?} / GPU {:?} (expected ±{tau:.3e} per step)", fc[bottom], fg[bottom]); - // The tangential force is the shear: a cancellation of the links' exchanges to 1e-6 per link, - // summed from the stored deviations alone (adding the weights back in single precision left - // it 1e-4 of itself on an RTX 4090; 3e-6 on lavapipe this way). The normal force is the - // pressure. - for g in [bottom, top] { - assert!((fc[g][0] - fg[g][0]).abs() <= 5e-5 * tau.abs(), "group {g} shear: {} vs {}", fc[g][0], fg[g][0]); - assert!((fc[g][1] - fg[g][1]).abs() <= 1e-6 * fc[g][1].abs(), "group {g} normal: {} vs {}", fc[g][1], fg[g][1]); + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + let (h, nu, uw, n) = (1.0, 0.1, 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; + d.faces.y_max = FaceCondition::MovingWall { velocity: [uw, 0.0, 0.0] }; + let m = model(1.0, nu, d, step(3.0 * h * h / nu, 0.5 * dx * dx / (3.0 * nu), None)); + let (setup, flow) = prepare(&m, &[]); + let Some((mut cpu, mut gpu)) = both(&setup, flow, streaming) else { continue }; + let (fc, fg) = (Backend::forces(&mut cpu).unwrap(), gpu.forces().unwrap()); + let bottom = BoxFace::ALL.iter().position(|f| f.name() == "y_min").unwrap(); + let top = BoxFace::ALL.iter().position(|f| f.name() == "y_max").unwrap(); + let tau = nu * uw / h / setup.force_scale() * dx * dx; + println!("Couette: wall force CPU {:?} / GPU {:?} (expected ±{tau:.3e} per step)", fc[bottom], fg[bottom]); + // The tangential force is the shear: a cancellation of the links' exchanges to 1e-6 per link, + // summed from the stored deviations alone (adding the weights back in single precision left + // it 1e-4 of itself on an RTX 4090; 3e-6 on lavapipe this way). The normal force is the + // pressure. + for g in [bottom, top] { + println!("Couette group {g}: CPU {:.9e}, GPU {:.9e}, shear parity {:.9e}", fc[g][0], fg[g][0], (fc[g][0] - fg[g][0]).abs() / tau.abs()); + assert!((fc[g][0] - fg[g][0]).abs() <= 5e-5 * tau.abs(), "group {g} shear: {} vs {}", fc[g][0], fg[g][0]); + assert!((fc[g][1] - fg[g][1]).abs() <= 1e-6 * fc[g][1].abs(), "group {g} normal: {} vs {}", fc[g][1], fg[g][1]); + } + let (du, drho) = field_difference(&mut cpu, &mut gpu); + assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); } - let (du, drho) = field_difference(&mut cpu, &mut gpu); - assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); } #[test] fn taylor_green_decays_like_the_cpu_and_conserves_mass() { - let (nu, u0, l, n) = (0.1, 0.1, 2.0 * std::f64::consts::PI, 32usize); - let dx = l / n as f64; - let mut d = domain([0.0; 3], [l, l, dx], dx); - for f in BoxFace::ALL { - *d.faces.get_mut(f) = PERIODIC; + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + let (nu, u0, l, n) = (0.1, 0.1, 2.0 * std::f64::consts::PI, 32usize); + let dx = l / n as f64; + let mut d = domain([0.0; 3], [l, l, dx], dx); + for f in BoxFace::ALL { + *d.faces.get_mut(f) = PERIODIC; + } + let dt = 0.3 * dx * dx / (3.0 * nu); + let half_life = std::f64::consts::LN_2 / (2.0 * nu); + let m = model(1.0, nu, d, step(half_life, dt, None)); + let (setup, mut flow) = prepare(&m, &[]); + let grid = *flow.grid(); + let to_lattice = 1.0 / setup.velocity_scale(); + flow.set_equilibrium(|c| { + let [x, y, _] = grid.centre(c); + let p = -u0 * u0 / 4.0 * ((2.0 * x).cos() + (2.0 * y).cos()); + (setup.lattice_density(p), [-x.cos() * y.sin() * u0 * to_lattice, x.sin() * y.cos() * u0 * to_lattice, 0.0]) + }); + let Some(mut gpu) = gpu(&flow, streaming) else { continue }; + let mass0 = gpu.mass().unwrap(); + let mut cpu = flow; + cpu.advance(setup.steps).unwrap(); + gpu.advance(setup.steps).unwrap(); + let drift = (gpu.mass().unwrap() - mass0).abs() / mass0; + let (du, drho) = field_difference(&mut cpu, &mut gpu); + let f = gpu.fields().unwrap(); + let (mut num, mut den) = (0.0, 0.0); + for c in 0..grid.len() { + let [x, y, _] = grid.centre(c); + let mode = [-x.cos() * y.sin(), x.sin() * y.cos()]; + num += f.velocity[c][0] * mode[0] + f.velocity[c][1] * mode[1]; + den += mode[0] * mode[0] + mode[1] * mode[1]; + } + let amplitude = num / den * setup.velocity_scale(); + let exact = u0 * (-2.0 * nu * setup.steps as f64 * setup.time_step).exp(); + let err = (amplitude / exact - 1.0).abs(); + println!("Taylor–Green on the GPU: amplitude error {err:.2e}; CPU/GPU velocity {du:.2e}, density {drho:.2e}; mass drift {drift:.2e}"); + // The amplitude error is the scheme's at 32 cells (the CPU test shows second order). + assert!(err < 5e-3, "{err}"); + assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); + assert!(drift < 1e-6, "{drift}"); } - let dt = 0.3 * dx * dx / (3.0 * nu); - let half_life = std::f64::consts::LN_2 / (2.0 * nu); - let m = model(1.0, nu, d, step(half_life, dt, None)); - let (setup, mut flow) = prepare(&m, &[]); - let grid = *flow.grid(); - let to_lattice = 1.0 / setup.velocity_scale(); - flow.set_equilibrium(|c| { - let [x, y, _] = grid.centre(c); - let p = -u0 * u0 / 4.0 * ((2.0 * x).cos() + (2.0 * y).cos()); - (setup.lattice_density(p), [-x.cos() * y.sin() * u0 * to_lattice, x.sin() * y.cos() * u0 * to_lattice, 0.0]) - }); - let Some(mut gpu) = gpu(&flow) else { return }; - let mass0 = gpu.mass().unwrap(); - let mut cpu = flow; - cpu.advance(setup.steps).unwrap(); - gpu.advance(setup.steps).unwrap(); - let drift = (gpu.mass().unwrap() - mass0).abs() / mass0; - let (du, drho) = field_difference(&mut cpu, &mut gpu); - let f = gpu.fields().unwrap(); - let (mut num, mut den) = (0.0, 0.0); - for c in 0..grid.len() { - let [x, y, _] = grid.centre(c); - let mode = [-x.cos() * y.sin(), x.sin() * y.cos()]; - num += f.velocity[c][0] * mode[0] + f.velocity[c][1] * mode[1]; - den += mode[0] * mode[0] + mode[1] * mode[1]; - } - let amplitude = num / den * setup.velocity_scale(); - let exact = u0 * (-2.0 * nu * setup.steps as f64 * setup.time_step).exp(); - let err = (amplitude / exact - 1.0).abs(); - println!("Taylor–Green on the GPU: amplitude error {err:.2e}; CPU/GPU velocity {du:.2e}, density {drho:.2e}; mass drift {drift:.2e}"); - // The amplitude error is the scheme's at 32 cells (the CPU test shows second order). - assert!(err < 5e-3, "{err}"); - assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); - assert!(drift < 1e-6, "{drift}"); } #[test] fn flow_past_a_sphere_matches_the_cpu_fields_and_forces() { - let l = 24usize; - let side = l as f64; - let mut d = domain([0.0; 3], [side; 3], 1.0); - for f in BoxFace::ALL { - *d.faces.get_mut(f) = PERIODIC; + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + let l = 24usize; + let side = l as f64; + let mut d = domain([0.0; 3], [side; 3], 1.0); + for f in BoxFace::ALL { + *d.faces.get_mut(f) = PERIODIC; + } + d.acceleration = [1e-5, 0.0, 0.0]; + d.obstacles = vec!["Ball".into()]; + let ball = obstacle("Ball", sphere([side / 2.0; 3], 6.0, 3), &["skin"]); + let m = model( + 1.0, + 1.0 / 6.0, + d, + Procedure::Fluid { domain: "Flow".into(), end_time: None, steps: Some(300), time_step: Some(1.0), frames: 1, magic: None }, + ); + let (setup, flow) = prepare(&m, std::slice::from_ref(&ball)); + let Some((mut cpu, mut gpu)) = both(&setup, flow, streaming) else { continue }; + let (du, drho) = field_difference(&mut cpu, &mut gpu); + let (fc, fg) = (Backend::forces(&mut cpu).unwrap(), gpu.forces().unwrap()); + let skin = cpu.groups().iter().position(|g| g.name == "skin").unwrap(); + let scale = fc[skin].iter().map(|v| v.abs()).fold(0.0, f64::max); + let dforce = (0..3).map(|a| (fc[skin][a] - fg[skin][a]).abs()).fold(0.0, f64::max) / scale; + println!("Sphere: CPU/GPU velocity {du:.2e}, density {drho:.2e}; drag CPU {:?} GPU {:?} ({dforce:.2e})", fc[skin], fg[skin]); + assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); + assert!(dforce < 5e-5, "{dforce}"); } - d.acceleration = [1e-5, 0.0, 0.0]; - d.obstacles = vec!["Ball".into()]; - let ball = obstacle("Ball", sphere([side / 2.0; 3], 6.0, 3), &["skin"]); - let m = model( - 1.0, - 1.0 / 6.0, - d, - Procedure::Fluid { domain: "Flow".into(), end_time: None, steps: Some(300), time_step: Some(1.0), frames: 1, magic: None }, - ); - let (setup, flow) = prepare(&m, std::slice::from_ref(&ball)); - let Some((mut cpu, mut gpu)) = both(&setup, flow) else { return }; - let (du, drho) = field_difference(&mut cpu, &mut gpu); - let (fc, fg) = (Backend::forces(&mut cpu).unwrap(), gpu.forces().unwrap()); - let skin = cpu.groups().iter().position(|g| g.name == "skin").unwrap(); - let scale = fc[skin].iter().map(|v| v.abs()).fold(0.0, f64::max); - let dforce = (0..3).map(|a| (fc[skin][a] - fg[skin][a]).abs()).fold(0.0, f64::max) / scale; - println!("Sphere: CPU/GPU velocity {du:.2e}, density {drho:.2e}; drag CPU {:?} GPU {:?} ({dforce:.2e})", fc[skin], fg[skin]); - assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); - assert!(dforce < 5e-5, "{dforce}"); } #[test] fn zou_he_channel_matches_the_cpu() { - let (h, l, umax, nu) = (1.0, 4.0, 0.1, 0.05); - let dx = h / 10.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 m = model(1.0, nu, d, step(2.0 * l / umax, 0.05 * dx / umax, None)); - let (setup, flow) = prepare(&m, &[]); - let Some((mut cpu, mut gpu)) = both(&setup, flow) else { return }; - let (du, drho) = field_difference(&mut cpu, &mut gpu); - println!("Zou–He channel: CPU/GPU velocity {du:.2e}, density {drho:.2e}"); - assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + let (h, l, umax, nu) = (1.0, 4.0, 0.1, 0.05); + let dx = h / 10.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 m = model(1.0, nu, d, step(2.0 * l / umax, 0.05 * dx / umax, None)); + let (setup, flow) = prepare(&m, &[]); + let Some((mut cpu, mut gpu)) = both(&setup, flow, streaming) else { continue }; + let (du, drho) = field_difference(&mut cpu, &mut gpu); + println!("Zou–He channel: CPU/GPU velocity {du:.2e}, density {drho:.2e}"); + assert!(du < 2e-5 && drho < 1e-6, "{du} {drho}"); + } } // ---- Heat ---- @@ -252,42 +268,45 @@ fn heat_difference(cpu: &mut Flow, gpu: &mut GpuFlow) -> (f64, f64) { /// 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()); + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + 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, streaming) else { continue }; + 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}"); } - 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 @@ -295,26 +314,29 @@ fn conduction_between_plates_matches_the_analytic_profile_and_the_cpu() { /// 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}"); + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + 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, streaming) else { continue }; + 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 @@ -322,27 +344,108 @@ fn a_heated_channel_matches_the_cpu() { /// 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}"); + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + 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, streaming) else { continue }; + 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}"); + } +} + +/// Readbacks and force caches must describe the same state on both AA orientations, +/// including calls that cross a 64-step submission boundary. The CPU defines the result. +fn compare_incremental(mut cpu: Flow, heated: bool, streaming: Streaming) { + let Some(mut gpu) = gpu(&cpu, streaming) else { return }; + let mut total = 0; + for steps in [1, 1, 61, 1, 1, 64, 65] { + cpu.advance(steps).unwrap(); + gpu.advance(steps).unwrap(); + total += steps; + assert_eq!(gpu.steps_done(), total); + let (du, drho) = field_difference(&mut cpu, &mut gpu); + println!("incremental step {total}: velocity {du:.9e}, density {drho:.9e}"); + assert!(du < 2e-5 && drho < 1e-6, "step {total}: {du} {drho}"); + let mass = Backend::mass(&mut cpu).unwrap(); + assert!((gpu.mass().unwrap() - mass).abs() < 1e-6 * mass.abs()); + let a = Backend::forces(&mut cpu).unwrap(); + let b = gpu.forces().unwrap(); + let scale = a.iter().flatten().map(|v| v.abs()).fold(0.0, f64::max).max(1e-30); + let error = a.iter().flatten().zip(b.iter().flatten()).map(|(a, b)| (a - b).abs()).fold(0.0, f64::max); + assert!(error < 5e-5 * scale, "step {total}: force {error} / {scale}"); + gpu.advance(0).unwrap(); + assert_eq!(gpu.forces().unwrap(), b, "a zero-step advance preserves the force cache"); + if heated { + let (dt, dh) = heat_difference(&mut cpu, &mut gpu); + assert!(dt < 2e-5 && dh < 1e-4, "step {total}: heat {dt} {dh}"); + } + } +} + +#[test] +fn odd_even_and_chunk_readbacks_match_the_cpu_with_heat_open_and_slip_faces() { + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + let mut d = domain([0.0; 3], [16.0, 8.0, 1.0], 1.0); + d.faces.x_min = FaceCondition::Velocity { velocity: [0.01, 0.0, 0.0], profile: Profile::Uniform }; + d.faces.x_max = FaceCondition::Pressure { pressure: 0.0 }; + d.faces.y_max = FaceCondition::Slip {}; + d.faces.z_min = PERIODIC; + d.faces.z_max = PERIODIC; + let mut t = FluidThermal { initial_temperature: 0.5, ..FluidThermal::default() }; + t.faces.x_min = ThermalCondition::Temperature { value: 1.0 }; + t.faces.y_min = ThermalCondition::Flux { value: 0.001 }; + t.faces.y_max = ThermalCondition::Temperature { value: 0.0 }; + d.thermal = Some(t); + let m = hot_model( + 1.0, + 0.1, + 0.05, + d, + Procedure::Fluid { domain: "Flow".into(), end_time: None, steps: Some(194), time_step: Some(1.0), frames: 1, magic: None }, + ); + let (_, flow) = prepare(&m, &[]); + compare_incremental(flow, true, streaming); + } +} + +#[test] +fn odd_even_and_chunk_readbacks_match_the_cpu_around_an_obstacle() { + // Both layouts: the default depends on the device (see `Streaming`). + for streaming in [Streaming::InPlace, Streaming::TwoBuffers] { + let mut d = domain([0.0; 3], [12.0; 3], 1.0); + for face in BoxFace::ALL { + *d.faces.get_mut(face) = PERIODIC; + } + d.acceleration = [1e-5, 0.0, 0.0]; + d.obstacles = vec!["Ball".into()]; + let ball = obstacle("Ball", sphere([6.0; 3], 3.0, 2), &["skin"]); + let m = model( + 1.0, + 1.0 / 6.0, + d, + Procedure::Fluid { domain: "Flow".into(), end_time: None, steps: Some(194), time_step: Some(1.0), frames: 1, magic: None }, + ); + let (_, flow) = prepare(&m, &[ball]); + compare_incremental(flow, false, streaming); + } } diff --git a/docs/architecture.md b/docs/architecture.md index 3d6703a..cb72822 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, 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-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, one population buffer streamed in place (the AA pattern, `docs/lbm-gpu-performance.md`), 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 | diff --git a/docs/lbm-gpu-performance.md b/docs/lbm-gpu-performance.md new file mode 100644 index 0000000..938f605 --- /dev/null +++ b/docs/lbm-gpu-performance.md @@ -0,0 +1,49 @@ +# GPU lattice Boltzmann: streaming layout and performance + +The momentum populations stream in one of two layouts, chosen per flow with +`fieldcraft_lbm_gpu::Options::streaming`: + +- **In place** (the AA pattern: Bailey, Myre, Walsh, Lilja and Saar, "Accelerating lattice + Boltzmann fluid flow simulations using graphics processors", ICPP 2009; Wittmann, Zeiser, Hager + and Wellein, Comput. Math. Appl. 65 (2013) 924–935, §6). One population buffer. Even steps read + population `q` of cell `x` at `(q, x)` and write the collided one at `(opposite(q), x)`; odd steps + read it at `(opposite(q), x − c_q)` and write it at `(q, x + c_q)`. In both, a cell writes exactly + the addresses it read, so cells never race, and streaming happens inside the collision: 152 + population bytes per cell update. +- **Two buffers** (A–B). Collide in place, then pull into the other buffer: 304 bytes per update. + +The wall, reflection and open-face rules are the same per-link rules in both. In place, a gather +pass reads every link's post-collision source into a scratch value before a scatter pass writes, +so no link overwrites another's source; the momentum exchange is taken in the gather of the step +whose forces are read back (a submission's last). The temperature lattice is two-buffer in both. + +`Streaming::Auto`, what a solve uses, streams in place: it was faster on every device measured +(the ignored throughput test, million cell updates per second, means of three runs; the 4090 row +is from the first version of this layout, before its addressing was made cheaper): + +| device | case | in place | two buffers | +|---|---|---|---| +| RTX 4090, Vulkan | periodic 128³ | 3897 | 2578 | +| RTX 4090, Vulkan | sphere 128³ | 3572 | 2441 | +| GB10 (unified memory), Vulkan | periodic 128³ | 1224 | 656 | +| GB10 (unified memory), Vulkan | sphere 128³ | 1146 | 630 | +| M4 Max, Metal | periodic 128³ | 1364 | 937 | +| M4 Max, Metal | sphere 128³ | 1256 | 872 | + +On the 4090 and the GB10 the two-buffer step is bound by memory bandwidth (~0.8 TB/s and ~0.2 of +the GB10's 0.27), so halving the bytes pays directly; on the M4 Max it gains once the in-place +addressing is cheap (a cell's coordinates once, neighbours wrapped by compare and add). With +128³ cells in place also needs one 152 MiB population buffer less. 64³ fits in cache, where the +two layouts are closer and vary more between runs. + +## Measuring + +``` +cargo test --release -p fieldcraft-lbm-gpu --test throughput -- --ignored --nocapture +``` + +runs each case in both layouts on the device `FIELDCRAFT_GPU` picks, marks the one `Auto` chooses +and prints how long creating the pipelines took (on Metal this reads near zero: the driver +compiles lazily or from its cache, so the cost lands in the first dispatches instead). With `FIELDCRAFT_LBM_GPU_PROFILE=1` it also prints +each kernel's GPU time (`Options::profile`, timestamp queries; the throughput is then not +representative, as every dispatch gets its own pass).