Nested Angular Grids (for adaptive quadrature)

[PaulWAyers]

We may need to manually generated nested grids. One can generate icosahedral nested grids by subdividing an icosahedral triangulation. The sequence would be $12, 48=4 \times 12, 192 = 4 \times 48, 768 = 4 \times 192, 3072 = 4 \times 768, \ldots$. Once the sequence is found, one has to solve for the weights that ensure that as many spherical harmonics as possible are integrated correctly. There are $(\ell+1)^2$ spherical harmonics with angular momentum less than or equal to $\ell$, so (at a minimum) the aforementioned grids have degree $3,6,13,27,55,\ldots$. In general, due to symmetry, one should be able to do a bit better than this.

I asked ChatGPT to make a script to generate the nested grid points, and got the following (not checked); the example at the end is for 192 points. The angular degrees would probably be updated based on these grids. We'd have angular grids:
(48, 192, 768) fine
(192, 768, 3072) crazy big
We'd need to figure out which radial grids are a good match for these angular grids.

import numpy as np

def icosahedron_vertices():
    # Golden ratio
    phi = (1 + np.sqrt(5)) / 2

    # Normalize the length to 1
    a = 1 / np.sqrt(1 + phi**2)
    b = phi * a

    # Vertices of an icosahedron
    vertices = [
        (-a,  b,  0), ( a,  b,  0), (-a, -b,  0), ( a, -b,  0),
        ( 0, -a,  b), ( 0,  a,  b), ( 0, -a, -b), ( 0,  a, -b),
        ( b,  0, -a), ( b,  0,  a), (-b,  0, -a), (-b,  0,  a)
    ]
    return np.array(vertices)

def subdivide(vertices, faces):
    # New vertex set, including old vertices
    new_vertices = vertices.tolist()
    vertex_map = {}

    # Helper function to find the midpoint and index
    def midpoint_index(v1, v2):
        # Check if we have already calculated this
        edge = tuple(sorted([v1, v2]))
        if edge in vertex_map:
            return vertex_map[edge]

        # Calculate midpoint and normalize it to the sphere
        midpoint = (vertices[v1] + vertices[v2]) / 2
        midpoint /= np.linalg.norm(midpoint)
        
        # Store the new vertex index
        index = len(new_vertices)
        new_vertices.append(midpoint)
        vertex_map[edge] = index
        return index

    new_faces = []

    for v1, v2, v3 in faces:
        # Get the midpoints
        a = midpoint_index(v1, v2)
        b = midpoint_index(v2, v3)
        c = midpoint_index(v3, v1)

        # Create four new faces
        new_faces.extend([
            [v1, a, c],
            [v2, b, a],
            [v3, c, b],
            [a, b, c]
        ])

    return np.array(new_vertices), np.array(new_faces)

def generate_icosphere(subdivisions):
    # Initial icosahedron vertices and faces
    vertices = icosahedron_vertices()
    faces = [
        [0, 11, 5], [0, 5, 1], [0, 1, 7], [0, 7, 10], [0, 10, 11],
        [1, 5, 9], [5, 11, 4], [11, 10, 2], [10, 7, 6], [7, 1, 8],
        [3, 9, 4], [3, 4, 2], [3, 2, 6], [3, 6, 8], [3, 8, 9],
        [4, 9, 5], [2, 4, 11], [6, 2, 10], [8, 6, 7], [9, 8, 1]
    ]

    # Subdivide the triangles the number of times specified
    for _ in range(subdivisions):
        vertices, faces = subdivide(vertices, faces)

    return vertices

# Generate vertices for a subdivided icosahedron with 2 subdivisions
icosphere_vertices = generate_icosphere(2)
icosphere_vertices

[izxle]

We (Anuar and I) have been working on this problem and have implemented a notebook that follows the approach outlined in this issue. Our implementation generates a radial grid and iterates over the radial points, adding the angular grid using an icosphere. It then evaluates the integral for a given function, increasing the order of the icosphere adaptively until the integral converges (based on a specified threshold). Once convergence is reached at a given radial point, the process moves to the next radial point.

A key limitation of our current approach is that it recalculates everything on each iteration when increasing the order of the icosphere. This is inefficient, as new points could be computed incrementally rather than recalculating the entire integral. However, this incremental update only works when integrating functions like the density. When applying this method to compute the potential, updating only the new points would lead to inconsistencies, as the Laplacian depends on the total number of points. In such cases, recalculating the integral over the entire updated grid is necessary to obtain consistent results.

We've submitted our implementation as a pull request in the nested-angular-grid branch, including the .ipynb notebook in the examples directory. We look forward to any feedback or suggestions for further improvements!

[PaulWAyers]

Thanks. I will discuss with @marco-2023 and we will get back to you. We may request a meeting, depending.

[PaulWAyers]

Update: I talked to @marco-2023 yesterday so sometime late this week or next week you should hear more.

[azhar-ikhtiarudin]

Hi @PaulWAyers @marco-2023, I have run and studied the implementation by @izxle in this pull request, and I appreciate the clear tutorial and step-by-step explanations in the notebook. I have a few questions regarding the angular iteration:

1. Is it possible to iteratively refine the angular grid based on the error difference between the previous and the refined grid, without calculating the true value from a reference grid? That is, if the difference between two consecutive steps is below a certain threshold, we no longer need to refine the grid.

2. Is it possible to locally refine each face of the icosahedral grid independently? Specifically, for each triangular face, we compare the difference between the face before and after increasing the number of vertices. If this difference exceeds a certain threshold, we refine only that specific region, independently for each icosahedral face.

3. Would this approach be feasible or practical? Or could it introduce other issues, such as difficulties in obtaining the weights of each point or problems when integrating across all points?

I would like to clarify this first, but if feasible, may I try implementing it? If this approach is not possible, are there any alternative methods to improve the implementation?

Thank you!

[PaulWAyers]

1. It is definitely feasible (and preferable) to sequentially expand the grid until the angular integral converges. In some cases, a certain integral may need to be used at the beginning to "establish" the grid, after which all successive integrals would use the same grid. This is important mostly in cases where you want to use the same grid for a set of integrals. In such cases, the "first" integral (which selects the angular grids) should be (one of) the trickiest/most difficult in your set.
2. It's a good idea. there are cases you would not want to do this (because one wishes to benefit from symmetry) but it would be a possible strategy. The issue would be the need to assess each face's accuracy and to work with the alternative triangulations that appear one one bisects one or two edges of a triangular face, but not all three.

[izxle]

Thanks for your response. We’ll wait for further feedback and are happy to discuss any details if needed.

In the meantime, we’re working on additional analyses to compare the results of our adaptive method with preset grids, as well as evaluating differences between grids optimized for density vs. those optimized for potential. We expect this to provide more insight into the strengths and limitations of this approach. Looking forward to hearing your thoughts!