fixup. blaze stuff

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

@@ -7,6 +7,7 @@
 #include <array>
 #include <optional>
 #include <pup.h>
+#include <type_traits>
 #include <vector>
 
 #include "DataStructures/Blaze/IntegerPow.hpp"
@@ -19,6 +20,7 @@
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeString.hpp"
 #include "Utilities/Math.hpp"
+#include "Utilities/StdArrayHelpers.hpp"
 #include "Utilities/StdHelpers.hpp"
 
 namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
@@ -93,13 +95,41 @@ DataVector SphereTransition::operator()(
     const std::array<DataVector, 3>& source_coords,
     const std::optional<size_t>& one_over_radius_power) const {
   DataVector result = source_coords[0];
-  // Go point by point
+#ifdef SPECTRE_DEBUG
+  // Do checks point by point
   for (size_t i = 0; i < source_coords[0].size(); i++) {
     result[i] = (*this)(std::array{source_coords[0][i], source_coords[1][i],
                                    source_coords[2][i]},
                         one_over_radius_power);
   }
+#endif
+
+// If we are in debug, just return what we already computed. Otherwise compute
+// using blaze for better performance
+#ifdef SPECTRE_DEBUG
   return result;
+#else
+  const DataVector mag = magnitude(source_coords);
+
+  if (interior_) {
+    return inverse_cube_r_min_ *
+           (one_over_radius_power.value_or(0_st) < 3
+                ? DataVector{integer_pow(
+                      mag, static_cast<int>(
+                               2 - one_over_radius_power.value_or(0_st)))}
+                : 1.0 /
+                      integer_pow(mag, static_cast<int>(
+                                           one_over_radius_power.value() - 2)));
+  }
+
+  result = (a_ * mag + b_);
+  // Avoid roundoff
+  result = blaze::clamp(result, 0.0, 1.0);
+
+  return result /
+         integer_pow(
+             mag, static_cast<int>(1 + one_over_radius_power.value_or(0_st)));
+#endif
 }
 
 std::optional<double> SphereTransition::original_radius_over_radius(
@@ -154,7 +184,9 @@ std::optional<double> SphereTransition::original_radius_over_radius(
 
         return A + B;
       }
-    } else if (equal_within_roundoff(mag, r_min_)) {
+    } else if (equal_within_roundoff(mag + radial_distortion, r_min_)) {
+      // IntMid: Since the point is on the boundary, we can't tell if this is
+      // the interior or not, but either map will work.
       return std::optional{r_min_ / (r_min_ - radial_distortion)};
     }
   }
