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 <array>
       8             : #include <cstddef>
       9             : #include <limits>
      10             : #include <set>
      11             : #include <variant>
      12             : #include <vector>
      13             : 
      14             : #include "DataStructures/DataBox/DataBox.hpp"
      15             : #include "DataStructures/Transpose.hpp"
      16             : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
      17             : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
      18             : #include "Options/String.hpp"
      19             : #include "Parallel/GlobalCache.hpp"
      20             : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
      21             : #include "ParallelAlgorithms/Interpolation/Protocols/ComputeTargetPoints.hpp"
      22             : #include "ParallelAlgorithms/Interpolation/Tags.hpp"
      23             : #include "ParallelAlgorithms/Interpolation/Targets/AngularOrdering.hpp"
      24             : #include "Utilities/Algorithm.hpp"
      25             : #include "Utilities/ErrorHandling/Assert.hpp"
      26             : #include "Utilities/PrettyType.hpp"
      27             : #include "Utilities/ProtocolHelpers.hpp"
      28             : 
      29             : /// \cond
      30             : class DataVector;
      31             : namespace PUP {
      32             : class er;
      33             : }  // namespace PUP
      34             : namespace db {
      35             : template <typename TagsList>
      36             : class DataBox;
      37             : }  // namespace db
      38             : namespace intrp {
      39             : namespace Tags {
      40             : template <typename TemporalId>
      41             : struct TemporalIds;
      42             : }  // namespace Tags
      43             : }  // namespace intrp
      44             : /// \endcond
      45             : 
      46             : namespace intrp {
      47             : 
      48             : namespace OptionHolders {
      49             : /// A series of concentric spherical surfaces.
      50             : ///
      51             : /// \details The parameter `LMax` sets the number of collocation points on
      52             : /// each spherical surface equal to `(l_max + 1) * (2 * l_max + 1)`. The
      53             : /// parameter `AngularOrdering` encodes the collocation ordering. For example,
      54             : /// the apparent horizon finder relies on spherepack routines that require
      55             : /// `Strahlkorper` for `AngularOrdering`, and using these surfaces for a CCE
      56             : /// worldtube requires `Cce` for `AngularOrdering`.
      57           1 : struct Sphere {
      58           0 :   struct LMax {
      59           0 :     using type = size_t;
      60           0 :     static constexpr Options::String help = {
      61             :         "The number of collocation points on each sphere will be equal to "
      62             :         "`(l_max + 1) * (2 * l_max + 1)`"};
      63             :   };
      64           0 :   struct Center {
      65           0 :     using type = std::array<double, 3>;
      66           0 :     static constexpr Options::String help = {"Center of every sphere"};
      67             :   };
      68           0 :   struct Radius {
      69           0 :     using type = std::variant<double, std::vector<double>>;
      70           0 :     static constexpr Options::String help = {"Radius of the sphere(s)"};
      71             :   };
      72           0 :   struct AngularOrdering {
      73           0 :     using type = intrp::AngularOrdering;
      74           0 :     static constexpr Options::String help = {
      75             :         "Chooses theta,phi ordering in 2d array"};
      76             :   };
      77           0 :   using options = tmpl::list<LMax, Center, Radius, AngularOrdering>;
      78           0 :   static constexpr Options::String help = {
      79             :       "An arbitrary number of spherical surface."};
      80           0 :   Sphere(const size_t l_max_in, const std::array<double, 3> center_in,
      81             :          const typename Radius::type& radius_in,
      82             :          intrp::AngularOrdering angular_ordering_in,
      83             :          const Options::Context& context = {});
      84             : 
      85           0 :   Sphere() = default;
      86             : 
      87             :   // NOLINTNEXTLINE(google-runtime-references)
      88           0 :   void pup(PUP::er& p);
      89             : 
      90           0 :   size_t l_max{0};
      91           0 :   std::array<double, 3> center{std::numeric_limits<double>::signaling_NaN()};
      92           0 :   std::set<double> radii;
      93           0 :   intrp::AngularOrdering angular_ordering;
      94             : };
      95             : 
      96           0 : bool operator==(const Sphere& lhs, const Sphere& rhs);
      97           0 : bool operator!=(const Sphere& lhs, const Sphere& rhs);
      98             : 
      99             : }  // namespace OptionHolders
     100             : 
     101             : namespace OptionTags {
     102             : template <typename InterpolationTargetTag>
     103           0 : struct Sphere {
     104           0 :   using type = OptionHolders::Sphere;
     105           0 :   static constexpr Options::String help{
     106             :       "Options for interpolation onto a sphere(s)."};
     107           0 :   static std::string name() {
     108             :     return pretty_type::name<InterpolationTargetTag>();
     109             :   }
     110           0 :   using group = InterpolationTargets;
     111             : };
     112             : }  // namespace OptionTags
     113             : 
     114             : namespace Tags {
     115             : template <typename InterpolationTargetTag>
     116           0 : struct Sphere : db::SimpleTag {
     117           0 :   using type = OptionHolders::Sphere;
     118           0 :   using option_tags = tmpl::list<OptionTags::Sphere<InterpolationTargetTag>>;
     119             : 
     120           0 :   static constexpr bool pass_metavariables = false;
     121           0 :   static OptionHolders::Sphere create_from_options(
     122             :       const OptionHolders::Sphere& option) {
     123             :     return option;
     124             :   }
     125             : };
     126             : 
     127             : template <typename Frame>
     128           0 : struct AllCoords : db::SimpleTag {
     129           0 :   using type = tnsr::I<DataVector, 3, Frame>;
     130             : };
     131             : }  // namespace Tags
     132             : 
     133             : namespace TargetPoints {
     134             : /// \brief Computes points on spherical surfaces.
     135             : ///
     136             : /// Conforms to the intrp::protocols::ComputeTargetPoints protocol
     137             : ///
     138             : /// For requirements on InterpolationTargetTag, see
     139             : /// intrp::protocols::InterpolationTargetTag
     140             : template <typename InterpolationTargetTag, typename Frame>
     141           1 : struct Sphere : tt::ConformsTo<intrp::protocols::ComputeTargetPoints> {
     142           0 :   using const_global_cache_tags =
     143             :       tmpl::list<Tags::Sphere<InterpolationTargetTag>>;
     144           0 :   using is_sequential = std::false_type;
     145           0 :   using frame = Frame;
     146             : 
     147           0 :   using simple_tags =
     148             :       tmpl::list<ylm::Tags::Strahlkorper<Frame>, Tags::AllCoords<Frame>>;
     149           0 :   using compute_tags = typename ylm::Tags::compute_items_tags<Frame>;
     150             : 
     151             :   template <typename DbTags, typename Metavariables>
     152           0 :   static void initialize(const gsl::not_null<db::DataBox<DbTags>*> box,
     153             :                          const Parallel::GlobalCache<Metavariables>& cache) {
     154             :     const auto& sphere =
     155             :         Parallel::get<Tags::Sphere<InterpolationTargetTag>>(cache);
     156             :     const size_t l_max = sphere.l_max;
     157             :     const auto& radii = sphere.radii;
     158             : 
     159             :     // Total number of points is number of points for one sphere times the
     160             :     // number of spheres we use.
     161             :     const size_t num_points = radii.size() * (l_max + 1) * (2 * l_max + 1);
     162             : 
     163             :     tnsr::I<DataVector, 3, Frame> all_coords{num_points};
     164             : 
     165             :     size_t index = 0;
     166             :     for (const double radius : radii) {
     167             :       ylm::Strahlkorper<Frame> strahlkorper(
     168             :           l_max, l_max, DataVector{(l_max + 1) * (2 * l_max + 1), radius},
     169             :           sphere.center);
     170             : 
     171             :       db::mutate<ylm::Tags::Strahlkorper<Frame>>(
     172             :           [&strahlkorper](const gsl::not_null<ylm::Strahlkorper<Frame>*>
     173             :                               local_strahlkorper) {
     174             :             *local_strahlkorper = std::move(strahlkorper);
     175             :           },
     176             :           box);
     177             : 
     178             :       // This copy is ok because it's just in initialization
     179             :       auto coords = db::get<ylm::Tags::CartesianCoords<Frame>>(*box);
     180             : 
     181             :       // If the angular ordering is Strahlkorper then we don't have to do
     182             :       // anything to the coords because they are already in the right order
     183             :       if (sphere.angular_ordering == intrp::AngularOrdering::Cce) {
     184             :         const auto physical_extents =
     185             :             strahlkorper.ylm_spherepack().physical_extents();
     186             :         auto transposed_coords =
     187             :             tnsr::I<DataVector, 3, Frame>(get<0>(coords).size());
     188             : 
     189             :         for (size_t i = 0; i < 3; ++i) {
     190             :           transpose(make_not_null(&transposed_coords.get(i)), coords.get(i),
     191             :                     physical_extents[0], physical_extents[1]);
     192             :         }
     193             :         coords = std::move(transposed_coords);
     194             :       }
     195             : 
     196             :       const size_t tmp_index = index;
     197             :       for (size_t i = 0; i < 3; ++i) {
     198             :         for (size_t local_index = 0; local_index < coords.get(i).size();
     199             :              local_index++, index++) {
     200             :           all_coords.get(i)[index] = coords.get(i)[local_index];
     201             :         }
     202             :         if (i != 2) {
     203             :           index = tmp_index;
     204             :         }
     205             :       }
     206             :     }
     207             : 
     208             :     // If this fails, there is a bug. Can't really test it
     209             :     // LCOV_EXCL_START
     210             :     ASSERT(index == all_coords.get(0).size(),
     211             :            "Didn't initialize points of Sphere target correctly. index = "
     212             :                << index << " all_coords.size() = " << all_coords.get(0).size());
     213             :     // LCOV_EXCL_STOP
     214             : 
     215             :     Initialization::mutate_assign<tmpl::list<Tags::AllCoords<Frame>>>(
     216             :         box, std::move(all_coords));
     217             :   }
     218             : 
     219             :   template <typename Metavariables, typename DbTags, typename TemporalId>
     220           0 :   static tnsr::I<DataVector, 3, Frame> points(
     221             :       const db::DataBox<DbTags>& box,
     222             :       const tmpl::type_<Metavariables>& metavariables_v,
     223             :       const TemporalId& /*temporal_id*/) {
     224             :     return points(box, metavariables_v);
     225             :   }
     226             :   template <typename Metavariables, typename DbTags>
     227           0 :   static tnsr::I<DataVector, 3, Frame> points(
     228             :       const db::DataBox<DbTags>& box,
     229             :       const tmpl::type_<Metavariables>& /*meta*/) {
     230             :     return db::get<Tags::AllCoords<Frame>>(box);
     231             :   }
     232             : };
     233             : }  // namespace TargetPoints
     234             : }  // namespace intrp