6 #include <boost/iterator/zip_iterator.hpp>
9 #include "DataStructures/ComplexDataVector.hpp"
10 #include "DataStructures/DataVector.hpp"
13 #include "Evolution/Systems/Cce/BoundaryData.hpp"
14 #include "Evolution/Systems/Cce/WorldtubeBufferUpdater.hpp"
15 #include "Helpers/Evolution/Systems/Cce/WriteToWorldtubeH5.hpp"
16 #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
17 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
18 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/Phi.hpp"
19 #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/Pi.hpp"
20 #include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp"
25 namespace TestHelpers {
26 template <
typename... Structure>
28 const Tensor<DataVector, Structure...>& nodal_data,
size_t l_max) noexcept {
33 for (
size_t i = 0; i < nodal_data.size(); ++i) {
35 goldberg_modal_data[i] =
40 return goldberg_modal_data;
43 template <
typename... Structure>
45 const Tensor<DataVector, Structure...>& nodal_data,
52 for (
size_t i = 0; i < nodal_data.size(); ++i) {
54 libsharp_modal_data[i] =
57 return libsharp_modal_data;
60 template <
typename AnalyticSolution>
61 void create_fake_time_varying_gh_nodal_data(
65 const AnalyticSolution& solution,
const double extraction_radius,
66 const double amplitude,
const double frequency,
const double time,
67 const size_t l_max) noexcept {
68 const size_t number_of_angular_points =
74 Spectral::Swsh::ComplexRepresentation::Interleaved>(l_max);
75 for (
const auto collocation_point : collocation) {
77 extraction_radius * (1.0 + amplitude *
sin(frequency * time)) *
78 sin(collocation_point.theta) *
cos(collocation_point.phi);
80 extraction_radius * (1.0 + amplitude *
sin(frequency * time)) *
81 sin(collocation_point.theta) *
sin(collocation_point.phi);
83 extraction_radius * (1.0 + amplitude *
sin(frequency * time)) *
84 cos(collocation_point.theta);
87 const auto kerr_schild_variables = solution.variables(
91 get<gr::Tags::Lapse<DataVector>>(kerr_schild_variables);
93 get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(kerr_schild_variables);
94 const auto& d_lapse = get<gr::Solutions::KerrSchild::DerivLapse<DataVector>>(
95 kerr_schild_variables);
97 const auto&
shift = get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(
98 kerr_schild_variables);
99 const auto& dt_shift =
100 get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
101 kerr_schild_variables);
102 const auto& d_shift = get<gr::Solutions::KerrSchild::DerivShift<DataVector>>(
103 kerr_schild_variables);
106 get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
107 kerr_schild_variables);
108 const auto& dt_spatial_metric =
get<
110 kerr_schild_variables);
111 const auto& d_spatial_metric =
112 get<gr::Solutions::KerrSchild::DerivSpatialMetric<DataVector>>(
113 kerr_schild_variables);
119 dt_spatial_metric, *
phi);
122 template <
typename AnalyticSolution>
123 void create_fake_time_varying_modal_data(
125 spatial_metric_coefficients,
127 dt_spatial_metric_coefficients,
129 dr_spatial_metric_coefficients,
130 const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> shift_coefficients,
131 const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> dt_shift_coefficients,
132 const gsl::not_null<tnsr::I<ComplexModalVector, 3>*> dr_shift_coefficients,
136 const AnalyticSolution& solution,
const double extraction_radius,
137 const double amplitude,
const double frequency,
const double time,
138 const size_t l_max,
const bool convert_to_goldberg =
true,
139 const bool apply_normalization_bug =
false) noexcept {
140 const size_t number_of_angular_points =
146 Spectral::Swsh::ComplexRepresentation::Interleaved>(l_max);
147 for (
const auto collocation_point : collocation) {
149 extraction_radius * (1.0 + amplitude *
sin(frequency * time)) *
150 sin(collocation_point.theta) *
cos(collocation_point.phi);
152 extraction_radius * (1.0 + amplitude *
sin(frequency * time)) *
153 sin(collocation_point.theta) *
sin(collocation_point.phi);
155 extraction_radius * (1.0 + amplitude *
sin(frequency * time)) *
156 cos(collocation_point.theta);
159 const auto kerr_schild_variables = solution.variables(
163 get<gr::Tags::Lapse<DataVector>>(kerr_schild_variables);
165 get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(kerr_schild_variables);
166 const auto& d_lapse = get<gr::Solutions::KerrSchild::DerivLapse<DataVector>>(
167 kerr_schild_variables);
169 const auto&
shift = get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(
170 kerr_schild_variables);
171 const auto& dt_shift =
172 get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
173 kerr_schild_variables);
174 const auto& d_shift = get<gr::Solutions::KerrSchild::DerivShift<DataVector>>(
175 kerr_schild_variables);
178 get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
179 kerr_schild_variables);
180 const auto& dt_spatial_metric =
get<
182 kerr_schild_variables);
183 const auto& d_spatial_metric =
184 get<gr::Solutions::KerrSchild::DerivSpatialMetric<DataVector>>(
185 kerr_schild_variables);
187 DataVector normalization_factor{number_of_angular_points, 1.0};
188 if (apply_normalization_bug) {
189 normalization_factor = 0.0;
190 const auto inverse_spatial_metric =
192 for (
size_t i = 0; i < 3; ++i) {
193 for (
size_t j = 0; j < 3; ++j) {
194 normalization_factor +=
197 square(extraction_radius *
198 (1.0 + amplitude *
sin(frequency * time)));
201 normalization_factor =
sqrt(normalization_factor);
208 (extraction_radius * normalization_factor);
209 tnsr::I<DataVector, 3> dr_shift{number_of_angular_points};
210 for (
size_t i = 0; i < 3; ++i) {
214 (extraction_radius * normalization_factor);
216 tnsr::ii<DataVector, 3> dr_spatial_metric{number_of_angular_points};
217 for (
size_t i = 0; i < 3; ++i) {
218 for (
size_t j = i; j < 3; ++j) {
219 dr_spatial_metric.get(i, j) =
223 (extraction_radius * normalization_factor);
227 if (convert_to_goldberg) {
228 *lapse_coefficients =
229 TestHelpers::tensor_to_goldberg_coefficients(
lapse, l_max);
230 *dt_lapse_coefficients =
231 TestHelpers::tensor_to_goldberg_coefficients(dt_lapse, l_max);
232 *dr_lapse_coefficients =
233 TestHelpers::tensor_to_goldberg_coefficients(dr_lapse, l_max);
235 *shift_coefficients =
236 TestHelpers::tensor_to_goldberg_coefficients(
shift, l_max);
237 *dt_shift_coefficients =
238 TestHelpers::tensor_to_goldberg_coefficients(dt_shift, l_max);
239 *dr_shift_coefficients =
240 TestHelpers::tensor_to_goldberg_coefficients(dr_shift, l_max);
242 *spatial_metric_coefficients =
243 TestHelpers::tensor_to_goldberg_coefficients(
spatial_metric, l_max);
244 *dt_spatial_metric_coefficients =
245 TestHelpers::tensor_to_goldberg_coefficients(dt_spatial_metric, l_max);
246 *dr_spatial_metric_coefficients =
247 TestHelpers::tensor_to_goldberg_coefficients(dr_spatial_metric, l_max);
249 *lapse_coefficients =
250 TestHelpers::tensor_to_libsharp_coefficients(
lapse, l_max);
251 *dt_lapse_coefficients =
252 TestHelpers::tensor_to_libsharp_coefficients(dt_lapse, l_max);
253 *dr_lapse_coefficients =
254 TestHelpers::tensor_to_libsharp_coefficients(dr_lapse, l_max);
256 *shift_coefficients =
257 TestHelpers::tensor_to_libsharp_coefficients(
shift, l_max);
258 *dt_shift_coefficients =
259 TestHelpers::tensor_to_libsharp_coefficients(dt_shift, l_max);
260 *dr_shift_coefficients =
261 TestHelpers::tensor_to_libsharp_coefficients(dr_shift, l_max);
263 *spatial_metric_coefficients =
264 TestHelpers::tensor_to_libsharp_coefficients(
spatial_metric, l_max);
265 *dt_spatial_metric_coefficients =
266 TestHelpers::tensor_to_libsharp_coefficients(dt_spatial_metric, l_max);
267 *dr_spatial_metric_coefficients =
268 TestHelpers::tensor_to_libsharp_coefficients(dr_spatial_metric, l_max);
272 template <
typename AnalyticSolution>
273 void write_test_file(
const AnalyticSolution& solution,
274 const std::string& filename,
const double target_time,
275 const double extraction_radius,
const double frequency,
276 const double amplitude,
const size_t l_max) noexcept {
277 const size_t goldberg_size =
square(l_max + 1);
278 tnsr::ii<ComplexModalVector, 3> spatial_metric_coefficients{goldberg_size};
279 tnsr::ii<ComplexModalVector, 3> dt_spatial_metric_coefficients{goldberg_size};
280 tnsr::ii<ComplexModalVector, 3> dr_spatial_metric_coefficients{goldberg_size};
281 tnsr::I<ComplexModalVector, 3> shift_coefficients{goldberg_size};
282 tnsr::I<ComplexModalVector, 3> dt_shift_coefficients{goldberg_size};
283 tnsr::I<ComplexModalVector, 3> dr_shift_coefficients{goldberg_size};
294 TestHelpers::WorldtubeModeRecorder recorder{filename, l_max};
295 for (
size_t t = 0; t < 30; ++t) {
296 const double time = 0.1 * t + target_time - 1.5;
297 TestHelpers::create_fake_time_varying_modal_data(
306 make_not_null(&dr_lapse_coefficients), solution, extraction_radius,
307 amplitude, frequency, time, l_max);
308 for (
size_t i = 0; i < 3; ++i) {
309 for (
size_t j = i; j < 3; ++j) {
310 recorder.append_worldtube_mode_data(
311 detail::dataset_name_for_component(
"/g", i, j), time,
312 spatial_metric_coefficients.get(i, j), l_max);
313 recorder.append_worldtube_mode_data(
314 detail::dataset_name_for_component(
"/Drg", i, j), time,
315 dr_spatial_metric_coefficients.get(i, j), l_max);
316 recorder.append_worldtube_mode_data(
317 detail::dataset_name_for_component(
"/Dtg", i, j), time,
318 dt_spatial_metric_coefficients.get(i, j), l_max);
320 recorder.append_worldtube_mode_data(
321 detail::dataset_name_for_component(
"/Shift", i), time,
322 shift_coefficients.get(i), l_max);
323 recorder.append_worldtube_mode_data(
324 detail::dataset_name_for_component(
"/DrShift", i), time,
325 dr_shift_coefficients.get(i), l_max);
326 recorder.append_worldtube_mode_data(
327 detail::dataset_name_for_component(
"/DtShift", i), time,
328 dt_shift_coefficients.get(i), l_max);
330 recorder.append_worldtube_mode_data(
331 detail::dataset_name_for_component(
"/Lapse"), time,
332 get(lapse_coefficients), l_max);
333 recorder.append_worldtube_mode_data(
334 detail::dataset_name_for_component(
"/DrLapse"), time,
335 get(dr_lapse_coefficients), l_max);
336 recorder.append_worldtube_mode_data(
337 detail::dataset_name_for_component(
"/DtLapse"), time,
338 get(dt_lapse_coefficients), l_max);