Change weight in ZonalJet()

[miniufo]

This is related #518. Changing the weight have reduced the unbalance. But we may still need better algorithm for cumsum along the meridional direction.

[milankl]

Awesome thanks Yu-Kun (is that your preferred first name?)!! Just looking at this I'm actually wondering whether we could just obtain $\eta$ from $u$ in the same way as we do in the documentation: https://speedyweather.github.io/SpeedyWeather.jl/dev/gradients/#Geostrophy

So use $u, v=0$ calculate vorticity $\zeta = \nabla \times (u,v)$, then $\nabla^{-2}(\tfrac{f\zeta}{g}) = \eta$. It's more costly because of the transforms, but maybe it'll already give more balanced initial conditions?

[miniufo]

Awesome thanks Yu-Kun (is that your preferred first name?)

Yes, Yu-Kun is my first name.

I'm actually wondering whether we could just obtain $\eta$ from $u$ in the same way as we do in the documentation

Yes, that's a workaround. But this $\eta$ is the geostrophically balanced one, not the gradient-wind-balanced one. You need to take into account $u^2$ term to be more balanced. This is not a problem in the ocean where $u$ is about 1 m/s. In the atmospehre, $u$ can be 100 m/s so that $u^2$ term is also important.

[miniufo]

Checks are not successful, the same as previously. Not sure how to deal with it.

[miniufo]

So your workaround could be like this (I did some math). The balanced momentum equation is:

\frac{\partial \mathbf u}{\partial t} = 0 = - \zeta_a \hat{\mathbf u} - \nabla\left(\mathbf u^2/2 + g\eta\right)

where $\hat{\mathbf u}=\left(-v, u\right)$. Taking $\nabla\cdot$ on both sides gives the steady divergence equation:

0 = - \nabla\cdot \left(\zeta_a \hat{\mathbf u}\right) - \nabla^2\left(\mathbf u^2/2 + g\eta\right)

After cleaning up, we have:

\nabla^2\left(g\eta\right)= - \nabla\cdot \left(\zeta_a \hat{\mathbf u}\right)  -\nabla^2\left(\mathbf u^2/2\right) 

which is ready for the inversion to get $\eta$.

[milankl]

Thanks, that looks great! I just implemented that, hope I got it all right, please check if you don't mind! In the end it's a few more transforms than I was hoping for (5 vs 2) but it looks stable, e.g. without perturbation T63 the jet remains stable for 100 days (FullGaussianGrid, reduced grids don't seem to be as stable presumably due to transform errors that are different from ring to ring)

with default perturbation after 6 days it looks like Galewsky's jet

$\eta$ initial conditions are now (blue) ~158m higher south of the jet and ~928.2m lower north of it (zero mean due to the inverse Laplacian), previously we had 0-1525m (yellow)

[miniufo]

That's great! I will give a try on this, because it eliminates most of the gravity waves bouncing over the globe.

I suddenly notice that, if I want to put two jets, one in northern hemispehre while the other in the southern one, how can I combine the jets into a single initial condition through using ZonalJet() (possibly different jet parameters)? Maybe a ZonalJet!() is needed?

[milankl]

That's great! I will give a try on this, because it eliminates most of the gravity waves bouncing over the globe.

Awesome, thanks! You can see that prior to this pull request the jet was balanced by increasing $\eta$ south of the jet, the new version would change $\eta$ so that the mean is zero, which I think I prefer because that way the initial conditions don't increase mass/volume and just inject kinetic and potential energy. But maybe you have another opinion on this?

I suddenly notice that, if I want to put two jets, one in northern hemispehre while the other in the southern one, how can I combine the jets into a single initial condition through using ZonalJet() (possibly different jet parameters)? Maybe a ZonalJet!() is needed?

I see two possibilites

• We make every parameter in ZonalJet() a vector, so there's one jet per element in the input vectors, the good thing is that this only affects the calculation of u. Once we have the field of u ($v = 0$ always I would still require because it's a zonal jet) the other calculations are actually all the same and don't have to be repeated.

