Merge pull request sxs-collaboration#6404 from nikwit/worldtube-gammas

Add compute tags for constraint gammas in worldtube evolution
kidder · Dec 13, 2024 · b6fc141 · b6fc141
2 parents cbe1362 + c55a6ed
commit b6fc141
Show file tree

Hide file tree

Showing 4 changed files with 97 additions and 2 deletions.
diff --git a/src/Evolution/Executables/CurvedScalarWave/EvolveWorldtubeCurvedScalarWave.hpp b/src/Evolution/Executables/CurvedScalarWave/EvolveWorldtubeCurvedScalarWave.hpp
@@ -292,8 +292,6 @@ struct EvolutionMetavars {
       Initialization::Actions::AddSimpleTags<
           CurvedScalarWave::Worldtube::Initialization::
               InitializeCurrentIteration,
-          CurvedScalarWave::Worldtube::Initialization::
-              InitializeConstraintDampingGammas<volume_dim>,
           CurvedScalarWave::Initialization::InitializeEvolvedVariables<
               volume_dim>>,
       Initialization::Actions::AddComputeTags<
@@ -302,6 +300,8 @@ struct EvolutionMetavars {
       Initialization::Actions::AddComputeTags<tmpl::list<
           CurvedScalarWave::Worldtube::Tags::ParticlePositionVelocityCompute<
               volume_dim>,
+          CurvedScalarWave::Worldtube::Tags::ConstraintGamma1Compute,
+          CurvedScalarWave::Worldtube::Tags::ConstraintGamma2Compute,
           CurvedScalarWave::Worldtube::Tags::FaceCoordinatesCompute<
               volume_dim, Frame::Inertial, true>,
           CurvedScalarWave::Worldtube::Tags::FaceCoordinatesCompute<

diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.cpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.cpp
@@ -319,6 +319,31 @@ void FaceQuantitiesCompute::function(
       face_inv_jacobian, direction.value());
 }
 
+void ConstraintGamma1Compute::function(
+    gsl::not_null<Scalar<DataVector>*> gamma1,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& coords) {
+  get(*gamma1).destructive_resize(get<0>(coords).size());
+  get(*gamma1) = 0.;
+}
+
+void ConstraintGamma2Compute::function(
+    gsl::not_null<Scalar<DataVector>*> gamma2,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& coords,
+    const std::array<tnsr::I<double, 3, Frame::Inertial>, 2>& pos_vel) {
+  get(*gamma2).destructive_resize(get<0>(coords).size());
+  const auto& pos = pos_vel[0];
+
+  // re-use allocation
+  auto& centered_radii = get(*gamma2);
+  centered_radii = sqrt(square(get<0>(coords) - get<0>(pos)) +
+                        square(get<1>(coords) - get<1>(pos)) +
+                        square(get<2>(coords) - get<2>(pos)));
+  const double amplitude = 10.;
+  const double sigma = 1e-1;
+  const double constant = 1e-3;
+  get(*gamma2) = amplitude * exp(-square(sigma * centered_radii)) + constant;
+}
+
 template struct BackgroundQuantitiesCompute<3>;
 template struct EvolvedParticlePositionVelocityCompute<3>;
 template struct GeodesicAccelerationCompute<3>;

diff --git a/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp b/src/Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp
@@ -839,5 +839,40 @@ struct dtPsi0 : db::SimpleTag {
   using type = Scalar<DataVector>;
 };
 
+/*!
+ * \brief Sets Gamma1 to zero throughout the domain. The equations are given in
+ * Initialization::InitializeConstraintDampingGammas.
+ */
+struct ConstraintGamma1Compute : CurvedScalarWave::Tags::ConstraintGamma1,
+                                 db::ComputeTag {
+  static constexpr size_t Dim = 3;
+  using base = CurvedScalarWave::Tags::ConstraintGamma1;
+  using return_type = Scalar<DataVector>;
+  using argument_tags =
+      tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>>;
+  static void function(gsl::not_null<Scalar<DataVector>*> gamma1,
+                       const tnsr::I<DataVector, Dim, Frame::Inertial>& coords);
+};
+
+/*!
+ * \brief Sets Gamma2 to a Gaussian that falls off to a constant value centered
+ * on the position of the particle. This was found to be necessary for a stable
+ * evolution. The equations are given in
+ * Initialization::InitializeConstraintDampingGammas.
+ */
+struct ConstraintGamma2Compute : CurvedScalarWave::Tags::ConstraintGamma2,
+                                 db::ComputeTag {
+  static constexpr size_t Dim = 3;
+  using base = CurvedScalarWave::Tags::ConstraintGamma2;
+  using return_type = Scalar<DataVector>;
+  using argument_tags =
+      tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
+                 ParticlePositionVelocity<Dim>>;
+  static void function(
+      gsl::not_null<Scalar<DataVector>*> gamma2,
+      const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
+      const std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>& pos_vel);
+};
+
 }  // namespace Tags
 }  // namespace CurvedScalarWave::Worldtube
diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Tags.cpp
@@ -381,6 +381,40 @@ void test_evolved_particle_position_velocity_compute() {
   }
 }
 
+void test_constraint_gammas_compute() {
+  static constexpr size_t Dim = 3;
+  MAKE_GENERATOR(gen);
+  const std::uniform_real_distribution<> dist(-1., 1.);
+  const auto pos = make_with_random_values<tnsr::I<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const auto vel = make_with_random_values<tnsr::I<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+
+  const size_t coord_size = 100;
+  const auto coords = make_with_random_values<tnsr::I<DataVector, Dim>>(
+      make_not_null(&gen), dist, DataVector(coord_size));
+
+  auto box = db::create<
+      db::AddSimpleTags<Tags::ParticlePositionVelocity<Dim>,
+                        domain::Tags::Coordinates<Dim, Frame::Inertial>>,
+      db::AddComputeTags<Tags::ConstraintGamma1Compute,
+                         Tags::ConstraintGamma2Compute>>(
+      std::array<tnsr::I<double, Dim>, 2>{pos, vel}, coords);
+
+  const auto& gamma1 = db::get<CurvedScalarWave::Tags::ConstraintGamma1>(box);
+  CHECK(get(gamma1) == DataVector(coord_size, 0.));
+
+  auto centered_coords = coords;
+  for (size_t i = 0; i < 3; ++i) {
+    centered_coords.get(i) -= pos.get(i);
+  }
+  const auto centered_radii = magnitude(centered_coords);
+  const auto gamma2_expected =
+      10. * exp(-square(0.1 * get(centered_radii))) + 1e-3;
+  const auto& gamma2 = db::get<CurvedScalarWave::Tags::ConstraintGamma2>(box);
+  CHECK_ITERABLE_APPROX(gamma2_expected, get(gamma2));
+}
+
 void test_geodesic_acceleration_compute() {
   static constexpr size_t Dim = 3;
   MAKE_GENERATOR(gen);
@@ -810,5 +844,6 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CurvedScalarWave.Worldtube.Tags",
   test_face_quantities_compute();
   test_puncture_field();
   test_check_input_file();
+  test_constraint_gammas_compute();
 }
 }  // namespace CurvedScalarWave::Worldtube