Properly extend shape map into interior region #6477

knelli2 · 2025-02-07T00:35:32Z

This is only allowed when the inner surface is spherical and when we are using the "regular" transition (i.e. not reversed).

Proposed changes

Fixes #6451.

A previous PR (#6359) attempted to extend the shape map into the excision, but did so incorrectly. This PR implements this correctly. In order to do this, the meaning of the transition function had to be changed slightly. Previously, the shape map was

$$ \tilde{r} = r \left(1 - \frac{f(r, \theta, \phi)}{r} \sum_{\ell,m} \lambda_{\ell m}Y_{\ell m}(\theta, \phi)\right) $$

and the transition function was $f(r, \theta, \phi)$. The functional form of the shape map hasn't changed, but now it reads

$$ \tilde{r} = r \left(1 - G(r, \theta, \phi) \sum_{\ell,m} \lambda_{\ell m}Y_{\ell m}(\theta, \phi)\right) $$ $$ G(r, \theta, \phi) = \frac{f(r, \theta, \phi)}{r} $$

where $G(r, \theta, \phi)$ is the transition. The $1/r$ factor was just absorbed into the transition function because the full $G$ is different inside the excision, not just $f$.

Also, the transition function inside the excisions is

$$ G(r, \theta, \phi) = \frac{r^2}{r_{\text{min}}^3} $$

which is different than it is in SpEC. The transition defined in SpEC causes the jacobian to be undefined at the center of the shape map, which means the full frame velocity of all the maps can't be computed at the center of the shape map. This new transition inside the excision allows for this and has the correct behavior.

Upgrade instructions

Code review checklist

The code is documented and the documentation renders correctly. Run
make doc to generate the documentation locally into BUILD_DIR/docs/html.
Then open index.html.
The code follows the stylistic and code quality guidelines listed in the
code review guide.
The PR lists upgrade instructions and is labeled bugfix or
new feature if appropriate.

Further comments

knelli2 · 2025-02-11T18:37:31Z

Ok this is ready for review (maybe @nilsvu ?). Just an FYI for @nikwit that this adds the shape map to the ExcisionSpheres, but if you're only evaluating it at the center of the excision, then this shouldn't behave any differently than before when the shape map wasn't included.

nilsdeppe · 2025-02-19T15:26:55Z

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

+                              interpolation_info);
+  }
+
+  auto transition_func_over_radius =


So this is G or G/r? Because my understanding from the dox change is this is G/r

You are correct, this is G/r. I'll add a small comment here to clarify

nilsdeppe · 2025-02-19T15:48:07Z

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Wedge.hpp

@@ -243,26 +276,21 @@ namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
 * \parblock
 *
 * \note If \p reverse is true, then the value multiplying $\Sigma$ in the
- * numerator is now $-|\vec x_0 - \vec P|$ and in the denomintor $\Sigma$ picks
+ * numerator is now $-|\vec x_0 - \vec P|$ and in the denominator $\Sigma$ picks


