2dPointVortex

2dPointVortex is a dependency-light C++20 solver for two-dimensional point-vortex dynamics. It supports the infinite plane, a plane periodic in one direction, a square periodic box, and a circular disk; CPU/OpenMP, MPI, and NVIDIA CUDA executables use the same input format, integrators, output files, and checkpoints.

The solver uses direct all-pairs velocity sums: O(N^2) except for the doubly periodic O(N^2 M) image sum. It is a clear numerical reference and a practical tool for small-to-medium simulations; it is not a tree-code or FMM implementation.

Quick start

The CPU example needs only CMake 3.20+ and a C++20 compiler. From the repository root:

cmake -S . -B build/cpu -DCMAKE_BUILD_TYPE=Release \
  -DPOINT_VORTEX_MPI=OFF -DPOINT_VORTEX_CUDA=OFF
cmake --build build/cpu --parallel
ctest --test-dir build/cpu --output-on-failure
./build/cpu/point_vortex_cpu examples/quickstart.params

This advances a two-vortex infinite-plane case to time 0.1. It writes:

File	Contents
runs/quickstart/trajectory.csv	Saved positions, circulations, and velocities
runs/quickstart/diagnostics.csv	Invariants and conservation drift
runs/quickstart/checkpoints/	Restart checkpoints

Every simulation has one managed run directory. Missing directories are created automatically. Existing solver output is protected; use a new runDirectory for another experiment, or set overwriteRun true to replace the managed output in that directory.

params.txt is a larger periodic example with 400 vortices, adaptive integration, and dipole removal/reinjection.

Features

• Four geometries: infinite, periodic_x, periodic, and disk
• Fixed-step classical RK4 and adaptive Dormand-Prince 5(4) integration
• CPU serial/OpenMP, MPI, and CUDA velocity backends
• Text initial conditions, geometry-aware generator, CSV output, and restart checkpoints
• Optional close dipole removal, with reinjection in bounded periodic and disk domains
• CMake and Make builds, numerical tests, and analysis/movie tools

Repository layout

.
├── src/                  Solver, kernels, integrators, I/O, and backend implementations
├── initial_conditions/   Initial-condition generator and its guide
├── examples/             Small parameter files, including the quick start
├── tests/                C++ unit/audit and Python integration/backend tests
├── scripts/
│   ├── vortices.ipynb    Vortex-configuration plotting notebook
│   ├── diagnostics.ipynb Invariants, drift, and dipole-event notebook
│   ├── movie_vortices.py Vortex-configuration MP4/GIF renderer
│   ├── point_vortex_plotting.py  Shared readers and plotting helpers
│   ├── analysis/         Jupyter notebook for diagnostics and configuration figures
│   └── movie/            Legacy streaming CSV-to-MP4 renderer
├── runs/                 Generated managed run output (ignored by Git)
├── CMakeLists.txt        Primary cross-platform build configuration
├── Makefile              Lightweight alternative build workflow
└── params.txt            Full periodic-run example

src/ at a glance

Area	Files	Responsibility
Program driver	main.cpp, print.cpp	Loads a run, schedules output, and reports diagnostics
Model and diagnostics	vortex.h, compute.cpp/.h	Vortex storage, velocity kernels, Hamiltonian, and invariants
Time integration	timestep.cpp/.h	RK4 and adaptive DOPRI5 stepping
Execution backends	backend_cpu.cpp, backend_mpi.cpp, backend_cuda.cu, backend_common.cpp, backend.h	CPU/OpenMP, MPI, CUDA, and shared backend interface
Configuration and input	params.h, read.cpp/.h	Parameter parsing, validation, and initial-condition loading
Events and restart	dipole.cpp/.h, checkpoint.cpp/.h	Dipole handling and versioned checkpoint I/O
Benchmark	benchmark.cpp	Standalone infinite-plane kernel benchmark

The velocity kernel is deliberately separate from the timestepper, so a new geometry or faster kernel can be added without rewriting the integrators.

Build

CMake (recommended)

The basic build attempts optional MPI and CUDA targets when their toolchains are installed; otherwise it still builds the CPU executable.

cmake -S . -B build/release -DCMAKE_BUILD_TYPE=Release
cmake --build build/release --parallel
ctest --test-dir build/release --output-on-failure

Use a CPU-only build when MPI or CUDA should not be detected:

