Conserved quantity verifying

[mini-DONG]

Hi, thanks a lot for your help at #391 . Now, I have run a shallow water model like:

spectral_grid = SpectralGrid(trunc=511,nlev=1)
initial_conditions = StartFromRest()
output = OutputWriter(spectral_grid,ShallowWater,
                      output_vars=[:u,:v,:pres,:vor,:div])
model = ShallowWaterModel(;spectral_grid,initial_conditions,output=output)
simulation = initialize!(model)

with initial conditions:

simulation.prognostic_variables.layers[1].timesteps[1].vor .= vor/spectral_grid.radius
simulation.prognostic_variables.layers[1].timesteps[1].div .= div/spectral_grid.radius
simulation.prognostic_variables.surface.timesteps[1].pres .= etas

I have run 1 year, and the first output (i.e. initial conditions) and last output is like: Like #323, I also want to verify the conserved quantity. I calculate the total thickness with python codes:

h = ds.pres + 8500
cos_phi =  np.cos(np.radians(ds.lat))
H = (h*cos_phi).sum(['lat','lon'])

The plot of H: H is growing fast, which doesn't seem reasonable. Or maybe I should use other way to calculate?

[milankl]

If I read your y axis correctly then you have

julia> h1 = 6200+6.3875e9
6.3875062e9

julia> h0 = 5500+6.3875e9
6.3875055e9

julia> h1/h0
1.0000001095889468

It looks like it's growing "fast" but actually it's conserved to 1e-7 which I find very reasonable given there are still several sources of error

• It's single precision so more than 1e-8 you cannot expect anyway
• By default you are using the OctahedralGaussianGrid, but the output will be on the corresponding full grid, that's the FullGaussianGrid. The default interpolation is not conservative.
• you sum the array, but on the FullGaussianGrid even with coslat scaling the grid cells are not of even size because of the Gaussian latitudes. You use a high resolution T511, but this error may still be not negligible.

Do you just want to be sure that volume is conserved? Then that's :point_up: probably good enough, if you want to try and close the mass budget as best as possible I'd start by those points.

[mini-DONG]

Oh! I see! Thank you very much! I made a mistake, because I just look at the trend. If I want to

try and close the mass budget as best as possible

How can I adjust, especially with the grid? maybe I should the even size lon-lat grid.

In addition, at #391, when making the initial conditions:

Get your initial conditions onto one of the grids of our supported grids, e.g. julia> u = randn(FullGaussianGrid,24) julia> v = randn(FullGaussianGrid,24)

How can I determine this number 24 which is the size of the array. I think this have something with the trunc:

spectral_grid = SpectralGrid(trunc=31,nlev=1)

So, how can I decide the array size by trunc-number? In other words, how do I know the the array size of initial conditions in diffrernt trunc-number?

[mini-DONG]

I am interested in the energy budget, so I also want to know the default viscosity coefficient and how to modify it. I read the output's parameters.txt which seem don't mention the viscosity.

[mini-DONG]

For the other conserved quantity, such as momentum, angular momentum, I calculated them with python codes:

h = ds.pres + 8500
cos_phi =  np.cos(np.radians(ds.lat))
momentum =(h*cos_phi*ds.u).sum(['lat','lon'])
ang_m = (h*(cos_phi**2)*ds.u).sum(['lat','lon'])

The plot is: Is it reasonable to lose momentum all the time?

Moreover, I find cycle in the parameter.txt:

daily_cycle::Bool = true; seasonal_cycle::Bool = true;

Is this the default setting? Does this mean that the model has external forcing? And how to close it?

[milankl]

What is the mathematical formula of the conservations you want to use? Once you multiply h with u I'd definitely include orography. And what happened to hv?

There's no forcing by default, you can check that with model.forcing.

I am interested in the energy budget, so I also want to know the default viscosity coefficient and how to modify it. I read the output's parameters.txt which seem don't mention the viscosity.

There's no viscosity in the simulation. There's hyper diffusion which removes energy at the smallest scales to keep to model stable. Read the documentation https://speedyweather.github.io/SpeedyWeather.jl/dev/barotropic/#diffusion for more information. While the operator is dimensionless, the associated time scale is by default

