Some Lagrangian particles like to travel Southwest on an immersed lat-lon grid + hydrostatic model #3917

ali-ramadhan · 2024-11-12T00:10:28Z

So I'm trying to advect some particles in a HydrostaticFreeSurfaceModel with an immersed LatitudeLongitudeGrid.

I noticed that some particles travel really fast in the southwest direction and get stuck at the southwestern corner. Depending on the dynamics and the immersed boundary, tons of particles move this way or very few do. The particles behave fine once you remove the immersed boundary.

I was able to reduce the issue down to a somewhat lengthy MWE. I add an immersed seamount and some particles (with a zero coefficient of restitution) and impose a surface momentum flux to get them moving. I use WENO to stabilize the simulation. It takes a few days for the issue to show up then I plot the particles at high temporal resolution so their trajectories are clear.

I'm digging into the code to see what the issue might, but if anyone has any ideas I'd appreciate it!

I suspect it has something to do with the particles bouncing off bathymetry incorrectly (well, they shouldn't be bouncing at all with restitution = 0). Will try running on the CPU and with a small coefficient of restitution to see what happens.

MWE:

using Printf
using Distributions
using JLD2
using CairoMakie

using Oceananigans
using Oceananigans.Units

using Oceananigans.Architectures: on_architecture
using Oceananigans.OutputReaders: iterations_from_file

# Grid setup

arch = GPU()

Nλ, Nφ, Nz = 100, 100, 32

Lλ = Lφ = 1
H = 100

underlying_grid = LatitudeLongitudeGrid(
    arch;
    topology = (Bounded, Bounded, Bounded),
    size = (Nλ, Nφ, Nz),
    longitude = (-Lλ, Lλ),
    latitude = (-Lφ, Lφ),
    z = (-H, 0),
    halo = (5, 5, 5)
)

height = H/2
width = Lλ / 3
mount(λ, φ) = height * exp(-λ^2 / 2width^2) * exp(-φ^2 / 2width^2)
bottom(λ, φ) = -H + mount(λ, φ)

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom))

# Boundary conditions

v_bcs = FieldBoundaryConditions(
      top = FluxBoundaryCondition(-1e-5)
)

boundary_conditions = (; v = v_bcs)

# Particle setup

Np = 100000 # lots of particles so the issue is clear

ε = 0.1
λ_particles = 0.0
φ_particles = 0.5

xs = on_architecture(arch, rand(Uniform(-0.1, 0.1),Np))
ys = on_architecture(arch, rand(Uniform(-0.5, 0.5), Np))
zs = on_architecture(arch, rand(Uniform(-25, -1), Np))

particles = LagrangianParticles(; x=xs, y=ys, z=zs, restitution = 0.0)

# Model setup

model = HydrostaticFreeSurfaceModel(;
    grid,
    boundary_conditions,
    particles,
    momentum_advection = WENOVectorInvariant(; order=5)
)

u₀(λ, φ, z) = 1e-3 * randn()
v₀(λ, φ, z) = 1e-3 * randn()
η₀(λ, φ, z) = 1e-3 * mount(λ, φ)

set!(model, u=u₀, v=v₀, η=η₀)

# Simulation setup

simulation = Simulation(model; Δt=1second, stop_time=8days)

wizard = TimeStepWizard(cfl = 0.1, min_Δt = 0.1second, max_change = 1.1)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

progress(sim) = @printf(
    "iteration: %d, time: %s, Δt: %s, U_max=(%.2e, %.2e, %.2e)\n",
    iteration(simulation),
    prettytime(time(simulation)),
    prettytime(simulation.Δt),
    maximum(abs, model.velocities.u),
    maximum(abs, model.velocities.v),
    maximum(abs, model.velocities.w)
)

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

simulation.output_writers[:particles_jld2] =
    JLD2OutputWriter(
        model,
        (; particles = model.particles);
        filename = "debug_particles",
        schedule = TimeInterval(10minutes),
        overwrite_existing = true
    )

run!(simulation)

# Plotting

