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]

[xkykai]

I've updated the code to ParticleVelocities and changed the function signature to

@inline particle_u_velocity(x, y, z, u_fluid, particles, p, advective_velocity::ParticleVelocities, grid, clock, Δt, model_fields) = advective_velocity.u(x, y, z, u_fluid, particles, p, grid, clock, Δt, model_fields)

I also swapped around advect_lagrangian_particles! and dynamics, so now the particles are advected before the dynamics function is applied. I found this to be easier when implementing an example with the drag.

I've implemented a simple example with a drag in the form of $\frac{d \boldsymbol{v}}{dt} = \frac{C_d}{\tau}(\boldsymbol{u} - \boldsymbol{v})$.

[xkykai]

Upon thinking about it I think perhaps it makes more sense to run dynamics first before advect_lagrangian_particles!, one example of which might be that particles sink depending on the radius of the particle, which changes with time. In such a case it is perhaps better to evolve the particle radius, then compute the sinking velocity given the new radius before advective it.

Also provided a draft example of how one could set up a problem where the particle sinks with a drag in the form of $\frac{d \boldsymbol{v}}{dt} = \frac{C_d}{\tau}(\boldsymbol{u} - \boldsymbol{v})$.
Note: in the calculation the velocity of the particle itself needs to be tracked. This is done in u_particle, v_particle, and w_particle in particles.properties where particles::LagrangianParticles. The particle velocities are computed and updated in the dynamics step, then ParticleVelocities only has functions that access the particle properties to grab the particle velocity. It is slightly clunky but unless we keep track of the particle velocities right out of the box and update them during the advect_lagrangian_particles! step, this is the way I could think of.

Since particle velocities are not updated when the particle is bounced, it will not work if the particles bounce from the boundaries back into the interior during the advection step, but for doubly-periodic domian and sinking particles it might not be very important.

[glwagner]

In such a case it is perhaps better to evolve the particle radius, then compute the sinking velocity given the new radius before advective it.

Are you sure? This could mean that the advecting velocity is extracted at time-step n but we are using a radius from n+1.

[xkykai]

In such a case it is perhaps better to evolve the particle radius, then compute the sinking velocity given the new radius before advective it.

Are you sure? This could mean that the advecting velocity is extracted at time-step n but we are using a radius from n+1.

I think that depends on how one defines the advective velocity eg. if it is defined as an explicit function of time then it is better that dynamics act first to compute the radius at time $t$ before it is being advected to time $t + \Delta t$.

[glwagner]

In such a case it is perhaps better to evolve the particle radius, then compute the sinking velocity given the new radius before advective it.

Are you sure? This could mean that the advecting velocity is extracted at time-step n but we are using a radius from n+1.

I think that depends on how one defines the advective velocity eg. if it is defined as an explicit function of time then it is better that dynamics act first to compute the radius at time t before it is being advected to time t+Δt.

But in that case, when a time-step is complete, the radius has been computed at time t but the particle position is updated to time t + Δt. This means that if the user asks for output, the two are inconsistent with one another.

[glwagner]

Actually, I think we need to change the time-stepping loop:

Oceananigans.jl/src/TimeSteppers/quasi_adams_bashforth_2.jl

Lines 93 to 100 in 3e26503

	ab2_step!(model, Δt, χ) # full step for tracers, fractional step for velocities.
	calculate_pressure_correction!(model, Δt)
	@apply_regionally correct_velocities_and_store_tendecies!(model, Δt)
	tick!(model.clock, Δt)
	update_state!(model, callbacks; compute_tendencies)
	step_lagrangian_particles!(model, Δt)

We should step particles forward before stepping the velocity field?

It also seems that we need to separate particle advection from "dynamics", and compute the dynamics once the time-stepping is complete (probably it makes sense to compute dynamics within update_state!?)

[xkykai]

In such a case it is perhaps better to evolve the particle radius, then compute the sinking velocity given the new radius before advective it.

Are you sure? This could mean that the advecting velocity is extracted at time-step n but we are using a radius from n+1.

I think that depends on how one defines the advective velocity eg. if it is defined as an explicit function of time then it is better that dynamics act first to compute the radius at time t before it is being advected to time t+Δt.

But in that case, when a time-step is complete, the radius has been computed at time t but the particle position is updated to time t + Δt. This means that if the user asks for output, the two are inconsistent with one another.

Makes sense. advect_lagrangian_particles! first before dynamics it is!

[xkykai]

We should step particles forward before stepping the velocity field?

It also seems that we need to separate particle advection from "dynamics", and compute the dynamics once the time-stepping is complete (probably it makes sense to compute dynamics within update_state!?)

Agreed. @simone-silvestri what do you think?