Particles and dye of a simulation¶

Tracer particles and a numerical dye that live inside a native simulation: they advance with each step of sim.run, on the CPU or the GPU, and write files in the cn.io layout.

Tracers¶

import numpy as np

import conatus as cn

box = cn.Box(x=cn.periodic(4.0), y=cn.walls(-1, 1))
m = cn.Model(box, name="plane Couette", nu=0.1)
m.velocity("u"), m.pressure("p")
m.equations("dt(u) + grad(p) - nu*lap(u) = -u@grad(u)", "div(u) = 0")
m.at("bottom", u=(-1, 0))
m.at("top", u=(1, 0))
sim = m.simulate(N=(32, 32), scheme="RK443", dt=0.01)
sim.initial(u=lambda x, y: (y, 0 * x))
tracers = sim.particles(np.array([[0.0, -1.0], [2.0, 0.3]]), path="tracers.h5", every=50,
                        velocity="u", at_walls="refuse")
sim.run(until=3.0)
print(tracers.explain())
X = tracers.positions(unwrapped=True)       # (P, d), coordinate order

What is computed, and declared by explain():

Velocity at a tracer: the spectral series of the velocity evaluated exactly at the point (all resolved modes: cos/sin along periodic axes, Chebyshev \(T_n\) across the walls), from the coefficients of the state; never an interpolation of grid values. The cost is \(P \times\) the number of modes per evaluation (four per step), in chunks of particles.
Time integration: one classical RK4 step per flow step. Its end velocities are the flow states before and after the step; the midpoint velocity is interpolated in time through the last four flow states. Order 4 in a steady flow; in an unsteady one, never better than the order of the flow scheme itself, whose time error the flow states carry (3 with RK443), and the first step of each run interpolates linearly between two states (also third order).
Periodic axes: positions are kept unwrapped (positions(unwrapped=True), for dispersion) and also given wrapped into \([a, a + L)\).
Walls: a step that ends outside the walls raises a NumericalError naming the particle (at_walls="refuse", the default) or puts it back on the wall and counts it (at_walls="clip", wall_clips(), in the file and in explain()). A tracer on a no-slip wall stays on it; a tracer on a moving wall moves with it. RK4 stage positions are not clipped.
File (path=, every=): position, position_unwrapped, velocity with dimensions (t, particle, component) and wall_clips with (t, particle); cn.io.read(path).to_xarray().
Not covered (refused with the reason): curved and mapped geometries, sine/cosine walls, half lines and lines, compound intervals, complex-valued problems (periodic(L, complex=True)). Positions are not in the simulation’s checkpoints: after a restart, declare the tracers again (from the last row of the file to continue them).

Dye¶

import numpy as np

import conatus as cn

box = cn.Box(x=cn.periodic(4.0), y=cn.walls(-1, 1))
m = cn.Model(box, name="plane Couette with dye", nu=0.1)
m.velocity("u"), m.pressure("p")
m.equations("dt(u) + grad(p) - nu*lap(u) = -u@grad(u)", "div(u) = 0")
m.at("bottom", u=(-1, 0))
m.at("top", u=(1, 0))
m.dye("c", velocity="u", kappa=0.01)
sim = m.simulate(N=(32, 32), scheme="RK443", dt=0.01)
sim.initial(u=lambda x, y: (y, 0 * x), c=lambda x, y: np.exp(-((x - 2) ** 2 + y ** 2) / 0.3))
out = sim.dye("c", "dye.h5", every=50)
sim.run(until=3.0)
mass = cn.io.read("dye.h5").tasks["mass"]
print(out.explain())
print("relative mass drift:", abs(mass[-1] - mass[0]) / abs(mass[0]))

m.dye("c", velocity="u", kappa=...) declares one more unknown and its equation in conservative form, \(\partial_t c - \kappa \nabla^2 c = -\nabla\cdot(c\,\boldsymbol u)\), with the no-flux condition \(\partial_n c = 0\) on the walls when \(\kappa > 0\) (both printed by m.explain()). The dye is advanced by the simulation’s own scheme: diffusion implicit, advection explicit, the same dealiasing and the same time step, on the CPU or the GPU. kappa is a number or the name of a declared parameter. kappa=0 is pure advection: no diffusion, no filter, no hyperviscosity, no wall condition; the dye’s spectral tail then grows as its filaments thin, and the file records it.

sim.dye(name, path, every=) writes the dye (grid values; field=False for the diagnostics only) and, at each row: mass (\(\int c\,dV\) of the series), variance (\(\int (c - \bar c)^2\,dV\)), minimum and maximum of the grid values (under- and overshoots of an under-resolved dye), and tail (largest relative coefficient among the last two modes of each axis). Mass is conserved by the equation; what the discretization keeps is measured: to rounding in a periodic box, to the level of the tau correction of the no-flux condition with walls (reported in the validation table of the paper).

On the GPU¶

m.simulate(..., backend="cuda"): the tracers’ positions and the evaluation of the series stay on the device; positions() and the files copy to the host. The test tests/model/test_particles_gpu.py compares GPU and CPU runs. Not yet verified on GPU hardware: that test has not run yet, so treat the GPU path of tracers and dye as experimental until it has.

Two interfaces¶

	`sim.particles`, `sim.dye` (this page)	`cn.diagnostics.Particles`, `Dye` (Diagnostics and visualization)
Velocity	the running simulation, its exact spectral series	any callable, or grid snapshots with a declared interpolation
Time step	one RK4 step per flow step, past flow states for the midpoint	the caller’s `dt`, RK4 or RK2
Walls	`at_walls="refuse"` or `"clip"` (counted)	`boundary="raise"` or `"clip"` (counted)
Output	HDF5 file (`cn.io.read`) and `positions()`	arrays in memory

Use sim.particles while a native simulation runs; use the diagnostics module for a velocity you own or for data after a run.

See also: Diagnostics and visualization (tracers from snapshots or analytic velocities, inertial particles, streaklines, timelines, FTLE, mixing), tutorial 16. Lagrangian tracers and a numerical dye.