julia> model.horizontal_diffusion.time_scale
2.4

hours. You can translate that back into a diffusion coefficient. But note that by default we use a power 4 Laplacian so it'll have units of m^8/s!

[milankl]

So, how can I decide the array size by trunc-number? In other words, how do I know the the array size of initial conditions in diffrernt trunc-number?

It's all in the documentation!

• https://speedyweather.github.io/SpeedyWeather.jl/dev/spectral_transform/#Available-horizontal-resolutions
• https://speedyweather.github.io/SpeedyWeather.jl/dev/grids/#Matching-spectral-and-grid-resolution

[mini-DONG]

I use the following formula to calculate momentum and angular momentum:

$$ m = \rho_0 \iint uh dxdy = \rho_0 *r^2 \iint uh \cos(\phi) d\lambda d\phi $$

$$ \Lambda = \rho_0 \iint uhr\cos(\phi) dxdy = \rho_0 *r^3 \iint uh \cos^2(\phi) d\lambda d\phi $$

with python codes:

momentum =(h*cos_phi*ds.u).sum(['lat','lon'])
ang_m = (h*(cos_phi**2)*ds.u).sum(['lat','lon'])

In the codes, I don't multiply the $\rho_0*r^n$, because it is constant. I only use the $u$, so I think I calculate the zonal momentum. As for meridional momentum $hv$, I also calculate:

momentum_v =(h*cos_phi*ds.v).sum(['lat','lon'])

with the plot: The meridional momentum is around zero, so I only focus on zonal momentum. And it seems that my model does not include orography?

[milankl]

Orography is included by default in the shallow water model, which you can check with model.orography

[milankl]

I use the following formula to calculate momentum and angular momentum:

m=ρ0∬uhdxdy=ρ0∗r2∬uhcos⁡(ϕ)dλdϕ

Λ=ρ0∬uhrcos⁡(ϕ)dxdy=ρ0∗r3∬uhcos2⁡(ϕ)dλdϕ

Can you provide a reference for this? I need to check whether the assumptions in that derivation are also given for our ShallowWaterModel. Also please include orography $H_o$ in $h = \eta + H - H_o$, while you can leave that out in the conservation of volume case, you'll need it here due to the product. Just add :orography to the output_vars argument when constructing the OutputWriter,i.e. OutputWriter(;spectral_grid, output_vars=[:u,:v,:pres,:orography]

[mini-DONG]

I checked out the formula and I found that I made a mistake. Only absolute angular momentum is conserved on spherical coordinates. I run model again with setting no-orography and FullClenshawGrid:

spectral_grid = SpectralGrid(trunc=511,nlev=1,Grid=FullClenshawGrid)
initial_conditions = StartFromRest()
orography = NoOrography(spectral_grid)
output = OutputWriter(spectral_grid,ShallowWater,   
                      output_vars=[:u,:v,:pres,:vor,:div,:orography])
model = ShallowWaterModel(;spectral_grid,initial_conditions,orography,output=output)
simulation = initialize!(model)

And than, I calculate absolute angular momentum $m = ur + \omega r^2$ , where $r=a \cos \phi$, $a$ is earth radius. Global integral:

$$ \Lambda = \iint \left( ur+\omega r^2 \right) dxdy = a^3 \iint h \left(u \cos^2 \phi + \omega a \cos^3 \phi \right) d\lambda d\phi $$

with python codes:

a =  6.371e6
cos_phi =  np.cos(np.radians(ds.lat))
omega = (2*np.pi)/(24*3600)

ang_m = (h*(u*cos_phi**2+omega*a*cos_phi*3)).sum(['lat','lon'])

The plotting: Bottom friction affects angular momentum, I want to know the default bottom friction and how to setting. Besides, I want to change hyper diffusion to Laplace viscosity.

[milankl]

Quick note, I'd recommend to use directly the values from the model

julia> model.planet.rotation
7.29e-5

julia> model.spectral_grid.radius
6.371e6

