OTOM Logo

OTOM CFD WebGPU

Advanced Simulation, 3D Raymarching, & AI Orchestration

Table of Contents

1. System Overview & Architecture

OTOM CFD is a high-performance computational fluid dynamics solver running entirely in the browser via WebGPU. It is heavily inspired by the methodologies popularized by the WaterLily.jl Julia package, specifically bringing Cartesian-grid, boundary-data immersion methods to massive parallel compute shaders.

The software supports both 2D and 3D incompressible Navier-Stokes simulations, a dedicated 2D compressible Euler engine (HLLC Riemann solver) for shock-bearing high-speed flows, multi-phase flows, and high-fidelity volumetric raymarching—all orchestrated through an interactive GUI and an integrated LLM-powered AI assistant.

2. Comprehensive Parameter Guide

The simulation is driven by a uniform buffer passed to the GPU every frame. Parameters are split between physical and visual properties.

ParameterDomainDescription / Physical Meaning
Engine Solver2DIncompressible (projection Navier-Stokes, low speed) or Compressible (Euler/HLLC, shocks & high speed). See §4.
Grid Resolution2D / 3DNumber of computational cells (e.g., 910x512 in 2D, 64³ to 256³ in 3D). Dictates spatial frequency limits.
Flow Mode2D / 3DVelocity enforces a constant inlet speed. Pressure drives flow via boundary pressure gradients (\(\Delta P\)).
Flow Speed (\(U_\infty\))2D / 3DInlet velocity in grid cells per second (or m/s if SI units enabled).
Lattice dt (\(\Delta t\))2D / 3DTime advanced per sub-step. Governs the CFL condition (\(CFL = U \frac{\Delta t}{\Delta x}\)).
Physics Steps2D / 3DNumber of compute passes executed per visual render frame. Increases stability for high-speed flows.
Jacobi Itr2D / 3DIterations of the Poisson pressure solver. Higher values ensure stricter incompressibility (\(\nabla \cdot \mathbf{u} = 0\)).
Viscosity (\(\nu\))2D / 3DKinematic viscosity. Defines fluid thickness. Air ~0.015, Water ~0.001.
Visual ParameterDomainDescription
View Mode2D / 3DSelects the Eulerian field to display: Vorticity, Velocity, Pressure, Phase, or Schlieren (Pressure Gradient).
Gain / VScale2D / 3DMultiplicative factor applied before colormapping. Adjusts the sensitivity of the visualizer.
Theme2D / 3DColormap applied to the main viewport field: Seismic, Smoky Seismic, Jet, Hot, Viridis, Turbo, Magma, Inferno, Plasma. See below.
Isoline Density2D / 3DFrequency of contour lines superimposed on the field to visualize gradients.
Arrow Density2D / 3DSpacing of the vector field arrows. In 3D, this samples a specific 2D slice defined by Slice Depth.

Colormap Themes

Jet and Hot are visually similar names but genuinely different palettes: Jet is the classic MATLAB-style rainbow spectrum — dark blue → blue → cyan → green → yellow → red → dark red, spanning the field's full normalized range — while Hot is a black/red/yellow ramp driven by the signed field value rather than a spatial position. Both live in the same shared colormap() WGSL function, so any 2D or 3D view (main field, Force Display overlay, particle color-by-field) can use either. Force Theme (Aerodynamics tab, 2D and 3D) is a second, independent colormap selector for the pressure/force overlay — it can be set differently from the main Theme, e.g. Seismic for the flow field with Jet reserved for the Force Display, so the two visual channels never fight for the same color language.

Colormap Window, Auto Range & Colour Scale

Colormap Min/Max defines the window mapped onto the colormap's two ends — field values at Min land on the cold end, values at Max on the hot end. The historical fixed mapping is just the default window (\(-1\), \(1\)); many fields don't naturally sit in that range (an incompressible pressure field, for instance, floats at an arbitrary additive gauge offset with no fixed natural range at all), so a mismatched window can make a field with real structure look like a single flat colour.

Auto Range removes the guesswork: while enabled, the app periodically measures the field actually on screen (its 1st/99th percentile, robust against a single outlier cell) and fits Colormap Min/Max — and, in 3D, the render Threshold — to it automatically. The read-only Measured line shows the fitted percentiles alongside the field's true min/max/mean — useful on its own for telling a field that "looks flat" (genuinely tiny values) apart from one that's merely mis-scaled (a large mean with real but small variation riding on top of it, the exact signature of the pressure gauge-offset case above).

Colour Scale (Linear or Log) is a separate, complementary control: it doesn't change the window, it redistributes contrast within it. Log applies a signed, sign-preserving logarithmic curve after the Min/Max window is applied, so weak structure near the middle of the range becomes visible without clipping the strong values at the ends — the window's two ends still map to exactly \(\pm 1\) either way. Fit the window with Auto Range first, then switch to Log if fine detail is still hard to see.

Both 2D and 3D pressure sampling honour this pipeline identically, and in 3D compressible mode the window/scale apply to the same decoded velocity/pressure textures the shared volume renderer already uses — Auto Range and Log work in the Compressible engine exactly as they do in Incompressible.

Physical Calibration & SI Units 2D3D

By default, every number in the sidebar is a lattice unit — self-consistent relative to the grid, not inherently meters or seconds. Flipping SI Units Mode (2D: "Physical Calibration" folder; 3D: its own "Physical Calibration" folder, right below Simulation — the two are independent, so turning it on in one tab never silently recalibrates the other) reinterprets Flow speed/Viscosity as real physical quantities (m/s, m²/s) and converts them into the lattice units the solver actually runs on:

$$u_{lattice} = u_{SI} \cdot \frac{\Delta t_{SI}}{\Delta x}, \qquad \nu_{lattice} = \nu_{SI} \cdot \frac{\Delta t_{SI}}{\Delta x^2}, \qquad \Delta x = \frac{\text{Domain Width (m)}}{\text{Resolution}}$$

