Lagrangian particle tracking
Models can keep track of the location and properties of neutrally buoyant particles. Particles are advected with the flow field using forward Euler time-stepping at every model iteration.
Simple particles
If you just need to keep of particle locations $(x, y, z)$ then you can construct some Lagrangian particles using the regular LagrangianParticles
constructor
grid = RectilinearGrid(size=(10, 10, 10), extent=(1, 1, 1))
Nparticles = 10
x₀ = zeros(Nparticles)
y₀ = rand(Nparticles)
z₀ = -0.5 * ones(Nparticles)
lagrangian_particles = LagrangianParticles(x=x₀, y=y₀, z=z₀)
# output
10 LagrangianParticles with eltype Particle:
├── 3 properties: (:x, :y, :z)
├── particle-wall restitution coefficient: 1.0
├── 0 tracked fields: ()
└── dynamics: no_dynamics
then pass it to a model constructor
model = NonhydrostaticModel(grid=grid, particles=lagrangian_particles)
# output
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 10×10×10 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: Centered reconstruction order 2
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
├── coriolis: Nothing
└── particles: 10 LagrangianParticles with eltype Particle and properties (:x, :y, :z)
Remember to use CuArray
instead of regular Array
when storing particle locations and properties on the GPU.
Custom particles
If you want to keep track of custom properties, such as the species or DNA of a Lagrangian particle representing a microbe in an agent-based model, then you can create your own custom particle type and pass a StructArray
to the LagrangianParticles
constructor.
using Oceananigans
using StructArrays
struct LagrangianMicrobe{T, S, D}
x :: T
y :: T
z :: T
species :: S
dna :: D
end
Nparticles = 3
x₀ = zeros(Nparticles)
y₀ = rand(Nparticles)
z₀ = -0.5 * ones(Nparticles)
species = [:rock, :paper, :scissors]
dna = ["TATACCCC", "CCTAGGAC", "CGATTTAA"]
particles = StructArray{LagrangianMicrobe}((x₀, y₀, z₀, species, dna));
lagrangian_particles = LagrangianParticles(particles)
# output
3 LagrangianParticles with eltype LagrangianMicrobe:
├── 5 properties: (:x, :y, :z, :species, :dna)
├── particle-wall restitution coefficient: 1.0
├── 0 tracked fields: ()
└── dynamics: no_dynamics
Not all data types can be passed to GPU kernels. If you intend to advect particles on the GPU make sure particle properties consist of only simple data types. The symbols and strings in this example won't work on the GPU.
Writing particle properties to disk
Particle properties can be written to disk using JLD2 or NetCDF.
When writing to JLD2 you can pass model.particles
as part of the named tuple of outputs.
JLD2OutputWriter(model, (; particles=model.particles), filename="particles", schedule=TimeInterval(15))
JLD2OutputWriter scheduled on TimeInterval(15 seconds):
├── filepath: build/model_setup/particles.jld2
├── 1 outputs: particles
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 19.4 KiB
When writing to NetCDF you should write particles to a separate file as the NetCDF dimensions differ for particle trajectories. You can just pass model.particles
straight to NetCDFOutputWriter
:
NetCDFOutputWriter(model, model.particles, filename="particles.nc", schedule=TimeInterval(15))
NetCDFOutputWriter scheduled on TimeInterval(15 seconds):
├── filepath: build/model_setup/particles.nc
├── dimensions: particle_id(3), time(0)
├── 1 outputs: particles
└── array type: Array{Float64}
├── file_splitting: NoFileSplitting
└── file size: 10.2 KiB
NetCDF does not support arbitrary data types. If you need to write custom particle properties to disk that are not supported by NetCDF then you should use JLD2 (which should support almost any Julia data type).