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 <unordered_map>
       9             : 
      10             : #include "DataStructures/DataBox/DataBox.hpp"
      11             : #include "DataStructures/DataBox/Prefixes.hpp"
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "DataStructures/ModalVector.hpp"
      14             : #include "DataStructures/Variables.hpp"
      15             : #include "Domain/Structure/ElementId.hpp"
      16             : #include "Domain/Tags.hpp"
      17             : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Inboxes.hpp"
      18             : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
      19             : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
      20             : #include "NumericalAlgorithms/SphericalHarmonics/YlmToStf.hpp"
      21             : #include "Parallel/AlgorithmExecution.hpp"
      22             : #include "Parallel/GlobalCache.hpp"
      23             : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
      24             : #include "Time/TimeStepId.hpp"
      25             : #include "Utilities/ErrorHandling/Assert.hpp"
      26             : #include "Utilities/Gsl.hpp"
      27             : #include "Utilities/TMPL.hpp"
      28             : 
      29             : /// \cond
      30             : namespace Tags {
      31             : struct TimeStepId;
      32             : }  // namespace Tags
      33             : /// \endcond
      34             : 
      35             : namespace CurvedScalarWave::Worldtube::Actions {
      36             : /*!
      37             :  * \brief Adds up the spherical harmonic projections from the different elements
      38             :  * abutting the worldtube.
      39             :  *
      40             :  * \details This action currently assumes that there is no h-refinement
      41             :  * ocurring in the elements abutting the worldtubes. This could be accounted for
      42             :  * by checking that data from at least one element has been sent from each
      43             :  * abutting block and then using its `ElementId` to figure out the current
      44             :  * refinement level and therefore how many elements are expected to send data
      45             :  * for each block.
      46             :  *
      47             :  * DataBox:
      48             :  * - Uses:
      49             :  *    - `Worldtube::Tags::ExpansionOrder`
      50             :  *    - `Worldtube::Tags::ExcisionSphere`
      51             :  *    - `Worldtube::Tags::ElementFacesGridCoordinates`
      52             :  *    - `Tags::TimeStepId`
      53             :  * - Mutates:
      54             :  *    - `Stf::Tags::StfTensor<Tags::PsiWorldtube, 0, Dim, Frame::Inertial>`
      55             :  *    - `Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 0, Dim,
      56             :  *                                     Frame::Inertial>`
      57             :  *    - `Stf::Tags::StfTensor<Tags::PsiWorldtube, 1, Dim, Frame::Inertial>`
      58             :  *    - `Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 1, Dim,
      59             :  *                                     Frame::Inertial>`
      60             :  */
      61           1 : struct ReceiveElementData {
      62           0 :   static constexpr size_t Dim = 3;
      63           0 :   using tags_list = tmpl::list<CurvedScalarWave::Tags::Psi,
      64             :                                ::Tags::dt<CurvedScalarWave::Tags::Psi>>;
      65           0 :   using inbox_tags = tmpl::list<
      66             :       ::CurvedScalarWave::Worldtube::Tags::SphericalHarmonicsInbox<Dim>>;
      67           0 :   using simple_tags = tmpl::list<
      68             :       Stf::Tags::StfTensor<Tags::PsiWorldtube, 0, Dim, Frame::Inertial>,
      69             :       Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 0, Dim,
      70             :                            Frame::Inertial>,
      71             :       Stf::Tags::StfTensor<Tags::PsiWorldtube, 1, Dim, Frame::Inertial>,
      72             :       Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 1, Dim,
      73             :                            Frame::Inertial>>;
      74             : 
      75             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
      76             :             typename ArrayIndex, typename ActionList,
      77             :             typename ParallelComponent>
      78           0 :   static Parallel::iterable_action_return_t apply(
      79             :       db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
      80             :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
      81             :       const ArrayIndex& /*array_index*/, ActionList /*meta*/,
      82             :       const ParallelComponent* const /*meta*/) {
      83             :     const size_t expected_number_of_senders =
      84             :         db::get<Tags::ElementFacesGridCoordinates<Dim>>(box).size();
      85             :     const auto& time_step_id = db::get<::Tags::TimeStepId>(box);
      86             :     auto& inbox = tuples::get<Tags::SphericalHarmonicsInbox<Dim>>(inboxes);
      87             :     if (inbox.count(time_step_id) == 0 or
      88             :         inbox.at(time_step_id).size() < expected_number_of_senders) {
      89             :       return {Parallel::AlgorithmExecution::Retry, std::nullopt};
      90             :     }
      91             :     ASSERT(inbox.at(time_step_id).size() == expected_number_of_senders,
      92             :            "Expected data from "
      93             :                << expected_number_of_senders << " senders, but received "
      94             :                << inbox.at(time_step_id).size() << " for TimeStepId "
      95             :                << time_step_id);
      96             :     const size_t order = db::get<Tags::ExpansionOrder>(box);
      97             :     const size_t num_modes = (order + 1) * (order + 1);
      98             : 
      99             :     Variables<tags_list> external_ylm_coefs{num_modes, 0.};
     100             :     for (const auto& [_, element_ylm_coefs] : inbox.at(time_step_id)) {
     101             :       external_ylm_coefs += element_ylm_coefs;
     102             :     }
     103             :     const double wt_radius = db::get<Tags::ExcisionSphere<Dim>>(box).radius();
     104             :     external_ylm_coefs /= wt_radius * wt_radius;
     105             : 
     106             :     DataVector& psi_ylm_coefs =
     107             :         get(get<CurvedScalarWave::Tags::Psi>(external_ylm_coefs));
     108             :     DataVector& dt_psi_ylm_coefs =
     109             :         get(get<::Tags::dt<CurvedScalarWave::Tags::Psi>>(external_ylm_coefs));
     110             : 
     111             :     const ModalVector psi_ylm_l0(&psi_ylm_coefs[0], 1);
     112             :     const ModalVector dt_psi_ylm_l0(&dt_psi_ylm_coefs[0], 1);
     113             :     tnsr::i<double, Dim, Frame::Inertial> psi_stf_l1{};
     114             :     tnsr::i<double, Dim, Frame::Inertial> dt_psi_stf_l1{};
     115             :     if (order > 0) {
     116             :       ModalVector psi_ylm_l1(&psi_ylm_coefs[1], 3);
     117             :       ModalVector dt_psi_ylm_l1(&dt_psi_ylm_coefs[1], 3);
     118             :       psi_ylm_l1 /= wt_radius;
     119             :       dt_psi_ylm_l1 /= wt_radius;
     120             :       psi_stf_l1 = ylm::ylm_to_stf_1<Frame::Inertial>(psi_ylm_l1);
     121             :       dt_psi_stf_l1 = ylm::ylm_to_stf_1<Frame::Inertial>(dt_psi_ylm_l1);
     122             :     }
     123             : 
     124             :     ::Initialization::mutate_assign<simple_tags>(
     125             :         make_not_null(&box), ylm::ylm_to_stf_0(psi_ylm_l0),
     126             :         ylm::ylm_to_stf_0(dt_psi_ylm_l0), std::move(psi_stf_l1),
     127             :         std::move(dt_psi_stf_l1));
     128             :     inbox.erase(time_step_id);
     129             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     130             :   }
     131             : };
     132             : }  // namespace CurvedScalarWave::Worldtube::Actions