Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <optional>
       8             : #include <pup.h>
       9             : #include <string>
      10             : #include <unordered_map>
      11             : #include <utility>
      12             : #include <vector>
      13             : 
      14             : #include "DataStructures/DataBox/DataBox.hpp"
      15             : #include "DataStructures/DataBox/ObservationBox.hpp"
      16             : #include "DataStructures/DataBox/TagName.hpp"
      17             : #include "DataStructures/DataVector.hpp"
      18             : #include "DataStructures/Tensor/Tensor.hpp"
      19             : #include "Domain/BlockLogicalCoordinates.hpp"
      20             : #include "Domain/Creators/Tags/Domain.hpp"
      21             : #include "Domain/ElementLogicalCoordinates.hpp"
      22             : #include "Domain/Structure/ElementId.hpp"
      23             : #include "Domain/Tags.hpp"
      24             : #include "Elliptic/Systems/SelfForce/Scalar/AnalyticData/CircularOrbit.hpp"
      25             : #include "Elliptic/Systems/SelfForce/Scalar/Tags.hpp"
      26             : #include "Elliptic/Tags.hpp"
      27             : #include "IO/Observer/GetSectionObservationKey.hpp"
      28             : #include "IO/Observer/Helpers.hpp"
      29             : #include "IO/Observer/ObservationId.hpp"
      30             : #include "IO/Observer/ObserverComponent.hpp"
      31             : #include "IO/Observer/ReductionActions.hpp"
      32             : #include "IO/Observer/TypeOfObservation.hpp"
      33             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      34             : #include "Options/String.hpp"
      35             : #include "Parallel/ArrayIndex.hpp"
      36             : #include "Parallel/GlobalCache.hpp"
      37             : #include "Parallel/Invoke.hpp"
      38             : #include "Parallel/Local.hpp"
      39             : #include "Parallel/Reduction.hpp"
      40             : #include "Parallel/TypeTraits.hpp"
      41             : #include "ParallelAlgorithms/Events/Tags.hpp"
      42             : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
      43             : #include "Utilities/ErrorHandling/Assert.hpp"
      44             : #include "Utilities/ErrorHandling/Error.hpp"
      45             : #include "Utilities/Functional.hpp"
      46             : #include "Utilities/OptionalHelpers.hpp"
      47             : #include "Utilities/PrettyType.hpp"
      48             : #include "Utilities/Serialization/CharmPupable.hpp"
      49             : #include "Utilities/TMPL.hpp"
      50             : 
      51           0 : namespace ScalarSelfForce::Events {
      52             : 
      53             : namespace detail {
      54             : tnsr::i<std::complex<double>, 2> extract_self_force(
      55             :     const tnsr::I<double, 2, Frame::ElementLogical>& puncture_logical_coords,
      56             :     const AnalyticData::CircularOrbit& circular_orbit,
      57             :     const Scalar<ComplexDataVector>& field, const Mesh<2>& mesh,
      58             :     const InverseJacobian<DataVector, 2, Frame::ElementLogical,
      59             :                           Frame::Inertial>& inv_jacobian);
      60             : }  // namespace detail
      61             : 
      62             : /*!
      63             :  * \brief Observe the self-force at the position of the scalar charge.
      64             :  *
      65             :  * \details This event interpolates the scalar m-mode field and its derivatives
      66             :  * to the position of the scalar charge (puncture) and computes the m-mode
      67             :  * contribution to the self-force acting on the charge. The self-force is then
      68             :  * written to a reduction file.
      69             :  *
      70             :  * The self-force is computed in the Boyer-Lindquist `r` and `theta`
      71             :  * coordinates (for scalar charge $q=1$) following Eq. (3.10) in
      72             :  * \cite Osburn:2022bby :
      73             :  * \begin{align}
      74             :  *   F_r^m &= \partial_r (\Psi_m^R / r) = \frac{1}{r \alpha}
      75             :  *     \partial_{r_*} \Psi_m^R - \frac{1}{r^2} \Psi_m^R \\
      76             :  *   F_\theta^m &= \partial_\theta (\Psi_m^R / r) = -\frac{1}{r}
      77             :  *     \partial_{\cos\theta} \Psi_m^R
      78             :  * \end{align}
      79             :  * with $\alpha = 1 - 2 M r / (r^2 + a^2)$. For modes $m > 0$ the self-force
      80             :  * includes an additional factor of 2 and a complex rotation by
      81             :  * $m \Delta\phi = \frac{m a}{r_+ - r_-}\ln\frac{r - r_+}{r - r_-}$.
      82             :  *
      83             :  * These contributions can be summed over all m-modes (in postprocessing) to
      84             :  * form the total self-force $F_\mu^\mathrm{self} = \sum_m F_\mu^m$, where the
      85             :  * sum is truncated at some maximum m-mode number.
      86             :  */
      87             : template <typename ArraySectionIdTag = void>
      88           1 : class ObserveSelfForce : public Event {
      89             :  public:
      90           0 :   static constexpr Options::String help =
      91             :       "Observe the self-force at the position of the scalar charge.";
      92           0 :   using options = tmpl::list<>;
      93             : 
      94           0 :   ObserveSelfForce() = default;
      95             : 
      96           0 :   explicit ObserveSelfForce(CkMigrateMessage* msg) : Event(msg) {}
      97             :   using PUP::able::register_constructor;
      98           0 :   WRAPPED_PUPable_decl_template(ObserveSelfForce);  // NOLINT
      99             : 
     100           0 :   using compute_tags_for_observation_box = tmpl::list<>;
     101             : 
     102           0 :   using return_tags = tmpl::list<>;
     103           0 :   using argument_tags = tmpl::list<::Tags::ObservationBox>;
     104             : 
     105           0 :   using BackgroundTag =
     106             :       elliptic::Tags::Background<elliptic::analytic_data::Background>;
     107             : 
     108             :   template <typename ComputeTagsList, typename DataBoxType,
     109             :             typename Metavariables, typename ParallelComponent>
     110           0 :   void operator()(const ObservationBox<ComputeTagsList, DataBoxType>& box,
     111             :                   Parallel::GlobalCache<Metavariables>& cache,
     112             :                   const ElementId<2>& element_id,
     113             :                   const ParallelComponent* const /*meta*/,
     114             :                   const ObservationValue& observation_value) const {
     115             :     // Skip observation on elements that are not part of the section
     116             :     const std::optional<std::string> section_observation_key =
     117             :         observers::get_section_observation_key<ArraySectionIdTag>(box);
     118             :     if (not section_observation_key.has_value()) {
     119             :       return;
     120             :     }
     121             :     const auto& background = get<BackgroundTag>(box);
     122             :     const auto& circular_orbit =
     123             :         dynamic_cast<const AnalyticData::CircularOrbit&>(background);
     124             :     // Get element-logical coords of puncture
     125             :     const auto& domain = get<domain::Tags::Domain<2>>(box);
     126             :     const auto puncture_position = circular_orbit.puncture_position();
     127             :     const auto& block = domain.blocks()[element_id.block_id()];
     128             :     const auto block_logical_coords =
     129             :         block_logical_coordinates_single_point(puncture_position, block);
     130             :     if (not block_logical_coords.has_value()) {
     131             :       return;
     132             :     }
     133             :     const auto puncture_logical_coords =
     134             :         element_logical_coordinates(block_logical_coords.value(), element_id);
     135             :     if (not puncture_logical_coords.has_value()) {
     136             :       return;
     137             :     }
     138             :     // Extract self-force
     139             :     const auto& field = get<Tags::MMode>(box);
     140             :     const auto& mesh = get<domain::Tags::Mesh<2>>(box);
     141             :     const auto& inv_jacobian =
     142             :         get<domain::Tags::InverseJacobian<2, Frame::ElementLogical,
     143             :                                           Frame::Inertial>>(box);
     144             :     const auto self_force =
     145             :         detail::extract_self_force(puncture_logical_coords.value(),
     146             :                                    circular_orbit, field, mesh, inv_jacobian);
     147             :     // Write result to file
     148             :     auto& reduction_writer = Parallel::get_parallel_component<
     149             :         observers::ObserverWriter<Metavariables>>(cache);
     150             :     Parallel::threaded_action<
     151             :         observers::ThreadedActions::WriteReductionDataRow>(
     152             :         reduction_writer[0], std::string{"/SelfForce"},
     153             :         std::vector<std::string>{"IterationId", "NumberOfGridPoints",
     154             :                                  "Re(SelfForce_r)", "Im(SelfForce_r)",
     155             :                                  "Re(SelfForce_theta)", "Im(SelfForce_theta)"},
     156             :         std::make_tuple(observation_value.value, mesh.number_of_grid_points(),
     157             :                         get<0>(self_force).real(), get<0>(self_force).imag(),
     158             :                         get<1>(self_force).real(), get<1>(self_force).imag()));
     159             :   }
     160             : 
     161           0 :   using observation_registration_tags = tmpl::list<::Tags::DataBox>;
     162             : 
     163             :   template <typename DbTagsList>
     164             :   std::optional<
     165             :       std::pair<observers::TypeOfObservation, observers::ObservationKey>>
     166           0 :   get_observation_type_and_key_for_registration(
     167             :       const db::DataBox<DbTagsList>& box) const {
     168             :     const std::optional<std::string> section_observation_key =
     169             :         observers::get_section_observation_key<ArraySectionIdTag>(box);
     170             :     if (not section_observation_key.has_value()) {
     171             :       return std::nullopt;
     172             :     }
     173             :     return {
     174             :         {observers::TypeOfObservation::Reduction,
     175             :          observers::ObservationKey("ObserveSelfForce" +
     176             :                                    section_observation_key.value() + ".dat")}};
     177             :   }
     178             : 
     179           0 :   using is_ready_argument_tags = tmpl::list<>;
     180             : 
     181             :   template <typename Metavariables, typename ArrayIndex, typename Component>
     182           0 :   bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
     183             :                 const ArrayIndex& /*array_index*/,
     184             :                 const Component* const /*meta*/) const {
     185             :     return true;
     186             :   }
     187             : 
     188           1 :   bool needs_evolved_variables() const override { return false; }
     189             : };
     190             : 
     191             : /// \cond
     192             : // NOLINTBEGIN(cppcoreguidelines-avoid-non-const-global-variables)
     193             : template <typename ArraySectionIdTag>
     194             : PUP::able::PUP_ID ObserveSelfForce<ArraySectionIdTag>::my_PUP_ID = 0;
     195             : // NOLINTEND(cppcoreguidelines-avoid-non-const-global-variables)
     196             : /// \endcond
     197             : }  // namespace ScalarSelfForce::Events