Neutron star surface finder algorithm

[nilsvu]

Find the surface of a neutron star in numeric volume data. This is needed to adapt the computational grid to the stellar surface so the surface discontinuity doesn't spoil exponential convergence of our initial data solver. It's also useful to impose the regularity condition on the elliptic equation for the velocity potential, which we want to solve only within the star. Basic idea:

• Choose the desired angular resolution for the surface in terms of spherical harmonic l_max and m_max. Then the goal is to find the radius of the surface at each angular collocation point to construct a Strahlkorper representing the stellar surface (look at NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp).
• Implement a root finding algorithm that determines the radius of zero density or unit enthalpy along every radial "ray" through the angular collocation points.
• Bonus points for finding the isosurface of any value of the enthalpy / density / etc, so we can try to adapt the domain also to discontinuities in equations of state with phase transitions.

Literature: Read up on "surface-fitting coordinates" in the context of binary neutron star initial data. Also talk to me.

[PunJustice]

Working on this! :)

[isaaclegred]

Just a note on this, @yoonso0-0 and I have recently found that after some evolution the surface can get a bit "fuzzy" (like the density goes is small but not atmosphere near the surface of the star) so it will be really helpful to have an algorithm which can find along any enthalpy/density contour, +bonus points for bonus points

[nilsvu]

Here's the plan we came up with today (involved were @PunJustice @ffoucart @nilsdeppe):

• We specialize to blocks where one logical coordinate axis is aligned with the radial direction. This is needed to conform the domain boundary to the star surface anyway. This assumption should be asserted in the surface finder code.
• We find the surface in the specific enthalpy, though the algorithm should work for any variable and target value.
• In each element we do a 2d interpolation along the angular directions to those angular collocation points of the Strahlkorper that lie in that element. This gives us the data on the 1d radial LGL grid points of the element along each ray. Then we either use the 1d Legendre basis directly to find the root along each ray (see implementation in Numpy: https://github.com/numpy/numpy/blob/v1.24.0/numpy/polynomial/legendre.py), or run TOMS748 to find the root, Lagrange-interpolating in each TOMS748 step.
• Either we find the root simultaneously in all elements of the spherical shell and discard those elements where no root is found, or we cascade outwards and inwards from each element to its radial neighbors when no root is found. Either way we send the roots to a 'SurfaceFinder' singleton which builds the Strahlkorper and can adjust a FunctionOfTime for a shape map in the domain.

[nilsvu]

Suggested plan for the parallel algorithm / communication pattern:

• The surface finder is an array parallel component. Surface finders get instantiated by a list in the input file (one per neutron star in the simulation). This gives each surface finder an initial guess Strahlkorper from the input file. Each surface finder controls a FunctionOfTime in the global cache, which in turn controls a shape map in the domain.
• All elements in a specified set of blocks try to find the surface, and those that find it send the radii to the surface finder. Note that the surface will be close to a block boundary and probably oscillates around it, since we're deforming the domain to coincide with the surface. So we could optimize the algorithm to take this into account.
• Elements send surface radii to the surface finders, which update their Strahlkorper and then the coord maps.

[nilsdeppe]

I don't understand why we need a surface finder component. Can't the elements just send the surface they find directly to the non-const global cache?

[nilsvu]

We need a place to assemble the radii coming from different elements into one Strahlkorper. I think we need a component with a DataBox for that, to keep track of filled/missing Strahlkorper collocation points.

[nilsdeppe]

Ah okay. Could we just have one component rather than an array? Since all elements need to do a reduction over the entire domain, we need to have a vector of N CoM calculations anyway and each element decidse which one it contributes to. So naively I think a singleton should be the right answer that then just sends to the different control systems.