file = jldopen(simulation.output_writers[:particles_jld2].filepath)
file_iterations = iterations_from_file(file)
times = [file["timeseries/t/$i"] for i in file_iterations]
particles = [file["timeseries/particles/$i"] for i in file_iterations]
close(file)

Nt = length(file_iterations)

n = Observable(1)

title = @lift @sprintf("Particles @ %s", prettytime(times[$n]))

p_x = @lift particles[$n].x
p_y = @lift particles[$n].y

fig = Figure(size = (1000, 1000))

ax = Axis(fig[1, 1]; title = "particles", limits = ((-Lλ, Lλ), (-Lφ, Lφ)))
sc = scatter!(ax, p_x, p_y, color=:black, markersize=5)

fig[0, :] = Label(fig, title, fontsize=24, tellwidth=false)

n_start = 1 # round(Int, 0.8Nt)
record(fig, "debug_particles.mp4", n_start:Nt, framerate=10) do i
    @info "Plotting frame $i/$Nt..."
    n[] = i
end

Output:

debug_particles.mp4

The text was updated successfully, but these errors were encountered:

ali-ramadhan · 2024-11-12T03:39:48Z

Confirming that this also happens on the CPU, and on the GPU with restitution = 0.01.

glwagner · 2024-11-12T03:57:26Z

This is absurd! What is the vertical position of the particles in the movie above?

ali-ramadhan · 2024-11-12T04:05:59Z

The initial depth of the particles is Uniform(-25, -1), which should be well above the seamount which reaches up to -50 m.

EDIT: Ah sorry all the distributions below are from a followup simulation I ran with Uniform(-50, 0).

Although looking at the vertical distribution of the particles over time, at t = 16.5 hours it seems that some particles have teleported to the bottom of the domain.

And this is the same distribution at t = 6.938 days

🤔

ali-ramadhan · 2024-11-12T04:18:49Z

Hmmm, the particles teleport to the bottom before t = 10 minutes (the first time snapshot), perhaps even at iteration 1. Not sure if this is related to the Southwest travel issue, but it seems quite wrong.

t = 0:

t = 10 minutes:

glwagner · 2024-11-12T12:57:58Z

I'm wondering if that the vertical location of the particles is intrinsic to the error. The reason is that it seems that there is some part of the particle "state" that is corrupted after interaction with the boundary. Before the state is corrupted, the particles behave as expected. After the state is corrupted, the particles have spurious motion. For simple particles, the location of the particle is the only state variable. Since the horizontal position is not spurious (we are plotting it), the only remaining aspect of the particle state that can be corrupted is the vertical location.

This still leaves open why the particles would travel SW for wrong vertical position.

ali-ramadhan · 2024-11-12T16:13:55Z

I updated the MWE a bit and also plotted a meridional slice to look at vertical positions.

Looking at the meridional slice animation:

It seems that at iteration 1 a whole bunch of particles that were close to the top of the seamount teleported to the bottom left corner. Maybe these particles were inside the seamount. But then they should just be stuck there?

And particles are constantly moving to the bottom southwest corner (-x, -y, -z direction). Looks to me like it's particles that are close to the immersed boundary. So maybe another strong hint that it's the collision/bouncing logic?

MWE:

using Printf
using Distributions
using JLD2
using CairoMakie

using Oceananigans
using Oceananigans.Units

using Oceananigans.Architectures: on_architecture
using Oceananigans.OutputReaders: iterations_from_file
using Oceananigans.ImmersedBoundaries: mask_immersed_field!

# Grid setup

arch = GPU()

Nλ, Nφ, Nz = 100, 100, 32

Lλ = Lφ = 1
H = 100

underlying_grid = LatitudeLongitudeGrid(
    arch;
    topology = (Bounded, Bounded, Bounded),
    size = (Nλ, Nφ, Nz),
    longitude = (-Lλ, Lλ),
    latitude = (-Lφ, Lφ),
    z = (-H, 0),
    halo = (5, 5, 5)
)

