WaterLily.jl: an immersed boundary fluid solver in pure Julia

A Navier-Stokes solver written from scratch in Julia

WaterLily.jl is a fast and simple fluid simulator written in pure Julia, supported by libraries from the Julia scientific community, with the stated aim of accelerating and enhancing fluid simulations.

The method section is compact and specific. WaterLily solves the unsteady incompressible 2D or 3D Navier-Stokes equations on a Cartesian grid. The pressure Poisson equation is solved with a geometric multigrid method. Solid boundaries are modelled using the Boundary Data Immersion Method. The solver can run on serial CPU, multi-threaded CPU, or GPU backends.

Four choices define the project. Pure Julia, so there is no Python or C wrapper in the loop. Cartesian grid, so the geometry is described rather than meshed. Geometric multigrid, so the pressure solve scales with the grid rather than degrading. Boundary data immersion, so a moving body is a function of position and time rather than a remeshed surface.

The project publishes its research in the Computer Physics Communications journal, in a 2025 paper by G.D. Weymouth and B. Font titled WaterLily.jl: A differentiable and backend-agnostic Julia solver for incompressible viscous flow around dynamic bodies, at doi 10.1016/j.cpc.2025.109748, and the README asks that you cite it if you use the library for research. That paper is described as covering the main features of the solver along with benchmarking, validation and profiling results.

There is also a JuliaCon2024 talk linked from the README, and a separate WaterLily-Examples repository holding the larger set of examples, notebooks and tutorials.

Grid dimensions must be powers of two for the multigrid solver

The example defines a domain of n by m cells and the README makes a point that is easy to miss: n and m have been set to multiples of powers of 2, which is important when using the very fast geometric multi-grid solver.

That is a real constraint, not a stylistic note. A geometric multigrid solver coarsens the grid by factors of two, so a grid that is not a power of two multiple cannot be coarsened cleanly, and either the solver will not run on that shape or it will round down and leave a partially resolved region. Anyone porting a pipeline that generates arbitrary grid sizes needs to size the domain accordingly before anything else.

The same example shows how the domain and the flow parameters are expressed. The user sets the domain size, the boundary conditions, the fluid viscosity which determines the Reynolds number, and immerses solid obstacles:

using WaterLily
function circle(n,m;Re=100,U=1)
    # signed distance function to circle
    radius, center = m/8, m/2-1
    sdf(x,t) = √sum(abs2, x .- center) - radius

    Simulation((n,m),   # domain size
               (U,0),   # domain velocity (& velocity scale)
               2radius; # length scale
               ν=U*2radius/Re,     # fluid viscosity
               body=AutoBody(sdf)) # geometry
end

The Reynolds number is written as Re=U*radius/ν, where nu is the kinematic viscosity of the fluid, and the README notes explicitly that the Greek nu is U+03BD and not the Latin lowercase v, a detail that has cost people real debugging time when a symbol renders identically.

Bodies are signed distance functions, not meshes

The geometry model is the part of the design that makes the solver extensible. The circle in the example is defined by a signed distance function, and the AutoBody function uses automatic differentiation to infer the other geometric parameters of the body automatically.

The consequence stated in the README is that you replace the circle's distance function with any other and now you have the flow around something else, with the examples repository holding a donut and the Julia logo as such cases. For more complex geometries, ParametricBodies.jl defines a body using any parametric curve, such as a spline.

So there are three levels of geometry complexity, and each is a different amount of work for you. A signed distance function is one function of position and time. A parametric body handles shapes you can describe as a curve. A general triangle mesh is not what these entry points are built around, and the README does not describe a mesh based body in this section.

Automatic differentiation is doing real work here. The body needs a normal and a curvature estimate to impose the boundary condition, and deriving those from the distance function by hand is tedious and error prone, so the solver differentiates it for you.

One small discrepancy is worth noticing in the example itself. The prose says the circle is centered at (m/2,m/2) while the code sets center to m/2-1, so the shipped example and its description differ by one cell, which matters if you are checking a result against a published figure.

From a Simulation object to a force curve and a vorticity field

The example runs in steps. The function returns a Simulation, and you initialise it and advance it in time:

circ = circle(3*2^5,2^6)
sim_step!(circ)

After that, the flow fields are ordinary arrays, and you can plot them with the Julia plotting ecosystem:

