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