diff --git a/docs/src/spectral_transform.md b/docs/src/spectral_transform.md index 9825a6737..d6713099b 100644 --- a/docs/src/spectral_transform.md +++ b/docs/src/spectral_transform.md @@ -292,44 +292,14 @@ v &= +\frac{1}{R}\partial_\theta\Phi + \frac{1}{R\cos\theta}\partial_\lambda\Psi \end{aligned} ``` -Alternatively, we can use the velocities ``U = u\cos\theta, V = v\cos\theta``, which we do as the meridional gradient -for spherical harmonics is easier implemented with a ``\cos\theta``-scaling included, and because the divergence and -curl in spherical coordinates evaluates the meridional gradient with ``U,V`` and not ``u,v``. From ``u,v`` we can -return to ``\zeta, \mathcal{D}`` via - -```math -\begin{aligned} -\zeta &= \frac{1}{R\cos\theta}\partial_\lambda v - \frac{1}{R\cos\theta}\partial_\theta (u \cos\theta) \\ -\mathcal{D} &= \frac{1}{R\cos\theta}\partial_\lambda u + \frac{1}{R\cos\theta}\partial_\theta (v \cos\theta). -\end{aligned} -``` - -Equivalently, we have - -```math -\begin{aligned} -U &= -\frac{\cos\theta}{R}\partial_\theta\Psi + \frac{1}{R}\partial_\lambda\Phi \\ -V &= +\frac{\cos\theta}{R}\partial_\theta\Phi + \frac{1}{R}\partial_\lambda\Psi \\ -\zeta &= \frac{1}{R}\partial_\lambda \left( \frac{V}{\cos^2\theta} \right) - -\frac{\cos\theta}{R}\partial_\theta \left( \frac{U}{\cos^2\theta} \right) \\ -\mathcal{D} &= \frac{1}{R}\partial_\lambda \left( \frac{U}{\cos^2\theta} \right) + -\frac{\cos\theta}{R}\partial_\theta \left( \frac{V}{\cos^2\theta} \right). -\end{aligned} -``` - -which is a more convenient formulation because of the way how the [Meridional derivative](@ref) -is implemented with a recursion relation, actually computing ``\cos\theta \partial_\theta`` -rather than ``\partial_\theta`` directly. The remaining cosine scalings in -``(U,V)*\cos^{-2}\theta`` are done in grid-point space. -If one wanted to get back to ``\zeta, \mathcal{D}`` this is how it would be done, but -it is often more convenient to unscale ``U,V`` on the fly in the spectral transform -to obtain ``u,v`` and then divide again by ``\cos\theta`` when any gradient (or divergence or -curl) is taken. This is because other terms would need that single ``\cos\theta`` unscaling -too before a gradient is taken. How the operators ``\nabla, \nabla \times, \nabla \cdot`` can +How the operators ``\nabla, \nabla \times, \nabla \cdot`` can be implemented with spherical harmonics is presented in the following sections. - -Also note that SpeedyWeather.jl scales the equations with the radius `R` (see [Radius scaling](@ref scaling)) -such that the divisions by `R` drop out in this last formulation too. +However, note that the actually implemented operators differ slightly in their +scaling with respect to the radius ``R`` and the cosine of latitude ``\cos(\theta)``. +For further details see [Gradient operators](@ref) which describes those as implemented +in the [SpeedyTransforms](@ref) module. Also note that the equations in SpeedyWeather.jl +are scaled with the radius ``R^2`` (see [Radius scaling](@ref scaling)) +which turns most operators into non-dimensional operators on the unit sphere anyway. ### Zonal derivative diff --git a/docs/src/speedytransforms.md b/docs/src/speedytransforms.md index 47cea3e60..948adf267 100644 --- a/docs/src/speedytransforms.md +++ b/docs/src/speedytransforms.md @@ -11,6 +11,19 @@ for clarifications how to work with these. We will also not discuss mathematical of the [Spherical Harmonic Transform](@ref) here, but will focus on the usage of the SpeedyTransforms module. +The SpeedyTransforms module also implements the gradient operators +``\nabla, \nabla \cdot, \nabla \times, \nabla^2, \nabla^{-2}`` in spectral space. +Combined with the spectral transform, you could for example start with a velocity +field in grid-point space, transform to spectral, compute its divergence and transform +back to obtain the divergence in grid-point space. Examples are outlined in [Gradient operators](@ref). + +!!! info "SpeedyTransforms assumes a unit sphere" + The operators in SpeedyTransforms generally assume a sphere of radius ``R=1``. + For the transforms themselves that does not make a difference, but the gradient operators + `div`,`curl`,`∇`,`∇²`,`∇⁻²` omit the radius scaling. You will have to do this manually. + Also note that the meridional derivate expects a ``\cos^{-1}(\theta)`` scaling. + More in [Gradient operators](@ref). + ## Example transform Lets start with a simple transform. We could be `using SpeedyWeather` but to be more verbose @@ -249,6 +262,152 @@ to kilobytes SpectralTransform(spectral_grid,recompute_legendre=true) ``` +## Gradient operators + +`SpeedyTransforms` also includes many gradient operators to take derivatives in +spherical harmonics. These are in particular ``\nabla, \nabla \cdot, \nabla \times, +\nabla^2, \nabla^{-2}``. However, the actually implemented operators are, +in contrast to the mathematical [Derivatives in spherical coordinates](@ref) +due to reasons of scaling as follows. Let the implemented operators be +``\hat{\nabla}`` etc. + +```math +\hat{\nabla} A = \left(\frac{\partial A}{\partial \lambda}, \cos(\theta)\frac{\partial A}{\partial \theta} \right) = +R\cos(\theta)\nabla A +``` +So the zonal derivative omits the radius and the ``\cos^{-1}(\theta)`` scaling. +The meridional derivative adds a ``\cos(\theta)`` due to a recursion relation +being defined that way, which, however, is actually convenient because the whole +operator is therefore scaled by ``R\cos(\theta)``. The curl and divergence operators +expect the input velocity fields to be scaled by ``\cos^{-1}(\theta)``, i.e. + +```math +\begin{aligned} +\hat{\nabla} \cdot (\cos^{-1}(\theta)\mathbf{u}) &= \frac{\partial u}{\partial \lambda} + +\cos\theta\frac{\partial v}{\partial \theta} = R\nabla \cdot \mathbf{u}, \\ +\hat{\nabla} \times (\cos^{-1}(\theta)\mathbf{u}) &= \frac{\partial v}{\partial \lambda} - +\cos\theta\frac{\partial u}{\partial \theta} = R\nabla \times \mathbf{u}. +\end{aligned} +``` + +And the Laplace operators omit a ``R^2`` (radius ``R``) scaling, i.e. + +```math +\hat{\nabla}^{-2}A = \frac{1}{R^2}\nabla^{-2}A , \quad \hat{\nabla}^{2}A = R^2\nabla^{2}A +``` + +## Example: Geostrophy + +Now, we want to use the following example to illustrate +their use: We have ``u,v`` and want to calculate ``\eta`` in the shallow water +system from it following geostrophy. Analytically we have +```math +-fv = -g\partial_\lambda \eta, \quad fu = -g\partial_\theta \eta +``` +which becomes, if you take the divergence of these two equations +```math +\zeta = \frac{g}{f}\nabla^2 \eta +``` +Meaning that if we start with ``u,v`` we can obtain the relative vorticity +``\zeta`` and, using Coriolis parameter ``f`` and gravity ``g``, invert +the Laplace operator to obtain displacement ``\eta``. How to do this with +SpeedyTransforms? + +Let us start by generating some data +```@example speedytransforms +using SpeedyWeather + +spectral_grid = SpectralGrid(trunc=31,nlev=1) +forcing = SpeedyWeather.JetStreamForcing(spectral_grid) +drag = QuadraticDrag(spectral_grid) +model = ShallowWaterModel(;spectral_grid,forcing,drag) +model.feedback.verbose = false # hide +simulation = initialize!(model); +run!(simulation,n_days=30) +nothing # hide +``` + +Now pretend you only have `u,v` to get vorticity (which is actually the prognostic variable in the model, +so calculated anyway...). +```@example speedytransforms +u = simulation.diagnostic_variables.layers[1].grid_variables.u_grid +v = simulation.diagnostic_variables.layers[1].grid_variables.v_grid +vor = SpeedyTransforms.curl(u,v) / spectral_grid.radius +nothing # hide +``` +Here, `u,v` are the grid-point velocity fields, and the function `curl` takes in either +`LowerTriangularMatrix`s (no transform needed as all gradient operators act in spectral space), +or, as shown here, arrays of the same grid and size. In this case, the function actually +runs through the following steps +```julia +RingGrids.scale_coslat⁻¹!(u_grid) +RingGrids.scale_coslat⁻¹!(v_grid) + +S = SpectralTransform(u_grid,one_more_degree=true) +us = spectral(u_grid,S) +vs = spectral(v_grid,S) + +return curl(us,vs) +``` +(Copies of) the velocity fields are unscaled by the cosine of latitude (see above), +then transformed into spectral space, and the returned `vor` requires a manual division +by the radius. We always unscale vector fields by the cosine of latitude before any +`curl`, or `div` operation, as these omit those. + +## One more degree for spectral fields + +The `SpectralTransform` in general takes a `one_more_degree` keyword argument, +if otherwise the returned `LowerTriangularMatrix` would be of size 32x32, setting this +to true would return 33x32. The reason is that while most people would expect +square lower triangular matrices for a triangular spectral truncation, all vector quantities +always need one more degree (= one more row) because of a recursion relation in the +meridional gradient. So as we want to take the curl of `us,vs` here, they need this +additional degree, but in the returned lower triangular matrix this row is set to zero. + +!!! info "One more degree for vector quantities" + All gradient operators expect the input lower triangular matrices of shape ``N+1 \times N``. + This one more degree of the spherical harmonics is required for the meridional derivative. + Scalar quantities contain this degree too for size compatibility but they should not + make use of it. Use `spectral_truncation` to add or remove this degree manually. + +## Example: Geostrophy (continued) + +Now we transfer `vor` into grid-point space, but specify that we want it on the grid +that we also used in `spectral_grid`. The Coriolis parameter for a grid like `vor_grid` +is obtained, and we do the following for ``f\zeta/g``. + +```@example speedytransforms +vor_grid = gridded(vor,Grid=spectral_grid.Grid) +f = SpeedyWeather.coriolis(vor_grid) +fζ_g = spectral_grid.Grid(vor_grid .* f ./ model.planet.gravity) +nothing # hide +``` +The last line is a bit awkward for now, as the element-wise multiplication between +two grids escapes the grid and returns a vector that we wrap again into a grid. +We will fix that in future releases. Now we need to apply the inverse +Laplace operator to ``f\zeta/g`` which we do as follows + +```@example speedytransforms +fζ_g_spectral = spectral(fζ_g,one_more_degree=true); +η = SpeedyTransforms.∇⁻²(fζ_g_spectral) * spectral_grid.radius^2 +η_grid = gridded(η,Grid=spectral_grid.Grid) +nothing # hide +``` + +Note the manual scaling with the radius ``R^2`` here. We now compare the results +```@example speedytransforms +plot(η_grid) +``` +Which is the interface displacement assuming geostrophy. The actual interface +displacement contains also ageostrophy +```@example speedytransforms +plot(simulation.diagnostic_variables.surface.pres_grid) +``` +Strikingly similar! The remaining differences are the ageostrophic motions but +also note that the mean is off. This is because geostrophy only use/defines the gradient +of ``\eta`` not the absolute values itself. Our geostrophic ``\eta_g`` has by construction +a mean of zero (that is how we define the inverse Laplace operator) but the actual ``\eta`` +is some 1400m higher. ## Functions and type index