Add custom forcing velocity to particles

[simone-silvestri]

With @xkykai we were looking at ways to add a sinking velocity to particles to be added to the advecting velocity.

This PR introduces a framework to do this by adding a particle_forcing_u function.

The usage will be something like this:

struct SinkingParticle{T} <: AbstractParticle
    x :: T
    y :: T
    z :: T
    sinking_velocity :: T
end

import Oceananigans: particle_advective_forcing_w

@inline particle_advective_forcing_w(particles::SinkingParticle, p) = particles.sinking_velocity[p]

[glwagner]

Wasn't this already possible?

[glwagner]

This is not a horrible interface, except for the name "particle_advective_forcing_w"...

But this deviates from the other interfaces we provide.

It might make more sense to build an interface based on LagrangianParticles. You may want

struct ParticleAdvectionForcing
    x
    y
    z
    parameters
end

then a property called forcing to LagrangianParticles.

And a simple way to specify particle sinking.

Since the forcing should depend on particle, users also have the option to dispatch on it.

I wasn't sure what p was in the above. We just need forcing(particle) or forcing(particle, parameters) right? The particle contains all relevant information.

[simone-silvestri]

right, particle is a list of particles though. p is the index that we have to pass in. we can call it forcing but it is specifically an advective forcing (or velocity forcing) because it is added to the velocity

[jagoosw]

But this deviates from the other interfaces we provide.

I think we kind of have a difference in the interface between particles and body forces. For tracers etc we always launch a kernel over every cell, get the arguments for the forcing at that point, and pass that to forcing functions. But for particles, the only way to add forcings (without making custom structs) is through the custom_dynamics where we just give the users forcing function all of the particles together and leave them to e.g. launch a kernel over all of them.

It might be helpful to add this kind of forcing more similar to the first interface and ask users to define slip velocity functions with arbitrary arguments e.g. w(x, y, z, t, $\rho$, d), and then we make the inside of the lagrangian particle advection kernel pass the correct particle properties to the function.

I don't think it would be as simple to generalise particle forcings because we don't know what users want to do with the properties and in the current setup couldn't assume that they want every particle property integrated from some tendency function.

[simone-silvestri]

The problem with custom forcing is that the dynamics function is not implemented where the particle is actually advected (inside the particle advection kernel). We need a way to calculate a particle velocity that differs from the fluid velocity

for example, particles with inertia that follow $$\frac{d \boldsymbol{v}}{dt} = \frac{1}{\tau_p}(\boldsymbol{u} - \boldsymbol{v})$$ where $\boldsymbol{v}$ is the particle velocity and $\boldsymbol{u}$ is the fluid velocity can be implemented as:

struct InertialParticle{T} <: AbstractParticle
    x :: T
    y :: T
    z :: T
    u :: T
    v :: T
    w :: T
    particle_respose_time :: T
end

import Oceananigans: particle_velocity_u
import Oceananigans: particle_velocity_v
import Oceananigans: particle_velocity_w

@inline particle_velocity_u(u_fluid, particles::InertialParticle, p, Δt) = particles.u[p] + Δt / particles.particle_response_time[p] * (u_fluid -  particles.u[p])
@inline particle_velocity_v(v_fluid, particles::InertialParticle, p, Δt) = particles.v[p] + Δt / particles.particle_response_time[p] * (v_fluid -  particles.v[p])
@inline particle_velocity_w(w_fluid, particles::InertialParticle, p, Δt) = particles.w[p] + Δt / particles.particle_response_time[p] * (w_fluid -  particles.w[p])

[glwagner]

right, particle is a list of particles though. p is the index that we have to pass in. we can call it forcing but it is specifically an advective forcing (or velocity forcing) because it is added to the velocity

That's a detail and can be implemented regardless of the type of interface.

The difference is

1. Users extend an Oceananigans function to implement a forcing, OR
2. Users pass a function into ParticleAdvectiveForcing to implement a forcing

Option 1 is fine, no doubt. It's just that the rest of our user interfaces implements option 2. I think we should have a uniform interface across the code. Our original motivation for choosing option 2 is because we believed it would be easier to use for people who don't know Julia very well. People should be able to use Oceananigans without knowing what "multiple dispatch" is.

