Using interpolate() to track lagrangian particle velocity/temperature over time

[sebasti4nmc]

Hi, my name is Sebastian and I'm an undergraduate student working on a summer research project investigating upper ocean mixing under Dr. Amrapalli Garanaik (@amrapallig) at Oregon State University.

I've implemented lagrangian particle tracking into my Oceananigans simulation with custom particles that have temperature and velocity properties, and I'm trying to see how I can track those fields on each individual particle using an interpolation function that I can pass to the "dynamics" keyword in LagrangianParticles(). However, when I try to run my simulation, I keep getting errors related to unsafe_trunc() and method errors related to my tracking function track_fields!(). Any idea how to fix this? Here is the setup for my particles:

### Lagrangian particles
using Oceananigans: VelocityFields
using Oceananigans.Fields: interpolate

const Nparticles = 10

struct Particle{FT}
    x :: FT
    y :: FT
    z :: FT
    u :: FT
    v :: FT
    w :: FT
    T :: FT
end

# Initial velocities + temps
u₀, v₀, w₀ = [zeros(Nparticles) for i in 1:3]
T₀ = fill(20.0, Nparticles)

# Reference depths
const z₁ = 5 # near surface
const z₂ = 10 # depth of mixed layer
const z₃ = 20 # thermocline 1
const z₄ = 28 # thermocline 2
const z₅ = 64 # bottom
const ϵ = 1e-6

# Initial positions
x₀ = [32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ,]
y₀ = [32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ,]
z₀ = [-z₁+ϵ, -z₁+ϵ, -z₂+ϵ, -z₂+ϵ, -z₃+ϵ, -z₃+ϵ, -z₄+ϵ, -z₄+ϵ, -z₅+ϵ, -z₅+ϵ]

# StructArray with initial conditions
particle_array = StructArray{Particle}((x₀, y₀, z₀, u₀, v₀, w₀, T₀))

# Define tracer + velocity fields for particles
u, v, w = VelocityFields(grid)

T = CenterField(grid)


# Define function to interpolate fields at particle locations
function track_fields!(ps, model)
    u_field = model.velocities.u
    v_field = model.velocities.v
    w_field = model.velocities.w
    T_field = model.tracers.T

    for i in eachindex(ps.x)
        xᵢ, yᵢ, zᵢ = ps.x[i], ps.y[i], ps.z[i]

        ps.u[i] = interpolate((xᵢ, yᵢ, zᵢ), u_field)
        ps.v[i] = interpolate((xᵢ, yᵢ, zᵢ), v_field)
        ps.w[i] = interpolate((xᵢ, yᵢ, zᵢ), w_field)
        ps.T[i] = interpolate((xᵢ, yᵢ, zᵢ),T_field)
    end
end

lagrangian_particles = LagrangianParticles(particle_array;
    tracked_fields = (; u, v, w, T),
    dynamics = track_fields!)

[jagoosw]

Hi @sebasti4nmc,

You shouldn't need to write a custom dynamics for tracking fields, if you omit the dynamics the tracked fields should get updated automatically.

(You haven't include the line where you make the model, but just in case make sure that you pass the fields u, v, w, and T to the model like model = NonhydrostaticModel(; grid, velocities = (; u, v, w), tracers = T) otherwise they won't get updated).

[sebasti4nmc]

Thanks for your help. I have done so and it looks like my velocities are tracking properly. However, my temperature still doesn't seem to be working right. I had set it up so that the temperature of all particles starts out as 20 degrees, but when I read the data from my simulation, it looks like it starts at that initial temperature for the t=0 timestep then immediately drops to 0 for the rest of the simulation. I'll add my model constructor as well in case there's any small mistakes I might have missed

grid = RectilinearGrid(CPU(), size=(64, 64, 80),   # Nx, Ny, Nz
                              extent=(128, 128, 80)) # Lx, Ly, Lz

### Lagrangian particles
using Oceananigans: VelocityFields
using Oceananigans.Fields: interpolate

const Nparticles = 10

struct Particle{FT}
    x :: FT
    y :: FT
    z :: FT
    u :: FT
    v :: FT
    w :: FT
    T :: FT
end

# Initial velocities + temps
u₀, v₀, w₀ = [zeros(Nparticles) for i in 1:3]
T₀ = fill(20.0, Nparticles)

# Reference depths
const z₁ = 5 # near surface
const z₂ = 10 # depth of mixed layer
const z₃ = 20 # thermocline 1
const z₄ = 28 # thermocline 2
const z₅ = 64 # bottom
const ϵ = 1e-6

# Initial positions
x₀ = [32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ,]
y₀ = [32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ, 32+ϵ, 96+ϵ,]
z₀ = [-z₁+ϵ, -z₁+ϵ, -z₂+ϵ, -z₂+ϵ, -z₃+ϵ, -z₃+ϵ, -z₄+ϵ, -z₄+ϵ, -z₅+ϵ, -z₅+ϵ]

# StructArray with initial conditions
particle_array = StructArray{Particle}((x₀, y₀, z₀, u₀, v₀, w₀, T₀))

# Define tracer + velocity fields for particles
u, v, w = VelocityFields(grid)

T = CenterField(grid)

lagrangian_particles = LagrangianParticles(particle_array;
    tracked_fields = (; u, v, w, T))

model = NonhydrostaticModel(; grid,
                            velocities = (; u, v, w),
                            tracers = (:T, :S),
                            particles = lagrangian_particles)

[jagoosw]

Do you set the model temperature or the temperature field directly after you construct the model like set!(model, T = 20) or set!(T, 20)? Because if not then the default value is 0 and the particles are correctly tracking it no?

[sebasti4nmc]

Yes, I have boundary conditions that I set before constructing the model, and initial conditions that I set after constructing the model. I'll add those below. Let me know if there's anything else that would be helpful to know.

Qʰ = 5.0  # W m⁻², surface heat flux  
ρₒ = 1026.0 # kg m⁻³, average density at the surface of the world ocean
cᴾ = 3991.0 # J K⁻¹ kg⁻¹, typical heat capacity for seawater
Qᵀ = Qʰ / (ρₒ * cᴾ) # K m s⁻¹, surface temperature flux
dTdz = 0.01 # K m⁻¹
T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵀ),
                                bottom = GradientBoundaryCondition(dTdz))

# In the previous version I posted, I left out all the other kwargs for simplicity, but I'll put them all here. 
model = NonhydrostaticModel(; grid,
                            buoyancy = SeawaterBuoyancy(),
                            advection = UpwindBiased(order=5),
                            timestepper = :RungeKutta3,
                            velocities = (; u, v, w),
                            tracers = (:T, :S),
                            coriolis = FPlane(f = 1e-4),
                            closure = SmagorinskyLilly(),
                            stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                            boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs),
                            forcing = (;u = u_forcing, v = u_forcing, w = u_forcing),
                            particles = lagrangian_particles)

### Initial conditions
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise
const initial_mixed_layer_depth = 33.0 # m
stratification(z) = z < - initial_mixed_layer_depth ? 20 + dTdz * z : 20 + dTdz * (-initial_mixed_layer_depth)

Tᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(z) * dTdz * model.grid.Lz

set!(model, u=uᵢ, w=wᵢ, T=Tᵢ, S=35)

[jagoosw]

Oh sorry I missed this before, you're not passing the T field you made and put in the particles, it should be something like:

model = NonhydrostaticModel(; grid,
                              velocities = (; u, v, w),
                              tracers = (; T, S),
                              particles = lagrangian_particles)

so that the T and S fields you made are passed to the model.

[sebasti4nmc]

Ah I see, that fixed it! Thanks again for your help