jmert / associatedlegendrepolynomials.jl Goto Github PK

Hi we were looking at some of the packages functionality for the CLIMA project, but having the package not registered is kind of a blocker for more than experimental use. It would be great if it were added to the Julia Registry.

The normalization interface functions should be renamed from Plm_* to legendre_*.

I'm finding that if I'm at high multipoles and high elevation, at some point the legndre polynomials seem to maybe underflow, so it becomes all zeros?

using Legendre, PyPlot
θ = deg2rad(160)
matshow(log.(abs.(λlm(0:3000, 0:3000, cos(θ)))))

Am I doing anything wrong here? If not, is there any way to get this to work out there?

It can be annoying to always remember to allocate output arrays with the requisite lmax+1 due to the fact that quantities are specified on the closed interval [0, lmax]. In the context of #6, it also then becomes convenient and/or necessary for the user to hand a properly constructed offset array as well, and that may add to the annoyingness.

Consider whether a helper function should be provided to do vector and matrix-output allocations. (And maybe the hardest problem, what should such a function be named?)

Right now the "broadcasting" syntax (e.g. Plm.(0:lmax, 0:lmax, x)) allocates an output array of the necessary size, but the abuse of broadcasting like that is probably something that should be removed before a v1.0 is tagged.

Right now the use of fma() restricts the ability to calculate upon complex arguments, even though the functions are mathematically defined across the complex plane.

Figure out how to change the code and still maintain accuracy in the typical real-argument case.

Improvements/changes should be validated via a set of benchmarks before being committed. To do that, though, first a suite of benchmarks must be created and added.

• First level goal is to just provide a basic suite which can be run manually.
• Second level would be to add some performance assertions that compare pairs of benchmarks within the suite (such as the scalar calculation should be faster than a length-1 vector calculation).
• Third level is to provide a framework for running and summarizing the benchmarks across multiple architectures (i.e. compare personal laptop and desktop which have Intel vs AMD processors from different generations).

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Running the recurrence relations require some state be held across iterations, and in general this requires allocating some memory to store the state.

There should be a way of pre-allocating the state variables and allowing that temporary storage to be reused across invocations.

Hey Justin, just had a chat with @kburns on spin-weighted spherical harmonics, which

1. remove the cos scaling when calculating gradients (convenience and simplicity)
2. simplify meridional gradients where otherwise the gradient for $\ell$ depends on $\ell+1,\ell-1$

Downside is obviously that one needs to get their head around the theory behind it. The recursion relations for the polynomials change as outlined in this paper https://doi.org/10.1016/j.jcpx.2019.100013 (equation 26-28 particularly). But apart from that, not much actually changes when implementing the transform, which is neat.

Keaton and I thought about starting a project whereby SpeedyWeather could be using spin-weighted harmonics but as a starting point we would need the calculation of the polynomials to support spin-weighted ones. Would you be keen to get those implemented here, either by contributing directly or by providing guidance how this could be done with your package without interfering with its current features? Keaton has implemented these before for Dedalus and understands the maths. I'm happy to help with the implementation but I just learned about the theory ;)

Since the Associated Legendre Polynomials are naturally 0-indexed, we should take advantage of the OffsetArrays package and work with 0-indexed arrays.

Using @fastmath(sqrt) as a way to get a square root calculation without the domain (x < 0) check is a bit of an abuse of the macro. Instead, there should be an unchecked_sqrt function that skips the domain check without also annotating code with fast-math flags since that can introduce unwanted semantics.

A constraint is that we don't want to break the genericness of the implementation, so a reasonable fallback path needs to be indentified. For instance, right now the method called is actually Base.FastMath.sqrt_fast which may already have useful overloads from other number-type packages that we don't necessarily want to lose. Therefore some research across the package ecosystem is probably warranted.

More context in a Julia Discourse post and a related LLVM bug report