Implement horizontal derivatives at constant z for `MutableVerticalDiscretization` by glwagner · Pull Request #5101 · CliMA/Oceananigans.jl

glwagner · 2026-01-06T05:16:27Z

This PR updates the implementation of the generalized vertical discretization (what we now call MutableVerticalDiscretization) so that horizontal derivatives are interpreted as horizontal derivatives at constant z. I believe this is necessary for code to produce similar / same outcomes for any vertical coordinate. In other words, if we have a static vertical discretization, we assume that horizontal derivatives apply at constant z. So when we change vertical coordinates, we need to specify that derivatives still should apply at constant z.

The fact that this simplifies the code is demonstrated in how it allows us to remove a term from the tendency. Previously we had a term grid_slope_contribution_x... which part of the horizontal pressure gradient. Now that term is incorporated into ∂xᶠᶜᶜ.

As an aside, I think it is more efficient to evaluate dpdz directly, rather than recomputing buoyancy (which can be expensive, eg involving a 55 term formula).

I'm still not happy with how generalized coordinates are implemented, but I think additional changes should await a future PR and this is a good start.

simone-silvestri · 2026-01-06T10:59:03Z

looks, good, make sure that these horizontal derivatives are not used anywhere in the dynamics where there are not intended to, and that this will work for fields that do not depend on z like the free surface (I think with this formulation will probably fail with an out-of-bounds error)

simone-silvestri · 2026-01-06T11:22:14Z

Probably we need to define two-dimensional derivatives ∂xᶜᶜᵃ to use for two-dimensional fields

glwagner · 2026-01-06T19:20:49Z

looks, good, make sure that these horizontal derivatives are not used anywhere in the dynamics where there are not intended to

if there is anywhere in the dynamics that uses dx at constant r, doesn't that mean that those are currently incorrect for z* vertical coordinate?

For example, with this definition there is no "additional" term needed to compute the horizontal pressure gradient. I expect the same is true for any gradient.

glwagner · 2026-01-06T20:44:02Z

this will work for fields that do not depend on z like the free surface (I think with this formulation will probably fail with an out-of-bounds error)

for z-reduced fields, the operators like ∂zᶠᶠᶜ return 0. I don't think there is the possibility of an out of bounds error.

I think defining the reduced derivatives is probably out of scope for this PR. In practice, the third location in a horizontal derivative of a z-reduced field is ignored.

simone-silvestri · 2026-01-08T16:31:19Z

I think this is not always the case, because we should do the same with the time derivative, aka we need the time derivative at constant r, not at constant z (we are time stepping in r-coordinate), which would be a bit inconvenient to express without some simplifications in the material derivatives to compensate. The dx at constant r are there because they elide with the chain rule for the time derivative. I think it is well documented here, see for example the derivations between equation 11 and 13 (the equivalent for momentum is 17 to 19). The pressure gradient (20)-(21) is a simple term that does not simplify, so the extra term remains hanging and for this reason we can use the chain rule. Also the vertical velocity (8) is not the same quantity in r-coordinate and z-coordinate, which implies other simplifications.

glwagner · 2026-01-08T21:57:02Z

I think this is not always the case, because we should do the same with the time derivative, aka we need the time derivative at constant r, not at constant z (we are time stepping in r-coordinate), which would be a bit inconvenient to express without some simplifications in the material derivatives to compensate. The dx at constant r are there because they elide with the chain rule for the time derivative. I think it is well documented here, see for example the derivations between equation 11 and 13 (the equivalent for momentum is 17 to 19). The pressure gradient (20)-(21) is a simple term that does not simplify, so the extra term remains hanging and for this reason we can use the chain rule. Also the vertical velocity (8) is not the same quantity in r-coordinate and z-coordinate, which implies other simplifications.

I am possibly confused, but if I understand the point of the derivation there, it is to write down the correct equations in $\xi, \eta, r$ given the equations of motion in $x, y, z$. For example, if we compute the divergence with $\delta = \partial_x u + \partial_y v + \partial_z w$, then this computation remains a valid expression for the divergence in any coordinate system. Whereas, before this PR, $\delta$ is only correct with ZCoordinate and is incorrect in ZStarCoordinate --- as very clearly illustrated by the derivation you have linked. So what do you mean by saying "this is not always the case"?