cmake -S . -B build/cpu -DCMAKE_BUILD_TYPE=Release \
  -DPOINT_VORTEX_MPI=OFF -DPOINT_VORTEX_CUDA=OFF

Useful configuration options:

Option	Default	Effect
POINT_VORTEX_OPENMP	ON	Enable OpenMP when the compiler supports it
POINT_VORTEX_MPI	ON	Build point_vortex_mpi when MPI is found
POINT_VORTEX_CUDA	ON	Build point_vortex_cuda when CUDA is found
POINT_VORTEX_CUDA_ARCHITECTURES	empty	Optional CUDA target, e.g. -DPOINT_VORTEX_CUDA_ARCHITECTURES=89

Make

make                 # CPU/OpenMP executable and initial-condition generator
make test            # C++ numerical tests
make test-output     # Python output/restart integration tests
make mpi             # MPI executable
make cuda            # CUDA executable
make benchmark        # Kernel benchmark

Make outputs are under build/make/; CMake outputs are under the selected build directory. They are alternative ways to build the same source tree.

Run a simulation

All backends receive a parameter file as their optional first argument; with no argument, they use params.txt.

./build/release/point_vortex_cpu run.params
mpirun -n 4 ./build/release/point_vortex_mpi run.params
./build/release/point_vortex_cuda run.params

Start with the CPU backend when checking a new input. The CPU executable uses OpenMP when it was compiled with support and the population is sufficiently large. Set numThreads in the parameter file, or use OMP_NUM_THREADS; a positive numThreads takes precedence.

MPI distributes target-vortex calculations across ranks while retaining the complete source state on each rank. Only rank zero writes output. In a hybrid MPI/OpenMP run, choose ranks × threads to fit the available CPU cores:

OMP_NUM_THREADS=8 mpirun -n 2 ./build/release/point_vortex_mpi run.params

The CUDA backend requires an NVIDIA GPU, driver, and CUDA toolkit. It keeps vortex state and RK4/DOPRI5 stages on the GPU between output events, avoiding per-stage host/device transfers. The host synchronizes state only for output, diagnostics, checkpoints, and enabled dipole processing. This uses additional GPU memory for the integration-stage buffers.

Parameter files

Each non-empty line is key value; # starts a comment. Keys are case-sensitive. Paths cannot contain whitespace. Invalid or unknown settings stop the run rather than being ignored.

# Minimal two-vortex run
N 2
initialCondition dipole
boundaryCondition infinite
integrator rk4
timeStep 0.001
endTime 0.1
outputTime 0.02
runDirectory runs/my-two-vortex-case

Most-used settings:

Setting	Values / default	Purpose
initialCondition	ring	random, ring, single, dipole, or file
N	100	Population for random/ring; ignored for fixed, file, and checkpoint input
boundaryCondition	infinite	infinite, periodic_x, periodic, or disk
integrator	dopri5	rk4 or adaptive dopri5
timeStep, endTime	0.001, 1.0	Initial/fixed step and final simulation time
outputTime	0.1	Trajectory interval; final state is always saved
diagnosticsTime, checkpointTime	outputTime	Optional independent output intervals
coreRadius	0.0	Infinite-plane regularization radius
numThreads	0	OpenMP thread count; zero defers to the runtime
initialConditionFile	unset	Required when initialCondition file; contains x y circulation rows
dipoleRemoval	false	Enable dipole-removal events
dipoleRemovalDistance	0.01	Lower separation cutoff
dipoleRemovalInterval	0.0	Zero checks every accepted step; positive values set an independent time interval
dipoleRemovalUpper	false	Also remove closest-first matched dipoles above the upper distance
dipoleRemovalUpperDistance	1.0	Upper separation cutoff; must exceed dipoleRemovalDistance
dipoleReinjection	none	none, independent, or paired; bounded domains only
restartFile	unset	Checkpoint to restore; overrides the initial condition
runDirectory	runs/default	Self-contained output root for this simulation
overwriteRun	false	Replace this directory's managed solver output

For a plane periodic only in $x$, select periodic_x and set boxLengthX; $y$ remains unbounded and total circulation need not vanish. For a periodic box, set boxLengthX, boxLengthY (currently equal), and optionally periodicImageLayers; total circulation must be zero. For a disk, set diskRadius; every vortex must remain strictly inside it. absoluteTolerance, relativeTolerance, minimumTimeStep, and maximumTimeStep control adaptive DOPRI5. See the commented params.txt for every supported setting, including dipole removal and reinjection.

