global spectrum

[miniufo]

Hi, I am trying to analyze the energy at different scales from a set of simulations based on this nice model. So how could I do this on the spherical earth?

I know I could do FFT along the zonal direction, and make sure some scaling can ensure that the energy in the spectral space sum up to the total energy in the physical (gridded) space. But how to do this in the meridional direction? It is not periodic along the meridion and some window function should be applied if FFT is applied, which also modified the total energy to some extent. Maybe the spherical harmonics will help here.

Any advice will be much appreciated.

[milankl]

Check Power spectrum and the power_spectrum function! I'm not 100% sure that the scaling is correct to get equal variances / energy in either space following Parseval/Plancherel theorem but it would be good to check and then we can add a scale correction!

[miniufo]

Great! I will take a look at it and report back.

Just curious, Does the power spectrum here, based on spherical harmonics, only give a scalar spectrum, not zonal and meridional ones (sometimes one may need something like this)? However, those harmonics are composed by zonal FFT and meridional Legendre eigen-functions. So I am trying to figure out the conceptual picture.

[milankl]

You could take abs.(L) to get the amplitude of every spherical harmonic in a lower triangular matrix L that way you would have it split between the zonal harmonics (first column) all the way to the sectoral harmonics (the diagonal). In the end the power spectrum is just the average or sum of that (I'm not sure which one you'd need to match the total energy) but you could just try and figure it out easily with some random field I guess!

[miniufo]

I have try to get the spectra of u, v, and h during a simulation. I used the template of a callback below:

# add callbacks for power spectrum
function spectrum_diagnostics(diagn::DiagnosticVariables, model::AbstractModel)
    u = diagn.grid.u_grid[:, 1].data
    v = diagn.grid.v_grid[:, 1].data
    η = diagn.grid.pres_grid.data
    
    umap = FullClenshawGrid(u, input_as=Matrix) # is this efficient repeating three times???
    vmap = FullClenshawGrid(v, input_as=Matrix)
    ηmap = FullClenshawGrid(η, input_as=Matrix)
    
    uspec = SpeedyTransforms.power_spectrum(transform(umap))
    vspec = SpeedyTransforms.power_spectrum(transform(vmap))
    hspec = SpeedyTransforms.power_spectrum(transform(ηmap))

    return uspec, vspec, hspec
end

# define a SpectrumDiagnostics callback and the fields it needs
Base.@kwdef mutable struct SpectrumDiagnostics <: SpeedyWeather.AbstractCallback
    timestep_counter::Int = 0

    time::Vector{DateTime} = []
    uspec::Vector{Float64} = []  # u spectrum
    vspec::Vector{Float64} = []  # v spectrum
    hspec::Vector{Float64} = []  # η spectrum
end

# define how to initialize a SpectrumDiagnostics callback
function SpeedyWeather.initialize!(
    callback::SpectrumDiagnostics,
    progn::PrognosticVariables,
    diagn::DiagnosticVariables,
    model::AbstractModel,
)
    # replace with vector of correct length
    n = progn.clock.n_timesteps + 1    # +1 for initial conditions
    callback.time = zeros(DateTime, n)

    uspec, vspec, hspec = spectrum_diagnostics(diagn, model)
    
    callback.uspec = zeros(n) # how to get the right size???
    callback.vspec = zeros(n)
    callback.hspec = zeros(n)

    callback.time[1] = progn.clock.time
    callback.uspec[1] = uspec  # array assignment could be wrong???
    callback.vspec[1] = vspec  # 
    callback.hspec[1] = hspec  # 

    callback.timestep_counter = 1  # (re)set counter to 1

    return nothing
end

# define what a GlobalDiagnostics callback does on every time step
function SpeedyWeather.callback!(
    callback::SpectrumDiagnostics,
    progn::PrognosticVariables,
    diagn::DiagnosticVariables,
    model::AbstractModel,
)
    callback.timestep_counter += 1
    i = callback.timestep_counter

    uspec, vspec, hspec = spectrum_diagnostics(diagn, model)

    # store current time and diagnostics for timestep i
    callback.time[i] = progn.clock.time
    callback.uspec[1] = uspec
    callback.vspec[1] = vspec
    callback.hspec[1] = hspec
end

# define how to finish a GlobalDiagnostics callback after simulation finished
function SpeedyWeather.finalize!(
    callback::SpectrumDiagnostics,
    progn::PrognosticVariables,
    diagn::DiagnosticVariables,
    model::AbstractModel,
)
    n_timesteps = callback.timestep_counter

    # create a netCDF file in current path
    ds = NCDataset(joinpath(pwd(), "spectrum_diagnostics.nc"), "c")

    # save diagnostics variables within
    defDim(ds, "time"  , n_timesteps)
    defVar(ds, "time"  , callback.time , ("time",))
    defVar(ds, "uspec" , callback.uspec, ("time",))
    defVar(ds, "vspec" , callback.vspec, ("time",))
    defVar(ds, "hspec" , callback.hspec, ("time",))

    close(ds)

    return nothing
