Add/remove Lagrangian particles #3884

alesok · 2024-10-30T08:32:26Z

alesok
Oct 30, 2024

Hello,

I am trying to model methane bubbles using Lagrangian particles and need to remove particles when they reach the surface. When I remove a particle using

            bubble = LazyRow(particles.properties, n)
            bubble.z += wₚ * Δt
            if bubble.z > -0.01
                StructArrays.foreachfield(v -> deleteat!(v, n), particles.properties)
            end

I got and DivideError (see below)

Is there a proper way to delete particles on the surface and release new ones on the bottom during simulation?
One way around is to move particle from the surface to the bottom, but it is not a good approach. Ideally I want to run model with variable number of particles in the domain.

Thank you.
Alexander

ERROR: LoadError: DivideError: integer division error
Stacktrace:
  [1] div
    @ ./int.jl:295 [inlined]
  [2] div
    @ ./div.jl:310 [inlined]
  [3] div
    @ ./div.jl:355 [inlined]
  [4] fld
    @ ./div.jl:321 [inlined]
  [5] mod
    @ ./int.jl:287 [inlined]
  [6] #4
    @ ~/.julia/packages/KernelAbstractions/491pi/src/nditeration.jl:121 [inlined]
  [7] ntuple
    @ ./ntuple.jl:48 [inlined]
  [8] partition(ndrange::Tuple{Int64}, __workgroupsize::Tuple{Int64})
    @ KernelAbstractions.NDIteration ~/.julia/packages/KernelAbstractions/491pi/src/nditeration.jl:119
  [9] partition
    @ ~/.julia/packages/KernelAbstractions/491pi/src/KernelAbstractions.jl:640 [inlined]
 [10] launch_config
    @ ~/.julia/packages/KernelAbstractions/491pi/src/cpu.jl:87 [inlined]
 [11] (::KernelAbstractions.Kernel{…})(::StructVector{…}, ::Vararg{…}; ndrange::Nothing, workgroupsize::Nothing)
    @ KernelAbstractions ~/.julia/packages/KernelAbstractions/491pi/src/cpu.jl:40
 [12] (::KernelAbstractions.Kernel{…})(::StructVector{…}, ::Vararg{…})
    @ KernelAbstractions ~/.julia/packages/KernelAbstractions/491pi/src/cpu.jl:39
 [13] advect_lagrangian_particles!(particles::LagrangianParticles{…}, model::NonhydrostaticModel{…}, Δt::Float64)
    @ Oceananigans.Models.LagrangianParticleTracking ~/.julia/packages/Oceananigans/BJj1C/src/Models/LagrangianParticleTracking/lagrangian_particle_advection.jl:193
 [14] step_lagrangian_particles!
    @ ~/.julia/packages/Oceananigans/BJj1C/src/Models/LagrangianParticleTracking/LagrangianParticleTracking.jl:143 [inlined]
 [15] step_lagrangian_particles!
    @ ~/.julia/packages/Oceananigans/BJj1C/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl:112 [inlined]
 [16] time_step!(model::NonhydrostaticModel{…}, Δt::Float64; callbacks::Tuple{})
    @ Oceananigans.TimeSteppers ~/.julia/packages/Oceananigans/BJj1C/src/TimeSteppers/runge_kutta_3.jl:131
 [17] time_step!(sim::Simulation{…})
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/BJj1C/src/Simulations/run.jl:134
 [18] run!(sim::Simulation{…}; pickup::Bool)
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/BJj1C/src/Simulations/run.jl:97
 [19] run!(sim::Simulation{…})
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/BJj1C/src/Simulations/run.jl:85
 [20] top-level scope
    @ ~/MyFiles/Oceananigans_Himmerfjärden/src/cube.jl:142
 [21] include(fname::String)
    @ Main ./sysimg.jl:38
 [22] top-level scope
    @ REPL[1]:1
in expression starting at /Users/sokolov/MyFiles/Oceananigans_Himmerfjärden/src/cube.jl:142
Some type information was truncated. Use `show(err)` to see complete types.

Answered by ali-ramadhan

Oct 30, 2024

Hi @alesok! Lagrangian particle advection should work even as you add and remove particles by modifying the particle.properties struct array as you seem to be doing.

I suspect you're seeing this error because there are zero particles left, and KernelAbstractions.jl is trying to partition the amount of work to be done but ends up dividing by zero somewhere using the div function. You can reproduce the same error easily:

julia> div(5, 0)
ERROR: DivideError: integer division error
Stacktrace:
 [1] div(x::Int64, y::Int64)
   @ Base ./int.jl:295
 [2] top-level scope
   @ REPL[3]:1

