Make shape map exactly like the one in SpEC #5636

knelli2 · 2023-11-21T01:15:57Z

Proposed changes

This is a follow up to #5627.

In SpEC, the shape map is (for centered coordinates)

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

However, in SpECTRE, the shape map currently is

$$ x^i = \xi^i \left(1 - h(r, \theta, \phi) \sum_{\ell,m} \lambda_{\ell m}Y_{\ell m}(\theta, \phi)\right) $$

In SpEC, we have that $0 \leq f(r, \theta, \phi) \leq 1$, while in SpECTRE we have $0 \leq h(r, \theta, \phi) \leq 1$. This slight difference of a factor of $r$ is a fundamentally different coordinate mapping.

This PR edits the shape map to be identical to the one in SpEC and consequently edits both of the transition functions to be from 1 to 0.

The wedge transition map also has additional edits that fixes some issues it had before.

This PR also fixes #5635.

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 · 2023-11-21T01:36:11Z

src/Domain/CoordinateMaps/TimeDependent/Shape.hpp

+ * \link
+ * domain::CoordinateMaps::ShapeMapTransitionFunctions::ShapeMapTransitionFunction::original_radius_over_radius
+ * ShapeMapTransitionFunction::original_radius_over_radius \endlink.


Wasn't sure what to do about this. Even without the \link it would be over 80 characters.

This looks like a good case of a long line. There's a whitelist in tools/FileTestDefs.sh you can add to. You can either add this file or whitelist \link in general.

Whitelisting \link in general seems like the best approach and is helpful in the future as well.

wthrowe · 2023-11-21T02:24:04Z

src/Domain/CoordinateMaps/TimeDependent/Shape.hpp

+ * \link
+ * domain::CoordinateMaps::ShapeMapTransitionFunctions::ShapeMapTransitionFunction::original_radius_over_radius
+ * ShapeMapTransitionFunction::original_radius_over_radius \endlink.


This looks like a good case of a long line. There's a whitelist in tools/FileTestDefs.sh you can add to. You can either add this file or whitelist \link in general.

wthrowe · 2023-11-21T02:41:18Z

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

   * Evaluate the transition function at the Cartesian coordinates divided by
-   * the radius. All divisions over the radius which could potentially be zero
-   * divisions are left to the transition function which uses its knowledge
-   * of the domain to do this safely or throw an appropriate error.
+   * the radius to a power. All divisions over the radius which could
+   * potentially be zero divisions are left to the transition function which
+   * uses its knowledge of the domain to do this safely or throw an appropriate
+   * error.


The general feature isn't new, but it's being significantly expanded upon so I want to understand it.

What is this for? Looking at the math, the place a division-by-zero could potentially be analytically handled is in the shape map, not the transition function (and I don't think it actually gives a usable map even then). All the implementations are just dividing by r, without any of the checks mentioned here.

Yes, there is one place in the shape map where dividing by r could be taken care of, but my best guess as to why the transition function is responsible for checking the 1/r is because the transition function could be anything and so it may also need to account for a 1/r somewhere itself. And, I think, because the only place that 1/r shows up is with the transition function itself, it made sense to just have the transition function do all of the checks.

Both the Wedge and Sphere transition functions do check the magnitude of the point (in debug mode) before dividing by it. So these checks do exist.

You're right that I missed the checks in the transition functions, but why not move them to the shape map itself to avoid duplication? The shape map isn't valid at r=0, so there can't be any issues with dividing by r.

Ok yeah, I'm convinced this should be in the shape map now. If the transition function itself needs to divide by r it can just assume it's not zero since the shape map took care of it.

If we want to consider using the map in the opposite sense, to remove distortion when going inwards, then you could get a valid map at the origin where I think you might want to remove cancellations in the transition function, but I still don't think you'd want to do it with varying powers in the fetch functions. You'd want the interface to only provide the valid quantities.

wthrowe · 2023-11-21T02:50:00Z

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

+  } else {
+    return (outer_distance - radius) /
+           ((outer_distance - inner_surface_.distance(rotated_coords))) *
+           integer_pow(1.0 / radius, power);


You're already doing a division. Put the new factors inside the denominator instead of doing another. (Particularly since I'm not sure whether integer_pow(1.0 / radius, 3) does one division or three.)

I changed this in the commit after, (no clue why) but I'll move it to the shape map edit commit

wthrowe · 2023-11-21T02:58:36Z

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

