diff --git a/AGENTS.md b/AGENTS.md index 2b5fb7c..ca61b31 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -79,7 +79,7 @@ People trust ModelCraft with their models; a crash loses their work. **This outr JSON `run`) + tests. The UI, scripts, CLI and control channel all reach the same commands. - **Programmatic calls never open dialogs.** - **Layering** is enforced by `cargo xtask layers` (`docs/architecture.md`): L0 `model` → L1 - `geometry`, `mesh`, `fem`, `decks`, `render` → L2 solvers (`solver`; `explicit`, `lbm` later) → + `geometry`, `mesh`, `fem`, `decks`, `render` → L2 solvers (`solver`, `lbm`; `explicit` later) → L3 `fsi` (later) → L4 `doc` → L5 `io` → L6 `engine` → L7 `ui-egui`, `mcp` → apps. Nothing below L7 depends on egui/eframe/winit/wgpu/rfd. **Only `crates/geometry` may use `solidcraft-*`.** diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 1db11fa..ff0d3c1 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -26,8 +26,6 @@ is committed and listed here with its exports. | assets/icons/file.recover.svg, assets/icons/file.recover.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/file.recover.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/file.save_as.svg, assets/icons/file.save_as.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/file.save_as.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/file.save.svg, assets/icons/file.save.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/file.save.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | -| assets/icons/flow.boundary.svg, assets/icons/flow.boundary.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/flow.boundary.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | -| assets/icons/flow.domain.svg, assets/icons/flow.domain.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/flow.domain.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/geometry.bodies.svg, assets/icons/geometry.bodies.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/geometry.bodies.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/geometry.datum.svg, assets/icons/geometry.datum.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/geometry.datum.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/geometry.import.mesh.svg, assets/icons/geometry.import.mesh.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/geometry.import.mesh.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | @@ -76,6 +74,8 @@ is committed and listed here with its exports. | assets/icons/model.boundary_condition.create.svg, assets/icons/model.boundary_condition.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.boundary_condition.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/model.check.svg, assets/icons/model.check.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.check.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/model.coupling.create.svg, assets/icons/model.coupling.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.coupling.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | +| assets/icons/model.fluid_domain.boundary.svg, assets/icons/model.fluid_domain.boundary.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.fluid_domain.boundary.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | +| assets/icons/model.fluid_domain.create.svg, assets/icons/model.fluid_domain.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.fluid_domain.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/model.info.svg, assets/icons/model.info.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.info.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/model.load.create.svg, assets/icons/model.load.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.load.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/model.material.create.svg, assets/icons/model.material.create.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/model.material.create.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | @@ -111,7 +111,7 @@ is committed and listed here with its exports. | assets/icons/results.xy.svg, assets/icons/results.xy.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/results.xy.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/solve.blast.svg, assets/icons/solve.blast.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/solve.blast.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/solve.explicit.svg, assets/icons/solve.explicit.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/solve.explicit.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | -| assets/icons/solve.flow.svg, assets/icons/solve.flow.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/solve.flow.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | +| assets/icons/solve.fluid.svg, assets/icons/solve.fluid.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/solve.fluid.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/solve.frequency.svg, assets/icons/solve.frequency.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/solve.frequency.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/solve.fsi.svg, assets/icons/solve.fsi.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/solve.fsi.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | | assets/icons/solve.static.svg, assets/icons/solve.static.vectorcraft | ModelCraft contributors | Original toolbar icon drawn with VectorCraft from `assets/icons/src/solve.static.json` (`cargo xtask icons`) | MIT OR Apache-2.0 | diff --git a/Cargo.lock b/Cargo.lock index f564384..9ebcd3e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1884,6 +1884,7 @@ dependencies = [ "modelcraft-fem", "modelcraft-geometry", "modelcraft-io", + "modelcraft-lbm", "modelcraft-mesh", "modelcraft-model", "modelcraft-solver", @@ -1940,6 +1941,15 @@ dependencies = [ "zip", ] +[[package]] +name = "modelcraft-lbm" +version = "0.1.0" +dependencies = [ + "modelcraft-mesh", + "modelcraft-model", + "thiserror 2.0.21", +] + [[package]] name = "modelcraft-mesh" version = "0.1.0" diff --git a/Cargo.toml b/Cargo.toml index cd9bf5c..4604860 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -20,6 +20,7 @@ modelcraft-io = { path = "crates/io" } modelcraft-fem = { path = "crates/fem" } modelcraft-solver = { path = "crates/solver" } modelcraft-explicit = { path = "crates/explicit" } +modelcraft-lbm = { path = "crates/lbm" } modelcraft-engine = { path = "crates/engine" } modelcraft-testkit = { path = "crates/testkit" } modelcraft-render = { path = "crates/render" } @@ -75,6 +76,11 @@ opt-level = 2 [profile.dev] opt-level = 1 +# The lattice Boltzmann kernels are compute loops over millions of cells; their verification +# runs (tests) need them optimised. +[profile.dev.package.modelcraft-lbm] +opt-level = 3 + # Tests and lints (`cargo xtask ci`, GitHub CI). Debug assertions are off because truck's # debug-only checks reject valid geometry (SolidCraft builds its tests the same way). [profile.ci] @@ -88,6 +94,9 @@ incremental = false [profile.ci.package."*"] opt-level = 2 +[profile.ci.package.modelcraft-lbm] +opt-level = 3 + [profile.release] opt-level = 3 lto = "thin" diff --git a/ROADMAP.md b/ROADMAP.md index 632869a..ef21b64 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -30,7 +30,7 @@ hex meshes today and it gives the best elements where stresses are read. | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned | | M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `modelcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (objective rates, J2, Johnson–Cook, damage, erosion); contact (global search, penalty with friction, rigid walls); runs on hex and tet meshes | started | | M11 Equations of state and shells | Polynomial, ideal-gas, Mie–Grüneisen, JWL with programmed burn, tabulated; energy coupling; shock verification (Rankine–Hugoniot, Sod, Noh, Sedov); Belytschko–Lin–Tsay shells with self-contact on quad meshes | planned | -| M12 CFD: lattice Boltzmann | D3Q19/D3Q27 on grids voxelised from the face-tagged surface; TRT, regularised and cumulant collision; interpolated walls, inlets and outlets; LES; octree refinement shared with the hex mesher; wgpu backend (deterministic per device, checked against the CPU) | planned | +| M12 CFD: lattice Boltzmann | D3Q19/D3Q27 on grids voxelised from the face-tagged surface; TRT, regularised and cumulant collision; interpolated walls, inlets and outlets; LES; octree refinement shared with the hex mesher; wgpu backend (deterministic per device, checked against the CPU) (CPU core: done — `modelcraft-lbm` with D3Q19 TRT and Guo forcing, halfway bounce-back and moving walls, slip and periodic faces, Zou–He velocity and pressure faces, exact voxelisation of obstacle bodies, momentum-exchange forces per face; fluid materials, fluid domains and fluid steps in the model; `model.fluid_domain.*`, `solve.fluid`; verified by Poiseuille channel and pipe, Couette, Taylor–Green decay, sphere-array Stokes drag and mass conservation; Schäfer–Turek 2D-1 as a benchmark) | started | | M13 Fluid–structure interaction | Immersed-boundary coupling of the lattice Boltzmann fluid with the explicit or implicit structure; interface force transfer; sub-cycling; Aitken then interface quasi-Newton coupling; free-surface flow | planned | | M14 Blast and impact fluids | EOS-driven fluids in the explicit solver: ALE multi-material or SPH, with fluid–structure contact | planned | | Continuous: verification and validation | Every capability checked against theory (closed forms), experiments (published measurements) and published benchmarks (NAFEMS, vendor benchmark and verification manuals such as Abaqus's), each case our own model in `tests/oracle/` with its cited reference and tolerance; `cargo xtask oracle` runs them and writes docs/oracle.md; once a case matches, our result is frozen as a fixture so fidelity cannot regress (docs/verification.md) | started | diff --git a/assets/icons/flow.boundary.vectorcraft b/assets/icons/flow.boundary.vectorcraft deleted file mode 100644 index 930f96c..0000000 Binary files a/assets/icons/flow.boundary.vectorcraft and /dev/null differ diff --git a/assets/icons/flow.domain.vectorcraft b/assets/icons/flow.domain.vectorcraft deleted file mode 100644 index 60bc6ad..0000000 Binary files a/assets/icons/flow.domain.vectorcraft and /dev/null differ diff --git a/assets/icons/flow.boundary.svg b/assets/icons/model.fluid_domain.boundary.svg similarity index 94% rename from assets/icons/flow.boundary.svg rename to assets/icons/model.fluid_domain.boundary.svg index 62755a1..73a9ed0 100644 --- a/assets/icons/flow.boundary.svg +++ b/assets/icons/model.fluid_domain.boundary.svg @@ -1,6 +1,6 @@ - flow.boundary + model.fluid_domain.boundary diff --git a/assets/icons/model.fluid_domain.boundary.vectorcraft b/assets/icons/model.fluid_domain.boundary.vectorcraft new file mode 100644 index 0000000..ed7a8bc Binary files /dev/null and b/assets/icons/model.fluid_domain.boundary.vectorcraft differ diff --git a/assets/icons/flow.domain.svg b/assets/icons/model.fluid_domain.create.svg similarity index 94% rename from assets/icons/flow.domain.svg rename to assets/icons/model.fluid_domain.create.svg index 7c362ae..f907cf3 100644 --- a/assets/icons/flow.domain.svg +++ b/assets/icons/model.fluid_domain.create.svg @@ -1,6 +1,6 @@ - flow.domain + model.fluid_domain.create diff --git a/assets/icons/model.fluid_domain.create.vectorcraft b/assets/icons/model.fluid_domain.create.vectorcraft new file mode 100644 index 0000000..58d931f Binary files /dev/null and b/assets/icons/model.fluid_domain.create.vectorcraft differ diff --git a/assets/icons/solve.flow.vectorcraft b/assets/icons/solve.flow.vectorcraft deleted file mode 100644 index f4ae396..0000000 Binary files a/assets/icons/solve.flow.vectorcraft and /dev/null differ diff --git a/assets/icons/solve.flow.svg b/assets/icons/solve.fluid.svg similarity index 96% rename from assets/icons/solve.flow.svg rename to assets/icons/solve.fluid.svg index 2de5b8b..3dcfb00 100644 --- a/assets/icons/solve.flow.svg +++ b/assets/icons/solve.fluid.svg @@ -1,6 +1,6 @@ - solve.flow + solve.fluid diff --git a/assets/icons/solve.fluid.vectorcraft b/assets/icons/solve.fluid.vectorcraft new file mode 100644 index 0000000..9b0f537 Binary files /dev/null and b/assets/icons/solve.fluid.vectorcraft differ diff --git a/assets/icons/src/flow.boundary.json b/assets/icons/src/model.fluid_domain.boundary.json similarity index 100% rename from assets/icons/src/flow.boundary.json rename to assets/icons/src/model.fluid_domain.boundary.json diff --git a/assets/icons/src/flow.domain.json b/assets/icons/src/model.fluid_domain.create.json similarity index 100% rename from assets/icons/src/flow.domain.json rename to assets/icons/src/model.fluid_domain.create.json diff --git a/assets/icons/src/solve.flow.json b/assets/icons/src/solve.fluid.json similarity index 100% rename from assets/icons/src/solve.flow.json rename to assets/icons/src/solve.fluid.json diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 835cf19..004173b 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -104,6 +104,7 @@ pub fn write(model: &Model) -> Result { /// [`write`] for a model already resolved. pub fn write_resolved(model: &Model, r: &Resolved) -> Result { + crate::refuse_fluid(FORMAT, model)?; let by_type = elements_by_type(model, r)?; check_expressible(model)?; let mut d = Deck::default(); @@ -757,6 +758,7 @@ impl Deck { for step in steps { self.comment(&format!("Step {}", step.name)); match step.procedure { + Procedure::Fluid { .. } => return Err(crate::fluid_step(FORMAT, &step.name)), Procedure::Static {} => { self.line(&format!("*STEP, NAME={}, NLGEOM=NO", step.name)); self.line("*STATIC"); diff --git a/crates/decks/src/lib.rs b/crates/decks/src/lib.rs index cb3ba32..bbad58e 100644 --- a/crates/decks/src/lib.rs +++ b/crates/decks/src/lib.rs @@ -83,6 +83,19 @@ pub fn with_steps(model: &Model, names: &[String]) -> Result { Ok(out) } +/// A fluid step: solved by ModelCraft's lattice Boltzmann solver, never written to a deck. +pub(crate) fn fluid_step(format: &'static str, step: &str) -> Error { + Error::Unsupported { format, what: format!("fluid step `{step}` (solve it with solve.fluid; write the structural steps only)") } +} + +/// Every writer's first check: the model has no fluid step. +pub(crate) fn refuse_fluid(format: &'static str, model: &Model) -> Result<()> { + match model.steps().iter().find(|s| s.procedure.fluid_domain().is_some()) { + Some(step) => Err(fluid_step(format, &step.name)), + None => Ok(()), + } +} + /// Every deck format, one row per writer. pub const FORMATS: &[Format] = &[ Format { name: "abaqus", extension: "inp", writer: abaqus::write_resolved }, diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 1f2c8ed..5c27750 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -126,6 +126,7 @@ pub(crate) fn write_resolved(model: &Model, r: &Resolved) -> Result { } fn write_with(model: &Model, r: &Resolved) -> Result { + crate::refuse_fluid(FORMAT, model)?; let estimate = (model.nodes().len() + model.reference_points().len()) * 130 + model.elements().len() * 260; if estimate > MAX_DECK_BYTES { return Err(Error::Limit(format!("the deck would be about {estimate} bytes; at most {MAX_DECK_BYTES} are written"))); @@ -566,6 +567,7 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste d.header(&["imflag", "dt0"]); d.card([Field::from(1), Field::from(1.0)]); match step.procedure { + Procedure::Fluid { .. } => return Err(crate::fluid_step(FORMAT, &step.name)), Procedure::Static {} => { d.line("*CONTROL_IMPLICIT_SOLUTION"); d.header(&["nsolvr"]); diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index 7fe8226..5b86fb9 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -60,6 +60,7 @@ pub fn write(model: &Model) -> Result { /// [`write`] for a model already resolved. pub fn write_resolved(model: &Model, r: &Resolved) -> Result { + crate::refuse_fluid(FORMAT, model)?; let mut deck = Deck::new(model, r, model.steps().iter().collect()); deck.plan()?; deck.emit() diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs index 1057c2e..3e4a380 100644 --- a/crates/decks/tests/abaqus.rs +++ b/crates/decks/tests/abaqus.rs @@ -15,7 +15,7 @@ use modelcraft_model::{ // Models fn steel() -> Material { - Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9) } + Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None } } fn faces(list: &[(u32, u8)]) -> Set { @@ -91,7 +91,7 @@ fn plate() -> Model { m.add_set("top", faces(&[(1, 1), (2, 1), (3, 1), (4, 1)])).unwrap(); m.add_set("edge", Set { nodes: vec![1, 4, 7], ..Set::default() }).unwrap(); m.add_set("mid", Set { nodes: vec![2, 5, 8], ..Set::default() }).unwrap(); - let aluminium = Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9) }; + let aluminium = Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9), fluid: None }; m.add_material("Aluminium", aluminium).unwrap(); let kind = SectionKind::Shell { thickness: 2.0, offset: 0.5, points: 5 }; let formulation = Formulation { integration: Integration::Reduced, ..Formulation::default() }; @@ -523,7 +523,8 @@ fn largest_labels_longest_numbers_and_names_fit_a_line() { m.set_element(modelcraft_model::MAX_LABEL, ElementKind::Hex20, (0..20).map(|i| base + i).collect()).unwrap(); let long = format!("A{}", "x".repeat(modelcraft_model::MAX_NAME_LEN - 1)); m.add_set(&long, Set { elements: vec![modelcraft_model::MAX_LABEL], nodes: (0..20).map(|i| base + i).collect(), ..Set::default() }).unwrap(); - m.add_material(&long, Material { elastic: Some(Elastic::Isotropic { young: ugly.abs(), poisson: 0.3 }), density: Some(1.0) }).unwrap(); + m.add_material(&long, Material { elastic: Some(Elastic::Isotropic { young: ugly.abs(), poisson: 0.3 }), density: Some(1.0), fluid: None }) + .unwrap(); m.add_section(&long, Section { set: long.clone(), material: long.clone(), kind: SectionKind::Solid {}, formulation: Formulation::default() }) .unwrap(); m.add_reference_point(&long, ReferencePoint { xyz: [ugly; 3], label: 0 }).unwrap(); @@ -639,7 +640,7 @@ fn what_abaqus_cannot_express_is_named() { let mut ortho = plate(); let carbon = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] }; - ortho.replace_material("Aluminium", Material { elastic: Some(carbon), density: Some(1.6e-9) }).unwrap(); + ortho.replace_material("Aluminium", Material { elastic: Some(carbon), density: Some(1.6e-9), fluid: None }).unwrap(); assert!(unsupported(&ortho).contains("orthotropic material `Aluminium` in a shell section")); let mut loose = frequency(); @@ -657,7 +658,7 @@ fn what_abaqus_cannot_express_is_named() { fn orthotropic_solids_beams_and_rods_are_written() { let mut m = cantilever(); let carbon = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] }; - m.replace_material("Steel", Material { elastic: Some(carbon), density: None }).unwrap(); + m.replace_material("Steel", Material { elastic: Some(carbon), density: None, fluid: None }).unwrap(); let deck = abaqus::write(&m).unwrap(); assert!(deck.contains("*ELASTIC, TYPE=ENGINEERING CONSTANTS\n140000.0, 10000.0, 10000.0, 0.3, 0.3, 0.4, 5000.0, 5000.0\n3800.0\n"), "{deck}"); diff --git a/crates/decks/tests/lsdyna.rs b/crates/decks/tests/lsdyna.rs index ae82486..62dc6d9 100644 --- a/crates/decks/tests/lsdyna.rs +++ b/crates/decks/tests/lsdyna.rs @@ -13,7 +13,7 @@ use modelcraft_model::*; // Our own models. fn steel() -> Material { - Material { elastic: Some(Elastic::Isotropic { young: 210000.0, poisson: 0.3 }), density: Some(7.85e-9) } + Material { elastic: Some(Elastic::Isotropic { young: 210000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None } } fn step(name: &str, procedure: Procedure, bcs: &[&str], loads: &[&str]) -> Step { @@ -72,7 +72,7 @@ fn cantilever_modes() -> Model { let mut m = Model::default(); cantilever_mesh(&mut m); let ortho = Elastic::EngineeringConstants { young: [140000.0, 10000.0, 10000.0], poisson: [0.3, 0.3, 0.45], shear: [5000.0, 5000.0, 3500.0] }; - m.add_material("Ply", Material { elastic: Some(ortho), density: Some(1.6e-9) }).unwrap(); + m.add_material("Ply", Material { elastic: Some(ortho), density: Some(1.6e-9), fluid: None }).unwrap(); let f = Formulation { incompatible_modes: true, ..Formulation::default() }; m.add_section("beam", Section { set: "all".into(), material: "Ply".into(), kind: SectionKind::Solid {}, formulation: f }).unwrap(); m.add_reference_point("RP", ReferencePoint { xyz: [5.0, 0.5, 0.5], label: 0 }).unwrap(); @@ -155,7 +155,8 @@ fn solids() -> Model { let h: Vec = h.iter().map(|p| node(&mut m, *p)).collect(); m.set_element(4, ElementKind::Hex20, h).unwrap(); - m.add_material("Aluminium", Material { elastic: Some(Elastic::Isotropic { young: 70e9, poisson: 0.33 }), density: Some(2700.0) }).unwrap(); + m.add_material("Aluminium", Material { elastic: Some(Elastic::Isotropic { young: 70e9, poisson: 0.33 }), density: Some(2700.0), fluid: None }) + .unwrap(); for (e, name) in [(1, "tet4"), (2, "tet10"), (3, "wedge"), (4, "hex20")] { m.add_set(name, Set { elements: vec![e], ..Set::default() }).unwrap(); let f = Formulation { hybrid: e == 1, ..Formulation::default() }; @@ -657,7 +658,7 @@ fn hostile_numbers_stay_within_their_fields() { m.set_node(2, [f64::MIN_POSITIVE, -f64::MIN_POSITIVE, 5e-324]).unwrap(); m.replace_material( "Steel", - Material { elastic: Some(Elastic::Isotropic { young: 1e300, poisson: 0.499_999_999_999_999_9 }), density: Some(1e-300) }, + Material { elastic: Some(Elastic::Isotropic { young: 1e300, poisson: 0.499_999_999_999_999_9 }), density: Some(1e-300), fluid: None }, ) .unwrap(); m.replace_load( diff --git a/crates/decks/tests/nastran.rs b/crates/decks/tests/nastran.rs index 3b8017b..b2856ff 100644 --- a/crates/decks/tests/nastran.rs +++ b/crates/decks/tests/nastran.rs @@ -131,7 +131,7 @@ impl Read { // Our own models. fn steel() -> Material { - Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9) } + Material { elastic: Some(Elastic::Isotropic { young: 210_000.0, poisson: 0.3 }), density: Some(7.85e-9), fluid: None } } /// A bar of `n` unit hex8 along X; node (i, y, z) = 1 + 4i + 2y + z. @@ -188,7 +188,11 @@ fn plate() -> Model { m.add_set("all", Set { elements: vec![1, 2, 3, 4], ..Set::default() }).unwrap(); m.add_set("edge", Set { nodes: vec![1, 4, 7], ..Set::default() }).unwrap(); m.add_set("top", Set { faces: (1..=4).map(|element| ElementFace { element, face: 1 }).collect(), ..Set::default() }).unwrap(); - m.add_material("Aluminium", Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9) }).unwrap(); + m.add_material( + "Aluminium", + Material { elastic: Some(Elastic::Isotropic { young: 70_000.0, poisson: 0.33 }), density: Some(2.7e-9), fluid: None }, + ) + .unwrap(); let shell = SectionKind::Shell { thickness: 2.0, offset: 0.0, points: 5 }; m.add_section("skin", Section { set: "all".into(), material: "Aluminium".into(), kind: shell, formulation: Formulation::default() }).unwrap(); m.add_boundary_condition("clamp", BoundaryCondition::Fixed { target: Target::Set("edge".into()) }).unwrap(); @@ -517,6 +521,7 @@ fn beams_rods_orthotropy_couplings_and_prescribed_motion() { let ortho = Material { elastic: Some(Elastic::EngineeringConstants { young: [140e3, 10e3, 12e3], poisson: [0.3, 0.28, 0.4], shear: [5e3, 4.5e3, 3.8e3] }), density: Some(1.6e-9), + fluid: None, }; m.add_material("Carbon", ortho.clone()).unwrap(); m.replace_section( @@ -671,7 +676,7 @@ fn what_nastran_cannot_express_is_named() { // Orthotropic beams, tilted orthotropic shells, one-point shells. let mut m = plate(); let ortho = Elastic::EngineeringConstants { young: [140e3, 10e3, 10e3], poisson: [0.3, 0.3, 0.4], shear: [5e3, 5e3, 3.8e3] }; - m.replace_material("Aluminium", Material { elastic: Some(ortho), density: Some(1.6e-9) }).unwrap(); + m.replace_material("Aluminium", Material { elastic: Some(ortho), density: Some(1.6e-9), fluid: None }).unwrap(); m.set_node(9, [10.0, 10.0, 0.5]).unwrap(); assert!(unsupported(&m).contains("planes of constant Z")); m.set_node(9, [10.0, 10.0, 0.0]).unwrap(); @@ -690,7 +695,7 @@ fn what_nastran_cannot_express_is_named() { let mut m = plate(); m.set_element(9, ElementKind::Rod2, vec![1, 9]).unwrap(); m.add_set("bars", Set { elements: vec![9], ..Set::default() }).unwrap(); - m.add_material("Fibre", Material { elastic: Some(ortho), density: Some(1.6e-9) }).unwrap(); + m.add_material("Fibre", Material { elastic: Some(ortho), density: Some(1.6e-9), fluid: None }).unwrap(); m.add_section( "bars", Section { set: "bars".into(), material: "Fibre".into(), kind: SectionKind::Rod { area: 1.0 }, formulation: Formulation::default() }, diff --git a/crates/engine/Cargo.toml b/crates/engine/Cargo.toml index 378f51e..dc95539 100644 --- a/crates/engine/Cargo.toml +++ b/crates/engine/Cargo.toml @@ -17,6 +17,7 @@ modelcraft-decks = { workspace = true } modelcraft-solver = { workspace = true } modelcraft-explicit = { workspace = true } modelcraft-fem = { workspace = true } +modelcraft-lbm = { workspace = true } serde = { workspace = true } serde_json = { workspace = true } thiserror = { workspace = true } diff --git a/crates/engine/src/catalogue.rs b/crates/engine/src/catalogue.rs index 5362fa5..c9fb222 100644 --- a/crates/engine/src/catalogue.rs +++ b/crates/engine/src/catalogue.rs @@ -261,8 +261,21 @@ pub static CATALOGUE: &[Capability] = &[ "Force, moment, pressure, traction, gravity or body force", &["model.load.edit", "model.load.delete"], ), - soon("flow.domain", "Fluid Domain", Loads, "Flow", "M12", "The fluid grid voxelised from the face-tagged surface"), - soon("flow.boundary", "Inlet & Outlet", Loads, "Flow", "M12", "Velocity inlets, pressure outlets and walls"), + ready_more( + "model.fluid_domain.create", + "Fluid Domain", + Loads, + "Flow", + "The box of fluid cells, its fluid, and the obstacle bodies voxelised from their face-tagged surfaces", + &["model.fluid_domain.edit", "model.fluid_domain.delete"], + ), + ready( + "model.fluid_domain.boundary", + "Inlet & Outlet", + Loads, + "Flow", + "Velocity inlets, pressure outlets, walls, moving and slip walls, periodic faces", + ), // ---- Interactions ---- ready_more( "model.reference_point.create", @@ -306,7 +319,7 @@ pub static CATALOGUE: &[Capability] = &[ "Transient dynamics of an explicit step by central differences (contact and erosion coming in M10)", ), soon("solve.blast", "Blast & Impact", Solve, "Explicit", "M14", "EOS-driven fluids with fluid–structure contact (ALE or SPH)"), - soon("solve.flow", "Lattice Boltzmann", Solve, "Fluids", "M12", "Lattice Boltzmann flow on the CPU or the GPU"), + ready("solve.fluid", "Lattice Boltzmann", Solve, "Fluids", "Lattice Boltzmann flow of a fluid step on the CPU, with wall forces"), soon("solve.fsi", "Fluid–Structure", Solve, "Fluids", "M13", "Coupled fluid and structure solve"), // ---- Results ---- ready_more("results.deformed", "Deformed", Results, "Display", "The deformed shape of a solved step, to scale", &["results.clear"]), diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index ed2b3c6..8541962 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -1,11 +1,14 @@ //! Model definition commands: create, edit and delete sets, materials, sections, reference -//! points, couplings, boundary conditions, loads and steps, and the step activity table. +//! points, couplings, boundary conditions, loads, fluid domains and steps, and the step activity +//! table. //! //! Each definition's parameters are its JSON form in the model (the same form a `.modelcraft` //! file stores), plus its `name`. `create` fails when the name is taken, `edit` replaces an //! existing definition, `delete` refuses while another definition refers to it. -use modelcraft_model::{BoundaryCondition, Coupling, Load, Material, Model, ReferencePoint, Section, Set, Step}; +use modelcraft_model::{ + BoundaryCondition, Coupling, FaceCondition, FluidBoundary, FluidDomain, Load, Material, Model, ReferencePoint, Section, Set, Step, +}; use serde::de::DeserializeOwned; use serde_json::{Value, json}; @@ -77,7 +80,7 @@ definition_commands!( Material, "material", "Material", - r#"{"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3} | {"type": "engineering_constants", "young": [E1, E2, E3], "poisson": [ν12, ν13, ν23], "shear": [G12, G13, G23]}, "density": 7.85e-9} — in the model's unit system (model.units), never converted"#, + r#"{"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3} | {"type": "engineering_constants", "young": [E1, E2, E3], "poisson": [ν12, ν13, ν23], "shear": [G12, G13, G23]}, "density": 7.85e-9, "fluid": {"type": "newtonian", "viscosity": 1.0e-9} (dynamic viscosity; a fluid needs a density)} — in the model's unit system (model.units), never converted"#, add_material, replace_material, remove_material, @@ -156,10 +159,40 @@ definition_commands!( LOAD_DELETE ); +definition_commands!( + "fluid domain", + FluidDomain, + "fluid_domain", + "Fluid Domain", + r#"{"name": "Flow", "material": "Water" (a density and fluid behaviour), "min": [x, y, z], "max": [x, y, z], "spacing": 1.0 (cubic cells; every extent a whole number of spacings), "faces": {"x_min": {"type": "velocity", "velocity": [u, v, w], "profile": "uniform|parabolic"}, "x_max": {"type": "pressure", "pressure": 0}, "y_min": {"type": "wall"} | {"type": "moving_wall", "velocity": [u, v, w] (in the face's plane)} | {"type": "slip"} | {"type": "periodic"} (in opposite pairs), …} (default walls), "obstacles": ["Body1"] (solid bodies of the design or imported surfaces, by name), "surfaces": [{"face": {"body": "Body1", "face": "F2:top"}, "condition": {"type": "moving_wall", "velocity": [u, v, w]}}] (obstacle faces are walls unless listed), "initial_velocity": [u, v, w], "acceleration": [x, y, z] (gravity or a driving force per mass), "tolerance": 0.1 (obstacle triangulation; default a tenth of the spacing)}"#, + add_fluid_domain, + replace_fluid_domain, + remove_fluid_domain, + FLUID_DOMAIN_CREATE, + FLUID_DOMAIN_EDIT, + FLUID_DOMAIN_DELETE +); + +pub const FLUID_BOUNDARY: CommandSpec = CommandSpec { + id: "model.fluid_domain.boundary", + label: "Inlet, Outlet or Wall", + params: r#"{"domain": "Flow", "at": "x_min" (a box face) | {"body": "Body1", "face": "F2:top"} (an obstacle face), "condition": {"type": "velocity", "velocity": [u, v, w], "profile": "uniform|parabolic"} | {"type": "pressure", "pressure": 0} | {"type": "wall"} | {"type": "moving_wall", "velocity": [u, v, w]} | {"type": "slip"} | {"type": "periodic"}} — sets one boundary's condition (obstacle faces take walls and moving walls)"#, + kind: CommandKind::Model, + run: fluid_boundary, +}; + +fn fluid_boundary(s: &mut Session, p: &Value) -> Result { + let domain: &str = req(p, "domain")?; + let at: FluidBoundary = req(p, "at")?; + let condition: FaceCondition = req(p, "condition")?; + s.model.set_fluid_boundary(domain, at, condition)?; + Ok(json!({"domain": domain})) +} + pub const STEP_CREATE: CommandSpec = CommandSpec { id: "model.step.create", label: "New Step", - params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal"}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, + params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "fluid", "domain": "Flow", "end_time": 1.0 | "steps": 1000, "time_step": (default: chosen for Mach 0.087 and τ ≤ 1), "frames": 10, "magic": 0.1875 (the two-relaxation-time Λ)} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal"}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, kind: CommandKind::Model, run: step_create, }; @@ -207,7 +240,7 @@ fn step_edit(s: &mut Session, p: &Value) -> Result { let n = name(p)?; let mut step = s.model.step(n).cloned().ok_or_else(|| EngineError::Params(format!("no step `{n}`")))?; if let Some(procedure) = opt(p, "procedure")? { - step.procedure = procedure; + step.set_procedure(procedure); } if let Some(o) = opt(p, "output")? { step.output = o; diff --git a/crates/engine/src/cmd/mod.rs b/crates/engine/src/cmd/mod.rs index aa4bdfb..8ed73e4 100644 --- a/crates/engine/src/cmd/mod.rs +++ b/crates/engine/src/cmd/mod.rs @@ -57,6 +57,10 @@ pub static ALL: &[CommandSpec] = &[ define::LOAD_CREATE, define::LOAD_EDIT, define::LOAD_DELETE, + define::FLUID_DOMAIN_CREATE, + define::FLUID_DOMAIN_EDIT, + define::FLUID_DOMAIN_DELETE, + define::FLUID_BOUNDARY, define::STEP_CREATE, define::STEP_EDIT, define::STEP_DELETE, @@ -75,6 +79,7 @@ pub static ALL: &[CommandSpec] = &[ solve::FREQUENCY, solve::EXPLICIT, solve::HISTORY, + solve::FLUID, solve::PROBE, results::DEFORMED, results::CONTOUR, diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index e63146d..5e9990f 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -1,11 +1,14 @@ -//! `solve.*`: run the model's steps with the implicit or the explicit solver, keep the results -//! in the session, and read them back. +//! `solve.*`: run the model's steps with the implicit, the explicit or the lattice Boltzmann +//! solver, keep the results in the session, and read them back. use modelcraft_explicit::{ExplicitResult, solve_explicit}; +use modelcraft_lbm::{FluidResult, Obstacle}; +use modelcraft_mesh::surface::{BoundarySurface, SurfaceOptions}; +use modelcraft_model::FluidDomain; use modelcraft_solver::{FrequencyOptions, Limits, Location, MassKind, StepResult as ImplicitResult, solve_frequency, solve_static}; use serde_json::{Value, json}; -use super::{labels, opt, req}; +use super::{labels, opt, points, req}; use crate::{CommandKind, CommandSpec, EngineError, Result, Session, StepResult}; pub const STATIC: CommandSpec = CommandSpec { @@ -40,10 +43,18 @@ pub const HISTORY: CommandSpec = CommandSpec { run: history, }; +pub const FLUID: CommandSpec = CommandSpec { + id: "solve.fluid", + label: "Solve Fluid Step", + params: r#"{"step": "Flow", "max_cells": 16000000, "max_steps": 10000000, "max_cell_updates": 1e12, "max_frame_bytes": 2e9} — runs a fluid step with the lattice Boltzmann solver (the only fluid step when `step` is left out): obstacles are triangulated from the open design or taken from imported surfaces; returns the time step, relaxation time, Mach number, mass and the final force on every wall; the caps may only be lowered"#, + kind: CommandKind::Model, + run: run_fluid, +}; + pub const PROBE: CommandSpec = CommandSpec { id: "solve.probe", label: "Probe Results", - params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1, "frame": 0} — the stored result at nodes: displacement, reaction, stress and von Mises (static), the shape of mode `mode` (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default); a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments)"#, + params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1, "frame": 0} — the stored result at nodes: displacement, reaction, stress and von Mises (static), the shape of mode `mode` (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default); a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity and pressure in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories"#, kind: CommandKind::Model, run: probe, }; @@ -51,6 +62,9 @@ pub const PROBE: CommandSpec = CommandSpec { /// Most nodes one probe lists. const MAX_PROBE: usize = 10_000; +/// Most force histories one fluid probe lists. +const MAX_PROBE_FORCES: usize = 100; + /// Most history samples one `solve.history` returns. const MAX_HISTORY_SAMPLES: usize = 10_000; @@ -206,6 +220,121 @@ fn history(s: &mut Session, p: &Value) -> Result { })) } +/// The surfaces of a fluid domain's obstacles: bodies of the open design (triangulated at the +/// domain's tolerance) or imported surfaces, by name. +fn obstacles(s: &Session, domain: &FluidDomain) -> Result> { + let tolerance = domain.tolerance.unwrap_or(domain.spacing / 10.0); + let mut out = Vec::with_capacity(domain.obstacles.len()); + for body in &domain.obstacles { + let (mesh, faces) = if let Some(shape) = s.design.as_deref().and_then(|d| d.body(body)) { + (std::borrow::Cow::Owned(shape.tessellate(tolerance)?), shape.faces().iter().map(|f| f.name.clone()).collect()) + } else if let Some(t) = s.surfaces.iter().find(|t| &t.name == body) { + // Imported faces have no names: the solver calls them face0, face1, …. + (std::borrow::Cow::Borrowed(&t.mesh), Vec::new()) + } else { + return Err(EngineError::Params(format!("obstacle `{body}` is neither a body of the open design nor an imported surface"))); + }; + let options = SurfaceOptions { weld_distance: mesh.weld_distance, ..SurfaceOptions::default() }; + let surface = BoundarySurface::new(&mesh.positions, &mesh.triangles, &mesh.tri_face, &options) + .map_err(|e| EngineError::Params(format!("obstacle `{body}`: its surface cannot bound a solid: {e}")))?; + out.push(Obstacle { body: body.clone(), surface, faces }); + } + Ok(out) +} + +fn run_fluid(s: &mut Session, p: &Value) -> Result { + let defaults = modelcraft_lbm::Limits::default(); + let limits = modelcraft_lbm::Limits { + max_cells: cap(p, "max_cells", defaults.max_cells)?, + max_steps: cap(p, "max_steps", defaults.max_steps)?, + max_cell_updates: cap(p, "max_cell_updates", defaults.max_cell_updates)?, + max_frame_bytes: cap(p, "max_frame_bytes", defaults.max_frame_bytes)?, + }; + let step = opt::(p, "step")?; + let (_, domain) = modelcraft_lbm::fluid_step(&s.model, step.as_deref())?; + let obstacles = obstacles(s, domain)?; + let r = modelcraft_lbm::solve_fluid(&s.model, step.as_deref(), &obstacles, &limits)?; + let summary = fluid_summary(&r); + s.store_result(StepResult::Fluid(r)); + Ok(summary) +} + +fn fluid_summary(r: &FluidResult) -> Value { + let speed = |v: &[f32; 3]| (f64::from(v[0]).powi(2) + f64::from(v[1]).powi(2) + f64::from(v[2]).powi(2)).sqrt(); + let last = r.frames.last(); + let max_speed = last.map_or(0.0, |f| f.velocity.iter().map(speed).fold(0.0, f64::max)); + // The mean over the last quarter of the history smooths the acoustic noise of the + // weakly compressible scheme. + let tail = r.history_times.len().div_ceil(4).max(1); + let forces: Vec = r + .forces + .iter() + .filter(|f| f.force.iter().any(|v| *v != [0.0; 3])) + .map(|f| { + let recent = f.force.get(f.force.len().saturating_sub(tail)..).unwrap_or_default(); + let n = recent.len().max(1) as f64; + let mean = [0, 1, 2].map(|a| recent.iter().map(|v| v[a]).sum::() / n); + json!({"name": f.name, "body": f.body, "final": f.force.last(), "mean": mean}) + }) + .collect(); + let (m0, m1) = (r.frames.first().map_or(0.0, |f| f.mass), last.map_or(0.0, |f| f.mass)); + json!({ + "step": r.step, + "cells": r.stats.cells, + "fluid_cells": r.stats.fluid_cells, + "steps": r.stats.steps, + "time_step": r.stats.time_step, + "end_time": r.stats.end_time, + "tau": r.stats.tau_plus, + "tau_minus": r.stats.tau_minus, + "magic": r.stats.magic, + "velocity_scale": r.stats.velocity, + "lattice_velocity": r.stats.lattice_velocity, + "mach": r.stats.mach, + "frames": r.frames.iter().map(|f| f.time).collect::>(), + "max_speed": max_speed, + "mass": {"first": m0, "last": m1, "relative_change": if m0 != 0.0 { (m1 - m0) / m0 } else { 0.0 }}, + "forces": forces, + "warnings": r.warnings, + }) +} + +fn probe_fluid(r: &FluidResult, p: &Value) -> Result { + if let Some(names) = opt::>(p, "forces")? { + if names.len() > MAX_PROBE_FORCES { + return Err(EngineError::Params(format!("`forces` lists at most {MAX_PROBE_FORCES} walls"))); + } + let out: Vec = names + .iter() + .map(|n| { + let f = + r.forces.iter().find(|f| &f.name == n || f.body.as_ref().is_some_and(|b| format!("{b}/{}", f.name) == *n)).ok_or_else(|| { + EngineError::Params(format!("no wall `{n}` (box faces x_min … z_max, or obstacle faces by name or body/name)")) + })?; + Ok(json!({"name": f.name, "body": f.body, "force": f.force})) + }) + .collect::>()?; + return Ok(json!({"step": r.step, "times": r.history_times, "forces": out})); + } + let pts = points(p, "points", MAX_PROBE)?; + let frame = match opt::(p, "frame")? { + Some(k) if k < r.frames.len() => k, + Some(_) => return Err(EngineError::Params(format!("`frame` must be from 0 to {}", r.frames.len().saturating_sub(1)))), + None => r.frames.len().saturating_sub(1), + }; + let time = r.frames.get(frame).map(|f| f.time); + let out: Vec = pts + .iter() + .map(|&pt| match r.probe(frame, pt) { + Some((cell, velocity, pressure)) => { + json!({"point": pt, "cell": cell, "solid": r.solid.get(cell), "velocity": velocity, "pressure": pressure}) + } + None => json!({"point": pt, "cell": null}), + }) + .collect(); + Ok(json!({"step": r.step, "frame": frame, "time": time, "points": out})) +} + fn probe(s: &mut Session, p: &Value) -> Result { let step = opt::(p, "step")?.ok_or_else(|| EngineError::Params("`step` is required".into()))?; let result = s.result(&step).ok_or_else(|| { @@ -215,9 +344,13 @@ fn probe(s: &mut Session, p: &Value) -> Result { format!("step `{step}` has no result; solve it first") }) })?; + if let StepResult::Fluid(r) = result.as_ref() { + return probe_fluid(r, p); + } let nodes = labels(p, "nodes", MAX_PROBE)?; let result = match result.as_ref() { StepResult::Explicit(r) => return probe_frame(r, &step, &nodes, p), + StepResult::Fluid(r) => return probe_fluid(r, p), StepResult::Implicit(r) => r, }; let locate = |label: u32| { @@ -409,6 +542,146 @@ pub(crate) mod tests { } } + /// Couette flow between y = 0 and y = 1 (eight cells), the top wall sliding at 0.1 along x, + /// periodic in x and z: material, domain and fluid step through commands. + pub(crate) fn couette(s: &mut Session) { + let periodic = json!({"type": "periodic"}); + let script = json!([ + {"command": "model.material.create", "params": {"name": "Oil", "density": 1.0, "fluid": {"type": "newtonian", "viscosity": 0.1}}}, + {"command": "model.fluid_domain.create", "params": {"name": "Gap", "material": "Oil", "min": [0, 0, 0], "max": [0.125, 1, 0.125], "spacing": 0.125, + "faces": {"x_min": periodic, "x_max": periodic, "z_min": periodic, "z_max": periodic}}}, + {"command": "model.fluid_domain.boundary", "params": {"domain": "Gap", "at": "y_max", "condition": {"type": "moving_wall", "velocity": [0.1, 0, 0]}}}, + {"command": "model.step.create", "params": {"name": "Shear", "procedure": {"type": "fluid", "domain": "Gap", "end_time": 50, "frames": 4}}}, + ]); + s.run_script(&script).unwrap(); + } + + #[test] + fn solve_a_fluid_step_and_probe_its_fields_and_wall_forces() { + let mut s = Session::default(); + couette(&mut s); + let journal = s.journal.len(); + let r = s.execute("solve.fluid", &json!({})).unwrap(); + assert_eq!(r["step"], "Shear"); + assert_eq!(r["frames"].as_array().unwrap().len(), 5, "{r}"); + assert!((r["tau"].as_f64().unwrap() - 1.0).abs() < 1e-3, "the automatic time step keeps τ ≤ 1: {r}"); + assert!(r["mass"]["relative_change"].as_f64().unwrap().abs() < 1e-12, "{r}"); + assert_eq!(s.journal.len(), journal, "solving does not change the model"); + // The wall shear μ U / H over the wall's area 0.125². + let bottom = r["forces"].as_array().unwrap().iter().find(|f| f["name"] == "y_min").unwrap(); + let expected = 0.1 * 0.1 / 1.0 * 0.125 * 0.125; + assert!((bottom["mean"][0].as_f64().unwrap() / expected - 1.0).abs() < 1e-3, "{bottom}"); + + let probe = s.execute("solve.probe", &json!({"step": "Shear", "points": [[0.06, 0.5625, 0.06], [5, 5, 5]]})).unwrap(); + let u = probe["points"][0]["velocity"][0].as_f64().unwrap(); + assert!((u - 0.05625).abs() < 1e-5, "{probe}"); + assert!(probe["points"][1]["cell"].is_null()); + let first = s.execute("solve.probe", &json!({"step": "Shear", "points": [[0.06, 0.5625, 0.06]], "frame": 0})).unwrap(); + assert_eq!(first["points"][0]["velocity"][0], 0.0, "the initial frame is at rest"); + let h = s.execute("solve.probe", &json!({"step": "Shear", "forces": ["y_min", "y_max"]})).unwrap(); + assert_eq!(h["forces"].as_array().unwrap().len(), 2); + assert_eq!(h["times"].as_array().unwrap().len(), h["forces"][0]["force"].as_array().unwrap().len()); + + // An edit makes the result stale; the step solves again. + s.execute( + "model.fluid_domain.boundary", + &json!({"domain": "Gap", "at": "y_max", "condition": {"type": "moving_wall", "velocity": [0.2, 0, 0]}}), + ) + .unwrap(); + assert!(s.execute("solve.probe", &json!({"step": "Shear", "points": [[0, 0, 0]]})).is_err()); + // The domain cannot go while a step uses it; the material cannot lose its fluid behaviour. + assert!(s.execute("model.fluid_domain.delete", &json!({"name": "Gap"})).is_err()); + assert!(s.execute("model.material.edit", &json!({"name": "Oil", "density": 1.0})).is_err()); + // A structural step is not a fluid one, and a fluid step is not written to decks. + assert!(s.execute("solve.static", &json!({"step": "Shear"})).is_err()); + for format in ["abaqus", "nastran", "lsdyna"] { + let e = s.execute("deck.write", &json!({"format": format, "path": "/no/such/dir/x.deck"})).unwrap_err().to_string(); + assert!(e.contains("fluid step"), "{format}: {e}"); + } + } + + #[test] + fn a_fluid_domain_takes_its_obstacles_from_the_open_design() { + use modelcraft_geometry::sample::{Plate, plate_design}; + let dir = std::env::temp_dir().join(format!("modelcraft-fluid-{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let path = dir.join("plate.solidcraft"); + std::fs::write(&path, plate_design(&Plate::default()).unwrap()).unwrap(); + let mut s = Session::default(); + s.execute("geometry.open", &json!({"path": path.to_str().unwrap()})).unwrap(); + let _ = std::fs::remove_dir_all(&dir); + let bodies = s.execute("geometry.bodies", &json!({})).unwrap(); + let [lo, hi] = [0, 1].map(|k| bodies[0]["bbox"][k].as_array().unwrap().iter().map(|v| v.as_f64().unwrap()).collect::>()); + // A box 20 mm larger than the plate on every side, 5 mm cells: a viscous liquid (ν = 100 mm²/s) flowing along x. + let min: Vec = lo.iter().map(|v| (v / 5.0).floor() * 5.0 - 20.0).collect(); + let max: Vec = hi.iter().map(|v| (v / 5.0).ceil() * 5.0 + 20.0).collect(); + let script = json!([ + {"command": "model.material.create", "params": {"name": "Water", "density": 1e-9, "fluid": {"type": "newtonian", "viscosity": 1e-7}}}, + {"command": "model.fluid_domain.create", "params": {"name": "Tank", "material": "Water", "min": min, "max": max, "spacing": 5, "obstacles": ["Body1"], + "faces": {"x_min": {"type": "velocity", "velocity": [10, 0, 0]}, "x_max": {"type": "pressure", "pressure": 0}}}}, + {"command": "model.step.create", "params": {"name": "Flow", "procedure": {"type": "fluid", "domain": "Tank", "steps": 60, "frames": 2}}}, + ]); + s.run_script(&script).unwrap(); + let r = s.execute("solve.fluid", &json!({"step": "Flow"})).unwrap(); + assert!(r["fluid_cells"].as_u64().unwrap() < r["cells"].as_u64().unwrap(), "{r}"); + let on_plate: Vec<&Value> = r["forces"].as_array().unwrap().iter().filter(|f| f["body"] == "Body1").collect(); + assert!(!on_plate.is_empty(), "{r}"); + let faces = bodies[0]["faces"].as_array().unwrap(); + for f in &on_plate { + assert!(faces.iter().any(|g| g["name"] == f["name"]), "forces are reported per named face: {f}"); + } + let drag: f64 = on_plate.iter().map(|f| f["final"][0].as_f64().unwrap()).sum(); + assert!(drag > 0.0, "the flow pushes the plate downstream: {r}"); + // An obstacle that is not in the design. + s.execute( + "model.fluid_domain.edit", + &json!({"name": "Tank", "material": "Water", "min": min, "max": max, "spacing": 5, "obstacles": ["Nope"]}), + ) + .unwrap(); + assert!(s.execute("solve.fluid", &json!({})).is_err()); + } + + #[test] + fn solve_fluid_hostile_params_never_panic() { + let mut s = Session::default(); + couette(&mut s); + s.execute("solve.fluid", &json!({"step": "Shear"})).unwrap(); + let hostile = [ + json!(null), + json!({"step": 7}), + json!({"step": "Nope"}), + json!({"step": "Shear", "max_cells": 0}), + json!({"step": "Shear", "max_cells": 1e300}), + json!({"step": "Shear", "max_cells": 3}), + json!({"step": "Shear", "max_steps": 3}), + json!({"step": "Shear", "max_steps": -1}), + json!({"step": "Shear", "max_cell_updates": 10}), + json!({"step": "Shear", "max_cell_updates": 0.5}), + json!({"step": "Shear", "max_frame_bytes": 16}), + json!({"step": "Shear", "max_frame_bytes": "lots"}), + json!({"step": "Shear", "points": [[1e308, -1e308, 0]]}), + json!({"step": "Shear", "points": [[f64::MAX, 0, 0]], "frame": 1e300}), + json!({"step": "Shear", "points": [[0, 0, 0]], "frame": 99}), + json!({"step": "Shear", "points": null}), + json!({"step": "Shear", "points": vec![[0, 0, 0]; 20_000]}), + json!({"step": "Shear", "forces": ["nope"]}), + json!({"step": "Shear", "forces": vec!["y_min"; 1000]}), + json!({"step": "Shear", "forces": "y_min"}), + json!({"step": "Shear", "nodes": [1]}), + ]; + for id in ["solve.fluid", "solve.probe", "solve.static", "solve.frequency"] { + for p in &hostile { + let mut t = s.clone(); + let before = t.model.clone(); + let results = t.solved_steps().count(); + if t.execute(id, p).is_err() { + assert_eq!(t.model, before); + assert_eq!(t.solved_steps().count(), results, "{id} {p}"); + } + } + } + } + #[test] fn solve_and_probe_a_coupled_reference_point() { // The defined bar: a force on a reference point kinematically coupled to the free end. diff --git a/crates/engine/src/lib.rs b/crates/engine/src/lib.rs index ebe507e..f7a829e 100644 --- a/crates/engine/src/lib.rs +++ b/crates/engine/src/lib.rs @@ -50,6 +50,8 @@ pub enum EngineError { Solver(#[from] modelcraft_solver::Error), #[error(transparent)] Explicit(#[from] modelcraft_explicit::Error), + #[error(transparent)] + Fluid(#[from] modelcraft_lbm::Error), /// A command panicked; the model is unchanged. #[error("internal error in {command}: {message}")] Internal { command: String, message: String }, @@ -57,11 +59,12 @@ pub enum EngineError { pub type Result = std::result::Result; -/// The result of a solved step, from the implicit or the explicit solver. +/// The result of a solved step, from the implicit, the explicit or the lattice Boltzmann solver. #[derive(Debug, Clone, PartialEq)] pub enum StepResult { Implicit(modelcraft_solver::StepResult), Explicit(ExplicitResult), + Fluid(modelcraft_lbm::FluidResult), } impl StepResult { @@ -69,6 +72,7 @@ impl StepResult { match self { StepResult::Implicit(r) => r.step(), StepResult::Explicit(r) => &r.step, + StepResult::Fluid(r) => &r.step, } } @@ -77,6 +81,8 @@ impl StepResult { match self { StepResult::Implicit(r) => r.nodes(), StepResult::Explicit(r) => &r.nodes, + // A fluid result is per cell, not per node. + StepResult::Fluid(_) => &[], } } } @@ -357,6 +363,19 @@ mod tests { json!({"step": "Bend", "loads": ["nope"], "active": "yes"}), json!({"step": "Bend", "boundary_conditions": [], "loads": []}), json!({"system": "furlongs"}), + // Fluid domains, boundaries and steps. + json!({"name": "A", "material": "Steel", "min": [0, 0, 0], "max": [1e308, 1, 1], "spacing": 1e-300}), + json!({"name": "A", "material": "Steel", "min": [0, 0, 0], "max": [1, 1, 1], "spacing": 0.5, "faces": {"x_min": {"type": "periodic"}}}), + json!({"name": "A", "material": "nope", "min": [0, 0, 0], "max": [1, 1, 1], "spacing": 0.5, "obstacles": ["", ""]}), + json!({"name": "A", "density": 1, "fluid": {"type": "newtonian", "viscosity": -1}}), + json!({"name": "A", "fluid": {"type": "newtonian", "viscosity": 1}}), + json!({"domain": "A", "at": "x_mid", "condition": {"type": "wall"}}), + json!({"domain": "A", "at": {"body": "B", "face": "f"}, "condition": {"type": "pressure", "pressure": 1e308}}), + json!({"domain": null, "at": null, "condition": null}), + json!({"name": "A", "procedure": {"type": "fluid", "domain": "A", "steps": 1e300}}), + json!({"name": "A", "procedure": {"type": "fluid", "domain": "A", "end_time": -1, "steps": 3}}), + json!({"name": "Bend", "procedure": {"type": "fluid", "domain": "nope", "end_time": 1}}), + json!({"step": "Bend", "max_cells": 1e300, "max_steps": -3}), json!({"geometry": {"body": "Bar"}, "elements": [1, 99], "faces": [{"element": 1, "face": 200}]}), json!({"geometry": {"body": "Bar", "vertex": "v"}, "nodes": [1], "edges": [{"element": 1, "edge": 0}]}), ]; diff --git a/crates/engine/src/results.rs b/crates/engine/src/results.rs index 19fd378..1212844 100644 --- a/crates/engine/src/results.rs +++ b/crates/engine/src/results.rs @@ -20,6 +20,8 @@ pub enum Kind<'a> { Static(&'a StaticResult), Frequency(&'a FrequencyResult), Explicit(&'a ExplicitResult), + /// Per cell, not per node: no nodal fields to display yet. + Fluid(&'a modelcraft_lbm::FluidResult), } impl StepResult { @@ -28,6 +30,7 @@ impl StepResult { StepResult::Implicit(modelcraft_solver::StepResult::Static(r)) => Kind::Static(r), StepResult::Implicit(modelcraft_solver::StepResult::Frequency(r)) => Kind::Frequency(r), StepResult::Explicit(r) => Kind::Explicit(r), + StepResult::Fluid(r) => Kind::Fluid(r), } } } @@ -146,8 +149,10 @@ impl Field { Kind::Static(_) => (true, false), Kind::Frequency(_) => (false, false), Kind::Explicit(e) => e.frames.last().map_or((false, false), |f| (!f.stress.is_empty(), !f.velocity.is_empty())), + // Cell fields have their own display (M12 streamlines); no nodal field. + Kind::Fluid(_) => return None.into_iter().flatten(), }; - Field::ALL.into_iter().filter(move |f| if *f == Field::Velocity { velocity } else { stress || !f.is_stress() }) + Some(Field::ALL.into_iter().filter(move |f| if *f == Field::Velocity { velocity } else { stress || !f.is_stress() })).into_iter().flatten() } } @@ -222,6 +227,7 @@ pub fn displacements(r: &StepResult, pick: Pick) -> Result<&[[f64; 3]]> { } Ok(&f.displacement) } + Kind::Fluid(f) => Err(EngineError::Params(format!("step `{}` is a fluid step: its fields are per cell (solve.probe)", f.step))), } } @@ -233,7 +239,7 @@ pub fn field_values(model: &Model, r: &StepResult, field: Field, pick: Pick) -> let stresses = || -> Result> { match r.kind() { Kind::Static(s) => Ok(std::borrow::Cow::Borrowed(s.stress.as_slice())), - Kind::Frequency(_) => Err(unavailable()), + Kind::Frequency(_) | Kind::Fluid(_) => Err(unavailable()), Kind::Explicit(e) => { let f = frame(e, pick)?; if f.stress.is_empty() { @@ -255,7 +261,7 @@ pub fn field_values(model: &Model, r: &StepResult, field: Field, pick: Pick) -> let f = frame(e, pick)?; if f.velocity.is_empty() { Err(unavailable()) } else { Ok(f.velocity.iter().map(|v| magnitude(*v)).collect()) } } - Kind::Static(_) | Kind::Frequency(_) => Err(unavailable()), + Kind::Static(_) | Kind::Frequency(_) | Kind::Fluid(_) => Err(unavailable()), }, Field::VonMises => each(&modelcraft_fem::von_mises), Field::MaxPrincipal => principal(0), @@ -338,6 +344,9 @@ impl Session { out["time"] = json!(f.time); } } + Kind::Fluid(f) => { + out["frames"] = f.frames.iter().map(|f| f.time).collect(); + } Kind::Static(_) => {} } out diff --git a/crates/explicit/tests/verification.rs b/crates/explicit/tests/verification.rs index 3746a79..99a2e3b 100644 --- a/crates/explicit/tests/verification.rs +++ b/crates/explicit/tests/verification.rs @@ -77,7 +77,8 @@ impl Bar { } } m.add_set("end_face", Set { faces, ..Set::default() }).unwrap(); - m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: Some(DENSITY) }).unwrap(); + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: Some(DENSITY), fluid: None }) + .unwrap(); let integration = if kind == ElementKind::Hex8 { Integration::Reduced } else { Integration::Full }; let formulation = Formulation { integration, ..Formulation::default() }; m.add_section("bar", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap(); @@ -333,10 +334,13 @@ fn caps_and_unsupported_models_are_errors() { assert!(solve_explicit(&full.model, None, &Limits::default()).unwrap_err().to_string().contains("not supported yet")); let mut ortho = Bar::new(ElementKind::Tet4, 4, ExplicitDynamic::new(end)); let constants = Elastic::EngineeringConstants { young: [YOUNG; 3], poisson: [0.0; 3], shear: [YOUNG / 2.0; 3] }; - ortho.model.replace_material("Mat", Material { elastic: Some(constants), density: Some(DENSITY) }).unwrap(); + ortho.model.replace_material("Mat", Material { elastic: Some(constants), density: Some(DENSITY), fluid: None }).unwrap(); assert!(solve_explicit(&ortho.model, None, &Limits::default()).unwrap_err().to_string().contains("orthotropic")); let mut light = Bar::new(ElementKind::Tet4, 4, ExplicitDynamic::new(end)); - light.model.replace_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: None }).unwrap(); + light + .model + .replace_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: YOUNG, poisson: 0.0 }), density: None, fluid: None }) + .unwrap(); assert!(solve_explicit(&light.model, None, &Limits::default()).is_err()); // A node behind its wall at the start. let mut behind = ExplicitDynamic::new(end); diff --git a/crates/lbm/Cargo.toml b/crates/lbm/Cargo.toml new file mode 100644 index 0000000..eb8bdbd --- /dev/null +++ b/crates/lbm/Cargo.toml @@ -0,0 +1,14 @@ +[package] +name = "modelcraft-lbm" +description = "ModelCraft's lattice Boltzmann flow solver: D3Q19 two-relaxation-time collision, bounce-back walls, Zou–He inlets and outlets, exact voxelisation of face-tagged surfaces, momentum-exchange wall forces" +version.workspace = true +edition.workspace = true +license.workspace = true + +[lints] +workspace = true + +[dependencies] +modelcraft-model = { workspace = true } +modelcraft-mesh = { workspace = true } +thiserror = { workspace = true } diff --git a/crates/lbm/src/flow.rs b/crates/lbm/src/flow.rs new file mode 100644 index 0000000..f4540d1 --- /dev/null +++ b/crates/lbm/src/flow.rs @@ -0,0 +1,598 @@ +//! The lattice Boltzmann state and one time step: collide, stream, open boundaries. +//! +//! # Scheme +//! - **Lattice and equilibrium:** D3Q19 with the second-order equilibrium +//! (Qian, d'Humières & Lallemand 1992; Krüger et al. 2017, ch. 3). +//! - **Collision:** two relaxation times (Ginzburg, Verhaeghe & d'Humières 2008): the symmetric +//! part of each opposite pair relaxes with `1/τ⁺` (which sets the viscosity), the antisymmetric +//! part with `1/τ⁻`, related by the parameter `Λ = (τ⁺ − ½)(τ⁻ − ½)`. With `Λ = 3/16` +//! halfway bounce-back puts the wall exactly halfway between nodes and steady solutions do not +//! depend on the viscosity. +//! - **Body force:** Guo, Zheng & Shi 2002 (PRE 65:046308): the velocity is +//! `(Σ f c + F/2) / ρ` and the forcing term is split into its symmetric part, relaxed with +//! `1 − 1/(2τ⁺)`, and its antisymmetric part, with `1 − 1/(2τ⁻)` (the TRT form of Guo's +//! term, Krüger et al. 2017, ch. 10). +//! - **Walls:** halfway bounce-back at every link from a fluid cell to a solid cell or across a +//! wall face of the box; moving walls add Ladd's term `6 w_q ρ₀ c_q·u_w` (Ladd 1994, JFM +//! 271:285). Slip faces reflect populations specularly. +//! - **Inlets and outlets:** Zou & He 1997 (Phys. Fluids 9:1591) on the first layer of cells: the +//! unknown incoming populations are the opposite ones plus the equilibrium difference +//! (bounce-back of the non-equilibrium part), the density (velocity inlet) or the normal velocity +//! (pressure outlet) follow from mass and normal momentum, and the transverse momentum is +//! corrected on the diagonal unknowns as M. Hecht and J. Harting do for D3Q19 (*Implementation +//! of on-site velocity boundary conditions for D3Q19 lattice Boltzmann simulations*, J. Stat. +//! Mech. 2010, P01018). The prescribed velocity is the one with the half force of Guo's scheme. +//! - **Forces:** momentum exchange on every wall link (Ladd 1994; R. Mei, D. Yu, W. Shyy, L.-S. +//! Luo, *Force evaluation in the lattice Boltzmann method involving curved geometry*, PRE +//! 65:041203, 2002): a population `f*` leaving a fluid cell toward the wall and the one `f'` +//! coming back transfer `c_out f* − c_in f'` to the wall. +//! +//! # Layout and determinism +//! Two arrays of populations, structure of arrays (`f[q n + cell]`), A–B pattern. Streaming +//! pulls: every population of a cell comes from one source (a neighbour, a periodic image, the +//! cell itself through a wall, or an open boundary). Cells whose 18 neighbours are all fluid +//! cells inside the box stream by one shifted copy per direction; the others follow per-link +//! rules. Every loop runs in a fixed order with no parallel reductions, so the result is the +//! same bit for bit on every platform. + +use modelcraft_model::{BoxFace, FaceCondition, Profile}; + +use crate::lattice::{C, Q, W, c, equilibrium, moving_wall, opposite}; +use crate::setup::{MAX_LATTICE_SPEED, Obstacle, Setup}; +use crate::voxel::{Grid, TriangleTree, voxelise}; +use crate::{Error, Result}; + +/// Where a population of a cell next to a wall or box face comes from. +#[derive(Debug, Clone, Copy, PartialEq)] +enum Rule { + /// The same direction from another cell (a neighbour or a periodic image). + Pull { from: u32 }, + /// Bounce-back at a wall of force group `group`; `du` is the moving-wall term. + Bounce { group: u32, du: f64 }, + /// Specular reflection at a slip face: direction `q` from cell `from`. + Reflect { from: u32, q: u8, group: u32 }, + /// From outside an inlet or outlet: set by the open boundary after streaming. + Open, +} + +/// A population, direction `q` of `cell`, and where it comes from. +#[derive(Debug, Clone, Copy)] +struct Link { + cell: u32, + q: u8, + rule: Rule, +} + +/// An inlet or outlet cell. +#[derive(Debug, Clone, Copy)] +struct OpenCell { + cell: u32, + axis: usize, + /// +1 at a min face (the inward normal points up the axis), −1 at a max face. + inward: f64, + target: Target, +} + +#[derive(Debug, Clone, Copy)] +enum Target { + /// The velocity (lattice units). + Velocity([f64; 3]), + /// The density. + Density(f64), +} + +/// A force group: a box face or a face of an obstacle. +#[derive(Debug, Clone, PartialEq)] +pub struct Group { + /// `x_min` … for box faces; the face name for obstacle faces. + pub name: String, + /// The obstacle body, for obstacle faces. + pub body: Option, +} + +/// The flow state. +#[derive(Debug, Clone)] +pub struct Flow { + grid: Grid, + n: usize, + /// The distance between two directions' arrays in `f` and `g` (at least `n`). + stride: usize, + fluid: Vec, + f: Vec, + g: Vec, + /// The links that do not stream by the shifted copies, in cell order. + links: Vec, + open: Vec, + omega_plus: f64, + omega_minus: f64, + acceleration: [f64; 3], + /// The momentum transferred to each group in the last step (lattice units). + forces: Vec<[f64; 3]>, + groups: Vec, + /// Per direction, the index offset of the neighbour a population streams from. + offset: [isize; Q], + fluid_cells: usize, + steps_done: u64, + /// Obstacles that cover no cell centre (smaller than the spacing, or outside the box). + empty_obstacles: Vec, +} + +const NONE: u32 = u32::MAX; + +/// A vector allocated without aborting when memory runs out. +fn filled(len: usize, value: T, what: &str) -> Result> { + let mut v = Vec::new(); + v.try_reserve_exact(len).map_err(|_| Error::Limit(format!("not enough memory for the {what} ({len} entries)")))?; + v.resize(len, value); + Ok(v) +} + +impl Flow { + /// Build the grid (voxelising the obstacles), the streaming rules and the initial state (at + /// rest or the domain's initial velocity, density 1). + pub fn new(setup: &Setup, obstacles: &[Obstacle]) -> Result { + let grid = setup.grid; + let [nx, ny, nz] = grid.cells; + let n = grid.len(); + if n == 0 || n >= NONE as usize { + return Err(Error::Limit(format!("{n} cells cannot be addressed"))); + } + // Each direction's array starts an odd number of cache lines (8 doubles) after the + // previous one, so the 19 arrays never fall into the same cache sets (a power-of-two grid + // would otherwise put them all in one). + let lines = n.div_ceil(8); + let stride = 8 * (lines + 1 - lines % 2); + let q_n = stride.checked_mul(Q).ok_or_else(|| Error::Limit("the grid is too large".into()))?; + + // Solid cells, obstacle by obstacle (the first obstacle that covers a cell owns it). + let mut fluid = filled(n, true, "cell flags")?; + let mut owner = filled(n, NONE, "cell owners")?; + let mut trees = Vec::with_capacity(setup.obstacles.len()); + let mut empty_obstacles = Vec::new(); + // Force groups: the six box faces, then every face of every obstacle. + let mut groups: Vec = BoxFace::ALL.iter().map(|f| Group { name: f.name().to_string(), body: None }).collect(); + let mut base = Vec::with_capacity(setup.obstacles.len()); + for (k, faces) in setup.obstacles.iter().enumerate() { + let o = obstacles.get(faces.index).ok_or_else(|| Error::NotFound(format!("obstacle {}", faces.index)))?; + let mut inside = filled(n, false, "obstacle cells")?; + if voxelise(&o.surface, &grid, &mut inside)? == 0 { + empty_obstacles.push(o.body.clone()); + } + for (c, &solid) in inside.iter().enumerate() { + if solid && owner[c] == NONE { + owner[c] = k as u32; + fluid[c] = false; + } + } + trees.push((TriangleTree::new(&o.surface), o.surface.tri_face())); + base.push(groups.len() as u32); + groups.extend(faces.names.iter().map(|name| Group { name: name.clone(), body: Some(o.body.clone()) })); + } + + let domain = &setup.domain; + let periodic = [domain.faces.periodic(0), domain.faces.periodic(1), domain.faces.periodic(2)]; + let wall_velocity = |face: BoxFace| match *domain.faces.get(face) { + FaceCondition::MovingWall { velocity } => Some(setup.to_lattice(velocity)), + FaceCondition::Wall {} => Some([0.0; 3]), + _ => None, + }; + let is_slip = |face: BoxFace| matches!(domain.faces.get(face), FaceCondition::Slip {}); + let dims = [nx as i64, ny as i64, nz as i64]; + let index = |p: [i64; 3]| grid.index([p[0] as usize, p[1] as usize, p[2] as usize]) as u32; + let mut offset = [0isize; Q]; + for (q, o) in offset.iter_mut().enumerate() { + let cq = C[q]; + *o = cq[0] as isize + nx as isize * (cq[1] as isize + ny as isize * cq[2] as isize); + } + + // The group and moving-wall term of a link into fluid cell `cell` (direction q) from the + // obstacle cell `solid`: the face of the triangle nearest the link's midpoint. + let obstacle_link = |cell: usize, solid: u32, q: usize| -> Rule { + let k = owner.get(solid as usize).copied().unwrap_or(NONE) as usize; + let (Some((tree, tags)), Some(faces), Some(&first)) = (trees.get(k), setup.obstacles.get(k), base.get(k)) else { + return Rule::Bounce { group: 0, du: 0.0 }; + }; + let (centre, cq) = (grid.centre(cell), c(q)); + let mid = [0, 1, 2].map(|a| centre[a] - 0.5 * cq[a] * grid.spacing); + let tag = tree.nearest(mid).and_then(|t| tags.get(t as usize).copied()).unwrap_or(0); + let u = faces.velocities.get(tag as usize).copied().unwrap_or([0.0; 3]); + Rule::Bounce { group: first + tag, du: moving_wall(q, u) } + }; + + // The links that do not stream from the neighbour the shifted copies read. + let mut links: Vec = Vec::new(); + let mut fluid_cells = 0; + for cell in 0..n { + if !fluid[cell] { + continue; + } + fluid_cells += 1; + let ijk = grid.coords(cell); + let p = [ijk[0] as i64, ijk[1] as i64, ijk[2] as i64]; + for q in 1..Q { + let cq = C[q]; + let mut s = [p[0] - i64::from(cq[0]), p[1] - i64::from(cq[1]), p[2] - i64::from(cq[2])]; + let mut wrapped = false; + let mut crossed: [Option; 3] = [None; 3]; + for a in 0..3 { + if s[a] < 0 || s[a] >= dims[a] { + if periodic[a] { + s[a] = s[a].rem_euclid(dims[a]); + wrapped = true; + } else { + crossed[a] = Some(BoxFace::of(a, s[a] >= dims[a])); + } + } + } + let crossed_faces = || crossed.iter().flatten().copied(); + let rule = match crossed_faces().next() { + None => { + let from = index(s); + match (fluid.get(from as usize).copied().unwrap_or(false), wrapped) { + (true, false) => continue, + (true, true) => Rule::Pull { from }, + (false, _) => obstacle_link(cell, from, q), + } + } + Some(_) => { + if let Some(face) = crossed_faces().find(|f| wall_velocity(*f).is_some()) { + Rule::Bounce { group: face.index() as u32, du: moving_wall(q, wall_velocity(face).unwrap_or([0.0; 3])) } + } else if let Some(face) = crossed_faces().find(|f| is_slip(*f)) { + // The mirror source: the same cell layer, the tangential components shifted. + let a = face.axis(); + let mut m = s; + m[a] = p[a]; + let mut mirrored = C[q]; + mirrored[a] = -mirrored[a]; + let q2 = C.iter().position(|v| *v == mirrored).unwrap_or(opposite(q)); + let single = crossed_faces().count() == 1; + let inside = (0..3).all(|b| m[b] >= 0 && m[b] < dims[b]); + match inside.then(|| index(m)) { + Some(from) if single && fluid.get(from as usize).copied().unwrap_or(false) => { + Rule::Reflect { from, q: q2 as u8, group: face.index() as u32 } + } + _ => Rule::Bounce { group: face.index() as u32, du: 0.0 }, + } + } else { + Rule::Open + } + } + }; + links.push(Link { cell: cell as u32, q: q as u8, rule }); + } + } + + // Inlets and outlets: the fluid cells of each open face's layer. + let mut open = Vec::new(); + for face in BoxFace::ALL { + let a = face.axis(); + let layer = if face.is_max() { grid.cells[a] - 1 } else { 0 }; + let (t1, t2) = ((a + 1) % 3, (a + 2) % 3); + for j in 0..grid.cells[t2] { + for i in 0..grid.cells[t1] { + let mut ijk = [0usize; 3]; + ijk[a] = layer; + ijk[t1] = i; + ijk[t2] = j; + let cell = grid.index(ijk); + if !fluid[cell] { + continue; + } + let target = match *domain.faces.get(face) { + FaceCondition::Velocity { velocity, profile } => { + let mut scale = 1.0; + if profile == Profile::Parabolic { + for t in [t1, t2].into_iter().filter(|t| !periodic[*t]) { + let s = (ijk[t] as f64 + 0.5) / grid.cells[t] as f64; + scale *= 4.0 * s * (1.0 - s); + } + } + Target::Velocity(setup.to_lattice(velocity).map(|v| v * scale)) + } + FaceCondition::Pressure { pressure } => Target::Density(setup.lattice_density(pressure)), + _ => break, + }; + open.push(OpenCell { cell: cell as u32, axis: a, inward: if face.is_max() { -1.0 } else { 1.0 }, target }); + } + } + } + + // The initial state: equilibrium at density 1 and the initial velocity. + let mut f = filled(q_n, 0.0, "populations")?; + let g = filled(q_n, 0.0, "populations")?; + for cell in 0..n { + for q in 0..Q { + f[q * stride + cell] = if fluid[cell] { equilibrium(q, 1.0, setup.initial_velocity) } else { 0.0 }; + } + } + Ok(Flow { + grid, + n, + stride, + fluid, + f, + g, + links, + open, + omega_plus: 1.0 / setup.tau_plus, + omega_minus: 1.0 / setup.tau_minus, + acceleration: setup.acceleration, + forces: vec![[0.0; 3]; groups.len()], + groups, + offset, + fluid_cells, + steps_done: 0, + empty_obstacles, + }) + } + + pub fn grid(&self) -> &Grid { + &self.grid + } + + pub fn fluid_cells(&self) -> usize { + self.fluid_cells + } + + pub fn is_fluid(&self, cell: usize) -> bool { + self.fluid.get(cell).copied().unwrap_or(false) + } + + pub fn groups(&self) -> &[Group] { + &self.groups + } + + /// The momentum each group received in the last step (lattice units: force per step). + pub fn forces(&self) -> &[[f64; 3]] { + &self.forces + } + + /// Obstacles that cover no cell centre. + pub fn empty_obstacles(&self) -> &[String] { + &self.empty_obstacles + } + + pub fn steps_done(&self) -> u64 { + self.steps_done + } + + /// Set every fluid cell to the equilibrium of `state(cell) = (density, velocity)` in lattice + /// units (velocity as defined with the half force). + pub fn set_equilibrium(&mut self, state: impl Fn(usize) -> (f64, [f64; 3])) { + let n = self.n; + let s = self.stride; + for cell in 0..n { + if !self.fluid[cell] { + continue; + } + let (rho, u) = state(cell); + let bare = [0, 1, 2].map(|a| u[a] - 0.5 * self.acceleration[a]); + for q in 0..Q { + self.f[q * s + cell] = equilibrium(q, rho, bare); + } + } + } + + /// The density and velocity of a cell (lattice units), or `None` for a solid cell. + pub fn moments(&self, cell: usize) -> Option<(f64, [f64; 3])> { + if !self.is_fluid(cell) { + return None; + } + let s = self.stride; + let (mut rho, mut j) = (0.0, [0.0; 3]); + for q in 0..Q { + let fq = self.f.get(q * s + cell).copied().unwrap_or(0.0); + let cq = c(q); + rho += fq; + j[0] += fq * cq[0]; + j[1] += fq * cq[1]; + j[2] += fq * cq[2]; + } + let u = [0, 1, 2].map(|a| (j[a] + 0.5 * rho * self.acceleration[a]) / rho); + Some((rho, u)) + } + + /// The total mass of the fluid cells (lattice units: the sum of every population). + pub fn mass(&self) -> f64 { + let n = self.n; + let s = self.stride; + let mut m = 0.0; + for cell in (0..n).filter(|c| self.fluid[*c]) { + for q in 0..Q { + m += self.f[q * s + cell]; + } + } + m + } + + /// One time step. Fails with [`Error::Unstable`] when a density stops being a positive + /// number or a lattice speed exceeds [`MAX_LATTICE_SPEED`]. + pub fn step(&mut self) -> Result<()> { + let max_u2 = self.collide(); + if !(max_u2 <= MAX_LATTICE_SPEED * MAX_LATTICE_SPEED) { + return Err(Error::Unstable(format!( + "the flow went unstable at time step {}: the lattice speed reached {:.3} (limit {MAX_LATTICE_SPEED}) or a density stopped being a positive number; lower the time step, refine the spacing or check the boundary conditions", + self.steps_done + 1, + max_u2.sqrt() + ))); + } + self.stream(); + self.open_boundaries(); + std::mem::swap(&mut self.f, &mut self.g); + self.steps_done += 1; + Ok(()) + } + + /// Two-relaxation-time collision with Guo forcing, in place. Returns the largest squared + /// lattice speed, or NaN when a density is not positive. + fn collide(&mut self) -> f64 { + let n = self.n; + let s = self.stride; + let (wp, wm) = (self.omega_plus, self.omega_minus); + let (sp, sm) = (1.0 - 0.5 * wp, 1.0 - 0.5 * wm); + let a = self.acceleration; + let forced = a != [0.0; 3]; + let mut max_u2: f64 = 0.0; + let mut fq = [0.0; Q]; + for cell in 0..n { + if !self.fluid[cell] { + continue; + } + let (mut rho, mut j) = (0.0, [0.0; 3]); + for q in 0..Q { + let v = self.f[q * s + cell]; + fq[q] = v; + let cq = C[q]; + rho += v; + j[0] += v * f64::from(cq[0]); + j[1] += v * f64::from(cq[1]); + j[2] += v * f64::from(cq[2]); + } + if !(rho > 0.0 && rho.is_finite()) { + return f64::NAN; + } + let force = [rho * a[0], rho * a[1], rho * a[2]]; + let u = [(j[0] + 0.5 * force[0]) / rho, (j[1] + 0.5 * force[1]) / rho, (j[2] + 0.5 * force[2]) / rho]; + let uu = u[0] * u[0] + u[1] * u[1] + u[2] * u[2]; + max_u2 = max_u2.max(uu); + let uf = u[0] * force[0] + u[1] * force[1] + u[2] * force[2]; + // Rest population: symmetric only. + let eq0 = W[0] * rho * (1.0 - 1.5 * uu); + let mut out0 = fq[0] - wp * (fq[0] - eq0); + if forced { + out0 += sp * W[0] * (-3.0 * uf); + } + self.f[cell] = out0; + for q in (1..Q).step_by(2) { + let o = q + 1; + let cq = C[q]; + let cv = [f64::from(cq[0]), f64::from(cq[1]), f64::from(cq[2])]; + let cu = cv[0] * u[0] + cv[1] * u[1] + cv[2] * u[2]; + let w = W[q]; + let eq_plus = w * rho * (1.0 + 4.5 * cu * cu - 1.5 * uu); + let eq_minus = w * rho * 3.0 * cu; + let plus = 0.5 * (fq[q] + fq[o]); + let minus = 0.5 * (fq[q] - fq[o]); + let relax_plus = wp * (plus - eq_plus); + let relax_minus = wm * (minus - eq_minus); + let (mut s_plus, mut s_minus) = (0.0, 0.0); + if forced { + let cf = cv[0] * force[0] + cv[1] * force[1] + cv[2] * force[2]; + s_plus = sp * w * (9.0 * cu * cf - 3.0 * uf); + s_minus = sm * w * 3.0 * cf; + } + self.f[q * s + cell] = fq[q] - relax_plus - relax_minus + s_plus + s_minus; + self.f[o * s + cell] = fq[o] - relax_plus + relax_minus + s_plus - s_minus; + } + } + max_u2 + } + + /// Pull streaming into `g`, with the wall rules and the momentum exchange. + fn stream(&mut self) { + let n = self.n; + let s = self.stride; + // Every direction as one shifted copy: right for cells whose neighbours are all fluid + // cells in the box; the others are overwritten by their rules below. + for q in 0..Q { + let off = self.offset[q]; + let (src, dst) = (&self.f[q * s..q * s + n], &mut self.g[q * s..q * s + n]); + if off >= 0 { + let o = off as usize; + if o < n { + dst[o..].copy_from_slice(&src[..n - o]); + } + } else { + let o = off.unsigned_abs(); + if o < n { + dst[..n - o].copy_from_slice(&src[o..]); + } + } + } + self.forces.fill([0.0; 3]); + for link in &self.links { + let (cell, q) = (link.cell as usize, link.q as usize); + let cq = c(q); + match link.rule { + Rule::Pull { from } => self.g[q * s + cell] = self.f[q * s + from as usize], + Rule::Bounce { group, du } => { + let out = self.f[opposite(q) * s + cell]; + let back = out + du; + self.g[q * s + cell] = back; + if let Some(force) = self.forces.get_mut(group as usize) { + // c_out f* − c_in f' with c_out = −c_q. + for a in 0..3 { + force[a] -= cq[a] * (out + back); + } + } + } + Rule::Reflect { from, q: q2, group } => { + let q2 = q2 as usize; + let out = self.f[q2 * s + from as usize]; + self.g[q * s + cell] = out; + let c2 = c(q2); + if let Some(force) = self.forces.get_mut(group as usize) { + for a in 0..3 { + force[a] += c2[a] * out - cq[a] * out; + } + } + } + // A placeholder the open boundary replaces. + Rule::Open => self.g[q * s + cell] = self.f[opposite(q) * s + cell], + } + } + } + + /// Zou–He velocity and pressure boundaries on the streamed populations. + fn open_boundaries(&mut self) { + let s = self.stride; + let a_half = self.acceleration.map(|v| 0.5 * v); + for oc in &self.open { + let cell = oc.cell as usize; + let axis = oc.axis; + let sign = oc.inward; + let ci = |q: usize| f64::from(C[q][axis]); + // Tangential (c_n = 0) and outgoing (c_n = −s, known) sums. + let (mut tangential, mut outgoing) = (0.0, 0.0); + for q in 0..Q { + let v = self.g[q * s + cell]; + let cn = ci(q) * sign; + if cn == 0.0 { + tangential += v; + } else if cn < 0.0 { + outgoing += v; + } + } + // The bare velocity (without the half force) the populations must carry. + let (rho, bare) = match oc.target { + Target::Velocity(u) => { + let bare = [u[0] - a_half[0], u[1] - a_half[1], u[2] - a_half[2]]; + let un = sign * bare[axis]; + ((tangential + 2.0 * outgoing) / (1.0 - un), bare) + } + Target::Density(rho) => { + let un = 1.0 - (tangential + 2.0 * outgoing) / rho; + let mut bare = [-a_half[0], -a_half[1], -a_half[2]]; + bare[axis] = sign * un; + (rho, bare) + } + }; + // Unknowns: bounce-back of the non-equilibrium part. + for q in (1..Q).filter(|q| ci(*q) * sign > 0.0) { + let cq = c(q); + let cu = cq[0] * bare[0] + cq[1] * bare[1] + cq[2] * bare[2]; + self.g[q * s + cell] = self.g[opposite(q) * s + cell] + 6.0 * W[q] * rho * cu; + } + // Transverse momentum: correct the two diagonal unknowns along each tangential axis. + for t in (0..3).filter(|t| *t != axis) { + let mut m = -rho * bare[t]; + for q in 0..Q { + m += self.g[q * s + cell] * f64::from(C[q][t]); + } + for q in (1..Q).filter(|q| ci(*q) * sign > 0.0 && C[*q][t] != 0) { + self.g[q * s + cell] -= 0.5 * m * f64::from(C[q][t]); + } + } + } + } +} diff --git a/crates/lbm/src/lattice.rs b/crates/lbm/src/lattice.rs new file mode 100644 index 0000000..7ef6648 --- /dev/null +++ b/crates/lbm/src/lattice.rs @@ -0,0 +1,141 @@ +//! The D3Q19 velocity set (Qian, d'Humières & Lallemand 1992; Krüger et al. 2017, ch. 3). +//! +//! Directions come in opposite pairs: 1–2, 3–4, …, 17–18, so the opposite of an odd `q` is +//! `q + 1` and of an even `q` (not 0) is `q − 1`. + +/// The number of velocities. +pub const Q: usize = 19; + +/// The lattice velocities `c_q` in lattice units. +pub const C: [[i32; 3]; Q] = [ + [0, 0, 0], + [1, 0, 0], + [-1, 0, 0], + [0, 1, 0], + [0, -1, 0], + [0, 0, 1], + [0, 0, -1], + [1, 1, 0], + [-1, -1, 0], + [1, -1, 0], + [-1, 1, 0], + [1, 0, 1], + [-1, 0, -1], + [1, 0, -1], + [-1, 0, 1], + [0, 1, 1], + [0, -1, -1], + [0, 1, -1], + [0, -1, 1], +]; + +/// The weights `w_q`: 1/3 at rest, 1/18 along the axes, 1/36 along the face diagonals. +pub const W: [f64; Q] = [ + 1.0 / 3.0, + 1.0 / 18.0, + 1.0 / 18.0, + 1.0 / 18.0, + 1.0 / 18.0, + 1.0 / 18.0, + 1.0 / 18.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, +]; + +/// The squared speed of sound in lattice units, `c_s² = 1/3`. +pub const CS2: f64 = 1.0 / 3.0; + +/// The direction opposite `q`. +#[inline(always)] +pub const fn opposite(q: usize) -> usize { + if q == 0 { + 0 + } else if q % 2 == 1 { + q + 1 + } else { + q - 1 + } +} + +/// `c_q` as floating-point numbers. +#[inline(always)] +pub fn c(q: usize) -> [f64; 3] { + let v = C.get(q).copied().unwrap_or([0; 3]); + [f64::from(v[0]), f64::from(v[1]), f64::from(v[2])] +} + +/// The second-order equilibrium `w_q ρ (1 + 3 c·u + 9/2 (c·u)² − 3/2 u²)`. +#[inline(always)] +pub fn equilibrium(q: usize, rho: f64, u: [f64; 3]) -> f64 { + let c = c(q); + let cu = c[0] * u[0] + c[1] * u[1] + c[2] * u[2]; + let uu = u[0] * u[0] + u[1] * u[1] + u[2] * u[2]; + W.get(q).copied().unwrap_or(0.0) * rho * (1.0 + 3.0 * cu + 4.5 * cu * cu - 1.5 * uu) +} + +/// Ladd's moving-wall term for a population bounced back into direction `q` off a wall moving +/// with `u` (lattice units, reference density 1): `2 w_q c_q·u / c_s²` (Ladd 1994). +#[inline(always)] +pub fn moving_wall(q: usize, u: [f64; 3]) -> f64 { + let c = c(q); + 2.0 * W.get(q).copied().unwrap_or(0.0) * (c[0] * u[0] + c[1] * u[1] + c[2] * u[2]) / CS2 +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn moments_of_the_velocity_set() { + // Σ w = 1, Σ w c = 0, Σ w c_a c_b = c_s² δ_ab, and opposite directions pair up. + let mut m0 = 0.0; + let mut m1 = [0.0; 3]; + let mut m2 = [[0.0; 3]; 3]; + for q in 0..Q { + let c = c(q); + m0 += W[q]; + for a in 0..3 { + m1[a] += W[q] * c[a]; + for b in 0..3 { + m2[a][b] += W[q] * c[a] * c[b]; + } + } + let o = opposite(q); + assert_eq!(opposite(o), q); + assert_eq!([C[o][0], C[o][1], C[o][2]], [-C[q][0], -C[q][1], -C[q][2]]); + assert_eq!(W[o], W[q]); + } + assert!((m0 - 1.0f64).abs() < 1e-15); + for a in 0..3 { + assert!(m1[a].abs() < 1e-15); + for b in 0..3 { + let want = if a == b { CS2 } else { 0.0 }; + assert!((m2[a][b] - want).abs() < 1e-15); + } + } + // The equilibrium carries the density and momentum it is built from. + let (rho, u) = (1.03, [0.02, -0.01, 0.005]); + let (mut r, mut j) = (0.0, [0.0; 3]); + for q in 0..Q { + let f = equilibrium(q, rho, u); + r += f; + for a in 0..3 { + j[a] += f * c(q)[a]; + } + } + assert!((r - rho).abs() < 1e-15); + for a in 0..3 { + assert!((j[a] - rho * u[a]).abs() < 1e-15); + } + } +} diff --git a/crates/lbm/src/lib.rs b/crates/lbm/src/lib.rs new file mode 100644 index 0000000..e45a5ea --- /dev/null +++ b/crates/lbm/src/lib.rs @@ -0,0 +1,107 @@ +//! ModelCraft's lattice Boltzmann flow solver (CPU). +//! +//! It solves a fluid step of the model ([`modelcraft_model::Procedure::Fluid`]) in its fluid +//! domain: weakly compressible, isothermal, Newtonian flow on a uniform grid of cubic cells that +//! fills the domain's box, less the cells inside the obstacle bodies. +//! +//! - [`voxel`]: which cells the obstacles take, decided with exact predicates from the +//! face-tagged boundary surface, and which face each wall link meets. +//! - [`setup`]: the step in lattice units, the time step, and the stability limits. +//! - [`flow`]: the D3Q19 two-relaxation-time scheme with Guo forcing, halfway bounce-back +//! walls (moving walls after Ladd), slip faces, Zou–He inlets and outlets, and the +//! momentum-exchange force on every wall. +//! - [`solve`]: the time loop, output frames and force histories in model units. +//! +//! # Sources +//! - T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, *The Lattice +//! Boltzmann Method: Principles and Practice*, Springer 2017 (velocity sets, equilibria, units, +//! stream–collide, bounce-back, forcing, boundary conditions). +//! - Y. H. Qian, D. d'Humières, P. Lallemand, *Lattice BGK models for Navier–Stokes equation*, +//! EPL 17(6):479, 1992 (D3Q19). +//! - I. Ginzburg, F. Verhaeghe, D. d'Humières, *Two-relaxation-time lattice Boltzmann scheme: +//! about parametrization, velocity, pressure and mixed boundary conditions*, Commun. Comput. +//! Phys. 3(2):427–478, 2008 (TRT and the parameter Λ). +//! - Z. Guo, C. Zheng, B. Shi, *Discrete lattice effects on the forcing term in the lattice +//! Boltzmann method*, PRE 65:046308, 2002 (body force). +//! - A. J. C. Ladd, *Numerical simulations of particulate suspensions via a discretized Boltzmann +//! equation. Part 1*, JFM 271:285–309, 1994 (bounce-back, moving walls, momentum exchange). +//! - R. Mei, D. Yu, W. Shyy, L.-S. Luo, *Force evaluation in the lattice Boltzmann method +//! involving curved geometry*, PRE 65:041203, 2002 (momentum exchange). +//! - Q. Zou, X. He, *On pressure and velocity boundary conditions for the lattice Boltzmann BGK +//! model*, Phys. Fluids 9(6):1591–1598, 1997; M. Hecht, J. Harting, J. Stat. Mech. 2010, +//! P01018 (the D3Q19 transverse correction). +//! - P. Lallemand, L.-S. Luo, *Theory of the lattice Boltzmann method: dispersion, dissipation, +//! isotropy, Galilean invariance, and stability*, PRE 61:6546, 2000 (stability limits). +//! - H. Edelsbrunner, E. P. Mücke, *Simulation of Simplicity*, ACM TOG 9(1):66–104, 1990; +//! J. R. Shewchuk 1997 (exact predicates, through `modelcraft_mesh::predicates`). +//! +//! No code of any lattice Boltzmann package was read; everything is written from these papers and +//! the book. +//! +//! # Determinism +//! Double precision, fixed loop orders, no threads and no parallel reductions: the same model +//! gives the same result bit for bit on every platform, wasm included. +#![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] +#![forbid(unsafe_code)] +// Population arrays are indexed by cell and direction counts fixed when they are allocated. +#![allow(clippy::needless_range_loop)] + +pub mod flow; +pub mod lattice; +pub mod setup; +pub mod solve; +pub mod voxel; + +pub use flow::{Flow, Group}; +pub use setup::{Obstacle, Setup, fluid_step}; +pub use solve::{FluidResult, FluidStats, ForceHistory, Frame, run, solve_fluid}; +pub use voxel::Grid; + +#[derive(Debug, Clone, PartialEq, thiserror::Error)] +pub enum Error { + #[error("invalid input: {0}")] + Invalid(String), + #[error("not found: {0}")] + NotFound(String), + #[error("limit reached: {0}")] + Limit(String), + /// The run diverged. + #[error("unstable: {0}")] + Unstable(String), + #[error(transparent)] + Model(#[from] modelcraft_model::Error), + #[error(transparent)] + Mesh(#[from] modelcraft_mesh::Error), +} + +pub type Result = std::result::Result; + +/// Caps on a run's size and work; a run beyond any of them is refused before it starts. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct Limits { + pub max_cells: u64, + pub max_steps: u64, + /// Fluid cells times time steps. + pub max_cell_updates: u64, + /// The memory the output frames may take. + pub max_frame_bytes: u64, +} + +impl Limits { + /// 16 million cells: about 5 GB of populations in double precision. + pub const MAX_CELLS: u64 = 16_000_000; + pub const MAX_STEPS: u64 = 10_000_000; + pub const MAX_CELL_UPDATES: u64 = 1_000_000_000_000; + pub const MAX_FRAME_BYTES: u64 = 2_000_000_000; +} + +impl Default for Limits { + fn default() -> Self { + Limits { + max_cells: Limits::MAX_CELLS, + max_steps: Limits::MAX_STEPS, + max_cell_updates: Limits::MAX_CELL_UPDATES, + max_frame_bytes: Limits::MAX_FRAME_BYTES, + } + } +} diff --git a/crates/lbm/src/setup.rs b/crates/lbm/src/setup.rs new file mode 100644 index 0000000..922de05 --- /dev/null +++ b/crates/lbm/src/setup.rs @@ -0,0 +1,334 @@ +//! From a fluid step of the model to lattice units: the grid, the time step, the relaxation +//! times, and the checks that keep the scheme in its stable and accurate range. +//! +//! # Units (Krüger et al. 2017, ch. 7) +//! Lattice units take the spacing `Δx`, the time step `Δt` and the fluid's density `ρ₀` as 1. A +//! velocity converts by `Δx/Δt`, a kinematic viscosity by `Δx²/Δt`, an acceleration by `Δx/Δt²`, +//! a pressure by `ρ₀ Δx²/Δt²` and a force by `ρ₀ Δx⁴/Δt²`. The lattice viscosity fixes the +//! symmetric relaxation time, `ν = c_s² (τ⁺ − ½)` with `c_s² = 1/3`. +//! +//! # Time step +//! Unless the step gives one, `Δt` is the largest that keeps the characteristic lattice velocity +//! at [`TARGET_LATTICE_VELOCITY`] (Mach 0.087: the weakly compressible error is O(Ma²)) and +//! `τ⁺` at most [`AUTO_TAU_MAX`] (beyond it the scheme only gets less accurate in time). The +//! characteristic velocity is the largest of the inlet and wall velocities, the initial velocity, +//! `√(2 Δp / ρ)` across pressure faces and, for a driving acceleration `a`, the smaller of the +//! viscous plane-channel peak `a H² / (8 ν)` (H the box's smallest extent) and the free-fall speed +//! `√(2 a L)` (L its largest). +//! +//! # Limits +//! A run is refused, naming the numbers and what to change, when the lattice Mach number of the +//! characteristic velocity exceeds [`MAX_MACH`] (the scheme is weakly compressible and loses +//! stability toward Ma ≈ 0.5), when `τ⁺ < `[`MIN_TAU`] (the grid cannot resolve this viscosity at +//! this velocity; linear stability analysis, Lallemand & Luo 2000, shows the stable velocity range +//! shrinks to nothing as τ → ½), or when the lattice acceleration exceeds [`MAX_LATTICE_ACCELERATION`]. +//! While running, a lattice speed above [`MAX_LATTICE_SPEED`] or a non-finite or non-positive +//! density stops the run with [`crate::Error::Unstable`]. + +use modelcraft_mesh::surface::BoundarySurface; +use modelcraft_model::{BoxFace, FaceCondition, FluidDomain, Model, Procedure, SurfaceCondition}; + +use crate::voxel::Grid; +use crate::{Error, Limits, Result}; + +/// The characteristic lattice velocity the automatic time step aims at. +pub const TARGET_LATTICE_VELOCITY: f64 = 0.05; +/// The largest symmetric relaxation time the automatic time step allows. +pub const AUTO_TAU_MAX: f64 = 1.0; +/// The largest lattice Mach number `u Δt / (Δx c_s)` of the characteristic velocity. +pub const MAX_MACH: f64 = 0.3; +/// The smallest symmetric relaxation time τ⁺. +pub const MIN_TAU: f64 = 0.505; +/// The largest lattice acceleration `a Δt² / Δx`. +pub const MAX_LATTICE_ACCELERATION: f64 = 0.01; +/// A lattice speed above this during a run means it went unstable (Mach 0.87). +pub const MAX_LATTICE_SPEED: f64 = 0.5; +/// The default two-relaxation-time parameter Λ: bounce-back walls exactly halfway between cells +/// (Ginzburg, Verhaeghe & d'Humières 2008). +pub const DEFAULT_MAGIC: f64 = 3.0 / 16.0; + +/// The surface of an obstacle body and the persistent name of each of its faces (by the +/// triangles' face tags); prepared by the caller from the geometry. +#[derive(Debug, Clone)] +pub struct Obstacle { + pub body: String, + pub surface: BoundarySurface, + /// Face names by tag; tags without a name are called `face0`, `face1`, …. + pub faces: Vec, +} + +/// An obstacle of the domain as the run sees it. +#[derive(Debug, Clone, PartialEq)] +pub struct ObstacleFaces { + /// The obstacle's index in the caller's obstacle list. + pub index: usize, + /// Every face tag's name. + pub names: Vec, + /// Every face tag's wall velocity in lattice units (zero for walls at rest). + pub velocities: Vec<[f64; 3]>, +} + +/// A fluid step in lattice units. +#[derive(Debug, Clone, PartialEq)] +pub struct Setup { + pub step: String, + pub domain: FluidDomain, + pub grid: Grid, + /// ρ₀ and the kinematic viscosity ν = μ / ρ₀, in model units. + pub density: f64, + pub viscosity: f64, + pub time_step: f64, + pub steps: u64, + pub frames: u32, + /// τ⁺ and τ⁻ (Λ = (τ⁺ − ½)(τ⁻ − ½)). + pub tau_plus: f64, + pub tau_minus: f64, + pub magic: f64, + /// The characteristic velocity in model and lattice units, and its lattice Mach number. + pub velocity: f64, + pub lattice_velocity: f64, + pub mach: f64, + /// The body acceleration and initial velocity in lattice units. + pub acceleration: [f64; 3], + pub initial_velocity: [f64; 3], + /// The domain's obstacles, in its order. + pub obstacles: Vec, + pub warnings: Vec, +} + +fn norm(v: [f64; 3]) -> f64 { + (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt() +} + +/// The fluid step named `step` (or the only fluid step when `step` is `None`) and its domain. +pub fn fluid_step<'m>(model: &'m Model, step: Option<&str>) -> Result<(&'m modelcraft_model::Step, &'m FluidDomain)> { + let s = model.select_step(step, |p| p.fluid_domain().is_some(), "fluid")?; + let name = s.procedure.fluid_domain().unwrap_or_default(); + let domain = model.fluid_domain(name).ok_or_else(|| Error::NotFound(format!("fluid domain `{name}`")))?; + Ok((s, domain)) +} + +impl Setup { + /// Convert a fluid step, checking it against the stability limits and `limits`. + pub fn new(model: &Model, step: Option<&str>, obstacles: &[Obstacle], limits: &Limits) -> Result { + // The model checked the domain and its material when they were defined. + let (s, domain) = fluid_step(model, step)?; + let domain = domain.clone(); + let Procedure::Fluid { end_time, steps, time_step, frames, magic, .. } = &s.procedure else { + return Err(Error::Invalid(format!("step `{}` is not a fluid step", s.name))); + }; + let material = model.material(&domain.material); + let (Some(density), Some(fluid)) = (material.and_then(|m| m.density), material.and_then(|m| m.fluid)) else { + return Err(Error::Invalid(format!("material `{}` needs a density and fluid behaviour", domain.material))); + }; + let viscosity = fluid.viscosity() / density; + if !(viscosity > 0.0 && viscosity.is_finite()) { + return Err(Error::Invalid(format!("the kinematic viscosity μ/ρ of material `{}` is {viscosity}", domain.material))); + } + + // The grid. + let counts = domain.cells()?; + let mut cells = [0usize; 3]; + let mut total: u64 = 1; + for a in 0..3 { + cells[a] = usize::try_from(counts[a]).map_err(|_| Error::Limit("the grid is too large".into()))?; + total = total.saturating_mul(counts[a]); + } + if total > limits.max_cells { + return Err(Error::Limit(format!( + "the domain has {total} cells ({} × {} × {}), above the cap of {}; coarsen the spacing or raise `max_cells`", + cells[0], cells[1], cells[2], limits.max_cells + ))); + } + for face in BoxFace::ALL { + let open = matches!(domain.faces.get(face), FaceCondition::Velocity { .. } | FaceCondition::Pressure { .. }); + if open && cells[face.axis()] < 2 { + return Err(Error::Invalid(format!("box face {} is an inlet or outlet, so the box needs at least 2 cells across it", face.name()))); + } + } + let grid = Grid { origin: domain.min, spacing: domain.spacing, cells }; + let h = domain.spacing; + + // The characteristic velocity. + let mut u_char: f64 = norm(domain.initial_velocity); + let mut pressures = Vec::new(); + for face in BoxFace::ALL { + match *domain.faces.get(face) { + FaceCondition::MovingWall { velocity } | FaceCondition::Velocity { velocity, .. } => u_char = u_char.max(norm(velocity)), + FaceCondition::Pressure { pressure } => pressures.push(pressure), + FaceCondition::Wall {} | FaceCondition::Slip {} | FaceCondition::Periodic {} => {} + } + } + for sb in &domain.surfaces { + if let SurfaceCondition::MovingWall { velocity } = sb.condition { + u_char = u_char.max(norm(velocity)); + } + } + let lo = pressures.iter().copied().fold(f64::INFINITY, f64::min); + let hi = pressures.iter().copied().fold(f64::NEG_INFINITY, f64::max); + if hi > lo { + u_char = u_char.max((2.0 * (hi - lo) / density).sqrt()); + } + let a = norm(domain.acceleration); + if a > 0.0 { + let extents: Vec = (0..3).map(|k| domain.max[k] - domain.min[k]).collect(); + let small = extents.iter().copied().fold(f64::INFINITY, f64::min); + let large = extents.iter().copied().fold(0.0, f64::max); + u_char = u_char.max((a * small * small / (8.0 * viscosity)).min((2.0 * a * large).sqrt())); + } + if !u_char.is_finite() { + return Err(Error::Invalid("the flow's velocity scale is not a finite number".into())); + } + + // The time step and the number of steps. + let diffusive = (AUTO_TAU_MAX - 0.5) * h * h / (3.0 * viscosity); + let auto = if u_char > 0.0 { (TARGET_LATTICE_VELOCITY * h / u_char).min(diffusive) } else { diffusive }; + let (dt, n) = match (*time_step, *end_time, *steps) { + (Some(dt), Some(t), _) => { + let n = (t / dt).round(); + if !(n >= 1.0) { + return Err(Error::Invalid(format!("the end time {t} is shorter than the time step {dt}"))); + } + (dt, n) + } + (Some(dt), None, Some(n)) => (dt, n as f64), + (None, Some(t), _) => { + let n = (t / auto).ceil().max(1.0); + (t / n, n) + } + (None, None, Some(n)) => (auto, n as f64), + (_, None, None) => return Err(Error::Invalid(format!("step `{}` needs `end_time` or `steps`", s.name))), + }; + if !(dt > 0.0 && dt.is_finite()) { + return Err(Error::Invalid(format!("the time step {dt} is not a positive number"))); + } + if !(n <= limits.max_steps as f64) { + return Err(Error::Limit(format!( + "the step needs {n} time steps of {dt}, above the cap of {}; shorten it or raise `max_steps`", + limits.max_steps + ))); + } + let n = n as u64; + let updates = (total as f64) * (n as f64); + if !(updates <= limits.max_cell_updates as f64) { + return Err(Error::Limit(format!( + "{total} cells for {n} time steps is {updates:e} cell updates, above the cap of {:e}; coarsen, shorten or raise `max_cell_updates`", + limits.max_cell_updates as f64 + ))); + } + + // Lattice parameters and the stability limits. + let nu = viscosity * dt / (h * h); + let tau_plus = 3.0 * nu + 0.5; + let lattice_velocity = u_char * dt / h; + let mach = lattice_velocity / crate::lattice::CS2.sqrt(); + if mach > MAX_MACH { + return Err(Error::Invalid(format!( + "the lattice Mach number is {mach:.3} (velocity {u_char} with time step {dt} and spacing {h}), above {MAX_MACH}: lower the time step" + ))); + } + if !(tau_plus >= MIN_TAU) { + return Err(Error::Invalid(format!( + "the relaxation time τ = {tau_plus:.5} is below {MIN_TAU}: the grid cannot resolve viscosity {viscosity} at this time step (cell Reynolds number {:.1}); refine the spacing or raise the time step", + if nu > 0.0 { lattice_velocity / nu } else { f64::INFINITY } + ))); + } + let magic = magic.unwrap_or(DEFAULT_MAGIC); + let tau_minus = 0.5 + magic / (tau_plus - 0.5); + let acceleration = domain.acceleration.map(|c| c * dt * dt / h); + if norm(acceleration) > MAX_LATTICE_ACCELERATION { + return Err(Error::Invalid(format!( + "the lattice acceleration is {:.3e}, above {MAX_LATTICE_ACCELERATION}: lower the time step", + norm(acceleration) + ))); + } + + // Obstacles, their face names and their faces' wall velocities. + let mut domain_obstacles = Vec::with_capacity(domain.obstacles.len()); + let mut warnings = Vec::new(); + for body in &domain.obstacles { + let (index, o) = obstacles + .iter() + .enumerate() + .find(|(_, o)| &o.body == body) + .ok_or_else(|| Error::NotFound(format!("obstacle body `{body}` has no surface (open the design or import its surface)")))?; + let tags = o.surface.tri_face().iter().copied().max().map_or(1, |m| m as usize + 1).max(o.faces.len()); + let names: Vec = (0..tags).map(|t| o.faces.get(t).cloned().unwrap_or_else(|| format!("face{t}"))).collect(); + let mut velocities = vec![[0.0; 3]; tags]; + for sb in domain.surfaces.iter().filter(|sb| &sb.face.body == body) { + let mut hit = false; + for (tag, name) in names.iter().enumerate() { + if sb.face.selects(&modelcraft_model::GeometryRef::face(body.as_str(), name.as_str())?) { + hit = true; + if let (Some(slot), SurfaceCondition::MovingWall { velocity }) = (velocities.get_mut(tag), sb.condition) { + *slot = velocity.map(|c| c * dt / h); + } + } + } + if !hit { + return Err(Error::NotFound(format!("the {} is not a face of the obstacle's surface", sb.face))); + } + } + domain_obstacles.push(ObstacleFaces { index, names, velocities }); + } + if mach > 0.1 + 1e-12 { + warnings.push(format!( + "lattice Mach number {mach:.3}: the compressibility error, of order Ma² ≈ {:.1}%, may be visible", + 100.0 * mach * mach + )); + } + if tau_plus < 0.51 { + warnings.push(format!("τ = {tau_plus:.4} is close to ½: under-resolved flows may go unstable; a finer spacing is safer")); + } + + Ok(Setup { + step: s.name.clone(), + grid, + density, + viscosity, + time_step: dt, + steps: n, + frames: *frames, + tau_plus, + tau_minus, + magic, + velocity: u_char, + lattice_velocity, + mach, + acceleration, + initial_velocity: domain.initial_velocity.map(|c| c * dt / h), + obstacles: domain_obstacles, + warnings, + domain, + }) + } + + /// Model units per lattice unit of velocity, pressure and force. + pub fn velocity_scale(&self) -> f64 { + self.grid.spacing / self.time_step + } + + pub fn pressure_scale(&self) -> f64 { + self.density * self.grid.spacing * self.grid.spacing / (self.time_step * self.time_step) + } + + pub fn force_scale(&self) -> f64 { + self.pressure_scale() * self.grid.spacing * self.grid.spacing + } + + /// A velocity in lattice units. + pub fn to_lattice(&self, v: [f64; 3]) -> [f64; 3] { + let k = 1.0 / self.velocity_scale(); + v.map(|c| c * k) + } + + /// The lattice density of a gauge pressure: `1 + p / (c_s² ρ₀ Δx²/Δt²)`. + pub fn lattice_density(&self, pressure: f64) -> f64 { + 1.0 + pressure / (crate::lattice::CS2 * self.pressure_scale()) + } + + /// The gauge pressure of a lattice density (the inverse of [`Setup::lattice_density`]). + pub fn gauge_pressure(&self, rho: f64) -> f64 { + (rho - 1.0) * crate::lattice::CS2 * self.pressure_scale() + } +} diff --git a/crates/lbm/src/solve.rs b/crates/lbm/src/solve.rs new file mode 100644 index 0000000..e209155 --- /dev/null +++ b/crates/lbm/src/solve.rs @@ -0,0 +1,190 @@ +//! Running a fluid step: the time loop, output frames and force histories, in model units. + +use modelcraft_model::Model; + +use crate::flow::Flow; +use crate::setup::{Obstacle, Setup}; +use crate::voxel::Grid; +use crate::{Error, Limits, Result}; + +/// Most samples a force history keeps; longer runs record every k-th step. +pub const MAX_HISTORY: u64 = 10_000; + +/// A relative density variation above this at the end of a run is reported. +pub const COMPRESSIBILITY_WARNING: f64 = 0.01; + +/// The velocity and pressure field at one time, per cell (index `i + nx (j + ny k)`; zero in +/// solid cells). Stored in single precision: output for display and probing. +#[derive(Debug, Clone, PartialEq)] +pub struct Frame { + pub time: f64, + pub step: u64, + pub velocity: Vec<[f32; 3]>, + /// Gauge pressure. + pub pressure: Vec, + /// The total fluid mass (density times volume). + pub mass: f64, +} + +/// The force the fluid exerts on one wall over time. +#[derive(Debug, Clone, PartialEq)] +pub struct ForceHistory { + /// `x_min` … for box faces, the face's persistent name for obstacle faces. + pub name: String, + /// The obstacle body, for obstacle faces. + pub body: Option, + /// One value per [`FluidResult::history_times`]. + pub force: Vec<[f64; 3]>, +} + +/// Run statistics. +#[derive(Debug, Clone, PartialEq)] +pub struct FluidStats { + pub cells: usize, + pub fluid_cells: usize, + pub steps: u64, + pub time_step: f64, + pub end_time: f64, + pub tau_plus: f64, + pub tau_minus: f64, + pub magic: f64, + /// The characteristic velocity (model units), its lattice value and its lattice Mach number. + pub velocity: f64, + pub lattice_velocity: f64, + pub mach: f64, + /// Fluid cells times steps. + pub cell_updates: u64, +} + +/// The result of a fluid step. +#[derive(Debug, Clone, PartialEq)] +pub struct FluidResult { + pub step: String, + pub grid: Grid, + /// Per cell: not fluid (inside an obstacle). + pub solid: Vec, + /// The initial field and `frames` more, evenly spaced to the end. + pub frames: Vec