Particle API

Configuration for particle initialization and advection.

Key parameters:

• Spatial domain: xmin/max, ymin/max, z_level for initial particle placement
• Particle count: nxparticles × nyparticles
• Physics options: useybjw (vertical velocity), use3dadvection
• Timing control: particleadvectime - when to start advecting particles
• Interpolation: interpolation_method for velocity interpolation
• Boundaries: periodicx/y, reflectz for boundary conditions
• I/O: saveinterval and maxsave_points for trajectory output

Advanced timing control:

• particleadvectime=0.0: Start advecting immediately (default)
• particleadvectime>0.0: Keep particles stationary until this time
• Useful for letting flow field develop before particle release
• Enables study of transient vs established flow patterns

Note: Time integration uses Euler method for simplicity and stability.

Particle state including positions, global IDs, velocities, and trajectory history.

GLM Position Decomposition

Following the GLM framework, we track both:

• Mean positions (x, y, z): The Lagrangian-mean trajectory X(t) advected by QG flow
• Wave displacements (xix, xiy): The oscillatory displacement ξ from the wave field

The physical (wiggly) position is: xphys = x + xix, yphys = y + xiy

Wave Displacement (from PDF equations)

The wave displacement is computed from the wave velocity amplitude LA: ξx + iξy = Re{(LA / (-if)) × e^{-ift}} [eq. (6)]

where LA uses the YBJ operator L (NOT L⁺): L = ∂/∂z(f²/N² ∂/∂z) [eq. (4)] L⁺ = L - k_h²/4 (YBJ+)

Since B = L⁺A and L = L⁺ + kh²/4: LA = (L⁺ + kh²/4)A = B + (k_h²/4)A

Main particle tracker that handles both serial and parallel execution.

GLM Framework

The tracker advects mean positions X(t) using QG velocities and separately tracks wave displacement ξ for reconstructing physical positions x = X + ξ.

particles_in_box([T=Float32], z_level; x_max, y_max, x_min=0, y_min=0, nx=10, ny=10, kwargs...)

Create 2D particle distribution at a fixed z-level. Default precision is Float32 for memory efficiency.

Arguments

• T: Optional type parameter (Float32 or Float64). Default: Float32
• z_level: The z-level (depth) for all particles
• x_max, y_max: Maximum domain bounds (REQUIRED)
• x_min, y_min: Minimum bounds (default: 0.0)
• nx, ny: Number of particles in each direction (default: 10 each)

Examples

# Default Float32 precision (recommended for memory efficiency)
config = particles_in_box(500.0; x_max=G.Lx, y_max=G.Ly, nx=20, ny=20)

# Explicit Float64 if higher precision needed
config = particles_in_box(Float64, 500.0; x_max=G.Lx, y_max=G.Ly, nx=20, ny=20)

particles_in_circle(z; center, radius, n, pattern, kwargs...)

Create particles distributed in a circular disk at a fixed z-level.

Arguments

• z: Vertical level where particles are placed
• center: (x, y) center of circle (default: (π, π))
• radius: Radius of the circular region (default: 1.0)
• n: Number of particles (default: 100)
• pattern: Distribution pattern (default: :sunflower)
  • :sunflower - Fibonacci spiral (very uniform, recommended)
  • :rings - Concentric rings
  • :random - Uniform random within disk

Example

# 100 particles in a circle at z = -π/2
config = particles_in_circle(-π/2; radius=1.0, n=100)

# Custom center and larger circle
config = particles_in_circle(-1.0; center=(2.0, 2.0), radius=2.0, n=200)

# Random distribution
config = particles_in_circle(-π/4; radius=0.5, n=50, pattern=:random)

particles_in_grid_3d(; x_max, y_max, z_max, nx, ny, nz, x_min=0, y_min=0, z_min=nothing, kwargs...)

Create particles uniformly distributed in a 3D rectangular grid.

Arguments

• x_max, y_max, z_max: Domain bounds (REQUIRED - use G.Lx, G.Ly, G.Lz)
• x_min, y_min: Minimum bounds (default: 0.0)
• z_min: Minimum z (default: nothing → uses full depth with z ∈ [-z_max, 0])
• nx, ny, nz: Number of particles in each direction

Example

# 1000 particles in a 10×10×10 3D grid
config = particles_in_grid_3d(x_max=G.Lx, y_max=G.Ly, z_max=G.Lz, nx=10, ny=10, nz=10)

# Custom subdomain
config = particles_in_grid_3d(x_max=250e3, y_max=250e3, z_max=-500.0, z_min=-2000.0, nx=8, ny=8, nz=5)

particles_in_layers(z_levels; x_max, y_max, nx, ny, x_min=0, y_min=0, kwargs...)

Create particles distributed in 2D grids at multiple z-levels.

Arguments

• z_levels: Vector of z-levels where particles are placed
• x_max, y_max: Domain bounds (REQUIRED - use G.Lx, G.Ly)
• x_min, y_min: Minimum bounds (default: 0.0)
• nx, ny: Number of particles per level in x and y (default: 10 each)

Example

# 3 layers at depths 1000m, 2000m, 3000m with 10×10 particles each
config = particles_in_layers([-1000.0, -2000.0, -3000.0]; x_max=G.Lx, y_max=G.Ly, nx=10, ny=10)

# Custom subdomain with 5 particles per side at each layer
config = particles_in_layers([-500.0, -1000.0, -1500.0]; x_max=250e3, y_max=250e3, nx=5, ny=5)

particles_random_3d(n; x_max, y_max, z_max, x_min=0, y_min=0, z_min=nothing, seed=1234, kwargs...)

Create randomly distributed particles in a 3D volume.