-  const auto transition_func_over_radius =
-      transition_func_->map_over_radius(centered_coords);
+  transition_func_over_square_radius =
+      transition_func_->map_over_radius(centered_coords, 2);


This re-use is not an optimization. It doesn't avoid an allocation and may prevent some form of copy elision.

wthrowe · 2023-11-21T02:59:09Z

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

+  auto& transition_func_over_cube_radius = get<1>(theta_phis);
+
+  transition_func_over_cube_radius =
+      transition_func_->map_over_radius(centered_coords, 3);


Again, not an optimization.

wthrowe · 2023-11-21T03:02:19Z

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

+  // re-use allocation
+  auto& transition_func_over_radius = get<1>(theta_phis);
+  transition_func_over_radius =
+      transition_func_->map_over_radius(centered_coords);


Again, not an optimization.

wthrowe · 2023-11-21T03:10:56Z

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

+            "The Wedge transition map was called with coordinates outside of "
+            "its wedge. The "
+            << (axis.value() == 0_st ? "x" : (axis.value() == 1_st ? "y" : "z"))
+            << "-coordinate shouldn't be 0 but it is.");


Print the coordinate or the index or something.

knelli2 · 2023-11-22T18:30:17Z

@wthrowe Pushed some fixups. I'll make the whitelist commit the first commit when I squash. I'll also squash the new changes into the appropriate commits as well since I edited the Wedge transition in both original commits.

wthrowe

Looks good. Squash (and fix clang-tidy).

wthrowe · 2023-11-27T18:22:41Z

src/Domain/CoordinateMaps/TimeDependent/Shape.cpp

+T Shape::check_and_compute_one_over_radius(
+    const std::array<T, 3>& centered_coords) const {
+#ifdef SPECTRE_DEBUG
+  const T radius = magnitude(centered_coords);


[optional] Move up one line and return 1.0 / radius.

wthrowe

Looks good. GitHub Actions is having problems (in general, not just for us) at the moment, but they usually resolve that sort of thing quickly.

knelli2 added the priority critical for progress label Nov 21, 2023

knelli2 added this to the BBH inspiral reaches common apparent horizon milestone Nov 21, 2023

knelli2 commented Nov 21, 2023

View reviewed changes

wthrowe requested changes Nov 21, 2023

View reviewed changes

wthrowe reviewed Nov 27, 2023

View reviewed changes

knelli2 added 2 commits November 27, 2023 11:03

Add \link to whitelist of long lines

d84c420

Fix failing shape map test

d318838

knelli2 force-pushed the spec_shape_map branch from 49387dc to b568435 Compare November 27, 2023 20:23

wthrowe previously approved these changes Nov 27, 2023

View reviewed changes

knelli2 added 2 commits November 27, 2023 14:10

Put one over r factor into shape map, not transition funcs

9a78117

Rewrite wedge transition cleaner

4f801a5

knelli2 dismissed wthrowe’s stale review via 4f801a5 November 27, 2023 22:11

knelli2 force-pushed the spec_shape_map branch from b568435 to 4f801a5 Compare November 27, 2023 22:11

wthrowe approved these changes Nov 28, 2023

View reviewed changes

wthrowe merged commit 64dda8f into sxs-collaboration:develop Nov 28, 2023
21 of 22 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Make shape map exactly like the one in SpEC #5636

Make shape map exactly like the one in SpEC #5636

knelli2 commented Nov 21, 2023

knelli2 Nov 21, 2023

wthrowe Nov 21, 2023

knelli2 Nov 21, 2023

wthrowe Nov 21, 2023

wthrowe Nov 21, 2023

knelli2 Nov 21, 2023

wthrowe Nov 21, 2023

knelli2 Nov 21, 2023

wthrowe Nov 21, 2023

wthrowe Nov 21, 2023

knelli2 Nov 21, 2023

wthrowe Nov 21, 2023

wthrowe Nov 21, 2023

wthrowe Nov 21, 2023

wthrowe Nov 21, 2023

knelli2 commented Nov 22, 2023

wthrowe left a comment

wthrowe Nov 27, 2023

wthrowe left a comment

Make shape map exactly like the one in SpEC #5636

Make shape map exactly like the one in SpEC #5636

Conversation

knelli2 commented Nov 21, 2023

Proposed changes

Upgrade instructions

Code review checklist

Further comments

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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

knelli2 commented Nov 22, 2023

wthrowe left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wthrowe left a comment

Choose a reason for hiding this comment