using Plots
u = circ.flow.u[:,:,1] # first component is x
contourf(u') # transpose the array for the plot

The README's interpretation of that first frame is instructive: the velocity within the circle is zero, the velocity far from the circle is one, and there are accelerated and decelerated regions around the circle, and since sim_step! has only taken a single time step the initial flow looks similar to potential flow because the viscous boundary layer has not separated yet.

For quantitative work, a set of flow metric functions lives in src/Metrics.jl, and the example uses pressure_force to measure the force on the immersed body over a time range:

function get_forces!(sim,t)
    sim_step!(sim,t,remeasure=false)
    force = WaterLily.pressure_force(sim)
    force./(0.5sim.L*sim.U^2) # scale the forces!
end

# Simulate through the time range and get forces
time = 1:0.1:50 # time scale is sim.L/sim.U
forces = [get_forces!(circ,t) for t in time];

The vorticity field is computed with the curl of the velocity and plotted with flood, a convenience function in WaterLily for 2D flood plots:

# Use curl(velocity) to compute vorticity `inside` the domain
ω = zeros(size(u));
@inside ω[I] = WaterLily.curl(3,I,circ.flow.u)*circ.L/circ.U

# Plot it using WaterLily's Plots Extension
flood(ω,clims = (-10,10),border=:none)

That last call goes through WaterLily's Plots extension, which is why the repository has an ext/ directory, so the plotting hooks ship as a package extension rather than a hard dependency.

One codebase for serial CPU, threads and GPU

Backend portability comes from KernelAbstractions.jl. The README states that WaterLily uses KernelAbstractions.jl to multi-thread on CPU and run on GPU backends, and the implementation method and speed-up are documented in the project documentation, in a section the truncated README points to rather than repeats.

The structure that makes this work is that the numerical kernels are written once against the KernelAbstractions interface. The same source then runs on a serial CPU, on multiple CPU threads, or on a GPU, without separate code paths per backend and without the user choosing a backend in the simulation setup.

The repository layout supports that design. There is a src/ directory for the solver, an ext/ directory for the package extension that adds the plotting helpers, a test/ directory, a docs/ directory, and a Project.toml declaring the dependencies. The build system is Buildkite, with a .buildkite/ directory at the top level, alongside GitHub Actions workflows for continuous integration and a codecov badge for coverage on the master branch.

For a solver whose inner loops are the whole cost, this matters more than it would in an application. If your code runs on a laptop CPU during development and on a GPU for a parameter sweep, the portability is already handled, and the sweep does not need a second implementation.

Licence, release cadence and the gaps in the documentation

The licence field for the project is filed as a custom licence that GitHub cannot classify, and the terms are in the LICENSE.md file at the repository root, so read that file before you build on the code or redistribute it.

The release cadence is fast and the repository is not archived. The last push was on 2026-10-19, and the recent tags are v1.9.2 from 2026-10-12, v1.9.1 from 2026-10-02 and v1.9.0 from 2026-09-30. Three patch releases in three weeks is the profile of a solver that is being used and adjusted rather than a library that sits still, and it also means you should pin a version rather than track the master branch.

What the documentation does not cover is as important as what it does. The README does not list the supported Julia versions beyond the Project.toml, does not describe the numerical tolerances or time stepping options in detail, does not document how the solver handles moving bodies beyond the time dependence of the distance function, and does not give a stability limit in terms of a time step. The truncated README also points at a section for the implementation method and speed-up that is not reproduced here, which is where the performance characterisation lives.

The paper is the place to look for validation. The README says the 2025 paper describes the main features of the solver and provides benchmarking, validation and profiling results, and if you are evaluating this for research that needs to be defensible, that document rather than the README is the artefact to cite and check.

Editorial conclusion

WaterLily.jl suits a numerics person who wants a readable, differentiable and backend agnostic incompressible flow solver they can extend, and it suits a Julia team that wants fluid simulation without leaving the language. It is the wrong choice if you need a general purpose CFD package with validated turbulence models, if your grid cannot be sized to powers of two, or if you need a solver with an existing track record in production engineering. Before adopting it, read the 2025 Computer Physics Communications paper, which the README describes as covering benchmarking, validation and profiling results, and check that the multigrid constraint on your grid dimensions holds for your case.

Frequently asked questions

Official sources

1. Issues
2. README
3. Releases
4. WaterLily-jl/WaterLily.jl on GitHub