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 : #include <vector>
10 :
11 : #include "ControlSystem/Component.hpp"
12 : #include "ControlSystem/Protocols/Measurement.hpp"
13 : #include "ControlSystem/Protocols/Submeasurement.hpp"
14 : #include "ControlSystem/RunCallbacks.hpp"
15 : #include "DataStructures/LinkedMessageId.hpp"
16 : #include "DataStructures/Tensor/TypeAliases.hpp"
17 : #include "Domain/FunctionsOfTime/Tags.hpp"
18 : #include "Domain/Structure/ObjectLabel.hpp"
19 : #include "IO/Logging/Verbosity.hpp"
20 : #include "IO/Observer/ObserverComponent.hpp"
21 : #include "IO/Observer/ReductionActions.hpp"
22 : #include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp"
23 : #include "Parallel/GlobalCache.hpp"
24 : #include "Parallel/Invoke.hpp"
25 : #include "Parallel/Reduction.hpp"
26 : #include "ParallelAlgorithms/Actions/FunctionsOfTimeAreReady.hpp"
27 : #include "ParallelAlgorithms/Events/Tags.hpp"
28 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
29 : #include "Utilities/GetOutput.hpp"
30 : #include "Utilities/ProtocolHelpers.hpp"
31 : #include "Utilities/Serialization/CharmPupable.hpp"
32 : #include "Utilities/StdArrayHelpers.hpp"
33 : #include "Utilities/TMPL.hpp"
34 :
35 : /// \cond
36 : class DataVector;
37 : template <size_t VolumeDim>
38 : class ElementId;
39 : template <size_t Dim>
40 : class Mesh;
41 : namespace control_system::Tags {
42 : struct WriteDataToDisk;
43 : } // namespace control_system::Tags
44 : namespace Parallel {
45 : template <typename Metavariables>
46 : class GlobalCache;
47 : } // namespace Parallel
48 : namespace domain::Tags {
49 : template <size_t Dim>
50 : struct Mesh;
51 : template <size_t Dim, typename Frame>
52 : struct Coordinates;
53 : } // namespace domain::Tags
54 : namespace grmhd::ValenciaDivClean::Tags {
55 : struct TildeD;
56 : } // namespace grmhd::ValenciaDivClean::Tags
57 : namespace Tags {
58 : struct PreviousTriggerTime;
59 : } // namespace Tags
60 : /// \endcond
61 :
62 : namespace control_system {
63 0 : namespace measurements {
64 : /// \cond
65 : template <typename ControlSystems>
66 : struct PostReductionSendBNSStarCentersToControlSystem;
67 : /// \endcond
68 :
69 : /*!
70 : * \brief Factored Center of Mass calculation (for easier testing)
71 : *
72 : * \details This function computes the integral of tildeD (assumed to be the
73 : * conservative baryon density in the inertial frame), as well as its first
74 : * moment in the grid frame. The integrals are limited to \f$x>0\f$ (label A) or
75 : * \f$x<0\f$ (label B).
76 : *
77 : * \param mass_a Integral of tildeD (x > 0)
78 : * \param mass_b Integral of tildeD (x < 0)
79 : * \param first_moment_A First moment of integral of tildeD (x > 0)
80 : * \param first_moment_B First moment of integral of tildeD (x < 0)
81 : * \param mesh The mesh
82 : * \param inv_det_jacobian The inverse determinant of the jacobian of the map
83 : * between logical and inertial coordinates \param tilde_d TildeD on the mesh
84 : * \param x_grid The coordinates in the grid frame
85 : */
86 1 : void center_of_mass_integral_on_element(
87 : const gsl::not_null<double*> mass_a, const gsl::not_null<double*> mass_b,
88 : const gsl::not_null<std::array<double, 3>*> first_moment_A,
89 : const gsl::not_null<std::array<double, 3>*> first_moment_B,
90 : const Mesh<3>& mesh, const Scalar<DataVector>& inv_det_jacobian,
91 : const Scalar<DataVector>& tilde_d,
92 : const tnsr::I<DataVector, 3, Frame::Grid>& x_grid);
93 : } // namespace measurements
94 :
95 : /*!
96 : * \brief An `::Event` that computes the center of mass for $x > 0$ and $x < 0$
97 : * where $x$ is the in the `Frame::Grid`.
98 : *
99 : * \details See
100 : * `control_system::measurements::center_of_mass_integral_on_element` for the
101 : * calculation of the CoM. This event then does a reduction and calls
102 : * `control_system::PostReductionSendBNSStarCentersToControlSystem` as a post
103 : * reduction callback.
104 : *
105 : * \tparam ControlSystems `tmpl::list` of all control systems that use this
106 : * event.
107 : */
108 : template <typename ControlSystems>
109 1 : class BNSEvent : public ::Event {
110 : public:
111 : /// \cond
112 : // LCOV_EXCL_START
113 : explicit BNSEvent(CkMigrateMessage* /*unused*/) {}
114 : using PUP::able::register_constructor;
115 : WRAPPED_PUPable_decl_template(BNSEvent); // NOLINT
116 : // LCOV_EXCL_STOP
117 : /// \endcond
118 :
119 : // This event is created during control system initialization, not
120 : // from the input file.
121 0 : static constexpr bool factory_creatable = false;
122 0 : BNSEvent() = default;
123 :
124 0 : using compute_tags_for_observation_box = tmpl::list<>;
125 :
126 0 : using return_tags = tmpl::list<>;
127 0 : using argument_tags =
128 : tmpl::list<::Tags::Time, ::Tags::PreviousTriggerTime,
129 : ::Events::Tags::ObserverMesh<3>,
130 : ::Events::Tags::ObserverDetInvJacobian<Frame::ElementLogical,
131 : Frame::Inertial>,
132 : grmhd::ValenciaDivClean::Tags::TildeD,
133 : ::Events::Tags::ObserverCoordinates<3, Frame::Grid>>;
134 :
135 : template <typename Metavariables, typename ArrayIndex,
136 : typename ParallelComponent>
137 0 : void operator()(const double time, const std::optional<double>& previous_time,
138 : const Mesh<3>& mesh,
139 : const Scalar<DataVector>& inv_det_jacobian,
140 : const Scalar<DataVector>& tilde_d,
141 : const tnsr::I<DataVector, 3, Frame::Grid>& x_grid,
142 : Parallel::GlobalCache<Metavariables>& cache,
143 : const ArrayIndex& array_index,
144 : const ParallelComponent* const /*meta*/,
145 : const ObservationValue& /*observation_value*/) const {
146 : const LinkedMessageId<double> measurement_id{time, previous_time};
147 :
148 : // Initialize integrals and perform local calculations
149 : double mass_a = 0.;
150 : double mass_b = 0.;
151 : std::array<double, 3> first_moment_a = {0., 0., 0.};
152 : std::array<double, 3> first_moment_b = {0., 0., 0.};
153 : measurements::center_of_mass_integral_on_element(
154 : &mass_a, &mass_b, &first_moment_a, &first_moment_b, mesh,
155 : inv_det_jacobian, tilde_d, x_grid);
156 :
157 : // Reduction
158 : auto my_proxy =
159 : Parallel::get_parallel_component<ParallelComponent>(cache)[array_index];
160 : // We need a place to run RunCallback on... this does not need to be
161 : // the control system using the CoM data.
162 : auto& reduction_target_proxy = Parallel::get_parallel_component<
163 : ControlComponent<Metavariables, tmpl::front<ControlSystems>>>(cache);
164 : Parallel::ReductionData<
165 : Parallel::ReductionDatum<LinkedMessageId<double>, funcl::AssertEqual<>>,
166 : Parallel::ReductionDatum<double, funcl::Plus<>>,
167 : Parallel::ReductionDatum<double, funcl::Plus<>>,
168 : Parallel::ReductionDatum<std::array<double, 3>, funcl::Plus<>>,
169 : Parallel::ReductionDatum<std::array<double, 3>, funcl::Plus<>>>
170 : reduction_data{measurement_id, mass_a, mass_b, first_moment_a,
171 : first_moment_b};
172 : Parallel::contribute_to_reduction<
173 : measurements::PostReductionSendBNSStarCentersToControlSystem<
174 : ControlSystems>>(std::move(reduction_data), my_proxy,
175 : reduction_target_proxy);
176 : }
177 :
178 0 : using is_ready_argument_tags = tmpl::list<>;
179 :
180 : template <typename Metavariables, typename ArrayIndex, typename Component>
181 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
182 : const ArrayIndex& /*array_index*/,
183 : const Component* const /*component*/) const {
184 : return true;
185 : }
186 :
187 1 : bool needs_evolved_variables() const override { return true; }
188 : };
189 :
190 : /// \cond
191 : template <typename ControlSystems>
192 : PUP::able::PUP_ID BNSEvent<ControlSystems>::my_PUP_ID = 0; // NOLINT
193 : /// \endcond
194 :
195 : namespace measurements {
196 0 : namespace Tags {
197 : /// \ingroup DataBoxTagsGroup
198 : /// \ingroup ControlSystemGroup
199 : /// DataBox tag for location of neutron star center (or more accurately, center
200 : /// of mass of the matter in the x>0 (label A) or x<0 (label B) region, in grid
201 : /// coordinates).
202 : template <::domain::ObjectLabel Center>
203 1 : struct NeutronStarCenter : db::SimpleTag {
204 0 : using type = std::array<double, 3>;
205 : };
206 : } // namespace Tags
207 :
208 : /*!
209 : * \brief Measurement providing the location of the center of mass of the matter
210 : * in the \f$x>0\f$ and \f$x<0\f$ regions (assumed to correspond to the center
211 : * of mass of the two neutron stars in a BNS merger).
212 : *
213 : * \details We use Events::Tags::ObserverXXX for tags that might need to be
214 : * retrieved from either the Subcell or DG grid.
215 : */
216 1 : struct BothNSCenters : tt::ConformsTo<protocols::Measurement> {
217 0 : struct FindTwoCenters : tt::ConformsTo<protocols::Submeasurement> {
218 0 : static std::string name() { return "BothNSCenters::FindTwoCenters"; }
219 : /// Unused tag needed to conform to the submeasurement protocol.
220 : template <typename ControlSystems>
221 1 : using interpolation_target_tag = void;
222 :
223 : template <typename ControlSystems>
224 0 : using event = BNSEvent<ControlSystems>;
225 : };
226 : /// List of submeasurements used by this measurement -- only FindTwoCenters
227 : /// here.
228 1 : using submeasurements = tmpl::list<FindTwoCenters>;
229 : };
230 :
231 : /*!
232 : * \brief Simple action called after reduction of the center of mass data.
233 : *
234 : * \details `mass_a`, `mass_b`, `first_moment_a`, and `first_moment_b` will
235 : * contain the reduced data for the integral of the density (and its first
236 : * moment) in the x>=0 (A label) and x<0 (B label) regions. This action
237 : * calculates the center of mass in each region, and sends the result to the
238 : * control system.
239 : *
240 : * If the `control_system::Tags::WriteDataToDisk` tag is true, then this will
241 : * also write the centers of mass to the `/ControlSystems/BnsCenters` subfile of
242 : * the reductions h5 file. The columns of the file are
243 : *
244 : * - %Time
245 : * - Center_A_x
246 : * - Center_A_y
247 : * - Center_A_z
248 : * - Center_B_x
249 : * - Center_B_y
250 : * - Center_B_z
251 : */
252 : template <typename ControlSystems>
253 1 : struct PostReductionSendBNSStarCentersToControlSystem {
254 : template <typename ParallelComponent, typename DbTags, typename Metavariables,
255 : typename ArrayIndex>
256 0 : static void apply(db::DataBox<DbTags>& /*box*/,
257 : Parallel::GlobalCache<Metavariables>& cache,
258 : const ArrayIndex& /*array_index*/,
259 : const LinkedMessageId<double>& measurement_id,
260 : const double& mass_a, const double& mass_b,
261 : const std::array<double, 3>& first_moment_a,
262 : const std::array<double, 3>& first_moment_b) {
263 : // Function called after reduction of the CoM data.
264 : // Calculate CoM from integrals
265 : std::array<double, 3> center_a = first_moment_a / mass_a;
266 : std::array<double, 3> center_b = first_moment_b / mass_b;
267 : const auto center_databox = db::create<
268 : db::AddSimpleTags<Tags::NeutronStarCenter<::domain::ObjectLabel::A>,
269 : Tags::NeutronStarCenter<::domain::ObjectLabel::B>>>(
270 : center_a, center_b);
271 : // Send results to the control system(s)
272 : RunCallbacks<BothNSCenters::FindTwoCenters, ControlSystems>::apply(
273 : center_databox, cache, measurement_id);
274 :
275 : if (Parallel::get<control_system::Tags::WriteDataToDisk>(cache)) {
276 : std::vector<double> data_to_write{
277 : measurement_id.id, center_a[0], center_a[1], center_a[2],
278 : center_b[0], center_b[1], center_b[2]};
279 :
280 : auto& writer_proxy = Parallel::get_parallel_component<
281 : observers::ObserverWriter<Metavariables>>(cache);
282 :
283 : Parallel::threaded_action<
284 : observers::ThreadedActions::WriteReductionDataRow>(
285 : // Node 0 is always the writer
286 : writer_proxy[0], subfile_path_, legend_,
287 : std::make_tuple(std::move(data_to_write)));
288 : }
289 : }
290 :
291 : private:
292 0 : const static inline std::vector<std::string> legend_{
293 : "Time", "Center_A_x", "Center_A_y", "Center_A_z",
294 : "Center_B_x", "Center_B_y", "Center_B_z"};
295 0 : const static inline std::string subfile_path_{"/ControlSystems/BnsCenters"};
296 : };
297 : } // namespace measurements
298 : } // namespace control_system
|
