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::Grid>`
      55             :  *    - `Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 0, Dim,
      56             :  *                                     Frame::Grid>`
      57             :  *    - `Stf::Tags::StfTensor<Tags::PsiWorldtube, 1, Dim, Frame::Grid>`
      58             :  *    - `Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 1, Dim,
      59             :  *                                     Frame::Grid>`
      60             :  *    - `Stf::Tags::StfTensor<Tags::PsiWorldtube, 2, Dim, Frame::Grid>`
      61             :  *    - `Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 2, Dim,
      62             :  *                                     Frame::Grid>`
      63             :  */
      64           1 : struct ReceiveElementData {
      65           0 :   static constexpr size_t Dim = 3;
      66           0 :   using tags_list = tmpl::list<CurvedScalarWave::Tags::Psi,
      67             :                                ::Tags::dt<CurvedScalarWave::Tags::Psi>>;
      68           0 :   using inbox_tags = tmpl::list<
      69             :       ::CurvedScalarWave::Worldtube::Tags::SphericalHarmonicsInbox<Dim>>;
      70           0 :   using simple_tags = tmpl::list<
      71             :       Stf::Tags::StfTensor<Tags::PsiWorldtube, 0, Dim, Frame::Grid>,
      72             :       Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 0, Dim, Frame::Grid>,
      73             :       Stf::Tags::StfTensor<Tags::PsiWorldtube, 1, Dim, Frame::Grid>,
      74             :       Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 1, Dim, Frame::Grid>,
      75             :       Stf::Tags::StfTensor<Tags::PsiWorldtube, 2, Dim, Frame::Grid>,
      76             :       Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 2, Dim,
      77             :                            Frame::Grid>>;
      78             : 
      79             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
      80             :             typename ArrayIndex, typename ActionList,
      81             :             typename ParallelComponent>
      82           0 :   static Parallel::iterable_action_return_t apply(
      83             :       db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
      84             :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
      85             :       const ArrayIndex& /*array_index*/, ActionList /*meta*/,
      86             :       const ParallelComponent* const /*meta*/) {
      87             :     const size_t expected_number_of_senders =
      88             :         db::get<Tags::ElementFacesGridCoordinates<Dim>>(box).size();
      89             :     const auto& time_step_id = db::get<::Tags::TimeStepId>(box);
      90             :     auto& inbox = tuples::get<Tags::SphericalHarmonicsInbox<Dim>>(inboxes);
      91             :     if (inbox.count(time_step_id) == 0 or
      92             :         inbox.at(time_step_id).size() < expected_number_of_senders) {
      93             :       return {Parallel::AlgorithmExecution::Retry, std::nullopt};
      94             :     }
      95             :     ASSERT(inbox.at(time_step_id).size() == expected_number_of_senders,
      96             :            "Expected data from "
      97             :                << expected_number_of_senders << " senders, but received "
      98             :                << inbox.at(time_step_id).size() << " for TimeStepId "
      99             :                << time_step_id);
     100             :     const size_t order = db::get<Tags::ExpansionOrder>(box);
     101             :     const size_t num_modes = (order + 1) * (order + 1);
     102             : 
     103             :     Variables<tags_list> external_ylm_coefs{num_modes, 0.};
     104             :     for (const auto& [_, element_ylm_coefs] : inbox.at(time_step_id)) {
     105             :       external_ylm_coefs += element_ylm_coefs;
     106             :     }
     107             :     const double wt_radius = db::get<Tags::ExcisionSphere<Dim>>(box).radius();
     108             :     external_ylm_coefs /= wt_radius * wt_radius;
     109             : 
     110             :     DataVector& psi_ylm_coefs =
     111             :         get(get<CurvedScalarWave::Tags::Psi>(external_ylm_coefs));
     112             :     DataVector& dt_psi_ylm_coefs =
     113             :         get(get<::Tags::dt<CurvedScalarWave::Tags::Psi>>(external_ylm_coefs));
     114             : 
     115             :     const ModalVector psi_ylm_l0(&psi_ylm_coefs[0], 1);
     116             :     const ModalVector dt_psi_ylm_l0(&dt_psi_ylm_coefs[0], 1);
     117             :     tnsr::i<double, Dim, Frame::Grid> psi_stf_l1{};
     118             :     tnsr::i<double, Dim, Frame::Grid> dt_psi_stf_l1{};
     119             :     tnsr::ii<double, Dim, Frame::Grid> psi_stf_l2{};
     120             :     tnsr::ii<double, Dim, Frame::Grid> dt_psi_stf_l2{};
     121             :     if (order > 0) {
     122             :       ModalVector psi_ylm_l1(&psi_ylm_coefs[1], 3);
     123             :       ModalVector dt_psi_ylm_l1(&dt_psi_ylm_coefs[1], 3);
     124             :       psi_ylm_l1 /= wt_radius;
     125             :       dt_psi_ylm_l1 /= wt_radius;
     126             :       psi_stf_l1 = ylm::ylm_to_stf_1<Frame::Grid>(psi_ylm_l1);
     127             :       dt_psi_stf_l1 = ylm::ylm_to_stf_1<Frame::Grid>(dt_psi_ylm_l1);
     128             :       if (order > 1) {
     129             :         ModalVector psi_ylm_l2(&psi_ylm_coefs[4], 5);
     130             :         ModalVector dt_psi_ylm_l2(&dt_psi_ylm_coefs[4], 5);
     131             :         psi_ylm_l2 /= wt_radius * wt_radius;
     132             :         dt_psi_ylm_l2 /= wt_radius * wt_radius;
     133             :         psi_stf_l2 = ylm::ylm_to_stf_2<Frame::Grid>(psi_ylm_l2);
     134             :         dt_psi_stf_l2 = ylm::ylm_to_stf_2<Frame::Grid>(dt_psi_ylm_l2);
     135             :       }
     136             :     }
     137             : 
     138             :     ::Initialization::mutate_assign<simple_tags>(
     139             :         make_not_null(&box), ylm::ylm_to_stf_0(psi_ylm_l0),
     140             :         ylm::ylm_to_stf_0(dt_psi_ylm_l0), std::move(psi_stf_l1),
     141             :         std::move(dt_psi_stf_l1), std::move(psi_stf_l2),
     142             :         std::move(dt_psi_stf_l2));
     143             :     inbox.erase(time_step_id);
     144             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     145             :   }
     146             : };
     147             : }  // namespace CurvedScalarWave::Worldtube::Actions