• An AddZonalJet! (I suggest this name to make clear you're not just writing the initial conditions to be one zonal jet) which recalculates a new u, takes the curl adds that vorticity to the existing vorticity (they're additive). But as I see this because of the $u^2$ term, the $\eta$ is not additive, which means we need to overwrite $\eta$. That's not a problem for the jet, but the perturbation of it, which is added in $\eta$. So I don't see an easy way to retain the information of the $\eta$ perturbation from jet 1 if jet 2 gets added.

[miniufo]

You can see that prior to this pull request the jet was balanced by increasing south of the jet, the new version would change so that the mean is zero, which I think I prefer because that way the initial conditions don't increase mass/volume and just inject kinetic and potential energy. But maybe you have another opinion on this?

That perfect! I did not notice this. The original cumsum always starts from 0, so it may increase or decrease the total mass afterward, depending on the sign of u. The new way is basically the inversion of Poisson equation. The solution is up to a constant. It looks like the constant is determined internally by spherical harmonics so that the global integration of the solution is zero. I have not checked this. But it looks right.

I see two possibilites

I would then prefer the first one. Yes, not all the fields are addictive. So adding u is more intuitive for ZonalJet(), and inverting everything else is then straightforwrd. So I guess the API would be

# I may use python list-style here
initial_condition = ZonalJet(latitude=[-50. 60], umax=[60, 40])

which gives a jet of 60 m/s at 50S and the other jet of 40 m/s at 60N.

[milankl]

The solution is up to a constant.

Yes, which is zero for all Laplace inversions, we have set this here

https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/00a2e2791cb0d20063fdcb391275f2cc984d2588/src/SpeedyTransforms/spectral_transform.jl#L193-L196

So if you wanted another integration constant (e.g. 1) you could always do spectral_transform.eigenvalues.eigenvalues⁻¹[1] = 1 although it's probably not advised to do that with a spectral transform that's also used in the model integration.

So I guess the API would be

Yes, (the syntax for python lists yields in Julia a vector, that's what I had anticipated), the ZonalJet would have fields that are ::Vector and we just need a generator function that also accepts ::Number (i.e. scalars) so that ZonalJet(latitude=10) is still possible which then gets translated to [10.0] meaning there's one jet and there's a loop inside the function that simply loops over all jets in ZonalJet. Alternatively we can split this: ZonalJet as before ZonalJets is for several ones, but that's probably more complicated than clear or helpful ...

[milankl]

I'm actually wondering whether something like this should live in its own repository similar to how we did it with https://github.com/SpeedyWeather/StochasticStir.jl this repository would then have SpeedyWeather as dependency and you could define ZonalJets in there. This would also nicely showcase how you can define your own initial conditions externally and pass them on to model = ShallowWaterModel(;spectral_grid, initial_conditions=ZonalJets())

[milankl]

Documentation deploy is hanging because

┌ Info: Deployment criteria for deploying preview build from GitHub Actions:
│ - ✔ ENV["GITHUB_REPOSITORY"]="SpeedyWeather/SpeedyWeather.jl" occurs in repo="github.com/SpeedyWeather/SpeedyWeather.jl.git"
│ - ✔ ENV["GITHUB_REF"] corresponds to a PR number
│ - ✘ PR originates from the same repository
│ - ✔ `push_preview` keyword argument to depl

so it's just because you're pushing from a fork, I'll merge to have the corrected zonal jet and I feel how to do the multiple jets is another question not in this PR

[milankl]

@miniufo I wanted to showcase how specific model components/setups (in this case initial conditions) can easily live externally and can be developed in their own respective packages having SpeedyWeather as a dependency, so I just went ahead and put the idea of several/multiple/many zonal jets following Galewsky into its own repository

https://github.com/SpeedyWeather/ManyZonalJets.jl

You find some videos in there. I invited you to be part of our SpeedyWeather organization, so you can also go ahead makes changes to that repository (let me know if anything doesn't work with access rights), or you can create your own model components in separate packages. Having them under the SpeedyWeather umbrella just helps for visibility and other developers here (incl me) having full access in case we update things that require updates in those separate packages too -- but otherwise they can be independent and don't add too much code to SpeedyWeather.jl itself.

[miniufo]

That's great! I just joined this group. The video looks great and this is exactly what I need. I shall take look at the details.