where Domain Width (m) and Real-world Time Step (s) are the two calibration inputs you provide, and Lattice dt is fixed to 1.0 internally once SI mode is on (the real timestep is already folded into the conversion above). In 3D, Domain Width (m) is calibrated for Length (X) = 1.0 (the domain's default ±1 world extent) and automatically rescales if you change the Length (X) slider, so resizing the non-cubic control volume doesn't silently invalidate your calibration. Height (Y)/Depth (Z) can still differ from Length(X) (the domain doesn't have to be a cube) — each axis then has a physically different cell size, which the solver's diffusion/gradient stencils don't individually correct for (they operate uniformly in normalized grid-cell space); this only matters if your domain is significantly non-cubic and you need diffusion to be quantitatively exact along the stretched axis, not just qualitatively correct.

Sanity-checking your inputs: both the 2D and 3D Physical Calibration folders show three read-only numbers — Cell Size (mm) (the physical size of one grid cell), Effective CFL (the Courant number of the actual lattice velocity that reaches the GPU, in either units mode — named "Effective" rather than "CFL Number" to avoid confusion with the Compressible engine's own CFL Number parameter, which is a different thing: a user-set target for the acoustic timestep, not a live diagnostic), and Est. Reynolds (computed from Flow speed, Viscosity, and a characteristic length derived from Domain Width — named "Est." rather than "Reynolds (Re)" to avoid confusion with the Aerodynamics folder's own Reynolds display, which only updates while "Live Analysis" is on; this one is always live). If Effective CFL climbs too high, a warning banner appears above the canvas — lower Flow speed, raise Domain Width (m), or shrink Real-world Time Step. A CFL warning means the fluid is moving more than about one grid cell per step, which produces visibly wrong (oscillating, exploding, or NaN) results — this is the first thing to check before trusting any result once SI units are on.

Turning SI Units Mode on auto-corrects Viscosity for the selected Fluid Material (unless it's set to "Custom"): the lattice-tuned viscosity values (e.g. Air's default 0.015) were chosen to look right as abstract lattice numbers, not as real physical properties — reinterpreting that same number as literal m²/s the instant SI mode turns on would imply an object roughly 1,000× more viscous than real air. To prevent a jarring change in flow behaviour, the app substitutes the real 20°C kinematic viscosity for the current material (Air ≈ 1.48×10⁻⁵ m²/s, Water ≈ 1.0×10⁻⁶, Engine Oil ≈ 1.2×10⁻⁴, Glycerin ≈ 1.18×10⁻³ m²/s) and widens the Viscosity slider's range/precision to represent these much smaller numbers usefully. Toggling SI Units Mode back off restores the original lattice-tuned value. Switching Fluid Material while SI is on picks from the real-value table instead of the lattice one.

Why the flow can still look "strange" right after enabling SI units even with this correction: real physical parameters at a coarse simulation grid very easily imply a Reynolds number in the hundreds of thousands or millions — e.g. real air, a modest few-m/s breeze, and a few-meter domain is already Re in the millions, which no 64³–256³ grid can resolve without looking chaotic and rapidly-changing, LES turbulence modelling or not. This is not a bug: it's the same "under-resolved turbulence" limitation every CFD tool has when asked to simulate a too-high-Re flow on too-coarse a grid. The Est. Reynolds readout crossing roughly 5,000 while SI Units Mode is on triggers a "HIGH REYNOLDS" banner explaining this and suggesting concrete fixes — raise Viscosity (a physically thicker fluid, or accept a coarser numerical/eddy viscosity), lower Flow speed, or shrink the domain (Domain Width and Real-world Time Step alone won't fix this if changed proportionally to each other, since that leaves the implied physical scenario — and therefore the Reynolds number — unchanged).

The Ground Truth tab (§10) remains the recommended way to validate results once your CFL is healthy and Reynolds number is in a resolvable range: it compares measured Cd/Cl against textbook correlations for canonical shapes at the current Reynolds regime.

3. Core Numerics & Theory Solver

The core fluid solver uses a staggered Cartesian grid with a fractional-step (Projection) method.

Advection & Projection

1. Advection: The velocity field is self-advected using a semi-Lagrangian back-tracing scheme. To reduce numerical dissipation, we utilize high-order filtering.
2. Divergence: The divergence of the intermediate velocity field (\(\mathbf{u}^*\)) is calculated.
3. Pressure Poisson: A Jacobi iterative solver resolves the pressure field needed to counteract the divergence.
4. Gradient Subtraction: The velocity is corrected: \(\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t \nabla p\).

Boundary Data Immersion Method (BDIM)

Instead of conforming a mesh to an obstacle, OTOM CFD uses Brinkman Penalization. Solid objects are defined by exact analytical Signed Distance Fields (SDFs). The velocity inside the SDF (\(d \le 0\)) is forced to match the object's velocity, allowing for arbitrary, moving, and overlapping geometries at near-zero computational cost.

4. Compressible Flow Solver Euler / HLLC

Alongside the default incompressible projection solver, OTOM CFD provides a dedicated compressible engine (set Engine Solver = Compressible) for high-speed, transonic, and supersonic flows where density variation and shock waves are physically significant. Instead of enforcing incompressibility, it solves the 2D compressible Euler equations in conservative finite-volume form, so it naturally captures shocks, expansion fans, and contact discontinuities.

Governing Equations

The conserved state stored in every cell is \(\mathbf{Q} = [\rho,\ \rho u,\ \rho v,\ E]^T\) (density, x/y momentum, total energy), evolved by the hyperbolic conservation law:

$$\frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial \mathbf{F}(\mathbf{Q})}{\partial x} + \frac{\partial \mathbf{G}(\mathbf{Q})}{\partial y} = 0$$

The system is closed by the ideal-gas equation of state, where \(\gamma\) is the ratio of specific heats and \(a\) the local speed of sound:

$$p = (\gamma - 1)\left(E - \tfrac{1}{2}\rho(u^2 + v^2)\right), \qquad a = \sqrt{\frac{\gamma p}{\rho}}$$

Godunov Finite-Volume Scheme & HLLC Riemann Solver

Each cell is advanced by the net flux across its four faces. The flux at every interface is obtained by solving a local Riemann problem with the HLLC (Harten–Lax–van Leer–Contact) approximate solver, which restores the central contact/shear wave that the simpler HLL solver smears out—essential for clean slip lines and wakes:

$$\mathbf{Q}^{n+1}_{i,j} = \mathbf{Q}^{n}_{i,j} - \frac{\Delta t}{\Delta x}\!\left(\mathbf{F}_{i+\frac12} - \mathbf{F}_{i-\frac12}\right) - \frac{\Delta t}{\Delta y}\!\left(\mathbf{G}_{j+\frac12} - \mathbf{G}_{j-\frac12}\right)$$

The full stencil—four face fluxes plus the conservative update—is fused into a single WebGPU compute pass for robustness and throughput. The scheme is first-order in space and time; a positivity guard floors density and internal energy each step to guarantee \(\rho>0\) and \(p>0\), preventing the NaN blow-up that otherwise occurs across strong shocks.

Boundary & Obstacle Treatment

  • Inlet (left): Dirichlet free-stream state fixed by the user's Inlet Mach, Inlet Density, and Inlet Pressure.
  • Outlet (right): Zero-gradient (transmissive) extrapolation.
  • Top / bottom: Inviscid slip walls (\(v = 0\)).
  • Obstacles: Treated as inviscid reflecting slip walls. The wall-face flux carries pressure only, \(\mathbf{F}_{wall} = [\,0,\ p\,n_x,\ p\,n_y,\ 0\,]^T\), which produces detached bow shocks ahead of blunt bodies and oblique shocks off wedges.

Stability, Parameters & Visualization

The explicit time step is bounded by the acoustic CFL condition \(\Delta t < \Delta x / (|u| + a)\). If a strong shock destabilizes the field, lower Lattice dt or raise Physics Steps. The engine exposes four parameters (active only in Compressible mode):

ParameterRangePhysical Meaning
Gamma (\(\gamma\))1.1 – 2.0Ratio of specific heats (1.4 for diatomic air).
Inlet Mach0 – 5.0Free-stream Mach number \(M = u/a\). Values > 1 are supersonic.
Inlet Density0.1 – 5.0Free-stream density \(\rho_\infty\).
Inlet Pressure0.1 – 5.0Free-stream static pressure \(p_\infty\).

Recommended view modes: Schlieren (density gradient \(|\nabla \rho|\), the classic shock-visualization technique), Velocity (rendered as Mach number), or Pressure. Lagrangian particles and streamlines are advected directly by the reconstructed velocity \(\mathbf{u} = (\rho u,\ \rho v)/\rho\), so they stream through and bend around shocks just like the field itself.

5. Turbulence & Fidelity LES

Large Eddy Simulation (LES)

The Smagorinsky model calculates sub-grid scale (SGS) stresses. It derives an "eddy viscosity" \(\nu_t\) from the resolved strain rate tensor \(\bar{S}_{ij}\), dissipating energy at the smallest resolvable scales.

$$\nu_{t} = (C_s \Delta)^2 \sqrt{2 \bar{S}_{ij} \bar{S}_{ij}}$$

