The Problem

A sloshing tank is a canonical problem in computational fluid dynamics. Water in a rectangular container, shaken interactively, rendered in real time. It shows up in aerospace (fuel tanks on rockets), naval engineering (liquid cargo ships), and civil engineering (water towers under seismic load). The core challenge is simulating an incompressible fluid with a free surface that can slosh, break, and splash — while keeping the physics quantitatively correct.

The governing equations are the incompressible Navier-Stokes equations:

The second equation — the divergence-free condition — is the incompressibility constraint. It says fluid volume is conserved everywhere. Enforcing it numerically is the hard part.

For validation, the relevant benchmark is the first-mode sloshing frequency. Linear wave theory predicts:

where L is the tank length and h is the water depth. Any simulation that reproduces this frequency to within a few percent is doing real physics.

The approach: give Claude Code the goal and one technical constraint — MAC grid — then describe what you observe and let it diagnose. Everything else was left to Claude.

What Claude Proposed

The one constraint I gave upfront was the grid type: MAC (Marker-And-Cell). The MAC grid places pressure at cell centers and velocity components at cell faces, which eliminates the checkerboard pressure oscillations that plague collocated grids. It’s the standard choice for free-surface flows, and naming it upfront told Claude something about the level of rigor expected.

Given that constraint and the goal, Claude designed the full method stack without further input:

The fractional-step method works in three stages per timestep:

1. Advect the velocity field (semi-Lagrangian: trace each grid point backward along the flow, interpolate velocity at the departure point)
2. Add body forces (gravity, tank acceleration)
3. Project onto the divergence-free subspace by solving a pressure Poisson equation and subtracting the gradient

The VOF field — a scalar between 0 and 1 indicating water fraction per cell — is advected separately using an operator-split scheme with the Superbee limiter, which is the most compressive TVD limiter: it keeps interfaces sharp while remaining monotone.

The project structure:

sloshing-2d/
├── CMakeLists.txt
├── src/
│   ├── grid.hpp/cpp       # MAC staggered grid
│   ├── solver.hpp/cpp     # Navier-Stokes solver
│   ├── vof.hpp/cpp        # VOF transport
│   ├── renderer.hpp/cpp   # OpenGL visualization
│   └── simulation.hpp/cpp # Orchestrator
└── tests/
    ├── test_grid.cpp      # 7 test cases
    ├── test_solver.cpp    # 6 test cases
    ├── test_vof.cpp       # 5 test cases
    └── test_physics.cpp   # 2 physics validation tests

All 20 tests passed on the first build. The simulation ran and showed water in the tank.

The Bug Hunt

The initial build worked in the sense that water appeared and moved. But watching it for a few seconds revealed problems. This is where the experiment’s approach kicked in: rather than reading code or diagnosing algorithms myself, I described what was visible and let Claude diagnose.

The bugs did not appear one at a time, and they didn’t resolve cleanly. Volume loss and a fuzzy surface showed up together. Some fixes were partially effective and left residuals. Wall artifacts persisted through multiple attempts. The simulation only reached a state suitable for experiments when a final fix resolved the outstanding issues simultaneously.

Volume loss and a fuzzy surface

Without pressing any keys, the water level slowly dropped. Once forces were applied and the surface started moving, the volume loss accelerated and a second problem appeared simultaneously: the free surface smeared out over 5–10 cells with no sharp boundary.

Claude’s diagnosis for the volume loss: the operator-split VOF advection needed the Weymouth-Yue divergence correction. In 2D, \nabla \cdot \mathbf{u} = 0, but when advection is split into separate x and y sweeps, each 1D sweep sees a non-zero apparent divergence. Without correction, this causes artificial compression or expansion of the volume fraction. The fix divides each sweep’s flux by:

A second contributor: the pressure solver was being zeroed every timestep and given only 80 SOR iterations to reconstruct the hydrostatic profile from scratch. It wasn’t converging, leaving residual velocity divergence. Warm-starting — reusing the previous timestep’s pressure as the initial guess — fixed this. In near-equilibrium, the pressure barely changes between steps, so warm-starting effectively provides hundreds of accumulated iterations for free.

For the fuzzy surface, Claude replaced the first-order upwind VOF flux with a Superbee TVD limiter:

static double superbee(double r) {
    return std::max(0.0, std::max(std::min(2.0*r, 1.0), std::min(r, 2.0)));
}

This sharpened the interface to roughly 2 cells wide. A smoothstep in the renderer made it look even crisper without touching the physics:

f = f * f * (3.0 - 2.0 * f); // Hermite smoothstep — visual only

The prompt that drove this iteration:

“The water is still loosing volume, also when there is no force applied. The surface looks much better, only slightly fuzzy now.”

Volume loss remained partially present and continued to be addressed in later iterations.

Air patches in the wall cells

Light-colored patches appeared inside the solid boundary cells, particularly in corners — they looked like air intrusion rather than ghost water, though Claude read the screenshot the other way. The fix — initializing boundary cells to VOF=0 and calling enforce_vof_boundaries() after each advection step — didn’t fully take on the first attempt and had to be revisited in a second session before the artifacts disappeared.

Patterns forming in still water

Even with no keyboard input, subtle color variations appeared in the bulk fluid — structure that shouldn’t exist in a static tank.

