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 <string>
       8             : #include <tuple>
       9             : 
      10             : #include "DataStructures/DataBox/DataBox.hpp"
      11             : #include "DataStructures/DataBox/Prefixes.hpp"
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      14             : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
      15             : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
      16             : #include "IO/Observer/Helpers.hpp"
      17             : #include "IO/Observer/ObservationId.hpp"
      18             : #include "IO/Observer/ObserverComponent.hpp"
      19             : #include "IO/Observer/ReductionActions.hpp"
      20             : #include "IO/Observer/TypeOfObservation.hpp"
      21             : #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
      22             : #include "Parallel/GlobalCache.hpp"
      23             : #include "Parallel/Reduction.hpp"
      24             : #include "Time/Tags/TimeStepId.hpp"
      25             : #include "Utilities/ConstantExpressions.hpp"
      26             : #include "Utilities/TMPL.hpp"
      27             : #include "Utilities/TaggedTuple.hpp"
      28             : 
      29             : /// \cond
      30             : namespace Tags {
      31             : struct Time;
      32             : }  // namespace Tags
      33             : /// \endcond
      34             : 
      35             : namespace CurvedScalarWave::Worldtube::Actions {
      36             : 
      37             : /*!
      38             :  * \brief When Tags::ObserveCoefficientsTrigger is triggered, write the
      39             :  * coefficients of the Taylor expansion of the regular field as well as the
      40             :  * current particle's position, velocity and acceleration to file. A brief
      41             :  * summary is printed to stdout.
      42             :  */
      43           1 : struct ObserveWorldtubeSolution {
      44           0 :   using reduction_data = Parallel::ReductionData<
      45             :       Parallel::ReductionDatum<double, funcl::AssertEqual<>>,
      46             :       Parallel::ReductionDatum<std::vector<double>, funcl::AssertEqual<>>>;
      47             : 
      48           0 :   using observed_reduction_data_tags =
      49             :       observers::make_reduction_data_tags<tmpl::list<reduction_data>>;
      50           0 :   static constexpr size_t Dim = 3;
      51             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
      52             :             typename ArrayIndex, typename ActionList,
      53             :             typename ParallelComponent>
      54           0 :   static Parallel::iterable_action_return_t apply(
      55             :       db::DataBox<DbTagsList>& box,
      56             :       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
      57             :       Parallel::GlobalCache<Metavariables>& cache,
      58             :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
      59             :       const ParallelComponent* const /*meta*/) {
      60             :     if (db::get<Tags::ObserveCoefficientsTrigger>(box).is_triggered(box) and
      61             :         db::get<::Tags::TimeStepId>(box).substep() == 0) {
      62             :       const auto& inertial_particle_position =
      63             :           db::get<Tags::EvolvedPosition<Dim>>(box);
      64             :       const auto& particle_velocity = db::get<Tags::EvolvedVelocity<Dim>>(box);
      65             :       const auto& particle_acceleration =
      66             :           db::get<::Tags::dt<Tags::EvolvedVelocity<Dim>>>(box);
      67             :       const size_t expansion_order = db::get<Tags::ExpansionOrder>(box);
      68             :       const auto& psi_monopole = db::get<
      69             :           Stf::Tags::StfTensor<Tags::PsiWorldtube, 0, Dim, Frame::Inertial>>(
      70             :           box);
      71             :       const auto& dt_psi_monopole =
      72             :           db::get<Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 0, Dim,
      73             :                                        Frame::Inertial>>(box);
      74             : 
      75             :       // number of components in Taylor series
      76             :       const size_t num_coefs = ((expansion_order + 3) * (expansion_order + 2) *
      77             :                                 (expansion_order + 1)) /
      78             :                                6;
      79             :       // offset everything by the 3 * 3 components for position, velocity and
      80             :       // acceleration
      81             :       const size_t pos_offset = 9;
      82             :       std::vector<double> psi_coefs(2 * num_coefs + pos_offset);
      83             :       for (size_t i = 0; i < 3; ++i) {
      84             :         psi_coefs[i] = inertial_particle_position.get(i)[0];
      85             :         psi_coefs[i + 3] = particle_velocity.get(i)[0];
      86             :         psi_coefs[i + 6] = particle_acceleration.get(i)[0];
      87             :       }
      88             :       psi_coefs[0 + pos_offset] = get(psi_monopole);
      89             :       psi_coefs[num_coefs + pos_offset] = get(dt_psi_monopole);
      90             :       if (expansion_order > 0) {
      91             :         const auto& psi_dipole = db::get<
      92             :             Stf::Tags::StfTensor<Tags::PsiWorldtube, 1, Dim, Frame::Inertial>>(
      93             :             box);
      94             :         const auto& dt_psi_dipole =
      95             :             db::get<Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 1, Dim,
      96             :                                          Frame::Inertial>>(box);
      97             :         for (size_t i = 0; i < Dim; ++i) {
      98             :           psi_coefs[1 + i + pos_offset] = psi_dipole.get(i);
      99             :           psi_coefs[num_coefs + 1 + i + pos_offset] = dt_psi_dipole.get(i);
     100             :         }
     101             :       }
     102             :       const auto legend = [&expansion_order]() -> std::vector<std::string> {
     103             :         switch (expansion_order) {
     104             :           case (0):
     105             :             return {"Time",           "Position_x",     "Position_y",
     106             :                     "Position_z",     "Velocity_x",     "Velocity_y",
     107             :                     "Velocity_z",     "Acceleration_x", "Acceleration_y",
     108             :                     "Acceleration_z", "Psi0",           "dtPsi0"};
     109             :             break;
     110             :           case (1):
     111             :             return {"Time",           "Position_x",     "Position_y",
     112             :                     "Position_z",     "Velocity_x",     "Velocity_y",
     113             :                     "Velocity_z",     "Acceleration_x", "Acceleration_y",
     114             :                     "Acceleration_z", "Psi0",           "Psix",
     115             :                     "Psiy",           "Psiz",           "dtPsi0",
     116             :                     "dtPsix",         "dtPsiy",         "dtPsiz"};
     117             :             break;
     118             :           default:
     119             :             ERROR("requested invalid expansion order");
     120             :         }
     121             :       }();
     122             :       const auto current_time = db::get<::Tags::Time>(box);
     123             :       if (Parallel::get<Tags::Verbosity>(cache) >= ::Verbosity::Quiet) {
     124             :         Parallel::printf(
     125             :             "Time: %.16f, field value: %.16f, worldtube radius: %.16f, orbital "
     126             :             "radius: %.16f\n",
     127             :             current_time, psi_coefs[9], db::get<Tags::WorldtubeRadius>(box),
     128             :             get(magnitude(inertial_particle_position))[0]);
     129             :       }
     130             : 
     131             :       const auto observation_id =
     132             :           observers::ObservationId(current_time, "/Worldtube");
     133             :       auto& reduction_writer = Parallel::get_parallel_component<
     134             :           observers::ObserverWriter<Metavariables>>(cache);
     135             : 
     136             :       Parallel::threaded_action<
     137             :           observers::ThreadedActions::WriteReductionDataRow>(
     138             :           reduction_writer[0], std::string{"/PsiTaylorCoefs"}, legend,
     139             :           std::make_tuple(current_time, psi_coefs));
     140             :     }
     141             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     142             :   }
     143             : };
     144             : }  // namespace CurvedScalarWave::Worldtube::Actions