Numerical Codes

A collection of numerical software for studying nonlinear waves, quantum fluids, and two-dimensional turbulence from continuous fields resolved in Fourier space to interacting systems of discrete vortices.

Each code is written in modern C++ and includes reproducible examples, diagnostics, restart support, and plotting and movie scripts. The repositories contain full installation and usage guides.

Isolated soliton field from the one-dimensional Schrödinger–Helmholtz model — Relaxation of a nonlinear Schrödinger soliton into a soliton bound state in the one-dimensional Schrödinger–Helmholtz model.

Nonlinear optical waves

1D Schrödinger–Helmholtz

A pseudo-spectral solver for the nonlocal 1D Schrödinger–Helmholtz equation and two closely related limits. It is designed to explore nonlocal optical wave turbulence and coherent solitonic structures in a periodic one-dimensional domain. The equation models the propagation of light, along the z axis, through a nematic liquid crystal cell.

Governing equation

i∂_zψ = c∂_xxψ+ gψ(1 − β∂_xx)⁻¹|ψ|²+ μψ

Here ψ(x,z) is a complex wave field and the inverse Helmholtz operator smooths the nonlinear density response over a scale set by β. Configurable forcing and large- and small-scale damping may also be included.

A single model setting also selects the long-wave approximation (1 + β∂_xx)|ψ|² or the local cubic nonlinear Schrödinger equation. All variants use a fully dealiased pseudo-spectral nonlinear evaluation and fixed-step exponential or integrating-factor time integrators.

Three optical wave models
Stochastic or single-mode forcing
FFTW and optional OpenMP

View the GitHub repository

Gross–Pitaevskii density field containing several quantised vortices — Gross–Pitaevskii density field showing Bogoliubov waves and quantised vortices.

Quantum fluids

2D Gross–Pitaevskii

A parallel pseudo-spectral solver for a complex 2D Gross–Pitaevskii equation on a doubly periodic domain. The model describes weakly interacting Bose gases and superfluids, as well as 2D nonlinear optical propagation defined through a Kerr nonlinearity. The model resolves both dispersive waves and quantised vortices within a single field equation.

Conservative governing equation

i∂_tψ = c∆ψ+ g|ψ|²ψ + μψ

The complex field ψ(x,t) contains the fluid density |ψ|² and phase. Optional Ginzburg–Landau damping, spectral dissipation, and external forcing extend the equation to driven, non-equilibrium systems.

The cubic nonlinearity is evaluated with a two-pass 3/2 dealiasing rule. CPU/OpenMP, hybrid MPI/OpenMP, and NVIDIA CUDA executables share the same model, inputs, output, and checkpoint format, making it possible to move a calculation between platforms.

Focussing/defocussing nonlinearity
Fully dealiased
CPU, MPI, and CUDA backends

View the GitHub repository

Complex vortical structures in a two-dimensional turbulent flow — Vorticity field of a two-dimensional Navier–Stokes turbulent flow.

Classical turbulence

2D Navier–Stokes

A parallel pseudo-spectral solver for incompressible 2D Navier-Stokes flow in vorticity form. It supports forced and decaying turbulence, ordinary or scale-selective dissipation, and an optional beta-plane term for idealised geophysical flows.

Governing equation

∂_tω + u · ∇ω + β∂_xψ= −ν(−∆)^pω− α(−∆)^qω + F

ω = ∆ψ,u = (−∂_yψ, ∂_xψ)

The vorticity ω determines the streamfunction ψ and incompressible velocity u. The coefficients ν and α control small- and large-scale dissipation, while F injects vorticity into selected Fourier modes.

Nonlinear advection is evaluated on a fully padded 3/2-rule grid. The solver records energy and enstrophy, their spectra and fluxes, and runs through matching CPU/OpenMP, hybrid MPI/OpenMP, and NVIDIA CUDA implementations.

Beta-plane dynamics
Energy and enstrophy diagnostics
CPU, MPI, and CUDA backends

View the GitHub repository

Positive and negative point vortices distributed across a periodic domain — A turbulent point-vortex flow displaying clustering of same-signed vortices.

Vortex dynamics

2D point vortices

A direct N-body solver for two-dimensional point-vortex dynamics. Instead of resolving a continuous velocity field, it follows idealised vortices with fixed circulations as they advect one another through the Biot–Savart interaction.

Equations of motion in the infinite plane

dx_idt= −12π∑_j≠i Γ_jy_i−y_jr_ij²+ε²dy_idt= 12π∑_j≠i Γ_jx_i−x_jr_ij²+ε²

Vortex i has position (x_i, y_i) and circulation Γ_i; ε optionally regularises close encounters. Geometry-specific kernels provide singly periodic, doubly periodic, and circular-disk domains as well as the infinite plane shown here.

The solver uses direct all-pairs velocity sums with adaptive Dormand–Prince time integration. CPU/OpenMP, MPI, and CUDA backends use identical input, trajectory, diagnostics, and checkpoint formats. The solver includes bespoke forcing and dissipation algorithms to study out-of-equilibrium point vortex motion.

Four domain geometries
Adaptive timestepping
CPU, MPI, and CUDA backends

View the GitHub repository