Vorticity Confinement

To counteract the inherent diffusion of the semi-Lagrangian advection, Vorticity Confinement applies a body force that pushes fluid towards local vorticity maxima, keeping vortices sharp and persistent.

$$\mathbf{f}_{conf} = \epsilon \Delta x (\mathbf{N} \times \omega), \quad \mathbf{N} = \frac{\nabla |\omega|}{|\nabla |\omega||}$$

6. Multi-Phase & Control Volumes

OTOM CFD supports multi-phase simulations (e.g., air and water interaction) by advecting a scalar phase field \(\phi \in [0,1]\).

7. 3D Simulation & Volumetric Rendering WebGPU 3D

The 3D solver extends the 2D logic into a texture_3d environment. Rendering is achieved via a single-pass Raymarching fragment shader.

Physically Based Raymarching

The camera casts rays through the domain bounding box. The shader performs two primary operations per ray:

  1. Sphere Tracing: Finds intersections with 3D solid obstacles using the get_sdf3d function. If hit, it calculates surface normals and applies PBR lighting (Metallic/Roughness/Ambient).
  2. Volume Accumulation: If the ray misses the solid (or before it hits), it steps through the 3D texture, sampling the chosen field (e.g., Vorticity). Color and alpha are accumulated via alpha blending (\(C = C_{src}\alpha + C_{dst}(1-\alpha)\)).

Advanced 3D Features

  • Isosurface Mode: Instead of fuzzy clouds, the volume renderer searches for a specific threshold value, calculates the 3D gradient normal of the fluid field, and renders it as a solid, shaded surface complete with volumetric soft shadows. The solid obstacle body is always composited as an opaque backing layer behind the isosurface, so the geometry remains visible even when high field values (e.g. stagnation zone) exist upstream of the surface.
  • Cutting Planes: Users can define a plane (\(\mathbf{n} \cdot \mathbf{x} - d = 0\)) to clip the volume or strictly render a 2D cross-section slice inside the 3D domain.

Obstacle Spin & Local Axis Gizmo 3D

Every 3D obstacle — the primary analytical shape (Sphere/Box/Cylinder/Wing/Extrude 2D), the STL mesh, and each of up to 4 secondary "Extra Obstacles" — can spin continuously about all 3 axes, the 3D analogue of the 2D Torque slider. Spin X/Y/Z (deg/frame) is a rotation rate: each rendered frame, the obstacle's current absolute orientation angle (Rot X/Y/Z, directly editable — set it to pose a fixed starting tilt with spin left at zero, e.g. pitching a blade before it starts turning) advances by the rate on each axis independently (wrapped to \([-180°,180°]\)), then that angle is fed into the SDF evaluation. Setting a rate to 0 simply holds the current orientation. Spin About Local Axes (per obstacle) applies the three rates in the obstacle's own rotated frame instead of world axes, so a tilted shape keeps spinning about its own axis instead of the components mixing across world axes as it tilts.

Mechanically, the rotation is applied as the exact inverse of the STL transform's own rotation (stl_inv_rotate, \(R^{-1} = R_x(-r_x)R_y(-r_y)R_z(-r_z)\)) to the query point before the primitive's local-space SDF runs, and the resulting distance is otherwise unaffected — so spin composes correctly with the Wing's own Angle-of-Attack slider (AoA still tilts the chord inside the wing's own frame, after the outer spin has already oriented the whole obstacle in world space). Because get_sdf3d/get_sdf_extra are the single shared functions consumed by rendering, particle/streamline collision, and the aerodynamic force integral, a spinning obstacle is mechanically consistent everywhere at once — no separate plumbing was needed for each consumer.

Show Local Axes (a checkbox per obstacle — one for the primary obstacle, one for STL if selected, one for each Extra Obstacle) draws that obstacle's local X (red) / Y (green) / Z (blue) axis as a short line from its centre, rotated by the same spin, so the user can always see which way "local X" currently points as the shape spins.

Compressible engine caveat: spin for the primary obstacle and Extra Obstacles is only applied when Engine Solver is Incompressible. The Compressible (HLLC) solver's own wall boundary condition uses a separate, unrotated copy of the SDF for its compute kernels; rotating the shared render-path SDF without also rotating that copy would make the rendered obstacle visually disagree with what the solver treats as solid, so rotation is forced back to identity while Compressible is active. STL spin is unaffected by this caveat and works in both engines, since the STL transform was already fully mirrored into the compressible solver's SDF copy before this feature existed.

Spin drives real flow — it is not a cosmetic rotation

A spinning obstacle's solid cells carry the rigid-body velocity \(\boldsymbol{\omega} \times \mathbf{r}\) (about the obstacle's own centre) instead of the historical zero — a genuine moving no-slip wall, exactly like the physical surface of a real spinning object. This is what lets a rotating shape actually drag the surrounding fluid, shed vorticity, and pump flow, rather than presenting fluid with a stationary wall that merely changes orientation between frames. The same mechanism drives the swirl/downwash behind a spinning box, sphere, or cylinder, and the suction/pressure-side split across a genuinely pitched blade (see below).

Blade pitch matters for thrust. A symmetric shape spinning about an axis through its own centre — a Box, a Sphere, an unpitched flat plate — mostly produces swirl and a weak centrifugal low-pressure core along the spin axis, not a net thrust-generating pressure differential between its faces. Real thrust needs the blade angled relative to its own rotation plane, exactly like a real propeller's twist: the Wing shape's AoA provides this directly, or an STL mesh with real blade pitch built into its geometry. Measure it with Torque (Aerodynamics tab, below) alongside the existing Drag/Lift: Drag/Lift give the thrust a spinning, pitched blade produces; Torque gives what it cost to produce it (shaft power is Torque × angular speed) — together they're what let a user read a spinning blade's aerodynamic efficiency off the simulation.

Primitive Sizing 3D

Beyond the uniform Scale slider, the primary Box and Cylinder shapes have independent per-axis controls, shown only when that shape is selected: Box Size X/Y/Z (ratios, default 1.0 each, multiplying the historical fixed 0.25 half-extent per axis) let a Box be stretched into a slab or a beam rather than staying a cube; Cylinder Radius / Cylinder Height (ratios, default 1.0, multiplying the historical fixed \(r=0.2\), \(h=0.4\)) do the same for a Cylinder. Both compose with Scale, which still resizes the whole result uniformly afterward.

Extrude 2D takes its footprint from whatever is currently painted on the flat 2D canvas (Draw Mode) — it is not a parametric shape, and the 2D canvas is the actual UI for its outline. Its Extrusion Length slider controls the extruded half-length along Z (ratio, default 1.0 = the historical fixed 0.85 half-span), independent of Scale (which would otherwise also zoom the 2D footprint along with the length).

Extra Obstacles (3D) & STL Instancing

The "Extra Obstacles (3D)" folder lets a scene hold up to 4 secondary obstacles in addition to the primary one, each with an independent Shape, position, scale, spin, and local-axis toggle (added/removed via + Add Obstacle / Remove). The Shape choices are Sphere, Box, Cylinder, Wing, and STL.

STL as an Extra Obstacle is an instance, not a second upload. There is only one voxelized STL SDF texture in the whole scene — the mesh loaded via the primary "Obstacle" folder's Upload STL button. Choosing STL for an Extra Obstacle samples that same voxel data at this slot's own independent position, scale, and rotation, exactly the same technique the primary obstacle itself uses to place the mesh. This is useful for repeating one part at several locations (e.g. an array of identical pillars, or a duplicate of a car body at a different offset) without the cost of a second voxelization pass or a second GPU texture. If no STL file has been uploaded yet, an STL-typed Extra Obstacle simply renders as nothing until one is.

