Skip to content

(0.95.13) Use active-weighted rather than naive reconstruction for f- and beta-plane Coriolis terms#4112

Merged
glwagner merged 29 commits intomainfrom
glw/active-weighted-interp
Feb 25, 2025
Merged

(0.95.13) Use active-weighted rather than naive reconstruction for f- and beta-plane Coriolis terms#4112
glwagner merged 29 commits intomainfrom
glw/active-weighted-interp

Conversation

@glwagner
Copy link
Copy Markdown
Member

@glwagner glwagner commented Feb 19, 2025 •
edited
Loading

Resolves #4111. Thanks for the bugfind @kenflat2!

Basically this implements boundary-aware reconstruction for Coriolis terms. Reconstruction is required for Coriolis terms on a staggered grid. Previously, reconstruction was not boundary-aware.

Note, I am not sure about ConstantCartesianCoriolis or NonTraditionalBetaPlane. Probably there is some work to do there, but reconstruction is different for those so I left it out of this PR.

Here's my updated script:

using Oceananigans
using Oceananigans.Units
using GLMakie
using Printf

Lx = 2000kilometers
Ly = 1000kilometers
Lz = 4kilometers

grid = RectilinearGrid(size = (128, 64, 64),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = (-Lz, 0),
                       halo = (4, 4, 4),
                       topology = (Periodic, Bounded, Bounded))

ridge(x, y) = -Lz + 1000meters * exp(-((x - Lx / 2)^2) / (100000kilometers^2))
ridge_grid = ImmersedBoundaryGrid(grid, GridFittedBottom(ridge))

coriolis = FPlane(f=-1e-4)
model = HydrostaticFreeSurfaceModel(; grid = ridge_grid, coriolis, momentum_advection = WENO(order=5))

uᵢ(x, y, z) = 0.1 * exp(z / 500)
set!(model, u=uᵢ)
simulation = Simulation(model, Δt=20minutes, stop_iteration=100)

wall_clock = Ref(time_ns())

function progress(sim)
    elapsed = 1e-9 * (time_ns() - wall_clock[])

    msg = @sprintf("iteration: %d, time: %s, wall time: %s, max|u|: (%6.3e, %6.3e, %6.3e), m s⁻¹\n",
                   iteration(sim), prettytime(sim), prettytime(elapsed),
                   maximum(abs, sim.model.velocities.u),
                   maximum(abs, sim.model.velocities.v),
                   maximum(abs, sim.model.velocities.w))

    wall_clock[] = time_ns()
    @info msg

    return nothing
end

add_callback!(simulation, progress, IterationInterval(10))

run!(simulation)

u, v, w = model.velocities
heatmap(view(v, :, 32, :))
current_figure()

on main

image

this PR

image

I'm not sure how to test this, ideas welcome.

@glwagner glwagner changed the title Use active-weighted reconstruction for f- and beta-plane Coriolis (0.95.13) Use active-weighted reconstruction for f- and beta-plane Coriolis Feb 19, 2025
simone-silvestri
simone-silvestri reviewed Feb 20, 2025
View reviewed changes
@simone-silvestri
Copy link
Copy Markdown
Collaborator

simone-silvestri commented Feb 20, 2025

Now that the scope for active_cell interpolation has expanded out of use for a specific type (ActiveCellEnstrophyConserving), we should probably move these functions in Operators (and add some tests) so that we can use them for other methods that require interpolation across immersed boundaries

@simone-silvestri
Copy link
Copy Markdown
Collaborator

simone-silvestri commented Feb 20, 2025

Looks like masking inactive nodes is always the numerically correct implementation and it is not a user choice, so I think we can remove ActiveCellEnsthrophyScheme and leave EnstrophyScheme where we mask inactive cells.

@glwagner
Copy link
Copy Markdown
Member Author

glwagner commented Feb 20, 2025

@simone-silvestri feel free to make these changes

@glwagner
Copy link
Copy Markdown
Member Author

glwagner commented Feb 20, 2025

@simone-silvestri feel free to make these changes

Nevermind I did it.

@glwagner
Copy link
Copy Markdown
Member Author

glwagner commented Feb 20, 2025

Now that the scope for active_cell interpolation has expanded out of use for a specific type (ActiveCellEnstrophyConserving), we should probably move these functions in Operators (and add some tests) so that we can use them for other methods that require interpolation across immersed boundaries

Not sure what you have in mind for test or something more general but feel free to add it.

Note that we've received feedback that metaprogramming makes it hard to understand Oceananigans code/errors and slows down development. So I think we should avoid it if we don't need it. Not sure how many of the "active-weighted" operators are needed.

@simone-silvestri
Copy link
Copy Markdown
Collaborator

simone-silvestri commented Feb 20, 2025

I was thinking to use it for sea ice, so I think we can add only the horizontal interpolations for now and add other ones as we see fit

@navidcy
Copy link
Copy Markdown
Member

navidcy commented Feb 21, 2025

Hm... I can't understand neither from this PR nor from the issue why Coriolis needs to be reconstructed.

@glwagner
Copy link
Copy Markdown
Member Author

glwagner commented Feb 21, 2025

Hm... I can't understand neither from this PR nor from the issue why Coriolis needs to be reconstructed.

The reconstruction is required because of staggering. For example computing f * v at f, c, c requires a 4-point stencil.

This PR does not change whether we reconstruct or not. Instead, it changes how we reconstruct by omitting inactive cells from the reconstruction.

As for why, the evidence is given by the empirical test. I'm not offering a theoretical explanation which would require more work...

@navidcy
Copy link
Copy Markdown
Member

navidcy commented Feb 21, 2025

Oh ok. I get it now. I wasn’t asking for deeper why.

@glwagner
Copy link
Copy Markdown
Member Author

glwagner commented Feb 21, 2025

Oh ok. I get it now. I wasn’t asking for deeper why.

updated the top-level description

@glwagner glwagner changed the title (0.95.13) Use active-weighted reconstruction for f- and beta-plane Coriolis (0.95.13) Use active-weighted rather than naive reconstruction for f- and beta-plane Coriolis terms Feb 21, 2025
glwagner
glwagner commented Feb 24, 2025
View reviewed changes
@glwagner glwagner merged commit 583dcd2 into main Feb 25, 2025
50 checks passed
@glwagner glwagner deleted the glw/active-weighted-interp branch February 25, 2025 01:01
@simone-silvestri simone-silvestri mentioned this pull request Mar 19, 2026
Merged
4 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@simone-silvestri simone-silvestri simone-silvestri approved these changes

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

Issue with Noise In Vorticity on Cartesian Grid with Immersed Bottom Boundary

3 participants

@glwagner @simone-silvestri @navidcy