Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 : #pragma once
4 :
5 : #include <cmath>
6 : #include <cstddef>
7 :
8 : #include "DataStructures/DataBox/AsAccess.hpp"
9 : #include "DataStructures/DataBox/DataBox.hpp"
10 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
11 : #include "DataStructures/DataBox/Prefixes.hpp"
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/TaggedContainers.hpp"
14 : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "DataStructures/Variables.hpp"
17 : #include "DataStructures/VectorImpl.hpp"
18 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
19 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
20 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
21 : #include "Evolution/Systems/Ccz4/BoundaryConditions/BoundaryCondition.hpp"
22 : #include "Evolution/Systems/Ccz4/BoundaryConditions/Sommerfeld.hpp"
23 : #include "Evolution/Systems/Ccz4/Christoffel.hpp"
24 : #include "Evolution/Systems/Ccz4/DerivChristoffel.hpp"
25 : #include "Evolution/Systems/Ccz4/DerivLapse.hpp"
26 : #include "Evolution/Systems/Ccz4/DerivZ4Constraint.hpp"
27 : #include "Evolution/Systems/Ccz4/FiniteDifference/BoundaryConditionGhostData.hpp"
28 : #include "Evolution/Systems/Ccz4/FiniteDifference/Derivatives.hpp"
29 : #include "Evolution/Systems/Ccz4/FiniteDifference/Tags.hpp"
30 : #include "Evolution/Systems/Ccz4/Ricci.hpp"
31 : #include "Evolution/Systems/Ccz4/RicciScalarPlusDivergenceZ4Constraint.hpp"
32 : #include "Evolution/Systems/Ccz4/Tags.hpp"
33 : #include "Evolution/Systems/Ccz4/TagsDeclarations.hpp"
34 : #include "Evolution/Systems/Ccz4/TempTags.hpp"
35 : #include "Evolution/Systems/Ccz4/TimeDerivative.hpp"
36 : #include "Evolution/Systems/Ccz4/Z4Constraint.hpp"
37 : #include "Utilities/CallWithDynamicType.hpp"
38 : #include "Utilities/ConstantExpressions.hpp"
39 : #include "Utilities/ContainerHelpers.hpp"
40 : #include "Utilities/ErrorHandling/Assert.hpp"
41 : #include "Utilities/ErrorHandling/Error.hpp"
42 : #include "Utilities/Gsl.hpp"
43 : #include "Utilities/TMPL.hpp"
44 :
45 : /*!
46 : * \brief The namespace for evolving the second-order Ccz4 system.
47 : * Spatial derivatives are computed using 4-th order finite
48 : * differencing. Currently this system only works in 3D.
49 : */
50 : namespace Ccz4::fd {
51 0 : const size_t Dim = 3;
52 :
53 : namespace detail {
54 : // Calculate the time derivative of the evolved variables for the second-order
55 : // Ccz4 system. There is quite some overlap between this apply() funcion
56 : // and the apply() function in the first-order Ccz4 system. However,
57 : // it is not straightforward to directly reuse the first-order
58 : // apply() function as the function signatures differ significantly.
59 : template <size_t Dim>
60 : static void apply(
61 : // LHS time derivatives of evolved variables: eq 4(a) - 4(i)
62 : const gsl::not_null<tnsr::ii<DataVector, Dim>*> dt_conformal_spatial_metric,
63 : const gsl::not_null<Scalar<DataVector>*> dt_lapse,
64 : const gsl::not_null<tnsr::I<DataVector, Dim>*> dt_shift,
65 : const gsl::not_null<Scalar<DataVector>*> dt_conformal_factor,
66 : const gsl::not_null<tnsr::ii<DataVector, Dim>*> dt_a_tilde,
67 : const gsl::not_null<Scalar<DataVector>*> dt_trace_extrinsic_curvature,
68 : const gsl::not_null<Scalar<DataVector>*> dt_theta,
69 : const gsl::not_null<tnsr::I<DataVector, Dim>*> dt_gamma_hat,
70 : const gsl::not_null<tnsr::I<DataVector, Dim>*> dt_b,
71 :
72 : // quantities we need for computing eq 4, 13-27
73 : const gsl::not_null<Scalar<DataVector>*> conformal_factor_squared,
74 : const gsl::not_null<Scalar<DataVector>*> det_conformal_spatial_metric,
75 : const gsl::not_null<tnsr::II<DataVector, Dim>*>
76 : inv_conformal_spatial_metric,
77 : const gsl::not_null<tnsr::II<DataVector, Dim>*> inv_spatial_metric,
78 : const gsl::not_null<tnsr::II<DataVector, Dim>*> inv_a_tilde,
79 : // temporary expressions
80 : const gsl::not_null<tnsr::ij<DataVector, Dim>*> a_tilde_times_field_b,
81 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
82 : a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde,
83 : const gsl::not_null<Scalar<DataVector>*> contracted_field_b,
84 : const gsl::not_null<tnsr::ijK<DataVector, Dim>*> symmetrized_d_field_b,
85 : const gsl::not_null<tnsr::i<DataVector, Dim>*>
86 : contracted_symmetrized_d_field_b,
87 : const gsl::not_null<tnsr::i<DataVector, Dim>*> field_d_up_times_a_tilde,
88 : const gsl::not_null<tnsr::I<DataVector, Dim>*>
89 : contracted_field_d_up, // temp for eq 18 -20
90 : const gsl::not_null<Scalar<DataVector>*>
91 : half_conformal_factor_squared, // temp for eq 25
92 : const gsl::not_null<tnsr::ij<DataVector, Dim>*>
93 : conformal_metric_times_field_b,
94 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
95 : conformal_metric_times_trace_a_tilde,
96 : const gsl::not_null<tnsr::i<DataVector, Dim>*>
97 : inv_conformal_metric_times_d_a_tilde,
98 : const gsl::not_null<tnsr::I<DataVector, Dim>*>
99 : gamma_hat_minus_contracted_conformal_christoffel,
100 : const gsl::not_null<tnsr::iJ<DataVector, Dim>*>
101 : d_gamma_hat_minus_contracted_conformal_christoffel,
102 : const gsl::not_null<tnsr::i<DataVector, Dim>*>
103 : contracted_christoffel_second_kind, // temp for eq 18 -20
104 : const gsl::not_null<tnsr::ij<DataVector, Dim>*>
105 : contracted_d_conformal_christoffel_difference, // temp for eq 18 -20
106 : const gsl::not_null<Scalar<DataVector>*> k_minus_2_theta_c,
107 : const gsl::not_null<Scalar<DataVector>*> k_minus_k0_minus_2_theta_c,
108 : const gsl::not_null<tnsr::ii<DataVector, Dim>*> lapse_times_a_tilde,
109 : const gsl::not_null<tnsr::ijj<DataVector, Dim>*> lapse_times_d_a_tilde,
110 : const gsl::not_null<tnsr::i<DataVector, Dim>*> lapse_times_field_a,
111 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
112 : lapse_times_conformal_spatial_metric,
113 : const gsl::not_null<Scalar<DataVector>*> lapse_times_slicing_condition,
114 : const gsl::not_null<Scalar<DataVector>*>
115 : lapse_times_ricci_scalar_plus_divergence_z4_constraint,
116 : const gsl::not_null<tnsr::I<DataVector, Dim>*> shift_times_deriv_gamma_hat,
117 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
118 : inv_tau_times_conformal_metric,
119 : // expressions and identities needed for evolution equations: eq 13 - 27
120 : const gsl::not_null<Scalar<DataVector>*> trace_a_tilde, // eq 13
121 : const gsl::not_null<tnsr::iJJ<DataVector, Dim>*> field_d_up, // eq 14
122 : const gsl::not_null<tnsr::Ijj<DataVector, Dim>*>
123 : conformal_christoffel_second_kind, // eq 15
124 : const gsl::not_null<tnsr::iJkk<DataVector, Dim>*>
125 : d_conformal_christoffel_second_kind, // eq 16
126 : const gsl::not_null<tnsr::Ijj<DataVector, Dim>*>
127 : christoffel_second_kind, // eq 17
128 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
129 : spatial_ricci_tensor, // eq 18 - 20
130 : const gsl::not_null<tnsr::ij<DataVector, Dim>*> grad_grad_lapse, // eq 21
131 : const gsl::not_null<Scalar<DataVector>*> divergence_lapse, // eq 22
132 : const gsl::not_null<tnsr::I<DataVector, Dim>*>
133 : contracted_conformal_christoffel_second_kind, // eq 23
134 : const gsl::not_null<tnsr::iJ<DataVector, Dim>*>
135 : d_contracted_conformal_christoffel_second_kind, // eq 24
136 : const gsl::not_null<tnsr::i<DataVector, Dim>*>
137 : spatial_z4_constraint, // eq 25
138 : const gsl::not_null<tnsr::I<DataVector, Dim>*>
139 : upper_spatial_z4_constraint, // eq 25
140 : const gsl::not_null<tnsr::ij<DataVector, Dim>*>
141 : grad_spatial_z4_constraint, // eq 26
142 : const gsl::not_null<Scalar<DataVector>*>
143 : ricci_scalar_plus_divergence_z4_constraint, // eq 27
144 :
145 : // fixed params for SO-CCZ4
146 : // c = 1.0, cleaning_speed = 1.0
147 : // one_over_relaxation_time = 0.0
148 : const double c, const double cleaning_speed, // e in the paper
149 : const double one_over_relaxation_time, // \tau^{-1}
150 :
151 : // free params for SO-CCZ4
152 : const Scalar<DataVector>& eta, const double f, const double kappa_1,
153 : const double kappa_2, const double kappa_3, const Scalar<DataVector>& k_0,
154 :
155 : // evolved variables
156 : const tnsr::ii<DataVector, Dim>& conformal_spatial_metric,
157 : const Scalar<DataVector>& lapse, const tnsr::I<DataVector, Dim>& shift,
158 : const Scalar<DataVector>& conformal_factor,
159 : const tnsr::ii<DataVector, Dim>& a_tilde,
160 : const Scalar<DataVector>& trace_extrinsic_curvature,
161 : const Scalar<DataVector>& theta, const tnsr::I<DataVector, Dim>& gamma_hat,
162 : const tnsr::I<DataVector, Dim>& b,
163 :
164 : // auxilliary fields and their derivatives
165 : const tnsr::i<DataVector, Dim>& field_a,
166 : const tnsr::iJ<DataVector, Dim>& field_b,
167 : const tnsr::ijj<DataVector, Dim>& field_d,
168 : const tnsr::i<DataVector, Dim>& field_p,
169 : const tnsr::ii<DataVector, Dim>& d_field_a,
170 : const tnsr::iiJ<DataVector, Dim>& d_field_b,
171 : const tnsr::iijj<DataVector, Dim>& d_field_d,
172 : const tnsr::ii<DataVector, Dim>& d_field_p,
173 :
174 : // spatial derivatives of other evolved variables
175 : const tnsr::ijj<DataVector, Dim>& d_a_tilde,
176 : const tnsr::i<DataVector, Dim>& d_trace_extrinsic_curvature,
177 : const tnsr::i<DataVector, Dim>& d_theta,
178 : const tnsr::iJ<DataVector, Dim>& d_gamma_hat,
179 : const tnsr::iJ<DataVector, Dim>& d_b, const bool shifting_shift,
180 : const bool evolve_lapse_and_shift) {
181 : constexpr double one_third = 1.0 / 3.0;
182 : // quantities we need for computing eq 4, 13 - 27
183 :
184 : determinant_and_inverse(det_conformal_spatial_metric,
185 : inv_conformal_spatial_metric,
186 : conformal_spatial_metric);
187 :
188 : get(*conformal_factor_squared) = square(get(conformal_factor));
189 :
190 : ::tenex::evaluate<ti::I, ti::J>(
191 : inv_spatial_metric, (*conformal_factor_squared)() *
192 : (*inv_conformal_spatial_metric)(ti::I, ti::J));
193 :
194 : ::tenex::evaluate<ti::I, ti::J>(
195 : inv_a_tilde, a_tilde(ti::k, ti::l) *
196 : (*inv_conformal_spatial_metric)(ti::I, ti::K) *
197 : (*inv_conformal_spatial_metric)(ti::J, ti::L));
198 :
199 : ASSERT(min(get(lapse)) > 0.0,
200 : "The lapse must be positive when using 1+log slicing but is: "
201 : << get(lapse));
202 :
203 : // slicing_condition and d_slicing_condition is not used in SO-CCZ4
204 : // \alpha g(\alpha) == 2
205 : *lapse_times_slicing_condition =
206 : make_with_value<Scalar<DataVector>>(lapse, 2.0);
207 :
208 : // expressions and identities needed for evolution equations: eq 13 - 27
209 :
210 : // eq 13
211 : ::tenex::evaluate(
212 : trace_a_tilde,
213 : (*inv_conformal_spatial_metric)(ti::I, ti::J) * a_tilde(ti::i, ti::j));
214 :
215 : // from eq 14: field_d_up is the D_k^{ij} tensor
216 : ::tenex::evaluate<ti::k, ti::I, ti::J>(
217 : field_d_up, (*inv_conformal_spatial_metric)(ti::I, ti::N) *
218 : (*inv_conformal_spatial_metric)(ti::M, ti::J) *
219 : field_d(ti::k, ti::n, ti::m));
220 :
221 : // eq 15
222 : ::Ccz4::conformal_christoffel_second_kind(conformal_christoffel_second_kind,
223 : *inv_conformal_spatial_metric,
224 : field_d);
225 :
226 : // eq 16
227 : ::Ccz4::deriv_conformal_christoffel_second_kind(
228 : d_conformal_christoffel_second_kind, *inv_conformal_spatial_metric,
229 : field_d, d_field_d, *field_d_up);
230 :
231 : // eq 17
232 : ::Ccz4::christoffel_second_kind(christoffel_second_kind,
233 : conformal_spatial_metric,
234 : *inv_conformal_spatial_metric, field_p,
235 : *conformal_christoffel_second_kind);
236 :
237 : // temporary expressions needed for eq 18 - 20
238 : ::tenex::evaluate<ti::l>(contracted_christoffel_second_kind,
239 : (*christoffel_second_kind)(ti::M, ti::l, ti::m));
240 :
241 : // comment for the future: we should probably ensure the traces are taken
242 : // before computing the differences as the off-diagonal terms are not
243 : // needed
244 : ::tenex::evaluate<ti::i, ti::j>(
245 : contracted_d_conformal_christoffel_difference,
246 : (*d_conformal_christoffel_second_kind)(ti::m, ti::M, ti::i, ti::j) -
247 : (*d_conformal_christoffel_second_kind)(ti::j, ti::M, ti::i, ti::m));
248 :
249 : ::tenex::evaluate<ti::L>(contracted_field_d_up,
250 : (*field_d_up)(ti::m, ti::M, ti::L));
251 :
252 : // eq 18 - 20
253 : ::Ccz4::spatial_ricci_tensor(
254 : spatial_ricci_tensor, *christoffel_second_kind,
255 : *contracted_christoffel_second_kind,
256 : *contracted_d_conformal_christoffel_difference, conformal_spatial_metric,
257 : *inv_conformal_spatial_metric, field_d, *field_d_up,
258 : *contracted_field_d_up, field_p, d_field_p);
259 :
260 : // eq 21
261 : ::Ccz4::grad_grad_lapse(grad_grad_lapse, lapse, *christoffel_second_kind,
262 : field_a, d_field_a);
263 :
264 : // eq 22
265 : ::Ccz4::divergence_lapse(divergence_lapse, *conformal_factor_squared,
266 : *inv_conformal_spatial_metric, *grad_grad_lapse);
267 :
268 : // eq 23
269 : ::Ccz4::contracted_conformal_christoffel_second_kind(
270 : contracted_conformal_christoffel_second_kind,
271 : *inv_conformal_spatial_metric, *conformal_christoffel_second_kind);
272 :
273 : // eq 24
274 : ::Ccz4::deriv_contracted_conformal_christoffel_second_kind(
275 : d_contracted_conformal_christoffel_second_kind,
276 : *inv_conformal_spatial_metric, *field_d_up,
277 : *conformal_christoffel_second_kind, *d_conformal_christoffel_second_kind);
278 :
279 : // temp for eq 25
280 : ::tenex::evaluate<ti::I>(
281 : gamma_hat_minus_contracted_conformal_christoffel,
282 : gamma_hat(ti::I) -
283 : (*contracted_conformal_christoffel_second_kind)(ti::I));
284 :
285 : // eq 25
286 : ::Ccz4::spatial_z4_constraint(
287 : spatial_z4_constraint, conformal_spatial_metric,
288 : *gamma_hat_minus_contracted_conformal_christoffel);
289 :
290 : // temp for eq 25
291 : ::tenex::evaluate(half_conformal_factor_squared,
292 : 0.5 * (*conformal_factor_squared)());
293 :
294 : // eq 25
295 : ::Ccz4::upper_spatial_z4_constraint(
296 : upper_spatial_z4_constraint, *half_conformal_factor_squared,
297 : *gamma_hat_minus_contracted_conformal_christoffel);
298 :
299 : // temp for eq 26
300 : ::tenex::evaluate<ti::i, ti::L>(
301 : d_gamma_hat_minus_contracted_conformal_christoffel,
302 : d_gamma_hat(ti::i, ti::L) -
303 : (*d_contracted_conformal_christoffel_second_kind)(ti::i, ti::L));
304 :
305 : // eq 26
306 : ::Ccz4::grad_spatial_z4_constraint(
307 : grad_spatial_z4_constraint, *spatial_z4_constraint,
308 : conformal_spatial_metric, *christoffel_second_kind, field_d,
309 : *gamma_hat_minus_contracted_conformal_christoffel,
310 : *d_gamma_hat_minus_contracted_conformal_christoffel);
311 :
312 : // eq 27
313 : ::Ccz4::ricci_scalar_plus_divergence_z4_constraint(
314 : ricci_scalar_plus_divergence_z4_constraint, *conformal_factor_squared,
315 : *inv_conformal_spatial_metric, *spatial_ricci_tensor,
316 : *grad_spatial_z4_constraint);
317 :
318 : // temporary expressions not already computed above
319 :
320 : ::tenex::evaluate(contracted_field_b, field_b(ti::k, ti::K));
321 :
322 : ::tenex::evaluate<ti::k, ti::j, ti::I>(
323 : symmetrized_d_field_b,
324 : 0.5 * (d_field_b(ti::k, ti::j, ti::I) + d_field_b(ti::j, ti::k, ti::I)));
325 :
326 : ::tenex::evaluate<ti::k>(contracted_symmetrized_d_field_b,
327 : (*symmetrized_d_field_b)(ti::k, ti::i, ti::I));
328 :
329 : ::tenex::evaluate<ti::k>(
330 : field_d_up_times_a_tilde,
331 : (*field_d_up)(ti::k, ti::I, ti::J) * a_tilde(ti::i, ti::j));
332 :
333 : ::tenex::evaluate<ti::i, ti::j>(
334 : conformal_metric_times_field_b,
335 : conformal_spatial_metric(ti::k, ti::i) * field_b(ti::j, ti::K));
336 :
337 : ::tenex::evaluate<ti::i, ti::j>(
338 : conformal_metric_times_trace_a_tilde,
339 : conformal_spatial_metric(ti::i, ti::j) * (*trace_a_tilde)());
340 :
341 : ::tenex::evaluate<ti::k>(inv_conformal_metric_times_d_a_tilde,
342 : (*inv_conformal_spatial_metric)(ti::I, ti::J) *
343 : d_a_tilde(ti::k, ti::i, ti::j));
344 :
345 : ::tenex::evaluate<ti::i, ti::j>(
346 : a_tilde_times_field_b, a_tilde(ti::k, ti::i) * field_b(ti::j, ti::K));
347 :
348 : ::tenex::evaluate<ti::i, ti::j>(
349 : a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde,
350 : a_tilde(ti::i, ti::j) -
351 : one_third * (*conformal_metric_times_trace_a_tilde)(ti::i, ti::j));
352 :
353 : ::tenex::evaluate(k_minus_2_theta_c,
354 : trace_extrinsic_curvature() - 2.0 * c * theta());
355 :
356 : ::tenex::evaluate(k_minus_k0_minus_2_theta_c, (*k_minus_2_theta_c)() - k_0());
357 :
358 : ::tenex::evaluate<ti::i, ti::j>(lapse_times_a_tilde,
359 : (lapse)() * a_tilde(ti::i, ti::j));
360 :
361 : tenex::evaluate<ti::k, ti::i, ti::j>(
362 : lapse_times_d_a_tilde, (lapse)() * d_a_tilde(ti::k, ti::i, ti::j));
363 :
364 : ::tenex::evaluate<ti::k>(lapse_times_field_a, (lapse)() * field_a(ti::k));
365 :
366 : ::tenex::evaluate<ti::i, ti::j>(
367 : lapse_times_conformal_spatial_metric,
368 : (lapse)() * conformal_spatial_metric(ti::i, ti::j));
369 :
370 : ::tenex::evaluate(
371 : lapse_times_ricci_scalar_plus_divergence_z4_constraint,
372 : (lapse)() * (*ricci_scalar_plus_divergence_z4_constraint)());
373 :
374 : ::tenex::evaluate<ti::I>(shift_times_deriv_gamma_hat,
375 : shift(ti::K) * d_gamma_hat(ti::k, ti::I));
376 :
377 : ::tenex::evaluate<ti::i, ti::j>(
378 : inv_tau_times_conformal_metric,
379 : one_over_relaxation_time * conformal_spatial_metric(ti::i, ti::j));
380 :
381 : // time derivative computation: eq 4(a) - 4(i)
382 : // The way we evaluate the following may seem weird compared
383 : // to 4(a) - 4(i). This is because we try to reuse the first-order
384 : // Ccz4 codes as much as possible. Hence, the following is written
385 : // based on 12a-12i, with s=1, and red terms canceled.
386 :
387 : // eq 12a : time derivative of the conformal spatial metric
388 : // this reduces to eq 4a for our choice of \tau^{-1}=0.
389 : ::tenex::evaluate<ti::i, ti::j>(
390 : dt_conformal_spatial_metric,
391 : 2.0 * shift(ti::K) * field_d(ti::k, ti::i, ti::j) +
392 : (*conformal_metric_times_field_b)(ti::i, ti::j) +
393 : (*conformal_metric_times_field_b)(ti::j, ti::i) -
394 : 2.0 * one_third * conformal_spatial_metric(ti::i, ti::j) *
395 : (*contracted_field_b)() -
396 : 2.0 * (lapse)() *
397 : (*a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde)(
398 : ti::i, ti::j) -
399 : (*inv_tau_times_conformal_metric)(ti::i, ti::j) *
400 : ((*det_conformal_spatial_metric)() - 1.0));
401 :
402 : if (evolve_lapse_and_shift) {
403 : // time derivative of the lapse
404 : ::tenex::evaluate<>(dt_lapse, (shift(ti::K) * field_a(ti::k) -
405 : (*lapse_times_slicing_condition)() *
406 : (*k_minus_k0_minus_2_theta_c)()) *
407 : lapse());
408 : // time derivative of the shift
409 : ::tenex::evaluate<ti::I>(dt_shift, f * b(ti::I));
410 : if (shifting_shift) {
411 : ::tenex::update<ti::I>(
412 : dt_shift, (*dt_shift)(ti::I) + shift(ti::K) * field_b(ti::k, ti::I));
413 : }
414 : } else {
415 : *dt_lapse = make_with_value<Scalar<DataVector>>(lapse, 0.0);
416 : *dt_shift = make_with_value<tnsr::I<DataVector, Dim>>(shift, 0.0);
417 : }
418 :
419 : // time derivative of the the conformal factor
420 : ::tenex::evaluate(dt_conformal_factor,
421 : (shift(ti::K) * field_p(ti::k) +
422 : one_third * ((lapse)() * trace_extrinsic_curvature() -
423 : (*contracted_field_b)())) *
424 : conformal_factor());
425 :
426 : // time derivative of the trace-free part of the extrinsic curvature
427 : ::tenex::evaluate<ti::i, ti::j>(
428 : dt_a_tilde,
429 : shift(ti::K) * d_a_tilde(ti::k, ti::i, ti::j) +
430 : (*conformal_factor_squared)() *
431 : ((lapse)() * ((*spatial_ricci_tensor)(ti::i, ti::j) +
432 : (*grad_spatial_z4_constraint)(ti::i, ti::j) +
433 : (*grad_spatial_z4_constraint)(ti::j, ti::i)) -
434 : (*grad_grad_lapse)(ti::i, ti::j)) -
435 : one_third * conformal_spatial_metric(ti::i, ti::j) *
436 : ((*lapse_times_ricci_scalar_plus_divergence_z4_constraint)() -
437 : (*divergence_lapse)()) +
438 : (*a_tilde_times_field_b)(ti::i, ti::j) +
439 : (*a_tilde_times_field_b)(ti::j, ti::i) -
440 : 2.0 * one_third * a_tilde(ti::i, ti::j) * (*contracted_field_b)() +
441 : (*lapse_times_a_tilde)(ti::i, ti::j) * (*k_minus_2_theta_c)() -
442 : 2.0 * (*lapse_times_a_tilde)(ti::i, ti::l) *
443 : (*inv_conformal_spatial_metric)(ti::L, ti::M) *
444 : a_tilde(ti::m, ti::j) -
445 : (*inv_tau_times_conformal_metric)(ti::i, ti::j) * (*trace_a_tilde)());
446 :
447 : // time derivative of the trace of the extrinsic curvature
448 : ::tenex::evaluate(
449 : dt_trace_extrinsic_curvature,
450 : shift(ti::K) * d_trace_extrinsic_curvature(ti::k) -
451 : (*divergence_lapse)() +
452 : (*lapse_times_ricci_scalar_plus_divergence_z4_constraint)() +
453 : (lapse)() * (trace_extrinsic_curvature() * (*k_minus_2_theta_c)() -
454 : 3.0 * kappa_1 * (1.0 + kappa_2) * theta()));
455 :
456 : // time derivative of the projection of the Z4 four-vector along
457 : // the normal direction
458 : ::tenex::evaluate(
459 : dt_theta,
460 : shift(ti::K) * d_theta(ti::k) +
461 : (lapse)() *
462 : (0.5 * square(cleaning_speed) *
463 : ((*ricci_scalar_plus_divergence_z4_constraint)() +
464 : 2.0 * one_third * square(trace_extrinsic_curvature()) -
465 : a_tilde(ti::i, ti::j) * (*inv_a_tilde)(ti::I, ti::J)) -
466 : c * theta() * trace_extrinsic_curvature() -
467 : (*upper_spatial_z4_constraint)(ti::I)*field_a(ti::i) -
468 : kappa_1 * (2.0 + kappa_2) * theta()));
469 :
470 : // time derivative \hat{\Gamma}^i
471 : // first, compute terms without s
472 : ::tenex::evaluate<ti::I>(
473 : dt_gamma_hat,
474 : // terms without lapse nor s
475 : (*shift_times_deriv_gamma_hat)(ti::I) +
476 : 2.0 * one_third *
477 : (*contracted_conformal_christoffel_second_kind)(ti::I) *
478 : (*contracted_field_b)() -
479 : (*contracted_conformal_christoffel_second_kind)(ti::K)*field_b(
480 : ti::k, ti::I) +
481 : 2.0 * kappa_3 * (*spatial_z4_constraint)(ti::j) *
482 : (2.0 * one_third * (*inv_conformal_spatial_metric)(ti::I, ti::J) *
483 : (*contracted_field_b)() -
484 : (*inv_conformal_spatial_metric)(ti::J, ti::K) *
485 : field_b(ti::k, ti::I)) +
486 : // terms with lapse but not s
487 : 2.0 * (lapse)() *
488 : (-2.0 * one_third *
489 : (*inv_conformal_spatial_metric)(ti::I, ti::J) *
490 : d_trace_extrinsic_curvature(ti::j) +
491 : (*inv_conformal_spatial_metric)(ti::K, ti::I) * d_theta(ti::k) +
492 : (*conformal_christoffel_second_kind)(ti::I, ti::j, ti::k) *
493 : (*inv_a_tilde)(ti::J, ti::K) -
494 : 3.0 * (*inv_a_tilde)(ti::I, ti::J) * field_p(ti::j) -
495 : (*inv_conformal_spatial_metric)(ti::K, ti::I) *
496 : (theta() * field_a(ti::k) +
497 : 2.0 * one_third * trace_extrinsic_curvature() *
498 : (*spatial_z4_constraint)(ti::k)) -
499 : (*inv_a_tilde)(ti::I, ti::J) * field_a(ti::j) -
500 : kappa_1 * (*inv_conformal_spatial_metric)(ti::I, ti::J) *
501 : (*spatial_z4_constraint)(ti::j)));
502 : // We add the following since s=1 (Gamma-driver gauge) is assumed in SoCcz4
503 : ::tenex::update<ti::I>(dt_gamma_hat,
504 : (*dt_gamma_hat)(ti::I) +
505 : // red terms should cancel
506 : // terms with s but not lapse
507 : (*inv_conformal_spatial_metric)(ti::K, ti::L) *
508 : (*symmetrized_d_field_b)(ti::k, ti::l, ti::I) +
509 : one_third *
510 : (*inv_conformal_spatial_metric)(ti::I, ti::K) *
511 : (*contracted_symmetrized_d_field_b)(ti::k));
512 :
513 : // time derivative b^i
514 : if (evolve_lapse_and_shift) {
515 : ::tenex::evaluate<ti::I>(dt_b, (*dt_gamma_hat)(ti::I)-eta() * b(ti::I));
516 : if (shifting_shift) {
517 : ::tenex::update<ti::I>(
518 : dt_b, (*dt_b)(ti::I) + shift(ti::K) * (d_b(ti::k, ti::I) -
519 : d_gamma_hat(ti::k, ti::I)));
520 : }
521 : } else {
522 : (*dt_b).get(0) = 0.0;
523 : (*dt_b).get(1) = 0.0;
524 : (*dt_b).get(2) = 0.0;
525 : }
526 :
527 : // Note that, we do not need to evolve the auxiliary variables in SO-CCZ4.
528 : }
529 : } // namespace detail
530 :
531 : /*!
532 : * \brief Compute the RHS of the second-order CCZ4 formulation of Einstein's
533 : * equations \cite Dumbser2017okk with finite differencing. Also
534 : * update the time derivative at outermost interior points in a sphere domain
535 : * if Sommerfeld boundary conditions are applied.
536 : *
537 : * \details The evolution equations are equations 4(a) - 4(i) of
538 : * \cite Dumbser2017okk Equations 13 - 27 define identities
539 : * used in the evolution equations.
540 : *
541 : * \note Different from the first-order system, the lapse
542 : * and the conformal factor are evolved instead of their natural logs.
543 : * The four auxiliary varialbels \f$A_i\f$, \f$B_k{}^{i}\f$,
544 : * \f$D_{kij}\f$, and \f$P_i\f$ are NOT evolved.
545 : */
546 1 : struct SoTimeDerivative {
547 : template <typename DbTagsList>
548 0 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box) {
549 : // compute the first spatial derivatives of the evolved variables
550 : using evolved_vars_tag = typename System::variables_tag;
551 : using gradients_tags = typename System::gradients_tags;
552 :
553 : // only 4-th order accurate second derivatives have been implemented
554 : // To keep the same order of accuracy we use the same order for
555 : // both first and second derivatives
556 : constexpr size_t fd_order = 4;
557 : const auto& evolved_vars = db::get<evolved_vars_tag>(*box);
558 : const Mesh<Dim>& subcell_mesh =
559 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(*box);
560 : const size_t num_pts = subcell_mesh.number_of_grid_points();
561 : using deriv_var_tag = db::wrap_tags_in<::Tags::deriv, gradients_tags,
562 : tmpl::size_t<Dim>, Frame::Inertial>;
563 : Variables<deriv_var_tag> cell_centered_Ccz4_derivs{num_pts};
564 : const auto& cell_centered_logical_to_inertial_inv_jacobian =
565 : db::get<evolution::dg::subcell::fd::Tags::
566 : InverseJacobianLogicalToInertial<Dim>>(*box);
567 :
568 : constexpr bool subcell_enabled_at_external_boundary =
569 : std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
570 : *box))>::SubcellOptions::subcell_enabled_at_external_boundary;
571 :
572 : const Element<3>& element = db::get<domain::Tags::Element<3>>(*box);
573 :
574 : const Ccz4::fd::Reconstructor& recons =
575 : db::get<Ccz4::fd::Tags::Reconstructor>(*box);
576 :
577 : // If the element has external boundaries and subcell is enabled for
578 : // boundary elements, compute FD ghost data with a given boundary condition.
579 : if constexpr (subcell_enabled_at_external_boundary) {
580 : if (not element.external_boundaries().empty()) {
581 : fd::BoundaryConditionGhostData::apply(box, element, recons);
582 : }
583 : }
584 :
585 : const bool evolve_lapse_and_shift =
586 : get<::Ccz4::fd::Tags::EvolveLapseAndShift>(*box);
587 :
588 : ::Ccz4::fd::spacetime_derivatives(
589 : make_not_null(&cell_centered_Ccz4_derivs), evolved_vars,
590 : db::get<evolution::dg::subcell::Tags::GhostDataForReconstruction<Dim>>(
591 : *box),
592 : fd_order, subcell_mesh, cell_centered_logical_to_inertial_inv_jacobian);
593 :
594 : // calculate the four auxiliary fields in eq. 6
595 : // auxiliary variables NOT evolved in SO-CCZ4
596 : const auto& d_lapse =
597 : get<::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<Dim>,
598 : Frame::Inertial>>(cell_centered_Ccz4_derivs);
599 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(evolved_vars);
600 : const auto field_a = ::tenex::evaluate<ti::i>(d_lapse(ti::i) / lapse());
601 :
602 : const auto& field_b =
603 : get<::Tags::deriv<gr::Tags::Shift<DataVector, Dim>, tmpl::size_t<Dim>,
604 : Frame::Inertial>>(cell_centered_Ccz4_derivs);
605 :
606 : const auto& d_spatial_conformal_metric =
607 : get<::Tags::deriv<::Ccz4::Tags::ConformalMetric<DataVector, Dim>,
608 : tmpl::size_t<Dim>, Frame::Inertial>>(
609 : cell_centered_Ccz4_derivs);
610 : tnsr::ijj<DataVector, Dim> field_d;
611 : ::tenex::evaluate<ti::i, ti::j, ti::k>(
612 : make_not_null(&field_d),
613 : 0.5 * d_spatial_conformal_metric(ti::i, ti::j, ti::k));
614 :
615 : const auto& conformal_factor =
616 : get<::Ccz4::Tags::ConformalFactor<DataVector>>(evolved_vars);
617 : const auto& d_conformal_factor =
618 : get<::Tags::deriv<::Ccz4::Tags::ConformalFactor<DataVector>,
619 : tmpl::size_t<Dim>, Frame::Inertial>>(
620 : cell_centered_Ccz4_derivs);
621 : const auto field_p = ::tenex::evaluate<ti::i>(d_conformal_factor(ti::i) /
622 : conformal_factor());
623 :
624 : // compute second derivatives of the evolved variables
625 : Variables<db::wrap_tags_in<::Tags::second_deriv, gradients_tags,
626 : tmpl::size_t<Dim>, Frame::Inertial>>
627 : cell_centered_Ccz4_second_derivs{num_pts};
628 :
629 : Ccz4::fd::second_spacetime_derivatives(
630 : make_not_null(&cell_centered_Ccz4_second_derivs), evolved_vars,
631 : db::get<evolution::dg::subcell::Tags::GhostDataForReconstruction<Dim>>(
632 : *box),
633 : fd_order, subcell_mesh, cell_centered_logical_to_inertial_inv_jacobian);
634 :
635 : // compute spatial derivative of the four auxiliary fields
636 : const auto& d_d_lapse =
637 : get<::Tags::second_deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<Dim>,
638 : Frame::Inertial>>(
639 : cell_centered_Ccz4_second_derivs);
640 : tnsr::ii<DataVector, Dim> d_field_a;
641 : ::tenex::evaluate<ti::i, ti::j>(
642 : make_not_null(&d_field_a),
643 : (d_d_lapse(ti::i, ti::j) - d_lapse(ti::i) * d_lapse(ti::j) / lapse()) /
644 : lapse());
645 :
646 : const auto& d_field_b =
647 : get<::Tags::second_deriv<gr::Tags::Shift<DataVector, Dim>,
648 : tmpl::size_t<Dim>, Frame::Inertial>>(
649 : cell_centered_Ccz4_second_derivs);
650 :
651 : const auto& d_d_conformal_metric =
652 : get<::Tags::second_deriv<::Ccz4::Tags::ConformalMetric<DataVector, Dim>,
653 : tmpl::size_t<Dim>, Frame::Inertial>>(
654 : cell_centered_Ccz4_second_derivs);
655 :
656 : tnsr::iijj<DataVector, Dim> d_field_d;
657 : ::tenex::evaluate<ti::i, ti::j, ti::k, ti::l>(
658 : make_not_null(&d_field_d),
659 : 0.5 * d_d_conformal_metric(ti::i, ti::j, ti::k, ti::l));
660 :
661 : const auto& d_d_conformal_factor =
662 : get<::Tags::second_deriv<::Ccz4::Tags::ConformalFactor<DataVector>,
663 : tmpl::size_t<Dim>, Frame::Inertial>>(
664 : cell_centered_Ccz4_second_derivs);
665 :
666 : tnsr::ii<DataVector, Dim> d_field_p;
667 : ::tenex::evaluate<ti::i, ti::j>(
668 : make_not_null(&d_field_p),
669 : (d_d_conformal_factor(ti::i, ti::j) - d_conformal_factor(ti::i) *
670 : d_conformal_factor(ti::j) /
671 : conformal_factor()) /
672 : conformal_factor());
673 :
674 : // intialize containers to be supplied in the SO-CCZ4 TimeDerivative.cpp
675 : // apply() function quantities we need for computing eq 4, 13 - 27
676 : using TempVars = Variables<tmpl::list<
677 : ::Ccz4::Tags::ConformalFactorSquared<DataVector>,
678 : ::Ccz4::Tags::DetConformalSpatialMetric<DataVector>,
679 : ::Ccz4::Tags::InverseConformalMetric<DataVector, Dim>,
680 : gr::Tags::InverseSpatialMetric<DataVector, Dim>,
681 : ::Ccz4::Tags::InvATilde<DataVector, Dim>,
682 : ::Ccz4::Tags::ATildeTimesFieldB<DataVector, Dim>,
683 : ::Ccz4::Tags::ATildeMinusOneThirdConformalMetricTimesTraceATilde<
684 : DataVector, Dim>,
685 : ::Ccz4::Tags::ContractedFieldB<DataVector>,
686 : ::Ccz4::Tags::SymmetrizedDerivFieldB<DataVector, Dim>,
687 : ::Ccz4::Tags::ContractedSymmetrizedDerivFieldB<DataVector, Dim>,
688 : ::Ccz4::Tags::FieldDUpTimesATilde<DataVector, Dim>,
689 : ::Ccz4::Tags::ContractedFieldDUp<DataVector, Dim>,
690 : ::Ccz4::Tags::HalfConformalFactorSquared<DataVector>,
691 : ::Ccz4::Tags::ConformalMetricTimesFieldB<DataVector, Dim>,
692 : ::Ccz4::Tags::ConformalMetricTimesTraceATilde<DataVector, Dim>,
693 : ::Ccz4::Tags::InverseConformalMetricTimesDerivATilde<DataVector, Dim>,
694 : ::Ccz4::Tags::GammaHatMinusContractedConformalChristoffel<DataVector,
695 : Dim>,
696 : ::Ccz4::Tags::DerivGammaHatMinusContractedConformalChristoffel<
697 : DataVector, Dim>,
698 : ::Ccz4::Tags::ContractedChristoffelSecondKind<DataVector, Dim>,
699 : ::Ccz4::Tags::ContractedDerivConformalChristoffelDifference<DataVector,
700 : Dim>,
701 : ::Ccz4::Tags::KMinus2ThetaC<DataVector>,
702 : ::Ccz4::Tags::KMinusK0Minus2ThetaC<DataVector>,
703 : ::Ccz4::Tags::LapseTimesATilde<DataVector, Dim>,
704 : ::Ccz4::Tags::LapseTimesDerivATilde<DataVector, Dim>,
705 : ::Ccz4::Tags::LapseTimesFieldA<DataVector, Dim>,
706 : ::Ccz4::Tags::LapseTimesConformalMetric<DataVector, Dim>,
707 : ::Ccz4::Tags::LapseTimesSlicingCondition<DataVector>,
708 : ::Ccz4::Tags::LapseTimesRicciScalarPlus2DivergenceZ4Constraint<
709 : DataVector>,
710 : ::Ccz4::Tags::ShiftTimesDerivGammaHat<DataVector, Dim>,
711 : ::Ccz4::Tags::InverseTauTimesConformalMetric<DataVector, Dim>,
712 : ::Ccz4::Tags::TraceATilde<DataVector>,
713 : ::Ccz4::Tags::FieldDUp<DataVector, Dim>,
714 : ::Ccz4::Tags::ConformalChristoffelSecondKind<DataVector, Dim>,
715 : ::Ccz4::Tags::DerivConformalChristoffelSecondKind<DataVector, Dim>,
716 : ::Ccz4::Tags::ChristoffelSecondKind<DataVector, Dim>,
717 : ::Ccz4::Tags::SpatialRicciTensor<DataVector, Dim>,
718 : ::Ccz4::Tags::GradGradLapse<DataVector, Dim>,
719 : ::Ccz4::Tags::DivergenceLapse<DataVector>,
720 : ::Ccz4::Tags::ContractedConformalChristoffelSecondKind<DataVector, Dim>,
721 : ::Ccz4::Tags::DerivContractedConformalChristoffelSecondKind<DataVector,
722 : Dim>,
723 : ::Ccz4::Tags::SpatialZ4Constraint<DataVector, Dim>,
724 : ::Ccz4::Tags::GradSpatialZ4Constraint<DataVector, Dim>,
725 : ::Ccz4::Tags::RicciScalarPlusDivergenceZ4Constraint<DataVector>>>;
726 :
727 : TempVars temp_vars(num_pts);
728 : tnsr::I<DataVector, Dim> upper_spatial_z4_constraint(num_pts);
729 :
730 : // free params
731 : const double c = 1.0; // c = 1.0 in SO-CCZ4
732 : const double cleaning_speed = 1.0; // e in the paper; e = 1.0 for SO-CCZ4
733 : const Scalar<DataVector>& eta = get<::Ccz4::Tags::Eta<DataVector>>(*box);
734 : const double f = Ccz4::fd::System::f;
735 : const Scalar<DataVector>& k_0 = get<::Ccz4::Tags::K0<DataVector>>(*box);
736 : const double kappa_1 = get<::Ccz4::Tags::Kappa1>(*box);
737 : const double kappa_2 = get<::Ccz4::Tags::Kappa2>(*box);
738 : const double kappa_3 = get<::Ccz4::Tags::Kappa3>(*box);
739 : const double one_over_relaxation_time = 0.0; // \tau^{-1} = 0 in SO-CCZ4
740 : const bool shifting_shift = Ccz4::fd::System::shifting_shift;
741 :
742 : // we assume the databox already has tags corresponding to dt of the evolved
743 : // variables
744 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, evolved_vars_tag>;
745 :
746 : db::mutate<dt_variables_tag,
747 : ::Ccz4::Tags::SpatialZ4ConstraintUp<DataVector, Dim>>(
748 : [&](const auto dt_vars_ptr,
749 : const auto upper_spatial_z4_constraint_ptr) {
750 : // resize here
751 : dt_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
752 : auto& [conformal_factor_squared, det_conformal_spatial_metric,
753 : inv_conformal_spatial_metric, inv_spatial_metric, inv_a_tilde,
754 : a_tilde_times_field_b,
755 : a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde,
756 : contracted_field_b, symmetrized_d_field_b,
757 : contracted_symmetrized_d_field_b, field_d_up_times_a_tilde,
758 : contracted_field_d_up, half_conformal_factor_squared,
759 : conformal_metric_times_field_b,
760 : conformal_metric_times_trace_a_tilde,
761 : inv_conformal_metric_times_d_a_tilde,
762 : gamma_hat_minus_contracted_conformal_christoffel,
763 : d_gamma_hat_minus_contracted_conformal_christoffel,
764 : contracted_christoffel_second_kind,
765 : contracted_d_conformal_christoffel_difference,
766 : k_minus_2_theta_c, k_minus_k0_minus_2_theta_c,
767 : lapse_times_a_tilde, lapse_times_d_a_tilde,
768 : lapse_times_field_a, lapse_times_conformal_spatial_metric,
769 : lapse_times_slicing_condition,
770 : lapse_times_ricci_scalar_plus_divergence_z4_constraint,
771 : shift_times_deriv_gamma_hat, inv_tau_times_conformal_metric,
772 : trace_a_tilde, field_d_up, conformal_christoffel_second_kind,
773 : d_conformal_christoffel_second_kind, christoffel_second_kind,
774 : spatial_ricci_tensor, grad_grad_lapse, divergence_lapse,
775 : contracted_conformal_christoffel_second_kind,
776 : d_contracted_conformal_christoffel_second_kind,
777 : spatial_z4_constraint, grad_spatial_z4_constraint,
778 : ricci_scalar_plus_divergence_z4_constraint] = temp_vars;
779 : detail::apply(
780 : // LHS time derivatives of evolved variables: eq 4a - 4i
781 : make_not_null(
782 : &get<::Tags::dt<
783 : ::Ccz4::Tags::ConformalMetric<DataVector, Dim>>>(
784 : *dt_vars_ptr)), // eq 4a
785 : make_not_null(&get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(
786 : *dt_vars_ptr)), // eq 4g
787 : make_not_null(&get<::Tags::dt<gr::Tags::Shift<DataVector, Dim>>>(
788 : *dt_vars_ptr)), // eq 4h
789 : make_not_null(
790 : &get<::Tags::dt<::Ccz4::Tags::ConformalFactor<DataVector>>>(
791 : *dt_vars_ptr)), // eq 4c
792 : make_not_null(
793 : &get<::Tags::dt<::Ccz4::Tags::ATilde<DataVector, Dim>>>(
794 : *dt_vars_ptr)), // eq 4b
795 : make_not_null(&get<::Tags::dt<
796 : gr::Tags::TraceExtrinsicCurvature<DataVector>>>(
797 : *dt_vars_ptr)), // eq 4d
798 : make_not_null(&get<::Tags::dt<::Ccz4::Tags::Theta<DataVector>>>(
799 : *dt_vars_ptr)), // eq 4e
800 : make_not_null(
801 : &get<::Tags::dt<::Ccz4::Tags::GammaHat<DataVector, Dim>>>(
802 : *dt_vars_ptr)), // eq 4f
803 : make_not_null(
804 : &get<::Tags::dt<
805 : ::Ccz4::Tags::AuxiliaryShiftB<DataVector, Dim>>>(
806 : *dt_vars_ptr)), // eq 4i
807 :
808 : // quantities we need for computing eq 4, 13 - 27
809 : make_not_null(&conformal_factor_squared),
810 : make_not_null(&det_conformal_spatial_metric),
811 : make_not_null(&inv_conformal_spatial_metric),
812 : make_not_null(&inv_spatial_metric), make_not_null(&inv_a_tilde),
813 : // temporary expressions
814 : make_not_null(&a_tilde_times_field_b),
815 : make_not_null(
816 : &a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde),
817 : make_not_null(&contracted_field_b),
818 : make_not_null(&symmetrized_d_field_b),
819 : make_not_null(&contracted_symmetrized_d_field_b),
820 : make_not_null(&field_d_up_times_a_tilde),
821 : make_not_null(&contracted_field_d_up), // temp for eq 18 -20
822 : make_not_null(&half_conformal_factor_squared), // temp for eq 25
823 : make_not_null(&conformal_metric_times_field_b),
824 : make_not_null(&conformal_metric_times_trace_a_tilde),
825 : make_not_null(&inv_conformal_metric_times_d_a_tilde),
826 : make_not_null(&gamma_hat_minus_contracted_conformal_christoffel),
827 : make_not_null(
828 : &d_gamma_hat_minus_contracted_conformal_christoffel),
829 : make_not_null(
830 : &contracted_christoffel_second_kind), // temp for eq 18 -20
831 : make_not_null(
832 : &contracted_d_conformal_christoffel_difference), // temp for
833 : // eq 18 -20
834 : make_not_null(&k_minus_2_theta_c),
835 : make_not_null(&k_minus_k0_minus_2_theta_c),
836 : make_not_null(&lapse_times_a_tilde),
837 : make_not_null(&lapse_times_d_a_tilde),
838 : make_not_null(&lapse_times_field_a),
839 : make_not_null(&lapse_times_conformal_spatial_metric),
840 : make_not_null(&lapse_times_slicing_condition),
841 : make_not_null(
842 : &lapse_times_ricci_scalar_plus_divergence_z4_constraint),
843 : make_not_null(&shift_times_deriv_gamma_hat),
844 : make_not_null(&inv_tau_times_conformal_metric),
845 : // expressions and identities needed for evolution equations: eq
846 : // 13
847 : // - 27
848 : make_not_null(&trace_a_tilde), // eq 13
849 : make_not_null(&field_d_up), // eq 14
850 : make_not_null(&conformal_christoffel_second_kind), // eq 15
851 : make_not_null(&d_conformal_christoffel_second_kind), // eq 16
852 : make_not_null(&christoffel_second_kind), // eq 17
853 : make_not_null(&spatial_ricci_tensor), // eq 18 - 20
854 : make_not_null(&grad_grad_lapse), // eq 21
855 : make_not_null(&divergence_lapse), // eq 22
856 : make_not_null(
857 : &contracted_conformal_christoffel_second_kind), // eq 23
858 : make_not_null(
859 : &d_contracted_conformal_christoffel_second_kind), // eq 24
860 : make_not_null(&spatial_z4_constraint), // eq 25
861 : make_not_null(&upper_spatial_z4_constraint), // eq 25
862 : make_not_null(&grad_spatial_z4_constraint), // eq 26
863 : make_not_null(
864 : &ricci_scalar_plus_divergence_z4_constraint), // eq 27
865 : // fixed params
866 : c, cleaning_speed, // e in the paper
867 : one_over_relaxation_time, // \tau^{-1}
868 : // free params
869 : eta, f, kappa_1, kappa_2, kappa_3, k_0,
870 : // evolved variables
871 : get<::Ccz4::Tags::ConformalMetric<DataVector, Dim>>(evolved_vars),
872 : get<gr::Tags::Lapse<DataVector>>(evolved_vars),
873 : get<gr::Tags::Shift<DataVector, Dim>>(evolved_vars),
874 : get<::Ccz4::Tags::ConformalFactor<DataVector>>(evolved_vars),
875 : get<::Ccz4::Tags::ATilde<DataVector, Dim>>(evolved_vars),
876 : get<gr::Tags::TraceExtrinsicCurvature<DataVector>>(evolved_vars),
877 : get<::Ccz4::Tags::Theta<DataVector>>(evolved_vars),
878 : get<::Ccz4::Tags::GammaHat<DataVector, Dim>>(evolved_vars),
879 : get<::Ccz4::Tags::AuxiliaryShiftB<DataVector, Dim>>(evolved_vars),
880 :
881 : field_a, // auxiliary variables NOT evolved in SO-CCZ4
882 : field_b, field_d, field_p,
883 : d_field_a, // spatial derivative of auxiliary variables
884 : d_field_b, d_field_d, d_field_p,
885 :
886 : // spatial derivatives of other evolved variables
887 : get<::Tags::deriv<::Ccz4::Tags::ATilde<DataVector, Dim>,
888 : tmpl::size_t<Dim>, Frame::Inertial>>(
889 : cell_centered_Ccz4_derivs),
890 : get<::Tags::deriv<gr::Tags::TraceExtrinsicCurvature<DataVector>,
891 : tmpl::size_t<Dim>, Frame::Inertial>>(
892 : cell_centered_Ccz4_derivs),
893 : get<::Tags::deriv<::Ccz4::Tags::Theta<DataVector>,
894 : tmpl::size_t<Dim>, Frame::Inertial>>(
895 : cell_centered_Ccz4_derivs),
896 : get<::Tags::deriv<::Ccz4::Tags::GammaHat<DataVector, Dim>,
897 : tmpl::size_t<Dim>, Frame::Inertial>>(
898 : cell_centered_Ccz4_derivs),
899 : get<::Tags::deriv<::Ccz4::Tags::AuxiliaryShiftB<DataVector, Dim>,
900 : tmpl::size_t<Dim>, Frame::Inertial>>(
901 : cell_centered_Ccz4_derivs),
902 : shifting_shift, evolve_lapse_and_shift);
903 :
904 : *upper_spatial_z4_constraint_ptr =
905 : std::move(upper_spatial_z4_constraint);
906 : },
907 : box);
908 :
909 : // handling Sommerfeld BCs
910 : // WARNING: We assume a complete sphere domain (all wedges)
911 : if constexpr (not subcell_enabled_at_external_boundary) {
912 : return;
913 : }
914 : if (element.external_boundaries().empty()) {
915 : return;
916 : }
917 :
918 : const auto& external_bcs =
919 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(*box).at(
920 : element.id().block_id());
921 : for (const auto& direction : element.external_boundaries()) {
922 : ASSERT(direction == Direction<Dim>::upper_zeta(),
923 : "external bc direction is not upper_zeta but " << direction);
924 : // this is potentially expensive; we should add a boolean flag for
925 : // time-dependent bc. This should be done in a future PR.
926 : if (dynamic_cast<const Ccz4::BoundaryConditions::Sommerfeld*>(
927 : external_bcs.at(direction).get()) == nullptr) {
928 : return; // no need to check more, as we do not allow mixed BCs
929 : }
930 : }
931 :
932 : db::mutate<dt_variables_tag>(
933 : [&](const auto dt_vars_ptr, const auto& inertial_coords,
934 : const auto& evolved_var) {
935 : const size_t num_face_pts =
936 : subcell_mesh.extents(0) * subcell_mesh.extents(1);
937 :
938 : ASSERT(inertial_coords.get(0).size() ==
939 : subcell_mesh.number_of_grid_points() and
940 : inertial_coords.get(1).size() ==
941 : subcell_mesh.number_of_grid_points() and
942 : inertial_coords.get(2).size() ==
943 : subcell_mesh.number_of_grid_points(),
944 : "Inertial coords size ["
945 : << inertial_coords.get(0).size()
946 : << ", " << inertial_coords.get(1).size()
947 : << ", " << inertial_coords.get(2).size()
948 : << "] does not match number of subcell points "
949 : << subcell_mesh.number_of_grid_points());
950 :
951 : std::array<DataVector, Dim> outermost_inertial_coords{};
952 : for (size_t i = 0; i < Dim; ++i) {
953 : make_const_view<DataVector>(
954 : make_not_null(&outermost_inertial_coords.at(i)),
955 : inertial_coords.get(i),
956 : (subcell_mesh.extents(2) - 1) * num_face_pts, num_face_pts);
957 : }
958 : const DataVector outermost_radial_coords =
959 : sqrt(square(outermost_inertial_coords[0]) +
960 : square(outermost_inertial_coords[1]) +
961 : square(outermost_inertial_coords[2]));
962 :
963 : // modify the time derivatives at the outermost interior points
964 : // per Sommerfeld BCs
965 : tmpl::for_each<Ccz4::fd::System::variables_tag_list>(
966 : [&]<typename Tag>(tmpl::type_<Tag> /*meta*/) {
967 : const auto& var = get<Tag>(evolved_var);
968 : const auto& d_var =
969 : get<::Tags::deriv<Tag, tmpl::size_t<Dim>, Frame::Inertial>>(
970 : cell_centered_Ccz4_derivs);
971 : auto& dt_var = get<::Tags::dt<Tag>>(*dt_vars_ptr);
972 :
973 : for (size_t tensor_index = 0; tensor_index < Tag::type::size();
974 : ++tensor_index) {
975 : DataVector outermost_var{};
976 : make_const_view<DataVector>(
977 : make_not_null(&outermost_var), var[tensor_index],
978 : (subcell_mesh.extents(2) - 1) * num_face_pts,
979 : num_face_pts);
980 :
981 : std::array<DataVector, Dim> outermost_d_var{};
982 : for (size_t i = 0; i < Dim; ++i) {
983 : make_const_view<DataVector>(
984 : make_not_null(&outermost_d_var.at(i)),
985 : d_var[Dim * tensor_index + i],
986 : (subcell_mesh.extents(2) - 1) * num_face_pts,
987 : num_face_pts);
988 : }
989 :
990 : const DataVector outermost_dr_var =
991 : (outermost_inertial_coords[0] * outermost_d_var[0] +
992 : outermost_inertial_coords[1] * outermost_d_var[1] +
993 : outermost_inertial_coords[2] * outermost_d_var[2]) /
994 : outermost_radial_coords;
995 : DataVector outermost_dt_var{};
996 : outermost_dt_var.set_data_ref(
997 : dt_var[tensor_index].data() +
998 : (subcell_mesh.extents(2) - 1) * num_face_pts,
999 : num_face_pts);
1000 : outermost_dt_var =
1001 : outermost_var * (-1.0 / outermost_radial_coords) -
1002 : outermost_dr_var;
1003 : }
1004 : });
1005 : },
1006 : box,
1007 : db::get<evolution::dg::subcell::Tags::Coordinates<
1008 : Dim, Frame::Inertial>>(*box),
1009 : db::get<evolved_vars_tag>(*box));
1010 : }
1011 : };
1012 : } // namespace Ccz4::fd
|