We can revisit this. If we think that extending forcing function is a better pattern, we should change it throughout the code.

[glwagner]

I think we kind of have a difference in the interface between particles and body forces. For tracers etc we always launch a kernel over every cell, get the arguments for the forcing at that point, and pass that to forcing functions. But for particles, the only way to add forcings (without making custom structs) is through the custom_dynamics where we just give the users forcing function all of the particles together and leave them to e.g. launch a kernel over all of them.

Agree, but what point are you trying to make? I don't understand.

It might be helpful to add this kind of forcing more similar to the first interface and ask users to define slip velocity functions with arbitrary arguments e.g. w(x, y, z, t, $\rho$, d), and then we make the inside of the lagrangian particle advection kernel pass the correct particle properties to the function.

This is basically what I was suggesting.

I don't think it would be as simple to generalise particle forcings because we don't know what users want to do with the properties and in the current setup couldn't assume that they want every particle property integrated from some tendency function.

I don't understand what you're trying to say. Can you elaborate or give an example?

[glwagner]

The problem with custom forcing is that the dynamics function is not implemented where the particle is actually advected (inside the particle advection kernel). We need a way to calculate a particle velocity that differs from the fluid velocity

for example, particles with inertia that follow dvdt=1τp(u−v) where v is the particle velocity and u is the fluid velocity can be implemented as:

First of all, I think it's great to implement a better interface for forcing particles. So I'm on board with this.

You could probably implement inertia (and many other things) by storing the previous particle velocity (or whatever properties you needed for correct dynamics) inside the particle struct. So I'm pretty sure you could also implement these dynamics with the current interface.

As I said though, I agree it's a good idea to make this kind of thing easier.

[glwagner]

The problem with custom forcing is that the dynamics function is not implemented where the particle is actually advected (inside the particle advection kernel). We need a way to calculate a particle velocity that differs from the fluid velocity for example, particles with inertia that follow dvdt=1τp(u−v) where v is the particle velocity and u is the fluid velocity can be implemented as:

First of all, I think it's great to implement a better interface for forcing particles. So I'm on board with this. You could probably implement inertia (and many other things) by storing the previous particle velocity (or whatever properties you needed for correct dynamics) inside the particle struct. So I'm pretty sure you could also implement these dynamics with the current interface.

As I said though, I agree it's a good idea to make this kind of thing easier.

Actually I realized that dynamics is called before advection, which makes it more annoying:

https://github.com/CliMA/Oceananigans.jl/blob/dd893745e50fd7aef945f765df929d2a7510db6e/src/Models/LagrangianParticleTracking/LagrangianParticleTracking.jl#L123-L132

That makes it a little trickier. It's still possible, since you can compute a "fake particle update", and then add the inertia term, before the "real update" that occurs in advect_lagrangian_particles!. But it's complicated.

Why do we call dynamics before advection? Should we switch the order of those two? It also seems like we should update particle properties after we advect the particles.

[simone-silvestri]

We could think about a ParticleAdvectiveForcing or also a ParticleVelocityTendency function that specifies the rhs of $$\frac{d \boldsymbol{v}}{dt} = \boldsymbol{G}$$

[simone-silvestri]

What is the use of the dynamics function?

ok I see, if we want to remove or add particles or other non-advective phenomena. I think, yes it should be after advection

[xkykai]

We could think about a ParticleAdvectiveForcing or also a ParticleVelocityTendency function that specifies the rhs of dvdt=G

One example that this might be tricky for this approach is where the particles have a constant sinking velocity on top of the fluid velocity.

But I suppose if dynamics is run after advection this can be done through dynamics, by calling another advect_lagrangian_particles! in dynamics with the relevant velocities.

[glwagner]

What is the use of the dynamics function?

ok I see, if we want to remove or add particles or other non-advective phenomena. I think, yes it should be after advection

It's a catch all that can be used to implement a wide variety of physics. Basically a callback.

[glwagner]

We could think about a ParticleAdvectiveForcing or also a ParticleVelocityTendency function that specifies the rhs of

What is the difference between this and what we have been discussing?

[glwagner]

We could think about a ParticleAdvectiveForcing or also a ParticleVelocityTendency function that specifies the rhs of dvdt=G