end

spectrum_diagnostics_recorder = SpectrumDiagnostics()
add!(model.callbacks, :diagnostics_recorder => spectrum_diagnostics_recorder)

The power spectrum codes are copied from here. I am not sure about:

1. Are the codes efficient if I need to repeat it three times to get spectra of u, v, and h?
2. How to store the results (I am not sure how to allocate the correct size of the vector to store the results).
3. Do I need to take abs.(L), as suggested, before writing out the results?

After this, I should have gridded fields and spectra, which makes it possible to validate Parseval/Plancherel theorem

[milankl]

I think you have two options

1. compute the spectra on the fly with a callback as you did above (although I would change a few things, more below)
2. just write the netcdf output and compute the spectra from there given that you want to compare it to the variance/energy in the gridded fields anyway?

For 1)
Given that you have both spectral (in the prognostic variables) and grid space (in the diagnostic variables) available on every time step just calculate the spectra from the spectral data, that's super cheap but transforming from grid to spectral and then compute the spectra is just undoing the spectral->grid transform that's happening on every time step anyway.

julia> progn
PrognosticVariables{Float32, Array}
├ vor:   T31, 1-layer, 2-steps LowerTriangularArray{Float32}
├ div:   T31, 1-layer, 2-steps LowerTriangularArray{Float32}
├ temp:  T31, 1-layer, 2-steps LowerTriangularArray{Float32}
├ humid: T31, 1-layer, 2-steps LowerTriangularArray{Float32}
├ pres:  T31, 1-layer, 2-steps LowerTriangularArray{Float32}
├ random_pattern: T31, 1-layer LowerTriangularArray{Float32}
├┐ocean: PrognosticVariablesOcean{Float32}
│├ sea_surface_temperature:  48-ring OctahedralGaussianGrid{Float32}
│└ sea_ice_concentration:    48-ring OctahedralGaussianGrid{Float32}
├┐land:  PrognosticVariablesLand{Float32}
│├ land_surface_temperature: 48-ring OctahedralGaussianGrid{Float32}
│├ snow_depth: 48-ring OctahedralGaussianGrid{Float32}
│├ soil_moisture_layer1:     48-ring OctahedralGaussianGrid{Float32}
│└ soil_moisture_layer2:     48-ring OctahedralGaussianGrid{Float32}
├ particles: 0-element Vector{Particle{Float32}}
├ scale: 1.0
└ clock: 2000-01-01T00:00:00

julia> SpeedyTransforms.power_spectrum(progn.pres[1])
32-element Vector{Float32}:
      1.3960999
 226390.88
...
      0.24390882
      0.21223953

This would be the power spectrum of $$\eta$$, you would need to subtract the power spectrum of orography $$hb$$ because $$h = h_0 + \eta - h_b$$ but that's already in model.orography.geopot_surf (and constant) as the surface geopotential is $$\Phi = gh_b$$. But you are right that u, v isn't available in spectral space so if you want to compute those spectra you'd need to do transform those anyway. Hence, I think it's just easier to work on the netCDF data:

For 2)

run with run!(simulation, output=true) then

ds = NCDataset("run_0001/output.nc")    # or other run id, you can also use `model.output.id` (="0001" here), i.e. "run_$(model.output.id)/output.nc"

# load u data, use .var to ignore possible "missing", use [:, :, 1, :] to remove singleton dimension for vertical
u = ds["u"].var[:, :, 1, :]

# wrap into `FullGaussianGrid`, this doesn't copy any data! = it's efficient.
u_grid = FullGaussianGrid(u, input_as=Matrix);

# transform to spectral space
u_spec = transform(u_grid)

# you can also create and reuse a spectral transform object
S = SpectralTransform(u_grid)
u_spec = transform(u_grid, S)

# now compute the spectra, one for each time step in the netcdf file
spectra = power_spectrum(u_spec)

which returns a lmax x ntimesteps matrix with one spectrum per timestep. This however only works with #618 without that you'd need to manually loop over the timesteps like so SpeedyTransforms.power_spectrum(u_spec[:, 1]) each returning a vector that's a column in the above.

[miniufo]

Didn't realize that the prognostic variables are spectral and diagnostic ones are gridded. I think your options 2) is more flexiable for the purpose, just hoping that the write-out and read-in do not lose any information as the on-the-fly calculation.

