Merge pull request #3921 from knelli2/obs_ah_centers

Add ObserveCenters post horizon find callback

src/ApparentHorizons/CMakeLists.txt

@@ -23,6 +23,7 @@ spectre_target_headers(
   ComputeItems.hpp
   FastFlow.hpp
   ObjectLabel.hpp
+  ObserveCenters.hpp
   StrahlkorperGr.hpp
   StrahlkorperInDifferentFrame.hpp
   Tags.hpp

src/ApparentHorizons/ObserveCenters.hpp

@@ -0,0 +1,93 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <tuple>
+
+#include "DataStructures/DataBox/DataBox.hpp"
+#include "IO/Observer/ObserverComponent.hpp"
+#include "IO/Observer/ReductionActions.hpp"
+#include "Parallel/GlobalCache.hpp"
+#include "Parallel/Invoke.hpp"
+#include "Parallel/ParallelComponentHelpers.hpp"
+#include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
+#include "Utilities/PrettyType.hpp"
+
+/// \cond
+namespace StrahlkorperTags {
+template <typename Frame>
+struct Strahlkorper;
+}  // namespace StrahlkorperTags
+namespace Frame {
+struct Grid;
+struct Inertial;
+}  // namespace Frame
+/// \endcond
+
+namespace ah {
+namespace callbacks {
+/*!
+ * \brief Writes the center of an apparent horizon to disk in both the
+ * Frame::Grid frame and Frame::Inertial frame. Intended to be used in the
+ * `post_horizon_find_callbacks` list of an InterpolationTargetTag.
+ *
+ * The centers will be written to a subfile with the name
+ * `/ApparentHorizons/TargetName_Centers.dat` where `TargetName` is the
+ * pretty_type::name of the InterpolationTargetTag template parameter.
+ *
+ * The columns of the dat file are:
+ * - %Time
+ * - GridCenter_x
+ * - GridCenter_y
+ * - GridCenter_z
+ * - InertialCenter_x
+ * - InertialCenter_y
+ * - InertialCenter_z
+ *
+ * \note Requires StrahlkorperTags::Strahlkorper<Frame::Grid> and
+ * StrahlkorperTags::Strahlkorper<Frame::Inertial> be in the DataBox of the
+ * InterpolationTarget.
+ */
+template <typename InterpolationTargetTag>
+struct ObserveCenters {
+  template <typename DbTags, typename Metavariables, typename TemporalId>
+  static void apply(const db::DataBox<DbTags>& box,
+                    Parallel::GlobalCache<Metavariables>& cache,
+                    const TemporalId& temporal_id) {
+    using GridTag = StrahlkorperTags::Strahlkorper<Frame::Grid>;
+    using InertialTag = StrahlkorperTags::Strahlkorper<Frame::Inertial>;
+
+    const auto& grid_horizon = db::get<GridTag>(box);
+    const auto& inertial_horizon = db::get<InertialTag>(box);
+
+    const std::array<double, 3> grid_center = grid_horizon.physical_center();
+    const std::array<double, 3> inertial_center =
+        inertial_horizon.physical_center();
+
+    // time, grid_x, grid_y, grid_z, inertial_x, inertial_y, inertial_z
+    const auto center_tuple = std::make_tuple(
+        intrp::InterpolationTarget_detail::get_temporal_id_value(temporal_id),
+        grid_center[0], grid_center[1], grid_center[2], inertial_center[0],
+        inertial_center[1], inertial_center[2]);
+
+    auto& observer_writer_proxy = Parallel::get_parallel_component<
+        observers::ObserverWriter<Metavariables>>(cache);
+
+    const std::string subfile_path{"/ApparentHorizons/" +
+                                   pretty_type::name<InterpolationTargetTag>() +
+                                   "_Centers"};
+
+    Parallel::threaded_action<
+        observers::ThreadedActions::WriteReductionDataRow>(
+        // Node 0 is always the writer
+        observer_writer_proxy[0], subfile_path, legend_, center_tuple);
+  }
+
+ private:
+  const static inline std::vector<std::string> legend_{
+      {"Time", "GridCenter_x", "GridCenter_y", "GridCenter_z",
+       "InertialCenter_x", "InertialCenter_y", "InertialCenter_z"}};
+};
+}  // namespace callbacks
+}  // namespace ah

src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp

@@ -9,6 +9,7 @@
 #include "ApparentHorizons/ComputeHorizonVolumeQuantities.hpp"
 #include "ApparentHorizons/ComputeHorizonVolumeQuantities.tpp"
 #include "ApparentHorizons/ComputeItems.hpp"
+#include "ApparentHorizons/ObserveCenters.hpp"
 #include "ApparentHorizons/Tags.hpp"
 #include "DataStructures/DataBox/PrefixHelpers.hpp"
 #include "DataStructures/DataBox/Tag.hpp"
@@ -269,7 +270,8 @@ struct EvolutionMetavars {
         intrp::callbacks::ErrorOnFailedApparentHorizon;
     using post_horizon_find_callbacks = tmpl::list<
         intrp::callbacks::ObserveTimeSeriesOnSurface<tags_to_observe, AhA>,
-        intrp::callbacks::ObserveSurfaceData<surface_tags_to_observe, AhA>>;
+        intrp::callbacks::ObserveSurfaceData<surface_tags_to_observe, AhA>,
+        ah::callbacks::ObserveCenters<AhA>>;
   };
 
   struct AhB : tt::ConformsTo<intrp::protocols::InterpolationTargetTag> {
@@ -288,7 +290,8 @@ struct EvolutionMetavars {
         intrp::callbacks::ErrorOnFailedApparentHorizon;
     using post_horizon_find_callbacks = tmpl::list<
         intrp::callbacks::ObserveTimeSeriesOnSurface<tags_to_observe, AhB>,
-        intrp::callbacks::ObserveSurfaceData<surface_tags_to_observe, AhB>>;
+        intrp::callbacks::ObserveSurfaceData<surface_tags_to_observe, AhB>,
+        ah::callbacks::ObserveCenters<AhB>>;
   };
 
   using interpolation_target_tags = tmpl::list<AhA, AhB>;

tests/Unit/ApparentHorizons/CMakeLists.txt

@@ -9,6 +9,7 @@ set(LIBRARY_SOURCES
   Test_ComputeItems.cpp
   Test_FastFlow.cpp
   Test_ObjectLabel.cpp
+  Test_ObserveCenters.cpp
   Test_StrahlkorperGr.cpp
   Test_StrahlkorperInDifferentFrame.cpp
   Test_Tags.cpp
@@ -29,6 +30,7 @@ target_link_libraries(
   GeneralRelativitySolutions
   LinearOperators
   Logging
+  ObserverHelpers
   Options
   ParallelInterpolation
   RootFinding

tests/Unit/ApparentHorizons/Test_ObserveCenters.cpp

@@ -0,0 +1,144 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <array>
+#include <cstddef>
+#include <random>
+#include <string>
+
+#include "ApparentHorizons/ObserveCenters.hpp"
+#include "DataStructures/DataBox/DataBox.hpp"
+#include "DataStructures/Matrix.hpp"
+#include "Framework/ActionTesting.hpp"
+#include "Helpers/DataStructures/MakeWithRandomValues.hpp"
+#include "Helpers/IO/Observers/MockWriteReductionDataRow.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
+#include "Parallel/GlobalCache.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/PrettyType.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace {
+
+// This is needed to template ObserveCenters, but it doesn't have to be an
+// actual InterpolationTargetTag. This is only used for the name of the subfile
+// written
+struct AhA {};
+
+struct TestMetavars {
+  enum class Phase { Initialization, WriteData, Exit };
+
+  using component_list =
+      tmpl::list<TestHelpers::observers::MockObserverWriter<TestMetavars>>;
+};
+
+using FoTPtr = std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>;
+
+SPECTRE_TEST_CASE("Unit.ApparentHorizons.ObserveCenters",
+                  "[Unit][ApparentHorizons]") {
+  using observer_writer =
+      TestHelpers::observers::MockObserverWriter<TestMetavars>;
+  MAKE_GENERATOR(gen);
+  std::uniform_real_distribution<double> center_dist{-10.0, 10.0};
+
+  // set up runner and stuff
+  ActionTesting::MockRuntimeSystem<TestMetavars> runner{{}};
+  runner.set_phase(TestMetavars::Phase::Initialization);
+  ActionTesting::emplace_nodegroup_component_and_initialize<observer_writer>(
+      make_not_null(&runner), {});
+  auto& cache = ActionTesting::cache<observer_writer>(runner, 0);
+
+  runner.set_phase(TestMetavars::Phase::WriteData);
+
+  const auto make_center = [&gen, &center_dist]() -> std::array<double, 3> {
+    return make_with_random_values<std::array<double, 3>>(
+        make_not_null(&gen), center_dist, std::array<double, 3>{});
+  };
+
+  // Lists of centers so we can check that the correct centers were written
+  std::vector<std::array<double, 3>> grid_centers{};
+  std::vector<std::array<double, 3>> inertial_centers{};
+
+  db::DataBox<tmpl::list<StrahlkorperTags::Strahlkorper<Frame::Grid>,
+                         StrahlkorperTags::Strahlkorper<Frame::Inertial>>>
+      box{};
+
+  const auto update_stored_centers = [&make_center, &grid_centers,
+                                      &inertial_centers, &box]() {
+    const auto grid_center = make_center();
+    const auto inertial_center = make_center();
+
+    grid_centers.push_back(grid_center);
+    inertial_centers.push_back(inertial_center);
+
+    db::mutate<StrahlkorperTags::Strahlkorper<Frame::Grid>,
+               StrahlkorperTags::Strahlkorper<Frame::Inertial>>(
+        make_not_null(&box),
+        [&grid_center, &inertial_center](
+            gsl::not_null<Strahlkorper<Frame::Grid>*> box_grid_horizon,
+            gsl::not_null<Strahlkorper<Frame::Inertial>*>
+                box_inertial_horizon) {
+          // We only care about the centers of the strahlkorper so the
+          // ell and radius are arbitrary
+          *box_grid_horizon = Strahlkorper<Frame::Grid>{2, 1.0, grid_center};
+          *box_inertial_horizon =
+              Strahlkorper<Frame::Inertial>{2, 1.0, inertial_center};
+        });
+  };
+
+  // times to write
+  const std::vector<double> times{0.0, 0.1, 0.2, 0.3, 0.4, 0.5};
+
+  // write some data
+  for (size_t i = 0; i < times.size(); i++) {
+    update_stored_centers();
+
+    ah::callbacks::ObserveCenters<AhA>::apply(box, cache, times[i]);
+
+    size_t num_threaded_actions =
+        ActionTesting::number_of_queued_threaded_actions<observer_writer>(
+            runner, 0);
+    CHECK(num_threaded_actions == 1);
+    ActionTesting::invoke_queued_threaded_action<observer_writer>(
+        make_not_null(&runner), 0);
+  }
+
+  // These must be the same as in ObserveCenters
+  const std::vector<std::string> compare_legend{
+      {"Time", "GridCenter_x", "GridCenter_y", "GridCenter_z",
+       "InertialCenter_x", "InertialCenter_y", "InertialCenter_z"}};
+  const std::string subfile_name =
+      "/ApparentHorizons/" + pretty_type::name<AhA>() + "_Centers";
+
+  auto& h5_file = ActionTesting::get_databox_tag<
+      observer_writer, TestHelpers::observers::MockReductionFileTag>(runner, 0);
+  const auto& dataset = h5_file.get_dat(subfile_name);
+  const Matrix data = dataset.get_data();
+  const std::vector<std::string>& legend = dataset.get_legend();
+
+  // Check legend is correct
+  for (size_t i = 0; i < legend.size(); i++) {
+    CHECK(legend[i] == compare_legend[i]);
+  }
+
+  // Check proper number of times were written
+  CHECK(data.rows() == times.size());
+
+  // Check centers
+  for (size_t i = 0; i < times.size(); i++) {
+    CHECK(data(i, 0) == times[i]);
+
+    const std::array<double, 3>& grid_center = grid_centers[i];
+    const std::array<double, 3>& inertial_center = inertial_centers[i];
+    for (size_t j = 0; j < grid_center.size(); j++) {
+      // Grid center is columns 2-4
+      CHECK(data(i, j + 1) == gsl::at(grid_center, j));
+      // Inertial center is columns 5-7
+      CHECK(data(i, j + 4) == gsl::at(inertial_center, j));
+    }
+  }
+}
+}  // namespace