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