Consistency between NetCDF output grid and internal grid?

[maximilian-gelbrecht]

I might misunderstand the way the netCDF is written, but at least for a FullGaussianGrid I expected this to work:

SpeedyWeather.run_speedy(Float32, trunc=ntrunc, σ_levels_half=[0.,0.4,0.6,1.], n_days=5, nlev=3, output=true)

P, D, M = SpeedyWeather.initialize_speedy(NF, initial_conditions=:restart, restart_id=1) # output dt made to fit this ntrunc to save every time step

vor = ncread("run0001/output.nc","vor")

S = M.spectral_transform
G = M.geometry 
Grid = G.Grid

Grid(vor[:,:,1,1]) # yields a 8192-element FullGaussianGrid 

gridded(P.layers[1].leapfrog[1].vor, S) # yields a 4608-element FullGaussianGrid

So half the grid points. What am I missing there?
This also results in things like spectral(Grid(vor[:,:,1,1], S) not working. It would be nice to an easy way to load the netCDF data and then work with the transforms etc from the library.

[milankl]

What is ntrunc? Because initialize_speedy() hasn't trunc specified so it will use T31 by default. If you want to use the spectral transform from initialize_speedy() then you'd need to pass on the same spectral resolution as keyword argument. However, note you can generally also do

map = Grid(vor[:,:,1,1])
S = SpectralTransform(map)    # create a spectral transform based on the size/resolution of map
spectral(map,S)    # use that spectral transform struct, allocate memory, do transform
spectral(map)     # recompute a spectral transform struct, allocate memory, then do transform

P,D,M = initialize_speedy()
gridded(P.layers[1].leapfrog[1].vor)    # create a spectral transform based on the size of the LowerTriangularMatrix vor
gridded(P.layers[1].leapfrog[1].vor,M.spectral_transform)    # reuse the precomputed spectral transform

I've tried to implement all convenient methods for gridded,spectral,SpectralTransform, but do let me know if something is missing. In general, we also can interpolate/coarsen the resolution in spectral space like so

julia> alms = LowerTriangularMatrix(randn(ComplexF32,65,64))    # T63 
julia> P,D,M = initialize_speedy()    # spectral transform for T31
julia> gridded(alms,M.spectral_transform)
ERROR: BoundsError
Stacktrace:
 [1] gridded!(map::...

julia> almst = spectral_truncation(alms,32,31);    # truncate T63 to T31
julia> gridded(almst,M.spectral_transform)
4608-element FullGaussianGrid{Float64}:
   1.0241511954278444
   1.305803380654694
   1.6143697638450374
   1.9524110446149656
   2.3229442295681544
...

[milankl]

On that note, I believe the error message could be specified, in your situation something like "ERROR: BoundsError, Spectral resolutions do not match, T31 spectral transform with T63 data."

[maximilian-gelbrecht]

ntrunc=42 I forgot that.

So, yeah, I can do

map = Grid(vor[:,:,1,1])
S = SpectralTransform(map)    # create a spectral transform based on the size/resolution of map
spectral(map,S)    # use that spectral transform struct, allocate memory, do transform

and put a spectral_truncation in the middle to get the same result.

But I am still confused why the netcdf output is a different grid, when I do look at the routines that I seemed to me that it just saves the full grid directly?

[milankl]

The netcdf output shouldn't be on a different grid, but just because you say initial_conditions=:restart it doesn't automatically use the resolution from the restart file. Maybe we want to change that, but at the moment it just does an interpolation/truncation onto the new spectral resolution here
https://github.com/milankl/SpeedyWeather.jl/blob/5a91ba8d0e90cd80a0441bb8823f7dbc2b97ad47/src/initial_conditions.jl#L78-L92
I did that because in general I wanted to be able to spin up the model at lower resolution and then ramp up to a higher resolution. I do see that one may expect with initial_conditions=:restart that also the spectral truncation is taken over unless further specified. We should probably do that, but it's easier said than done because Parameters is currently immutable, and we'd need to read in restart.jld2 first, before we set trunc in Parameters.

[milankl]

But I am still confused why the netcdf output is a different grid, when I do look at the routines that I seemed to me that it just saves the full grid directly?

Yes, this is the logic we discussed in #134.

• Full grids should be saved directly as full grids in netcdf
• reduced grids should be
  i) either transformed onto their corresponding full grid (e.g. FullGaussianGrid - OctahedralGaussianGrid),
  ii) or resorted into a matrix (which however, currently is a dead end and only useful for visualisation, see discussion in that issue)

However note that with initial_conditions=:restart data is not read from output.nc, but from restart.jld2 in the same folder. The former stores all gridded data at regular time steps, the latter only stores the prognostic variables in spectral space at the last time step. I've split this up because .nc is what people expect for model output, but jld2 also supports output directly in custom types, meaning literally the whole PrognosticVariables struct is stored in there.

[maximilian-gelbrecht]

The netcdf output shouldn't be on a different grid, but just because you say initial_conditions=:restart it doesn't automatically use the resolution from the restart file. Maybe we want to change that, but at the moment it just does an interpolation/truncation onto the new spectral resolution here

Ah, I see, this is what I did not think about. I just assumed restart continues exactly the same model setup. Thanks for the quick reply

[milankl]

No, we could change that though. We'd need

• to output also the Parameters struct
• if initial_conditions=:restart then load the restart.jld2 file first in initialize_speedy (before Parameters(...) is called
• use keyword arguments from restart file if not explicitly specified in initialize_speedy(kwargs...) such that restarting with different parameters is still possible

[maximilian-gelbrecht]

The solution I went for myself as a temporary fix is to just save the kwargs as a NamedTuple separately. I am fine with that.

[milankl]

But that's hardly reproducible whenever we change the default kwargs for Parameters. As in, while the kwargs may not specify a keyword argument that same argument may be different in the version used to create the restart.jld2 than in the version that uses the restart file. restart.jld2 already contains a field version for that purpose, maybe we can extend that funcationality and throw a warning if the major versions don't align?

[maximilian-gelbrecht]

Yeah, that would make sense

[milankl]

Btw, with #226 merged, you can now output onto any (full) grid at any resolution (output_grid, output_NF are defined in default_parameters.jl, output_nlat_half could be easily added).

[milankl]

Closing this as I believe it's solved, but definitely also outdated as we now output on arbitrary full grids via interpolation (#233, #226). We still have the situation that restarting from file doesn't automatically copy over the spectral resolution, which is more a feature than a bug because you may want to change the resolution at restart.

Feel free to reopen if this is still causing problems.