The omega you use is slightly off because the planet actually rotates in 24 hours more than around itself. This is because in 24 hours it also travels 1/365.25 around the sun. Check https://en.wikipedia.org/wiki/Sidereal_time Normally I wouldn't bother but because you are trying to close some budget you should use what's actually used inside the model. Bottom friction/drag by default is

julia> model.drag
SpeedyWeather.NoDrag{Float32} <: AbstractDrag

which you should always check as default setups may change from version to version.

[milankl]

For a normal Laplacian do power=1 (default is 4)

julia> horizontal_diffusion = HyperDiffusion(spectral_grid,power=1,time_scale=24)
HyperDiffusion{Float32}
├ trunc::Int64 = 31
├ nlev::Int64 = 8
├ power::Float64 = 1.0
├ time_scale::Float64 = 24.0
├ resolution_scaling::Float64 = 0.5
├ power_stratosphere::Float64 = 2.0
├ tapering_σ::Float64 = 0.2
├ adaptive::Bool = true
├ vor_max::Float64 = 0.0001
├ adaptive_strength::Float64 = 2.0
└── arrays: ∇²ⁿ_2D, ∇²ⁿ_2D_implicit, ∇²ⁿ, ∇²ⁿ_implicit

julia> model = ShallowWaterModel(;spectral_grid,horizontal_diffusion)

