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/TaggedTuple.hpp"
14 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
15 : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
16 : #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
17 : #include "IO/Observer/Helpers.hpp"
18 : #include "IO/Observer/ObservationId.hpp"
19 : #include "IO/Observer/ObserverComponent.hpp"
20 : #include "IO/Observer/ReductionActions.hpp"
21 : #include "IO/Observer/TypeOfObservation.hpp"
22 : #include "NumericalAlgorithms/Strahlkorper/Tags.hpp"
23 : #include "Parallel/GlobalCache.hpp"
24 : #include "Parallel/Reduction.hpp"
25 : #include "Time/Tags/TimeStepId.hpp"
26 : #include "Utilities/ConstantExpressions.hpp"
27 : #include "Utilities/TMPL.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
|