I'm not sure whether to classify this as a bug. Oceananigans could do nothing if there are zero particles. If thi…

View full answer

ali-ramadhan · 2024-10-30T13:03:55Z

ali-ramadhan
Oct 30, 2024
Maintainer

Hi @alesok! Lagrangian particle advection should work even as you add and remove particles by modifying the particle.properties struct array as you seem to be doing.

I suspect you're seeing this error because there are zero particles left, and KernelAbstractions.jl is trying to partition the amount of work to be done but ends up dividing by zero somewhere using the div function. You can reproduce the same error easily:

julia> div(5, 0)
ERROR: DivideError: integer division error
Stacktrace:
 [1] div(x::Int64, y::Int64)
   @ Base ./int.jl:295
 [2] top-level scope
   @ REPL[3]:1

I'm not sure whether to classify this as a bug. Oceananigans could do nothing if there are zero particles. If this is something you're interested in working on, pull requests are very welcome!

The quick and easy solution might be to ensure that there's always at least one particle around to be advected.

PS: Hope it's okay that I modified your post to format your code and stack trace using triple backticks (```) so it is easier to read.

8 replies

glwagner Oct 31, 2024
Maintainer

Huh true, so the question is why this kernel seems to slip through

alesok Oct 31, 2024
Author

So I can run a model with 0 particles for a while and then can start releasing particles every say 60 seconds from the bottom and need to take care not to remove the very last particle. Good for me.
Thank you

simone-silvestri Oct 31, 2024
Maintainer

Huh true, so the question is why this kernel seems to slip through

Looks like particle advection uses another syntax 😅 .

Oceananigans.jl/src/Models/LagrangianParticleTracking/lagrangian_particle_advection.jl

Lines 192 to 193 in f2a8fb3

    
           advect_particles_kernel! = _advect_particles!(device(arch), workgroup, worksize) 
        
           advect_particles_kernel!(particles.properties, particles.restitution, model.grid, Δt, total_velocities(model))

I guess using the launch! function here might clear this issue

glwagner Oct 31, 2024
Maintainer

This is probably a legacy from a time when launch! was fairly inflexible.

glwagner Oct 31, 2024
Maintainer

So I can run a model with 0 particles for a while and then can start releasing particles every say 60 seconds from the bottom and need to take care not to remove the very last particle. Good for me.
Thank you

For sure @alesok there is no question in my mind that this feature is important, and that not supporting is legitimately a bug. Users should not have to stress about how many particles they have. 0 is a number.

alesok · 2024-10-31T18:58:46Z

alesok
Oct 31, 2024
Author

On 31 Oct 2024, at 19:43, Gregory L. Wagner ***@***.***> wrote: So I can run a model with 0 particles for a while and then can start releasing particles every say 60 seconds from the bottom and need to take care not to remove the very last particle. Good for me. Thank you For sure @alesok <https://github.com/alesok> there is no question in my mind that this feature is important, and that not supporting is legitimately a bug. Users should not have to stress about how many particles they have. 0 is a number.