You will also need to increase the time_scale from 2.4 hours (default for power 4) to something larger because otherwise the model is crazy diffusive. But you'll need to play around with it, I don't know what to use, I just picked 24 hours out of the hat here. Spectral models typically don't use a simple Laplacian because it's too diffusive (and there's no additional cost in using higher orders).

[mini-DONG]

Oh, thanks, I will use the the values from the model directly. And I check the bottom drag:

model.drag
SpeedyWeather.NoDrag{Float32} <: AbstractDrag

I think this means no drag. So why is the absolute angular momentum decreasing in the first two months?

[miniufo]

You will also need to increase the time_scale from 2.4 hours (default for power 4) to something larger because otherwise the model is crazy diffusive. But you'll need to play around with it, I don't know what to use, I just picked 24 hours out of the hat here. Spectral models typically don't use a simple Laplacian because it's too diffusive (and there's no additional cost in using higher orders).

Hi @milankl, I am not very familiar with spectral method but interested here. Standard viscous term, $\nu \nabla^2 u$ where $\nu$ is the viscosity, does not have a explicit timescale. So where does the timescale come from? Also there are resolution_scaling, power_stratosphere, tapering_σ etc which are not explicit in the expression. Could you please point me to the references to these?

I and @mini-DONG just want to close the budget as best as we could, so we now choose the simplest form of viscous dissipation. Another thing is that if we choose to output at daily frequency, does the default model output is instantansous or daily-averaged one? I guess this may also affect the budget.

[milankl]

On time scale, see https://speedyweather.github.io/SpeedyWeather.jl/dev/barotropic/#Normalization-of-diffusion Our Laplace operators (regardless which power) are normalized by the largest eigenvalue $l{max}(l{max} + 1)$ so for trunc=511 as you have used before, it is normalised by 511*512. The diffusion coefficient becomes with that scaling a time scale, and we only use that time scale (next to the power) the control the strength of the diffusion.

power_stratosphere and tapering_sigma are irrelevant for the ShallowWaterModel because there's only one level. You can ignore those as they are somewhat experimental features I have been playing around with, hence they are not documented. resolution_scaling currently only applies for the primitive equation model, so you can forget about that too.

[milankl]

I and @mini-DONG just want to close the budget as best as we could, so we now choose the simplest form of viscous dissipation.

Yes, but it doesn't really make a difference whether you use time scale times power 1,2,3 or 4 Laplacian on your side, does it? I don't know how you want to calculate that term, but to be consistent with the model you'll need to transform into spectral space and apply it there anyway, in which case you also run into the problem that the diffusion inside the model is applied implicitly, check https://speedyweather.github.io/SpeedyWeather.jl/dev/barotropic/#diffusion

Another thing is that if we choose to output at daily frequency, does the default model output is instantansous or daily-averaged one? I guess this may also affect the budget.

Output is always instantaneous, the only exception is precipiation which is accumulated - but you do not use that anyway. Non-instantaneous values would give you a hard time anyway as any non-linear metric wouldn't close I guess: For example in your formulation you have h times u, which is not the same as h averaged over some time times u averaged over some time.

[milankl]

So why is the absolute angular momentum decreasing in the first two months?

I don't know. For these kind of budgets I would always try to replicate as best as possible what's done in the model. With angular momentum I do not have enough experience to say what's affecting that budget and what is not. We would need to play around with it. There's many numerics inside the model that make this a non-trivial exercise

• the transform comes with an error, choosing an exact grid and the dealiasing factor sufficiently high can minimize this, search for 'exactness' in https://speedyweather.github.io/SpeedyWeather.jl/dev/grids/ (we currently do not have a specific section on this)
• the time integration is not necessarily conserving because of time filters we need to use, so you could check how much changes when you reduce the time step
• the semi-implicit integration has its tweaks that may be important on how you calculate budgets
• ...

[mini-DONG]

Hi, I ran a shallow-water simulation earlier, and now I want to add forcing on that. If I think of the shallow-water as an ocean, and I want to put wind stress on the upper surface, what should I do?

[milankl]

Hey @mini-DONG, sorry I never responded to that. We have a JetStreamForcing that acts like a windstress on top of your shallow water layer

help?> JetStreamForcing
search: JetStreamForcing

  Forcing term for the Barotropic or ShallowWaterModel with an idealised jet
  stream similar to the initial conditions from Galewsky, 2004, but mirrored
  for both hemispheres.

    •  nlat::Int64: Number of latitude rings

    •  latitude::Any: jet latitude [˚N]

    •  width::Any: jet width [˚], default ≈ 19.29˚

    •  speed::Any: jet speed scale [m/s]

    •  time_scale::Second: time scale [days]

    •  amplitude::Vector: precomputed amplitude vector [m/s²]

julia> forcing = JetStreamForcing(spectral_grid)
JetStreamForcing{Float32} <: Main.SpeedyWeather.AbstractForcing
├ nlat::Int64 = 256
├ latitude::Float32 = 45.0
├ width::Float32 = 19.285715
├ speed::Float32 = 85.0
├ time_scale::Second = 2592000 seconds
└── arrays: amplitude

This part in the documentation gives you an idea of how to set this up https://speedyweather.github.io/SpeedyWeather.jl/stable/setups/#Polar-jet-streams-in-shallow-water

[miniufo]

Hi @milankl, I just move back to this issue, where previously I have tried some schemes to see several conserved quantities like:

• total mass $M=\iint h dA$;
• total absolute angular momentum $\Lambda = \iint ur + \Omega r^2 dA$ ($r=a \cos\phi$);
• total circulation $C = \iint q h dA$ (should be zero);
• total enstrophy $\iint q^2/2 h dA$, and
• total energy $\iint h\mathbf{u}^2/2 + gh^2/2dA$.

I and @mini-DONG found that the evaluation of these integrations depend on the grids and schemes. So I would like to know:

1. Is there any intrinsic way to do these integration or diagnostic-calculation from spectral space, which is expected to be most accurate as some kind of reference;
2. If we 'have to' evaluate all these diagnostics on a certain grid, which one is the 'best' one (close to those reference obtained in spectral space)?

We need to close the energy budget using speedy simulation. So we want to know how well speedy is doing this (not perfect for every metric). Considering previous discussions at here, there could be several other things affecting the accuracy. Please remind us if you feel it is necessary.

[milankl]

On 1: The linear quantities you can calculate in spectral space, e.g. total mass. Linear here means that any variable that depends on space can be added but not multiplied with another variable that also depends on space. E.g. $u$, $cu$ ($c$ is a global constant like gravity), $u+v$ are fine, but $u$ is not (element-wise multiplication). This includes things like $u\cos(\theta)$ ($\theta$ being latitude). In your case for total mass you'd have

$$ \iint h dA = \iint \eta + H - H_b dA$$

where layer thickness $H$ is time-independent, as it the orography $H_b$. So in order to verify conservation of total mass you can just check that the integral of $\eta$ doesn't change over time. Evey spectral field represented as LowerTriangularMatrix has its mean stored in the first value, i.e. the $l=m=0$ mode. So the global mean interface displacement $\eta$ of default initial conditions is

spectral_grid = SpectralGrid(trunc=31, nlev=1)
model = ShallowWaterModel(;spectral_grid)
simulation = initialize!(model)

Now the l=m=0 of $\eta$ is in pres[1]

julia> simulation.prognostic_variables.surface.timesteps[1].pres[1]
4634.5415f0 + 0.0f0im

It's imaginary part is always zero, 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

julia> a = model.spectral_transform.norm_sphere
3.5449078f0

julia> real(simulation.prognostic_variables.surface.timesteps[1].pres[1]) / a
1307.38f0

So with the chosen initial conditions $\int \eta dA = \bar{\eta}$ is 1307 m times the surface of the sphere. Now run the simulation and check whether this changes

julia> run!(simulation, period=Day(10));
Weather is speedy: 100%|████████████████████████████████████████████| Time: 0:00:00 (12.73 millenia/day)

julia> real(simulation.prognostic_variables.surface.timesteps[1].pres[1]) / a
1307.38f0

And in fact it cannot, because $\partial_t \eta = -\nabla \cdot (\mathbf{u} h)$ is a divergence of a flux with zero mean and calculating the divergence in spherical harmonics always sets the $l=m=0$ to zero exactly, so timestepping here is always with a zero tendency.

[milankl]

On 2: The Gaussian and Clenshaw-Curtis grids are exact, meaning that the transform has only floating-point errors but no transform error. For Gaussian grids that happens at a dealiasing factor of $>=2$, with the Clenshaw-Curtis grids it's $>=3$. So I'd suggest you do use

spectral_grid = SpectralGrid(Grid=FullClenshawGrid, dealiasing=3)

The FullClenshawGrid is a regular longitude-latitude grid, which likely makes your analysis easier. With Gaussian grids, the latitudes are not equally spaced because they are Gaussian latitudes, meaning you'll get another complication when calculating global integrals. Instead of calculating the global integral by hand you can also just transform a grid back into spectral space and then use the $l=m=0$ as shown above for the global average. You find more information on how to do this here https://speedyweather.github.io/SpeedyWeather.jl/dev/speedytransforms/

So for $\iint \frac{1}{2} (u^2 + v^2) dA$ you would do $u^2, v^2$ in grid-point space as element-wise multiplication then transform the result to spectral space and check again the $l=m=0$ mode as shown above.

julia> u = simulation.diagnostic_variables.layers[1].grid_variables.u_grid
julia> v = simulation.diagnostic_variables.layers[1].grid_variables.v_grid
julia> KE = zero(u)               # allocate grid of same size
julia> @. KE = 1/2 * (u^2 + v^2)  # with @. everything element-wise
julia> real(spectral(KE)[1])/a    # to spectral space, get l=m=0 and normalise
160.22296f0

[milankl]

Wrapped into a function

julia> function total_energy(u,v,η,model)
           E = zero(u)
           H = model.atmosphere.layer_thickness
           Hb = model.orography.orography
           g = model.planet.gravity
           @. E = 1/2 * ((η + H - Hb)*(u^2 + v^2) + g*(η + H - Hb)^2)
           real(spectral(E)[1])/model.spectral_transform.norm_sphere
       end

You'd have initially

julia> simulation = initialize!(model)
julia> run!(simulation, period=Day(0));    # run for 0 days to propagate spectral initial conditions to grid

julia> u = simulation.diagnostic_variables.layers[1].grid_variables.u_grid;   # flat copy
julia> v = simulation.diagnostic_variables.layers[1].grid_variables.v_grid;
julia> η = simulation.diagnostic_variables.surface.pres_grid;

julia> total_energy(u, v, η, model)
4.5406352f8

then after 100days

julia> run!(simulation, period=Day(100));
Weather is speedy: 100%|████████████████████████████████████████████| Time: 0:00:01 (14.91 millenia/day)

julia> total_energy(u, v, η, model)   # u,v, eta still point to the same mutable array
4.5268032f8

The model has lost about 0.3% in total energy. I find that acceptable, and I wasn't expecting anything more than that to be honest.

[milankl]

Just repeated that over 1000 days, I get

[miniufo]

That's awesome! I am trying now but a tech problem is how to calculate absolute vorticiy. I have:

• relative vorticity $\zeta$ = simulation.diagnostic_variables.layers[1].grid_variables.vor_grid (2D variable)
• Coriolis parameter $f$ = model.coriolis.f (1D variable)

it is not possible to sum them up $\zeta_a = \zeta + f$ as their dimensions are not the same. I've already got used to python-xarray's broadcast capability. So is there an easy way to do the summation here for two arrays of different shapes?

[milankl]

By default we use reduced grids that cannot be represented as a matrix, hence also numpy's broadcasting (Julia's broadcasting is similar) wouldn't help you. You can use a full grid instead or interpolate (as the output always is) but if you want to do it on the reduced grid avoiding other errors you'd do

ζ = simulation.diagnostic_variables.layers[1].grid_variables.vor_grid
ζ ./= model.spectral_grid.radius    # unscale radius-scaling
f = SpeedyWeather.coriolis(ζ)     # create f on that grid

η = simulation.diagnostic_variables.surface.pres_grid
h = zero(η)
H = model.atmosphere.layer_thickness
Hb = model.orography.orography
@. h = η + H - Hb

q = zero(ζ)
@. q = (f + ζ) / h

• By default, accessing vorticity or divergence in the diagnostic variables will be scaled by the radius, I realise now that for these use cases we should unscale them when the simulation is over. I'll create a pull request and tag you in it once that's done. For now, do ./= model.spectral_grid.radius to unscale.
• you see I always create a new variable like h = zero(η) and then work with in-place operations through @., this guarantees that your variable stays a ::AbstractGrid and doesn't become a vector (which is just the underlying vector but therefore you lose the information what grid you were originally on
• Given a grid::AbstractGrid you can do f = SpeedyWeather.coriolis(ζ) to get the Coriolis parameter on that grid. We should export this function, we currently don't, sorry.

Note also that if you stay within Julia and work directly with our grids you can always do plot(grid) to get a quick unicodeplot of a given variable grid, e.g.

This works because the underlying vector is wrapped into an object (here OctahedralGaussianGrid for which geographical locations of the vector elements are unambiguously defined. Maybe that's helpful.

[miniufo]

Great, that works for me.

By default we use reduced grids that cannot be represented as a matrix, hence also numpy's broadcasting (Julia's broadcasting is similar) wouldn't help you. You can use a full grid instead or interpolate (as the output always is)

What's the difference between reduced and full grids? Does the output always use interpolation so that there are unavoidable errors (e.g., cannot get exact total energy as in the spectral grid) even the Gaussian and Clenshaw-Curtis grids you suggested?

I also found that the outputs do not have grid points at poles (+90 or -90). Does that affect the calculation of total (potential) energy?

[milankl]

Read on grids here https://speedyweather.github.io/SpeedyWeather.jl/dev/grids/#Implemented-grids

Output is by default on the corresponding full grid of whatever reduced grid the simulation uses and using the same resolution parameter nlat_half. You can check this with

julia> RingGrids.full_grid(OctahedralGaussianGrid)
FullGaussianGrid

So when using reduced grids, the output is interpolated, for full grids it's not (unless you specify a different resolution). Read here https://speedyweather.github.io/SpeedyWeather.jl/dev/output/#Example-2:-Output-onto-a-higher/lower-resolution-grid

No grid in SpeedyWeather (and in any other spectral model afaik) has a grid point on the pole. It's not needed for an exact transform and, in fact, its FFT would be just the identity. If you do need a point on the pole you can average the pole-most ring, e.g.

# create a random grid of given type (octahedral Gaussian) and resolution parameter nlat_half = 4
grid = rand(OctahedralGaussianGrid, 4)
rings = RingGrids.eachring(grid)    # index range per ring
first_ring = rings[1]    # ring around north pole
last_ring = rings[end]  # ring around south pole

using Statistics
north_pole = mean(grid[first_ring])
south_pole = mean(grid[last_ring])

for full grids you may intuitively think of doing mean(grid[1,:]) (or [end,:]) but this doesn't work with reduced grids. The above workflow works for all grids.

To answer your question: If you calculate the global integral via the spectral transform as I did it, no it doesn't matter that we don't have a point on the pole. If you calculate the integral based on grid points by simply adding them up with area weights then you'll have a much larger error from numerical integration anyway (because that's not a good quadrature rule https://en.wikipedia.org/wiki/Numerical_integration). For the Gaussian grids we use the Gaussian quadrature (hence the name) https://en.wikipedia.org/wiki/Gaussian_quadrature which is exact (rounding errors only) if you use at leas the default dealiasing = 2. Hence, I generally recommend using the spectral transform for this kind of stuff as it's more precise and also how the model does it. I know it's probably less familiar to people that work a lot with gridded model output, but that's why I'm grateful for you to raise these issues so that we can work together for a better documentation of these methods and interfaces.

On that note, I'm imagining a Analyzing a simulation section in the documentation where we could document things like

• Calculating the total mass
• Calculating the potential vorticity
• Calculating total energy

as we discussed them here. If you'd like to contribute I'm happy to co-create a pull request to the documentation based on everything you've learned through this issue. We have already a similar section on Geostrophy which shows how to use the gradient operators directly from the model on the model grids https://speedyweather.github.io/SpeedyWeather.jl/dev/speedytransforms/#Example:-Geostrophy

[miniufo]

Great. I just learn a lot on grids here. So even if I prescribe the full grid for SpectralGrid at the beginning, the model in the internal will still use its corresponding reduced grid to do the integration, right? I guess this would increase the memory efficiency without reducing the accuracy. 'More' grid points are useless as they are spectral truncated.

I see your concern on the summation over grid points. But the offline diagnoses would be more popular in practice. I and @mini-DONG are trying to see how well the offline calculation is close to the online calculation. We can put all these into a doc section once we figure these out clearly. You can point us to an outline or a template (probably the Geostrophy example page), so that we can fill in step by step.

[milankl]

So even if I prescribe the full grid for SpectralGrid at the beginning, the model in the internal will still use its corresponding reduced grid to do the integration, right?

No no. You get what you ask for. SpectralGrid(Grid=SomeGrid) will use that SomeGrid internally. It's just that the default output would be the full grid-equivalent of SomeGrid. Check here

julia> spectral_grid = SpectralGrid(trunc=31, nlev=5)
SpectralGrid:
├ Spectral:   T31 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid:       48-ring OctahedralGaussianGrid{Float32}, 3168 grid points
├ Resolution: 401km (average)
└ Vertical:   5-level SigmaCoordinates

julia> output = OutputWriter(spectral_grid, ShallowWater)
julia> output.output_Grid
FullGaussianGrid

julia> spectral_grid.Grid
OctahedralGaussianGrid

if you try outputting directly on the OctahedralGaussianGrid you'll get an error that relates to the fact that that grid cannot be represented as a matrix:

julia> output = OutputWriter(spectral_grid, ShallowWater, output_Grid=OctahedralGaussianGrid);
ERROR: MethodError: no method matching get_nlon(::Type{OctahedralGaussianGrid}, ::Int64)

Closest candidates are:
  get_nlon(::Type{<:SpeedyWeather.RingGrids.AbstractFullGrid}, ::Integer)

suggesting you to use a full grid instead. (Maybe that error message could be clearer).

But the offline diagnoses would be more popular in practice.

I know because people don't expect that you can do these things interactively with a model! I think that's a massive problem in the climate modelling community that we have such a strong divide between the model and its analysis. The model is too often a stiff blackbox with code that you cannot reuse, and so everyone invents their own workflow based on the data the model outputs. What a missed opportunity! Models should be libraries that allow you to do all these things directly! That's at least what I'm working towards with SpeedyWeather!!

We can put all these into a doc section once we figure these out clearly. You can point us to an outline or a template (probably the Geostrophy example page), so that we can fill in step by step.

I would create a pull request but then you cannot add to it. Can you create a dummy pull request? Then I'll add the outline in the docs and you can fill in the rest. That would be very much appreciated!