Handle Spherical SL solutions for large harmonic orders

[The9Cat]

What

The inner and outer solutions to the Laplace equation scale as $r^l$ and $r^{-1-1}$ respectively. These very strong radial dependencies break the convergence of the SL solver in the tails. The solution is to restrict the evaluation grid of the SL solver for large r, replacing the tabulated values with zero.

Implementation

• The updated code uses the full radial range for $l<L_{switch}$ and uses $r_{min} = r_{map} 10^{-\alpha/l}$ and $r_{max} = r_{map} 10^{\alpha/l}$ for $l\ge L_{switch}$.
• The constants were selected to produce good grids up to $l_{max}=512$. The resulting values are $L_{switch}=32$ and $\alpha=100$.

Tests

• I used slcheck to verify that that SLGridSph generates orthogonal bases in both scalar and MPI modes
• I spot checked high l-order eigenfunctions to make sure that the choice of $\alpha$ produced inner and outer boundaries many orders of magnitudes into the tails.
• This strategy will rely on the user picking an appropriate value for $r_{map}$, since the heuristic choice for the boundary points depends on this value. The exact value is not critical, but $r_{map}$ should be in the heart of the nodes for the basis.

Notes

• I did not change the cache version value since all previous caches will remain valid.
• The heuristic change only affects the SL solver calls in compute_table so the risk of breaking existing behavior is low.
• Also, I recommend increasing ncylr from the default value of 2000 to something proportionately higher. E.g. at least 10000 or more for $l_{max}=256$.
• More tests are welcome, of course, but this seems like a fairly safe set of changes.

[Copilot]

Pull Request Overview

This pull request implements a solution to handle large harmonic orders in spherical Sturm-Liouville (SL) solutions by restricting the evaluation grid for large harmonic orders where the radial dependencies ($r^l$ and $r^{-1-l}$) break solver convergence. The implementation adds heuristic boundary constraints for $l \geq L_{switch}=32$ using $r_{min/max} = r_{map} \cdot 10^{\pm\alpha/l}$ with $\alpha=100$.

• Introduces adaptive grid boundaries for SL solver based on harmonic order to prevent convergence issues
• Refactors variable names for clarity (LMAX→LMAXFID, NMAX→NMAXFID, NORDER→NMAX)
• Adds improved spherical harmonic evaluation with prenormalized recursion for numerical stability
• Enhances debug output and MPI-aware console messaging

Reviewed Changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 6 comments.

Show a summary per file

File	Description
utils/ICs/check_coefs2.cc	Major refactor of parameter names, improved MPI output handling, and debug features
utils/ICs/DiskEval.cc	Added MPI-aware console output and improved thread handling
utils/ICs/DiskEval.H	Implemented prenormalized spherical harmonic evaluation for numerical stability
include/SLGridMP2.H	Added static constants for grid boundary heuristics
exputil/SLGridMP2.cc	Core implementation of adaptive grid boundaries and enhanced SL solver diagnostics

[The9Cat]

I appended a bunch of comment changes to this commit for the sole purpose of checking how they are rendered in the Doxygen+Breathe+Sphinx generation of the manual. I probably should have made a separate PR for this, sorry about that.

Anyway, seems good to me. The lists of YAM configuration parameters is not particular front and center, but they are easy enough to find. A bunch of classes already have documented YAML keys, btw.

[The9Cat]

I think we can go ahead and merge this one. Or perhaps there are some additional tests that we need?