Exports more gradient operators from SpeedyTransforms

[milankl]

@miniufo as discussed some convience functions to compute the gradient of a LowerTriangularMatrix or AbstractGrid. Also other operators ∇ ∇² ∇! ∇⁻² ∇²! ∇⁻²! are now all exported. Gradient can be used as

# create some data with wavenumbers 0,1,2,3,4
L = randn(LowerTriangularMatrix{ComplexF32}, 33, 32)
spectral_truncation!(L, 5)

# now assume you have that data in grid-point space as G
G = gridded(L)

# calculate the gradient (optionally pass on a spectral transform as second argument for performance)
dGdx, dGdy = ∇(G)

as always, radius scaling is omitted, so you'd need to /radius manually. Results in

More methods in the docstring

help?> ∇
"∇" can be typed by \del<tab>

search: ∇ ∇² ∇! ∇⁻² ∇²! ∇⁻²!

  ∇(
      grid::Main.SpeedyWeather.RingGrids.AbstractGrid,
      S::SpectralTransform
  ) -> Tuple{Any, Any}

  The zonal and meridional gradient of grid. Transform to spectral space, takes the gradient
  and unscales the 1/coslat scaling in the gradient, but omits the 1/radius scaling that has to
  be applied manually. Makes use of an existing spectral transform S.

  ─────────────────────────────────────────────────────────────────────────────────────────────

  ∇(
      grid::Main.SpeedyWeather.RingGrids.AbstractGrid
  ) -> Tuple{Any, Any}

  The zonal and meridional gradient of grid. Transform to spectral space, takes the gradient
  and unscales the 1/coslat scaling in the gradient, but omits the 1/radius scaling that has to
  be applied manually.

  ─────────────────────────────────────────────────────────────────────────────────────────────

  ∇(p::LowerTriangularMatrix)

  The zonal and meridional gradient of p. Precomputes a SpectralTransform S. The
  1/radius-scaling is omitted, as is the 1/coslat scaling.

  ─────────────────────────────────────────────────────────────────────────────────────────────

  ∇(
      p::LowerTriangularMatrix,
      S::SpectralTransform
  ) -> Tuple{LowerTriangularMatrix{Complex{NF}} where NF<:AbstractFloat, LowerTriangularMatrix{Complex{NF}} where NF<:AbstractFloat}

  The zonal and meridional gradient of p using an existing SpectralTransform S. The
  1/radius-scaling is omitted, as is the 1/coslat scaling.

[miniufo]

That's great! Now it looks very close to the math expression.

[milankl]

Yes that's intended but I'm still struggling to make our functions consistent with respect to radius and or coslat scaling. Maybe you have ideas/opinions on that, for context

• None of the internal gradient operators divide or multiply by the radius, all act on the unit sphere as most of these scalings cancel out in the dynamical core
• as you can only efficiently apply a coslat scaling in grid-point space this is also omitted for internal operators to avoid back and forth tranforms just for this. There's therefore only deliberate places where we have to unscale coslats and we chain them together for efficiency
• the operators for users (like you) are convenience functions around those internal functions to allocate the result or precompute the transform

Now for the user-facing functions I keep facing questions of

• do we automatically unscale the coslat if we can without a transform (yes?)
• do we never scale by the radius and completely leave this to the user to be consistent? (yes?)
• when do we return a spectral field when do we return a grid point field? Grid in grid out? Spectral in spectral out? Make this a keyword argument? Or just leave any final transform to the user?

Feel free to let us know what you would find most intuitive and clear or what you think a general user would expect?

[miniufo]

Yes, I see your "struggling" too. When I design some gradient APIs in python, I need to add a kwarg as

grdx, grdy = gradient(v, Rearth=6371e3)

so that users could change the default radius of earth (like scaling or not here, they may choose a slightly value of Rearth).

In such function-style APIs, one needs this "extra" kwarg in every function (gradient, laplacian, divergence etc.), which is bit annoying. So a solution would be define an Application class (OOP-style), and remember users' choice at the start when Application is initialized. I do this sometimes. I also tried to define some internal function with long list of args so that it is flexible for developers (e.g., grdx, grdy, grdz = __gradient(v, dims=['x', 'y', 'x'], Rearth=6371e3)), and then exposed them to users with pre-defined default parameters (e.g., grdx, grdy = hgrad(v, Rearth=6371e3), which uses __gradient) so that many details are hidden.

I feel the case here is a bit similar. For users with gridded data, I feel like hidding all the scaling is better, but also allow users to change Rearth, as users do not need to know the details of the spectral transform. For (advanced) users dealing with spectral data, hidding all scaling is not possible, as cos-scaling is intrinsic in spectral method (for velocity, and users need to know this). So let users take this responsibility should be fine. Although this cos/radius scaling should not be an overhead for diagnostic purpose, it is much better to store the pre-computed spectral coefficient for faster calculation (this seems to fit the OOP-style).

As one of the users, I currently prefer:

1. Grid in/grid out, without lossing the information/accuracy;
2. Scaling, and other details of spectral transform, being transparent to me;
3. APIs being as simple as possible (not very long list of arguments), and close to maths;

Glad you share your thoughts, but I feel like never find a perfect solution (╥﹏╥).

[milankl]

Awesome thanks for the input, I've put this example into the docs: https://speedyweather.github.io/SpeedyWeather.jl/previews/PR523/gradients/

• All methods exposed to the user (the allocating ones without the in-place !) take now an optional keyword argument radius so that dGdx, dGdy = ∇(G, radius=6.371e6) would compute the gradient on Earth in 1/meters. I don't like the idea of having a somewhat global radius one can set as it also violates the (Julia / functional programming) philosophy to avoid globals as much as possible to have functions that only depend on their arguments, not on some "spooky action at a distance". By default the radius is 1.

• Grid in, grid out and then the coslat-scaling is also directly implemented, when providing spectral fields as inputs the scaling isn't done as the gradient operators (all acting on spectral fields) would require another two transforms to account for scaling. I think that's a consistent logic we can probably stick to (grid = with scaling, spectral = without).

• I've (hopefully all correct) added these specifics to every docstring too, e.g. ?∇, ?divergence etc.