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

[ali-ramadhan]

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:

https://github.com/user-attachments/assets/77d17470-2dae-44ad-a38b-a026cce1f02e

[ali-ramadhan]

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

[glwagner]

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

[ali-ramadhan]

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

:thinking:

[ali-ramadhan]

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]

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]

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

https://github.com/user-attachments/assets/b835c970-a76e-406d-a10c-b21aa824c72b

https://github.com/user-attachments/assets/2bff08ab-3398-4178-8719-17ba06f5ff93

[ali-ramadhan]

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

[glwagner]

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

[ali-ramadhan]

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]

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]

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]

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]

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

https://github.com/user-attachments/assets/2ffb3b2f-01e3-4479-9486-cba5519eb23f

https://github.com/user-attachments/assets/c142c911-3040-40d4-adbf-125daf61c875

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

[ali-ramadhan]

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]

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]

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

[ali-ramadhan]

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.

https://github.com/user-attachments/assets/d575d04e-4c29-43e3-aacb-c0f0d592fe17

https://github.com/user-attachments/assets/4c03866d-7d34-481c-955c-37168f15dc04

[simone-silvestri]

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

[simone-silvestri]

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

https://github.com/user-attachments/assets/0ad263c2-f42e-4a25-9966-321e772d3b78