height = H/2
width = Lλ / 3
mount(λ, φ) = height * exp(-λ^2 / 2width^2) * exp(-φ^2 / 2width^2)
bottom(λ, φ) = -H + mount(λ, φ)

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom))

# Boundary conditions

v_bcs = FieldBoundaryConditions(
      top = FluxBoundaryCondition(-1e-5)
)

boundary_conditions = (; v = v_bcs)

# Particle setup

Np = 100000 # lots of particles so the issue is clear

ε = 0.1
λ_particles = 0.0
φ_particles = 0.5

xs = on_architecture(arch, rand(Uniform(-0.1, 0.1),Np))
ys = on_architecture(arch, rand(Uniform(-0.5, 0.5), Np))
zs = on_architecture(arch, rand(Uniform(-50, 0), Np))

particles = LagrangianParticles(; x=xs, y=ys, z=zs, restitution=0.0)

# Model setup

model = HydrostaticFreeSurfaceModel(;
    grid,
    boundary_conditions,
    particles,
    momentum_advection = WENOVectorInvariant(; order=5)
)

u₀(λ, φ, z) = 1e-3 * randn()
v₀(λ, φ, z) = 1e-3 * randn()
η₀(λ, φ, z) = 1e-3 * mount(λ, φ)

set!(model, u=u₀, v=v₀, η=η₀)

# Simulation setup

simulation = Simulation(model; Δt=1second, stop_time=8days)

wizard = TimeStepWizard(cfl = 0.1, min_Δt = 0.1second, max_change = 1.1)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

progress(sim) = @printf(
    "iteration: %d, time: %s, Δt: %s, U_max=(%.2e, %.2e, %.2e)\n",
    iteration(simulation),
    prettytime(time(simulation)),
    prettytime(simulation.Δt),
    maximum(abs, model.velocities.u),
    maximum(abs, model.velocities.v),
    maximum(abs, model.velocities.w)
)

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

simulation.output_writers[:particles_jld2] =
    JLD2OutputWriter(
        model,
        (; particles = model.particles);
        filename = "debug_particles",
        schedule = TimeInterval(10minutes),
        overwrite_existing = true
    )

run!(simulation)

# Plotting

file = jldopen(simulation.output_writers[:particles_jld2].filepath)
file_iterations = iterations_from_file(file)
times = [file["timeseries/t/$i"] for i in file_iterations]
particles = [file["timeseries/particles/$i"] for i in file_iterations]
close(file)

Nt = length(file_iterations)

# Surface plot

n = Observable(1)

title = @lift @sprintf("Particles @ %s", prettytime(times[$n]))

p_x = @lift particles[$n].x
p_y = @lift particles[$n].y

fig = Figure(size = (600, 600))

ax = Axis(fig[1, 1]; title = "particles", limits = ((-Lλ, Lλ), (-Lφ, Lφ)))
sc = scatter!(ax, p_x, p_y, color=:black, markersize=5)

fig[0, :] = Label(fig, title, fontsize=24, tellwidth=false)

n_start = 1 # round(Int, 0.8Nt)
record(fig, "debug_particles_surface.mp4", n_start:Nt, framerate=10) do i
    @info "Plotting frame $i/$Nt..."
    n[] = i
end

# Meridional slice

xu, _, zu = nodes(model.velocities.u)

mask_immersed_field!(model.velocities.u, grid, (Center(), Center(), Center()), NaN)

n = Observable(1)

title = @lift @sprintf("Particles @ %s", prettytime(times[$n]))

p_y = @lift particles[$n].y
p_z = @lift particles[$n].z

fig = Figure(size = (600, 600))

ax = Axis(fig[1, 1]; title = "particles", limits = ((-Lφ, Lφ), (-H, 0)))
hm = heatmap!(ax, xu, zu, interior(model.velocities.u)[:, round(Int, Nφ/2), :] |> Array, colormap=:balance, colorrange=(-1e5, 1e5), nan_color=:gray)
sc = scatter!(ax, p_y, p_z, color=:black, markersize=5)

fig[0, :] = Label(fig, title, fontsize=24, tellwidth=false)

