# WaterfallCFD

A fast incompressible CFD solver for external aerodynamics that turns an STL into a drag number, a pressure field and a wake video on a desktop CPU, with a local web UI.

## Overview

WaterfallCFD solves the incompressible Navier-Stokes equations on a uniform staggered grid with an immersed boundary. It is 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, and run again. It is not a substitute for a wall-resolved commercial solve. The places where it is wrong are measured and written down in VALIDATION_REPORT.md rather than left for the user to find. The intended workflow is A/B comparison of geometries at identical settings and identical resolution: rank designs, do not quote coefficients.

## How it works

### Governing equations

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

Incompressible, isothermal, Newtonian, constant density.

### Numerical choices

| 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 per 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; an 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. The compute kernels are written in Numba; they are the solver, not an optional speed-up.

### Axis convention

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

### Code 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)
legacy/             the previous ray-based estimator, kept for reference
```

### Validation

The benchmark suite covers exact checks (Poisson, freestream, divergence), a laminar Blasius boundary layer against the similarity solution, a circular cylinder (Cd, Strouhal number, recirculation length) and a square cylinder for sharp-edged separation. Measured results:

- Velocity and pressure fields: Blasius profile L2 error 0.7-1.3%
- Wake structure, separation and recirculation length converge into the reference range
- Shedding dynamics: Strouhal number within 2.6% of reference
- Drag: treat as +/-10%

### Legacy ray estimator

`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. It could not run any of the benchmarks. It is kept only as a record of what the numbers in old `summary.txt` files meant.

## Status

The current revision is a complete rewrite as a CPU Navier-Stokes LES solver with a local web UI. Validation confirms the flow field against benchmarks and puts drag within +/-10%. Current revision: **o2.0A**.

## Build and use

### Install

Create a virtual environment and install the dependencies:

```bash
python -m venv .venv
source .venv/bin/activate        # Windows: .venv\Scripts\activate
python -m pip install numpy scipy trimesh matplotlib rtree numba scikit-image imageio imageio-ffmpeg pyyaml fastapi "uvicorn[standard]" python-multipart
```

Numba is required.

### The app

```bash
python ui.py
```

A browser opens at `http://127.0.0.1:8010`. The app listens on loopback only, so nothing is exposed to the network.

- **New run**: drag an STL onto the page. The app reads the mesh, guesses the flow axis, and draws side, plan and front previews with the flow direction marked. Enter the part's real length in metres and the scale factor is worked out automatically. Runtime and cell count update live as settings change.
- **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 the run produced, all downloadable.

**Always check the orientation preview before running.** The flow-axis guess assumes the longest dimension faces the flow. That 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.

Pre-flight only (geometry validity, orientation preview, domain bounds, blockage, cell count and estimated wall-clock time). Always run this first:

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

Full run, writing `summary.json` (Cd, Cl, forces, diagnostics) and a directory of frames:

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

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

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

### Configuration

See `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 settings from the command line without editing the file:

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

### Validation benchmarks

```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
```

## Known limits

- Force magnitude is not grid-converged: on a Re=40 cylinder Cd moves away from a limit under refinement (1.610, 1.622, 1.779 at 16/24/32 cells per diameter) while the flow field converges. Treat drag as +/-10% and never compare cases run at different resolutions, because the bias is resolution dependent.
- Absolute vehicle Cd is out of envelope: the Model S case returns 0.557 against a published 0.208. Subgrid viscosity on an affordable grid reaches 1.2e4 times molecular, putting the effective Reynolds number near 850 instead of 1e7.
- Streamlined bodies are the worst case: at 138 mm cells the smooth pressure-recovery taper becomes a staircase drowned in eddy viscosity, and the body behaves bluff. Every passing benchmark is a bluff body.
- There is no wall model: the first cell off the body sits far outside the viscous sublayer at vehicle Reynolds numbers, so friction drag is not physical.
- The grid is uniform with no local refinement near the body, a direct consequence of the FFT pressure solve and the reason resolution is expensive.
- An incorrect axis_order still produces a plausible drag number instead of an error, so orientation must be checked by the user.

## Revisions

- **o2.0A** (current): Complete rewrite as a CPU incompressible Navier-Stokes LES solver with immersed boundary, DCT pressure solve and local web UI; validated flow field, drag within +/-10%.
- **c1.1A** (closed): Final ray-casting estimator using Numba, sector sampling and a calibrated Cd; superseded because the ray method contained no fluid mechanics.
- **c1.0A** (closed): First ray-casting quick CFD drag estimator.

Revision codes read o (open) or c (closed), then major.minor, then the branch letter; A is the main line.

## License

PolyForm Noncommercial License 1.0.0. See [LICENSE](Legal/LICENSE).