Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cmath>
       8             : #include <limits>
       9             : #include <tuple>
      10             : 
      11             : #include "DataStructures/DataBox/DataBox.hpp"
      12             : #include "IO/Observer/ObserverComponent.hpp"
      13             : #include "IO/Observer/ReductionActions.hpp"
      14             : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
      15             : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
      16             : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
      17             : #include "Parallel/GlobalCache.hpp"
      18             : #include "Parallel/Invoke.hpp"
      19             : #include "Parallel/ParallelComponentHelpers.hpp"
      20             : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
      21             : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
      22             : #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
      23             : #include "Utilities/Gsl.hpp"
      24             : #include "Utilities/PrettyType.hpp"
      25             : 
      26             : /// \cond
      27             : namespace Frame {
      28             : struct Grid;
      29             : struct Distorted;
      30             : struct Inertial;
      31             : }  // namespace Frame
      32             : /// \endcond
      33             : 
      34             : namespace ah {
      35             : namespace callbacks {
      36             : /*!
      37             :  * \brief Writes the center of an apparent horizon to disk in both the
      38             :  * `Frame` template parameter frame and Frame::Inertial frame. Intended to be
      39             :  * used in the `post_horizon_find_callbacks` list of an InterpolationTargetTag.
      40             :  *
      41             :  * The centers will be written to a subfile with the name
      42             :  * `/ApparentHorizons/TargetName_Centers.dat` where `TargetName` is the
      43             :  * pretty_type::name of the InterpolationTargetTag template parameter.
      44             :  *
      45             :  * The columns of the dat file are:
      46             :  * - %Time
      47             :  * - GridCenter_x
      48             :  * - GridCenter_y
      49             :  * - GridCenter_z
      50             :  * - InertialCenter_x
      51             :  * - InertialCenter_y
      52             :  * - InertialCenter_z
      53             :  *
      54             :  * The `Frame` template parameter must be either `::Frame::Grid` or
      55             :  * `::Frame::Distorted`. Even though the template parameter can be
      56             :  * `::Frame::Distorted`, we still write `GridCenter_?` because the centers of
      57             :  * the objects are the same in the Grid and Distorted frames.
      58             :  *
      59             :  * \note Requires ylm::Tags::Strahlkorper<Frame>
      60             :  * and ylm::Tags::CartesianCoords<Frame::Inertial> and
      61             :  * ylm::Tags::EuclideanAreaElement<Frame> to be in the DataBox
      62             :  * of the InterpolationTarget.
      63             :  */
      64             : template <typename InterpolationTargetTag, typename Frame>
      65           1 : struct ObserveCenters {
      66             :   // Note that we don't add a const_global_cache_tags type alias here with
      67             :   // ah::Tags::ObserveCenters because we want to use the base tag and so must be
      68             :   // agnostic to how the tag was added to the cache. We do this so anything that
      69             :   // uses ObserveCenters can control when it gets printed.
      70             : 
      71             :   template <typename DbTags, typename Metavariables, typename TemporalId>
      72           0 :   static void apply(const db::DataBox<DbTags>& box,
      73             :                     Parallel::GlobalCache<Metavariables>& cache,
      74             :                     const TemporalId& temporal_id) {
      75             :     static_assert(std::is_same_v<Frame, ::Frame::Grid> or
      76             :                       std::is_same_v<Frame, ::Frame::Distorted>,
      77             :                   "Frame must be either Grid or Distorted.");
      78             :     using HorizonTag = ylm::Tags::Strahlkorper<Frame>;
      79             :     using CoordsTag = ylm::Tags::CartesianCoords<::Frame::Inertial>;
      80             :     using EuclideanAreaElementTag = ylm::Tags::EuclideanAreaElement<Frame>;
      81             :     static_assert(db::tag_is_retrievable_v<HorizonTag, db::DataBox<DbTags>>,
      82             :                   "DataBox must contain ylm::Tags::Strahlkorper<Frame>");
      83             :     static_assert(db::tag_is_retrievable_v<CoordsTag, db::DataBox<DbTags>>,
      84             :                   "DataBox must contain "
      85             :                   "ylm::Tags::CartesianCoords<Frame::Inertial>");
      86             :     static_assert(
      87             :         db::tag_is_retrievable_v<EuclideanAreaElementTag, db::DataBox<DbTags>>,
      88             :         "DataBox must contain ylm::Tags::EuclideanAreaElement<Frame>");
      89             : 
      90             :     // Only print the centers if we want to.
      91             :     if (not Parallel::get<ah::Tags::ObserveCentersBase>(cache)) {
      92             :       return;
      93             :     }
      94             : 
      95             :     const auto& grid_horizon = db::get<HorizonTag>(box);
      96             :     const std::array<double, 3> grid_center = grid_horizon.physical_center();
      97             : 
      98             :     // computes the inertial center to be the average value of the
      99             :     // inertial coordinates over the Strahlkorper, where the average is
     100             :     // computed by a surface integral.
     101             :     // Note that we use the Euclidean area element here, since we are trying
     102             :     // to find a geometric center of a surface without regard to GR.
     103             :     const auto& inertial_coords = db::get<CoordsTag>(box);
     104             :     const auto& area_element = db::get<EuclideanAreaElementTag>(box);
     105             :     std::array<double, 3> inertial_center =
     106             :         make_array<3>(std::numeric_limits<double>::signaling_NaN());
     107             :     auto integrand = make_with_value<Scalar<DataVector>>(get(area_element), 0.);
     108             :     const double denominator = grid_horizon.ylm_spherepack().definite_integral(
     109             :         get(area_element).data());
     110             :     for (size_t i = 0; i < 3; ++i) {
     111             :       get(integrand) = get(area_element) * inertial_coords.get(i);
     112             :       gsl::at(inertial_center, i) =
     113             :           grid_horizon.ylm_spherepack().definite_integral(
     114             :               get(integrand).data()) /
     115             :           denominator;
     116             :     }
     117             : 
     118             :     // time, grid_x, grid_y, grid_z, inertial_x, inertial_y, inertial_z
     119             :     const auto center_tuple = std::make_tuple(
     120             :         intrp::InterpolationTarget_detail::get_temporal_id_value(temporal_id),
     121             :         grid_center[0], grid_center[1], grid_center[2], inertial_center[0],
     122             :         inertial_center[1], inertial_center[2]);
     123             : 
     124             :     auto& observer_writer_proxy = Parallel::get_parallel_component<
     125             :         observers::ObserverWriter<Metavariables>>(cache);
     126             : 
     127             :     const std::string subfile_path{"/ApparentHorizons/" +
     128             :                                    pretty_type::name<InterpolationTargetTag>() +
     129             :                                    "_Centers"};
     130             : 
     131             :     Parallel::threaded_action<
     132             :         observers::ThreadedActions::WriteReductionDataRow>(
     133             :         // Node 0 is always the writer
     134             :         observer_writer_proxy[0], subfile_path, legend_, center_tuple);
     135             :   }
     136             : 
     137             :  private:
     138           0 :   const static inline std::vector<std::string> legend_{
     139             :       {"Time", "GridCenter_x", "GridCenter_y", "GridCenter_z",
     140             :        "InertialCenter_x", "InertialCenter_y", "InertialCenter_z"}};
     141             : };
     142             : }  // namespace callbacks
     143             : }  // namespace ah