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 :
9 : #include "DataStructures/DataBox/Tag.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Domain/Tags.hpp"
12 : #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Gauges.hpp"
13 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
14 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
15 : #include "Utilities/ContainerHelpers.hpp"
16 : #include "Utilities/TMPL.hpp"
17 :
18 : /// \cond
19 : class DataVector;
20 : namespace gsl {
21 : template <class T>
22 : class not_null;
23 : } // namespace gsl
24 : /// \endcond
25 :
26 : namespace gh {
27 : namespace gauges {
28 : /*!
29 : * \brief Damped harmonic gauge source function and its spacetime derivative.
30 : *
31 : * \details The gauge condition has been taken from \cite Szilagyi2009qz and
32 : * \cite Deppe2018uye. We provide both a "rollon" version
33 : * (`damped_harmonic_rollon`), and a "non-rollon" version (`damped_harmonic`).
34 : * In the non-rollon version the rollon function \f$R(t)=1\f$.
35 : *
36 : * \warning Only the non-rollon version can be used with a moving mesh.
37 : *
38 : * The covariant form of the source function \f$H_a\f$ is written as:
39 : *
40 : * \f{align*}
41 : * H_a := [1 - R(t)] H_a^\mathrm{init} +
42 : * [\mu_{L1} \mathrm{log}(\sqrt{\gamma}/\alpha) + \mu_{L2}
43 : * \mathrm{log}(1/\alpha)] n_a
44 : * - \mu_S \gamma_{ai} \beta^i / \alpha
45 : * \f}
46 : *
47 : * where \f$\alpha, \beta^k\f$ are the lapse and shift respectively, \f$n_a\f$
48 : * is the unit normal one-form to the spatial slice, and \f$\gamma_{ab}\f$ is
49 : * the spatial metric (obtained by projecting the spacetime metric onto the
50 : * 3-slice, i.e. \f$\gamma_{ab} = g_{ab} + n_a n_b\f$). The prefactors are:
51 : *
52 : * \f{align*}
53 : * \mu_{L1} &= A_{L1} R(t) W(x^i) \mathrm{log}(\sqrt{\gamma}/\alpha)^{e_{L1}},
54 : * \\
55 : * \mu_{L2} &= A_{L2} R(t) W(x^i) \mathrm{log}(1/\alpha)^{e_{L2}}, \\
56 : * \mu_{S} &= A_{S} R(t) W(x^i) \mathrm{log}(\sqrt{\gamma}/\alpha)^{e_{S}},
57 : * \f}
58 : *
59 : * temporal roll-on function \f$ R(t)\f$ is:
60 : *
61 : * \f{align*}
62 : * \begin{array}{ll}
63 : * R(t) & = 0, & t< t_{0} \\
64 : * & = 1 - \exp[-((t - t_{0})/ \sigma_t)^4], & t\geq t_{0} \\
65 : * \end{array}
66 : * \f}
67 : *
68 : * and the spatial weight function is:
69 : *
70 : * \f{align*}
71 : * W(x^i) = \exp[-(r/\sigma_r)^2].
72 : * \f}
73 : *
74 : * This weight function can be written with multiple constant factors in the
75 : * exponent in literature \cite Deppe2018uye, but we absorb them all into
76 : * \f$ \sigma_r\f$ here. The coordinate \f$ r\f$ is the Euclidean radius
77 : * in Inertial coordinates.
78 : *
79 : * Note that for the last three terms in \f$H_a\f$ (with \f$ X = \{L1, L2, S\}
80 : * \f$):
81 : * - Amplitude factors \f$ A_{X} \f$ are taken as input here as `amp_coef_X`
82 : * - Exponents \f$ e_X\f$ are taken as input here as `exp_X`.
83 : * - Spatial weight function \f$W\f$ is specified completely by
84 : * \f$\sigma_r\f$, which is taken as input here as `sigma_r`.
85 : *
86 : * Also computes spacetime derivatives, i.e. \f$\partial_a H_b\f$, of the damped
87 : * harmonic source function H. Using notation from damped_harmonic_h(), we
88 : * rewrite the same as:
89 : *
90 : * \f{align*}
91 : * \partial_a H_b =& \partial_a T_1 + \partial_a T_2 + \partial_a T_3, \\
92 : * H_a =& T_1 + T_2 + T_3,
93 : * \f}
94 : *
95 : * where:
96 : *
97 : * \f{align*}
98 : * T_1 =& [1 - R(t)] H_a^\mathrm{init}, \\
99 : * T_2 =& [\mu_{L1} \mathrm{log}(\sqrt{\gamma}/\alpha) + \mu_{L2}
100 : * \mathrm{log}(1/\alpha)] n_a,
101 : * \\
102 : * T_3 =& - \mu_S \gamma_{ai} \beta^i / \alpha.
103 : * \f}
104 : *
105 : * Derivation:
106 : *
107 : * \f$\blacksquare\f$ For \f$ T_1 \f$, the derivatives are:
108 : * \f{align*}
109 : * \partial_a T_1 = (1 - R(t))
110 : * \partial_a H_b^\mathrm{init}
111 : * - H_b^\mathrm{init} \partial_a R.
112 : * \f}
113 : *
114 : * \f$\blacksquare\f$ Write \f$ T_2 \equiv (\mu_1 + \mu_2) n_b \f$. Then:
115 : *
116 : * \f{align*}
117 : * \partial_a T_2 =& (\partial_a \mu_1 + \partial_a \mu_2) n_b \\
118 : * +& (\mu_1 + \mu_2) \partial_a n_b,
119 : * \f}
120 : *
121 : * where
122 : *
123 : * \f{align*}
124 : * \partial_a n_b =& \left(-\partial_a \alpha, 0, 0, 0\right) \\
125 : *
126 : * \partial_a \mu_1
127 : * =& \partial_a [A_{L1} R(t) W(x^i) \mathrm{log}(\sqrt{\gamma}/\alpha)^{e_{L1}
128 : * +
129 : * 1}], \\
130 : * =& A_{L1} R(t) W(x^i) \partial_a [\mathrm{log}(\sqrt{\gamma}/\alpha)^{e_{L1}
131 : * +
132 : * 1}] \\
133 : * +& A_{L1} \mathrm{log}(\sqrt{\gamma}/\alpha)^{e_{L1} + 1} \partial_a [R(t)
134 : * W(x^i)],\\
135 : *
136 : * \partial_a \mu_2
137 : * =& \partial_a [A_{L2} R(t) W(x^i) \mathrm{log}(1/\alpha)^{e_{L2} + 1}], \\
138 : * =& A_{L2} R(t) W(x^i) \partial_a [\mathrm{log}(1/\alpha)^{e_{L2} + 1}] \\
139 : * +& A_{L2} \mathrm{log}(1/\alpha)^{e_{L2} + 1} \partial_a [R(t) W(x^i)],
140 : * \f}
141 : *
142 : * where \f$\partial_a [R W] = \left(\partial_0 R(t), \partial_i
143 : * W(x^j)\right)\f$.
144 : *
145 : * \f$\blacksquare\f$ Finally, the derivatives of \f$ T_3 \f$ are:
146 : *
147 : * \f[
148 : * \partial_a T_3 = -\partial_a(\mu_S/\alpha) \gamma_{bi} \beta^i
149 : * -(\mu_S/\alpha) \partial_a(\gamma_{bi}) \beta^i
150 : * -(\mu_S/\alpha) \gamma_{bi}\partial_a \beta^i,
151 : * \f]
152 : *
153 : * where
154 : *
155 : * \f{align*}
156 : * \partial_a(\mu_S / \alpha) =& (1/\alpha)\partial_a \mu_S
157 : * - \frac{\mu_S}{\alpha^2}\partial_a \alpha,
158 : * \,\,\mathrm{and}\\ \partial_a \mu_S =& \partial_a [A_S R(t) W(x^i)
159 : * \mathrm{log}(\sqrt{\gamma}/\alpha)^{e_S}], \\
160 : * =& A_S R(t) W(x^i) \partial_a
161 : * [\mathrm{log}(\sqrt{\gamma}/\alpha)^{e_S}] \\
162 : * +& A_S \mathrm{log}(\sqrt{\gamma} / \alpha)^{e_S} \partial_a
163 : * [R(t) W(x^i)]. \f}
164 : */
165 : template <size_t SpatialDim, typename Frame>
166 1 : void damped_harmonic_rollon(
167 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*> gauge_h,
168 : gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame>*> d4_gauge_h,
169 : const tnsr::a<DataVector, SpatialDim, Frame>& gauge_h_init,
170 : const tnsr::ab<DataVector, SpatialDim, Frame>& dgauge_h_init,
171 : const Scalar<DataVector>& lapse,
172 : const tnsr::I<DataVector, SpatialDim, Frame>& shift,
173 : const Scalar<DataVector>& sqrt_det_spatial_metric,
174 : const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
175 : const tnsr::abb<DataVector, SpatialDim, Frame>& d4_spacetime_metric,
176 : const Scalar<DataVector>& half_pi_two_normals,
177 : const tnsr::i<DataVector, SpatialDim, Frame>& half_phi_two_normals,
178 : const tnsr::aa<DataVector, SpatialDim, Frame>& spacetime_metric,
179 : const tnsr::iaa<DataVector, SpatialDim, Frame>& phi, double time,
180 : const tnsr::I<DataVector, SpatialDim, Frame>& coords, double amp_coef_L1,
181 : double amp_coef_L2, double amp_coef_S, int exp_L1, int exp_L2, int exp_S,
182 : double rollon_start_time, double rollon_width, double sigma_r);
183 :
184 : /*!
185 : * \copydoc damped_harmonic_rollon()
186 : */
187 : template <size_t SpatialDim, typename Frame>
188 1 : void damped_harmonic(
189 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*> gauge_h,
190 : gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame>*> d4_gauge_h,
191 : const Scalar<DataVector>& lapse,
192 : const tnsr::I<DataVector, SpatialDim, Frame>& shift,
193 : const Scalar<DataVector>& sqrt_det_spatial_metric,
194 : const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
195 : const tnsr::abb<DataVector, SpatialDim, Frame>& d4_spacetime_metric,
196 : const Scalar<DataVector>& half_pi_two_normals,
197 : const tnsr::i<DataVector, SpatialDim, Frame>& half_phi_two_normals,
198 : const tnsr::aa<DataVector, SpatialDim, Frame>& spacetime_metric,
199 : const tnsr::iaa<DataVector, SpatialDim, Frame>& phi,
200 : const tnsr::I<DataVector, SpatialDim, Frame>& coords, double amp_coef_L1,
201 : double amp_coef_L2, double amp_coef_S, int exp_L1, int exp_L2, int exp_S,
202 : double sigma_r);
203 :
204 : /*!
205 : * \brief Impose damped harmonic gauge.
206 : *
207 : * \see `damped_harmonic()`
208 : */
209 1 : class DampedHarmonic final : public GaugeCondition {
210 : public:
211 : /// The width of the Gaussian for the spatial decay of the damped harmonic
212 : /// gauge.
213 1 : struct SpatialDecayWidth {
214 0 : using type = double;
215 0 : static constexpr Options::String help{
216 : "Spatial width (sigma_r) of weight function (W(x^i)) used in the "
217 : "damped harmonic gauge."};
218 : };
219 : /// The amplitudes for the L1, L2, and S terms, respectively, for the damped
220 : /// harmonic gauge.
221 1 : struct Amplitudes {
222 0 : using type = std::array<double, 3>;
223 0 : static constexpr Options::String help{
224 : "Amplitudes [A_{L1}, A_{L2}, A_{S}] for the damped harmonic gauge."};
225 : };
226 : /// The exponents for the L1, L2, and S terms, respectively, for the damped
227 : /// harmonic gauge.
228 1 : struct Exponents {
229 0 : using type = std::array<int, 3>;
230 0 : static constexpr Options::String help{
231 : "Exponents [e_{L1}, e_{L2}, e_{S}] for the damped harmonic gauge."};
232 : };
233 :
234 0 : static constexpr Options::String help{
235 : "Apply damped harmonic/damped wave gauge."};
236 :
237 0 : using options = tmpl::list<SpatialDecayWidth, Amplitudes, Exponents>;
238 :
239 0 : DampedHarmonic(double width, const std::array<double, 3>& amps,
240 : const std::array<int, 3>& exps);
241 :
242 0 : DampedHarmonic() = default;
243 0 : DampedHarmonic(const DampedHarmonic&) = default;
244 0 : DampedHarmonic& operator=(const DampedHarmonic&) = default;
245 0 : DampedHarmonic(DampedHarmonic&&) = default;
246 0 : DampedHarmonic& operator=(DampedHarmonic&&) = default;
247 0 : ~DampedHarmonic() override = default;
248 :
249 : /// \cond
250 : explicit DampedHarmonic(CkMigrateMessage* msg);
251 : using PUP::able::register_constructor;
252 : WRAPPED_PUPable_decl_template(DampedHarmonic); // NOLINT
253 : /// \endcond
254 :
255 0 : std::unique_ptr<GaugeCondition> get_clone() const override;
256 :
257 : template <size_t SpatialDim>
258 0 : void gauge_and_spacetime_derivative(
259 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame::Inertial>*> gauge_h,
260 : gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame::Inertial>*>
261 : d4_gauge_h,
262 : const Scalar<DataVector>& lapse,
263 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& shift,
264 : const Scalar<DataVector>& sqrt_det_spatial_metric,
265 : const tnsr::II<DataVector, SpatialDim, Frame::Inertial>&
266 : inverse_spatial_metric,
267 : const tnsr::abb<DataVector, SpatialDim, Frame::Inertial>&
268 : d4_spacetime_metric,
269 : const Scalar<DataVector>& half_pi_two_normals,
270 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
271 : half_phi_two_normals,
272 : const tnsr::aa<DataVector, SpatialDim, Frame::Inertial>& spacetime_metric,
273 : const tnsr::iaa<DataVector, SpatialDim, Frame::Inertial>& phi,
274 : double time,
275 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& inertial_coords)
276 : const;
277 :
278 : // NOLINTNEXTLINE(google-runtime-references)
279 0 : void pup(PUP::er& p) override;
280 :
281 : private:
282 0 : double spatial_decay_width_{std::numeric_limits<double>::signaling_NaN()};
283 0 : std::array<double, 3> amplitudes_{
284 : {std::numeric_limits<double>::signaling_NaN(),
285 : std::numeric_limits<double>::signaling_NaN(),
286 : std::numeric_limits<double>::signaling_NaN()}};
287 0 : std::array<int, 3> exponents_{{std::numeric_limits<int>::max(),
288 : std::numeric_limits<int>::max(),
289 : std::numeric_limits<int>::max()}};
290 : };
291 :
292 : namespace DampedHarmonicGauge_detail {
293 : // Used in the test
294 : double roll_on_function(double time, double t_start, double sigma_t);
295 :
296 : double time_deriv_of_roll_on_function(double time, double t_start,
297 : double sigma_t);
298 : } // namespace DampedHarmonicGauge_detail
299 : } // namespace gauges
300 : } // namespace gh
|
