add a :t- property to Lagrangian particles?

[ammtD]

Hello,

in relation to a discussion I started earlier, about forcing a tracer at a specific time in the model (see here, I wonder if one could add particles to a model at a specific time, and how to code that?
I am using the LagrangianParticles() option for 4800 particles with defined properties (:x, :y, :z) only.

[glwagner]

The particles are StructArray:

Oceananigans.jl/src/LagrangianParticleTracking/LagrangianParticleTracking.jl

Line 49 in 3322f18

particles = StructArray{Particle}((x, y, z))

so you should be able to add particles at any time. You can run a model for some time:

simulation.stop_time = 1
run!(simulation)

and then continue:

# add particles
simulation.stop_time += 1
run!(simulation)

You can also use Callback to add particles, I think.

I'm not immediately familiar with the syntax to do this because I don't use Lagrangian particles --- it'd be great if you post what you figure out to this discussion!

[glwagner]

The code I sent uses Particle as the particle type, but I think you can use any particle that's compatible with the particle advection scheme, etc.

It lookjs like you add to StructArray by pushing to the individual vectors:

julia> using StructArrays

julia> struct TestParticle{T}
       x :: T
       y :: T
       z :: T
       end

julia> a = StructArray{TestParticle}(([1.0], [2.0], [3.0]))
1-element StructArray(::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}) with eltype TestParticle:
 TestParticle{Float64}(1.0, 2.0, 3.0)

julia> push!(a.x, 1.1)
2-element Vector{Float64}:
 1.0
 1.1

julia> push!(a.y, 2.1)
2-element Vector{Float64}:
 2.0
 2.1

julia> push!(a.z, 3.1)
2-element Vector{Float64}:
 3.0
 3.1

julia> a
2-element StructArray(::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}) with eltype TestParticle:
 TestParticle{Float64}(1.0, 2.0, 3.0)
 TestParticle{Float64}(1.1, 2.1, 3.1)

[ammtD]

Thanks Greg,

I am going to give it a try. will have some learning to do as I haven't yet stopped and resumed a model before...

simulation.stop_time +=1 
run!(simulation) 

doesn't seem to work; the stop_time has not been updated

[glwagner]

Maybe give more information? This works for me (using the example from the README):

julia> using Oceananigans

julia> grid = RectilinearGrid(CPU(), size=(128, 128), x=(0, 2π), y=(0, 2π), topology=(Periodic, Periodic, Flat))
128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── Periodic x ∈ [-7.51279e-18, 6.28319) regularly spaced with Δx=0.0490874
├── Periodic y ∈ [-7.51279e-18, 6.28319) regularly spaced with Δy=0.0490874
└── Flat z

julia> model = NonhydrostaticModel(; grid, advection=WENO5())
┌ Warning: defaulting to uniform WENO scheme with Float64 precision, use WENO5(grid = grid) if this was not intended
└ @ Oceananigans.Advection ~/Projects/Oceananigans.jl/src/Advection/weno_fifth_order.jl:154
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
└── coriolis: Nothing

julia> ϵ(x, y, z) = 2rand() - 1
ϵ (generic function with 1 method)

julia> set!(model, u=ϵ, v=ϵ)

julia> simulation = Simulation(model; Δt=0.01, stop_time=4)
Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 10 ms
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN years
├── Stop time: 4 seconds
├── Stop iteration : Inf
├── Wall time limit: Inf
├── Callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entries

julia> simulation.stop_time
4.0

julia> run!(simulation)
[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (986.500 μs)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (10.564 ms).
[ Info: Simulation is stopping. Model time 4 seconds has hit or exceeded simulation stop time 4 seconds.

julia> simulation.stop_time
4.0

julia> simulation.stop_time += 4
8.0

julia> simulation.stop_time
8.0

julia> run!(simulation)
[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (837.166 μs)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (9.822 ms).
[ Info: Simulation is stopping. Model time 8 seconds has hit or exceeded simulation stop time 8 seconds.

julia> simulation.stop_time
8.0

julia> time(simulation)
8.0

[ammtD]

Thanks. Yes that works for me too. For the model were I have particles from the start I was using
simple_particles = LagrangianParticles(x=xp, y=yp, z=zp)

model = NonhydrostaticModel(; grid,
                            advection = WENO5(),
                            timestepper = :RungeKutta3,
                            particles = simple_particles,
                            closure = AnisotropicMinimumDissipation(),
                            ...) 

When I want to add the particles later I remove the line particles = ... from the model and run the simulation until the first stop-time.
Then I define a structure array a as in your earlier answer and use it in delayed_particles = LagrangianParticles(a).
And then I redefine model = (grid, ..., particles = delayed_particles)
adjust the simulation.stop_time, and run the simulation.

this now actually works! I don't know which step I skipped or made a mistake in earlier... thanks again.

[glwagner]

What you describe (creating a fake model from a named tuple?) might work right now but is very likely to break in the future so I wouldn't recommend continuing to use it. Instead, I think you should modify model.particles at the same time that you change simulation.stop_time. This is also probably easier / less code.

I think this recipe will achieve what you want:

1. Always create your model with the particles in it.
2. Run the model forward for t = t_first (eg Simulation(model, stop_time=t_first, ...))
3. Reset the particle positions to 0, or to whatever distribution you want. Ignore what happened to the particles before that. I'm not sure exactly what syntax this requires. Something like model.particles.x .= 0 ? (and the same for y and z). Unless if you have a custom particles type, then the syntax is whatever you chose.
4. Change simulation.stop_time.
5. Add whatever output writers you need for the next phase. If you had spin up writers, remove them.
6. Consider also resetting the model clock and iteration to 0: simulation.model.clock=0 simulation.model.time=0`.

Step 3 is the key one. Because you are completely free to modify your particle properties whenever you like (I think you can also add and remove particles at will), there's no need to build a model without particles.

It'd awesome to see some results here too! Sounds like cool work!