@@ -230,15 +262,30 @@ std::array<double, 3> SphereTransition::gradient(
 }
 std::array<DataVector, 3> SphereTransition::gradient(
     const std::array<DataVector, 3>& source_coords) const {
+#ifdef SPECTRE_DEBUG
   auto result = source_coords;
+  // Do checks point by point
   for (size_t i = 0; i < source_coords[0].size(); i++) {
     auto double_result = gradient(std::array{
         source_coords[0][i], source_coords[1][i], source_coords[2][i]});
     for (size_t j = 0; j < 3; j++) {
       gsl::at(result, j)[i] = gsl::at(double_result, j);
     }
   }
+#endif
+
+// If we are in debug, just return what we already computed. Otherwise compute
+// using blaze for better performance
+#ifdef SPECTRE_DEBUG
   return result;
+#else
+  if (interior_) {
+    return 2.0 * inverse_cube_r_min_ * source_coords;
+  } else {
+    return source_coords * (a_ / dot(source_coords, source_coords) -
+                            (*this)(source_coords, {2}));
+  }
+#endif
 }
 
 bool SphereTransition::operator==(

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

@@ -13,6 +13,7 @@
 #include <type_traits>
 #include <vector>
 
+#include "DataStructures/Blaze/IntegerPow.hpp"
 #include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
@@ -413,14 +414,66 @@ double Wedge::operator()(
 DataVector Wedge::operator()(
     const std::array<DataVector, 3>& source_coords,
     const std::optional<size_t>& one_over_radius_power) const {
+#ifdef SPECTRE_DEBUG
   DataVector result = source_coords[0];
-  // Go point by point to avoid the center point
+  // Do checks point by point
   for (size_t i = 0; i < source_coords[0].size(); i++) {
     result[i] = (*this)(std::array{source_coords[0][i], source_coords[1][i],
                                    source_coords[2][i]},
                         one_over_radius_power);
   }
+#endif
+
+// If we are in debug, just return what we already computed. Otherwise compute
+// using blaze for better performance
+#ifdef SPECTRE_DEBUG
   return result;
+#else
+  const DataVector centered_coords_magnitude = magnitude(source_coords);
+
+  if (axis_ == Axis::Interior) {
+    return 1.0 / cube(inner_surface_.radius) *
+           (one_over_radius_power.value_or(0_st) < 3
+                ? DataVector{integer_pow(
+                      centered_coords_magnitude,
+                      static_cast<int>(2 -
+                                       one_over_radius_power.value_or(0_st)))}
+                : 1.0 / integer_pow(centered_coords_magnitude,
+                                    static_cast<int>(
+                                        one_over_radius_power.value() - 2)));
+  }
+
+  const DataVector lambda =
+      compute_lambda(source_coords, centered_coords_magnitude, {axis_});
+
+  const DataVector inner_distance =
+      inner_surface_.sphericity == 1.0
+          ? make_with_value<DataVector>(centered_coords_magnitude,
+                                        inner_surface_.radius)
+          : magnitude(compute_inner_surface_vector(
+                source_coords, centered_coords_magnitude, {axis_}));
+
+  const std::array<DataVector, 3> outer_surface_vector =
+      compute_outer_surface_vector(source_coords, lambda);
+  const DataVector outer_distance = magnitude(outer_surface_vector);
+
+  DataVector linear_transition_func =
+      (outer_distance - centered_coords_magnitude) /
+      (outer_distance - inner_distance);
+
+  if (reverse_) {
+    linear_transition_func = 1.0 - linear_transition_func;
+  }
+
+  // Accounts for roundoff
+  linear_transition_func = blaze::clamp(linear_transition_func, 0.0, 1.0);
+
+  return linear_transition_func /
+         integer_pow(
+             centered_coords_magnitude,
+             static_cast<int>(1 + one_over_radius_power.value_or(0_st)));
+
+#endif
 }
 
 std::optional<double> Wedge::original_radius_over_radius(
@@ -580,56 +633,84 @@ std::optional<double> Wedge::original_radius_over_radius(
 
 std::array<double, 3> Wedge::gradient(
     const std::array<double, 3>& source_coords) const {
+  return gradient_impl<double>(source_coords);
+}
+std::array<DataVector, 3> Wedge::gradient(
+    const std::array<DataVector, 3>& source_coords) const {
+  return gradient_impl<DataVector>(source_coords);
+}
+
+template <typename T>
+std::array<T, 3> Wedge::gradient_impl(
+    const std::array<T, 3>& source_coords) const {
   // The source coords are centered
-  const double centered_coords_magnitude = magnitude(source_coords);
+  const T centered_coords_magnitude = magnitude(source_coords);
 
   // Short circuit if we are in the interior
   if (axis_ == Axis::Interior) {
     return 2.0 / cube(inner_surface_.radius) * source_coords;
   }
 
-  // First check if we are at the inner center to avoid dividing by zero below
-  if (UNLIKELY(equal_within_roundoff(centered_coords_magnitude, 0.0))) {
-    ERROR(
-        "The gradient of the wedge transition was called with a point that has "
-        "zero centered radius, but the axis isn't the 'Interior'. Point is "
-        << source_coords << ", axis is " << axis_);
+// First check if we are at the inner center to avoid dividing by zero below
+#ifdef SPECTRE_DEBUG
+  for (size_t i = 0; i < get_size(centered_coords_magnitude); i++) {
+    if (UNLIKELY(equal_within_roundoff(
+            get_element(centered_coords_magnitude, i), 0.0))) {
+      ERROR(
+          "The gradient of the wedge transition was called with a point that "
+          "has "
+          "zero centered radius, but the axis isn't the 'Interior'. Point is "
+          << inner_surface_.center << ", axis is " << axis_);
+    }
   }
+#endif
 
-  const double one_over_centered_coords_magnitude =
-      1.0 / centered_coords_magnitude;
+  const T one_over_centered_coords_magnitude = 1.0 / centered_coords_magnitude;
 
-  const double lambda =
+  const T lambda =
       compute_lambda(source_coords, centered_coords_magnitude, {axis_});
 
-  const double inner_distance =
+  const T inner_distance =
       inner_surface_.sphericity == 1.0
-          ? make_with_value<double>(lambda, inner_surface_.radius)
+          ? make_with_value<T>(lambda, inner_surface_.radius)
           : magnitude(compute_inner_surface_vector(
                 source_coords, centered_coords_magnitude, {axis_}));
 
-  const std::array<double, 3> outer_surface_vector =
+  const std::array<T, 3> outer_surface_vector =
       compute_outer_surface_vector(source_coords, lambda);
-  const double outer_distance = magnitude(outer_surface_vector);
+  const T outer_distance = magnitude(outer_surface_vector);
 
-  // Check if the point we were passed is within the transition region
-  if (centered_coords_magnitude < (1.0 - eps_) * inner_distance) {
-    ERROR("Wedge transition called with centered coordinate "
-          << source_coords << " (with radius " << centered_coords_magnitude
-          << ") which is inside the inner surface (" << inner_distance
-          << ") while the axis was not 'Interior' (" << axis_ << ")");
-  } else if (centered_coords_magnitude > (1.0 + eps_) * outer_distance) {
-    ERROR("Wedge transition called with centered coordinate "
-          << source_coords << " (with radius " << centered_coords_magnitude
-          << ") which is outside the outer surface (" << outer_distance
-          << ").");
+// Check if the point we were passed is within the transition region
+#ifdef SPECTRE_DEBUG
+  for (size_t i = 0; i < get_size(centered_coords_magnitude); i++) {
+    if (get_element(centered_coords_magnitude, i) <
+        (1.0 - eps_) * get_element(inner_distance, i)) {
+      ERROR("Wedge transition called with centered coordinate "
+            << (std::array{get_element(source_coords[0], i),
+                           get_element(source_coords[1], i),
+                           get_element(source_coords[2], i)})
+            << " (with radius " << get_element(centered_coords_magnitude, i)
+            << ") which is inside the inner surface ("
+            << get_element(inner_distance, i)
+            << ") while the axis was not 'Interior' (" << axis_ << ")");
+    } else if (get_element(centered_coords_magnitude, i) >
+               (1.0 + eps_) * get_element(outer_distance, i)) {
+      ERROR("Wedge transition called with centered coordinate "
+            << (std::array{get_element(source_coords[0], i),
+                           get_element(source_coords[1], i),
+                           get_element(source_coords[2], i)})
+            << " (with radius " << get_element(centered_coords_magnitude, i)
+            << ") which is outside the outer surface ("
+            << get_element(outer_distance, i) << ").");
+    }
   }
+#endif
 
   // This can only be called if the projection center is 0, otherwise this
   // formula won't work. And if the projection center isn't 0, then we require
   // that the sphericity of the inner surface is 1, so we ASSERT that here
   // as well
-  const auto inner_surface_gradient = [&]() -> std::array<double, 3> {
+  const auto inner_surface_gradient = [&]() -> std::array<T, 3> {
     using ::operator<<;
     ASSERT((projection_center_ == std::array{0.0, 0.0, 0.0}),
            "Should not be calculating the inner surface gradient when the "
@@ -643,9 +724,9 @@ std::array<double, 3> Wedge::gradient(
     const size_t axis_plus_one = (axis_idx + 1) % 3;
     const size_t axis_plus_two = (axis_idx + 2) % 3;
 
-    const double& axis_coord = gsl::at(source_coords, axis_idx);
+    const T& axis_coord = gsl::at(source_coords, axis_idx);
 
-    std::array<double, 3> grad = make_array<3>(
+    std::array<T, 3> grad = make_array<3, T>(
         inner_surface_.radius * (1.0 - inner_surface_.sphericity) / sqrt(3.0) *
         one_over_centered_coords_magnitude / abs(axis_coord));
 
@@ -661,12 +742,12 @@ std::array<double, 3> Wedge::gradient(
 
   // We always need to compute the outer gradient regardless of the projection
   // center or the outer sphericity
-  const std::array<double, 3> outer_gradient =
+  const std::array<T, 3> outer_gradient =
       lambda * source_coords * one_over_centered_coords_magnitude +
       compute_lambda_gradient(source_coords, centered_coords_magnitude) *
           centered_coords_magnitude;
 
-  const double one_over_distance_difference =
+  const T one_over_distance_difference =
       1.0 / (outer_distance - inner_distance);
 
   // Avoid allocating an array of 0 if the inner surface is a sphere
@@ -679,7 +760,7 @@ std::array<double, 3> Wedge::gradient(
                one_over_centered_coords_magnitude * (reverse_ ? -1.0 : 1.0) -
            source_coords * (*this)(source_coords, {2});
   } else {
-    const std::array<double, 3> inner_gradient = inner_surface_gradient();
+    const std::array<T, 3> inner_gradient = inner_surface_gradient();
 
     return (outer_gradient -
             source_coords * one_over_centered_coords_magnitude -
@@ -692,19 +773,6 @@ std::array<double, 3> Wedge::gradient(
   }
 }
 
-std::array<DataVector, 3> Wedge::gradient(
-    const std::array<DataVector, 3>& source_coords) const {
-  auto result = source_coords;
-  for (size_t i = 0; i < source_coords[0].size(); i++) {
-    const std::array<double, 3> double_result = gradient(std::array{
-        source_coords[0][i], source_coords[1][i], source_coords[2][i]});
-    for (size_t j = 0; j < 3; j++) {
-      gsl::at(result, j)[i] = gsl::at(double_result, j);
-    }
-  }
-  return result;
-}
-
 bool Wedge::operator==(const ShapeMapTransitionFunction& other) const {
   if (dynamic_cast<const Wedge*>(&other) == nullptr) {
     return false;

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

@@ -505,8 +505,7 @@ class Wedge final : public ShapeMapTransitionFunction {
 
  private:
   template <typename T>
-  T call_impl(const std::array<T, 3>& source_coords,
-              const std::optional<size_t>& one_over_radius_power) const;
+  std::array<T, 3> gradient_impl(const std::array<T, 3>& source_coords) const;
 
   // This is x_0 - P in the docs
   template <typename T>