I played with my own dataset:

ds = NCDataset("/Data/Speedy/SWMrun/run_T127/output.nc")

u = ds["u"].var[:, :, 1, :]

u_grid = FullGaussianGrid(u, input_as=Matrix) # u[x, y, t] is 512 * 255 * 201

but it showed:

130560×201 Matrix{Float32} cannot be used to create a 256-ring FullGaussianArray{Float32, 2, Matrix{Float32}}

Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] error_message(data::Matrix{Float32}, nlat_half::Int64, rings::Vector{UnitRange{Int64}}, G::Type, T::Type, N::Int64, A::Type)
   @ SpeedyWeather.RingGrids ~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/general.jl:107
 [3] FullGaussianArray
   @ [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/grids/full_gaussian.jl:19](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/grids/full_gaussian.jl#line=18) [inlined]
 [4] AbstractGridArray
   @ [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/general.jl:115](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/general.jl#line=114) [inlined]
 [5] (FullGaussianGrid)(data::Matrix{Float32}, input_as::Type{Vector})
   @ SpeedyWeather.RingGrids ~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/general.jl:124
 [6] AbstractFullGridArray
   @ [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl:29](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl#line=28) [inlined]
 [7] AbstractFullGridArray
   @ [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl:34](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl#line=33) [inlined]
 [8] (FullGaussianGrid)(M::Array{Float32, 3}; input_as::Type)
   @ SpeedyWeather.RingGrids [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl:29](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl#line=28)
 [9] top-level scope
   @ In[27]:5

I changed to:

ds = NCDataset("/Data/Speedy/SWMrun/run_T127/output.nc")

u = ds["u"].var[:, :, 1, :]

u_grid = FullGaussianGrid(u[:,:,1], input_as=Matrix) # u[x, y] is 512 * 255

but it showed:

130560-element Vector{Float32} cannot be used to create a 256-ring FullGaussianArray{Float32, 1, Vector{Float32}}

Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] error_message(data::Vector{Float32}, nlat_half::Int64, rings::Vector{UnitRange{Int64}}, G::Type, T::Type, N::Int64, A::Type)
   @ SpeedyWeather.RingGrids ~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/general.jl:107
 [3] FullGaussianArray
   @ [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/grids/full_gaussian.jl:19](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/grids/full_gaussian.jl#line=18) [inlined]
 [4] AbstractGridArray
   @ [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/general.jl:115](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/general.jl#line=114) [inlined]
 [5] (FullGaussianGrid)(data::Vector{Float32}, input_as::Type{Vector})
   @ SpeedyWeather.RingGrids ~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/general.jl:124
 [6] AbstractFullGridArray
   @ [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl:29](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl#line=28) [inlined]
 [7] AbstractFullGridArray
   @ [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl:34](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl#line=33) [inlined]
 [8] (FullGaussianGrid)(M::Matrix{Float32}; input_as::Type)
   @ SpeedyWeather.RingGrids [~/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl:29](http://192.168.1.105:8888/home/qianyk/.julia/packages/SpeedyWeather/7mNUb/src/RingGrids/full_grids.jl#line=28)
 [9] top-level scope
   @ In[28]:5

Does this relate to the version of my speedy?

[milankl]

It sounds like you're using one type of grid for output (maybe FullClenshawGrid?) but wrap it here into a FullGaussianGrid? Those have to be consistent. If you're worried about losing information but also in general, I recommend using the gaussian grids, e.g. with SpectralGrid(Grid=OctahedralGaussianGrid, ...) (which is the default) you'd output on the FullGaussianGrid by default and then wrap the netCDF variable into a FullGaussianGrid again!

[milankl]

I agree that the error message isn't the most informative. Maybe we should just have a function like SpeedyWeather.grid(::NCVariable) that takes any information in the netCDF file and does the wrapping back into a grid for you!?

[miniufo]

OK, you are right. I see the inconsistency. So the Gaussian grid is recommended. I just don't remember why I used FullChenshawGrid, with dealiasing=3. I need to go through previous posts. Maybe the FullClenshawGrid is information-keeping when output to an even lat-lon grid?

Then how to sumup all the energy on a Gaussian grid? For an equal lat-lon grid, I used sum((u^2+v^2)*dS/2) / sum(dS) where dS is the area of each grid point (this weight can also be $\cos\phi$). For an accurate integration on Gaussian grid, one may need some efforts (Gaussian quadrature?).

Maybe we should just have a function like SpeedyWeather.grid(::NCVariable) that takes any information in the netCDF file and does the wrapping back into a grid for you!?

Great, that's really thinking on the user side! But I am not sure how to do this. Maybe we could store the grid type as a global netcdf attribute.

[milankl]

Great, that's really thinking on the user side! But I am not sure how to do this. Maybe we could store the grid type as a global netcdf attribute.

Yes, that's technically all that's needed an then a convenience function around the functionality of NCDatasets.jl

I just don't remember why I used FullChenshawGrid, with dealiasing=3

you can do that too, the clenshaw grids are exact with dealiasing >=3. But that also makes the grid bigger and therefore increases the netCDF files as well as makes the model slower.

Then how to sumup all the energy on a Gaussian grid?

Best as before, transform to spectral space and use the $l=m=0$ harmonic divided by the norm. $dS$ should also work although the integration is not as exact. Here remember that the latitudes aren't equally spaced so you also need to take into account the latitude spacing like latd[i] - latd[i+1] etc

[miniufo]

you can do that too, the clenshaw grids are exact with dealiasing >=3. But that also makes the grid bigger and therefore increases the netCDF files as well as makes the model slower.

Can we use reduced grid during model integration and then output with a full grid so that we can both make the model faster and keep output as exact as possible?

We also tried to verify the Parseval/Plancherel theorem. We get the total energy from

• global integration over the FullChenshawGrid output $E_{grid}=\iint u^2+v^2 dA$;
• global integration of $u_{spec}$ (from Power spectrum) over wavenumber space $E_{spec1}=\int u_{spec}+v_{spec}dk$; and
• similar to above but $E_{spec2}=\int (u_{spec}+v_{spec})kdk$.

For a simulation of 30 days, we got the three time series as:

So there are a couple of things we are not sure:

• What is the unit of $u_{spec}$ (output from Power spectrum)?
• How to convert the wavenumber $k$ to wavelength $L$? Is it something like $L=2\pi R / k$?
• Do we need to multiply $k$ before integration over wavenumber space? It looks so from the figure. But this may drop the energy at $k=0$.
• $E_{grid}$ is quite similar in shape to $E_{spec2}$. To verify the Parseval/Plancherel theorem, $E_{spec2}$ should be multiplied by $R^2$, but there are still a O(1) constant (~2.5).

Do you have any suggestions here?

[mini-DONG]

Hi，it is very convenient to use the a method SpeedyTransforms.power_spectrum , and we analyzed the output from Shallow Water model of T682:

Regarding this figure above, we want to know how this spectrum compares with observations, that is, how to convert the wavenumber (T682) into the real wavelength (kilometer). I don't know if this is related to the effective resolution, but there are many ways to estimate the effective grid resolution of spectral models.
And, we should use which figure from panel a, b? It seems that panel b (power * wavenumber) fits the power law (-3) better.

[milankl]

How to convert the wavenumber k to wavelength L ? Is it something like L = 2 π R / k ?

Yes, look at the harmonics here for l=1, m=0, so wavenumber k=1, this is a wave that goes from positive (N) to negative (S) to positive (N) so exactly one wavenlength would be around the globe.

What is the unit ?

should be $m^3/s^2$, because $u$ has m/s, power means you square, see here

SpeedyWeather.jl/src/SpeedyTransforms/spectrum.jl

Line 25 in 28dbcac

spectrum[l, k] += 2*abs(spec[l, m, k])^2

it and then it should be divided by the "frequency dimesion" given we're computing the spectrum across space in meters, it's "frequency" is the inverse 1/m. if you do spectrum of u plus spectrum of v you essentially compute u^2 + v^2 per meter so that's a kinetic energy density. See also here. But please check that this is all correctly implemented. I haven't usede this functionality much yet and haven't tested it thoroughly therefore.

Do we need to multiply k before integration over wavenumber space? It looks so from the figure. But this may drop the energy at k = 0 .

Not sure what's needed to close Parseval's theorem, never attempeted that with spherical harmonics. But given above and because you integrate the spectral density (=energy per meter) that makes sense to me.

but there are still a O(1) constant (~2.5).

Note that we use a definition of the spherical harmonics such that the l=m=0 mode is the mean of the gridded field but mulitplied by $2\sqrt{pi}$ maybe that's missing?

[milankl]

Also relevant for you might be to check how the default HyperDiffusion compares to SpectralFilter, see #601. You'd use the latter with

horizontal_diffusion = SpectralFilter(spectral_grid)

and just pass that on to the model constructor. Both HyperDiffusion and SpectralFilter have a bunch of options you could play with. I'm saying this because it looks like your spectrum shoots up at the end, this is technically something the diffusion should take care of but for some reason here it doesn't. Depends also on forcing and drag as they can also have significant spectral fluxes, not sure what you use.