n_start = 1 # round(Int, 0.8Nt)
record(fig, "debug_particles_meridional.mp4", n_start:Nt, framerate=10) do i
    @info "Plotting frame $i/$Nt..."
    n[] = i
end

debug_particles_surface.mp4

debug_particles_meridional.mp4

ali-ramadhan · 2024-11-12T16:23:35Z

Don't think it would help debug that much more but would be cool to visualize this in 3D...

glwagner · 2024-11-12T20:34:50Z

They aren't teleporting right? They seem to take a few iterations to get there.

ali-ramadhan · 2024-11-12T22:09:27Z

The ones between frame 1 (t = 0) and frame 2 (t = 10 minutes) seem to be teleporting there. They could just be taking <10 minutes. But yeah the rest take a few time snapshots (10-40 minutes) to get there.

jagoosw · 2024-11-12T22:22:13Z

In the meridional slice, it looks like there is an offset between the seamount and the particles that start moving. I wonder if there is some issue interpolating something to the particles, or if the particles maybe incorrectly think they're in the boundary when they're not?

ali-ramadhan · 2024-11-12T23:42:21Z

Yeah that could also be the case. Maybe an off-by-one error?

I think it could be a plotting thing with faces and centers. I'm basically plotting the seamount as it looks in Center(), Center() but maybe I should also plot the exact bathymetry...

Although the particles should all be above z = -50 m which is the highest possible point of the seamount. Although the discretized seamount could be a bit higher I suppose.

simone-silvestri · 2024-11-13T13:34:13Z

It's an interpolation issue, and not by one, but by 2 or 3 cells! I have simplified your MWE a bit to show this:

using Oceananigans
using Oceananigans.Fields: fractional_indices, truncate_fractional_indices
using Oceananigans.Grids: ξnode, ηnode, rnode

# Grid setup
const c = Center()
const f = Face()

underlying_grid = LatitudeLongitudeGrid(;
    topology = (Bounded, Bounded, Bounded),
         size = (20, 20, 32),
    longitude = (-1, 1),
     latitude = (-1, 1),
            z = (-100, 0),
         halo = (5, 5, 5)
)

(x, y, z)  = X = (-0.082, 0.034, -49.9)
fi, fj, fk = fractional_indices(X, underlying_grid, c, c, c)
i, j, k    = truncate_fractional_indices(fi, fj, fk)

@show i, j, k

xi = ξnode(i, j, k, underlying_grid, f, f, f)
yi = ηnode(i, j, k, underlying_grid, f, f, f)
zi = rnode(i, j, k, underlying_grid, f, f, f)

xi⁻ = ξnode(i-1, j, k, underlying_grid, f, f, f)
yi⁻ = ηnode(i, j-1, k, underlying_grid, f, f, f)
zi⁻ = rnode(i, j, k-1, underlying_grid, f, f, f)

xi⁺ = ξnode(i+1, j, k, underlying_grid, f, f, f)
yi⁺ = ηnode(i, j+1, k, underlying_grid, f, f, f)
zi⁺ = rnode(i, j, k+1, underlying_grid, f, f, f)

@show x,  y,  z
@show xi, yi, zi
@show xi⁻, yi⁻, zi⁻
@show xi⁺, yi⁺, zi⁺

julia> @show x,  y,  z
(x, y, z) = (-0.082, 0.034, -49.9)
(-0.082, 0.034, -49.9)

julia> @show xi, yi, zi
(xi, yi, zi) = (-0.3, -0.2, -56.25)
(-0.3, -0.2, -56.25)

julia> @show xi⁻, yi⁻, zi⁻
(xi⁻, yi⁻, zi⁻) = (-0.4, -0.3, -59.375)
(-0.4, -0.3, -59.375)

julia> @show xi⁺, yi⁺, zi⁺
(xi⁺, yi⁺, zi⁺) = (-0.2, -0.1, -53.125)
(-0.2, -0.1, -53.125)

I ll open a PR to fix this.

simone-silvestri · 2024-11-13T14:54:29Z

