Merge pull request sxs-collaboration#5543 from acifajkya/CE_Norm_Gh

Add constraint energy normalization in GH

src/Evolution/Systems/GeneralizedHarmonic/Constraints.cpp

@@ -10,6 +10,7 @@
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/Christoffel.hpp"
 #include "Utilities/ConstantExpressions.hpp"
+#include "Utilities/ContainerHelpers.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeWithValue.hpp"
@@ -1325,6 +1326,60 @@ void constraint_energy(
     }
   }
 }
+
+template <typename DataType, size_t SpatialDim, typename Frame>
+Scalar<DataType> constraint_energy_normalization(
+    const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
+    const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
+    const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
+    const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
+    const Scalar<DataType>& sqrt_spatial_metric_determinant,
+    const double dimensional_constant) {
+  Scalar<DataType> energy_norm{get_size(get(sqrt_spatial_metric_determinant))};
+  constraint_energy_normalization<DataType, SpatialDim, Frame>(
+      make_not_null(&energy_norm), d_spacetime_metric, d_pi, d_phi,
+      inverse_spatial_metric, sqrt_spatial_metric_determinant,
+      dimensional_constant);
+  return energy_norm;
+}
+
+template <typename DataType, size_t SpatialDim, typename Frame>
+void constraint_energy_normalization(
+    const gsl::not_null<Scalar<DataType>*> energy_norm,
+    const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
+    const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
+    const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
+    const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
+    const Scalar<DataType>& sqrt_spatial_metric_determinant,
+    const double dimensional_constant) {
+  const double square_dimensional_constant = square(dimensional_constant);
+  get(*energy_norm) = 0.0;
+  for (size_t a = 0; a < SpatialDim + 1; ++a) {
+    for (size_t c = 0; c < SpatialDim + 1; ++c) {
+      for (size_t i = 0; i < SpatialDim; ++i) {
+        for (size_t j = 0; j < SpatialDim; ++j) {
+          get(*energy_norm) += inverse_spatial_metric.get(i, j) *
+                                   square_dimensional_constant *
+                                   d_spacetime_metric.get(i, a, c) *
+                                   d_spacetime_metric.get(j, a, c) +
+                               inverse_spatial_metric.get(i, j) *
+                                   d_pi.get(i, a, c) * d_pi.get(j, a, c);
+
+          for (size_t k = 0; k < SpatialDim; ++k) {
+            for (size_t l = 0; l < SpatialDim; ++l) {
+              get(*energy_norm) += inverse_spatial_metric.get(i, j) *
+                                   inverse_spatial_metric.get(k, l) *
+                                   d_phi.get(i, k, a, c) *
+                                   d_phi.get(j, l, a, c);
+            }
+          }
+        }
+      }
+    }
+  }
+  get(*energy_norm) *= get(sqrt_spatial_metric_determinant);
+}
+
 }  // namespace gh
 
 // Explicit Instantiations
@@ -1525,7 +1580,26 @@ void constraint_energy(
       double gauge_constraint_multiplier,                                     \
       double two_index_constraint_multiplier,                                 \
       double three_index_constraint_multiplier,                               \
-      double four_index_constraint_multiplier);
+      double four_index_constraint_multiplier);                               \
+  template Scalar<DTYPE(data)> gh::constraint_energy_normalization(           \
+      const tnsr::iaa<DTYPE(data), DIM(data), FRAME(data)>&                   \
+          d_spacetime_metric,                                                 \
+      const tnsr::iaa<DTYPE(data), DIM(data), FRAME(data)>& d_pi,             \
+      const tnsr::ijaa<DTYPE(data), DIM(data), FRAME(data)>& d_phi,           \
+      const tnsr::II<DTYPE(data), DIM(data), FRAME(data)>&                    \
+          inverse_spatial_metric,                                             \
+      const Scalar<DTYPE(data)>& sqrt_spatial_metric_determinant,             \
+      double dimensional_constant);                                           \
+  template void gh::constraint_energy_normalization(                          \
+      const gsl::not_null<Scalar<DTYPE(data)>*> energy_norm,                  \
+      const tnsr::iaa<DTYPE(data), DIM(data), FRAME(data)>&                   \
+          d_spacetime_metric,                                                 \
+      const tnsr::iaa<DTYPE(data), DIM(data), FRAME(data)>& d_pi,             \
+      const tnsr::ijaa<DTYPE(data), DIM(data), FRAME(data)>& d_phi,           \
+      const tnsr::II<DTYPE(data), DIM(data), FRAME(data)>&                    \
+          inverse_spatial_metric,                                             \
+      const Scalar<DTYPE(data)>& sqrt_spatial_metric_determinant,             \
+      double dimensional_constant);
 
 #define INSTANTIATE_ONLY_3D(_, data)                                       \
   template tnsr::iaa<DTYPE(data), DIM(data), FRAME(data)>                  \

