unbalanced initial condition for Galewsky Jet

[miniufo]

Following previous discussion, the initial mass and flow fields for Galewsky Jet are not in balance to each other. The ZonalJet() first prescibes a u-field, and then do the meridional integration $h(\phi) = \int_{90}^{\phi} f(u) d\phi'$ to get the mass field h.

I read some docs on Gauss-quadrature method, which approximates the integration by:

$$ h(\phi) = \int_{90}^{\phi} f(u) d\phi' \approx \sum_n f(u_i) w_i $$

The corresponding codes are:

# integration for layer thickness h / interface height η
w = weights[j]
η_sum += 2w*(radius*u_θ/gravity * (f + tan(θ)/radius*u_θ))

It is a bit strange to me that the weight is 2w rather than w. But after some tests (in T170 resolution), both 2w and w won't give the balanced mass field. I tried the naive way of integration and replace 2w by dlat as:

    ...
    latstr = π/2

    for (j, ring) in enumerate(eachring(u_grid, η_grid))
        θ = π/2 - colats[j]             # latitude in radians
        coslat⁻¹j = 1/cos(θ)
        f = 2rotation*sin(θ)

        # velocity per latitude
        if θ₀ < θ < θ₁
            u_θ = umax/eₙ*exp(1/(θ-θ₀)/(θ-θ₁))  # u as in Galewsky, 2004
        else
            u_θ = 0
        end

        # integration for layer thickness h / interface height η
        w = weights[j]
        dlat = latstr - θ
        latstr = θ
        η_sum += dlat*(radius*u_θ/gravity * (f + tan(θ)/radius*u_θ))
       ...
    end

then I get roughly the balanced mass field h. So I am pretty sure there is something wrong with the weight in the integration but cannot go further to check why it is wrong.

The Gauss-quadrature method seems to be best for definite integral from pole to pole (a reduction in meridional direction). Not so sure if it still make sense for a cumulative integration.

[milankl]

Yeah I probably didn't think too hard about how to do the meridional integration with Gaussian quadrature, will look into it.

[milankl]

I didn't find the time yet to look into this. Ideally we should do the integration it with Gaussian quadrature, but if it's currently not correctly implemented (my bad!) then that's obviously not very helpful. I'd be happy for you to create a pull request with any integration that yields a more balanced jet. We can just add a # TODO explaining what's still needed there, but any improvement would be appreciated. A few lines above you can see that the initial conditions are computed on a FullGaussianGrid with some given resolution, you can technically also change that if it leads to better balanced initial conditions, I'd just be reluctant if you wanted to dial up the resolution too much as it's a default that needs to be run frequently for continuous integration etc. But maybe a FullClenshawGrid with your way is already better? Any help would be awesome!!

For context, this is how we use the quadrature weights in the transforms

https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/7f999536df4d3680f9662c3866c3f4bdeddd41a6/src/SpeedyTransforms/spectral_transform.jl#L141-L146

the get_quadrature_weights(::Type{<:AbstractGrid}, nlat_half::Integer) function returns the quadrature weights on one hemisphere, Eq. included, over both hemispheres their sum is 2

julia> using SpeedyWeather

julia> RingGrids.get_quadrature_weights(FullGaussianGrid, 24)
24-element Vector{Float64}:
 0.0031533460523056585    # first ring j=1
 0.007327553901276272
 ⋮
 0.06446616443595003
 0.06473769681268386      # just north of Equator, j = nlat_half

julia> sum(RingGrids.get_quadrature_weights(FullGaussianGrid, 24))
0.9999999999999996

julia> sum(RingGrids.get_quadrature_weights(FullClenshawGrid, 24))
1.0322910836705603   # this includes the Equator that wouldn't be doubled counted for the southern hemisphere

[milankl]

The Gauss-quadrature method seems to be best for definite integral from pole to pole (a reduction in meridional direction). Not so sure if it still make sense for a cumulative integration.

Yes, maybe that's an issue... from Wikipedia

So maybe that's what's missing here when doing cumsum

[miniufo]

I can send a PR, in which the weight will be replaced by $\Delta \phi$. This is equivalent to the trapezoidal integration along the meridional direction. But the result looks better.

So far, I cannot figure out whether the cumsum using Gaussian quadrature is appropriate. My guess is that the "best" weights cannel some error when perform the integration from pole to pole. The intermediate results from cumsum will not have this kind of cannellation.