Equations of motion

Vortex $i$ has position $\boldsymbol r_i=(x_i,y_i)$ and constant circulation $\Gamma_i$. For every geometry the solver advances

$$\frac{d\boldsymbol r_i}{dt} =\sum_j\Gamma_j\,\boldsymbol K_\Omega (\boldsymbol r_i,\boldsymbol r_j), \qquad i=1,\ldots,N,$$

where the kernel $\boldsymbol K_\Omega$ is selected by boundaryCondition.

Infinite plane

For boundaryCondition infinite, the implemented regularized Biot–Savart equations are

$$\begin{aligned} \frac{dx_i}{dt} &=-\frac{1}{2\pi}\sum_{j\ne i}\Gamma_j \frac{y_i-y_j}{r_{ij}^2+\varepsilon^2},\\\ \frac{dy_i}{dt} &= \frac{1}{2\pi}\sum_{j\ne i}\Gamma_j \frac{x_i-x_j}{r_{ij}^2+\varepsilon^2},\\\ r_{ij}^2&=(x_i-x_j)^2+(y_i-y_j)^2. \end{aligned}$$

Here $\varepsilon=\texttt{coreRadius}$. Setting $\varepsilon=0$ gives the singular point-vortex model. The Hamiltonian reported by the infinite plane diagnostics is

$$H=-\frac{1}{4\pi}\sum_{i<j}\Gamma_i\Gamma_j \log\!\left(r_{ij}^2+\varepsilon^2\right).$$

Periodic in one direction

For boundaryCondition periodic_x, the $x$ direction has period $L=\texttt{boxLengthX}$ and the $y$ direction is unbounded. This is the singly periodic plane, topologically a cylinder. Define

$$X_{ij}=\frac{2\pi}{L}\operatorname{remainder}(x_i-x_j,L), \qquad Y_{ij}=\frac{2\pi}{L}(y_i-y_j).$$

Summing the infinite row of periodic images with $\sum_n(z+nL)^{-1}=(\pi/L)\cot(\pi z/L)$ gives the exact real-valued kernel

$$\begin{aligned} \frac{dx_i}{dt} &=-\frac{1}{2L}\sum_{j\ne i}\Gamma_j \frac{\sinh Y_{ij}}{\cosh Y_{ij}-\cos X_{ij}},\\\ \frac{dy_i}{dt} &= \frac{1}{2L}\sum_{j\ne i}\Gamma_j \frac{\sin X_{ij}}{\cosh Y_{ij}-\cos X_{ij}}. \end{aligned}$$

No image truncation or zero-total-circulation constraint is needed, so this kernel costs $O(N^2)$. The implementation uses scaled expressions at large $|Y_{ij}|$ to avoid hyperbolic-function overflow. Its reported Hamiltonian is

$$H=-\frac{1}{4\pi}\sum_{i<j}\Gamma_i\Gamma_j \log\!\left(\cosh Y_{ij}-\cos X_{ij}\right),$$

which is defined up to the usual circulation-dependent additive constant.

Square periodic box

For boundaryCondition periodic, the current implementation requires $L_x=L_y=L$ and $\sum_i\Gamma_i=0$. Define

$$\kappa=\frac{2\pi}{L}, \qquad X_{ij}=\kappa\,\operatorname{remainder}(x_i-x_j,L), \qquad Y_{ij}=\kappa\,\operatorname{remainder}(y_i-y_j,L).$$

The truncated Weiss–McWilliams image sum used by the code is

$$\begin{aligned} \frac{dx_i}{dt} &=-\frac{1}{2L}\sum_j\Gamma_j\sin Y_{ij} \sum_{n=-M}^{M} \frac{1}{\cosh(X_{ij}-2\pi n)-\cos Y_{ij}},\\\ \frac{dy_i}{dt} &= \frac{1}{2L}\sum_j\Gamma_j\sin X_{ij} \sum_{n=-M}^{M} \frac{1}{\cosh(Y_{ij}-2\pi n)-\cos X_{ij}}. \end{aligned}$$