The previous fifth shape option was Ahmed Body — a standard simplified-automotive-aerodynamics benchmark shape (a bluff body with a characteristic slanted rear window, introduced by S.R. Ahmed in 1984, widely used to validate drag-coefficient and wake-structure predictions against wind-tunnel data) — which has been replaced by STL instancing.

Non-Cubic Control Volume (Domain Resize) 3D

The 3D domain is no longer constrained to a unit cube. Three independent half-extent sliders — Length (X), Height (Y), and Depth (Z) — let users reshape the control volume to a rectangular cuboid. Default value 1.0 corresponds to the original ±1 extent in each axis. All internal coordinate transformations (world-to-UV, raymarching box, particle seeding, streamline reset bounds) scale with these parameters via three uniform floats (domainX, domainY, domainZ) pushed to the GPU every frame.

SliderRangeEffect
Length (X)0.5 – 4.0Stretches / compresses the streamwise (flow) direction.
Height (Y)0.5 – 4.0Stretches / compresses the vertical extent.
Depth (Z)0.5 – 4.0Stretches / compresses the spanwise direction.

Changing any dimension triggers a full simulation re-initialisation (initSim()) so that the velocity field and obstacle mask are consistent with the new geometry.

Velocity Perturbation Contour 3D Viz

In View Mode = Velocity, the 3D volume renderer displays the signed perturbation from free-stream rather than raw speed. The formula is:

$$ v_{display} = \left(\frac{|\mathbf{u}|}{U_\infty} - 1\right) \times G $$

where \(G\) is the user's Gain (VScale) slider × 10. The undisturbed free-stream maps to zero (transparent), regions faster than free-stream (e.g. flow accelerating around the sides) appear red, and the stagnation / wake (slower) regions appear blue. A bilateral threshold on \(|v_{display}|\) means both acceleration and deceleration contribute to opacity. To suppress numerical checkerboard noise (a typical artifact of collocated-grid pressure solvers) the velocity is sampled at the corner between 8 neighbouring cells so the GPU's trilinear filter averages them to zero for any ±δ alternating pattern — at no extra texture read cost.

8. Lagrangian Particles & Trails

Both 2D and 3D modes support massless Lagrangian particles advected by the velocity field.

Sub-Step Tunneling Prevention (3D) 3D

At high flow speeds a particle's displacement in one frame can exceed the obstacle diameter, causing it to pass through solid geometry without detection (tunneling). The 3D advection kernel solves this by splitting each full step into 4 sub-steps of \(\tfrac{1}{4} \Delta t_{vis}\) each. The SDF is tested after every sub-step; the particle is guaranteed to be caught within one cell-diameter of travel.

Obstacle Surface Interaction (3D) 3D

When a particle's sub-step would carry it inside an obstacle (SDF < 0.01), two strategies are applied in sequence rather than teleporting the particle back to the inlet:

  1. Tangential sliding: The sub-step velocity is projected onto the surface tangent plane using the SDF normal \(\hat{\mathbf{n}} = \nabla(\text{SDF})\) at the last-valid position:
    $$ \mathbf{v}_{tang} = \mathbf{v} - (\mathbf{v} \cdot \hat{\mathbf{n}})\,\hat{\mathbf{n}} $$
    The particle slides along the surface with this tangential component. If the slid position is still inside the SDF the particle is frozen at the last valid position for that sub-step.
  2. Persistence: Particles that become stuck on an obstacle surface remain visible at that position rather than disappearing. They resume motion once the flow carries them off the surface.

Field-Colored Particles 2D3D

Toggle Color by Field (in the Streamlines/Particles or 3D Particles folder) to color each particle and trail segment by the local flow field value at its current position, using the active colorbar theme.

  • Vorticity view: particles show local rotation intensity — strong vortices appear at the colormap extremes, irrotational flow is near-neutral.
  • Velocity view: particles range from blue (below free-stream speed) through neutral to red (above free-stream).
  • Pressure view: high-pressure stagnation regions vs. low-pressure suction regions are directly visible on each particle.

When disabled, particles render as white (dark theme) or black (light theme). The field value is sampled in the GPU vertex shader at zero extra texture-fetch cost, so performance impact is negligible.

Temperature (Buoyancy) Model 3D

An optional Temperature (Buoyancy) folder (above Performance in the 3D sidebar) adds Boussinesq natural-convection physics to the 3D incompressible solver — designed for greenhouse ventilation, stack-effect studies, and any scenario where a temperature difference between zones drives airflow.

Physics: the Boussinesq approximation

The solver remains fully incompressible. Density is assumed uniform everywhere except inside one vertical body-force term. The acceleration added each timestep is:

$$\mathbf{f}_{buoy} = g\beta\,T(\mathbf{x})\,\hat{y}$$

where \(T\) is the local temperature deviation from ambient (ambient = 0), \(\beta \approx 1/T_{abs} \approx 0.0034\,\mathrm{K}^{-1}\) for air, and \(g\beta\) is the single Buoyancy slider. This force is applied before the pressure-projection step, so the pressure Poisson solve produces the inside/outside pressure difference that drives flow through windows and vents — the physics is not prescribed but emerges from the equations, making the model self-consistent with any wind-driven inlet flow.

Temperature field and zone sources

A 3D scalar field \(T(\mathbf{x},t)\) is advected and diffused by the velocity field each substep (semi-Lagrangian advection + explicit Laplacian diffusion). Inlet planes always supply ambient air (\(T=0\)). A user-defined box zone injects heat in one of two modes:

  • Fixed Temperature (Dirichlet): the zone is held at the set overheat each step — analogous to a surface at constant temperature (warm soil, heated floor). Use Zone Temp (dT) to set the value.
  • Heat Flux (volumetric source): the zone temperature increases at a constant rate each step — analogous to solar gain or an electric heater continuously adding energy. Use Flux Rate. This mode has no direction: it is a W/m³-style volumetric source, not a directional beam. Direction appears only in the subsequent advection (carried by the wind), diffusion (spread isotropically), and buoyancy (warm air rises in +Y).

Negative values make a cooling zone (chilled wall, evaporative cooling) — the resulting cold air sinks.

Zone Temp units

Zone Temp is a normalized overheat above ambient, not a Celsius value. The user defines the mapping, for example "Zone Temp 1.0 = 20 °C overheat". The buoyancy force at each cell is Buoyancy × T_local (pointwise), so the slider acts on whatever temperature gradient the field has developed.

Buoyancy slider and the Richardson number

The Buoyancy slider encodes \(g\beta\) in simulation units. Its relationship to the physical Richardson number \(Ri = g\beta\Delta T\,H / U^2\) (the ratio of stack-effect strength to wind strength) is approximately:

$$\text{Buoyancy} \approx \frac{7 \cdot Ri}{T_{zone}}$$

Practical guidance for a greenhouse (height \(\approx 3\,\text{m}\), 10 °C overheat, Zone Temp = 0.6):

  • Breezy day (\(\geq\)3 m/s, \(Ri \approx 0.1\)) → Buoyancy 1–2 (wind dominates)
  • Light wind (\(\approx\)1 m/s, \(Ri \approx 1\)) → Buoyancy 8–15 (default 10; wind and stack comparable)
  • Near calm (\(\approx\)0.3–0.5 m/s, \(Ri \approx 4\text{–}10\)) → Buoyancy 40–50 (stack dominates)

Persistent field and the "Clear Heat" button

