Add InsideHorizon trigger for worldtube excision

src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt

@@ -48,3 +48,5 @@ target_link_libraries(
 
 add_subdirectory(SingletonActions)
 add_subdirectory(ElementActions)
+add_subdirectory(Triggers)
+

src/Evolution/Systems/CurvedScalarWave/Worldtube/Triggers/CMakeLists.txt

@@ -0,0 +1,16 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+
+spectre_target_sources(
+  ${LIBRARY}
+  PRIVATE
+  InsideHorizon.cpp
+  )
+
+spectre_target_headers(
+  ${LIBRARY}
+  INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
+  HEADERS
+  InsideHorizon.hpp
+)

src/Evolution/Systems/CurvedScalarWave/Worldtube/Triggers/InsideHorizon.cpp

@@ -0,0 +1,31 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/Triggers/InsideHorizon.hpp"
+
+#include <array>
+
+#include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/RadiusFunctions.hpp"
+
+namespace Triggers {
+
+bool InsideHorizon::operator()(
+    const std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>&
+        position_and_velocity,
+    const std::array<double, 4>& worldtube_radius_params) const {
+  const double orbit_radius = get(magnitude(position_and_velocity[0]));
+  // optimization that will almost always be true
+  if (LIKELY(orbit_radius > 1.99)) {
+    return false;
+  }
+  const double wt_radius_inertial =
+      CurvedScalarWave::Worldtube::smooth_broken_power_law(
+          orbit_radius, worldtube_radius_params[0], worldtube_radius_params[1],
+          worldtube_radius_params[2], worldtube_radius_params[3]);
+  return orbit_radius + wt_radius_inertial < 1.99;
+}
+
+PUP::able::PUP_ID InsideHorizon::my_PUP_ID = 0;  // NOLINT
+}  // namespace Triggers

src/Evolution/Systems/CurvedScalarWave/Worldtube/Triggers/InsideHorizon.hpp

@@ -0,0 +1,45 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <limits>
+#include <pup.h>
+
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
+#include "Utilities/Serialization/CharmPupable.hpp"
+
+namespace Triggers {
+
+/*!
+ * \brief This trigger is true when the worldtube is entirely within a
+ * coordinate sphere of radius 1.99 M centered on the origin in the inertial
+ * frame. This assumes a black hole mass of 1M.
+ */
+class InsideHorizon : public Trigger {
+ public:
+  /// \cond
+  InsideHorizon() = default;
+  explicit InsideHorizon(CkMigrateMessage* /*unused*/) {}
+  using PUP::able::register_constructor;
+  WRAPPED_PUPable_decl_template(InsideHorizon);  // NOLINT
+  /// \endcond
+
+  static constexpr Options::String help =
+      "Triggers if the worldtube is entirely inside a coordinate sphere with "
+      "radius 1.99 M centered on the origin in the inertial frame.";
+  static constexpr size_t Dim = 3;
+  using options = tmpl::list<>;
+  using argument_tags = tmpl::list<
+      CurvedScalarWave::Worldtube::Tags::ParticlePositionVelocity<Dim>,
+      CurvedScalarWave::Worldtube::Tags::WorldtubeRadiusParameters>;
+
+  bool operator()(const std::array<tnsr::I<double, Dim, Frame::Inertial>, 2>&
+                      position_and_velocity,
+                  const std::array<double, 4>& worldtube_radius_params) const;
+
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& /*p*/) override {}
+};
+}  // namespace Triggers

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt

@@ -10,6 +10,7 @@ set(LIBRARY_SOURCES
   Test_RadiusFunctions.cpp
   Test_SelfForce.cpp
   Test_Tags.cpp
+  Test_Triggers.cpp
   ElementActions/Test_SendToWorldtube.cpp
   ElementActions/Test_ReceiveWorldtubeData.cpp
   ElementActions/Test_InitializeConstraintGammas.cpp

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_Triggers.cpp

@@ -0,0 +1,51 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <cmath>
+
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/RadiusFunctions.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/Triggers/InsideHorizon.hpp"
+#include "Framework/TestHelpers.hpp"
+
+namespace CurvedScalarWave::Worldtube {
+namespace {
+void test_inside_horizon() {
+  const double exp = 1.5;
+  const double rinf_large = 0.1;
+  const double rb = 0.1;
+  const double delta = 0.01;
+
+  const std::array<double, 4> wt_radius_params_large{
+      {exp, rinf_large, rb, delta}};
+  // unused
+  const tnsr::I<double, 3> vel{{0., 0., 0.}};
+
+  // particle outside the horizon
+  const tnsr::I<double, 3> pos_outside{{0., 2., 0.}};
+  const std::array<tnsr::I<double, 3>, 2> pos_vel_outside{{pos_outside, vel}};
+  const Triggers::InsideHorizon inside_horizon_trigger{};
+  CHECK(inside_horizon_trigger(pos_vel_outside, wt_radius_params_large) ==
+        false);
+
+  // particle inside horizon but worldtube sticks out
+  const tnsr::I<double, 3> pos_inside{{1.95, 0., 0.}};
+  const std::array<tnsr::I<double, 3>, 2> pos_vel_inside{{pos_inside, vel}};
+  CHECK(inside_horizon_trigger(pos_vel_inside, wt_radius_params_large) ==
+        false);
+
+  // shrink worldtube radius to fit inside horizons
+  const double rinf_small = 0.01;
+  const std::array<double, 4> wt_radius_params_small{
+      {exp, rinf_small, rb, delta}};
+
+  CHECK(inside_horizon_trigger(pos_vel_inside, wt_radius_params_small) == true);
+}
+
+SPECTRE_TEST_CASE("Unit.Evolution.Systems.CurvedScalarWave.Worldtube.Triggers",
+                  "[Unit][Evolution]") {
+  test_inside_horizon();
+}
+}  // namespace
+}  // namespace CurvedScalarWave::Worldtube