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 <optional>
8 : #include <pup.h>
9 : #include <string>
10 : #include <unordered_map>
11 : #include <utility>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/DataBox/ObservationBox.hpp"
16 : #include "DataStructures/DataBox/TagName.hpp"
17 : #include "DataStructures/DataVector.hpp"
18 : #include "DataStructures/Tensor/Tensor.hpp"
19 : #include "Domain/Creators/Tags/Domain.hpp"
20 : #include "Domain/Structure/ElementId.hpp"
21 : #include "Domain/Tags.hpp"
22 : #include "Elliptic/Systems/SelfForce/Scalar/AnalyticData/CircularOrbit.hpp"
23 : #include "Elliptic/Systems/SelfForce/Scalar/Tags.hpp"
24 : #include "Elliptic/Tags.hpp"
25 : #include "IO/Observer/GetSectionObservationKey.hpp"
26 : #include "IO/Observer/Helpers.hpp"
27 : #include "IO/Observer/ObservationId.hpp"
28 : #include "IO/Observer/ObserverComponent.hpp"
29 : #include "IO/Observer/ReductionActions.hpp"
30 : #include "IO/Observer/TypeOfObservation.hpp"
31 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
32 : #include "Options/String.hpp"
33 : #include "Parallel/ArrayIndex.hpp"
34 : #include "Parallel/GlobalCache.hpp"
35 : #include "Parallel/Invoke.hpp"
36 : #include "Parallel/Local.hpp"
37 : #include "Parallel/Reduction.hpp"
38 : #include "Parallel/TypeTraits.hpp"
39 : #include "ParallelAlgorithms/Events/Tags.hpp"
40 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
41 : #include "Utilities/ErrorHandling/Assert.hpp"
42 : #include "Utilities/ErrorHandling/Error.hpp"
43 : #include "Utilities/Functional.hpp"
44 : #include "Utilities/OptionalHelpers.hpp"
45 : #include "Utilities/PrettyType.hpp"
46 : #include "Utilities/Serialization/CharmPupable.hpp"
47 : #include "Utilities/TMPL.hpp"
48 :
49 0 : namespace ScalarSelfForce::Events {
50 :
51 : namespace detail {
52 : std::optional<tnsr::i<std::complex<double>, 2>> extract_self_force(
53 : const Domain<2>& domain, const ElementId<2>& element_id,
54 : const AnalyticData::CircularOrbit& circular_orbit,
55 : const Scalar<ComplexDataVector>& field, const Mesh<2>& mesh,
56 : const InverseJacobian<DataVector, 2, Frame::ElementLogical,
57 : Frame::Inertial>& inv_jacobian);
58 : } // namespace detail
59 :
60 : /*!
61 : * \brief Observe the self-force at the position of the scalar charge.
62 : *
63 : * \details This event interpolates the scalar m-mode field and its derivatives
64 : * to the position of the scalar charge (puncture) and computes the m-mode
65 : * contribution to the self-force acting on the charge. The self-force is then
66 : * written to a reduction file. If multiple elements contain the puncture, their
67 : * contributions are averaged (they should only differ by truncation error).
68 : *
69 : * The self-force is computed in the Boyer-Lindquist `r` and `theta`
70 : * coordinates (for scalar charge $q=1$) following Eq. (3.10) in
71 : * \cite Osburn:2022bby :
72 : * \begin{align}
73 : * F_r^m &= \partial_r (\Psi_m^R / r) = \frac{1}{r \alpha}
74 : * \partial_{r_*} \Psi_m^R - \frac{1}{r^2} \Psi_m^R \\
75 : * F_\theta^m &= \partial_\theta (\Psi_m^R / r) = -\frac{1}{r}
76 : * \partial_{\cos\theta} \Psi_m^R
77 : * \end{align}
78 : * with $\alpha = 1 - 2 M r / (r^2 + a^2)$. For modes $m > 0$ the self-force
79 : * includes an additional factor of 2 and a complex rotation by
80 : * $m \Delta\phi = \frac{m a}{r_+ - r_-}\ln\frac{r - r_+}{r - r_-}$.
81 : *
82 : * These contributions can be summed over all m-modes (in postprocessing) to
83 : * form the total self-force $F_\mu^\mathrm{self} = \sum_m F_\mu^m$, where the
84 : * sum is truncated at some maximum m-mode number.
85 : */
86 : template <typename ArraySectionIdTag = void>
87 1 : class ObserveSelfForce : public Event {
88 : public:
89 0 : using ReductionData = Parallel::ReductionData<
90 : // Observation value
91 : Parallel::ReductionDatum<double, funcl::AssertEqual<>>,
92 : // Number of grid points
93 : Parallel::ReductionDatum<size_t, funcl::Plus<>>,
94 : // Num contributing elements
95 : Parallel::ReductionDatum<size_t, funcl::Plus<>>,
96 : // Re(SelfForce_r)
97 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
98 : std::index_sequence<2>>,
99 : // Im(SelfForce_r)
100 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
101 : std::index_sequence<2>>,
102 : // Re(SelfForce_theta)
103 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
104 : std::index_sequence<2>>,
105 : // Im(SelfForce_theta)
106 : Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
107 : std::index_sequence<2>>>;
108 :
109 0 : static constexpr Options::String help =
110 : "Observe the self-force at the position of the scalar charge.";
111 0 : using options = tmpl::list<>;
112 :
113 0 : ObserveSelfForce() = default;
114 :
115 0 : explicit ObserveSelfForce(CkMigrateMessage* msg) : Event(msg) {}
116 : using PUP::able::register_constructor;
117 0 : WRAPPED_PUPable_decl_template(ObserveSelfForce); // NOLINT
118 :
119 0 : using compute_tags_for_observation_box = tmpl::list<>;
120 0 : using observed_reduction_data_tags =
121 : observers::make_reduction_data_tags<tmpl::list<ReductionData>>;
122 :
123 0 : using return_tags = tmpl::list<>;
124 0 : using argument_tags = tmpl::list<::Tags::ObservationBox>;
125 :
126 0 : using BackgroundTag =
127 : elliptic::Tags::Background<elliptic::analytic_data::Background>;
128 :
129 : template <typename ComputeTagsList, typename DataBoxType,
130 : typename Metavariables, typename ParallelComponent>
131 0 : void operator()(const ObservationBox<ComputeTagsList, DataBoxType>& box,
132 : Parallel::GlobalCache<Metavariables>& cache,
133 : const ElementId<2>& element_id,
134 : const ParallelComponent* const /*meta*/,
135 : const ObservationValue& observation_value) const {
136 : // Skip observation on elements that are not part of the section
137 : const std::optional<std::string> section_observation_key =
138 : observers::get_section_observation_key<ArraySectionIdTag>(box);
139 : if (not section_observation_key.has_value()) {
140 : return;
141 : }
142 : // Extract self-force
143 : const auto& background = get<BackgroundTag>(box);
144 : const auto& circular_orbit =
145 : dynamic_cast<const AnalyticData::CircularOrbit&>(background);
146 : const auto& mesh = get<domain::Tags::Mesh<2>>(box);
147 : const size_t num_points = mesh.number_of_grid_points();
148 : const auto self_force = detail::extract_self_force(
149 : get<domain::Tags::Domain<2>>(box), element_id, circular_orbit,
150 : get<Tags::MMode>(box), mesh,
151 : get<domain::Tags::InverseJacobian<2, Frame::ElementLogical,
152 : Frame::Inertial>>(box));
153 : // Send data to reduction observer
154 : auto& local_observer = *Parallel::local_branch(
155 : Parallel::get_parallel_component<
156 : tmpl::conditional_t<Parallel::is_nodegroup_v<ParallelComponent>,
157 : observers::ObserverWriter<Metavariables>,
158 : observers::Observer<Metavariables>>>(cache));
159 : const std::string subfile_name = "SelfForce";
160 : observers::ObservationId observation_id{observation_value.value,
161 : subfile_name + ".dat"};
162 : Parallel::ArrayComponentId array_component_id{
163 : std::add_pointer_t<ParallelComponent>{nullptr},
164 : Parallel::ArrayIndex<ElementId<2>>(element_id)};
165 : ReductionData reduction_data =
166 : self_force.has_value()
167 : ? ReductionData{observation_value.value,
168 : num_points,
169 : 1_st,
170 : get<0>(*self_force).real(),
171 : get<0>(*self_force).imag(),
172 : get<1>(*self_force).real(),
173 : get<1>(*self_force).imag()}
174 : : ReductionData{
175 : observation_value.value, num_points, 0_st, 0., 0., 0., 0.};
176 : std::vector<std::string> legend{
177 : "IterationId", "NumGridPoints", "NumContributingElements",
178 : "Re(SelfForce_r)", "Im(SelfForce_r)", "Re(SelfForce_theta)",
179 : "Im(SelfForce_theta)"};
180 : if constexpr (Parallel::is_nodegroup_v<ParallelComponent>) {
181 : Parallel::threaded_action<
182 : observers::ThreadedActions::CollectReductionDataOnNode>(
183 : local_observer, std::move(observation_id),
184 : std::move(array_component_id), subfile_name + ".dat",
185 : std::move(legend), std::move(reduction_data));
186 : } else {
187 : Parallel::simple_action<observers::Actions::ContributeReductionData>(
188 : local_observer, std::move(observation_id),
189 : std::move(array_component_id), subfile_name + ".dat",
190 : std::move(legend), std::move(reduction_data));
191 : }
192 : }
193 :
194 0 : using observation_registration_tags = tmpl::list<::Tags::DataBox>;
195 :
196 : template <typename DbTagsList>
197 : std::optional<
198 : std::pair<observers::TypeOfObservation, observers::ObservationKey>>
199 0 : get_observation_type_and_key_for_registration(
200 : const db::DataBox<DbTagsList>& box) const {
201 : const std::optional<std::string> section_observation_key =
202 : observers::get_section_observation_key<ArraySectionIdTag>(box);
203 : if (not section_observation_key.has_value()) {
204 : return std::nullopt;
205 : }
206 : return {{observers::TypeOfObservation::Reduction,
207 : observers::ObservationKey(
208 : "SelfForce" + section_observation_key.value() + ".dat")}};
209 : }
210 :
211 0 : using is_ready_argument_tags = tmpl::list<>;
212 :
213 : template <typename Metavariables, typename ArrayIndex, typename Component>
214 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
215 : const ArrayIndex& /*array_index*/,
216 : const Component* const /*meta*/) const {
217 : return true;
218 : }
219 :
220 1 : bool needs_evolved_variables() const override { return false; }
221 : };
222 :
223 : /// \cond
224 : // NOLINTBEGIN(cppcoreguidelines-avoid-non-const-global-variables)
225 : template <typename ArraySectionIdTag>
226 : PUP::able::PUP_ID ObserveSelfForce<ArraySectionIdTag>::my_PUP_ID = 0;
227 : // NOLINTEND(cppcoreguidelines-avoid-non-const-global-variables)
228 : /// \endcond
229 : } // namespace ScalarSelfForce::Events
|