I am not sure the motivation for rewriting everything in terms of derivatives at constant $r$. Is it for computational efficiency, or numerical stability / accuracy?

This PR proposes to fix the meaning of derivatives such as $\partial_x$ as being defined at constant z. This means they ahve the same interpretation regardless of the choice of coordinate.

glwagner · 2026-01-19T04:01:05Z

@simone-silvestri do you approve?

simone-silvestri · 2026-01-19T15:41:13Z

I think windowed fields like the free surface would still error when being derived in x because it does not have the extra term in z?

I am not sure the motivation for rewriting everything in terms of derivatives at constant r. Is it for computational efficiency, or numerical stability / accuracy?

This PR proposes to fix the meaning of derivatives such as ∂x as being defined at constant z. This means they ahve the same interpretation regardless of the choice of coordinate.

This might be ok because we do not use the derivative operator for the advection, but if we want to use that approach, we should do the same for the time derivative, which is not really expressed in terms of operators. In fact we do not timestep at constant z but at constant r (the z changes between time steps in the same cell). So we need to reformulate the equations in r to make sure the time derivative is correct. Timestepping at constant z in a zstar coordinate would be a lagrangian approach (probably I am not 100% sure).

simone-silvestri · 2026-01-19T15:42:23Z

Now for diagnostics this approach is correct, but we need to make sure that dynamics are solved in r-coordinate, not in z-coordinate.

glwagner · 2026-01-19T16:24:32Z

This might be ok because we do not use the derivative operator for the advection, but if we want to use that approach, we should do the same for the time derivative, which is not really expressed in terms of operators

This is true, but arguably time information should be encapsulated in the coordinate, not in the discretization. Time stepping is the domain of models. At least this PR will fix static generalized vertical discretizations.

glwagner · 2026-01-19T16:25:39Z

Now for diagnostics this approach is correct, but we need to make sure that dynamics are solved in r-coordinate, not in z-coordinate.

This PR does not change the fact that dynamics must be implemented correctly (this was also true before this PR)

glwagner · 2026-01-19T16:28:58Z

I think windowed fields like the free surface would still error when being derived in x because it does not have the extra term in z?

I think you're right -- if we define dx as "derivative in x at constant z" then this is not a valid operator to apply to the free surface. We need dx_r (derivative in x at constant r).

simone-silvestri · 2026-01-19T16:39:10Z

Ok I guess this is a good compromise, for static generalized coordinates we get the discretization for free, for moving ones we need to make sure that the time derivative is correctly discretized

glwagner · 2026-01-19T17:40:25Z

What do you view this as a compromise between?

simone-silvestri · 2026-01-19T20:38:01Z

I guess the amount of work required to have a new coordinate working with the internals of a model and the correctness of the diagnostics associated with the output of the model.

glwagner · 2026-01-20T04:39:13Z

I guess the amount of work required to have a new coordinate working with the internals of a model and the correctness of the diagnostics associated with the output of the model.

How does this increase the amount of work required to get a coordinate working? It decreased the amount of work required for ZStar. What coordinates does this create a problem for?

glwagner · 2026-01-20T04:48:03Z

The primary motivation for this PR is to make it trivial to implement terrain-following static coordinates. This PR eliminates the need for special implementations of horizontal derivatives.

I do not understand why a generalizing abstraction wasn't useful for ZStarVerticalCoordinate. At this point, it may be true that we need to revert features implemented for ZStar in order to make this PR work. But that is not a deficiency of the abstraction. That is a problem with the original implementation of ZStar -- rather than finding a general abstraction that would work for both ZStar and other future coordinate implementations, a different, non-general route was taken that, if not reverted, would lead to code duplication. That is of course the reason to solve abstraction problems sooner rather than later. The central point of an abstraction is that it can be used to provide consistent mathematical representations of operations. In other codes, perhaps internal implementation of dynamics are sufficient. In our code, we need abstractions that can be used by both dynamics and users.

Note that it is completely unacceptable that diagnostics become meaningless when switching coordinates. Diagnostics are actually more important than dynamics implementation. A dynamics implementation is pointless if the code cannot be used for science.

simone-silvestri · 2026-01-20T10:51:04Z

I realize maybe "compromise" is not the correct term, the thing is that z-based horizontal derivatives do not solve the issue of implementing generic coordinates, that's just what I wanted to express.