The singular $j=i,n=0$ contribution is omitted. Here $L=\texttt{boxLengthX}=\texttt{boxLengthY}$ and $M=\texttt{periodicImageLayers}$. Larger $M$ retains more periodic image layers at greater $O(N^2M)$ cost. This is the Weiss–McWilliams square-torus construction.

Circular disk

For boundaryCondition disk, let $R=\texttt{diskRadius}$, require $|\boldsymbol r_i|&lt;R$, and define the inverse image

$$\boldsymbol r_j^*=\frac{R^2}{|\boldsymbol r_j|^2}\boldsymbol r_j.$$

With $\boldsymbol J(a,b)=(-b,a)$, the circle-theorem velocity is

$$\frac{d\boldsymbol r_i}{dt} =\frac{1}{2\pi}\boldsymbol J\!\left[ \sum_{j\ne i}\Gamma_j \frac{\boldsymbol r_i-\boldsymbol r_j} {|\boldsymbol r_i-\boldsymbol r_j|^2} -\sum_j\Gamma_j \frac{\boldsymbol r_i-\boldsymbol r_j^*} {|\boldsymbol r_i-\boldsymbol r_j^*|^2} \right].$$

The second sum includes the vortex's own opposite-sign image and enforces an impermeable circular wall. The implementation evaluates this term in an algebraically equivalent form that stays finite when a source is at the disk center.

Parameters and non-Hamiltonian events

Symbol	Parameter key	Meaning
$N$	N	Built-in initial vortex count
$\Gamma_i$	third initial-condition column	Circulation of vortex $i$
$\varepsilon$	coreRadius	Infinite-plane regularization radius
$L_x,L_y$	boxLengthX, boxLengthY	Periodic lengths; periodic_x uses only $L_x$
$M$	periodicImageLayers	Doubly periodic image-sum truncation
$R$	diskRadius	Circular-domain radius
$\Delta t,t_{\mathrm{end}}$	timeStep, endTime	Initial/fixed step and final time
tolerances	absoluteTolerance, relativeTolerance	Adaptive DOPRI5 error controls

There is no continuous forcing or viscous damping term in these ODEs. dipoleRemoval, dipoleRemovalDistance, and dipoleReinjection instead define discrete population events. dipoleRemovalInterval selects their schedule: zero checks after every accepted timestep, while a positive value uses an independent physical-time interval. The integrator lands exactly on scheduled removal times. These events can change circulation moments and the Hamiltonian; they are recorded in the diagnostics and should not be interpreted as part of the conservative point-vortex equations.

When dipoleRemovalUpper is enabled, the closest-first matching also removes opposite-sign pairs farther apart than dipoleRemovalUpperDistance. This models large-scale dissipation when like-signed vortices have clustered into separated regions. The same dipoleReinjection setting applies to lower- and upper-cutoff events. Diagnostics report their combined count in removed_pairs and the upper-cutoff subset in removed_upper_pairs.

Without removal events, circulation and Hamiltonian are reported for every geometry. The diagnostics additionally treat both linear impulses as relevant for infinite, periodic_x, and periodic, and angular impulse as relevant for infinite and disk. CSV files retain every invariant column so their schema remains identical across geometries; console and notebook displays select the geometry-appropriate subset.

Initial conditions

Built-in initial conditions are geometry-aware. random samples inside the periodic box or disk (and uses finite sampling extents in unbounded directions), while ring, single, and dipole use geometry-appropriate scales. Every generated state is checked against the selected domain.

Provide a plain text file with one x y circulation row per vortex (whitespace or commas are accepted), then select file and set initialConditionFile:

# initial.dat
-1.0, 0.0,  1.0
 1.0, 0.0, -1.0

initialCondition file
initialConditionFile initial.dat

File-loaded states are validated as well: periodic coordinates must lie in the fundamental box, disk coordinates must lie strictly inside the circle, and fully periodic states must have zero total circulation. Invalid input stops before the run begins.

Or build and use the generator:

cmake --build build/release --target point_vortex_initial --parallel
./build/release/point_vortex_initial \
  --geometry periodic --case random --count 400 --seed 20261376 \
  --box-length 2 --min-separation 0.01 --output runs/periodic_n400/initial_n400.dat

The generator supports all four geometries and the single, pair, dipole, ring, and random cases, records geometry metadata, and refuses incompatible solver settings. Full options and examples are in initial_conditions/README.md.

Output, run records, and restarting