One example that this might be tricky for this approach is where the particles have a constant sinking velocity on top of the fluid velocity.

A constant sinking velocity is easy to implement using dynamics, since it's constant everywhere and doesn't depend on the position / fluid velocity experienced by the particle.

[glwagner]

There used to be an example of using dynamics to implement particles that played "rock, paper, scissors"; ie near-neighbors were identified, and then made to interact in some way that changed particle properties.

There is scope for building out a user interface for more generic particle interactions of course. I'm not convinced that this is so important that we should spoon-feed users with an interface though. In principle, dynamics can do everything and advanced users can write their own code.

It depends on how important the uses get, and how clever we can be with an interface (so it is easy to maintain and doesn't screw things up for future development)

[xkykai]

We could think about a ParticleAdvectiveForcing or also a ParticleVelocityTendency function that specifies the rhs of dvdt=G

One example that this might be tricky for this approach is where the particles have a constant sinking velocity on top of the fluid velocity.

A constant sinking velocity is easy to implement using dynamics, since it's constant everywhere and doesn't depend on the position / fluid velocity experienced by the particle.

I don't think a constant sinking velocity can be implemented in dynamics in this setup where dynamics comes before advection. You could imagine a particle being advected with its sinking velocity at the dynamics step, then later on at the advect_lagrangian_particles! it'll be advected by the field velocity at a new location, which will be of a different value. I suppose this is fine if dynamics and advect_lagrangian_particles! are reversed.

However, I think that advection forcing should be run at the same time as advect_lagrangian_particles!. This would be important for cases where the additional advection is a function of its location for example.

[glwagner]

We could think about a ParticleAdvectiveForcing or also a ParticleVelocityTendency function that specifies the rhs of dvdt=G

One example that this might be tricky for this approach is where the particles have a constant sinking velocity on top of the fluid velocity.

A constant sinking velocity is easy to implement using dynamics, since it's constant everywhere and doesn't depend on the position / fluid velocity experienced by the particle.

I don't think a constant sinking velocity can be implemented in dynamics in this setup where dynamics comes before advection. You could imagine a particle being advected with its sinking velocity at the dynamics step, then later on at the advect_lagrangian_particles! it'll be advected by the field velocity at a new location, which will be of a different value. I suppose this is fine if dynamics and advect_lagrangian_particles! are reversed.

For a constant velocity, the only difference that switching the order of function calls makes is the way we interpret the initial condition or final state. For example, if we call dynamics after advection, but also shift the initial condition vertical positions by dt * w_sinking, then the result would be identical as calling dynamics before advection.

[glwagner]

However, I think that advection forcing should be run at the same time as advect_lagrangian_particles!. This would be important for cases where the additional advection is a function of its location for example.

Agreed

[xkykai]

I added a ParticleAdvectionForcing that behaves similarly to dynamics(particles, model, Δt) which allows user to specify any velocities it wishes to add to the model by passing in NonHydrostaticModel(..., advective_forcing = ParticleAdvectiveForcing(u=some_u, v=some_v, w=some_w)) . The arguments are functions that take (particles, model, Δt) as inputs.

When I am trying to test my code I ran into a problem:

ERROR: UndefVarError: `flattened_node` not defined

https://github.com/CliMA/Oceananigans.jl/blob/a73e845a7a7bce8e22e9453670c1c20b67dbfc3a/src/Models/LagrangianParticleTracking/lagrangian_particle_advection.jl#L92

@simone-silvestri I see that you added this and I couldn't find anywhere in the code where flattened_node is defined. What is it supposed to be?

[simone-silvestri]

I added a ParticleAdvectionForcing that behaves similarly to dynamics(particles, model, Δt) which allows user to specify any velocities it wishes to add to the model by passing in NonHydrostaticModel(..., advective_forcing = ParticleAdvectiveForcing(u=some_u, v=some_v, w=some_w)) . The arguments are functions that take (particles, model, Δt) as inputs.

When I am trying to test my code I ran into a problem:

ERROR: UndefVarError: `flattened_node` not defined

https://github.com/CliMA/Oceananigans.jl/blob/a73e845a7a7bce8e22e9453670c1c20b67dbfc3a/src/Models/LagrangianParticleTracking/lagrangian_particle_advection.jl#L92

@simone-silvestri I see that you added this and I couldn't find anywhere in the code where flattened_node is defined. What is it supposed to be?

This is from #3395. We should wait for that PR to be merged before merging this one

[simone-silvestri]

to test the code you can merge #3395 in this PR

[simone-silvestri]

I am still not super convinced about the ParticleForcing implementation. This would not work well in case we want to solve the tendency for the particle-velocity evolution equation.

For example, there is no non-convoluted way to implement a particle with drag where $\frac{d \boldsymbol{v}}{dt} = \frac{C_d}{\tau}(\boldsymbol{u} - \boldsymbol{v})$ with an "advecting forcing" (we would have to remove the fluid velocity to not count it twice).

In my opinion, I would propose to call it ParticleVelocities with fields u_velocity, v_velocity, and w_velocity and parameters which are functions that default to the fluid velocity if not specified (where fluid velocity is passed to the function) with more complicated implementations possible such as in the above case with drag

the actual implementation is very similar it's just a matter of including the fluid velocity in the advecting_velocity. So instead of

u = particle_u_velocity(i, j, k, u_fluid, advective_forcing.u, particles, p, grid, clock, Δt, model_fields, advective_forcing.u.parameters)

where

@inline particle_u_velocity(i, j, k, u_fluid, forcing::ParticleDiscreteForcing, particles, p, grid, clock, Δt, model_fields, parameters) = u_fluid + forcing(i, j, k, particles, p, grid, clock, Δt, model_fields, parameters) 

just write

u = advective_velocity.u(i, j, k, u_fluid, particles, p, grid, clock, Δt, model_fields, parameters)

where advective_velocity is a ParticleVelocities and by default (this is more of pseudocode but it gives an idea)

ParticleVelocities(; u = (i, j, k, u_fluid, args...) -> u_fluid, v = (i, j, k, v_fluid, args...) -> v_fluid, w = (i, j, k, w_fluid, args...) -> w_fluid)...

[glwagner]

What does i, j, k refer to?

[simone-silvestri]

in this case they would be the signature of the function as in the particle_advecting_forcing implementation

[glwagner]

But particles don't have an index i, j, k. Are these fractional indices?

[simone-silvestri]

these are if we want to make the particle velocity depend on other eulerian fields (like forcing for our prognostic variables). An example can be the buoyancy field for buoyant particles. If we decide that that role can be undertaken by the dynamics kernel, we can remove the i, j, k from the signature.

We should probably come up with a couple of examples to include in this PR that show the implementation of different particle dynamics. I would include:

• buoyant particles (sinking or rising with a density-dependent vertical velocity)
• particles with drag

[glwagner]

these are if we want to make the particle velocity depend on other eulerian fields (like forcing for our prognostic variables). An example can be the buoyancy field for buoyant particles. If we decide that that role can be undertaken by the dynamics kernel, we can remove the i, j, k from the signature.

We should probably come up with a couple of examples to include in this PR that show the implementation of different particle dynamics. I would include:

• buoyant particles (sinking or rising with a density-dependent vertical velocity)
• particles with drag

The discussion is fine but we first need to resolve something more fundamental: the particles do not have indices i, j, k. The only way to use Eulerian fields to determine particle forcing is to interpolate the Eulerian fields to the location of the particle. Obviously, that's how advecting particles works in the first place, by interpolating the velocity field to the location of the particle.

[simone-silvestri]

right, sorry I confused i, j, k with particle position, I meant we have to pass the particle position x, y, z

[glwagner]

right, sorry I confused i, j, k with particle position, I meant we have to pass the particle position x, y, z

Ah ok, that makes more sense!

The particle position is available via

x = particles[p].x

etc right?

We should implement a helper function,

position(p) = (p.x, p.y, p.z)

Note that particles are also annoying on Flat grids. One needs to beware the particle position in the Flat direction, but they shouldn't have to. The same issue will affect interpolate once we merge #3395

EDIT: even better would simply be to extend interpolate to work with the particle as an argument, eg something like

up = interpolate(particles[p], u, u_loc, grid)

I think ultimately it would be nice to push this under the hood with some nice user interface for ParticleVelocity or whatever