Nodes1D: Nodal points for Quadrature and Interpolation

[klappi-s]

Nodes1D

Adds a header-only ippl::nodes1d module for classical 1D quadrature (or interpolation) nodes and weights.

Functionalities

Includes: Gauss–Jacobi, Gauss–Legendre, Gauss–Lobatto (GLL), and Gauss–Chebyshev (closed form).
Also provides Affine maps for nodes/weights from [-1,1] onto other intervals.

Configurable backends for Root Finding:

• Golub–Welsch (default)
• dense Golub–Welsch
• Newton + Dekker-Brent fallback (https://en.wikipedia.org/wiki/Brent%27s_method)

Configuable Initial Guesses for Gauss Jacobi Newton Root finding:

• Asymptotic (default): Initial Newton guess for the k-th Jacobi root using a Tricomi-type large-n approximation.
• Chebyshev: Initial Newton guess for the k-th Jacobi root using a Chebyshev node. Good near alpha = beta = -1/2.
• StroudSecrest (LehrFEM like legacy): Initial Newton guess for the k-th Jacobi root using a LehrFEM-style heuristic using previously accepted ascending nodes.

Fallback Ladder for Newton Root Finders

computeGaussJacobi

If not all n nodes/roots are found by a newton method, it tries to find additional nodes, by rerunning with different initial guesses or the brent method.
LehrFEM

primary Guess: InitialGuessType(Asymptotic | Chebyshev | StroudSecrest)
→ Asymptotic (if not already primary)
→ Chebyshev (if not already primary)
→ Brent (scans samples of P_n on (-1,1) for sign changes; Brent on each bracket)
→ throw

computeGaussLobatto
No user InitialGuessType.

Endpoints $\pm 1$ ($j = 0$ and $j=n-1$) are fixed. Each interior node (roots of $P_{n-1}'$) is found from left to right ($j = 1,\ldots,n-2$) one after the other. If the j'th node is not found by the first method, the next method in the ladder is attempted and so on.

Asymptotic Jacobi(1,1)
→ Cosine (Chebyshev–Lobatto: -cos((j+1) π / (n-1)))
→ Brent (scans samples of P_{n-1}' on (lo, hi) for sign-change, then widen toward +1 if needed)
→ throw

Here lo is the previous accepted node $(j-1)$ and hi is the midpoint of the raw asymptotic guesses for interiors $j$ and $j+1$. GLL therefore never retries a later guess family for all interiors at once; it finishes one interior before starting the next.

Consumers

• Will be needed in upcoming FEM updates.
• Also replaces NUFFTUtilities.h; FFT/NUFFT/Correction.h now calls computeGaussLegendre from Nodes1D.

Tests:

unit_tests/Nodes1D/

Nodes1D.cpp:
Large test suite set which checks all root finding backends (except LehrFEM...) over the subset of data that is provided in
Nodes1DOracleData.h

Nodes1DOracleData.h:
File providing test data for nodes and weight finding for node samples (2, 3, 7, 10, 17, 64) + test data for the affine map functionality for the 7 Nodes test case.

create_Nodes1DOracleData.py:
Provides Python script which creates the file Nodes1DOracleData.h:. If needed we can produce a large set of test data.
by changing

NS = (2, 3, 7, 10, 17, 64)
AFFINE_NS = (7,)

and then Nodes1D.cpp test suite should still work without any adaptations.

[aaadelmann]

In general is a function is KOKKOS_INLINE_FUNCTION labled or not seams to me a bit random.

[klappi-s]

@aaadelmann

In general is a function is KOKKOS_INLINE_FUNCTION labled or not seams to me a bit random.

The Code Snippets which this submodule is based on are from the NUFFTUtilities and the Quadrature Code from FEM. Both of them do Node computation exclusively on Host side. Since Nodes computation is usually only run once during initialisation, with relatively small n (for FEM n<10 and for NUFFT n=100). So I guess for simplicity, there was reduced focus on performance.

All of computeGauss*** methods are therefor like before written as pure host functions. The few functions writtten with KOKKOS are small functions that "should" be inlined, but instead of inline KOKKOS_INLINE_FUNCTION was used.

So if you wan't I can either switch it all to inline, so it visibly is clear that this is all host side code. Or go over everything and optimise further with complete kokkos integration.

[klappi-s]

I will now additional wire up FEM/Quadrature to use, the Nodes1D module.