This is spurious currents: parasitic velocities at the free surface driven by the mismatch between full-strength gravity and the pressure boundary condition p = 0 at the free surface. The solver can’t perfectly cancel gravity there, leaving a residual velocity field.

Claude’s fix: scale gravity by the local VOF value at each velocity face. This eliminated the visible patterns.

The fix for spurious currents broke the physics

With VOF-weighted gravity in place, the sloshing frequency validation test reported 11.4% error — the simulated oscillation was substantially slower than the analytical prediction. This was the first indication anything was wrong with the physics; nothing before the validation test had signalled a problem. Claude was prepared to accept the number and move on. It took pushing back — pointing out that 11% is far outside the range CFD results are held to — before the cause was investigated.

The mechanism: scaling gravity by the local VOF value (between 0 and 1 at the interface) weakens the effective restoring force at the surface. Since sloshing frequency goes as \sqrt{g_\text{eff}}, even a modest reduction in effective gravity slows the oscillation measurably. Every unit test still passed. Only the physics validation test caught it.

The fix: revert to full-strength gravity, applied based on cell type rather than VOF value. Any velocity face adjacent to at least one FLUID cell gets full gravity. The warm-started pressure solver converges well enough that the spurious currents are manageable without the VOF weighting. After this change, sloshing frequency error dropped to 5.9% and the simulation reached a usable state.

Results

All 20 tests pass:

GridTests ........... Passed  (7 test cases)
SolverTests ......... Passed  (6 test cases)
VofTests ............ Passed  (5 test cases)
PhysicsTests ........ Passed  (2 test cases)

Hydrostatic pressure: RMS error < 0.01%. The pressure at every depth matches p = \rho g (H – y) to 4+ significant digits:

  y [m]   | p_sim [Pa]  | p_analytical [Pa] | error%
  0.023    |  5655.89    |  5656.08          | -0.003%
  0.289    |  3050.10    |  3050.30          | -0.006%
  0.555    |   444.33    |   444.52          | -0.043%

First-mode sloshing frequency:

The measured period is longer than analytical — the expected direction for semi-Lagrangian advection, which introduces numerical diffusion that acts as artificial viscosity and slows the oscillation. Published CFD results for similar schemes at moderate resolution report 2–7%.

The oscillation trace shows clean sinusoidal behavior with gradual amplitude decay:

t=0.00  xcm= 0.0168    (initial tilt)
t=0.54  xcm= 0.0139    (peak right)
t=1.08  xcm=-0.0312    (peak left)
t=1.61  xcm=-0.0185    (returning)
t=2.15  xcm= 0.0248    (peak right, slightly damped)
t=2.69  xcm= 0.0217    (continuing)
t=3.12  xcm=-0.0014    (crossing zero)

Lessons Learned

Describing what you see is a surprisingly effective bug report. Every bug in this session was caught by visual observation and described in plain language: water disappearing, fuzzy surface, blue patches in walls, patterns in still water. Claude consistently identified the numerical root cause from the symptom. The symptom space in CFD is constrained enough that this works — most visual artifacts have known causes in the literature.

Tests that pass aren’t enough. Physics tests are required. The VOF-weighted gravity bug passed all unit tests. Only the sloshing frequency validation — a test against an analytical solution — caught it. For scientific code, behavioral correctness and physical correctness are different things.

Fixes don’t work cleanly and problems don’t appear sequentially. Volume loss and fuzzy surface appeared together. Some fixes were partial. Wall artifacts took multiple attempts. The simulation reached a usable state when the final gravity fix resolved everything at once, not before.

At this resolution, there’s room left. At 128×64, the frequency error is 5.9%. Doubling resolution would likely bring it under 3%. The SOR pressure solver becomes a bottleneck above 256×128 — conjugate gradient or multigrid would scale better. Known limitations, reasonable starting points for a next iteration.

Update: Dambreak with Baffle

A dambreak comparison against OpenFOAM’s interFoam revealed a weakness: the wavefront was advancing at 28% of Ritter’s analytical speed. Published 2D CFD lands between 70% and 95%.

Changes

Velocity extrapolation. The semi-Lagrangian back-trace routinely lands in the air region above the free surface, where nothing sets the velocity. Every back-trace into air returns zero, damping surface motion each step. Fix: propagate fluid velocities into the air band before each advection step, weighted by the VOF gradient along the interface normal (Foster & Fedkiw 2001).

Geometric VOF (PLIC). The algebraic VOF scheme (Superbee TVD) smeared air into the water bulk. During oscillation, surface cells lose VOF to numerical diffusion; on flow reversal, depleted values ratchet inward. This is structural to diffuse-interface tracking — no limiter tuning fixes it. CICSAM eliminated the smearing but introduced diagonal flotsam, a known NVD artefact. Fix: PLIC — reconstruct a straight line segment per interface cell from a Youngs-estimated normal, compute fluxes as exact polygon intersections. Bulk cells stay at exactly 0 or 1 by construction.

Boundary conditions. The SOLID border layer and explicit enforce_boundary_conditions() were replaced by folding wall BCs into the pressure Poisson equation as non-homogeneous Neumann. This exposed a CSF surface tension bug near internal obstacles (baffle VOF = 0 was being read as a fake water-air interface).

Numbers

The Ritter ratio at 0.71 sits inside the published range for 2D CFD. Volume conservation is exact to machine precision with PLIC.

Source Code: https://github.com/ai-sci-computing/sloshing-2d