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/Structure/ObjectLabel.hpp"
19 : #include "IO/Observer/ObserverComponent.hpp"
20 : #include "IO/Observer/ReductionActions.hpp"
21 : #include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp"
22 : #include "Parallel/Invoke.hpp"
23 : #include "Parallel/Reduction.hpp"
24 : #include "ParallelAlgorithms/Events/Tags.hpp"
25 : #include "Utilities/ProtocolHelpers.hpp"
26 : #include "Utilities/StdArrayHelpers.hpp"
27 : #include "Utilities/TMPL.hpp"
28 :
29 : /// \cond
30 : class DataVector;
31 : template <size_t VolumeDim>
32 : class ElementId;
33 : template <size_t Dim>
34 : class Mesh;
35 : namespace control_system::Tags {
36 : struct WriteDataToDisk;
37 : } // namespace control_system::Tags
38 : namespace Parallel {
39 : template <typename Metavariables>
40 : class GlobalCache;
41 : } // namespace Parallel
42 : namespace domain::Tags {
43 : template <size_t Dim>
44 : struct Mesh;
45 : template <size_t Dim, typename Frame>
46 : struct Coordinates;
47 : } // namespace domain::Tags
48 : namespace grmhd::ValenciaDivClean::Tags {
49 : struct TildeD;
50 : }
51 : /// \endcond
52 :
53 0 : namespace control_system::measurements {
54 :
55 : template <typename ControlSystems>
56 : struct PostReductionSendBNSStarCentersToControlSystem;
57 :
58 0 : namespace Tags {
59 : /// \ingroup DataBoxTagsGroup
60 : /// \ingroup ControlSystemGroup
61 : /// DataBox tag for location of neutron star center (or more accurately, center
62 : /// of mass of the matter in the x>0 (label A) or x<0 (label B) region, in grid
63 : /// coordinates).
64 : template <::domain::ObjectLabel Center>
65 1 : struct NeutronStarCenter : db::SimpleTag {
66 0 : using type = std::array<double, 3>;
67 : };
68 : } // namespace Tags
69 :
70 : /*!
71 : * \brief Factored Center of Mass calculation (for easier testing)
72 : *
73 : * \details This function computes the integral of tildeD (assumed to be the
74 : * conservative baryon density in the inertial frame), as well as its first
75 : * moment in the grid frame. The integrals are limited to \f$x>0\f$ (label A) or
76 : * \f$x<0\f$ (label B).
77 : *
78 : * \param mass_a Integral of tildeD (x > 0)
79 : * \param mass_b Integral of tildeD (x < 0)
80 : * \param first_moment_A First moment of integral of tildeD (x > 0)
81 : * \param first_moment_B First moment of integral of tildeD (x < 0)
82 : * \param mesh The mesh
83 : * \param inv_det_jacobian The inverse determinant of the jacobian of the map
84 : * between logical and inertial coordinates \param tilde_d TildeD on the mesh
85 : * \param x_grid The coordinates in the grid frame
86 : */
87 1 : void center_of_mass_integral_on_element(
88 : const gsl::not_null<double*> mass_a, const gsl::not_null<double*> mass_b,
89 : const gsl::not_null<std::array<double, 3>*> first_moment_A,
90 : const gsl::not_null<std::array<double, 3>*> first_moment_B,
91 : const Mesh<3>& mesh, const Scalar<DataVector>& inv_det_jacobian,
92 : const Scalar<DataVector>& tilde_d,
93 : const tnsr::I<DataVector, 3, Frame::Grid>& x_grid);
94 :
95 : /*!
96 : * \brief Measurement providing the location of the center of mass of the matter
97 : * in the \f$x>0\f$ and \f$x<0\f$ regions (assumed to correspond to the center
98 : * of mass of the two neutron stars in a BNS merger).
99 : *
100 : * \details We use Events::Tags::ObserverXXX for tags that might need to be
101 : * retrieved from either the Subcell or DG grid.
102 : */
103 1 : struct BothNSCenters : tt::ConformsTo<protocols::Measurement> {
104 0 : struct FindTwoCenters : tt::ConformsTo<protocols::Submeasurement> {
105 0 : static std::string name() { return "BothNSCenters::FindTwoCenters"; }
106 : /// Unused tag needed to conform to the submeasurement protocol.
107 : template <typename ControlSystems>
108 1 : using interpolation_target_tag = void;
109 :
110 0 : using compute_tags_for_observation_box = tmpl::list<>;
111 :
112 : /// Tags for the arguments to the apply function.
113 1 : using argument_tags =
114 : tmpl::list<Events::Tags::ObserverMesh<3>,
115 : Events::Tags::ObserverDetInvJacobian<Frame::ElementLogical,
116 : Frame::Inertial>,
117 : grmhd::ValenciaDivClean::Tags::TildeD,
118 : Events::Tags::ObserverCoordinates<3, Frame::Grid>>;
119 :
120 : /// Calculate integrals needed for CoM computation on each element,
121 : /// then reduce the data.
122 : template <typename Metavariables, typename ParallelComponent,
123 : typename ControlSystems>
124 1 : static void apply(const Mesh<3>& mesh,
125 : const Scalar<DataVector>& inv_det_jacobian,
126 : const Scalar<DataVector>& tilde_d,
127 : const tnsr::I<DataVector, 3, Frame::Grid> x_grid,
128 : const LinkedMessageId<double>& measurement_id,
129 : Parallel::GlobalCache<Metavariables>& cache,
130 : const ElementId<3>& array_index,
131 : const ParallelComponent* const /*meta*/,
132 : ControlSystems /*meta*/) {
133 : // Initialize integrals and perform local calculations
134 : double mass_a = 0.;
135 : double mass_b = 0.;
136 : std::array<double, 3> first_moment_a = {0., 0., 0.};
137 : std::array<double, 3> first_moment_b = {0., 0., 0.};
138 : center_of_mass_integral_on_element(&mass_a, &mass_b, &first_moment_a,
139 : &first_moment_b, mesh,
140 : inv_det_jacobian, tilde_d, x_grid);
141 :
142 : // Reduction
143 : auto my_proxy = Parallel::get_parallel_component<ParallelComponent>(
144 : cache)[array_index];
145 : // We need a place to run RunCallback on... this does not need to be
146 : // the control system using the CoM data.
147 : auto& reduction_target_proxy = Parallel::get_parallel_component<
148 : ControlComponent<Metavariables, tmpl::front<ControlSystems>>>(cache);
149 : Parallel::ReductionData<
150 : Parallel::ReductionDatum<LinkedMessageId<double>,
151 : funcl::AssertEqual<>>,
152 : Parallel::ReductionDatum<double, funcl::Plus<>>,
153 : Parallel::ReductionDatum<double, funcl::Plus<>>,
154 : Parallel::ReductionDatum<std::array<double, 3>, funcl::Plus<>>,
155 : Parallel::ReductionDatum<std::array<double, 3>, funcl::Plus<>>>
156 : reduction_data{measurement_id, mass_a, mass_b, first_moment_a,
157 : first_moment_b};
158 : Parallel::contribute_to_reduction<
159 : PostReductionSendBNSStarCentersToControlSystem<ControlSystems>>(
160 : std::move(reduction_data), my_proxy, reduction_target_proxy);
161 : }
162 : };
163 : /// List of submeasurements used by this measurement -- only FindTwoCenters
164 : /// here.
165 1 : using submeasurements = tmpl::list<FindTwoCenters>;
166 : };
167 :
168 : /*!
169 : * \brief Simple action called after reduction of the center of mass data.
170 : *
171 : * \details `mass_a`, `mass_b`, `first_moment_a`, and `first_moment_b` will
172 : * contain the reduced data for the integral of the density (and its first
173 : * moment) in the x>=0 (A label) and x<0 (B label) regions. This action
174 : * calculates the center of mass in each region, and sends the result to the
175 : * control system.
176 : *
177 : * If the `control_system::Tags::WriteDataToDisk` tag is true, then this will
178 : * also write the centers of mass to the `/ControlSystems/BnsCenters` subfile of
179 : * the reductions h5 file. The columns of the file are
180 : *
181 : * - %Time
182 : * - Center_A_x
183 : * - Center_A_y
184 : * - Center_A_z
185 : * - Center_B_x
186 : * - Center_B_y
187 : * - Center_B_z
188 : */
189 : template <typename ControlSystems>
190 1 : struct PostReductionSendBNSStarCentersToControlSystem {
191 : template <typename ParallelComponent, typename DbTags, typename Metavariables,
192 : typename ArrayIndex>
193 0 : static void apply(db::DataBox<DbTags>& /*box*/,
194 : Parallel::GlobalCache<Metavariables>& cache,
195 : const ArrayIndex& /*array_index*/,
196 : const LinkedMessageId<double>& measurement_id,
197 : const double& mass_a, const double& mass_b,
198 : const std::array<double, 3>& first_moment_a,
199 : const std::array<double, 3>& first_moment_b) {
200 : // Function called after reduction of the CoM data.
201 : // Calculate CoM from integrals
202 : std::array<double, 3> center_a = first_moment_a / mass_a;
203 : std::array<double, 3> center_b = first_moment_b / mass_b;
204 : const auto center_databox = db::create<
205 : db::AddSimpleTags<Tags::NeutronStarCenter<::domain::ObjectLabel::A>,
206 : Tags::NeutronStarCenter<::domain::ObjectLabel::B>>>(
207 : center_a, center_b);
208 : // Send results to the control system(s)
209 : RunCallbacks<BothNSCenters::FindTwoCenters, ControlSystems>::apply(
210 : center_databox, cache, measurement_id);
211 :
212 : if (Parallel::get<control_system::Tags::WriteDataToDisk>(cache)) {
213 : std::vector<double> data_to_write{
214 : measurement_id.id, center_a[0], center_a[1], center_a[2],
215 : center_b[0], center_b[1], center_b[2]};
216 :
217 : auto& writer_proxy = Parallel::get_parallel_component<
218 : observers::ObserverWriter<Metavariables>>(cache);
219 :
220 : Parallel::threaded_action<
221 : observers::ThreadedActions::WriteReductionDataRow>(
222 : // Node 0 is always the writer
223 : writer_proxy[0], subfile_path_, legend_,
224 : std::make_tuple(std::move(data_to_write)));
225 : }
226 : }
227 :
228 : private:
229 0 : const static inline std::vector<std::string> legend_{
230 : "Time", "Center_A_x", "Center_A_y", "Center_A_z",
231 : "Center_B_x", "Center_B_y", "Center_B_z"};
232 0 : const static inline std::string subfile_path_{"/ControlSystems/BnsCenters"};
233 : };
234 :
235 : } // namespace control_system::measurements
|
