Quadrature grid has too many nodes for given L, M

[f0uriest]

for n in range(3,13):

    basis = ZernikePolynomial(n,n)
    grid1 = QuadratureGrid(n,n,0)
    grid2 = ConcentricGrid(n,n,0)

    print(f"{n} {grid1.num_nodes/basis.num_modes:.2f} {grid2.num_nodes/basis.num_modes:.2f}")

3  2.80 1.00
4  3.00 1.00
5  3.14 1.00
6  3.25 1.00
7  3.33 1.00
8  3.40 1.00
9  3.45 1.00
10 3.50 1.00
11 3.54 1.00
12 3.57 1.00

We expect the concentric grid to have as many nodes as the basis has modes (which is True), and the quadrature grid to have double that. However in practice the quadrature grid seems to have 3x-3.5x as many nodes.

[dpanici]

@dpanici double check this. by computing volume of torus or something, try lower L nodes, recheck old notes

[YigitElma]

I have some formula for total number of modes for Zernike Polynomials here. I should check them again but for now, it seems like this is the expected behavior (according to my formula).

[YigitElma]

And for Quadrature grid, we have (L+1)(2*M+1) nodes with N=0. So, the ratio is around 4.

[dpanici]

@dpanici double check this. by computing volume of torus or something, try lower L nodes, recheck old notes

or just integrate R

[dpanici]

Rory is right, here is finding the average R by integrating over a quad grid for DSHAPE, we hit the floor at L=M=14 but the default 0D grid is at the eq grid resolutions of L_grid=52, M_grid=26,

from desc.examples import get
from desc.grid import QuadratureGrid
import numpy as np
import matplotlib.pyplot as plt
eq = get("DSHAPE")
Rs = []
LMs = np.concatenate([np.arange(2,2*eq.L+9,4), np.array([eq.L_grid+10])])
for LM in LMs:
    grid = QuadratureGrid(L=LM,M=LM,N=0)
    data = eq.compute(["R","sqrt(g)"],grid=grid)
    num = np.sum(data["R"]*data["sqrt(g)"]*grid.weights)
    den = np.sum(data["sqrt(g)"]*grid.weights)
    R_avg = num/den
    Rs.append(R_avg)
    print(R_avg)
# compute at the default grid
grid = QuadratureGrid(eq.L_grid, eq.M_grid, eq.N_grid, eq.NFP)
data = eq.compute(["R","sqrt(g)"],grid=grid)
num = np.sum(data["R"]*data["sqrt(g)"]*grid.weights)
den = np.sum(data["sqrt(g)"]*grid.weights)
R_avg = num/den
print(R_avg)

plt.scatter(LMs,np.abs(np.array(Rs)-Rs[-1]))

plt.scatter(grid.L,np.abs(R_avg-Rs[-1]),marker="x",label=f"L={grid.L} M={grid.M} default 0D grid res")
plt.yscale("log")
plt.xlabel("L=M")
plt.ylabel("R0 - R0 at highest res")
plt.axvline(14,label="L=M=14 (should be enough to integrat the eq.L=26 eq.M=13)")

plt.legend()

including some other quantities, seems for easy quantities yes our defaul is too much, but for force error it is not yet at a floor

[dpanici]

@dpanici do this, make nodes created with L/2 not L