From eaa506189ba5764d607c453e5779e6c62985f18a Mon Sep 17 00:00:00 2001 From: Milan Date: Wed, 20 Mar 2024 13:44:15 -0400 Subject: [PATCH 1/5] Add analysis.md to docs --- docs/make.jl | 1 + docs/src/analysis.md | 233 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 234 insertions(+) create mode 100644 docs/src/analysis.md diff --git a/docs/make.jl b/docs/make.jl index 6736e42de..412fb21e0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,6 +28,7 @@ makedocs( "Running SpeedyWeather" => [ "How to run SpeedyWeather"=>"how_to_run_speedy.md", "Model setups"=>"setups.md", + "Analysis"=>"analysis.md", "Tree structure"=>"structure.md", "Particle advection"=>"particles.md", "NetCDF output"=>"output.md", diff --git a/docs/src/analysis.md b/docs/src/analysis.md new file mode 100644 index 000000000..efd473f5c --- /dev/null +++ b/docs/src/analysis.md @@ -0,0 +1,233 @@ +# Analysing a simulation + +While you can analyze a SpeedyWeather simulation through its [NetCDF output](@ref), +as most users will be used to with other models, you can also reuse a lot of +functionality from SpeedyWeather interactively for analysis. This makes +SpeedyWeather beyond being a model for the atmospheric general circulation +also a library with many functions for the analysis of simulations. +Often this also avoids the two language problem that you will face if you +run a simulation with a model in one language but then do the data +analysis in another, treating the model as a blackbox although it likely +has many of the functions you will need for analysis already defined. +With SpeedyWeather we try to avoid this and are working towards a +more unified approach in atmospheric modelling where simulation +and analysis are done interactively with the same library: SpeedyWeather.jl. + +## Advantages of online analysis + +Now you could run a SpeedyWeather simulation, and analyse the [NetCDF output](@ref) +but that comes with several issues related to accuracy + +- If you use a reduced grid for the simulation, then the output will (by default) be +interpolated on a full grid. This interpolation comes introduces an error. +- Computing integrals over gridded data on the sphere by weighting every grid point +according to its area is not the most accurate numerical integration. +- Computing gradients over gridded data comes with similar issues. While our +[RingGrids](@ref) are always equidistant in longitude, they are not necessarily +in latitude. + +The first point you can avoid by running a simulation on one of the full grids that +are implemented, see [SpectralGrid](@ref). But that also impacts the simulation +and for various reasons we don't run on full grids by default. +The second point you can address by defining a more advanced numerical integration +scheme, but that likely requires you to depend on external libraries and then, +well, you could also just depend on SpeedyWeather.jl directly, because we +have to do these computations internally anyway. Similar for the third point, +many gradients have to be computed on every time step and we do that with +spectral transforms to reduce the discretization error. + +The following contains a (hopefully growing) list of examples +of how a simulation can be analysed interactively. We call this +_online_ analysis because you are directly using the functionality from +the model as if it was a library. + +## Mass conservation + +In the absence of sources and sinks for the interface displacement ``\eta`` +in the [shallow water equations](@ref shallow_water_model), total mass +(or equivalently volume as density is constant) is conserved. +The total volume is defined as the integral of the dynamic layer thickness +``h = \eta + H - H_b`` (``H`` is the layer thickness at rest, ``H_b`` is orography) +over the surface ``A`` of the sphere + +```math +\iint h dA = \iint \eta dA + \iint H dA - \\int H_b dA +``` + +to check for conservation we want to assess that + +```math +\frac{\partial}{\partial t} \\int h dA = 0 +``` + +And because ``V = \iint H dA - \\int H_b dA``, the total volume at rest, +is a constant (``H`` is a global constant, the orography ``H_b`` does not change with time) +we can just check whether ``\\int \eta dA`` changes over time. +Instead of computing this integral in grid-point space, we use the spectral +``\eta`` whereby the coefficient of the first spherical harmonic (the ``l = m = 0`` mode, +see [Spherical Harmonic Transform](@ref) encodes the global average. + +```@example analysis +using SpeedyWeather +spectral_grid = SpectralGrid(trunc=31, nlev=1) +model = ShallowWaterModel(;spectral_grid) +simulation = initialize!(model) +``` + +Now we check ``\eta_{0,0}`` the ``l = m = 0`` coefficent of the inital conditions +of that simulation with + +```@example analysis +simulation.prognostic_variables.surface.timesteps[1].pres[1] +``` + +Its imaginary part is always zero (which is true for any zonal harmonic ``m=0`` as its +imaginary part would just unnecessarily rotate something zonally constant in zonal direction), +so you can `real` it. Also for spherical harmonic transforms there is a norm of the sphere +by which you have to divide to get your mean value in the original units + +```@example analysis +a = model.spectral_transform.norm_sphere # = 2√π = 3.5449078 +η_mean = real(simulation.prognostic_variables.surface.timesteps[1].pres[1]) / a +``` + +So the initial conditions in this simulation are such that the global mean interface displacement +is that value in meters. You would need to multiply by the area of the sphere +``4\pi r^2`` (radius ``r``) to get the actual integral from above, but because that doesn't +change either, we just want to check that `η_mean` doesn't change. Which is equivalent +to ``\partial_t \\iint \eta dA = 0`` and so volume conservation and because density is constant +also mass conservation. Let's check what happens after the simulation ran for some days + +```@example analysis +model.feedback.verbose = false # hide +run!(simulation, period=Day(10)) + +# now we check η_mean again +η_mean_later = real(simulation.prognostic_variables.surface.timesteps[1].pres[1]) / a +``` + +which is _exactly_ the same. So mass is conserved, woohoo. + +Insight from a numerical perspective: The tendency of ``\eta`` is +``\partial_t \eta = -\nabla \cdot (\mathbf{u} h)`` which is a divergence of a flux. +Calculating the divergence in spherical harmonics always sets the ``l=m=0`` mode to zero exactly +(the gradient of a periodic variable has zero mean) so the time integration here is always with +an exactly zero tendency. + +## Energy conservation + +The total energy in the shallow water equation is the sum of the kinetic energy +``\frac{1}{2}(u^2 + v^2)`` and the potential energy ``gh`` integrated over the +total volume (times ``h`` for the vertical then integrated over the sphere``\\int dA``). + +```math +\\iint \frac{h}{2}(u^2 + v^2) + gh^2 dA +``` + +In contrast to the [Mass conservation](@ref) which, with respect to the +spectral transform, is a linear calculation, here we need to multiply +variables, which is done in grid-point space. Then we can transform to spectral +space for the global integral. Let us define a `total_energy` function as + +```@example analysis +using SpeedyWeather +function total_energy(u, v, η, model) + + h = zero(u) + E = zero(u) # allocate grid variable + + H = model.atmosphere.layer_thickness + Hb = model.orography.orography + g = model.planet.gravity + + @. h = η + H - Hb + @. E = h/2*(u^2 + v^2) + g*h^2 + + # transform to spectral, take l=m=0 mode at [1] and normalize for mean + E_mean = real(spectral(E)[1]) / model.spectral_transform.norm_sphere +end +``` + +So at the current state of our simulation we have a total energy +(per square meter as we haven't multiplied by the surface area) + +```@example analysis +# flat copies for convenience +u = simulation.diagnostic_variables.layers[1].grid_variables.u_grid +v = simulation.diagnostic_variables.layers[1].grid_variables.v_grid +η = simulation.diagnostic_variables.surface.pres_grid + +TE = total_energy(u, v, η, model) +``` + +with units of ``m^3 s^{-2}`` (multiplying by surface area of the sphere +and density of the fluid would turn it into joule = ``kg m^2 s^{-2}``). +Now let us continue the simulation +```@example analysis +run!(simulation, period=Day(10)) + +# we don't need to reassign u, v, η as they were flat copies +# pointing directly to the diagnostic variables inside the simulation +# which got updated during run! +TE_later = total_energy(u, v, η, model) +``` +So the total energy has somewhat changed, it decreased to +```@example analysis +TE_later/TE +``` +of its previous value over 10 days. While technically energy should +be conserved in an unforced system, numerically this is rarely +exactly the case. We need some [Horizontal diffusion](@ref diffusion) +for numerical stability and also the time integration is dissipative +due to temporal filtering, see [Time integration](@ref leapfrog). + +## Potential vorticity + +Potential vorticity in the shallow water equations is defined as + +```math +q = \frac{f + \zeta}{h} +``` + +with ``f`` the Coriolis parameter, relative vorticity ``\zeta`` +and ``h`` the layer thickness as before. We can calculate this +conveniently directly on the model grid (whichever you chose) +as + +```@example analysis +# vorticity +ζ = simulation.diagnostic_variables.layers[1].grid_variables.vor_grid +f = coriolis(ζ) # create f on that grid + +# layer thickness +η = simulation.diagnostic_variables.surface.pres_grid +h = zero(η) +H = model.atmosphere.layer_thickness +Hb = model.orography.orography +@. h = η + H - Hb + +# potential vorticity +q = zero(ζ) +@. q = (f + ζ) / h +``` + +and we can compare the relative vorticity field to +```@example analysis +plot(ζ) +``` +the potential vorticity +```@example analysis +plot(q) +``` + +## Absolute angular momentum + +More to follow ... + +## Circulation + +More to follow ... + +## Enstrophy + +More to follow ... \ No newline at end of file From d3add65efd45ac8f9c6ba2e1e1815223cf6d23a6 Mon Sep 17 00:00:00 2001 From: Milan Date: Wed, 20 Mar 2024 14:11:37 -0400 Subject: [PATCH 2/5] formatting issues in analysis.md --- docs/src/analysis.md | 17 +++++++++++------ docs/src/barotropic.md | 32 +++++++++++++++++++------------- 2 files changed, 30 insertions(+), 19 deletions(-) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index efd473f5c..eac57fc80 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -20,8 +20,10 @@ but that comes with several issues related to accuracy - If you use a reduced grid for the simulation, then the output will (by default) be interpolated on a full grid. This interpolation comes introduces an error. + - Computing integrals over gridded data on the sphere by weighting every grid point according to its area is not the most accurate numerical integration. + - Computing gradients over gridded data comes with similar issues. While our [RingGrids](@ref) are always equidistant in longitude, they are not necessarily in latitude. @@ -51,21 +53,21 @@ The total volume is defined as the integral of the dynamic layer thickness over the surface ``A`` of the sphere ```math -\iint h dA = \iint \eta dA + \iint H dA - \\int H_b dA +\iint h dA = \iint \eta dA + \iint H dA - \iint H_b dA ``` to check for conservation we want to assess that ```math -\frac{\partial}{\partial t} \\int h dA = 0 +\frac{\partial}{\partial t} \iint h dA = 0 ``` -And because ``V = \iint H dA - \\int H_b dA``, the total volume at rest, +And because ``V = \iint H dA - \iint H_b dA``, the total volume at rest, is a constant (``H`` is a global constant, the orography ``H_b`` does not change with time) -we can just check whether ``\\int \eta dA`` changes over time. +we can just check whether ``\iint \eta dA`` changes over time. Instead of computing this integral in grid-point space, we use the spectral ``\eta`` whereby the coefficient of the first spherical harmonic (the ``l = m = 0`` mode, -see [Spherical Harmonic Transform](@ref) encodes the global average. +or wavenumber 0, see [Spherical Harmonic Transform](@ref)) encodes the global average. ```@example analysis using SpeedyWeather @@ -81,6 +83,8 @@ of that simulation with simulation.prognostic_variables.surface.timesteps[1].pres[1] ``` +`[1]` pulls the first element of the underlying [LowerTriangularMatrix](@ref lowertriangularmatrices) +which is the coefficient of the ``l = m = 0`` mode. Its imaginary part is always zero (which is true for any zonal harmonic ``m=0`` as its imaginary part would just unnecessarily rotate something zonally constant in zonal direction), so you can `real` it. Also for spherical harmonic transforms there is a norm of the sphere @@ -118,7 +122,7 @@ an exactly zero tendency. The total energy in the shallow water equation is the sum of the kinetic energy ``\frac{1}{2}(u^2 + v^2)`` and the potential energy ``gh`` integrated over the -total volume (times ``h`` for the vertical then integrated over the sphere``\\int dA``). +total volume (times ``h`` for the vertical then integrated over the sphere``\iint dA``). ```math \\iint \frac{h}{2}(u^2 + v^2) + gh^2 dA @@ -209,6 +213,7 @@ Hb = model.orography.orography # potential vorticity q = zero(ζ) @. q = (f + ζ) / h +nothing # hide ``` and we can compare the relative vorticity field to diff --git a/docs/src/barotropic.md b/docs/src/barotropic.md index 2ebc7db1e..4dd27681a 100644 --- a/docs/src/barotropic.md +++ b/docs/src/barotropic.md @@ -9,9 +9,10 @@ The dynamical core presented here to solve the barotropic vorticity equations la the idealized models with spectral dynamics developed at the Geophysical Fluid Dynamics Laboratory[^1]: A barotropic vorticity model[^2]. -Many concepts of the [Shallow water model](@ref shallow_water_model) and the [Primitive equation model](@ref primitive_equation_model) are similar, -such that for example [horizontal diffusion](@ref diffusion) and the [Time integration](@ref leapfrog) -are only explained here. +Many concepts of the [Shallow water model](@ref shallow_water_model) and the +[Primitive equation model](@ref primitive_equation_model) are similar, +such that for example [horizontal diffusion](@ref diffusion) and the +[Time integration](@ref leapfrog) are only explained here. ## Barotropic vorticity equation @@ -62,8 +63,8 @@ in the `BarotropicModel`, as outlined in the following section. ## Algorithm -We briefly outline the algorithm that SpeedyWeather.jl uses in order to integrate the barotropic vorticity equation. -As an initial step +We briefly outline the algorithm that SpeedyWeather.jl uses in order to integrate the barotropic +vorticity equation. As an initial step 0\. Start with initial conditions of ``\zeta_{lm}`` in spectral space and transform this model state to grid-point space: @@ -91,18 +92,22 @@ Now loop over ## [Horizontal diffusion](@id diffusion) -In SpeedyWeather.jl we use hyperdiffusion through an ``n``-th power Laplacian ``(-1)^{n+1}\nabla^{2n}`` (hyper when ``n>1``) which -can be implemented as a multiplication of the spectral coefficients ``\Psi_{lm}`` with ``(-l(l+1))^nR^{-2n}`` (see -spectral [Laplacian](@ref)) It is therefore computationally not more expensive to apply hyperdiffusion over diffusion -as the ``(-l(l+1))^nR^{-2n}`` can be precomputed. Note the sign change ``(-1)^{n+1}`` here is such that the dissipative nature of the diffusion operator is retained for ``n`` odd and even. +In SpeedyWeather.jl we use hyperdiffusion through an ``n``-th power Laplacian ``(-1)^{n+1}\nabla^{2n}`` +(hyper when ``n>1``) which +can be implemented as a multiplication of the spectral coefficients ``\Psi_{lm}`` with +``(-l(l+1))^nR^{-2n}`` (see spectral [Laplacian](@ref)). It is therefore computationally not more +expensive to apply hyperdiffusion over diffusion as the ``(-l(l+1))^nR^{-2n}`` can be precomputed. +Note the sign change ``(-1)^{n+1}`` here is such that the dissipative nature of the diffusion operator +is retained for ``n`` odd and even. In SpeedyWeather.jl the diffusion is applied _implicitly_. For that, consider a -[leapfrog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration) with time step ``\Delta t`` of variable ``\zeta`` to -obtain from time steps ``i-1`` and ``i``, the next time step ``i+1`` +[leapfrog scheme](https://en.wikipedia.org/wiki/Leapfrog_integration) with time step ``\Delta t`` of +variable ``\zeta`` to obtain from time steps ``i-1`` and ``i``, the next time step ``i+1`` ```math \zeta_{i+1} = \zeta_{i-1} + 2\Delta t d\zeta, ``` -with ``d\zeta`` being some tendency evaluated from ``\zeta_i``. Now we want to add a diffusion term ``(-1)^{n+1}\nu \nabla^{2n}\zeta`` +with ``d\zeta`` being some tendency evaluated from ``\zeta_i``. Now we want to add a diffusion term +``(-1)^{n+1}\nu \nabla^{2n}\zeta`` with coefficient ``\nu``, which however, is implicitly calculated from ``\zeta_{i+1}``, then ```math @@ -129,7 +134,8 @@ only an element-wise multiplication in spectral space. For every spectral harmon d\zeta \to D_\text{implicit}^{-1}(d\zeta + D_\text{explicit}\zeta_{i-1}). ``` Hence 2 multiplications and 1 subtraction with precomputed constants. -However, we will normalize the (hyper-)Laplacians as described in the following. This also will take care of the alternating sign such that the diffusion operation is dissipative regardless the power ``n``. +However, we will normalize the (hyper-)Laplacians as described in the following. This also will take care of +the alternating sign such that the diffusion operation is dissipative regardless the power ``n``. ### Normalization of diffusion From 083a86959fdec4ccc1067a83aff4760f2e623e24 Mon Sep 17 00:00:00 2001 From: Milan Date: Wed, 20 Mar 2024 14:15:17 -0400 Subject: [PATCH 3/5] More formatting issues in analysis.md --- docs/src/analysis.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index eac57fc80..2533d81f7 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -98,8 +98,8 @@ a = model.spectral_transform.norm_sphere # = 2√π = 3.5449078 So the initial conditions in this simulation are such that the global mean interface displacement is that value in meters. You would need to multiply by the area of the sphere ``4\pi r^2`` (radius ``r``) to get the actual integral from above, but because that doesn't -change either, we just want to check that `η_mean` doesn't change. Which is equivalent -to ``\partial_t \\iint \eta dA = 0`` and so volume conservation and because density is constant +change with time either, we just want to check that `η_mean` doesn't change with time. +Which is equivalent to ``\partial_t \iint \eta dA = 0`` and so volume conservation and because density is constant also mass conservation. Let's check what happens after the simulation ran for some days ```@example analysis @@ -125,7 +125,7 @@ The total energy in the shallow water equation is the sum of the kinetic energy total volume (times ``h`` for the vertical then integrated over the sphere``\iint dA``). ```math -\\iint \frac{h}{2}(u^2 + v^2) + gh^2 dA +\iint \frac{h}{2}(u^2 + v^2) + gh^2 dA ``` In contrast to the [Mass conservation](@ref) which, with respect to the From 349cb6c6f9789ab0839c6c50b1da3d2b8363551f Mon Sep 17 00:00:00 2001 From: Milan Date: Wed, 20 Mar 2024 14:21:13 -0400 Subject: [PATCH 4/5] Typos in land_sea_mask.md --- docs/src/land_sea_mask.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/src/land_sea_mask.md b/docs/src/land_sea_mask.md index ceded14ef..69b2d1a21 100644 --- a/docs/src/land_sea_mask.md +++ b/docs/src/land_sea_mask.md @@ -43,7 +43,7 @@ land_sea_mask = LandSeaMask(spectral_grid) which will automatically interpolate the land-sea mask onto grid and resolution as defined in `spectral_grid` at initialization. The actual mask is in -`land_sea_mask.land_sea_mask` and you can visualise it with +`land_sea_mask.mask` and you can visualise it with ```@example landseamask model = PrimitiveWetModel(;spectral_grid, land_sea_mask) @@ -180,4 +180,6 @@ plot(model.land_sea_mask.mask) ``` And the land-sea mask has succesfully been set to ocean everywhere at the start -of the 21st century. +of the 21st century. Note that while we added an `@info` line into the +`callback!` function, this is here not printed because of how the +Documenter works. If you execute this in the REPL you'll see it. From 6d9974189c52996cd2cda0838d143bb806604adc Mon Sep 17 00:00:00 2001 From: Milan Date: Wed, 20 Mar 2024 14:47:42 -0400 Subject: [PATCH 5/5] analysis.md list reformatted --- docs/src/analysis.md | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/docs/src/analysis.md b/docs/src/analysis.md index 2533d81f7..002da03c1 100644 --- a/docs/src/analysis.md +++ b/docs/src/analysis.md @@ -18,15 +18,9 @@ and analysis are done interactively with the same library: SpeedyWeather.jl. Now you could run a SpeedyWeather simulation, and analyse the [NetCDF output](@ref) but that comes with several issues related to accuracy -- If you use a reduced grid for the simulation, then the output will (by default) be -interpolated on a full grid. This interpolation comes introduces an error. - -- Computing integrals over gridded data on the sphere by weighting every grid point -according to its area is not the most accurate numerical integration. - -- Computing gradients over gridded data comes with similar issues. While our -[RingGrids](@ref) are always equidistant in longitude, they are not necessarily -in latitude. +- If you use a reduced grid for the simulation, then the output will (by default) be interpolated on a full grid. This interpolation comes introduces an error. +- Computing integrals over gridded data on the sphere by weighting every grid point according to its area is not the most accurate numerical integration. +- Computing gradients over gridded data comes with similar issues. While our [RingGrids](@ref) are always equidistant in longitude, they are not necessarily in latitude. The first point you can avoid by running a simulation on one of the full grids that are implemented, see [SpectralGrid](@ref). But that also impacts the simulation