Looking for functions to calculate vorticity and divergence

[winash12]

I am trying to migrate from https://www2.cisl.ucar.edu/resources/legacy/spherepack Spherepack because the package is no longer being maintained. I am in the atmospheric sciences and I am looking for functions to calculate vorticity , shear and divergence. I also need functions that transform from grid to spectral and from spectral to grid plus functions for gradients. Basically everything that SPHEREPACK provides. Is that possible ?

[MarkWieczorek]

These are not implemented natively in the shtools (i.e., the Fortran 95) code. We just added the option to support multiple backends, and ducc is the first one that was implemented. If ducc doesn't support this, could you let us know if there is another open source library that does this? (perhaps SHTns?) In principle, it is not too hard to write wrapper function so that other codes can be used inside of pyshtools.

[winash12]

@MarkWieczorek I did try SHTns but I have this issue with that library - https://bitbucket.org/nschaeff/shtns/issues/38/valueerror-nlat-must-be-even . It's been a year now and that has not been fixed. The issue is my latitudes can be even or odd. I do not create the data. I download it from a model forecast. It's not possible to remove the even latitude because that phenomena could affect downstream.

[mreineck]

Well, both Fortran shtools and ducc can calculate first derivatives, so it should be relatively straightforward to build the required functionality on top of that, correct?

[winash12]

@mreineck I do not know. Vorticity uses first derivatives yes so does divergence. But I have my data on a grid. I also need to convert it to a spectral data.

[mreineck]

Yes, the idea would be to go from your gridded data to spherical harmonics, and compute the maps of both first derivatives from that. With these, you should easily arrive at maps of divergence and vorticity.

What I'm saying is that computing this information probably will not require highly-optimized low-level functionality (that might be painful to implement), but it should be achievable using fairly basic operations involving SHGrids and SHCoeffs.

[winash12]

@mreineck Excellent. You have given me hope where none existed. In what language would one need to calculate vorticity and divergence ?

[mreineck]

I'm not sure in which language your application is written, but the easiest (and still performant) way would probably be to use Python and work directly on top of pyshtools. All the really compute-intensive work (i.e. the spherical harmonic transforms) will be executed under the hood in Fortran or C++ code.

[winash12]

I am quite flexible. I use xtensor (which is a C++ equivalent of numpy). They also have python bindings. I only need 4 functions to call vorticity and divergence and grid to spectral in my code.

[winash12]

@mreineck I looked at the vorticity code in SPHEREPACK. It is a non trivial piece of code to calculate vorticity on a spectral grid i.e. it is not a simple matter of taking derivatives as in a finite differences approach.

[mreineck]

Oh, I wasn't suggesting to use finite differences. My idea was to convert from grid to spherical harmonic coefficients first (using, say, ExpandDH), and from these then back to grids containing the first derivatives (using MakeGradientDH or similar). This approach is in principle exact (up to numerical noise of course) and does not involve approximations. With the partial derivatives on a grid, you should be able to compute the desired quantity on the same grid. I'm fairly certain that Spherepack will use a similar approach, but I didn't have the time yet to look more deeply into the code.

[winash12]

@mreineck Sounds very reasonable. I believe I will have some time in the first week of September to come up with a prototype (example that uses what you have just mentioned) and then we can take it from there. In that example we can output the value of vorticity from pyspharm and SPHEREPACK(pyspharm is the python wrapper around SPHEREPACK) and compare it with the value I get from SHTOOLS.

[winash12]

@mreineck https://github.com/jswhit/pyspharm/blob/master/Lib/spharm.py#L541 This appears to the be methodology and I think I can handle it. We first need to do a vector harmonic analysis. Is that possible with SHTools ? The u and v wind (zonal and meridional velocity) has to be subjected to a vector harmonic analysis. Then the output of that is sent to this Fortran code - ` do i=1,nt nmstrt = 0 do m=1,ntrunc+1 do n=m,ntrunc+1 nm = nmstrt + n - m + 1 divspec(nm,i) = -(sqrt(float(n)float(n-1))/rsphere)

• scalecmplx(br(m,n,i),bi(m,n,i)) vrtspec(nm,i) = (sqrt(float(n)float(n-1))/rsphere)*
• scale*cmplx(cr(m,n,i),ci(m,n,i)) enddo nmstrt = nmstrt + ntrunc - m + 2 enddo enddo`

Also one needs to check for even or odd latitudes.

[mreineck]

The u and v wind (zonal and meridional velocity) has to be subjected to a vector harmonic analysis

Not necessarily, I think (but I might be wrong). My thought was to apply the simple ExpandDH/MakeGridDH procedure first to u (resulting in du/dtheta and du/dphi) and then to v (resulting in dv/dtheta and dv/dphi), which is all you need for divergence and vorticity.

[winash12]

@mreineck - The partial derivatives in the spectral space only gives you the relative vorticity. What I am looking for is the absolute vorticity(relative vorticity + planetary vorticity) . For that you need to add the value of the Coriolis force(that is a function of latitude or the sin of the latitude).

Is it possible to do a vector harmonic analysis ?

[mreineck]

I don't really know all the SHTOOLS capabilities ... the low-level transforms I have been working on only deal with scalar quantities. (My own background is cosmology, where the term "vector harmonic analysis" isn't used at all, but we often perform analysis on spin-n quantities where spin-0 would be scalars and spin-1 are tangential vectors. The conventions seem to differ quite a lot between the different scientific areas, but there must be a way to translate between all the different representations.)

[MarkWieczorek]

shtools doesn't do vector harmonics. Is it possible to write the vorticity equation in terms of scalar harmonics?

At some point, I would like to implement vector harmonics (mostly for working with magnetic fields), but this is not high on my priorities now...

[winash12]

@MarkWieczorek The vorticity is a pseudo vector. So one can try to get away with scalar harmonics. But the vorticity is a function of a vector field. So I do not think the output would be comparable to the output of Spherepack without vector harmonics.

[mreineck]

OK, I'm fairly confident I can provide ducc0 vector harmonic transforms analogous to those in shtns for the next SHTOOLS release. I'm not entirely sure how this functionality can be presented in SHTOOLS's object-oriented interface, since such a call involves at least two different objects (e.g. coeff_spheroidal and coeff_toroidal). We would have to do something like

coeff_spheroidal.expandDH_vector(coeff_toroidal, ...)

which looks a bit asymmetric.

[winash12]

@mreineck Please let me know. I will be happy to try that release out for sure.