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 <limits>
8 : #include <memory>
9 : #include <pup.h>
10 : #include <type_traits>
11 :
12 : #include "DataStructures/CachedTempBuffer.hpp"
13 : #include "DataStructures/DataBox/Prefixes.hpp"
14 : #include "DataStructures/TaggedTuple.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "Elliptic/Systems/Xcts/Tags.hpp"
17 : #include "NumericalAlgorithms/LinearOperators/Divergence.hpp"
18 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
19 : #include "Options/String.hpp"
20 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
21 : #include "PointwiseFunctions/AnalyticSolutions/Xcts/CommonVariables.hpp"
22 : #include "PointwiseFunctions/AnalyticSolutions/Xcts/Flatness.hpp"
23 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
24 : #include "PointwiseFunctions/GeneralRelativity/Tags/Conformal.hpp"
25 : #include "PointwiseFunctions/Hydro/Tags.hpp"
26 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
27 : #include "Utilities/Gsl.hpp"
28 : #include "Utilities/PrettyType.hpp"
29 : #include "Utilities/Requires.hpp"
30 : #include "Utilities/Serialization/CharmPupable.hpp"
31 : #include "Utilities/TMPL.hpp"
32 :
33 : namespace Xcts::Solutions {
34 : namespace detail {
35 :
36 : template <typename DataType, size_t Dim>
37 : using gr_solution_vars =
38 : tmpl::list<gr::Tags::SpatialMetric<DataType, Dim>,
39 : gr::Tags::InverseSpatialMetric<DataType, Dim>,
40 : ::Tags::deriv<gr::Tags::SpatialMetric<DataType, Dim>,
41 : tmpl::size_t<Dim>, Frame::Inertial>,
42 : gr::Tags::Lapse<DataType>,
43 : ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<Dim>,
44 : Frame::Inertial>,
45 : gr::Tags::Shift<DataType, Dim>,
46 : ::Tags::deriv<gr::Tags::Shift<DataType, Dim>, tmpl::size_t<Dim>,
47 : Frame::Inertial>,
48 : gr::Tags::ExtrinsicCurvature<DataType, Dim>>;
49 :
50 : template <typename DataType>
51 : using WrappedGrVariablesCache =
52 : cached_temp_buffer_from_typelist<tmpl::push_back<
53 : common_tags<DataType>,
54 : hydro::Tags::MagneticFieldDotSpatialVelocity<DataType>,
55 : hydro::Tags::ComovingMagneticFieldSquared<DataType>,
56 : gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0>,
57 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>, 0>,
58 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, 3>, 0>>>;
59 :
60 : template <typename DataType, bool HasMhd>
61 : struct WrappedGrVariables
62 : : CommonVariables<DataType, WrappedGrVariablesCache<DataType>> {
63 : static constexpr size_t Dim = 3;
64 : using Cache = WrappedGrVariablesCache<DataType>;
65 : using Base = CommonVariables<DataType, WrappedGrVariablesCache<DataType>>;
66 : using Base::operator();
67 :
68 : WrappedGrVariables(
69 : std::optional<std::reference_wrapper<const Mesh<Dim>>> local_mesh,
70 : std::optional<std::reference_wrapper<const InverseJacobian<
71 : DataType, Dim, Frame::ElementLogical, Frame::Inertial>>>
72 : local_inv_jacobian,
73 : const tnsr::I<DataType, 3>& local_x,
74 : const tuples::tagged_tuple_from_typelist<gr_solution_vars<DataType, Dim>>&
75 : local_gr_solution,
76 : const tuples::tagged_tuple_from_typelist<hydro_tags<DataType>>&
77 : local_hydro_solution)
78 : : Base(std::move(local_mesh), std::move(local_inv_jacobian)),
79 : x(local_x),
80 : gr_solution(local_gr_solution),
81 : hydro_solution(local_hydro_solution) {}
82 :
83 : const tnsr::I<DataType, Dim>& x;
84 : const tuples::tagged_tuple_from_typelist<gr_solution_vars<DataType, Dim>>&
85 : gr_solution;
86 : const tuples::tagged_tuple_from_typelist<hydro_tags<DataType>>&
87 : hydro_solution;
88 :
89 : void operator()(
90 : gsl::not_null<tnsr::ii<DataType, Dim>*> conformal_metric,
91 : gsl::not_null<Cache*> cache,
92 : Xcts::Tags::ConformalMetric<DataType, Dim, Frame::Inertial> /*meta*/)
93 : const override;
94 : void operator()(gsl::not_null<tnsr::II<DataType, Dim>*> inv_conformal_metric,
95 : gsl::not_null<Cache*> cache,
96 : Xcts::Tags::InverseConformalMetric<
97 : DataType, Dim, Frame::Inertial> /*meta*/) const override;
98 : void operator()(
99 : gsl::not_null<tnsr::ijj<DataType, Dim>*> deriv_conformal_metric,
100 : gsl::not_null<Cache*> cache,
101 : ::Tags::deriv<Xcts::Tags::ConformalMetric<DataType, Dim, Frame::Inertial>,
102 : tmpl::size_t<Dim>, Frame::Inertial> /*meta*/)
103 : const override;
104 : void operator()(
105 : gsl::not_null<tnsr::ii<DataType, Dim>*> spatial_metric,
106 : gsl::not_null<Cache*> cache,
107 : gr::Tags::SpatialMetric<DataType, Dim> /*meta*/) const override;
108 : void operator()(
109 : gsl::not_null<tnsr::II<DataType, Dim>*> inv_spatial_metric,
110 : gsl::not_null<Cache*> cache,
111 : gr::Tags::InverseSpatialMetric<DataType, Dim> /*meta*/) const override;
112 : void operator()(
113 : gsl::not_null<tnsr::ijj<DataType, Dim>*> deriv_spatial_metric,
114 : gsl::not_null<Cache*> cache,
115 : ::Tags::deriv<gr::Tags::SpatialMetric<DataType, Dim>, tmpl::size_t<Dim>,
116 : Frame::Inertial> /*meta*/) const override;
117 : void operator()(
118 : gsl::not_null<Scalar<DataType>*> trace_extrinsic_curvature,
119 : gsl::not_null<Cache*> cache,
120 : gr::Tags::TraceExtrinsicCurvature<DataType> /*meta*/) const override;
121 : void operator()(
122 : gsl::not_null<Scalar<DataType>*> dt_trace_extrinsic_curvature,
123 : gsl::not_null<Cache*> cache,
124 : ::Tags::dt<gr::Tags::TraceExtrinsicCurvature<DataType>> /*meta*/)
125 : const override;
126 : void operator()(
127 : gsl::not_null<Scalar<DataType>*> conformal_factor,
128 : gsl::not_null<Cache*> cache,
129 : Xcts::Tags::ConformalFactor<DataType> /*meta*/) const override;
130 : void operator()(
131 : gsl::not_null<Scalar<DataType>*> conformal_factor_minus_one,
132 : gsl::not_null<Cache*> cache,
133 : Xcts::Tags::ConformalFactorMinusOne<DataType> /*meta*/) const override;
134 : void operator()(
135 : gsl::not_null<tnsr::i<DataType, Dim>*> conformal_factor_gradient,
136 : gsl::not_null<Cache*> cache,
137 : ::Tags::deriv<Xcts::Tags::ConformalFactorMinusOne<DataType>,
138 : tmpl::size_t<Dim>, Frame::Inertial> /*meta*/)
139 : const override;
140 : void operator()(gsl::not_null<Scalar<DataType>*> lapse,
141 : gsl::not_null<Cache*> cache,
142 : gr::Tags::Lapse<DataType> /*meta*/) const override;
143 : void operator()(
144 : gsl::not_null<Scalar<DataType>*> lapse_times_conformal_factor_minus_one,
145 : gsl::not_null<Cache*> cache,
146 : Xcts::Tags::LapseTimesConformalFactorMinusOne<DataType> /*meta*/)
147 : const override;
148 : void operator()(
149 : gsl::not_null<Scalar<DataType>*> lapse_times_conformal_factor,
150 : gsl::not_null<Cache*> cache,
151 : Xcts::Tags::LapseTimesConformalFactor<DataType> /*meta*/) const override;
152 : void operator()(
153 : gsl::not_null<tnsr::i<DataType, Dim>*>
154 : lapse_times_conformal_factor_gradient,
155 : gsl::not_null<Cache*> cache,
156 : ::Tags::deriv<Xcts::Tags::LapseTimesConformalFactorMinusOne<DataType>,
157 : tmpl::size_t<Dim>, Frame::Inertial> /*meta*/)
158 : const override;
159 : void operator()(
160 : gsl::not_null<tnsr::I<DataType, Dim>*> shift_background,
161 : gsl::not_null<Cache*> cache,
162 : Xcts::Tags::ShiftBackground<DataType, Dim, Frame::Inertial> /*meta*/)
163 : const override;
164 : void operator()(
165 : gsl::not_null<tnsr::iJ<DataType, 3, Frame::Inertial>*>
166 : deriv_shift_background,
167 : gsl::not_null<Cache*> cache,
168 : ::Tags::deriv<Xcts::Tags::ShiftBackground<DataType, 3, Frame::Inertial>,
169 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
170 : void operator()(gsl::not_null<tnsr::II<DataType, Dim, Frame::Inertial>*>
171 : longitudinal_shift_background_minus_dt_conformal_metric,
172 : gsl::not_null<Cache*> cache,
173 : Xcts::Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
174 : DataType, Dim, Frame::Inertial> /*meta*/) const override;
175 : void operator()(
176 : gsl::not_null<tnsr::I<DataType, Dim>*> shift_excess,
177 : gsl::not_null<Cache*> cache,
178 : Xcts::Tags::ShiftExcess<DataType, Dim, Frame::Inertial> /*meta*/)
179 : const override;
180 : void operator()(
181 : gsl::not_null<tnsr::iJ<DataType, Dim>*> deriv_shift_excess,
182 : gsl::not_null<Cache*> cache,
183 : ::Tags::deriv<Xcts::Tags::ShiftExcess<DataType, 3, Frame::Inertial>,
184 : tmpl::size_t<3>, Frame::Inertial> /*meta*/) const override;
185 : void operator()(
186 : gsl::not_null<tnsr::ii<DataType, 3>*> extrinsic_curvature,
187 : gsl::not_null<Cache*> cache,
188 : gr::Tags::ExtrinsicCurvature<DataType, 3> /*meta*/) const override;
189 : void operator()(
190 : gsl::not_null<Scalar<DataType>*> magnetic_field_dot_spatial_velocity,
191 : gsl::not_null<Cache*> cache,
192 : hydro::Tags::MagneticFieldDotSpatialVelocity<DataType> /*meta*/) const;
193 : void operator()(
194 : gsl::not_null<Scalar<DataType>*> comoving_magnetic_field_squared,
195 : gsl::not_null<Cache*> cache,
196 : hydro::Tags::ComovingMagneticFieldSquared<DataType> /*meta*/) const;
197 : void operator()(
198 : gsl::not_null<Scalar<DataType>*> energy_density,
199 : gsl::not_null<Cache*> cache,
200 : gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0> /*meta*/) const;
201 : void operator()(
202 : gsl::not_null<Scalar<DataType>*> stress_trace,
203 : gsl::not_null<Cache*> cache,
204 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>, 0> /*meta*/) const;
205 : void operator()(gsl::not_null<tnsr::I<DataType, Dim>*> momentum_density,
206 : gsl::not_null<Cache*> cache,
207 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, Dim>,
208 : 0> /*meta*/) const;
209 : };
210 :
211 : } // namespace detail
212 :
213 : /*!
214 : * \brief XCTS quantities for a solution of the Einstein equations
215 : *
216 : * This class computes all XCTS quantities from the `GrSolution`. To do so, it
217 : * chooses the conformal factor
218 : *
219 : * \f{equation}{
220 : * \psi = 1
221 : * \text{,}
222 : * \f}
223 : *
224 : * so the spatial metric of the `GrSolution` is used as conformal metric,
225 : * \f$\bar{\gamma}_{ij = \gamma_{ij}\f$. This is particularly useful for
226 : * superpositions, because it means that the superposed conformal metric of two
227 : * `WrappedGr` solutions is probably a good conformal background to solve for a
228 : * binary solution (see Xcts::AnalyticData::Binary).
229 : *
230 : * For example, when the `GrSolution` is `gr::Solutions::KerrSchild`, the
231 : * conformal metric is the spatial Kerr metric in Kerr-Schild coordinates and
232 : * \f$\psi = 1\f$. It is also possible to
233 : * choose a different \f$\psi\f$ so the solution is non-trivial in this
234 : * variable, though that is probably only useful for testing and currently not
235 : * implemented. It should be noted, however, that the combination of
236 : * \f$\psi=1\f$ and apparent-horizon boundary conditions poses a hard problem to
237 : * the nonlinear solver when starting at a flat initial guess. This is because
238 : * the strongly-nonlinear boundary-conditions couple the variables in such a way
239 : * that the solution is initially corrected away from \f$\psi=1\f$ and is then
240 : * unable to recover. A conformal-factor profile such as \f$\psi=1 +
241 : * \frac{M}{2r}\f$ (resembling isotropic coordinates) resolves this issue. In
242 : * production solves this is not an issue because we choose a much better
243 : * initial guess than flatness, such as a superposition of Kerr solutions for
244 : * black-hole binary initial data.
245 : *
246 : * \warning
247 : * The computation of the XCTS matter source terms (energy density $\rho$,
248 : * momentum density $S^i$, stress trace $S$) uses GR quantities (lapse $\alpha$,
249 : * shift $\beta^i$, spatial metric $\gamma_{ij}$), which means these GR
250 : * quantities are not treated dynamically in the source terms when solving the
251 : * XCTS equations. If the GR quantities satisfy the Einstein constraints (as is
252 : * the case if the `GrSolution` is actually a solution to the Einstein
253 : * equations), then the XCTS solve will reproduce the GR quantities given the
254 : * fixed sources computed here. However, if the GR quantities don't satisfy the
255 : * Einstein constraints (e.g. because a magnetic field was added to the solution
256 : * but ignored in the gravity sector, or because it is a hydrodynamic solution
257 : * on a fixed background metric) then the XCTS solution will depend on our
258 : * treatment of the source terms: fixing the source terms (the simple approach
259 : * taken here) means we're making a choice of $W$ and $u^i$. This is what
260 : * initial data codes usually do when they iterate back and forth between a
261 : * hydro solve and an XCTS solve (e.g. see \cite Tacik2016zal). Alternatively,
262 : * we could fix $v^i$ and compute $W$ and $u^i$ from $v^i$ and the dynamic
263 : * metric variables at every step in the XCTS solver algorithm. This requires
264 : * adding the source terms and their linearization to the XCTS equations, and
265 : * could be interesting to explore.
266 : *
267 : * \tparam GrSolution Any solution to the Einstein constraint equations
268 : * \tparam HasMhd Enable to compute matter source terms. Disable to set matter
269 : * source terms to zero.
270 : */
271 :
272 : template <typename GrSolution, bool HasMhd = false>
273 : class WrappedGr;
274 :
275 : template <typename GrSolution, bool HasMhd>
276 1 : class WrappedGr : public elliptic::analytic_data::AnalyticSolution {
277 : public:
278 0 : static constexpr size_t Dim = 3;
279 :
280 0 : using options = typename GrSolution::options;
281 0 : static constexpr Options::String help = GrSolution::help;
282 0 : static std::string name() { return pretty_type::name<GrSolution>(); }
283 :
284 0 : WrappedGr() = default;
285 0 : WrappedGr(const WrappedGr&) = default;
286 0 : WrappedGr& operator=(const WrappedGr&) = default;
287 0 : WrappedGr(WrappedGr&&) = default;
288 0 : WrappedGr& operator=(WrappedGr&&) = default;
289 0 : ~WrappedGr() = default;
290 :
291 : template <typename... Args,
292 : Requires<std::is_constructible_v<GrSolution, Args...>> = nullptr>
293 0 : explicit WrappedGr(Args&&... args)
294 : : gr_solution_(std::forward<Args>(args)...) {}
295 :
296 0 : const GrSolution& gr_solution() const { return gr_solution_; }
297 :
298 : /// \cond
299 : explicit WrappedGr(CkMigrateMessage* m)
300 : : elliptic::analytic_data::AnalyticSolution(m) {}
301 : using PUP::able::register_constructor;
302 : WRAPPED_PUPable_decl_template(WrappedGr);
303 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
304 : const override {
305 : return std::make_unique<WrappedGr>(*this);
306 : }
307 : /// \endcond
308 :
309 : template <typename DataType, typename... RequestedTags>
310 0 : tuples::TaggedTuple<RequestedTags...> variables(
311 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
312 : tmpl::list<RequestedTags...> /*meta*/) const {
313 : return variables_impl<DataType>(x, std::nullopt, std::nullopt,
314 : tmpl::list<RequestedTags...>{});
315 : }
316 :
317 : template <typename DataType, typename... RequestedTags>
318 0 : tuples::TaggedTuple<RequestedTags...> variables(
319 : const tnsr::I<DataType, 3, Frame::Inertial>& x, const Mesh<3>& mesh,
320 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
321 : Frame::Inertial>& inv_jacobian,
322 : tmpl::list<RequestedTags...> /*meta*/) const {
323 : return variables_impl<DataVector>(x, mesh, inv_jacobian,
324 : tmpl::list<RequestedTags...>{});
325 : }
326 :
327 0 : void pup(PUP::er& p) override {
328 : elliptic::analytic_data::AnalyticSolution::pup(p);
329 : p | gr_solution_;
330 : }
331 :
332 : private:
333 : template <typename DataType, typename... RequestedTags>
334 0 : tuples::TaggedTuple<RequestedTags...> variables_impl(
335 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
336 : std::optional<std::reference_wrapper<const Mesh<3>>> mesh,
337 : std::optional<std::reference_wrapper<const InverseJacobian<
338 : DataType, 3, Frame::ElementLogical, Frame::Inertial>>>
339 : inv_jacobian,
340 : tmpl::list<RequestedTags...> /*meta*/) const {
341 : tuples::tagged_tuple_from_typelist<detail::gr_solution_vars<DataType, Dim>>
342 : gr_solution;
343 : if constexpr (is_analytic_solution_v<GrSolution>) {
344 : gr_solution = gr_solution_.variables(
345 : x, std::numeric_limits<double>::signaling_NaN(),
346 : detail::gr_solution_vars<DataType, Dim>{});
347 : } else {
348 : gr_solution =
349 : gr_solution_.variables(x, detail::gr_solution_vars<DataType, Dim>{});
350 : }
351 : tuples::tagged_tuple_from_typelist<hydro_tags<DataType>> hydro_solution;
352 : if constexpr (HasMhd) {
353 : if constexpr (is_analytic_solution_v<GrSolution>) {
354 : hydro_solution = gr_solution_.variables(
355 : x, std::numeric_limits<double>::signaling_NaN(),
356 : hydro_tags<DataType>{});
357 : } else {
358 : hydro_solution = gr_solution_.variables(x, hydro_tags<DataType>{});
359 : }
360 : }
361 : using VarsComputer = detail::WrappedGrVariables<DataType, HasMhd>;
362 : const size_t num_points = get_size(*x.begin());
363 : typename VarsComputer::Cache cache{num_points};
364 : VarsComputer computer{mesh, inv_jacobian, x, gr_solution, hydro_solution};
365 : const auto get_var = [&cache, &computer, &hydro_solution, &x](auto tag_v) {
366 : using tag = std::decay_t<decltype(tag_v)>;
367 : if constexpr (tmpl::list_contains_v<hydro_tags<DataType>, tag>) {
368 : (void)cache;
369 : (void)computer;
370 : if constexpr (HasMhd) {
371 : (void)x;
372 : return get<tag>(hydro_solution);
373 : } else {
374 : (void)hydro_solution;
375 : return get<tag>(Flatness{}.variables(x, tmpl::list<tag>{}));
376 : }
377 : } else {
378 : (void)hydro_solution;
379 : (void)x;
380 : return cache.get_var(computer, tag{});
381 : }
382 : };
383 : return {get_var(RequestedTags{})...};
384 : }
385 :
386 0 : friend bool operator==(const WrappedGr<GrSolution, HasMhd>& lhs,
387 : const WrappedGr<GrSolution, HasMhd>& rhs) {
388 : return lhs.gr_solution_ == rhs.gr_solution_;
389 : }
390 :
391 0 : GrSolution gr_solution_;
392 : };
393 :
394 : template <typename GrSolution, bool HasMhd>
395 0 : inline bool operator!=(const WrappedGr<GrSolution, HasMhd>& lhs,
396 : const WrappedGr<GrSolution, HasMhd>& rhs) {
397 : return not(lhs == rhs);
398 : }
399 :
400 : template <typename GrMhdSolution>
401 0 : using WrappedGrMhd = WrappedGr<GrMhdSolution, true>;
402 :
403 : } // namespace Xcts::Solutions
|