For some reason particles = LagrangianParticles(StructArray{Bubble}(undef,0), restitution=1.0, dynamics=bubble_dynamics!, tracked_fields=(T=tracers.T, S=tracers.S) ) Gives ERROR: DivideError: integer division error (see below) But particles = LagrangianParticles(StructArray{Bubble}(undef,1), restitution=1.0, dynamics=bubble_dynamics!, tracked_fields=(T=tracers.T, S=tracers.S) ) Works as expected (despite undef) It is OK with me, seems to be something Thank you ``` [ Info: Initializing simulation... Iteration: 0000, time: 0 seconds, Δt: 11 seconds, max(|u|, |v|, |w|) = (1.5e-05, 1.6e-05, 1.2e-05) ms⁻¹, wall time: 0 seconds [ Info: ... simulation initialization complete (6.467 seconds) [ Info: Executing initial time step... ERROR: DivideError: integer division error Stacktrace: [1] div @ ./int.jl:295 [inlined] [2] div @ ./div.jl:310 [inlined] [3] div @ ./div.jl:355 [inlined] [4] fld @ ./div.jl:321 [inlined] [5] mod @ ./int.jl:287 [inlined] [6] #4 @ ~/.julia/packages/KernelAbstractions/491pi/src/nditeration.jl:121 [inlined] [7] ntuple @ ./ntuple.jl:48 [inlined] [8] partition(ndrange::Tuple{Int64}, __workgroupsize::Tuple{Int64}) @ KernelAbstractions.NDIteration ~/.julia/packages/KernelAbstractions/491pi/src/nditeration.jl:119 [9] partition @ ~/.julia/packages/KernelAbstractions/491pi/src/KernelAbstractions.jl:640 [inlined] [10] launch_config @ ~/.julia/packages/KernelAbstractions/491pi/src/cpu.jl:87 [inlined] [11] (::KernelAbstractions.Kernel{…})(::Vector{…}, ::Vararg{…}; ndrange::Nothing, workgroupsize::Nothing) @ KernelAbstractions ~/.julia/packages/KernelAbstractions/491pi/src/cpu.jl:40 [12] (::KernelAbstractions.Kernel{…})(::Vector{…}, ::Vararg{…}) @ KernelAbstractions ~/.julia/packages/KernelAbstractions/491pi/src/cpu.jl:39 [13] update_lagrangian_particle_properties!(particles::LagrangianParticles{…}, model::NonhydrostaticModel{…}, Δt::Float64) @ Oceananigans.Models.LagrangianParticleTracking ~/.julia/packages/Oceananigans/BJj1C/src/Models/LagrangianParticleTracking/update_lagrangian_particle_properties.jl:31 [14] step_lagrangian_particles!(particles::LagrangianParticles{…}, model::NonhydrostaticModel{…}, Δt::Float64) @ Oceananigans.Models.LagrangianParticleTracking ~/.julia/packages/Oceananigans/BJj1C/src/Models/LagrangianParticleTracking/LagrangianParticleTracking.jl:137 [15] step_lagrangian_particles! @ ~/.julia/packages/Oceananigans/BJj1C/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl:112 [inlined] [16] time_step!(model::NonhydrostaticModel{…}, Δt::Float64; callbacks::Tuple{}) @ Oceananigans.TimeSteppers ~/.julia/packages/Oceananigans/BJj1C/src/TimeSteppers/runge_kutta_3.jl:115 [17] time_step!(sim::Simulation{…}) @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/BJj1C/src/Simulations/run.jl:122 [18] run!(sim::Simulation{…}; pickup::Bool) @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/BJj1C/src/Simulations/run.jl:97 [19] run!(sim::Simulation{…}) @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/BJj1C/src/Simulations/run.jl:85 [20] top-level scope @ ~/MyFiles/Oceananigans_Himmerfjärden/src/cube.jl:114 Some type information was truncated. Use `show(err)` to see complete types. ```

…

— Reply to this email directly, view it on GitHub <#3884 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AH5SLJ5CEWQ7HNI4HTZ5QZ3Z6J25DAVCNFSM6AAAAABQ3QFXU2VHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCMJRGUZDKOI>. You are receiving this because you were mentioned.

2 replies

glwagner Oct 31, 2024
Maintainer

can you please format your code with triple backticks ("```")? It's hard to read otherwise

glwagner Oct 31, 2024
Maintainer

what is length(StructArray{Bubble}(undef,0))?

alesok · 2024-10-31T23:02:07Z

alesok
Oct 31, 2024
Author

On 31 Oct 2024, at 23:17, Gregory L. Wagner ***@***.***> wrote: what is length(StructArray{Bubble}(undef,0))?

``` julia> length(StructArray{Bubble}(undef,0)) 0 ```

11 replies

glwagner Nov 1, 2024
Maintainer

Thank you! I think there were some unrelated stuff in that script, like WENO and the buble dynamics? Let me know if this is a good minimal setup:

using Oceananigans

xp = Array{Float64}(undef, 0)
yp = Array{Float64}(undef, 0)
zp = Array{Float64}(undef, 0)

grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
particles = LagrangianParticles(x=xp, y=yp, z=zp)
model = NonhydrostaticModel(; grid, particles)
time_step!(model, 1)

glwagner Nov 1, 2024
Maintainer

Indeed it gives me

