Spectral gradient leakage?

[milankl]

When starting with a random vorticity field in spectral space (randn(Complex{NF}) for l,m in 2:25,2:25 (or other lmax,mmax) one should be able to

1. [spectral space] Obtain the stream function by inverting the spectral laplace operator
2. [spectral space] Obtain coslatu, coslatv from the zonal&meridional gradient of the stream function
3. [spectral->grid space] Obtain coslatu, coslatv in grid space with the spectral transform
4. [grid space] unscale coslatu, coslatv to obtain u,v without the coslat scaling
5. [grid->spectral space] Obtain u, v in spectral space via spectral transform
6. [spectral space] Calculate the zonal gradient of v and the meridional gradient of u
7. [spectral space] subtract the zonal gradient of v and the meridional gradient of u to obtain vorticity
8. [spectral->grid->spectral space] unscale vorticity with /coslat in grid space

And hence reobtain the original vorticity field, this works somehow:

However, it's unclear to me how exact this method should be, and whether there's another (simpler) loop one can do to check that all gradients work as they are supposed to. @white-alistair @maximilian-gelbrecht any idea?

[milankl]

In spectral space using T85 this "leakage" looks like

which could mean that something is wrong in the meridional gradient because of the visual leakage across degree's l (y-axis) of the spherical harmonics?

[maximilian-gelbrecht]

I saw similar differences/errors when I was trying to get my spectral transforms and derivatives right. This is just the test from test/spectral_gradients.jl, right?

I can have a look at it tomorrow, but I only implemented the meridional gradient with a pseudo-spectral approach and not the recursion relation. One possibility is also to check the gradients against the ones of the SPHEREPACK library. Not for a proper a unit test, but to be more certain that they are correct in principle.

[milankl]

This is just the test from test/spectral_gradients.jl, right?

I haven't pushed it yet to #60 as I'm still figuring out how to best formulate the test, but it's something like

NF = Float64
prog_vars,diag_vars,model_setup = initialize_speedy(NF,trunc=85)
(;spectral,geometry) = model_setup.geospectral
(;radius_earth) = model_setup.parameters

# use only one leapfrog index from vorticitiy
vor = view(prog_vars.vor,:,:,1,:)    

# some large scale initial conditions (zero mean with all l=0 modes being zero)
lmax,mmax = 50,50
vor[2:lmax,2:mmax,:] = randn(Complex{NF},lmax-1,mmax-1,
                                        model_setup.parameters.nlev)
SpeedyWeather.spectral_truncation!(vor,lmax,mmax)   # set upper triangle to 0

# convert spectral vorticity to spectral stream function and to spectral u,v and transform to u_grid, v_grid
SpeedyWeather.gridded!(diag_vars,prog_vars,model_setup)
(;u_grid, v_grid) = diag_vars.grid_variables    # retrieve u_grid, v_grid from struct
u = zero(vor)
v = zero(vor)

SpeedyWeather.scale_coslat!(v_grid,geometry)
SpeedyWeather.spectral!(u,u_grid,spectral)
SpeedyWeather.spectral!(v,v_grid,spectral)

# zonal gradient of v, meridional gradient of u for vorticity
(;coslat_u,coslat_v) = diag_vars.intermediate_variables
SpeedyWeather.gradient_longitude!(coslat_v,v,radius_earth)
SpeedyWeather.gradient_latitude!(coslat_u,u,radius_earth)

# subtract to obtain vorticity + unscale in grid space
vor2 = coslat_v - coslat_u
vor_grid2 = gridded(vor2[:,:,1])
SpeedyWeather.unscale_coslat!(vor_grid2,geometry)

# back to spectral for plotting
SpeedyWeather.spectral!(view(vor2,:,:,1),vor_grid2,spectral)

# plotting
fig,(ax1,ax2,ax3) = subplots(3,1)
q1 = ax1.imshow(abs.(vor[:,:,1]))
colorbar(q1,ax=ax1)
q2 = ax2.imshow(abs.(vor2[:,:,1]))
colorbar(q2,ax=ax2)
q3 = ax3.imshow(abs.(vor[:,:,1])-abs.(vor2[:,:,1]))
colorbar(q3,ax=ax3)
ax1.set_title("Vorticity, random",loc="left")
ax1.set_title("a",loc="right")
ax2.set_title("Vorticity->stream function->u,v->vorticity")
ax2.set_title("b",loc="right")
ax3.set_title("Error: a-b")

[milankl]

test is now in #main but not run https://github.com/milankl/SpeedyWeather.jl/blob/f21b69eae24f9be0b0667ebe381671d994b9b14f/test/spectral_gradients.jl#L52-L97

[milankl]

thanks to the tests in #75 which pass with https://github.com/milankl/SpeedyWeather.jl/pull/92/commits/a68b670099d9cfaafd1b60bcdce23ff64e5c3cef there should be no spectral gradient leakage anymore.