# waterfall

Fast incompressible CFD for external aerodynamics. Solves the Navier-Stokes
equations on a uniform staggered grid with an immersed boundary, sized so a
full vehicle case finishes inside 30 minutes on a desktop CPU.

The design target is **fast iteration on low-poly geometry**: drop in an STL,
get a drag number, a pressure field and a wake video, change the shape, run
again. It is not a substitute for a wall-resolved commercial solve, and the
places where it is wrong are measured and written down in
[VALIDATION_REPORT.md](VALIDATION_REPORT.md) rather than left for you to find.

## What it actually solves

```
du/dt + div(u u) = -grad(p)/rho + nu*lap(u) + f_ibm
div(u) = 0
```

Incompressible, isothermal, Newtonian, constant density.

| Piece | Choice | Why |
|---|---|---|
| Grid | Uniform staggered MAC | Lets the pressure solve be a direct transform |
| Pressure | DCT-II Poisson, all-Neumann | Direct, non-iterative. 57 ms at 2.4M cells |
| Time | Adams-Bashforth 2, one projection/step | 2nd order for one Poisson solve, not three |
| Convection | 2nd-order central, divergence form | Non-dissipative; upwinding would fake viscosity |
| Diffusion | Explicit 7-point Laplacian | Convective CFL binds first, so it is free |
| Geometry | Immersed boundary, fractional `eta` | No meshing step; STL change costs seconds |
| Turbulence | Vreman LES (default) | Required above Re ~1e5; benchmarks set `sgs: null` and run laminar |

The FFT pressure solver is the load-bearing decision. It is what makes the
runtime budget achievable, and it is why the grid must be uniform: there is
no local refinement near the body, and that is the method's main cost.

## Install

```bash
python -m venv .venv
source .venv/bin/activate        # Windows: .venv\Scripts\activate
```

```bash
python -m pip install numpy scipy trimesh matplotlib rtree numba scikit-image imageio imageio-ffmpeg pyyaml fastapi "uvicorn[standard]" python-multipart
```

Numba is required, not optional. The kernels are the solver.

## The app

Double-click **`waterfall-ui.bat`**, or:

```bash
python ui.py
```

A browser opens at `http://127.0.0.1:8010` (loopback only, so nothing is exposed
to your network).

- **New run**: drag an STL onto the page. It reads the mesh, guesses the flow
  axis, and draws side/plan/front previews with the flow direction marked.
  Type the part's real length in metres and the scale factor is worked out for
  you. The runtime and cell count update live as you change settings.
- **Running**: progress bar, ETA, step count, live divergence and force, a
  tail of the solver log, and a Stop button.
- **Projects**: every past run as a card with its Cd and a thumbnail. Open one
  for the result tiles, the animation, the geometry preview, the full log, and
  every file it produced, downloadable.

**Always look at the orientation preview before running.** The flow-axis guess
assumes the longest dimension faces the flow, which is right for a car and
wrong for a bumper, wing or splitter. A wrong axis does not raise an error; it
returns a confident, wrong drag number.

Results in the app carry their own caveats: if subgrid viscosity or resolution
made a number untrustworthy, the run says so next to the number.

## Run a case from the terminal

The Tesla case geometry (`geometries/tesla_model_s_2023.stl`) is a third-party
model and is not shipped; bring your own STL and point `geometry.path` at it.

```bash
python waterfall.py cases/tesla.yaml --check
```

`--check` runs the pre-flight only: geometry validity, orientation preview,
domain bounds, blockage, cell count, and an estimated wall-clock time. Always
run it first. The single most expensive mistake is an `axis_order` that puts
the model sideways to the flow, because that still produces a plausible drag
number.

```bash
python waterfall.py cases/tesla.yaml
```

Writes `summary.json` (Cd, Cl, forces, diagnostics) and a directory of frames.

```bash
python visualize_waterfall.py output/tesla/frames --layout dashboard --fps 30
```

Writes PNG frames and an MP4: coloured slice, streamlines, and a 3D
Q-criterion iso-surface of the wake.

## Configuration

See [cases/tesla.yaml](cases/tesla.yaml) for a fully commented example.