ERROR: LoadError: DivideError: integer division error
Stacktrace:
  [1] div
    @ ./int.jl:295 [inlined]
  [2] div
    @ ./div.jl:308 [inlined]
  [3] div
    @ ./div.jl:353 [inlined]
  [4] fld
    @ ./div.jl:319 [inlined]
  [5] mod(x::Int64, y::Int64)
    @ Base ./int.jl:287
  [6] #4
    @ ~/.julia/packages/KernelAbstractions/iW1Rw/src/nditeration.jl:121 [inlined]
  [7] ntuple
    @ ./ntuple.jl:48 [inlined]
  [8] partition
    @ ~/.julia/packages/KernelAbstractions/iW1Rw/src/nditeration.jl:119 [inlined]
  [9] partition
    @ ~/.julia/packages/KernelAbstractions/iW1Rw/src/KernelAbstractions.jl:640 [inlined]
 [10] launch_config
    @ ~/.julia/packages/KernelAbstractions/iW1Rw/src/cpu.jl:87 [inlined]
 [11] (::KernelAbstractions.Kernel{…})(::StructArrays.StructVector{…}, ::Vararg{…}; ndrange::Nothing, workgroupsize::Nothing)
    @ KernelAbstractions ~/.julia/packages/KernelAbstractions/iW1Rw/src/cpu.jl:40
 [12] (::KernelAbstractions.Kernel{…})(::StructArrays.StructVector{…}, ::Vararg{…})
    @ KernelAbstractions ~/.julia/packages/KernelAbstractions/iW1Rw/src/cpu.jl:39
 [13] advect_lagrangian_particles!(particles::LagrangianParticles{…}, model::NonhydrostaticModel{…}, Δt::Float64)
    @ Oceananigans.Models.LagrangianParticleTracking ~/Projects/Oceananigans.jl/src/Models/LagrangianParticleTracking/lagrangian_particle_advection.jl:193
 [14] step_lagrangian_particles!
    @ ~/Projects/Oceananigans.jl/src/Models/LagrangianParticleTracking/LagrangianParticleTracking.jl:143 [inlined]
 [15] step_lagrangian_particles!
    @ ~/Projects/Oceananigans.jl/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl:102 [inlined]
 [16] time_step!(model::NonhydrostaticModel{…}, Δt::Int64; callbacks::Vector{…})
    @ Oceananigans.TimeSteppers ~/Projects/Oceananigans.jl/src/TimeSteppers/runge_kutta_3.jl:115
 [17] time_step!(model::NonhydrostaticModel{…}, Δt::Int64)
    @ Oceananigans.TimeSteppers ~/Projects/Oceananigans.jl/src/TimeSteppers/runge_kutta_3.jl:81
 [18] top-level scope
    @ ~/Projects/Oceananigans.jl/mwe2.jl:10
 [19] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [20] top-level scope
    @ REPL[1]:1
in expression starting at /Users/gregorywagner/Projects/Oceananigans.jl/mwe2.jl:10
Some type information was truncated. Use `show(err)` to see complete types.

alesok Nov 1, 2024
Author

Yes, it is nice minimal setup to show problem with zero particles from the beginning run (initialisation with zero particles). I guess it is a bit other story, when zero particles situation occurred during model run. That is why particles dynamics is there.

glwagner Nov 1, 2024
Maintainer

ok, i used the above mwe to implement a test and i think things may work now: #3891

@alesok please let me know if the mwe covers your case (eg the bubble dynamics are not relevant to the error). You could also test the code on #3891

glwagner Nov 1, 2024
Maintainer

That is why particles dynamics is there.

Ah right I see. I'm not sure we can do anything on the Oceananigans side about that. But I will think about it. We do want to run dynamics even if there are no particles (because we may want to create a particle?) There might be some ways to make things a little smoother though

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add/remove Lagrangian particles #3884

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 3 comments 21 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Add/remove Lagrangian particles #3884

alesok Oct 30, 2024

Replies: 3 comments · 21 replies

ali-ramadhan Oct 30, 2024 Maintainer

glwagner Oct 31, 2024 Maintainer

alesok Oct 31, 2024 Author

simone-silvestri Oct 31, 2024 Maintainer

glwagner Oct 31, 2024 Maintainer

glwagner Oct 31, 2024 Maintainer

alesok Oct 31, 2024 Author

glwagner Oct 31, 2024 Maintainer

glwagner Oct 31, 2024 Maintainer

alesok Oct 31, 2024 Author

glwagner Nov 1, 2024 Maintainer

glwagner Nov 1, 2024 Maintainer

alesok Nov 1, 2024 Author

glwagner Nov 1, 2024 Maintainer

glwagner Nov 1, 2024 Maintainer

alesok
Oct 30, 2024

Replies: 3 comments 21 replies

ali-ramadhan
Oct 30, 2024
Maintainer

glwagner Oct 31, 2024
Maintainer

alesok Oct 31, 2024
Author

simone-silvestri Oct 31, 2024
Maintainer

glwagner Oct 31, 2024
Maintainer

glwagner Oct 31, 2024
Maintainer

alesok
Oct 31, 2024
Author

glwagner Oct 31, 2024
Maintainer

glwagner Oct 31, 2024
Maintainer

alesok
Oct 31, 2024
Author

glwagner Nov 1, 2024
Maintainer

glwagner Nov 1, 2024
Maintainer

alesok Nov 1, 2024
Author

glwagner Nov 1, 2024
Maintainer

glwagner Nov 1, 2024
Maintainer