In branch #3923 the MWE of @ali-ramadhan (albeit with Nx = 20 and Ny = 20 and on the CPU) spits out

debug_particles_surface.mp4

debug_particles_meridional.mp4

It looks like there are some grid effects I am not sure of, I ll continue working on the PR

ali-ramadhan · 2024-11-13T15:44:46Z

Wow nice investigative work @simone-silvestri! No more buggy particles!

In your new animations the paths of some particles does look a bit weird but I wonder how much of it is the weird flow field due to the contrived exampled.

simone-silvestri · 2024-11-13T16:09:17Z

Yeah, I am not sure about it...
Can it maybe be due to the larger cells in the horizontal? What if you try on the GPU with your fine grid?

ali-ramadhan · 2024-11-13T16:17:16Z

Running at 250x250x32 on the GPU. Will post back once the animations are made!

ali-ramadhan · 2024-11-13T17:11:17Z

Here's what I'm seeing now. I'm thinking I should plot v and w on the meridional animation to see if the particle movement makes sense.

3917_surface.mp4

3917_meridional.mp4

simone-silvestri · 2024-11-13T17:24:10Z

Ok there is definitely something wrong here. I will continue investigating.

simone-silvestri · 2024-11-13T18:57:55Z

I think the last iteration works. These are the results:

debug_particles_meridional.mp4

ali-ramadhan added the bug 🐞 Even a perfect program still has bugs label Nov 12, 2024

simone-silvestri linked a pull request Nov 13, 2024 that will close this issue

Fix bug in advect_particle #3923

Open

ali-ramadhan mentioned this issue Nov 13, 2024

Ability to set the coefficient of restitution for particles OceanBioME/OceanBioME.jl#232

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Some Lagrangian particles like to travel Southwest on an immersed lat-lon grid + hydrostatic model #3917

Some Lagrangian particles like to travel Southwest on an immersed lat-lon grid + hydrostatic model #3917

ali-ramadhan commented Nov 12, 2024 •

edited

Loading

ali-ramadhan commented Nov 12, 2024

glwagner commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024 •

edited

Loading

ali-ramadhan commented Nov 12, 2024

glwagner commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024

glwagner commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024 •

edited

Loading

jagoosw commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024

simone-silvestri commented Nov 13, 2024 •

edited

Loading

simone-silvestri commented Nov 13, 2024 •

edited

Loading

ali-ramadhan commented Nov 13, 2024

simone-silvestri commented Nov 13, 2024

ali-ramadhan commented Nov 13, 2024

ali-ramadhan commented Nov 13, 2024

simone-silvestri commented Nov 13, 2024

simone-silvestri commented Nov 13, 2024

Some Lagrangian particles like to travel Southwest on an immersed lat-lon grid + hydrostatic model #3917

Some Lagrangian particles like to travel Southwest on an immersed lat-lon grid + hydrostatic model #3917

Comments

ali-ramadhan commented Nov 12, 2024 • edited Loading

ali-ramadhan commented Nov 12, 2024

glwagner commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024 • edited Loading

ali-ramadhan commented Nov 12, 2024

glwagner commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024

glwagner commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024 • edited Loading

jagoosw commented Nov 12, 2024

ali-ramadhan commented Nov 12, 2024

simone-silvestri commented Nov 13, 2024 • edited Loading

simone-silvestri commented Nov 13, 2024 • edited Loading

ali-ramadhan commented Nov 13, 2024

simone-silvestri commented Nov 13, 2024

ali-ramadhan commented Nov 13, 2024

ali-ramadhan commented Nov 13, 2024

simone-silvestri commented Nov 13, 2024

simone-silvestri commented Nov 13, 2024

ali-ramadhan commented Nov 12, 2024 •

edited

Loading

ali-ramadhan commented Nov 12, 2024 •

edited

Loading

ali-ramadhan commented Nov 12, 2024 •

edited

Loading

simone-silvestri commented Nov 13, 2024 •

edited

Loading

simone-silvestri commented Nov 13, 2024 •

edited

Loading