demonator? 😅 (JK I'm just teasing!)

nilsdeppe · 2025-02-20T22:56:44Z

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp

+      if (multiple_real_roots) {
+        // Numerical Recipes (3rd edition) pg 228, eqs 5.6.11 - 5.6.12
+        const double theta = acos(R / sqrt(cube(Q)));
+        std::vector<double> roots{


why not use an array? or a static vector? remove_if returns an iterator to the end, so you can use it on an array.

I had assumed remove_if would resize the container, but I see I was incorrect. I think what I actually want is to use std::erase_if so the vector actually contains the positive roots.

nilsdeppe · 2025-02-21T02:42:30Z

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Wedge.cpp

+      if (multiple_real_roots) {
+        // Numerical Recipes (3rd edition) pg 228, eqs 5.6.11 - 5.6.12
+        const double theta = acos(R / sqrt(cube(Q)));
+        std::vector<double> roots{


same question about array

This is only allowed when the inner surface is spherical and when we are using the "regular" transition (i.e. not reversed).

nilsvu · 2025-03-06T00:39:46Z

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

+
+  // G / r
+  auto transition_func_over_radius =
+      transition_func_->operator()(centered_coords, {1});


(optional) don't need the curly braces I think

nilsvu · 2025-03-06T00:42:14Z

...main/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp

   */
  virtual double operator()(
-      const std::array<double, 3>& source_coords) const = 0;
+      const std::array<double, 3>& source_coords,
+      const std::optional<size_t>& one_over_radius_power) const = 0;


(optional) default to std::nullopt

Quick google search says having default parameters in pure virtual functions is not recommended so I'll leave it as is.

nilsvu · 2025-03-06T00:56:02Z

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp

+    const std::optional<size_t>& one_over_radius_power) const {
+  DataVector result = source_coords[0];
+  // Go point by point
+  for (size_t i = 0; i < source_coords[0].size(); i++) {


The performance of this might actually be relevant, so if it's easy enough I'd stick with Blaze for this computation instead of doing the loop yourself.

Well for each point, I'd still need to do the checks I do in the double implementation since this can be called with any set of points and to ensure everything is consistent which would require this to still be a loop.

Oh hide the checks behind #ifdef SPECTRE_DEBUG then their performance doesn't matter and you can do loops in them

Please do these like the gradient_impl function you did below, that's also how these functions were implemented before you split them into separate double and DataVector overloads in this PR. Otherwise the code can do something different in debug vs release, which seems very hard to debug 😅

nilsvu · 2025-03-06T01:01:53Z

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp

+  // a_ being positive is a sentinel for reversed.
+  // If we aren't reversed, check near or within r_min_.
+  if (a_ < 0.0) {
+    if (mag + radial_distortion < (1.0 - eps_) * r_min_) {


Why the + radial_distortion here?

I think you have to be careful that you compare r_min and r_max only to the original radius, not the target radius. E.g. in line 157 there's a comparison to the target radius, which I think could be a problem.

It's undoing the distortion assuming we are at the f=1 surface. Therefore if this condition is true, we are guaranteed to be within the excision and if not, we're outside of it. I ported this check over from SpEC.

Below on 157 should still be the target radius, but I forgot to add the radial distortion since I'm trying to see if we are exactly at the excision (because then the formula is simpler). In the Wedge it is correct

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp

src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Wedge.cpp

knelli2 · 2025-03-06T22:28:48Z

@nilsvu Pushed fixups

knelli2 added this to the BBH larger q & higher spins milestone Feb 7, 2025

knelli2 force-pushed the shape_map_fix branch from 82f1113 to a839da1 Compare February 7, 2025 01:04

knelli2 added the in progress Don't review, used for sharing code and getting feedback label Feb 9, 2025

knelli2 force-pushed the shape_map_fix branch from 03d5644 to 3e38d85 Compare February 11, 2025 18:33

knelli2 removed the in progress Don't review, used for sharing code and getting feedback label Feb 11, 2025

nilsdeppe requested changes Feb 21, 2025

View reviewed changes

knelli2 added 4 commits February 26, 2025 09:46

Properly extend shape map into interior region

a531092

This is only allowed when the inner surface is spherical and when we are using the "regular" transition (i.e. not reversed).

Add shape map to excision spheres

722cfd0

fixup

fd84067

fixup. epsilons

f13de35

knelli2 force-pushed the shape_map_fix branch from 3e38d85 to f13de35 Compare February 26, 2025 17:48

nilsvu requested changes Mar 6, 2025

View reviewed changes

fixup. blaze stuff

c023fc6

knelli2 force-pushed the shape_map_fix branch from 8e1bbfc to c023fc6 Compare March 6, 2025 22:33

fixup. blaze stuff

e3d0e6c

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Properly extend shape map into interior region #6477

Properly extend shape map into interior region #6477

knelli2 commented Feb 7, 2025 •

edited

Loading

knelli2 commented Feb 11, 2025

nilsdeppe Feb 19, 2025

knelli2 Feb 21, 2025

nilsdeppe Feb 19, 2025

nilsdeppe Feb 20, 2025

knelli2 Feb 21, 2025

nilsdeppe Feb 21, 2025

nilsvu Mar 6, 2025

nilsvu Mar 6, 2025

knelli2 Mar 6, 2025

nilsvu Mar 7, 2025

nilsvu Mar 6, 2025

knelli2 Mar 6, 2025

nilsvu Mar 6, 2025

nilsvu Mar 7, 2025 •

edited

Loading

knelli2 Mar 7, 2025

nilsvu Mar 6, 2025

knelli2 Mar 6, 2025

knelli2 commented Mar 6, 2025

Properly extend shape map into interior region #6477

Are you sure you want to change the base?

Properly extend shape map into interior region #6477

Conversation

knelli2 commented Feb 7, 2025 • edited Loading

Proposed changes

Upgrade instructions

Code review checklist

Further comments

knelli2 commented Feb 11, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nilsvu Mar 7, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

knelli2 commented Mar 6, 2025

knelli2 commented Feb 7, 2025 •

edited

Loading

nilsvu Mar 7, 2025 •

edited

Loading