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 <set>
      10             : #include <unordered_map>
      11             : #include <unordered_set>
      12             : #include <vector>
      13             : 
      14             : #include "DataStructures/DataVector.hpp"
      15             : #include "DataStructures/LinkedMessageId.hpp"
      16             : #include "DataStructures/Tensor/Tensor.hpp"
      17             : #include "DataStructures/Variables.hpp"
      18             : #include "Domain/BlockLogicalCoordinates.hpp"
      19             : #include "Domain/Structure/ElementId.hpp"
      20             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      21             : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
      22             : #include "Parallel/MultiReaderSpinlock.hpp"
      23             : #include "ParallelAlgorithms/ApparentHorizonFinder/Destination.hpp"
      24             : #include "ParallelAlgorithms/ApparentHorizonFinder/FastFlow.hpp"
      25             : #include "ParallelAlgorithms/ApparentHorizonFinder/HorizonAliases.hpp"
      26             : 
      27           0 : namespace ah::Storage {
      28             : /*!
      29             :  * \brief Holds the `ah::source_vars`, mesh, and other variables on the horizon
      30             :  * finder for a given element.
      31             :  */
      32             : template <typename Fr>
      33           1 : struct VolumeVariables {
      34             :   /*!
      35             :    * \brief The mesh that corresponds to the volume vars.
      36             :    */
      37           1 :   Mesh<3> mesh;
      38             : 
      39             :   /*!
      40             :    * \brief A `Variables` of the tensors in the volume that we need to
      41             :    * interpolate onto the horizon.
      42             :    */
      43             :   Variables<ah::vars_to_interpolate_to_target<3, Fr>>
      44           1 :       vars_to_interpolate_to_target{};
      45             : 
      46             :   // NOLINTNEXTLINE(google-runtime-references)
      47           0 :   void pup(PUP::er& p);
      48             : };
      49             : 
      50             : template <typename Fr>
      51           0 : bool operator==(const VolumeVariables<Fr>& lhs, const VolumeVariables<Fr>& rhs);
      52             : template <typename Fr>
      53           0 : bool operator!=(const VolumeVariables<Fr>& lhs, const VolumeVariables<Fr>& rhs);
      54             : 
      55             : /*!
      56             :  * \brief Holds the `ylm::Strahlkorper` and associated quantities for a single
      57             :  * `FastFlow` iteration.
      58             :  */
      59             : template <typename Fr>
      60           1 : struct Iteration {
      61             :   /*!
      62             :    * \brief The surface for this iteration
      63             :    */
      64           1 :   ylm::Strahlkorper<Fr> strahlkorper;
      65             :   /*!
      66             :    *  \brief Holds the list of all points (in block logical coordinates) that
      67             :    *  need to be interpolated onto.
      68             :    *
      69             :    * \details If an element of the vector is `std::nullopt`, then that point is
      70             :    * outside the domain.
      71             :    */
      72           1 :   std::optional<std::vector<BlockLogicalCoords<3>>> block_coord_holders;
      73             :   /*!
      74             :    * \brief Holds the interpolated `Variables` on the points in
      75             :    * `block_coord_holders`.
      76             :    *
      77             :    * \details The grid points inside are indexed according to
      78             :    * `block_coord_holders`.
      79             :    */
      80           1 :   Variables<ah::vars_to_interpolate_to_target<3, Fr>> interpolated_vars{};
      81             :   /*!
      82             :    * \brief Keeps track of the indices in `interpolated_vars` that have
      83             :    * already been interpolated to.
      84             :    */
      85           1 :   std::vector<bool> indices_interpolated_to_thus_far{};
      86             :   /*!
      87             :    * \brief Holds the element IDs of all elements that intersect with the
      88             :    * current iteration surface. Used to determine which elements will send
      89             :    * data for the next horizon find.
      90             :    */
      91           1 :   std::unordered_set<ElementId<3>> intersecting_element_ids{};
      92             :   /*!
      93             :    * \brief Offsets of newly interpolated points in the overall tensor (used as
      94             :    * memory buffer)
      95             :    */
      96           1 :   std::vector<size_t> offsets_of_newly_interpolated_points{};
      97             :   /*!
      98             :    * \brief Logical coordinates of newly interpolated points (used as memory
      99             :    * buffer)
     100             :    *
     101             :    * These `std::vector`s are used to reserve memory and then append points to
     102             :    * them as we find them in an element. The memory is reused for each element.
     103             :    * Then, a non-owning DataVector is created by pointing into this memory.
     104             :    * That's why this is a `std::array` of `std::vector`s, not vice versa.
     105             :    */
     106             :   std::array<std::vector<double>, 3>
     107           1 :       x_element_logical_of_newly_interpolated_points{};
     108             :   /*!
     109             :    * \brief Buffer for newly interpolated variables (used as memory buffer)
     110             :    */
     111           1 :   std::vector<double> newly_interpolated_vars_buffer{};
     112             : 
     113             :   /*!
     114             :    * \brief How many times we've tried to compute the coordinates for this
     115             :    * iteration.
     116             :    */
     117           1 :   size_t compute_coords_retries = 0;
     118             : 
     119             :   /*!
     120             :    * \brief Whether all points in `interpolated_vars` have been filled.
     121             :    */
     122           1 :   bool interpolation_is_complete() const;
     123             : 
     124           0 :   void reset_for_next_iteration();
     125             : 
     126             :   // NOLINTNEXTLINE(google-runtime-references)
     127           0 :   void pup(PUP::er& p);
     128             : };
     129             : 
     130             : template <typename Fr>
     131           0 : bool operator==(const Iteration<Fr>& lhs, const Iteration<Fr>& rhs);
     132             : template <typename Fr>
     133           0 : bool operator!=(const Iteration<Fr>& lhs, const Iteration<Fr>& rhs);
     134             : 
     135             : /*!
     136             :  * \brief Holds all data necessary for a single horizon find.
     137             :  *
     138             :  * \details This includes volume variables which persist for the entire horizon
     139             :  * find, and also interpolated variables that are updated for each iteration.
     140             :  */
     141             : template <typename Fr>
     142           1 : struct SingleTimeStorage {
     143             :   /*!
     144             :    * \brief Map between `ElementId`s and the volume variables from that element.
     145             :    */
     146           1 :   std::unordered_map<ElementId<3>, VolumeVariables<Fr>> all_volume_variables;
     147             :   /*!
     148             :    * \brief Elements in which we have found points to interpolate to in previous
     149             :    * iterations, to try first before searching all elements.
     150             :    *
     151             :    * This is not only a performance optimization, but also important for
     152             :    * robustness. If we try to interpolate from elements in a different order in
     153             :    * each iteration, then points that lie directly on element boundaries can
     154             :    * fluctuate in interpolated value, preventing convergence (see
     155             :    * https://github.com/sxs-collaboration/spectre/issues/3899).
     156             :    */
     157           1 :   std::vector<ElementId<3>> element_order{};
     158             :   /*!
     159             :    * \brief The `Iteration` data for the current fast flow iteration.
     160             :    */
     161           1 :   Iteration<Fr> current_iteration;
     162             :   /*!
     163             :    * \brief The previous iteration surface, used if interpolation fails.
     164             :    */
     165           1 :   ylm::Strahlkorper<Fr> previous_iteration_surface;
     166             :   /*!
     167             :    * \brief The `ah::Destination` for this horizon find.
     168             :    */
     169           1 :   Destination destination{};
     170             :   /*!
     171             :    * \brief Whether we have checked if the functions of time are up to date for
     172             :    * this horizon find.
     173             :    */
     174           1 :   bool time_is_ready = false;
     175             : 
     176             :   // NOLINTNEXTLINE(google-runtime-references)
     177           0 :   void pup(PUP::er& p);
     178             : };
     179             : 
     180             : template <typename Fr>
     181           0 : bool operator==(const SingleTimeStorage<Fr>& lhs,
     182             :                 const SingleTimeStorage<Fr>& rhs);
     183             : template <typename Fr>
     184           0 : bool operator!=(const SingleTimeStorage<Fr>& lhs,
     185             :                 const SingleTimeStorage<Fr>& rhs);
     186             : 
     187             : /*!
     188             :  * \brief The time and final surface for a previous horizon find.
     189             :  */
     190             : template <typename Fr>
     191           1 : struct PreviousSurface {
     192           0 :   PreviousSurface() = default;
     193           0 :   PreviousSurface(const LinkedMessageId<double>& time_in,
     194             :                   ylm::Strahlkorper<Fr> surface_in,
     195             :                   std::unordered_set<ElementId<3>> intersecting_element_ids_in);
     196             : 
     197           0 :   LinkedMessageId<double> time;
     198           0 :   ylm::Strahlkorper<Fr> surface;
     199           0 :   std::unordered_set<ElementId<3>> intersecting_element_ids;
     200             : 
     201             :   // NOLINTNEXTLINE(google-runtime-references)
     202           0 :   void pup(PUP::er& p);
     203             : };
     204             : 
     205             : template <typename Fr>
     206           0 : bool operator==(const PreviousSurface<Fr>& lhs, const PreviousSurface<Fr>& rhs);
     207             : template <typename Fr>
     208           0 : bool operator!=(const PreviousSurface<Fr>& lhs, const PreviousSurface<Fr>& rhs);
     209             : 
     210             : /*!
     211             :  * \brief Holds a previous surface and a lock that protects it.
     212             :  *
     213             :  * \details This is used to store and update a previous surface in the global
     214             :  * cache and allow multiple readers (elements) to access it simultaneously.
     215             :  */
     216             : template <typename Fr>
     217           1 : struct LockedPreviousSurface {
     218           0 :   std::optional<PreviousSurface<Fr>> surface;
     219             :   // Lock is mutable so it can be retrieved from the const global cache and put
     220             :   // in read-lock mode by elements.
     221             :   // NOLINTNEXTLINE(spectre-mutable)
     222           0 :   mutable Parallel::MultiReaderSpinlock lock;
     223             : 
     224           0 :   LockedPreviousSurface();
     225           0 :   explicit LockedPreviousSurface(const PreviousSurface<Fr>& rhs);
     226           0 :   LockedPreviousSurface(const LockedPreviousSurface& rhs);
     227           0 :   LockedPreviousSurface& operator=(const LockedPreviousSurface& rhs);
     228           0 :   LockedPreviousSurface(LockedPreviousSurface&& rhs);
     229           0 :   LockedPreviousSurface& operator=(LockedPreviousSurface&& rhs);
     230           0 :   ~LockedPreviousSurface() = default;
     231             : 
     232             :   // NOLINTNEXTLINE(google-runtime-references)
     233           0 :   void pup(PUP::er& p);
     234             : };
     235             : 
     236             : template <typename Fr>
     237           0 : bool operator==(const LockedPreviousSurface<Fr>& lhs,
     238             :                 const LockedPreviousSurface<Fr>& rhs);
     239             : template <typename Fr>
     240           0 : bool operator!=(const LockedPreviousSurface<Fr>& lhs,
     241             :                 const LockedPreviousSurface<Fr>& rhs);
     242             : 
     243             : }  // namespace ah::Storage