Non-hydrostatic curvilinear shell?

[YoffeG]

Hi all,

We wish to simulate a part of a spherical shell (e.g., within lat +-70) in a non-hydrostatic regime. Is there a readily available curvilinear grid framework to implement for such a setup?

Thanks in advance,

Gidi

[glwagner]

Hi @YoffeG do you mind if I convert this to a discussion? I think our conversation will be better served by that format.

It sounds like you want to use a LatitudeLongitudeGrid. However, while we have a conjugate gradient solver for such a grid, we have never tested it. Note that all nonhydrostatic work I am aware of uses RectilinearGrid and a direct FFT solver. There is some more experimental work with the CG solver for nonhydrostatic simulations in complex domains.

It could be a fun research project to try to get the CG solver working on LatitudeLongitudeGrid if you are up for it.

[glwagner]

@YoffeG I am on here every day 😂

MITgcm uses a CG solver with an analytical preconditioner as described by Marshall et al 1997. However, a CG solver with a multigrid preconditioner would likely be more performant.

I personally don't have plans to work on physics for nonhydrostatic deep atmosphere/oceans because I am working on developing a thin-shell hydrostatic ocean climate model at the moment. But Oceananigans is designed to facilitate rapid development by... anyone... so the correct question is "does anybody have plans to implement this?" Obviously I don't know the answer! Perhaps, @Yixiao-Zhang or @jm-c have plans, or might know if any friends have plans to look into this. Here are some additional comments:

I'm not sure what a full Coriolis force entails. If it uses the existing grid, it is probably pretty easy. For example, this is the implementation of FPlane coriolis force:

Oceananigans.jl/src/Coriolis/f_plane.jl

Line 44 in 2c955da

@inline x_f_cross_U(i, j, k, grid, coriolis::FPlane, U) = - coriolis.f * ℑxyᶠᶜᵃ(i, j, k, grid, U[2])

However if this is more than simply a new stencil calculation, but also requires one not to make the "thin shell approximation" when constructing the grid, then more substantial changes may be required. I am not sure how hard this will be, since much of the code was written to eventually support generalization to such a case. However, whether or not the design is successful can only be revealed by trying it.

In the interest of providing a complete description of what this would entail, I believe that 1) formulas to compute the vertically-varying horizontal spacings must be provided and 2) the functions that access grid metrics need to be modified, eg changing

Oceananigans.jl/src/Operators/spacings_and_areas_and_volumes.jl

Line 122 in 2c955da

@inline Δxᶜᶠᵃ(i, j, k, grid::LLG) = @inbounds grid.Δxᶜᶠᵃ[i, j]

[glwagner]

Out of curiosity...

using Oceananigans
using Oceananigans.Solvers: solve!, ConjugateGradientPoissonSolver
using GLMakie

grid = LatitudeLongitudeGrid(size=(120, 40, 5),
                             latitude = (-60, 60),
                             longitude = (0, 360),
                             z = (0, 1))

solver = ConjugateGradientPoissonSolver(grid)

r = solver.right_hand_side
ϕ = CenterField(grid)
set!(r, (λ, φ, z) -> exp(-(λ - 180)^2 / 50 - φ^2 / 50))

solve!(ϕ, solver.conjugate_gradient_solver, r)

fig = Figure(size=(800, 400))
axr = Axis(fig[1, 1], title="Source term, r")
axϕ = Axis(fig[1, 2], title="Poisson solution ϕ to ∇²ϕ = r")
heatmap!(axr, view(r, :, :, 1))
heatmap!(axϕ, view(ϕ, :, :, 1))

[glwagner]

Also @YoffeG I may have misunderstood what you said --- if you are merely talking about the "non-traditional" terms in a spherical Coriolis formulation for a NonhydrostaticModel, I think that is fairly trivial and could be implemented easily. It just means extending the HydrostaticSphericalCoriolis to include the extra vertical term, and our Coriolis infrastructure already generalizes to vertical terms. Happy to help with that if you want to try this.