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/DataBox/Tag.hpp"
16 : #include "DataStructures/LinkedMessageId.hpp"
17 : #include "DataStructures/Tensor/TypeAliases.hpp"
18 : #include "Domain/BlockLogicalCoordinates.hpp"
19 : #include "Domain/Creators/Tags/Domain.hpp"
20 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
21 : #include "Domain/FunctionsOfTime/Tags.hpp"
22 : #include "Domain/Structure/ObjectLabel.hpp"
23 : #include "IO/Logging/Verbosity.hpp"
24 : #include "IO/Observer/ObserverComponent.hpp"
25 : #include "IO/Observer/ReductionActions.hpp"
26 : #include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp"
27 : #include "Parallel/GlobalCache.hpp"
28 : #include "Parallel/Invoke.hpp"
29 : #include "Parallel/Reduction.hpp"
30 : #include "ParallelAlgorithms/Actions/FunctionsOfTimeAreReady.hpp"
31 : #include "ParallelAlgorithms/Events/Tags.hpp"
32 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
33 : #include "Utilities/Algorithm.hpp"
34 : #include "Utilities/GetOutput.hpp"
35 : #include "Utilities/ProtocolHelpers.hpp"
36 : #include "Utilities/Serialization/CharmPupable.hpp"
37 : #include "Utilities/StdArrayHelpers.hpp"
38 : #include "Utilities/TMPL.hpp"
39 :
40 : /// \cond
41 : class DataVector;
42 : template <size_t VolumeDim>
43 : class ElementId;
44 : template <size_t Dim>
45 : class Mesh;
46 : namespace control_system::Tags {
47 : struct WriteDataToDisk;
48 : } // namespace control_system::Tags
49 : namespace Parallel {
50 : template <typename Metavariables>
51 : class GlobalCache;
52 : } // namespace Parallel
53 : namespace domain::Tags {
54 : template <size_t Dim>
55 : struct Mesh;
56 : template <size_t Dim, typename Frame>
57 : struct Coordinates;
58 : } // namespace domain::Tags
59 : namespace grmhd::ValenciaDivClean::Tags {
60 : struct TildeD;
61 : } // namespace grmhd::ValenciaDivClean::Tags
62 : namespace Tags {
63 : struct PreviousTriggerTime;
64 : } // namespace Tags
65 : /// \endcond
66 :
67 : namespace control_system {
68 1 : namespace measurements {
69 : /// \cond
70 : template <typename ControlSystems>
71 : struct PostReductionSendBNSStarCentersToControlSystem;
72 : /// \endcond
73 :
74 : /*!
75 : * \brief Factored Center of Mass calculation (for easier testing)
76 : *
77 : * \details This function computes the integral of tildeD (assumed to be the
78 : * conservative baryon density in the inertial frame), as well as its first
79 : * moment in the grid frame. The integrals are limited to \f$x>0\f$ (label A) or
80 : * \f$x<0\f$ (label B).
81 : *
82 : * \param mass_a Integral of tildeD (x > 0)
83 : * \param mass_b Integral of tildeD (x < 0)
84 : * \param first_moment_A First moment of integral of tildeD (x > 0)
85 : * \param first_moment_B First moment of integral of tildeD (x < 0)
86 : * \param mesh The mesh
87 : * \param inv_det_jacobian The inverse determinant of the jacobian of the map
88 : * between logical and inertial coordinates \param tilde_d TildeD on the mesh
89 : * \param x_grid The coordinates in the grid frame
90 : */
91 1 : void center_of_mass_integral_on_element(
92 : const gsl::not_null<double*> mass_a, const gsl::not_null<double*> mass_b,
93 : const gsl::not_null<std::array<double, 3>*> first_moment_A,
94 : const gsl::not_null<std::array<double, 3>*> first_moment_B,
95 : const Mesh<3>& mesh, const Scalar<DataVector>& inv_det_jacobian,
96 : const Scalar<DataVector>& tilde_d,
97 : const tnsr::I<DataVector, 3, Frame::Grid>& x_grid);
98 : } // namespace measurements
99 :
100 : /*!
101 : * \brief An `::Event` that computes the center of mass for $x > 0$ and $x < 0$
102 : * where $x$ is the in the `Frame::Grid`.
103 : *
104 : * \details See
105 : * `control_system::measurements::grid_center_of_mass_integral_on_element` for
106 : * the calculation of the CoM. This event then does a reduction and calls
107 : * `control_system::PostReductionSendBNSStarCentersToControlSystem` as a post
108 : * reduction callback.
109 : *
110 : * \tparam ControlSystems `tmpl::list` of all control systems that use this
111 : * event.
112 : */
113 : template <typename ControlSystems>
114 1 : class BNSEvent : public ::Event {
115 : public:
116 : /// \cond
117 : // LCOV_EXCL_START
118 : explicit BNSEvent(CkMigrateMessage* /*unused*/) {}
119 : using PUP::able::register_constructor;
120 : WRAPPED_PUPable_decl_template(BNSEvent); // NOLINT
121 : // LCOV_EXCL_STOP
122 : /// \endcond
123 :
124 : // This event is created during control system initialization, not
125 : // from the input file.
126 0 : static constexpr bool factory_creatable = false;
127 0 : BNSEvent() = default;
128 :
129 0 : using compute_tags_for_observation_box = tmpl::list<>;
130 :
131 0 : using return_tags = tmpl::list<>;
132 0 : using argument_tags =
133 : tmpl::list<::Tags::Time, ::Tags::PreviousTriggerTime,
134 : ::Events::Tags::ObserverMesh<3>,
135 : ::Events::Tags::ObserverDetInvJacobian<Frame::ElementLogical,
136 : Frame::Inertial>,
137 : grmhd::ValenciaDivClean::Tags::TildeD,
138 : ::Events::Tags::ObserverCoordinates<3, Frame::Grid>>;
139 :
140 : template <typename Metavariables, typename ArrayIndex,
141 : typename ParallelComponent>
142 0 : void operator()(const double time, const std::optional<double>& previous_time,
143 : const Mesh<3>& mesh,
144 : const Scalar<DataVector>& inv_det_jacobian,
145 : const Scalar<DataVector>& tilde_d,
146 : const tnsr::I<DataVector, 3, Frame::Grid>& x_grid,
147 : Parallel::GlobalCache<Metavariables>& cache,
148 : const ArrayIndex& array_index,
149 : const ParallelComponent* const /*meta*/,
150 : const ObservationValue& /*observation_value*/) const {
151 : const LinkedMessageId<double> measurement_id{time, previous_time};
152 :
153 : // Initialize integrals and perform local calculations
154 : double mass_a = 0.;
155 : double mass_b = 0.;
156 : std::array<double, 3> first_moment_a = {0., 0., 0.};
157 : std::array<double, 3> first_moment_b = {0., 0., 0.};
158 : measurements::center_of_mass_integral_on_element(
159 : &mass_a, &mass_b, &first_moment_a, &first_moment_b, mesh,
160 : inv_det_jacobian, tilde_d, x_grid);
161 :
162 : // Reduction
163 : auto my_proxy =
164 : Parallel::get_parallel_component<ParallelComponent>(cache)[array_index];
165 : // We need a place to run RunCallback on... this does not need to be
166 : // the control system using the CoM data.
167 : auto& reduction_target_proxy = Parallel::get_parallel_component<
168 : ControlComponent<Metavariables, tmpl::front<ControlSystems>>>(cache);
169 : Parallel::ReductionData<
170 : Parallel::ReductionDatum<LinkedMessageId<double>, funcl::AssertEqual<>>,
171 : Parallel::ReductionDatum<double, funcl::Plus<>>,
172 : Parallel::ReductionDatum<double, funcl::Plus<>>,
173 : Parallel::ReductionDatum<std::array<double, 3>, funcl::Plus<>>,
174 : Parallel::ReductionDatum<std::array<double, 3>, funcl::Plus<>>>
175 : reduction_data{measurement_id, mass_a, mass_b, first_moment_a,
176 : first_moment_b};
177 : Parallel::contribute_to_reduction<
178 : measurements::PostReductionSendBNSStarCentersToControlSystem<
179 : ControlSystems>>(std::move(reduction_data), my_proxy,
180 : reduction_target_proxy);
181 : }
182 :
183 0 : using is_ready_argument_tags = tmpl::list<>;
184 :
185 : template <typename Metavariables, typename ArrayIndex, typename Component>
186 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
187 : const ArrayIndex& /*array_index*/,
188 : const Component* const /*component*/) const {
189 : return true;
190 : }
191 :
192 1 : bool needs_evolved_variables() const override { return true; }
193 : };
194 :
195 : /// \cond
196 : template <typename ControlSystems>
197 : PUP::able::PUP_ID BNSEvent<ControlSystems>::my_PUP_ID = 0; // NOLINT
198 : /// \endcond
199 :
200 : namespace measurements {
201 1 : namespace Tags {
202 : /// \ingroup DataBoxTagsGroup
203 : /// \ingroup ControlSystemGroup
204 : /// DataBox tag for location of neutron star center (or more accurately, center
205 : /// of mass of the matter in the x>0 (label A) or x<0 (label B) region, in grid
206 : /// coordinates).
207 : template <::domain::ObjectLabel Center, typename Fr>
208 1 : struct NeutronStarCenter : db::SimpleTag {
209 0 : using type = std::array<double, 3>;
210 : };
211 :
212 : /// \ingroup DataBoxTagsGroup
213 : /// \ingroup ControlSystemGroup
214 : /// DataBox tag for location of center of mass of all of the matter.
215 : template <typename Fr>
216 1 : struct SystemCenterOfMass : db::SimpleTag {
217 0 : using type = std::array<double, 3>;
218 : };
219 : } // namespace Tags
220 :
221 : /*!
222 : * \brief Measurement providing the location of the center of mass of the matter
223 : * in the \f$x>0\f$ and \f$x<0\f$ regions (assumed to correspond to the center
224 : * of mass of the two neutron stars in a BNS merger).
225 : *
226 : * \details We use Events::Tags::ObserverXXX for tags that might need to be
227 : * retrieved from either the Subcell or DG grid.
228 : */
229 1 : struct BothNSCenters : tt::ConformsTo<protocols::Measurement> {
230 0 : struct FindTwoCenters : tt::ConformsTo<protocols::Submeasurement> {
231 0 : static std::string name() { return "BothNSCenters::FindTwoCenters"; }
232 : /// Unused aliases needed to conform to the submeasurement protocol.
233 : template <typename ControlSystems>
234 1 : using interpolation_target_tag = void;
235 : template <typename ControlSystems>
236 0 : using horizon_metavars = void;
237 :
238 : template <typename ControlSystems>
239 0 : using event = BNSEvent<ControlSystems>;
240 : };
241 : /// List of submeasurements used by this measurement -- only FindTwoCenters
242 : /// here.
243 1 : using submeasurements = tmpl::list<FindTwoCenters>;
244 : };
245 :
246 : /*!
247 : * \brief Simple action called after reduction of the center of mass data.
248 : *
249 : * \details `mass_a`, `mass_b`, `first_moment_a`, and `first_moment_b` will
250 : * contain the reduced data for the integral of the density (and its first
251 : * moment) in the x>=0 (A label) and x<0 (B label) regions. This action
252 : * calculates the center of mass in each region, and sends the result to the
253 : * control system.
254 : *
255 : * If the `control_system::Tags::WriteDataToDisk` tag is true, then this will
256 : * also write the centers of mass to the `/ControlSystems/BnsCenters` subfile of
257 : * the reductions h5 file. The columns of the file are
258 : *
259 : * - %Time
260 : * - Center_A_x
261 : * - Center_A_y
262 : * - Center_A_z
263 : * - Center_B_x
264 : * - Center_B_y
265 : * - Center_B_z
266 : */
267 : template <typename ControlSystems>
268 1 : struct PostReductionSendBNSStarCentersToControlSystem {
269 : template <typename ParallelComponent, typename DbTags, typename Metavariables,
270 : typename ArrayIndex>
271 0 : static void apply(db::DataBox<DbTags>& /*box*/,
272 : Parallel::GlobalCache<Metavariables>& cache,
273 : const ArrayIndex& /*array_index*/,
274 : const LinkedMessageId<double>& measurement_id,
275 : const double& mass_a, const double& mass_b,
276 : const std::array<double, 3>& first_moment_a,
277 : const std::array<double, 3>& first_moment_b) {
278 : // Function called after reduction of the CoM data.
279 : // Calculate CoM from integrals
280 : std::array<double, 3> grid_center_a = first_moment_a / mass_a;
281 : std::array<double, 3> grid_center_b = first_moment_b / mass_b;
282 : std::array<double, 3> grid_center_total =
283 : (first_moment_a + first_moment_b) / (mass_a + mass_b);
284 :
285 : // To convert grid coords to inertial coords, we must find the block that
286 : // these coords are in and use that grid to inertial map
287 : const Domain<3>& domain = Parallel::get<domain::Tags::Domain<3>>(cache);
288 : const domain::FunctionsOfTimeMap& functions_of_time =
289 : Parallel::get<domain::Tags::FunctionsOfTime>(cache);
290 : tnsr::I<DataVector, 3, Frame::Grid> grid_tnsr_center{};
291 : for (size_t i = 0; i < 3; i++) {
292 : grid_tnsr_center.get(i) =
293 : DataVector{gsl::at(grid_center_a, i), gsl::at(grid_center_b, i),
294 : gsl::at(grid_center_total, i)};
295 : }
296 :
297 : const auto block_logical_coords = block_logical_coordinates(
298 : domain, grid_tnsr_center, measurement_id.id, functions_of_time);
299 :
300 : ASSERT(alg::all_of(block_logical_coords,
301 : [](const auto& logical_coord_holder) {
302 : return logical_coord_holder.has_value();
303 : }),
304 : "Grid centers of BNS ("
305 : << grid_center_a << ", " << grid_center_b << ", "
306 : << grid_center_total
307 : << ") could not be mapped to the logical frame.");
308 :
309 : const auto& blocks = domain.blocks();
310 :
311 : ASSERT(block_logical_coords.size() == 3,
312 : "There should be exactly 3 block logical coordinates for the two "
313 : "centers of the BNS plus the system COM, but instead there are "
314 : << block_logical_coords.size());
315 :
316 : std::array<double, 3> inertial_center_a{};
317 : std::array<double, 3> inertial_center_b{};
318 : std::array<double, 3> inertial_center_total{};
319 :
320 : for (size_t n = 0; n < block_logical_coords.size(); n++) {
321 : const auto& logical_coord_holder = block_logical_coords[n];
322 : const auto& block_id = logical_coord_holder.value().id;
323 :
324 : const auto& block = blocks[block_id.get_index()];
325 : const auto& grid_to_inertial_map =
326 : block.moving_mesh_grid_to_inertial_map();
327 :
328 : const auto inertial_center = grid_to_inertial_map(
329 : tnsr::I<double, 3, Frame::Grid>{
330 : n == 0 ? grid_center_a
331 : : (n == 1 ? grid_center_b : grid_center_total)},
332 : measurement_id.id, functions_of_time);
333 :
334 : auto& center_to_set =
335 : n == 0 ? inertial_center_a
336 : : (n == 1 ? inertial_center_b : inertial_center_total);
337 : for (size_t i = 0; i < 3; i++) {
338 : gsl::at(center_to_set, i) = inertial_center.get(i);
339 : }
340 : }
341 :
342 : const auto center_databox = db::create<db::AddSimpleTags<
343 : Tags::NeutronStarCenter<::domain::ObjectLabel::A, Frame::Grid>,
344 : Tags::NeutronStarCenter<::domain::ObjectLabel::B, Frame::Grid>,
345 : Tags::NeutronStarCenter<::domain::ObjectLabel::A, Frame::Inertial>,
346 : Tags::NeutronStarCenter<::domain::ObjectLabel::B, Frame::Inertial>,
347 : Tags::SystemCenterOfMass<Frame::Grid>,
348 : Tags::SystemCenterOfMass<Frame::Inertial>>>(
349 : grid_center_a, grid_center_b, inertial_center_a, inertial_center_b,
350 : grid_center_total, inertial_center_total);
351 :
352 : // Send results to the control system(s)
353 : RunCallbacks<BothNSCenters::FindTwoCenters, ControlSystems>::apply(
354 : center_databox, cache, measurement_id);
355 :
356 : if (Parallel::get<control_system::Tags::WriteDataToDisk>(cache)) {
357 : std::vector<double> grid_data_to_write{
358 : measurement_id.id, grid_center_a[0], grid_center_a[1],
359 : grid_center_a[2], grid_center_b[0], grid_center_b[1],
360 : grid_center_b[2], grid_center_total[0], grid_center_total[1],
361 : grid_center_total[2]};
362 :
363 : std::vector<double> inertial_data_to_write{measurement_id.id};
364 : inertial_data_to_write.insert(inertial_data_to_write.end(),
365 : inertial_center_a.begin(),
366 : inertial_center_a.end());
367 : inertial_data_to_write.insert(inertial_data_to_write.end(),
368 : inertial_center_b.begin(),
369 : inertial_center_b.end());
370 : inertial_data_to_write.insert(inertial_data_to_write.end(),
371 : inertial_center_total.begin(),
372 : inertial_center_total.end());
373 :
374 : auto& writer_proxy = Parallel::get_parallel_component<
375 : observers::ObserverWriter<Metavariables>>(cache);
376 :
377 : Parallel::threaded_action<
378 : observers::ThreadedActions::WriteReductionDataRow>(
379 : // Node 0 is always the writer
380 : writer_proxy[0], grid_subfile_path_, legend_,
381 : std::make_tuple(std::move(grid_data_to_write)));
382 : Parallel::threaded_action<
383 : observers::ThreadedActions::WriteReductionDataRow>(
384 : // Node 0 is always the writer
385 : writer_proxy[0], inertial_subfile_path_, legend_,
386 : std::make_tuple(std::move(inertial_data_to_write)));
387 : }
388 : }
389 :
390 : private:
391 0 : const static inline std::vector<std::string> legend_{
392 : "Time", "Center_A_x", "Center_A_y",
393 : "Center_A_z", "Center_B_x", "Center_B_y",
394 : "Center_B_z", "Center_System_x", "Center_System_y",
395 : "Center_System_z"};
396 0 : const static inline std::string grid_subfile_path_{
397 : "/ControlSystems/BnsGridCenters"};
398 0 : const static inline std::string inertial_subfile_path_{
399 : "/ControlSystems/BnsInertialCenters"};
400 : };
401 : } // namespace measurements
402 : } // namespace control_system
|