The temperature field is a GPU volume that persists across frames. Setting Zone Temp or Flux Rate to 0 stops injecting new heat but does not remove heat already in the domain — residual warmth floating in recirculation zones continues to drive buoyancy. This is physically correct behaviour (real greenhouse air does not instantly return to ambient when a heater is turned off). Press Clear Heat (back to ambient) to instantly zero the entire field. The field also resets automatically on sim re-init.

Unsteady behaviour is normal

Natural convection with residual heat is inherently unsteady — thermal plumes flicker and meander (like smoke above a candle). Only the time-averaged field is steady. Use the Steady Average toggle in the Visualization folder to observe the converged mean plume, exactly as a CFD engineer would report time-averaged natural-convection results. Setting Buoyancy = 0, or pressing Clear Heat while sources are off, produces a truly neutral (isothermal, no-buoyancy) flow.

Temperature view mode

Select Temperature in the Field dropdown (Visualization folder) to render the thermal plume directly. Ambient air (T = 0) is transparent; warmer regions render red, cooler regions render blue. The ghost box (Show Zone) draws the source zone as a translucent orange (heating) or blue (cooling) box — occluded by the STL obstacle — to help with positioning.

Ventilation Fan Model 3D

An optional Ventilation Fan folder (directly below Temperature (Buoyancy) in the 3D sidebar) lets the user place up to 4 independent fan zones — for cross-ventilation studies, forced-convection cooling, or modelling a spinning rotor's downwash without resolving actual blade geometry. Each fan has a Shape: Box (a straight-through vent) or Disk (an actuator-disk rotor/propeller).

Physics: two different actuator models

Box is a hard velocity clamp — the standard simplified-CFD treatment for a duct/desk/extractor fan: inside the zone, velocity is hard-set every step to a fixed directional jet,

$$\mathbf{u}(\mathbf{x}) = P \, \hat{d} \quad \text{for } \mathbf{x} \in \text{zone}$$

where \(P\) is the fan's Power (Jet Speed) and \(\hat{d}\) is its unit direction vector — a Dirichlet velocity boundary condition, the same technique the domain inlet itself uses. This behaves like a section of duct: the flow inside simply is that speed, which suits a straight-through vent but cannot, by itself, distinguish an upstream (suction) side from a downstream (discharge) side. Incompressible engine only.

Disk instead adds momentum to whatever air passes through it, every step:

$$\mathbf{u}(\mathbf{x}) \mathrel{+}= P \, \hat{d} \, \Delta t \quad \text{for } \mathbf{x} \in \text{disk}$$

