Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <optional>
7 : #include <pup.h>
8 : #include <string>
9 : #include <vector>
10 :
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/Tensor/Tensor.hpp"
13 : #include "Domain/Structure/ElementId.hpp"
14 : #include "Domain/Tags.hpp"
15 : #include "Domain/Tags/FaceNormal.hpp"
16 : #include "Domain/Tags/Faces.hpp"
17 : #include "Elliptic/Systems/Xcts/Tags.hpp"
18 : #include "IO/Observer/GetSectionObservationKey.hpp"
19 : #include "IO/Observer/ObservationId.hpp"
20 : #include "IO/Observer/ObserverComponent.hpp"
21 : #include "IO/Observer/ReductionActions.hpp"
22 : #include "IO/Observer/TypeOfObservation.hpp"
23 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
24 : #include "Options/String.hpp"
25 : #include "Parallel/ArrayIndex.hpp"
26 : #include "Parallel/GlobalCache.hpp"
27 : #include "Parallel/Reduction.hpp"
28 : #include "Parallel/TypeTraits.hpp"
29 : #include "ParallelAlgorithms/Events/Tags.hpp"
30 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
31 : #include "Utilities/Serialization/CharmPupable.hpp"
32 :
33 0 : namespace Events {
34 :
35 : /// @{
36 : /*!
37 : * \brief Computes the ADM integrals locally (within one element).
38 : *
39 : * To get the total ADM integrals, the results need to be summed over in a
40 : * reduction.
41 : */
42 1 : void local_adm_integrals(
43 : gsl::not_null<Scalar<double>*> adm_mass,
44 : gsl::not_null<tnsr::I<double, 3>*> adm_linear_momentum,
45 : gsl::not_null<Scalar<double>*> adm_angular_momentum_z,
46 : gsl::not_null<tnsr::I<double, 3>*> center_of_mass,
47 : const Scalar<DataVector>& conformal_factor,
48 : const Scalar<DataVector>& conformal_factor_minus_one,
49 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
50 : const tnsr::ii<DataVector, 3>& conformal_metric,
51 : const tnsr::II<DataVector, 3>& inv_conformal_metric,
52 : const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
53 : const tnsr::i<DataVector, 3>& conformal_christoffel_contracted,
54 : const tnsr::ii<DataVector, 3>& spatial_metric,
55 : const tnsr::II<DataVector, 3>& inv_spatial_metric,
56 : const tnsr::ii<DataVector, 3>& extrinsic_curvature,
57 : const Scalar<DataVector>& trace_extrinsic_curvature,
58 : const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
59 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
60 : Frame::Inertial>& inv_jacobian,
61 : const Mesh<3>& mesh, const Element<3>& element,
62 : const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals);
63 : /// @}
64 :
65 : /// @{
66 : /*!
67 : * \brief Observe ADM integrals after the XCTS solve.
68 : *
69 : * The surface integrals are taken over the outer boundary, which is defined as
70 : * the domain boundary in the upper logical zeta direction.
71 : *
72 : * Writes reduction quantities:
73 : * - Number of points in the domain
74 : * - ADM mass
75 : * - ADM linear momentum
76 : * - ADM angular momentum (z-component)
77 : * - Center of mass
78 : */
79 : template <typename ArraySectionIdTag = void>
80 1 : class ObserveAdmIntegrals : public Event {
81 : private:
82 0 : using ReductionData = Parallel::ReductionData<
83 : // Number of points
84 : Parallel::ReductionDatum<size_t, funcl::Plus<>>,
85 : // ADM Mass
86 : Parallel::ReductionDatum<double, funcl::Plus<>>,
87 : // ADM Linear Momentum (x-component)
88 : Parallel::ReductionDatum<double, funcl::Plus<>>,
89 : // ADM Linear Momentum (y-component)
90 : Parallel::ReductionDatum<double, funcl::Plus<>>,
91 : // ADM Linear Momentum (z-component)
92 : Parallel::ReductionDatum<double, funcl::Plus<>>,
93 : // ADM Angular Momentum (z-component)
94 : Parallel::ReductionDatum<double, funcl::Plus<>>,
95 : // Center of Mass (x-component)
96 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
97 : std::index_sequence<1>>,
98 : // Center of Mass (y-component)
99 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
100 : std::index_sequence<1>>,
101 : // Center of Mass (z-component)
102 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
103 : std::index_sequence<1>>>;
104 :
105 : public:
106 : /// \cond
107 : explicit ObserveAdmIntegrals(CkMigrateMessage* msg) : Event(msg) {}
108 : using PUP::able::register_constructor;
109 : WRAPPED_PUPable_decl_template(ObserveAdmIntegrals); // NOLINT
110 : /// \endcond
111 :
112 0 : using options = tmpl::list<>;
113 0 : static constexpr Options::String help =
114 : "Observe ADM integrals after the XCTS solve.\n"
115 : "\n"
116 : "Writes reduction quantities:\n"
117 : "- Number of points in the domain\n"
118 : "- ADM mass\n"
119 : "- ADM linear momentum\n"
120 : "- ADM angular momentum (z-component)\n"
121 : "- Center of mass";
122 :
123 0 : ObserveAdmIntegrals() = default;
124 :
125 0 : using observed_reduction_data_tags =
126 : observers::make_reduction_data_tags<tmpl::list<ReductionData>>;
127 :
128 0 : using compute_tags_for_observation_box = tmpl::list<>;
129 :
130 0 : using return_tags = tmpl::list<>;
131 :
132 0 : using argument_tags = tmpl::list<
133 : Xcts::Tags::ConformalFactor<DataVector>,
134 : Xcts::Tags::ConformalFactorMinusOne<DataVector>,
135 : ::Tags::deriv<Xcts::Tags::ConformalFactor<DataVector>, tmpl::size_t<3>,
136 : Frame::Inertial>,
137 : Xcts::Tags::ConformalMetric<DataVector, 3, Frame::Inertial>,
138 : Xcts::Tags::InverseConformalMetric<DataVector, 3, Frame::Inertial>,
139 : Xcts::Tags::ConformalChristoffelSecondKind<DataVector, 3,
140 : Frame::Inertial>,
141 : Xcts::Tags::ConformalChristoffelContracted<DataVector, 3,
142 : Frame::Inertial>,
143 : gr::Tags::SpatialMetric<DataVector, 3, Frame::Inertial>,
144 : gr::Tags::InverseSpatialMetric<DataVector, 3, Frame::Inertial>,
145 : gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame::Inertial>,
146 : gr::Tags::TraceExtrinsicCurvature<DataVector>,
147 : domain::Tags::Coordinates<3, Frame::Inertial>,
148 : domain::Tags::InverseJacobian<3, Frame::ElementLogical, Frame::Inertial>,
149 : domain::Tags::Mesh<3>, domain::Tags::Element<3>,
150 : domain::Tags::Faces<3, domain::Tags::FaceNormal<3>>,
151 : ::Tags::ObservationBox>;
152 :
153 : template <typename DataBoxType, typename ComputeTagsList,
154 : typename Metavariables, typename ArrayIndex,
155 : typename ParallelComponent>
156 0 : void operator()(
157 : const Scalar<DataVector>& conformal_factor,
158 : const Scalar<DataVector>& conformal_factor_minus_one,
159 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
160 : const tnsr::ii<DataVector, 3>& conformal_metric,
161 : const tnsr::II<DataVector, 3>& inv_conformal_metric,
162 : const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
163 : const tnsr::i<DataVector, 3>& conformal_christoffel_contracted,
164 : const tnsr::ii<DataVector, 3>& spatial_metric,
165 : const tnsr::II<DataVector, 3>& inv_spatial_metric,
166 : const tnsr::ii<DataVector, 3>& extrinsic_curvature,
167 : const Scalar<DataVector>& trace_extrinsic_curvature,
168 : const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
169 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
170 : Frame::Inertial>& inv_jacobian,
171 : const Mesh<3>& mesh, const Element<3>& element,
172 : const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals,
173 : const ObservationBox<DataBoxType, ComputeTagsList>& box,
174 : Parallel::GlobalCache<Metavariables>& cache,
175 : const ArrayIndex& array_index, const ParallelComponent* const /*meta*/,
176 : const ObservationValue& observation_value) const {
177 : // Skip observation on elements that are not part of a section
178 : const std::optional<std::string> section_observation_key =
179 : observers::get_section_observation_key<ArraySectionIdTag>(box);
180 : if (not section_observation_key.has_value()) {
181 : return;
182 : }
183 : const std::string subfile_path = subfile_path_ + *section_observation_key;
184 :
185 : Scalar<double> adm_mass;
186 : tnsr::I<double, 3> adm_linear_momentum;
187 : Scalar<double> adm_angular_momentum_z;
188 : tnsr::I<double, 3> center_of_mass;
189 : local_adm_integrals(
190 : make_not_null(&adm_mass), make_not_null(&adm_linear_momentum),
191 : make_not_null(&adm_angular_momentum_z), make_not_null(¢er_of_mass),
192 : conformal_factor, conformal_factor_minus_one, deriv_conformal_factor,
193 : conformal_metric, inv_conformal_metric,
194 : conformal_christoffel_second_kind, conformal_christoffel_contracted,
195 : spatial_metric, inv_spatial_metric, extrinsic_curvature,
196 : trace_extrinsic_curvature, inertial_coords, inv_jacobian, mesh, element,
197 : conformal_face_normals);
198 :
199 : // Save components of linear momentum as reduction data
200 : ReductionData reduction_data{
201 : mesh.number_of_grid_points(), get(adm_mass),
202 : get<0>(adm_linear_momentum), get<1>(adm_linear_momentum),
203 : get<2>(adm_linear_momentum), get(adm_angular_momentum_z),
204 : get<0>(center_of_mass), get<1>(center_of_mass),
205 : get<2>(center_of_mass)};
206 : std::vector<std::string> legend{
207 : "NumberOfPoints", "AdmMass",
208 : "AdmLinearMomentum_x", "AdmLinearMomentum_y",
209 : "AdmLinearMomentum_z", "AdmAngularMomentum_z",
210 : "CenterOfMass_x", "CenterOfMass_y",
211 : "CenterOfMass_z"};
212 :
213 : // Get information required for reduction
214 : auto& local_observer = *Parallel::local_branch(
215 : Parallel::get_parallel_component<
216 : tmpl::conditional_t<Parallel::is_nodegroup_v<ParallelComponent>,
217 : observers::ObserverWriter<Metavariables>,
218 : observers::Observer<Metavariables>>>(cache));
219 : observers::ObservationId observation_id{observation_value.value,
220 : subfile_path + ".dat"};
221 : Parallel::ArrayComponentId array_component_id{
222 : std::add_pointer_t<ParallelComponent>{nullptr},
223 : Parallel::ArrayIndex<ElementId<3>>(array_index)};
224 :
225 : // Send reduction action
226 : if constexpr (Parallel::is_nodegroup_v<ParallelComponent>) {
227 : Parallel::threaded_action<
228 : observers::ThreadedActions::CollectReductionDataOnNode>(
229 : local_observer, std::move(observation_id),
230 : std::move(array_component_id), subfile_path, std::move(legend),
231 : std::move(reduction_data));
232 : } else {
233 : Parallel::simple_action<observers::Actions::ContributeReductionData>(
234 : local_observer, std::move(observation_id),
235 : std::move(array_component_id), subfile_path, std::move(legend),
236 : std::move(reduction_data));
237 : }
238 : }
239 :
240 0 : using observation_registration_tags = tmpl::list<::Tags::DataBox>;
241 :
242 : template <typename DbTagsList>
243 : std::optional<
244 : std::pair<observers::TypeOfObservation, observers::ObservationKey>>
245 0 : get_observation_type_and_key_for_registration(
246 : const db::DataBox<DbTagsList>& box) const {
247 : const std::optional<std::string> section_observation_key =
248 : observers::get_section_observation_key<ArraySectionIdTag>(box);
249 : if (not section_observation_key.has_value()) {
250 : return std::nullopt;
251 : }
252 : return {{observers::TypeOfObservation::Reduction,
253 : observers::ObservationKey(
254 : subfile_path_ + section_observation_key.value() + ".dat")}};
255 : }
256 :
257 0 : using is_ready_argument_tags = tmpl::list<>;
258 :
259 : template <typename Metavariables, typename ArrayIndex, typename Component>
260 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
261 : const ArrayIndex& /*array_index*/,
262 : const Component* const /*meta*/) const {
263 : return true;
264 : }
265 :
266 1 : bool needs_evolved_variables() const override { return false; }
267 :
268 : // NOLINTNEXTLINE(google-runtime-references)
269 0 : void pup(PUP::er& p) override {
270 : Event::pup(p);
271 : p | subfile_path_;
272 : }
273 :
274 : private:
275 0 : std::string subfile_path_ = "/AdmIntegrals";
276 : };
277 : /// @}
278 :
279 : /// \cond
280 : template <typename ArraySectionIdTag>
281 : PUP::able::PUP_ID ObserveAdmIntegrals<ArraySectionIdTag>::my_PUP_ID =
282 : 0; // NOLINT
283 : /// \endcond
284 :
285 : } // namespace Events
|
