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