This single change — an increment rather than a replacement — is what makes a Disk fan behave like an actual rotor: because the zone no longer forces a fixed value, the pressure/continuity solve is free to draw replacement air in from behind the disk (suction) and build pressure up just ahead of it (discharge), a genuine upstream/downstream pressure difference — which is aerodynamic thrust — instead of a section of duct. Disk works in both engine modes: the compressible solver applies the equivalent source term directly to its conservative momentum and energy variables at the end of each Runge-Kutta stage (adding the matching kinetic-energy delta alongside momentum, so the derived pressure doesn't silently jump). Box has no clean equivalent as a source term in a conservative finite-volume scheme, so it remains incompressible-only.

Where two Box zones overlap, the higher-numbered fan's velocity wins (a simple sequential override). Where two Disk zones overlap, their momentum contributions sum — the physically correct combination for multiple momentum sources acting on the same air.

The zone itself is oriented — not just the velocity inside it

Both shapes' zones are built in a local frame derived directly from the direction vector \(\hat{d}\) (one axis along \(\hat{d}\), the other two spanning the plane perpendicular to it) — so tilting Direction Yaw/Pitch rotates the zone's own shape bodily, not just the velocity direction inside a fixed footprint. For a Box this means its own "outlet face" always faces wherever it's aimed, exactly like an oriented duct segment; for a Disk it means the disk's flat face is always perpendicular to its thrust axis, exactly like a real rotor plane. At the default direction (\(0°,0°\), blowing toward \(+X\)) this local frame reduces exactly to world \(X\)/\(Y\)/\(Z\), so an un-rotated fan looks identical to the original world-axis-aligned box.

Zone Size X/Y/Z follow that same local frame rather than world axes: Size X is always the half-extent along the thrust direction (thickness), and Size Y/Size Z are the two extents across it. For Disk, Size Y is the radius and Size Z is unused. Zone Center X/Y/Z and all sizes are normalized (\(-1\) to \(+1\) for center, \(0.02\) to \(1\) for size) and scaled by the domain half-extents, the same convention as the Temperature zone.

Multiple fans

Fans are added and removed independently via + Add Fan / Remove buttons, each opening its own sub-folder ("Fan 1", "Fan 2", ...) with its own shape, position, size, direction, and power — the same add/remove pattern as the Extra Obstacles (3D) folder. The Enable Fans toggle is a single master switch for all fans at once: turning it off zeroes every fan's effect on the next step without deleting any of their configurations, so settings are preserved if re-enabled.

Visualizing direction: the ghost zone and arrow

Show Zone + Arrow draws, for every active fan, a translucent teal ghost (a different colour from the Temperature zone's orange/blue, so the two are never confused) marking its true oriented extent — an oriented box for Box fans, an actual tilted cylinder (not a bounding box) for Disk fans, both computed with a proper ray intersection in the zone's own local frame so the ghost always matches what the solver is really doing — plus a shaft-and-arrowhead glyph spanning the zone and pointing exactly along the direction it blows. The glyph reuses the same rendering technique as the Aerodynamics tab's Force Vectors arrows.

No persistent field

Unlike the Temperature model, a fan has no memory: it is re-applied fresh every step, with nothing to advect or diffuse. There is therefore no "Clear Fan" button — toggling Enable Fans off removes the effect immediately on the next frame.

Modelling a drone/propeller blade's downwash

A Disk fan is the recommended way to model a spinning rotor's induced flow: real rotor RPM is far higher than this explicit-timestep solver could resolve for actual rotating blade geometry without an impractically small CFL-limited timestep and a finer grid than the 3D solver's 256³ cap — an actuator disk sidesteps this entirely by imposing the induced velocity directly, the same simplification real rotor-aerodynamics tools use. The alternative — a literally spinning, blade-pitched STL mesh (see Obstacle Spin below) — gives real blade-tip vortices and a measurable Torque, but is only realistic at low RPM.

9. Integrated AI Agent Architecture LLM

OTOM CFD includes a built-in CFDChatbot that allows users to alter the simulation using natural language. The architecture works as follows:

Function Calling & Parameter Mapping

The agent communicates with the Gemini API using strict JSON schemas (Function Calling). When a user asks "Make the fluid thicker" or "Add a large spinning wing," the LLM maps intent to two primary functions:

  • update_settings: Maps semantic concepts to internal variables (e.g., "thicker" \(\rightarrow\) viscosity: 0.05, "slower" \(\rightarrow\) inletVel: 5.0).
  • set_obstacles / modify_existing_obstacles: Translates geometry requests into SDF parameters (e.g., obsType: 'Engineering', engShape: 'NACA 0012', angle: 15, scale: 2.0).

The Chatbot class intercepts these tool calls and updates the Lil-GUI controllers, which in turn flush the updated data to the WebGPU uniform and storage buffers.

Array-valued 3D features without a dedicated action — the Ventilation Fan's fans3D list and the Extra Obstacles (3D) list — are set through the same update_settings call, passing the whole array as one of the "changes". Because update_settings replaces whichever keys it touches rather than merging into them, asking the agent to add one more fan requires it to resend every existing fan alongside the new one; it has no live read-back of the current scene and relies on the conversation history to know what is already there.

10. Aerodynamics & Analytics

The software calculates Drag (\(C_d\)) and Lift (\(C_l\)) coefficients by integrating pressure and viscous shear forces over the surface of the obstacles.

The Reynolds number is calculated dynamically based on the largest obstacle's chord length \(L\):

$$Re = \frac{U_\infty L}{\nu}$$

Analytics can be exported as structured text or baked into PDF reports with canvas screenshots using the GUI tools.

Torque — measuring a spinning blade's efficiency 3D

When the 3D obstacle is spinning (see Obstacle Spin above), the same surface integral that produces Drag/Lift also accumulates Torque about whichever spin axis is currently active — the shaft torque a motor would need to supply to hold that rotation rate:

$$\tau = \oint \big(\mathbf{r} \times (-p\,\mathbf{n})\big) \cdot \hat{\boldsymbol{\omega}} \; dA$$

where \(\mathbf{r}\) is measured from the obstacle's own centre and \(\hat{\boldsymbol{\omega}}\) is its current spin axis — exactly the same axis/centre solid_vel3d uses for the wall-velocity boundary condition, so the reading is about the same rotation actually driving the flow. Torque reads exactly 0 whenever the obstacle isn't spinning (there is no single shaft axis to report a torque about). Combined with Drag/Lift, this is what lets Cd/Cl (thrust) and Torque (cost) together describe a spinning blade's aerodynamic efficiency — shaft power is Torque × angular speed, though the app displays the raw torque rather than computing that product.

Surface Friction — the pressure/friction drag split 3D

The 3D surface integral now resolves the viscous component of the wall force alongside the pressure component. Wall shear stress is evaluated as

$$\boldsymbol{\tau} = \mu \frac{\partial \mathbf{u}_t}{\partial n} \approx \mu \,\big(\mathbf{u}_{\text{fluid}} - \mathbf{u}_{\text{wall}}\big)_t$$

sampling the fluid velocity one cell along the outward normal and differencing it against the wall's own velocity — the separation is exactly one cell and velocity is stored in cells/step, so the finite-difference denominator is unity. Using \(\mathbf{u}_{\text{wall}}\) rather than zero is what keeps this correct for a spinning obstacle: the shear a rotating blade feels is set by the fluid's velocity relative to the moving surface, which for a well-driven rotor is far smaller than the absolute velocity and would otherwise be grossly overestimated.

The Surface Friction sub-folder reports Cd — Pressure, Cd — Friction (which sum to the headline total \(C_d\)) and the area-weighted Mean Wall Shear. The ratio is a useful shape diagnostic: pressure drag dominates for bluff bodies, friction drag for slender bodies aligned with the flow. Friction also contributes to the reported Torque, which matters for rotors — omitting it would systematically under-report required shaft power.

Propulsion in real units — Thrust, Power, RPM 3D · SI

The Propulsion (SI) sub-folder converts the lattice-unit force and torque integrals into Newtons and Watts. The conversion is derived from the same relation the solver already uses for velocity in SI Units Mode (\(u_{\text{lat}} = u_{SI}\,\Delta t / \Delta x\)), so one lattice velocity unit is \(\Delta x/\Delta t\) m/s. Because the projection step computes \(\mathbf{u} \mathrel{-}= \nabla p\) with velocity in cells/step, the solver's pressure is a kinematic pressure \(p/\rho\) carrying units of velocity squared — so converting it to Pascals requires \(\rho\,(\Delta x/\Delta t)^2\) with no extra length factor:

$$F_{SI} = \rho \left(\frac{\Delta x}{\Delta t}\right)^{2} m_{\text{world}}^{2}\; F_{\text{lat}}, \qquad \tau_{SI} = F_{SI}\, m_{\text{world}}$$

where \(m_{\text{world}}\) is metres per world unit. RPM follows from the spin sliders (degrees per frame, converted through Physics Steps and Real-world Time Step), thrust is the force projected onto the spin axis when rotating (streamwise \(+X\) otherwise), and shaft power is \(P = |\tau|\,\omega\). Thrust / Power (N/W) is reported rather than a "figure of merit" because it needs no disk-area assumption and is the metric rotor designers actually compare.

These fields require SI Units Mode (and Live Analysis). Lattice numbers carry no physical scale whatsoever, so rather than showing a fabricated Newton value the fields state "needs SI Units Mode". The physical scale is set entirely by Domain Width (m) and Real-world Time Step, so those must describe your actual case for the output to mean anything.

Wind Noise — aeroacoustic proxy 3D

The Wind Noise folder (directly below Ventilation Fan) estimates flow-generated noise from the rate of change of surface pressure on the obstacle. This is the source term of Curle's acoustic analogy: at low Mach number the dominant noise mechanism is not the turbulence itself (a weak quadrupole, scaling as \(M^8\)) but the fluctuating pressure on solid surfaces, radiating as a dipole that scales as \(M^6\) — which is why a body in a flow is far louder than the same flow without it.

$$p'(r) \;\sim\; \frac{1}{4\pi c\, r}\sqrt{\;\int_S \left(\frac{\partial p}{\partial t}\right)^{2} A_c \; dS }$$

The integrand is accumulated once per solver substep — not on the ~2 Hz analysis timer, which is far too slow to resolve the unsteadiness that actually radiates. Noise Index (dB rel.) is a relative scale with an arbitrary zero point, always available and the right number for A/B comparisons. Est. SPL is an absolute estimate re \(20\,\mu\text{Pa}\) at the chosen Observer Dist, and requires SI Units Mode; it takes the surface correlation area \(A_c\) as one cell face, the standard resolved-scale assumption.

What this is not. There is no acoustic wave equation, no propagation delay, no directivity, no reflection or diffraction, and turbulence finer than the grid radiates nothing at all. The number is a source-strength indicator for ranking designs — genuinely useful for answering "is this shape quieter than that one?" — and must not be presented as a certifiable or measured sound-pressure level. Enabling the panel also costs one extra full-grid pass per substep plus a \(\text{res}_{3D}^3\) history texture pair, so it is off by default.

Ground Truth — Reference Comparison Truth tab

The Ground Truth tab (fact_check icon in the left sidebar) cross-checks the live simulation against textbook and empirical correlations for the active canonical shape. It auto-detects the current geometry — 2D Circle/Square/Rectangle/Triangle/Hexagon/Pentagon and Engineering shapes (NACA airfoils, Flat Plate, Ellipse, Tear Drop, Wedge), or 3D Sphere/Cube/Cylinder/Wing — and presents a Measured | Reference | Δ% table (colour-coded by error), a reference plot with the current operating point marked, and an engineering assessment.

Reference correlations include the circular cylinder fit \(C_d \approx 1 + 10\,Re^{-2/3}\) (White) with Strouhal \(St \approx 0.20\); the sphere Schiller–Naumann/Morrison drag law; the laminar flat plate \(C_d = 1.328/\sqrt{Re}\) (Blasius, parallel) or \(C_d \approx 1.98\) (normal); a generic sharp-edged bluff-body fit for prisms (e.g. square \(\approx 2.05\) face-on, \(\approx 1.55\) corner-on; triangle \(\approx 1.55\) apex-upstream, \(\approx 2.0\) base-upstream); and thin-airfoil theory \(C_l = 2\pi\sin\alpha\) closed by a parabolic drag polar \(C_d = C_{d0} + k\,C_l^2\). The reference plot is \(C_d\) vs \(Re\) (log axis) for bluff bodies, or the \(C_l\) vs \(\alpha\) lift line for airfoils.

Honesty of the measured column: in 3D the “Measured” values are a genuine GPU surface-force integral. In 2D they are presently an illustrative placeholder (marked with an asterisk) pending a dedicated 2D force kernel; only the Reference column is authoritative in 2D. Hand-drawn or imported (STL) geometry has no closed-form solution, and the panel states this rather than fabricating a number.

Engineering Judgement

Below the table, a scene-aware assessment translates the live observables into plain-language guidance with colour-coded markers (✓ good, ! caution, • info): the solver-vs-theory agreement, the Reynolds flow regime and expected wake behaviour (creeping → steady recirculation → laminar von Kármán street → subcritical → drag crisis), vortex-shedding side-force and vortex-induced-vibration warnings, airfoil stall (incidence beyond \(\sim\!12\text{--}15^\circ\)) and lift-to-drag ratio, and the compressibility regime (incompressible/subsonic/transonic/supersonic, including the Mach angle \(\mu = \arcsin(1/M)\)).

Local AI Scene Explanation & Chat WebLLM · in-browser

The Ground Truth tab includes an optional “Explain this scene” assistant that runs a small large-language model entirely in the browser on the user’s GPU via WebLLM (MLC). It produces a plain-language explanation of the current scene and then supports a multi-turn conversation for follow-up questions. The model is grounded on the structured scene data — the \(C_d\), \(C_l\), \(Re\), Mach, regime and judgement values come from the solver and reference correlations, not the model — so it narrates and reasons but does not invent coefficients.

A dropdown selects the local model (Llama 3.2 1B, Gemma 2 2B default, Qwen2.5 3B, Llama 3.2 3B, or DeepSeek-R1 Distill 7B). The chosen model is downloaded from the MLC CDN on first use and cached in-browser thereafter; everything runs locally, so no scene data leaves the device. On hardware without WebGPU (or if a model cannot be allocated), it transparently falls back to the cloud assistant. The simulation auto-pauses during generation so the language model and the renderer do not contend for the GPU.

11. AI Surrogate Models Experimental · ML

The AI tab (brain icon in the left sidebar, directly below Demo) provides a machine-learning workflow for training and running neural-network surrogate models that approximate the physics solver. A surrogate learns the solver's behaviour from data and, once trained, can advance the flow with a single network evaluation instead of the full iterative solve. This is offered as an experimental research feature: surrogates trade physical fidelity for speed, and are most interesting in 3D where the iterative pressure solve is the dominant cost.

Concept: Learning the Per-Step Residual

Rather than predicting the next flow field directly, the networks are trained to predict the residual — the change in the field over one solver step — which the simulation then integrates:

$$ \mathbf{q}_{t+\Delta t} = \mathbf{q}_{t} + \mathcal{N}_\theta(\mathbf{q}_{t}, \, \text{mask}, \, Re) $$

where \( \mathbf{q} \) is the stacked velocity/pressure state, \( \text{mask} \) is the obstacle occupancy, \( Re \) is the (normalised) Reynolds number, and \( \mathcal{N}_\theta \) is the trained network. Predicting the residual keeps the targets small and well-conditioned and improves step-to-step stability. All inputs and outputs are normalised internally by fixed reference scales so the network operates on \(\mathcal{O}(1)\) values.

Resolution lock: a network's weights encode the physical scale of the grid it trained on, so each model is bound to a single grid resolution. Applying a model automatically sets the simulation to that grid. To run at a different resolution you train a separate model — resolution is part of a model's identity, not a runtime knob.

2D Training — Train New Model

Trains a compact U-Net (three encoder/decoder levels with skip connections, ~tens of thousands of parameters) on the 2D incompressible solver using TensorFlow.js. The data pipeline drives the real GPU solver headlessly across randomised parameters, captures (state, next-state) pairs near quasi-steady state, and converts them to normalised residual targets. Controls:

  • Model name — identifier used for saving / loading.
  • Obstacle — Circle (the 2D v1 shape).
  • Re min / Re max, U_in min / U_in max — the Reynolds-number and inlet-velocity ranges sampled during data collection.
  • Samples — number of randomised simulation runs; each contributes several training pairs.
  • Epochs — passes over the collected dataset (Adam optimiser, mean-squared-error loss).
  • Status — live progress: Collecting…, then Epoch n — loss …, then Saved ….

2D inference uses a CPU round-trip (read the field off the GPU, evaluate the network in TensorFlow.js, write the result back). At small 2D grids the data volume is tiny, so this stays real-time.

3D Training — Train 3D Model

Trains a small fully-convolutional 3D CNN (a flat stack of 3×3×3 convolutions, 8 → 16 → 16 channels, plus a 1×1×1 output convolution — deliberately not a U-Net so the GPU inference path needs only one kernel type). The data pipeline runs the real 3D incompressible solver headlessly, reads back the velocity and pressure volumes, and bakes the obstacle's occupancy mask as an additional input channel. Controls:

  • Model name.
  • Grid (³) — cubic resolution: 32, 64, 128, 256, 512 (the same set the 3D solver offers, so an applied model always lands on a valid sim resolution).
  • Obstacle — a dropdown: choose a single shape (Sphere, Box, Cylinder, Wing) or All to train one model that generalises across every shape. Multi-shape generalisation is possible because the obstacle's signed-distance / occupancy mask is supplied to the network as a dedicated input channel.
  • Re / U_in ranges, Samples, Epochs, Status — as in the 2D trainer (3D uses tiny batch sizes because the volumes are large).

Input channels (6): \(v_x, v_y, v_z\), pressure, obstacle occupancy (SDF), normalised Reynolds number. Output channels (4): residual \(\Delta v_x, \Delta v_y, \Delta v_z, \Delta p\).

⚠ Memory & cost at large grids. Cost scales with the number of cells, which grows as the cube of the resolution. Grids 16–64 are practical; 128 is heavy; 256 and 512 are extreme and may exhaust GPU or CPU memory during data collection, dataset assembly, or training (a single 256³ activation tensor is on the order of 1 GB; 512³ is ~8× larger). These options are provided for completeness and experimentation — for a first run prefer Grid 16–32, one shape, ~8–12 Samples, ~10 Epochs, and increase only once it trains successfully.

3D Inference — GPU-Resident Forward Pass

Unlike the 2D path, the trained 3D network runs entirely on the GPU. Its weights are uploaded to GPU buffers and the forward pass executes as a chain of WebGPU compute shaders operating directly on the simulation's velocity/pressure volumes:

$$ \text{pack} \;\rightarrow\; \text{conv3d} \times 4 \;\rightarrow\; \text{write-back} $$

The pack stage gathers the velocity, pressure, and obstacle-mask channels into a packed input buffer; each conv3d stage performs one zero-padded ("same") convolution with ReLU; the write-back stage integrates the predicted residual and writes the new velocity and pressure back to the volume textures (zeroing velocity inside solid cells). Because there is no CPU round-trip, the surrogate can in principle outrun the iterative solver — whose 3D cost is dominated by the many Jacobi pressure-Poisson iterations per step. The advantage is largest at higher resolutions; at small grids the physics solver is already inexpensive, so the surrogate may show little or no speed-up.

Correctness self-test. When a 3D model is applied, a one-off check runs the GPU forward pass and compares it against the reference TensorFlow.js evaluation on the identical input, logging the maximum absolute difference to the browser console. A small value (≈\(10^{-3}\) or below) confirms the GPU kernels reproduce the trained network.

Managing & Applying Models — My Models / Apply Model

My Models lists every saved surrogate with its parameter ranges and final training loss; each row has a trash icon to delete it. Models (weights plus a metadata record) persist in the browser through IndexedDB, so they survive page reloads.

Apply Model selects an Active Model and toggles Enable AI Solver. Enabling a model swaps the surrogate in for the physics solver: a 2D model switches the simulation to 2D incompressible at the model's grid; a 3D model switches to 3D incompressible at the model's grid (open the 3D tab to view the result). Disabling the toggle — or selecting None — instantly restores the real physics solver.

Scope, Accuracy & Practical Tips

  • Experimental: surrogates approximate the solver and are not a substitute for it when accuracy matters. A small network trained on few samples can drift or become unstable over many steps; if the field degrades, disable the AI Solver to recover.
  • Improving quality: increase Samples (raises the achievable accuracy ceiling), increase Epochs (reaches that ceiling), and keep the Re/U_in ranges as narrow as your use-case allows (less variety means the network has less to generalise over). Watch the loss in Status: if it is still falling at the last epoch, add epochs; if it has plateaued, add samples or accept the limit.
  • Speed comparison: to judge whether the surrogate is actually faster, watch the FPS readout in the Performance folder with the AI Solver on versus off at the same grid.
  • Supported obstacles (3D training): Sphere, Box, Cylinder, Wing (the analytical primitives whose occupancy can be computed directly). STL meshes and extruded 2D profiles are not yet covered by the 3D trainer.

13. Fly-Through Camera Mode (3D) Navigation

Inspired by the CATIA V5 Fly Mode, the 3D view supports a first-person fly-through camera that lets you navigate inside the control volume — go forward, backward, turn left/right, look up/down — independently of the normal orbit/zoom interaction.

Activating & Deactivating

  • Activate: Click the middle mouse button (scroll wheel) while the 3D view is active. A HUD overlay appears with a crosshair and a blue directional arrow pointing in the flight direction.
  • Deactivate: Click the middle mouse button again, or press Esc.
  • On exit the orbit camera is automatically re-centred at the position where you stopped, so normal orbit/zoom continues smoothly from there.

Controls

InputAction
Mouse movementSteer — left/right yaws the view; up/down pitches. Moving the mouse to the left turns the camera left (intuitive first-person mapping).
Scroll wheel upDecrease forward speed.
Scroll wheel downIncrease forward speed.
Right mouse button (hold)Fly backward while held; release to resume forward flight.
Middle-click or EscExit fly mode.

HUD Overlay

While fly mode is active, a semi-transparent overlay is drawn on top of the 3D canvas:

  • A ring crosshair at the screen centre with four orthogonal tick marks.
  • A blue chevron arrow above the crosshair indicating the flight direction (always pointing forward into the scene).
  • A status bar at the bottom of the screen showing the current speed (e.g. FLY MODE | Speed: 0.20) and a reminder of the key controls.

Default Speed & Range

The default forward speed is 0.2 world units / second. The control volume spans ±1 unit per axis by default, so the camera crosses it in about 10 seconds at default speed. Speed is adjustable in real time via the scroll wheel over the range 0.05 – 10 world units/s. If the Pointer Lock API is available in the browser, the cursor is hidden during fly mode and movement is captured as raw deltas for smoother steering.

14. Session Persistence UX

All simulation settings and the current 3D camera state are automatically saved to localStorage under the key wl_session_v1. The session is restored when the page is reloaded — you continue from exactly where you left off.

What is Saved

  • All GUI-exposed parameters (inlet velocity, viscosity, obstacle type, domain size, rendering options, etc.).
  • The 3D camera state (azimuth θ, elevation φ, radius r).
  • Whether the 2D or 3D tab was active when the session ended.

Ephemeral display values (FPS, frame counter, live Cd/Cl/Reynolds, STL triangle count) are intentionally excluded so that restored values are always meaningful.

Save Triggers

  • On page unload (beforeunload event) — fires when you close the tab, navigate away, or refresh.
  • Every 30 seconds — a periodic autosave guards against browser crashes.

To start fresh, clear localStorage in your browser's developer tools (Application → Storage → Local Storage) and reload.

Backward-Compatible Loading

The GUI grows new per-item fields over time (e.g. the Spin X/Y/Z and Show Local Axes controls added to Extra Obstacles and Ventilation Fans). A session, saved scene file, or Demo written before such a field existed won't have it — restoring that old data verbatim would leave the corresponding slider with no value to bind to and fail to build the sidebar. Every load path (session restore, Load button, Demo, and the AI assistant's scene edits) therefore backfills any missing per-item field with its default before rebuilding the GUI, so older saves keep working after an update instead of failing to load.

Play/Pause and Reset All Settings

The bottom of both the 2D and 3D sidebars ends with two controls, below Load/Save:

  • Play/Pause — a single button (not a checkbox) that starts or stops the simulation loop. Since only one simulation runs at a time (whichever of the 2D/3D tabs is active), the 2D and 3D sidebars each have their own button, but both control the same underlying state and stay in sync with each other.
  • Reset All Settings — restores every setting, in both the 2D and 3D panels, back to the software's original defaults: obstacles, Extra Obstacles, Ventilation Fans, the Temperature/Buoyancy zone, colormap themes, engine mode, everything. This is distinct from RESET FIELDS / Reset 3D Simulation, which only re-seed the flow field and leave your configuration untouched. Because it discards the current setup with no undo, it asks for confirmation before proceeding.

12. Academic References & Theoretical Foundations

The numerical methods and boundary immersion techniques driving this WebGPU solver heavily draw upon academic research published in top-tier Elsevier journals. For further reading on the mathematics of WaterLily-style models, consult the following literature:

Weymouth, G. D., & Yue, D. K. P. (2011).
Boundary data immersion method for Cartesian-grid simulations of fluid-body interaction problems.
Journal of Computational Physics, 230(22), 6233-6247. [Elsevier]
Establishes the BDIM framework used for obstacle rendering and interaction without conformal meshing.
Maertens, A. P., & Weymouth, G. D. (2015).
Accurate Cartesian-grid simulations of near-body flows at intermediate Reynolds numbers.
Computer Methods in Applied Mechanics and Engineering, 283, 106-129. [Elsevier]
Expands on projection methods and pressure Poisson solvers over non-conformal grids, foundational for the Jacobi iterations used herein.
Laizet, S., & Lamballais, E. (2009).
High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy.
Journal of Computational Physics, 228(15), 5989-6015. [Elsevier]
Provides context on high-order spatial discretizations necessary for maintaining stability in Cartesian solvers at high Reynolds numbers.
Toro, E. F., Spruce, M., & Speares, W. (1994).
Restoration of the contact surface in the HLL-Riemann solver.
Shock Waves, 4(1), 25-34. [Springer]
Introduces the HLLC approximate Riemann solver used by the compressible Euler engine to capture shocks, contact surfaces, and shear waves.
Antoniadis, A. F., Drikakis, D., Farmakis, P. S., et al. (2022).
UCNS3D: An open-source high-order finite-volume unstructured CFD solver.
Computer Physics Communications, 279, 108453. [Elsevier]
High-order finite-volume methodology (MUSCL/WENO reconstruction with Riemann solvers) that informs the compressible engine's flux-reconstruction roadmap.
Smagorinsky, J. (1963).
General circulation experiments with the primitive equations.
Monthly Weather Review (Historical context; modern implementations reviewed widely in Elsevier's Computers & Fluids).
The basis for the Large Eddy Simulation (LES) sub-grid scale modeling implemented in the software's turbulence settings.