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 <deque>
       8             : #include <optional>
       9             : #include <string>
      10             : #include <type_traits>
      11             : #include <utility>
      12             : 
      13             : #include "DataStructures/DataBox/DataBox.hpp"
      14             : #include "DataStructures/DataVector.hpp"
      15             : #include "DataStructures/LinkedMessageId.hpp"
      16             : #include "DataStructures/Tensor/IndexType.hpp"
      17             : #include "DataStructures/Variables.hpp"
      18             : #include "DataStructures/VariablesTag.hpp"
      19             : #include "Domain/Creators/Tags/Domain.hpp"
      20             : #include "Domain/FunctionsOfTime/Tags.hpp"
      21             : #include "Domain/StrahlkorperTransformations.hpp"
      22             : #include "IO/Logging/Verbosity.hpp"
      23             : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
      24             : #include "NumericalAlgorithms/SphericalHarmonics/StrahlkorperFunctions.hpp"
      25             : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
      26             : #include "Parallel/GlobalCache.hpp"
      27             : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
      28             : #include "ParallelAlgorithms/ApparentHorizonFinder/HorizonAliases.hpp"
      29             : #include "ParallelAlgorithms/ApparentHorizonFinder/Storage.hpp"
      30             : #include "ParallelAlgorithms/ApparentHorizonFinder/Tags.hpp"
      31             : #include "Utilities/Gsl.hpp"
      32             : #include "Utilities/TMPL.hpp"
      33             : 
      34             : namespace ah {
      35             : 
      36             : namespace detail {
      37             : // Update the previous surface in the mutable global cache
      38             : template <typename Fr>
      39             : struct UpdatePreviousSurface {
      40             :   static void apply(const gsl::not_null<Storage::LockedPreviousSurface<Fr>*>
      41             :                         locked_previous_surface,
      42             :                     Storage::PreviousSurface<Fr> previous_surface) {
      43             :     auto& lock = locked_previous_surface->lock;
      44             :     lock.write_lock();
      45             :     locked_previous_surface->surface = std::move(previous_surface);
      46             :     lock.write_unlock();
      47             :   }
      48             : };
      49             : }  // namespace detail
      50             : 
      51             : /*!
      52             :  * \brief Invoke the callbacks specified in the `horizon_find_callbacks` alias
      53             :  * of the \p HorizonMetavars.
      54             :  *
      55             :  * \details Before invoking the callbacks, this function
      56             :  *
      57             :  * 1. Restricts the final interpolated variables from the $L_\mathrm{mesh}$ used
      58             :  *    for the FastFlow algorithm, to the actual $L$ of the Strahlkorper.
      59             :  * 2. Adds the current Strahlkorper to the `ah::Tags::PreviousSurfaces`.
      60             :  * 3. Copies the Strahlkorper, its time derivative, and the \p dependency into
      61             :  *    the box. Also possibly computes the Inertial coordinates of the final
      62             :  *    Strahlkorper and stores them in the box if the `frame` of the
      63             :  *    \p HorizonMetavars isn't the Inertial frame.
      64             :  */
      65             : template <typename HorizonMetavars, typename DbTags, typename Metavariables>
      66           1 : void invoke_callbacks(const gsl::not_null<db::DataBox<DbTags>*> box,
      67             :                       Parallel::GlobalCache<Metavariables>& cache,
      68             :                       const std::optional<std::string>& dependency,
      69             :                       const FastFlow::Status status) {
      70             :   using Fr = typename HorizonMetavars::frame;
      71             : 
      72             :   const auto& current_time = db::get<ah::Tags::CurrentTime>(*box).value();
      73             :   const auto& all_storage = db::get<ah::Tags::Storage<Fr>>(*box);
      74             :   const auto& current_time_storage = all_storage.at(current_time);
      75             :   const auto& current_iteration = current_time_storage.current_iteration;
      76             :   const auto& current_strahlkorper = current_iteration.strahlkorper;
      77             :   auto& previous_surfaces =
      78             :       db::get_mutable_reference<ah::Tags::PreviousSurfaces<Fr>>(box);
      79             :   const auto& fast_flow = db::get<ah::Tags::FastFlow>(*box);
      80             : 
      81             :   // Update the previous surfaces. We do this before the callbacks in case any
      82             :   // of the callbacks need the previous strahlkorpers with the current
      83             :   // strahlkorper already in it.
      84             :   previous_surfaces.emplace_front(current_time, current_strahlkorper,
      85             :                                   current_iteration.intersecting_element_ids);
      86             : 
      87             :   // Broadcast the previous surface to the mutable global cache so that elements
      88             :   // know whether they need to send data for the next horizon find. We do this
      89             :   // early so that elements get this information as soon as possible.
      90             :   Parallel::mutate<Tags::PreviousSurface<HorizonMetavars>,
      91             :                    detail::UpdatePreviousSurface<Fr>>(
      92             :       cache, previous_surfaces.front());
      93             : 
      94             :   // Remove old previous_strahlkorpers that are no longer relevant.
      95             :   const size_t num_previous_strahlkorpers = 3;
      96             :   while (previous_surfaces.size() > num_previous_strahlkorpers) {
      97             :     previous_surfaces.pop_back();
      98             :   }
      99             : 
     100             :   // The interpolated variables have been interpolated from the volume to
     101             :   // the points on the prolonged_strahlkorper, not to the points on the
     102             :   // actual strahlkorper. So here we do a restriction of these quantities
     103             :   // onto the actual strahlkorper.
     104             :   const auto& current_ylm = current_strahlkorper.ylm_spherepack();
     105             :   const size_t L_mesh = fast_flow.current_l_mesh(current_strahlkorper);
     106             :   const auto prolonged_strahlkorper =
     107             :       ylm::Strahlkorper<Fr>(L_mesh, L_mesh, current_strahlkorper);
     108             :   const auto& prolonged_ylm = prolonged_strahlkorper.ylm_spherepack();
     109             :   const auto& prolonged_interpolated_vars = current_iteration.interpolated_vars;
     110             :   db::mutate<::Tags::Variables<ah::vars_to_interpolate_to_target<3, Fr>>>(
     111             :       [&](const gsl::not_null<
     112             :           ::Variables<ah::vars_to_interpolate_to_target<3, Fr>>*>
     113             :               new_interpolated_vars) {
     114             :         new_interpolated_vars->initialize(current_ylm.physical_size());
     115             :         tmpl::for_each<ah::vars_to_interpolate_to_target<3, Fr>>(
     116             :             [&]<typename Tag>(tmpl::type_<Tag>) {
     117             :               const auto& prolonged_var = get<Tag>(prolonged_interpolated_vars);
     118             :               auto& new_var = get<Tag>(*new_interpolated_vars);
     119             :               for (size_t i = 0; i < prolonged_var.size(); i++) {
     120             :                 new_var[i] =
     121             :                     current_ylm.spec_to_phys(prolonged_ylm.prolong_or_restrict(
     122             :                         prolonged_ylm.phys_to_spec(prolonged_var[i]),
     123             :                         current_ylm));
     124             :               }
     125             :             });
     126             :       },
     127             :       box);
     128             : 
     129             :   db::mutate<ylm::Tags::Strahlkorper<Fr>, ylm::Tags::TimeDerivStrahlkorper<Fr>,
     130             :              ah::Tags::Dependency>(
     131             :       [&](const gsl::not_null<ylm::Strahlkorper<Fr>*> horizon,
     132             :           const gsl::not_null<ylm::Strahlkorper<Fr>*> time_deriv_of_horizon,
     133             :           const gsl::not_null<std::optional<std::string>*>
     134             :               callback_dependency) {
     135             :         *horizon = current_strahlkorper;
     136             : 
     137             :         // This is a hack to use ylm::time_deriv_of_strahlkorper function below
     138             :         std::deque<std::pair<double, ylm::Strahlkorper<Fr>>>
     139             :             previous_horizons{};
     140             :         for (const auto& previous_surface : previous_surfaces) {
     141             :           previous_horizons.emplace_back(previous_surface.time.id,
     142             :                                          previous_surface.surface);
     143             :         }
     144             : 
     145             :         *time_deriv_of_horizon = current_strahlkorper;
     146             :         ylm::time_deriv_of_strahlkorper(time_deriv_of_horizon,
     147             :                                         previous_horizons);
     148             : 
     149             :         *callback_dependency = dependency;
     150             :       },
     151             :       box);
     152             : 
     153             :   // Put inertial coords in the box if they aren't already there
     154             :   if constexpr (not std::is_same_v<Fr, Frame::Inertial>) {
     155             :     db::mutate<ylm::Tags::CartesianCoords<Frame::Inertial>>(
     156             :         [&](const gsl::not_null<tnsr::I<DataVector, 3>*> inertial_coords) {
     157             :           const auto& domain = Parallel::get<domain::Tags::Domain<3>>(cache);
     158             :           const auto& functions_of_time =
     159             :               Parallel::get<domain::Tags::FunctionsOfTime>(cache);
     160             :           strahlkorper_coords_in_different_frame(
     161             :               inertial_coords, current_strahlkorper, domain, functions_of_time,
     162             :               current_time.id);
     163             :         },
     164             :         box);
     165             :   }
     166             : 
     167             :   // Finally call callbacks
     168             :   tmpl::for_each<typename HorizonMetavars::horizon_find_callbacks>(
     169             :       [&]<typename Callback>(tmpl::type_<Callback>) {
     170             :         Callback::apply(*box, cache, status);
     171             :       });
     172             : }
     173             : }  // namespace ah