Arguments

• n: Number of particles
• x_max, y_max, z_max: Domain bounds (REQUIRED - use G.Lx, G.Ly, G.Lz)
• x_min, y_min: Minimum bounds (default: 0.0)
• z_min: Minimum z (default: nothing → uses full depth with z ∈ [-z_max, 0])
• seed: Random seed for reproducibility (default: 1234)

Example

# 500 random particles in the full domain
config = particles_random_3d(500; x_max=G.Lx, y_max=G.Ly, z_max=G.Lz)

# 1000 random particles in a subdomain
config = particles_random_3d(1000; x_max=250e3, y_max=250e3, z_max=-500.0, z_min=-2000.0)

particles_custom(positions; kwargs...)

Create particles at user-specified positions.

Arguments

• positions: Vector of (x, y, z) tuples

Example

# 4 particles at specific locations
config = particles_custom([(1.0, 1.0, 0.5), (2.0, 2.0, 1.0), (3.0, 1.5, 0.75), (1.5, 3.0, 1.25)])

initialize_particles!(tracker, config)

Initialize particles uniformly in specified region, handling both serial and parallel.

initialize_particles!(tracker, config3d)

Initialize particles using enhanced 3D configuration.

advect_particles!(tracker, state, grid, dt, current_time=nothing; params=nothing, N2_profile=nothing)

Advect particles using GLM (Generalized Lagrangian Mean) framework.

GLM Algorithm (following PDF Section 1.4)

1. Advect mean position X using QG velocity (Lagrangian-mean flow): X^{n+1} = X^n + Δt × u^L_QG(X^n, t^n)

2. Reconstruct wave displacement at (X^{n+1}, t^{n+1}): ξ^{n+1} = Re{(LA(X^{n+1}, t^{n+1}) / (-if)) × e^(-if t^{n+1})}

3. Physical position (available via particles.x + particles.xi_x, etc.): x^{n+1} = X^{n+1} + ξ^{n+1}

Parameters

• tracker: ParticleTracker instance
• state: Current fluid state (contains B, A fields for wave displacement)
• grid: Grid information
• dt: Time step
• current_time: Current simulation time (required for wave phase computation)
• params: Model parameters (QGParams). Required for f₀.
• N2_profile: Optional N²(z) profile for nonuniform stratification.

Note

The stored positions (particles.x, particles.y, particles.z) are the MEAN positions X. The wave displacements (particles.xix, particles.xiy) provide the oscillatory correction. The physical (wiggly) trajectory is x = X + ξ.

interpolate_velocity_at_position(x, y, z, tracker)

Interpolate velocity at particle position with advanced schemes and cross-domain capability.

write_particle_trajectories(filename, tracker; metadata=Dict())

Write complete particle trajectory history to NetCDF file.

read_particle_trajectories(filename) -> NamedTuple

Read particle trajectory history from NetCDF file.

Returns a NamedTuple with fields:

• x: Matrix of x positions (np × ntime)
• y: Matrix of y positions (np × ntime)
• z: Matrix of z positions (np × ntime)
• time: Vector of time values (ntime)
• particle_ids: Vector of particle identifiers (np)
• attributes: Dict of global attributes from the file

This is the inverse of write_particle_trajectories.

Example

# Write trajectories
write_particle_trajectories("particles.nc", tracker)

# Read them back
traj = read_particle_trajectories("particles.nc")
println("Number of particles: ", size(traj.x, 1))
println("Number of time steps: ", length(traj.time))
println("Initial x positions: ", traj.x[:, 1])

write_particle_snapshot(filename, tracker, time)

Write current particle positions to NetCDF file (single time snapshot).

write_particle_trajectories_by_zlevel(base_filename, tracker;
                                      z_tolerance=1e-6, metadata=Dict())

Write particle trajectories to separate files based on z-level.

This function groups particles by their initial z-level and saves each group to a separate NetCDF file. Useful for analyzing particles initialized at different depths independently.

Parameters:

• basefilename: Base name for output files (e.g., "particles" -> "particlesz1.23.nc")
• tracker: ParticleTracker instance with trajectory history
• z_tolerance: Tolerance for grouping particles by z-level (default: 1e-6)
• metadata: Additional metadata to include in all files

Returns: Dictionary mapping z-levels to filenames

Example:

# Initialize particles at multiple z-levels
config = create_layered_distribution(0.0, 2π, 0.0, 2π, [π/4, π/2, 3π/4], 4, 4)
tracker = ParticleTracker(config, grid, parallel_config)

# Run simulation...

# Save each z-level to separate file
files = write_particle_trajectories_by_zlevel("particles", tracker)
# Creates: particles_z0.785.nc, particles_z1.571.nc, particles_z2.356.nc

Available interpolation methods.

validate_particle_cfl(tracker, max_velocity, dt)

Check if timestep satisfies CFL condition for particle advection in parallel mode.

For Euler integration, particles move up to dt*max_velocity from their starting position. If this exceeds the halo region, interpolation will be inaccurate.

Returns true if timestep is safe, false if timestep may cause issues.

Warning

If this returns false, consider:

• Reducing dt
• Increasing halo_width (use higher-order interpolation which has wider halos)

Enhanced particle configuration supporting 3D distributions.

Domain bounds (xmax, ymax, zmax) are REQUIRED - no defaults. Use the Grid to get domain size: `xmax = G.Lx, ymax = G.Ly, zmax = G.Lz(depth). If you passzmin, thenzmaxis treated as a coordinate (typically ≤ 0). Single-level distributions can usezmin == zmax`.

Available particle distribution patterns.

initialize_particles_3d!(tracker, config3d)

Initialize particles using 3D configuration.