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 : * See `Xcts::adm_linear_momentum_surface_integrand` for details on the formula
43 : * for each integrand.
44 : */
45 1 : void local_adm_integrals(
46 : gsl::not_null<Scalar<double>*> adm_mass,
47 : gsl::not_null<tnsr::I<double, 3>*> adm_linear_momentum,
48 : gsl::not_null<tnsr::I<double, 3>*> center_of_mass,
49 : const Scalar<DataVector>& conformal_factor,
50 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
51 : const tnsr::ii<DataVector, 3>& conformal_metric,
52 : const tnsr::II<DataVector, 3>& inv_conformal_metric,
53 : const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
54 : const tnsr::i<DataVector, 3>& conformal_christoffel_contracted,
55 : const tnsr::ii<DataVector, 3>& spatial_metric,
56 : const tnsr::II<DataVector, 3>& inv_spatial_metric,
57 : const tnsr::ii<DataVector, 3>& extrinsic_curvature,
58 : const Scalar<DataVector>& trace_extrinsic_curvature,
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 : const DirectionMap<3, tnsr::I<DataVector, 3>>&
64 : conformal_face_normal_vectors);
65 : /// @}
66 :
67 : /// @{
68 : /*!
69 : * \brief Observe ADM integrals after the XCTS solve.
70 : *
71 : * The surface integrals are taken over the outer boundary, which is defined as
72 : * the domain boundary in the upper logical zeta direction.
73 : *
74 : * Writes reduction quantities:
75 : * - ADM mass
76 : * - Linear momentum
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 : // ADM Mass
84 : Parallel::ReductionDatum<double, funcl::Plus<>>,
85 : // ADM Linear Momentum (x-component)
86 : Parallel::ReductionDatum<double, funcl::Plus<>>,
87 : // ADM Linear Momentum (y-component)
88 : Parallel::ReductionDatum<double, funcl::Plus<>>,
89 : // ADM Linear Momentum (z-component)
90 : Parallel::ReductionDatum<double, funcl::Plus<>>,
91 : // Center of Mass (x-component)
92 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
93 : std::index_sequence<0>>,
94 : // Center of Mass (y-component)
95 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
96 : std::index_sequence<0>>,
97 : // Center of Mass (z-component)
98 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
99 : std::index_sequence<0>>>;
100 :
101 : public:
102 : /// \cond
103 : explicit ObserveAdmIntegrals(CkMigrateMessage* msg) : Event(msg) {}
104 : using PUP::able::register_constructor;
105 : WRAPPED_PUPable_decl_template(ObserveAdmIntegrals); // NOLINT
106 : /// \endcond
107 :
108 0 : using options = tmpl::list<>;
109 0 : static constexpr Options::String help =
110 : "Observe ADM integrals after the XCTS solve.\n"
111 : "\n"
112 : "Writes reduction quantities:\n"
113 : "- Linear momentum";
114 :
115 0 : ObserveAdmIntegrals() = default;
116 :
117 0 : using observed_reduction_data_tags =
118 : observers::make_reduction_data_tags<tmpl::list<ReductionData>>;
119 :
120 0 : using compute_tags_for_observation_box = tmpl::list<>;
121 :
122 0 : using return_tags = tmpl::list<>;
123 :
124 0 : using argument_tags = tmpl::list<
125 : Xcts::Tags::ConformalFactor<DataVector>,
126 : ::Tags::deriv<Xcts::Tags::ConformalFactor<DataVector>, tmpl::size_t<3>,
127 : Frame::Inertial>,
128 : Xcts::Tags::ConformalMetric<DataVector, 3, Frame::Inertial>,
129 : Xcts::Tags::InverseConformalMetric<DataVector, 3, Frame::Inertial>,
130 : Xcts::Tags::ConformalChristoffelSecondKind<DataVector, 3,
131 : Frame::Inertial>,
132 : Xcts::Tags::ConformalChristoffelContracted<DataVector, 3,
133 : Frame::Inertial>,
134 : gr::Tags::SpatialMetric<DataVector, 3, Frame::Inertial>,
135 : gr::Tags::InverseSpatialMetric<DataVector, 3, Frame::Inertial>,
136 : gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame::Inertial>,
137 : gr::Tags::TraceExtrinsicCurvature<DataVector>,
138 : domain::Tags::InverseJacobian<3, Frame::ElementLogical, Frame::Inertial>,
139 : domain::Tags::Mesh<3>, domain::Tags::Element<3>,
140 : domain::Tags::Faces<3, domain::Tags::FaceNormal<3>>,
141 : domain::Tags::Faces<3, domain::Tags::FaceNormalVector<3>>,
142 : ::Tags::ObservationBox>;
143 :
144 : template <typename DataBoxType, typename ComputeTagsList,
145 : typename Metavariables, typename ArrayIndex,
146 : typename ParallelComponent>
147 0 : void operator()(
148 : const Scalar<DataVector>& conformal_factor,
149 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
150 : const tnsr::ii<DataVector, 3>& conformal_metric,
151 : const tnsr::II<DataVector, 3>& inv_conformal_metric,
152 : const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
153 : const tnsr::i<DataVector, 3>& conformal_christoffel_contracted,
154 : const tnsr::ii<DataVector, 3>& spatial_metric,
155 : const tnsr::II<DataVector, 3>& inv_spatial_metric,
156 : const tnsr::ii<DataVector, 3>& extrinsic_curvature,
157 : const Scalar<DataVector>& trace_extrinsic_curvature,
158 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
159 : Frame::Inertial>& inv_jacobian,
160 : const Mesh<3>& mesh, const Element<3>& element,
161 : const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals,
162 : const DirectionMap<3, tnsr::I<DataVector, 3>>&
163 : conformal_face_normal_vectors,
164 : const ObservationBox<DataBoxType, ComputeTagsList>& box,
165 : Parallel::GlobalCache<Metavariables>& cache,
166 : const ArrayIndex& array_index, const ParallelComponent* const /*meta*/,
167 : const ObservationValue& observation_value) const {
168 : // Skip observation on elements that are not part of a section
169 : const std::optional<std::string> section_observation_key =
170 : observers::get_section_observation_key<ArraySectionIdTag>(box);
171 : if (not section_observation_key.has_value()) {
172 : return;
173 : }
174 : const std::string subfile_path = subfile_path_ + *section_observation_key;
175 :
176 : Scalar<double> adm_mass;
177 : tnsr::I<double, 3> adm_linear_momentum;
178 : tnsr::I<double, 3> center_of_mass;
179 : local_adm_integrals(
180 : make_not_null(&adm_mass), make_not_null(&adm_linear_momentum),
181 : make_not_null(¢er_of_mass), conformal_factor,
182 : deriv_conformal_factor, conformal_metric, inv_conformal_metric,
183 : conformal_christoffel_second_kind, conformal_christoffel_contracted,
184 : spatial_metric, inv_spatial_metric, extrinsic_curvature,
185 : trace_extrinsic_curvature, inv_jacobian, mesh, element,
186 : conformal_face_normals, conformal_face_normal_vectors);
187 :
188 : // Save components of linear momentum as reduction data
189 : ReductionData reduction_data{get(adm_mass),
190 : get<0>(adm_linear_momentum),
191 : get<1>(adm_linear_momentum),
192 : get<2>(adm_linear_momentum),
193 : get<0>(center_of_mass),
194 : get<1>(center_of_mass),
195 : get<2>(center_of_mass)};
196 : std::vector<std::string> legend{"AdmMass",
197 : "AdmLinearMomentum_x",
198 : "AdmLinearMomentum_y",
199 : "AdmLinearMomentum_z",
200 : "CenterOfMass_x",
201 : "CenterOfMass_y",
202 : "CenterOfMass_z"};
203 :
204 : // Get information required for reduction
205 : auto& local_observer = *Parallel::local_branch(
206 : Parallel::get_parallel_component<
207 : tmpl::conditional_t<Parallel::is_nodegroup_v<ParallelComponent>,
208 : observers::ObserverWriter<Metavariables>,
209 : observers::Observer<Metavariables>>>(cache));
210 : observers::ObservationId observation_id{observation_value.value,
211 : subfile_path + ".dat"};
212 : Parallel::ArrayComponentId array_component_id{
213 : std::add_pointer_t<ParallelComponent>{nullptr},
214 : Parallel::ArrayIndex<ElementId<3>>(array_index)};
215 :
216 : // Send reduction action
217 : if constexpr (Parallel::is_nodegroup_v<ParallelComponent>) {
218 : Parallel::threaded_action<
219 : observers::ThreadedActions::CollectReductionDataOnNode>(
220 : local_observer, std::move(observation_id),
221 : std::move(array_component_id), subfile_path, std::move(legend),
222 : std::move(reduction_data));
223 : } else {
224 : Parallel::simple_action<observers::Actions::ContributeReductionData>(
225 : local_observer, std::move(observation_id),
226 : std::move(array_component_id), subfile_path, std::move(legend),
227 : std::move(reduction_data));
228 : }
229 : }
230 :
231 0 : using observation_registration_tags = tmpl::list<::Tags::DataBox>;
232 :
233 : template <typename DbTagsList>
234 : std::optional<
235 : std::pair<observers::TypeOfObservation, observers::ObservationKey>>
236 0 : get_observation_type_and_key_for_registration(
237 : const db::DataBox<DbTagsList>& box) const {
238 : const std::optional<std::string> section_observation_key =
239 : observers::get_section_observation_key<ArraySectionIdTag>(box);
240 : if (not section_observation_key.has_value()) {
241 : return std::nullopt;
242 : }
243 : return {{observers::TypeOfObservation::Reduction,
244 : observers::ObservationKey(
245 : subfile_path_ + section_observation_key.value() + ".dat")}};
246 : }
247 :
248 0 : using is_ready_argument_tags = tmpl::list<>;
249 :
250 : template <typename Metavariables, typename ArrayIndex, typename Component>
251 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
252 : const ArrayIndex& /*array_index*/,
253 : const Component* const /*meta*/) const {
254 : return true;
255 : }
256 :
257 1 : bool needs_evolved_variables() const override { return false; }
258 :
259 : // NOLINTNEXTLINE(google-runtime-references)
260 0 : void pup(PUP::er& p) override {
261 : Event::pup(p);
262 : p | subfile_path_;
263 : }
264 :
265 : private:
266 0 : std::string subfile_path_ = "/AdmIntegrals";
267 : };
268 : /// @}
269 :
270 : /// \cond
271 : template <typename ArraySectionIdTag>
272 : PUP::able::PUP_ID ObserveAdmIntegrals<ArraySectionIdTag>::my_PUP_ID =
273 : 0; // NOLINT
274 : /// \endcond
275 :
276 : } // namespace Events
|
