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