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 <limits>
       8              : #include <memory>
       9              : #include <string>
      10              : #include <tuple>
      11              : #include <vector>
      12              : 
      13              : #include "Evolution/Systems/Cce/AnalyticSolutions/WorldtubeData.hpp"
      14              : #include "Evolution/Systems/Cce/BoundaryData.hpp"
      15              : #include "Evolution/Systems/Cce/Tags.hpp"
      16              : #include "Evolution/Systems/Cce/WorldtubeDataManager.hpp"
      17              : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
      18              : #include "IO/Observer/ObserverComponent.hpp"
      19              : #include "IO/Observer/ReductionActions.hpp"
      20              : #include "NumericalAlgorithms/Interpolation/SpanInterpolator.hpp"
      21              : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCoefficients.hpp"
      22              : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshTransform.hpp"
      23              : #include "Parallel/GlobalCache.hpp"
      24              : #include "Parallel/Invoke.hpp"
      25              : #include "Parallel/Local.hpp"
      26              : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      27              : #include "Utilities/MakeString.hpp"
      28              : #include "Utilities/Serialization/CharmPupable.hpp"
      29              : #include "Utilities/System/ParallelInfo.hpp"
      30              : 
      31              : namespace Cce {
      32            1 : namespace Tags {
      33              : /// \cond
      34              : struct ObservationLMax;
      35              : /// \endcond
      36              : }  // namespace Tags
      37              : 
      38              : /// A boundary data manager that constructs the desired boundary data into
      39              : /// the `Variables` from the data provided by the analytic solution.
      40            1 : class AnalyticBoundaryDataManager {
      41              :  public:
      42              :   // charm needs an empty constructor.
      43            0 :   AnalyticBoundaryDataManager() = default;
      44              : 
      45            0 :   AnalyticBoundaryDataManager(
      46              :       size_t l_max, double extraction_radius,
      47              :       std::unique_ptr<Solutions::WorldtubeData> generator,
      48              :       double du_dr_j_stencil_step);
      49              : 
      50              :   /*!
      51              :    * \brief Update the `boundary_data_variables` entries for all tags in
      52              :    * `Tags::characteristic_worldtube_boundary_tags` except `Du<Dr<BondiJ>>`
      53              :    * (see `populate_du_dr_j`) to the boundary data from the analytic solution
      54              :    * at  `time`.
      55              :    *
      56              :    * \details This class retrieves metric boundary data from the
      57              :    * `Cce::Solutions::WorldtubeData` derived class that represents an analytic
      58              :    * solution, then dispatches to `Cce::create_bondi_boundary_data()` to
      59              :    * construct the Bondi values into the provided `Variables`
      60              :    */
      61            1 :   bool populate_hypersurface_boundary_data(
      62              :       gsl::not_null<Variables<
      63              :           Tags::characteristic_worldtube_boundary_tags<Tags::BoundaryValue>>*>
      64              :           boundary_data_variables,
      65              :       double time) const;
      66              : 
      67              :   /// Use `observers::ThreadedActions::WriteSimpleData` to output the expected
      68              :   /// news at `time` from the analytic data to dataset `/expected_news.dat`.
      69              :   template <typename ParallelComponent, typename Metavariables>
      70            1 :   void write_news(Parallel::GlobalCache<Metavariables>& cache,
      71              :                   double time) const;
      72              : 
      73              :   /*!
      74              :    * \brief Update the `Du<Dr<BondiJ>>` entry of `boundary_data_variables` at
      75              :    * `time` if an interpolator was given with `set_du_dr_j_interpolator`.
      76              :    *
      77              :    * \details Samples `Dr<BondiJ>` from the analytic solution at
      78              :    * `2 * interpolator->required_number_of_points_before_and_after()` times
      79              :    * spaced by `Cce::OptionTags::DuDrJStencilStep` and centered on `time`, and
      80              :    * differentiates it with the interpolator.
      81              :    */
      82            1 :   void populate_du_dr_j(
      83              :       gsl::not_null<Variables<
      84              :           Tags::characteristic_worldtube_boundary_tags<Tags::BoundaryValue>>*>
      85              :           boundary_data_variables,
      86              :       double time) const;
      87              : 
      88              :   /// Set the interpolator used by `populate_du_dr_j`. Without one (the
      89              :   /// default, or after passing `nullptr`), `populate_du_dr_j` does nothing.
      90            1 :   void set_du_dr_j_interpolator(
      91              :       std::unique_ptr<intrp::SpanInterpolator> interpolator);
      92              : 
      93            0 :   double get_du_dr_j_stencil_step() const { return du_dr_j_stencil_step_; }
      94              : 
      95            0 :   size_t get_l_max() const { return l_max_; }
      96              : 
      97            0 :   const Solutions::WorldtubeData& get_generator() const { return *generator_; }
      98              : 
      99              :   /// Serialization for Charm++.
     100              :   // NOLINTNEXTLINE(google-runtime-references)
     101            1 :   void pup(PUP::er& p);
     102              : 
     103              :  private:
     104            0 :   size_t l_max_ = 0;
     105            0 :   std::unique_ptr<Solutions::WorldtubeData> generator_;
     106            0 :   double extraction_radius_ = std::numeric_limits<double>::signaling_NaN();
     107            0 :   std::unique_ptr<intrp::SpanInterpolator> du_dr_j_interpolator_;
     108            0 :   double du_dr_j_stencil_step_ = std::numeric_limits<double>::signaling_NaN();
     109              : };
     110              : 
     111              : template <typename ParallelComponent, typename Metavariables>
     112              : void AnalyticBoundaryDataManager::write_news(
     113              :     Parallel::GlobalCache<Metavariables>& cache, const double time) const {
     114              :   const auto news = get<Tags::News>(
     115              :       generator_->variables(l_max_, time, tmpl::list<Tags::News>{}));
     116              :   const size_t observation_l_max = Parallel::get<Tags::ObservationLMax>(cache);
     117              :   std::vector<double> data_to_write(2 * square(observation_l_max + 1) + 1);
     118              :   std::vector<std::string> file_legend;
     119              :   file_legend.reserve(2 * square(observation_l_max + 1) + 1);
     120              :   file_legend.emplace_back("time");
     121              :   for (int i = 0; i <= static_cast<int>(observation_l_max); ++i) {
     122              :     for (int j = -i; j <= i; ++j) {
     123              :       file_legend.push_back(MakeString{} << "Real Y_" << i << "," << j);
     124              :       file_legend.push_back(MakeString{} << "Imag Y_" << i << "," << j);
     125              :     }
     126              :   }
     127              :   const ComplexModalVector goldberg_modes =
     128              :       Spectral::Swsh::libsharp_to_goldberg_modes(
     129              :           Spectral::Swsh::swsh_transform(l_max_, 1, get(news)), l_max_)
     130              :           .data();
     131              :   data_to_write[0] = time;
     132              :   for (size_t i = 0; i < square(observation_l_max + 1); ++i) {
     133              :     data_to_write[2 * i + 1] = real(goldberg_modes[i]);
     134              :     data_to_write[2 * i + 2] = imag(goldberg_modes[i]);
     135              :   }
     136              :   auto observer_proxy = Parallel::get_parallel_component<
     137              :       observers::ObserverWriter<Metavariables>>(cache)[0];
     138              :   const std::string prefix =
     139              :       generator_->use_noninertial_news() ? "Noninertial_" : "";
     140              :   Parallel::threaded_action<observers::ThreadedActions::WriteReductionDataRow>(
     141              :       observer_proxy, "/Cce/News_" + prefix + "expected"s, file_legend,
     142              :       std::make_tuple(data_to_write));
     143              : }
     144              : }  // namespace Cce