Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <cmath>
       8             : #include <cstddef>
       9             : #include <pup.h>
      10             : #include <string>
      11             : 
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "Evolution/Systems/GeneralizedHarmonic/Bbh/CompletionCriteria.hpp"
      14             : #include "Evolution/Systems/GeneralizedHarmonic/Bbh/CompletionSingleton.hpp"
      15             : #include "Evolution/Systems/GeneralizedHarmonic/Constraints.hpp"
      16             : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
      17             : #include "Options/String.hpp"
      18             : #include "Parallel/ArrayCollection/IsDgElementCollection.hpp"
      19             : #include "Parallel/GlobalCache.hpp"
      20             : #include "Parallel/Reduction.hpp"
      21             : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
      22             : #include "Time/Tags/Time.hpp"
      23             : #include "Utilities/ErrorHandling/Error.hpp"
      24             : #include "Utilities/Serialization/CharmPupable.hpp"
      25             : #include "Utilities/TMPL.hpp"
      26             : 
      27           0 : namespace gh::bbh::Events {
      28             : /*!
      29             :  * \brief Computes per-element Linf norms of the generalized-harmonic
      30             :  * constraints, latching BBH completion criteria when thresholds are exceeded.
      31             :  *
      32             :  * \details This event is intended to run once common-horizon finding is
      33             :  * active (typically behind a separation-based trigger). It monitors the gauge
      34             :  * and three-index constraints, using a reduction to determine if their Linf
      35             :  * norms exceed completion thresholds.
      36             :  */
      37           1 : class CheckConstraintThresholds : public Event {
      38           0 :   using ReductionData = Parallel::ReductionData<
      39             :       Parallel::ReductionDatum<double, funcl::AssertEqual<>>,
      40             :       Parallel::ReductionDatum<double, funcl::Max<>>,
      41             :       Parallel::ReductionDatum<double, funcl::Max<>>>;
      42             : 
      43             :  public:
      44             :   /// \cond
      45             :   explicit CheckConstraintThresholds(CkMigrateMessage* /*unused*/) {}
      46             :   using PUP::able::register_constructor;
      47             :   WRAPPED_PUPable_decl_template(CheckConstraintThresholds);  // NOLINT
      48             :   /// \endcond
      49             : 
      50           0 :   using compute_tags_for_observation_box =
      51             :       tmpl::list<gh::Tags::GaugeConstraintCompute<3, Frame::Inertial>,
      52             :                  gh::Tags::ThreeIndexConstraintCompute<3, Frame::Inertial>>;
      53           0 :   using options = tmpl::list<>;
      54           0 :   static constexpr Options::String help =
      55             :       "Checks local Linf norms of constraints against BBH completion "
      56             :       "thresholds and forwards reduced maxima to the BBH completion singleton "
      57             :       "reduction callback.";
      58           0 :   static std::string name() { return "BbhCheckConstraintThresholds"; }
      59             : 
      60           0 :   CheckConstraintThresholds() = default;
      61             : 
      62           0 :   using return_tags = tmpl::list<>;
      63           0 :   using argument_tags = tmpl::list<
      64             :       ::Tags::Time, gh::Tags::GaugeConstraint<DataVector, 3, Frame::Inertial>,
      65             :       gh::Tags::ThreeIndexConstraint<DataVector, 3, Frame::Inertial>>;
      66             : 
      67             :   template <typename Metavariables, typename ArrayIndex, typename Component>
      68           0 :   void operator()(
      69             :       const double time,
      70             :       const tnsr::a<DataVector, 3, Frame::Inertial>& gauge_constraint,
      71             :       const tnsr::iaa<DataVector, 3, Frame::Inertial>& three_index_constraint,
      72             :       Parallel::GlobalCache<Metavariables>& cache,
      73             :       const ArrayIndex& array_index, const Component* const /*component*/,
      74             :       const ObservationValue& /*observation_value*/) const {
      75             :     const double local_gauge_linf = local_linf_norm(gauge_constraint);
      76             :     const double local_three_index_linf =
      77             :         local_linf_norm(three_index_constraint);
      78             :     if constexpr (Parallel::is_dg_element_collection_v<Component>) {
      79             :       ERROR(
      80             :           "BbhCheckConstraintThresholds currently requires array components "
      81             :           "(not DgElementCollection).");
      82             :     } else {
      83             :       const auto& self_proxy =
      84             :           Parallel::get_parallel_component<Component>(cache)[array_index];
      85             :       auto& reduction_target_proxy = Parallel::get_parallel_component<
      86             :           gh::bbh::CompletionSingleton<Metavariables>>(cache);
      87             :       Parallel::contribute_to_reduction<
      88             :           gh::bbh::Actions::ProcessConstraintMaxima>(
      89             :           ReductionData{time, local_gauge_linf, local_three_index_linf},
      90             :           self_proxy, reduction_target_proxy);
      91             :     }
      92             :   }
      93             : 
      94           0 :   using is_ready_argument_tags = tmpl::list<>;
      95             :   template <typename Metavariables, typename ArrayIndex, typename Component>
      96           0 :   bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
      97             :                 const ArrayIndex& /*array_index*/,
      98             :                 const Component* const /*meta*/) const {
      99             :     return true;
     100             :   }
     101             : 
     102           1 :   bool needs_evolved_variables() const override { return true; }
     103             : 
     104             :  private:
     105             :   template <typename TensorType>
     106           0 :   static double local_linf_norm(const TensorType& tensor) {
     107             :     double result = 0.0;
     108             :     for (size_t storage_index = 0; storage_index < tensor.size();
     109             :          ++storage_index) {
     110             :       const auto& component = tensor[storage_index];
     111             :       result = std::max(result, max(abs(component)));
     112             :     }
     113             :     return result;
     114             :   }
     115             : };
     116             : }  // namespace gh::bbh::Events