Output grids

[milankl]

Given that we have now the flexibility to use various grids, we need to think about how to output gridded data to netcdf. Several issues

• The octahedral grids cannot be reshaped/raveled into a matrix, although this is possible for the full grids and the healpix grids
• If the latitudes of the compute grid and the output grid don't match then the Legendre polynomials have to be recomputed. That's either expensive or doubles our memory footprint because two sets of Legendre polynomials have to be kept in memory

So I see several possibilities

1. We always output to FullGaussianGrid (currently), but that's expensive and considerably slows down the simulation
2. We output to the full grid that corresponds to the compute grid. E.g. OctahedralGaussianGrid-FullGaussianGrid, OctahedralClenshawGrid-FullClenshawGrid meaning we would need to define a full version to the HEALPix grid and create a S_output = SpectralTransform(S,output_grid=...) which changes parameters like nlon but points to the already precomputed Legendre polynomials
3. Like 2) but we output the HEALPix grids as a matrix by concatenating their basepixel squares.
4. We don't do some spectral transform onto the output grid and interpolate on the spatial grid instead, but that requires a interpolate!(grid::Grid,lat,lon) function for every Grid<:AbstractGrid, which, however, we may need eventually anyway if we want semi-Lagrangian advection at some point

I don't like 1) what we currently do, hence this issue. I do like 2) and 3) and think that's the simplest to implement. 4) is possible but sounds more like a long-term solution. I reckon it's faster than 2) but not than 3) (which is just memory reordering)

@hottad do you agree?

[hottad]

I agree. If the purpose of output is only visualization or casual post-processing, then the option 2) would be most preferable.
Implementing a full grid corresponding to each reduced grid (like octahedral or HEALPix) is quite easy, and maintenance of such code should not be too tedious.
Note that, for this, we do not need SpectralTransform(S,output_grid=...) (a full spherical harmonics transform) because doing direct Fourier transform on each reduced latitude circle and then doing inverse transform onto a full circle (which is essentially zonal interpolation) should suffice. If we do like this, then 2) may not be necessarily slower than 4).
For more serious applications beyond visualization, like coupling with data assimilation, we should provide capacity to dump the model state on the native grid.

[milankl]

Note to self: Data on an octahedral Clenshaw grid can also (for each of the 4 basepixels) be output like

julia> A
7×7 Matrix{Int64}:
  0   0   0   0   5  11  18
  0   0   0   4  10  17  24
  0   0   3   9  16  23  29
  0   2   8  15  22  28   0
  1   7  14  21  27   0   0
  6  13  20  26   0   0   0
 12  19  25   0   0   0   0

with the poles being 0s or NaNs or something else indicating missing. The octahedral Gaussian grid would need to double the diagonal though as there's no grid line directly on the equator

julia> B
8×7 Matrix{Int64}:
  0   0   0   0   5  11  18
  0   0   0   4  10  17  25
  0   0   3   9  16  24  31
  0   2   8  15  23  30  36
  1   7  14  22  29  35   0
  6  13  21  28  34   0   0
 12  20  27  33   0   0   0
 19  26  32   0   0   0   0

[milankl]

Orography on the octahedral Clenshaw-Curtis grid O512 directly output as matrix (no interpolation), following the idea from above, but rotating every quadrant such that the North Pole sits in the middle

• the matrix is always square (not true for Gaussian latitudes)
• At the poles there's a 4x4 grid point 45˚ rotated gap because the first latitude ring starts with 20 points in longitude
• at 0,90,180,270˚E there's a 90˚ turn for latitude rings, quite obvious from the map...

[hottad]

Wow, this is amazing!

[milankl]

With #140 one can now choose between 2 and 3 from above (at least for those grids where things are implemented).

Vorticity at T255 output directly on the octahedral clenshaw grid O256 (north pole view)

[milankl]

@hottad With #140 merged, I just ran a T1023 simulation using the octahedral Clenshaw-Curtis grid O1024 but output is regridded on the fly onto a full Clenshaw-Curtis grid with 4096x2047 grid points (hence reusing the Legendre polynomials, but replanning the FFTs).

Do you have a suggestions of how we should compare the grids? Like how would we determine whether HEALPix is better/worse than OctahedralClenshaw? Some idealised simulation and take, say, OctahedralGaussian as a reference? Or some power spectrum analysis?

[hottad]

Thank you @milankl for the great job. Those figures are mesmerizing!

I think the question of "ClenshawGrid better than HEALPix (or vice versa)?" cannot be answered purely by simulation results.
Clenshaw-Curtis should give identical results as the Gaussian (as long as the implementation is right), but HEALPix will not give the same results as the Gaussian because the quadrature is inexact, so judging superiority based on the closeness
of the solution to the reference solution from the Gaussian grid is not a fair comparison.

That said, comparison with the tried-and-truth Gaussian simulation is very useful to examine at which resolution HEALPix solution becomes reasonable. Spectral transform on HEALPix should become increasingly accurate as the resolution increases, and checking its convergence property by comparison to Gaussian grid is a good idea.

Besides comparison to Gaussian grid, we should check if there is any undesirable feature in the solutions, like

• noises near the both poles
• Gibbs ringing, particularly around steep orography (Alps, Himalayas, Andes etc)

Examination of convergence is perhaps best done with idealized test cases (Williamson's tests or Galewsky for shallow-water, and Jablonowski and Williamson for 3D primitive equation, for example). To check the presence of undesirable features, we should use more realistic set up with realistic orography and physics (if already available).

Power spectrum analysis is also interesting, but it is very sensitive to factors like choice of hyper diffusion, time step,... so interpretation is rather difficult.

If the simulation results are close enough across the grids, we can then decide our preference based on aspects other than accuracy, like computational cost, model behavior with reduced precision arithmetic, ease of data handling, ease of semi-Lagrangian advection implementation, coupling with other component models (ocean, sea-ice, atmospheric composition etc)...