src/Evolution/Systems/GeneralizedHarmonic/Constraints.hpp

@@ -551,6 +551,44 @@ void constraint_energy(
     double four_index_constraint_multiplier = 1.0);
 /// @}
 
+/// @{
+/*!
+ * \brief Computes the generalized-harmonic normalized constraint energy.
+ *
+ * \details Computes the generalized-harmonic normalized constraint energy
+ * integrand [Eq. (70) of \cite Lindblom2005qh with \f$m^{ab}=\delta^{ab}\f$],
+ * \f{eqnarray}{
+ * ||E||=\gamma^{ij}\delta^{ab}\delta^{cd}
+ * (\Lambda^{2}\partial_{i}g_{ac}\partial_{j}g_{bd}
+ * + \partial_{i}\Pi_{ac}\partial_{j}\Pi_{bd} +
+ * \gamma^{kl}\partial_{i}\Phi_{kac}\partial_{j}\Phi_{lbd})\gamma^{1/2}
+ * \f}
+ *
+ * Here \f$\gamma^{ij}\f$ is the inverse spatial metric, and \f$\Lambda^{2}\f$
+ * is the squared dimensional constant. Note
+ * that in this equation, spacetime indices \f$a,b\f$ are raised and lowered
+ * with the Kronecker delta.
+ */
+template <typename DataType, size_t SpatialDim, typename Frame>
+Scalar<DataType> constraint_energy_normalization(
+    const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
+    const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
+    const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
+    const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
+    const Scalar<DataType>& sqrt_spatial_metric_determinant,
+    double dimensional_constant);
+
+template <typename DataType, size_t SpatialDim, typename Frame>
+void constraint_energy_normalization(
+    gsl::not_null<Scalar<DataType>*> energy_norm,
+    const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
+    const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
+    const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
+    const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
+    const Scalar<DataType>& sqrt_spatial_metric_determinant,
+    double dimensional_constant);
+/// @}
+
 namespace Tags {
 /*!
  * \brief Compute item to get the gauge constraint for the

tests/Unit/Evolution/Systems/GeneralizedHarmonic/TestFunctions.py

@@ -1095,3 +1095,29 @@ def constraint_energy(
 
 
 # End test functions for constraint energy
+
+# Begin test function for normalized constraint energy
+
+
+def constraint_energy_normalization(
+    d_spacetime_metric,
+    d_pi,
+    d_phi,
+    inverse_spatial_metric,
+    sqrt_spatial_metric_determinant,
+    dimensional_constant,
+):
+    energy_norm = (
+        dimensional_constant**2
+        * np.einsum("iac,jac", d_spacetime_metric, d_spacetime_metric)
+        + np.einsum("iac,jac", d_pi, d_pi)
+        + np.einsum("kl,ikac,jlac", inverse_spatial_metric, d_phi, d_phi)
+    )
+    energy_norm = (
+        np.einsum("ij,ij", inverse_spatial_metric, energy_norm)
+        * sqrt_spatial_metric_determinant
+    )
+    return energy_norm
+
+
+# End test functions for normalized constraint energy

tests/Unit/Evolution/Systems/GeneralizedHarmonic/Test_Constraints.cpp

@@ -748,6 +748,52 @@ void test_constraint_energy_analytic(const Solution& solution,
       numerical_approx);
 }
 
+// Test the return-by-value normalized constraint energy function using random
+// values
+template <typename DataType, size_t SpatialDim, typename Frame>
+void test_constraint_energy_normalization_random(
+    const DataType& used_for_size) {
+  pypp::check_with_random_values<1>(
+      static_cast<Scalar<DataType> (*)(
+          const tnsr::iaa<DataType, SpatialDim, Frame>&,
+          const tnsr::iaa<DataType, SpatialDim, Frame>&,
+          const tnsr::ijaa<DataType, SpatialDim, Frame>&,
+          const tnsr::II<DataType, SpatialDim, Frame>&, const Scalar<DataType>&,
+          double)>(
+          &gh::constraint_energy_normalization<DataType, SpatialDim, Frame>),
+      "TestFunctions", "constraint_energy_normalization", {{{-1.0, 1.0}}},
+      used_for_size);
+}
+
+// Test the return-by-value normalized constraint energy function for Minkowski
+void test_constraint_energy_normalization_minkowski(
+    const DataVector& used_for_size) {
+  const auto d_spacetime_metric =
+      make_with_value<tnsr::iaa<DataVector, 3, Frame::Inertial>>(used_for_size,
+                                                                 0.0);
+  const auto d_pi = make_with_value<tnsr::iaa<DataVector, 3, Frame::Inertial>>(
+      used_for_size, 0.0);
+  const auto d_phi =
+      make_with_value<tnsr::ijaa<DataVector, 3, Frame::Inertial>>(used_for_size,
+                                                                  0.0);
+  auto inverse_spatial_metric =
+      make_with_value<tnsr::II<DataVector, 3, Frame::Inertial>>(used_for_size,
+                                                                0.0);
+  for (size_t i = 0; i < 3; ++i) {
+    inverse_spatial_metric.get(i, i) = 1.0;
+  }
+  const auto sqrt_spatial_metric_determinant =
+      make_with_value<Scalar<DataVector>>(used_for_size, 1.0);
+  const double dimensional_constant = 3.0;
+
+  // type result = func(args...);
+  const Scalar<DataVector> result = gh::constraint_energy_normalization(
+      d_spacetime_metric, d_pi, d_phi, inverse_spatial_metric,
+      sqrt_spatial_metric_determinant, dimensional_constant);
+
+  CHECK(get(result) == 0.0);
+}
+
 // Test compute items for various constraints via insertion and retrieval
 // in a databox
 template <typename Solution>
@@ -1287,6 +1333,36 @@ void constraint_energy() {
       std::numeric_limits<double>::signaling_NaN());
 }
 
+void constraint_energy_normalization() {
+  // Test the normalized constraint energy in minkowski
+  test_constraint_energy_normalization_minkowski(
+      DataVector(4, std::numeric_limits<double>::signaling_NaN()));
+
+  // Test the normalized constraint energy with random numbers
+  test_constraint_energy_normalization_random<DataVector, 1, Frame::Grid>(
+      DataVector(4, std::numeric_limits<double>::signaling_NaN()));
+  test_constraint_energy_normalization_random<DataVector, 1, Frame::Inertial>(
+      DataVector(4, std::numeric_limits<double>::signaling_NaN()));
+  test_constraint_energy_normalization_random<double, 1, Frame::Grid>(
+      std::numeric_limits<double>::signaling_NaN());
+  test_constraint_energy_normalization_random<double, 1, Frame::Inertial>(
+      std::numeric_limits<double>::signaling_NaN());
+
+  test_constraint_energy_normalization_random<DataVector, 2, Frame::Inertial>(
+      DataVector(4, std::numeric_limits<double>::signaling_NaN()));
+  test_constraint_energy_normalization_random<double, 2, Frame::Grid>(
+      std::numeric_limits<double>::signaling_NaN());
+  test_constraint_energy_normalization_random<double, 2, Frame::Inertial>(
+      std::numeric_limits<double>::signaling_NaN());
+
+  test_constraint_energy_normalization_random<DataVector, 3, Frame::Inertial>(
+      DataVector(4, std::numeric_limits<double>::signaling_NaN()));
+  test_constraint_energy_normalization_random<double, 3, Frame::Grid>(
+      std::numeric_limits<double>::signaling_NaN());
+  test_constraint_energy_normalization_random<double, 3, Frame::Inertial>(
+      std::numeric_limits<double>::signaling_NaN());
+}
+
 void test_compute_tags() {
   // Test the F constraint against Kerr Schild
   const double mass = 1.4;
@@ -1313,5 +1389,6 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.GeneralizedHarmonic.Constraints",
   four_index_constraint();
   f_constraint();
   constraint_energy();
+  constraint_energy_normalization();
   test_compute_tags();
 }