Every simulation writes to a self-contained run directory. runDirectory defaults to runs/default; give each experiment a descriptive directory:

runDirectory runs/periodic_n400

The solver creates this fixed layout, which keeps each experiment's outputs together.

runs/periodic_n400/
├── trajectory.csv
├── diagnostics.csv
├── checkpoints/
├── resolved_parameters.txt
└── segments/
    └── segment_00000001/resolved_parameters.txt

resolved_parameters.txt records the validated settings, resolved output paths, selected backend, and available runtime details (OpenMP threads, MPI ranks, or CUDA device). Every fresh run, restart, or branch gets a new numbered segment record, retaining the provenance of the invocation even when the top-level record is updated. A run directory that already contains solver output is rejected by default. Set overwriteRun true only when intentionally replacing its trajectory, diagnostics, checkpoints, and provenance records; unrelated files in that directory are not removed.

Every backend writes the same portable formats:

Output	Location	Notes
Trajectory	runDirectory/trajectory.csv	time,frame,index,x,y,circulation,u,v rows
Diagnostics	runDirectory/diagnostics.csv	Invariants, drift, and dipole-event counts
Checkpoints	runDirectory/checkpoints/checkpoint_*.dat	Versioned restart state

Trajectory, diagnostics, and checkpoint intervals are simulation time, not wall-clock time. The solver always saves the initial and final states.

To branch from a checkpoint, create a new parameter file with a new run directory:

restartFile runs/periodic_n400/checkpoints/checkpoint_00000005.dat
endTime 2.0
runDirectory runs/periodic_n400_branch

The geometry, integrator, core radius, and dipole settings must match the checkpoint. The restart begins new CSV files; it does not append to the source trajectory. frame is a monotonically increasing output-event identifier stored in trajectory, diagnostics, and checkpoints. It lets analysis join streams reliably when their independent schedules coincide; gaps in an individual CSV mean that stream was not scheduled at that event.

Analysis and movies

The plotting notebooks read a complete run directory, infer its geometry and box from resolved_parameters.txt, and write figures beneath that run by default. Open the notebooks and edit their clearly marked Configuration cells:

jupyter lab scripts/vortices.ipynb scripts/diagnostics.ipynb
python3 scripts/movie_vortices.py --run-dir runs/periodic_n400 \
  --output runs/periodic_n400/figures/vortices.mp4

The configuration notebook can override the geometry, the singly periodic length, square or rectangular periodic display box, disk radius, and explicit viewing limits. The movie exposes equivalent command-line options. Frame selection, GIF/MP4 output, diagnostics ranges, smoothing, dependencies, and more examples are documented in scripts/README.md. The original analysis notebook and movie entry point remain under scripts/analysis/ and scripts/movie/ for compatibility.

Test and validate changes

Run the CMake test suite after changes to solver code:

ctest --test-dir build/release --output-on-failure

Optional backend comparisons require the corresponding executable/toolchain:

python3 tests/backend_consistency.py \
  ./build/release/point_vortex_cpu ./build/release/point_vortex_cuda
python3 tests/backend_consistency.py \
  ./build/release/point_vortex_cpu mpirun -n 2 ./build/release/point_vortex_mpi

tests/tests.cpp covers core numerical behavior; tests/audit_tests.cpp targets numerical edge cases; tests/output_integration.py covers output and restart workflows. The analysis and movie smoke tests can be run with python3 tests/tooling_tests.py build/release when their Python dependencies are installed.

Limitations

• Velocity evaluation is direct O(N^2) (O(N^2 M) for the doubly periodic image sum); MPI still replicates the source arrays on each rank.
• CUDA keeps velocity evaluation and RK4/DOPRI5 integration on the device; diagnostics and file I/O remain host-side, and the stage buffers increase GPU-memory use.
• Doubly periodic dynamics currently requires a square, zero-net-circulation domain; periodic_x has neither restriction on the unbounded direction nor a neutrality requirement.
• Dipole reinjection is not defined for the unbounded infinite and periodic_x geometries.
• Core regularization is available only for the infinite plane.
• Singular encounters, disk-boundary violations, and non-finite states stop the run.

License and citation

Copyright (c) 2022–2026 Jason Laurie. This project is distributed under the BSD 3-Clause License. Third-party dependencies remain subject to their own license terms.

If this software contributes to research or a publication, please cite it using the metadata in CITATION.cff.