Fields interpolation on lagrangian particle trajectories

[weigwu]

Hi, I am using the LagrangianParticles() function in Oceananigans. I am wondering if there is a way to interpolate the fields (e.g. u, v, w, b) onto particle trajectories and save the fields for each particle.

Many thanks,
Roger

[jagoosw]

Hi Roger,

You could do this by making a custom particle, using tracked fields to set the particle properties, and saving the particles the normal way. For example:

using Oceananigans
using Oceananigans.Fields: VelocityFields
using StructArrays

struct TrackingParticle{FT}
    x :: FT
    y :: FT
    z :: FT
    u :: FT
    v :: FT
    w :: FT
    b :: FT
end

grid = RectilinearGrid(; size = (128, 128, 128), extent = (100, 100, 100))

n = 100

x₀, y₀, z₀, u₀, v₀, w₀, b₀ = [zeros(n) for i in 1:7]

particles = StructArray{TrackingParticle}((x₀, y₀, z₀, u₀, v₀, w₀, b₀))

u, v, w = VelocityFields(grid)

b = CenterField(grid)

lagrangian_particles = LagrangianParticles(particles, tracked_fields = (; u, v, w, b))

model = NonhydrostaticModel(; grid, velocities = (; u, v, w), particles = lagrangian_particles, buoyancy = BuoyancyTracer(), tracers = (; b))

simulation = Simulation(model, Δt = 1, stop_iteration = 100)

simulation.output_writers[:particles] = JLD2OutputWriter(model, (particles=model.particles,), filename="particles", schedule=IterationInterval(1))

[glwagner]

Also worth noting that, as far as I can tell, interpolate works with AbstractOperations and so it should be possible to use them as tracked fields (ie the tracked fields do not have to be precomputed).

[weigwu]

Hi Jago and Gregory, thanks a lot for providing an example. Do you have any idea in tracking a field involving the derivative of velocities (e.g., vorticity)? Can I still make a custom particle as the way Jago showed?

[jagoosw]

Are tracked_fields used in a kernel? Maybe it can be a Dict (which can be appended to after a model is created) rather than NamedTuple.

Would definitely also be useful to have a way to auto add prognostic and auxiliary model fields!

I don't think they're used in kernels. I don't know why we would need particles to have their own fields that can't just be part of the model. Wouldn't adding in that way require a materialisation step in the model setup also?

[jagoosw]

Hi Jago and Gregory, thanks a lot for providing an example. Do you have any idea in tracking a field involving the derivative of velocities (e.g., vorticity)? Can I still make a custom particle as the way Jago showed?

You can make an abstract operation and track that, e.g. ζ = ∂x(v) - ∂y(u), then add ζ as a property in the struct and add to tracked_fields.

[glwagner]

Note the tracked fields can be any abstract array, so fields, abstract operations, and any other array-like object that supports the 3D index i, j, k will work.

[jagoosw]

Are tracked_fields used in a kernel? Maybe it can be a Dict (which can be appended to after a model is created) rather than NamedTuple.

Would definitely also be useful to have a way to auto add prognostic and auxiliary model fields!

I don't think they're used in kernels. I don't know why we would need particles to have their own fields that can't just be part of the model. Wouldn't adding in that way require a materialisation step in the model setup also?

[glwagner]

I remember we had this discussion before also about auxiliary_fields --- those would also benefit from being Dict. One argument against is that model cannot be concretely typed. I just don't know if that kind of change would have unforseen negative future consequences --- for example, for automatic differentiation of time_step!(model, dt).

[glwagner]

I wonder --- should we fix how particles are implemented so that they can be built outside the model? What is the advantage of having them be a property of model? It seems like their dynamics could be handled via callbacks. This would add flexibility and might simplify things.

[jagoosw]

Hi Jago and Gregory, thanks a lot for providing an example. Do you have any idea in tracking a field involving the derivative of velocities (e.g., vorticity)? Can I still make a custom particle as the way Jago showed?

You can make an abstract operation and track that, e.g. ζ = ∂x(v) - ∂y(u), then add ζ as a property in the struct and add to tracked_fields.