```yaml
geometry:
  path: geometries/car.stl
  units: mm
  axis_order: yzx        # which STL axis becomes solver x (streamwise), y (up), z (span)
domain:
  cells_per_length: 40   # resolution, measured along the STREAMWISE extent
  upstream: 1.5          # domain padding, in body lengths
  downstream: 4.0
flow:
  speed: 30.0            # m/s
  viscosity: 1.5e-5      # kinematic, m^2/s
  rolling_road: true     # a fixed floor is a tunnel artefact
run:
  convective_times: 6.0
  output: output/car
```

Override from the command line without editing the file:

```bash
python waterfall.py cases/tesla.yaml --override domain.cells_per_length=24
```

## Axis convention

`x` streamwise (inlet to outlet), `y` vertical (ground at `y=0`, lift along
`+y`), `z` spanwise. `axis_order` permutes the STL's axes into that frame.

## Validation

```bash
python benchmarks/smoke.py        # exact checks: Poisson, freestream, divergence
python benchmarks/blasius.py      # laminar boundary layer vs similarity solution
python benchmarks/cylinder.py     # bluff body: Cd, Strouhal, recirculation
python benchmarks/bluff_body.py   # square cylinder, sharp-edged separation
```

Results and error analysis: [VALIDATION_REPORT.md](VALIDATION_REPORT.md).

## Known limits

Read these before trusting a number. They are measured, not assumed. See
[VALIDATION_REPORT.md](VALIDATION_REPORT.md) for the data behind each.

- **Force magnitude is not grid-converged.** Cd moves *away* from a limit under
  refinement (1.610 -> 1.622 -> 1.779 across 16/24/32 cells/D on a Re=40
  cylinder) while the flow field in the same runs converges properly. Treat
  drag as +/-10%, and never compare two cases run at different resolutions:
  the bias is resolution dependent, so it stops being common-mode.
- **Absolute vehicle Cd is out of envelope.** The Model S case returns 0.557
  against a published 0.208. The run is stable and divergence-free; subgrid
  viscosity on an affordable grid reaches 1.2e4 x molecular, putting the
  *effective* Reynolds number near 850 instead of 1e7. Not fixable by
  debugging; it is the arithmetic of LES at this cell size.
- **A streamlined body is the worst case.** Low-drag shapes earn it through
  smooth pressure recovery over a long taper. At 138 mm cells that taper is a
  staircase drowned in eddy viscosity, and the body behaves bluff. Every
  benchmark that passes here is a bluff body.
- **No wall model.** The first cell off the body sits far outside the viscous
  sublayer at vehicle Reynolds numbers. Do not read friction drag as physical.
- **Uniform grid.** No local refinement near the body. Direct consequence of
  the FFT pressure solve, and the reason resolution is expensive.

## What it is good for

- Velocity and pressure fields (Blasius profile L2 error 0.7-1.3%)
- Wake structure, separation, recirculation length (converges into reference range)
- Shedding dynamics (Strouhal within 2.6% of reference)
- **A/B comparison of geometries at identical settings and identical resolution**

That last one is the intended workflow: rank designs, do not quote coefficients.

## Repository layout

```
waterfall_core/     solver package
  grid.py           staggered MAC grid; nz=1 collapses to 2D
  poisson.py        DCT pressure solver
  kernels.py        numba convection, diffusion, divergence, IBM, derived fields
  boundary.py       boundary conditions
  solver.py         fractional-step time loop
  geometry.py       STL -> solid volume fraction
  shapes.py         analytic geometry for benchmarks
  sgs.py            Vreman subgrid model
  config.py         YAML config, domain sizing, pre-flight
  io.py             frame output (npz) and VTK export
  runs.py           run registry and subprocess supervision (UI)
benchmarks/         verification cases with published reference values
cases/              example configurations
ui.py               local web app (FastAPI)
ui/                 front end (no build step)
waterfall-ui.bat    double-click launcher for Windows
legacy/             the previous ray-based estimator, kept for reference
```

## legacy/

`legacy/ray_solver.py` is the original ray-based estimator. It contained no
fluid physics (no momentum equation, no continuity, no pressure-velocity
coupling, and no Reynolds-number dependence at all) and could not run any of
the benchmarks above. It is kept only as a record of what the numbers in old
`summary.txt` files meant. See VALIDATION_REPORT.md for the details.