I think it might happen that you need to explicitly avoid using these operators in many terms (which implies building new operators at constant r) to have correct dynamics. Most of the zstar formulation does not use these derivatives (w is not the vertical velocity, the time derivative is at constant r, and so on...), it is only a matter of the pressure gradient (and closures!). Maybe for terrain following coordinates this is all that is needed and that is great!

I agree with pushing forward this PR and I never claimed that diagnostics are not important, so I see this PR more as a positive addition for diagnostics rather than a "generalized abstraction" for implementing vertical coordinates.

We need to solve the problem of windowed fields like the free surface though: this PR at the moment is removing support for horizontal derivatives of the free surface in generalized coordinates (https://buildkite.com/clima/oceananigans/builds/28653#019bd4b7-503e-4e26-9f3f-1f52b336c3e5)

Indeed, we should probably have two sets of derivatives, one in constant r (just difference divided by spacing) and one in constant z.

simone-silvestri · 2026-01-20T12:34:41Z

We need to solve the problem of windowed fields like the free surface though: this PR at the moment is removing support > for horizontal derivatives of the free surface in generalized coordinates (https://buildkite.com/clima/oceananigans/builds/28653#019bd4b7-503e-4e26-9f3f-1f52b336c3e5)

I can give it a go?

glwagner · 2026-01-20T16:11:32Z

We need to solve the problem of windowed fields like the free surface though: this PR at the moment is removing support > for horizontal derivatives of the free surface in generalized coordinates (https://buildkite.com/clima/oceananigans/builds/28653#019bd4b7-503e-4e26-9f3f-1f52b336c3e5)

I can give it a go?

I am happy to tackle this, but also happy to accept help if you want.

I think it might happen that you need to explicitly avoid using these operators in many terms (which implies building new operators at constant r)

We should of course add dx_r and dy_r. These will be one-liners.

Most of the zstar formulation does not use these derivatives (w is not the vertical velocity, the time derivative is at constant r, and so on...), it is only a matter of the pressure gradient (and closures!).

Are you sure that is the optimal way to implement it? But indeed, these operators are also needed to compute horizontal buoyancy gradients and isopycnal slopes.

the thing is that z-based horizontal derivatives do not solve the issue of implementing generic coordinates, that's just what I wanted to express.

The critical question is not whether this is a complete solution but rather whether it represents a step towards a system of abstractions that is a complete solution.

We cannot solve issues with moving coordinates purely with spatial derivative operators. However, if we do open the support of static, terrain following coordinates with these operators, I would argue that this is the correct abstraction strategy.

simone-silvestri · 2026-01-21T11:08:25Z

Ok, I guess to push this PR, we can just substitute the derivative in the explicit_barotropic_pressure_*_gradient and then open another PR for the r-based derivatives.

glwagner · 2026-01-21T20:44:40Z

src/Models/HydrostaticFreeSurfaceModels/explicit_free_surface.jl


 @inline explicit_barotropic_pressure_y_gradient(i, j, k, grid, free_surface::ExplicitFreeSurface) =
-    free_surface.gravitational_acceleration * ∂yᶜᶠᶜ(i, j, grid.Nz+1, grid, free_surface.displacement)
+    free_surface.gravitational_acceleration * δyᶜᶠᶜ(i, j, grid.Nz+1, grid, free_surface.displacement) * Δyᶜᶠᶜ⁻¹(i, j, k, grid)


@simone-silvestri do you mind if I add operators for constant_r derivatives?

Ah no go ahead!

glwagner · 2026-01-24T20:29:15Z

Also @simone-silvestri did you see this comment?

As an aside, I think it is more efficient to evaluate dpdz directly, rather than recomputing buoyancy (which can be expensive, eg involving a 55 term formula).

glwagner · 2026-01-24T20:30:17Z

Also @simone-silvestri did you see this comment?

As an aside, I think it is more efficient to evaluate dpdz directly, rather than recomputing buoyancy (which can be expensive, eg involving a 55 term formula).

Resolve conflicts: - Project.toml: take main version (0.105.0), remove old extras/targets sections - hydrostatic_free_surface_tendency_kernel_functions.jl: keep removal of grid_slope_contribution terms (this PR's intent) while adopting the diffusivities→closure_fields rename from main Co-Authored-By: Claude Opus 4.6 <[email protected]>

Add explicit methods for Number arguments on AbstractMutableGrid and MutableImmersedGrid derivative operators to resolve method ambiguities with the base derivative_operators.jl Number methods. Co-Authored-By: Claude Opus 4.6 <[email protected]>

* main: Add fields with non-default indices to Field tutorial (#5308) Implement horizontal derivatives at constant z for `MutableVerticalDiscretization` (#5101)

briochemc · 2026-02-26T10:58:10Z

src/Operators/generalized_derivatives.jl

+@inline function ∂xᶠᶜᶜ(i, j, k, grid::AMG, ϕ)
+    ∂x_at_r = δxᶠᶜᶜ(i, j, k, grid, ϕ) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid)
+    ∂z_ϕ = ∂zᶠᶜᶜ(i, j, k, grid, ϕ)
+    ∂x_z = ∂x_zᶠᶜᶜ(i, j, k, grid)
+    return ∂x_at_r - ∂x_z * ∂z_ϕ
+end


Is that correct in terms of staggering?

When I look at this function signature I think that ϕ must be at a CCC location for its ∂xᶠᶜᶜ(ϕ) to be located at FCC. But inside the body, ϕ seems to be located at a FCF location for its ∂zᶠᶜᶜ(ϕ) to be at FCC:

@inline function ∂xᶠᶜᶜ(i, j, k, grid::AMG, ϕ) # <- ϕ at CCC location ∂x_at_r = δxᶠᶜᶜ(i, j, k, grid, ϕ) * Δx⁻¹ᶠᶜᶜ(i, j, k, grid) ∂z_ϕ = ∂zᶠᶜᶜ(i, j, k, grid, ϕ) # <- ϕ at FCF location? ∂x_z = ∂x_zᶠᶜᶜ(i, j, k, grid) return ∂x_at_r - ∂x_z * ∂z_ϕ end

AFAIU ∂zᶠᶜᶜ(ϕ) "siliently shifts" ϕ west and below by half an index (CCC->FCF). I understand this could be a good-enough approximation but if the function needs to be correct, shouldn't there be some CCC->FCF interpolation in there?

But maybe I am missing something in how these derivative and difference operators work? Please let me know!

cc @glwagner @simone-silvestri

(FWIW I noticed that because the transport matrix I built became structurally non-symmetric, which should not happen in my model setup.)

yeah this looks wrong.

I think we have to compute the z derivative and then do a 2D interpolation in x and z to obtain the z-derivative at fcc?

Will this become a bit too computationally expensive? (horizontal interpolations of vertical derivatives are terribly slow)

is there another solution?

I thought this might have been intentional for performance actually, but I didn't know it would be dramatic. Would it make any difference to do the horizontal interpolation first and then take the vertical derivatives?

simone-silvestri · 2026-03-03T16:20:29Z

I am testing a one degree resolution ocean model and it seems like with the new derivatives GM / Redi goes crazy, I will try again with #5339 to see if it is an issue of z-derivatives not at the correct position (hopefully). Otherwise I will try to distill a MWE and open an issue.

glwagner · 2026-03-03T18:56:58Z

I am testing a one degree resolution ocean model and it seems like with the new derivatives GM / Redi goes crazy, I will try again with #5339 to see if it is an issue of z-derivatives not at the correct position (hopefully). Otherwise I will try to distill a MWE and open an issue.

It is because there is a huge bug in the pressure gradient

simone-silvestri · 2026-03-04T08:36:09Z

the bug in the pressure gradient does not change much the solution. It's really GM which is completely unphysical:

The difference between these two free surfaces (evolved with atmospheric forcing for one year) is just adding GM and Redi. I think it might be fixed with #5339. However I am thinking whether we can just embed the chain rule in the slopes:
$$S_y = \partial y_r(b) / \partial z(b) - \partial y(z)$$
$$S_x = \partial x_r(b) / \partial z(b) - \partial x(z)$$

I think this will be much more efficient and avoid us all those interpolations which GM is full of already

implement horizontal derivatives at constant z

b60c56f

glwagner requested a review from simone-silvestri January 6, 2026 05:16

glwagner changed the title ~~Implement horizontal derivatives at constant z~~ Implement horizontal derivatives at constant z for MutableVerticalDiscretization Jan 6, 2026

Merge branch 'main' into glw/better-generalized-coords

815be01

glwagner and others added 3 commits January 6, 2026 13:44

fix tests

98d084d

remove stale import

345055b

Merge branch 'main' into glw/better-generalized-coords

2ac3c52

simone-silvestri mentioned this pull request Jan 12, 2026

(0.104.0) Reformulate hydrostatic model timestepping #4811

Merged

9 tasks

Merge main into glw/better-generalized-coords

abebc36

simone-silvestri approved these changes Jan 19, 2026

View reviewed changes

explicit derivative in eta

d192068

simone-silvestri added 4 commits January 21, 2026 12:18

reciprocal in front

af65b2f

fix more derivatives

6edb444

import correct operators

9933ed9

remove other derivatives

c59a0f9

glwagner commented Jan 21, 2026

View reviewed changes

glwagner added 6 commits January 21, 2026 12:52

add derivatives at constant r

baeeb25

Merge branch 'main' into glw/better-generalized-coords

f221ae1

fix

c913015

fix ambiguities

16a7b8a

bump to 0.104.3

3acd276

Merge branch 'main' into glw/better-generalized-coords

5f7b977

glwagner and others added 2 commits February 17, 2026 19:48

glwagner merged commit 957587e into main Feb 18, 2026
77 of 78 checks passed

glwagner deleted the glw/better-generalized-coords branch February 18, 2026 12:21

briochemc added a commit that referenced this pull request Feb 18, 2026

Merge branch 'main' into glw/prescribed-free-surface

a476387

* main: Add fields with non-default indices to Field tutorial (#5308) Implement horizontal derivatives at constant z for `MutableVerticalDiscretization` (#5101)

briochemc reviewed Feb 26, 2026

View reviewed changes

glwagner mentioned this pull request Feb 26, 2026

(0.105.3) Fix interpolation of terms in chain horizontal derivatives for MutableVerticalDiscretization #5339

Merged

simone-silvestri mentioned this pull request Mar 17, 2026

(0.106.2) Refactor coriolis and add Triad scheme #5383

Merged

4 tasks

Conversation

glwagner commented Jan 6, 2026

Uh oh!

simone-silvestri commented Jan 6, 2026

Uh oh!

simone-silvestri commented Jan 6, 2026

Uh oh!

glwagner commented Jan 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

glwagner commented Jan 6, 2026

Uh oh!

simone-silvestri commented Jan 8, 2026

Uh oh!

glwagner commented Jan 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

glwagner commented Jan 19, 2026

Uh oh!

simone-silvestri commented Jan 19, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

simone-silvestri commented Jan 19, 2026

Uh oh!

glwagner commented Jan 19, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

glwagner commented Jan 19, 2026

Uh oh!

glwagner commented Jan 19, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

simone-silvestri commented Jan 19, 2026

Uh oh!

glwagner commented Jan 19, 2026

Uh oh!

simone-silvestri commented Jan 19, 2026

Uh oh!

glwagner commented Jan 20, 2026

Uh oh!

glwagner commented Jan 20, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

simone-silvestri commented Jan 20, 2026

Uh oh!

simone-silvestri commented Jan 20, 2026

Uh oh!

glwagner commented Jan 20, 2026

Uh oh!

simone-silvestri commented Jan 21, 2026

Uh oh!

glwagner Jan 21, 2026

Choose a reason for hiding this comment

Uh oh!

simone-silvestri Jan 21, 2026

Choose a reason for hiding this comment

Uh oh!

glwagner commented Jan 24, 2026

Uh oh!

glwagner commented Jan 24, 2026

Uh oh!

Uh oh!

briochemc Feb 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

glwagner Feb 26, 2026

Choose a reason for hiding this comment

Uh oh!

simone-silvestri Feb 26, 2026

Choose a reason for hiding this comment

Uh oh!

glwagner Feb 26, 2026

Choose a reason for hiding this comment

Uh oh!

briochemc Feb 26, 2026

Choose a reason for hiding this comment

Uh oh!

simone-silvestri commented Mar 3, 2026

Uh oh!

glwagner commented Jan 6, 2026 •

edited

Loading

glwagner commented Jan 8, 2026 •

edited

Loading

simone-silvestri commented Jan 19, 2026 •

edited

Loading

glwagner commented Jan 19, 2026 •

edited

Loading

glwagner commented Jan 19, 2026 •

edited

Loading

glwagner commented Jan 20, 2026 •

edited

Loading

briochemc Feb 26, 2026 •

edited

Loading