Z4c: Make evolveA, evolveB evolution stationary for Kerr-Schild

Author: eschnett

Z4c/param.ccl

@@ -51,23 +51,29 @@ CCTK_REAL eta "eta" STEERABLE=always
 
 # Integrated ("second order") gauge conditions.
 #
-# With these off, the gauge is the usual first order one:
-#     d/dt alpha  = -alpha f_mu_L Khat
-#     d/dt beta^i = f_mu_S Gamt^i - eta beta^i
+# With these off, the gauge is the usual first order one, advected along the
+# shift:
+#     d/dt alpha  = -alpha f_mu_L Khat         + beta^i d_i alpha
+#     d/dt beta^i = f_mu_S Gamt^i - eta beta^i + beta^j d_j beta^i
 #
 # With them on, d/dt alpha and d/dt beta^i become evolved variables A and B^i
-# (group dtalphaG / dtbetaG), obeying
+# (groups dtalphaG / dtbetaG), the complete coordinate time derivatives:
 #     d/dt alpha = A
-#     d/dt A     = d/dt (-alpha f_mu_L Khat) - alphaDriver (A + alpha f_mu_L Khat)
-# and correspondingly for the shift. The driver terms relax A and B^i towards
-# the first order values; with a driver of zero and A, B^i initialised to
-# those values the first order gauge is reproduced exactly, which is what the
-# unit test checks.
+#     d/dt A     = d/dt (-alpha f_mu_L Khat + beta^i d_i alpha)
+#                  - alphaDriver (A + alpha f_mu_L Khat - beta^i d_i alpha)
+# and correspondingly for the shift, where d/dt on the right hand side is the
+# complete time derivative (d/dt alpha = A, d/dt beta^i = B^i, and
+# d/dt Khat = Kh_rhs + beta^i d_i Khat). The driver terms relax A and B^i
+# towards the first order values; with a driver of zero and A, B^i initialised
+# to those values the first order gauge is reproduced exactly.
 #
 # The point of the integrated form is that A and B^i can be initialised from
 # ADMBaseX::dtlapse and ADMBaseX::dtshift. For stationary initial data in
-# horizon-penetrating coordinates those are zero, so the gauge starts at rest
-# instead of being kicked by a nonzero -alpha f_mu_L Khat at t = 0.
+# horizon-penetrating coordinates those are zero, and since d/dt Khat and
+# d/dt Gamt^i vanish as well, every right hand side vanishes: the data stay
+# stationary up to truncation error instead of the gauge being kicked by a
+# nonzero -alpha f_mu_L Khat at t = 0. A nonzero driver re-introduces that
+# kick, since the first order values are not zero for such data.
 
 BOOLEAN evolveA "Evolve A = d/dt alpha instead of setting it algebraically"
 {
@@ -97,6 +103,48 @@ CCTK_REAL alphaG_floor "Floor for alphaG" STEERABLE=always
   (0:* :: ""
 } 1.0e-10
 
+# Clamping of extreme values. Kerr-Schild data with a "stuffed" interior
+# collapse there in finite coordinate time when the lapse does not collapse
+# with them (the stuffed shift makes Khat ~ 1/(alpha rtilde(0)), and
+# d/dt K ~ alpha K^2). Clamping bounds that garbage: it acts only where values
+# are extreme, which for such data is inside the horizon, where any
+# modification is acceptable; a healthy exterior never reaches the bounds.
+# Clamps introduce kinks, which the Kreiss-Oliger dissipation smooths.
+
+CCTK_REAL alpha_max "Ceiling for the lapse alpha; 0 disables" STEERABLE=always
+{
+  0.0 :: "disabled"
+  (0.0:* :: "clamp alpha to at most this value"
+} 0.0
+
+CCTK_REAL clamp_max "Clamp |Khat|, |Atilde_ij|, |Gamma_tilde^i|, |Theta|, |A| and |B^i| (all of dimension 1/M) to at most this value; 0 disables" STEERABLE=always
+{
+  0.0 :: "disabled"
+  (0.0:* :: "clamp to this value"
+} 0.0
+
+# Once the lapse has collapsed inside, only the shift-driven terms remain
+# active there; an unbounded shift then drives the conformal metric through
+# the Lie-derivative terms. These two bound that path.
+
+CCTK_REAL beta_max "Clamp |beta^i| (each component) to at most this value; 0 disables" STEERABLE=always
+{
+  0.0 :: "disabled"
+  (0.0:* :: "clamp to this value"
+} 0.0
+
+CCTK_REAL gammat_max "Clamp |gamma_tilde_ij - delta_ij| (each component) to at most this value, before and after det gamma_tilde = 1 is enforced; where delta + gamma_tilde is not positive definite it is reset to flat; 0 disables" STEERABLE=always
+{
+  0.0 :: "disabled"
+  (0.0:* :: "clamp to this value"
+} 0.0
+
+CCTK_REAL chi_max "Ceiling for chi (the physical metric is gamma_tilde / chi); 0 disables" STEERABLE=always
+{
+  0.0 :: "disabled"
+  (0.0:* :: "clamp chi to at most this value"
+} 0.0
+
 CCTK_REAL epsdiss "Dissipation coefficient <arXiv:gr-qc/0610128>" STEERABLE=always
 {
   0.0:* :: ""

Z4c/schedule.ccl

@@ -147,12 +147,24 @@ SCHEDULE Z4c_Enforce IN Z4c_PostStepGroup
   LANG: C
   READS: chi(interior)
   READS: gamma_tilde(interior)
+  READS: K_hat(interior)
   READS: A_tilde(interior)
+  READS: Gam_tilde(interior)
+  READS: Theta(interior)
   READS: alphaG(interior)
+  READS: betaG(interior)
+  READS: dtalphaG(interior)
+  READS: dtbetaG(interior)
   WRITES: chi(interior)
   WRITES: gamma_tilde(interior)
+  WRITES: K_hat(interior)
   WRITES: A_tilde(interior)
+  WRITES: Gam_tilde(interior)
+  WRITES: Theta(interior)
   WRITES: alphaG(interior)
+  WRITES: betaG(interior)
+  WRITES: dtalphaG(interior)
+  WRITES: dtbetaG(interior)
   SYNC: chi
   SYNC: gamma_tilde
   SYNC: K_hat

Z4c/src/adm.cxx

@@ -133,16 +133,16 @@ extern "C" void Z4c_ADM(CCTK_ARGUMENTS) {
   nvtxRangeEnd(range);
 #endif
 
-  // vars.dtalpha and vars.dtbeta are only the source terms of the gauge
-  // conditions. The conditions themselves are advected along the shift,
+  // Without evolveA, vars.dtalpha is only the source term of the lapse
+  // condition. The condition itself is advected along the shift,
   //
   //     d/dt alpha  = -alpha f_mu_L Khat + beta^i d_i alpha
   //
-  // and correspondingly for the shift, so ADMBaseX::dtlapse and dtshift are
-  // wrong by beta^i d_i alpha (a few percent near a black hole horizon)
-  // unless that term is added here as well. It is added for both gauges: with
-  // evolveA the stored value is A, and the lapse is still advected on top of
-  // it.
+  // and correspondingly for the shift, so ADMBaseX::dtlapse and dtshift would
+  // be wrong by beta^i d_i alpha (a few percent near a black hole horizon)
+  // unless that term is added here as well. With evolveA the stored value is
+  // A, which is the complete d/dt alpha (rhs.cxx does not advect alphaG on
+  // top of it), so nothing is added.
   //
   // Kreiss-Oliger dissipation is deliberately NOT added. It is applied to the
   // evolved right hand sides in rhs.cxx because it stabilises the discrete
@@ -155,10 +155,12 @@ extern "C" void Z4c_ADM(CCTK_ARGUMENTS) {
   // applies the configured CarpetX outer boundary condition; only with
   // CarpetX::boundary_* = "none" does the outer boundary layer keep the
   // source term alone.
-  apply_upwind(cctkGH, gf_alphaG1, gf_betaG1, gf_dtalp1);
+  if (!evolveA)
+    apply_upwind(cctkGH, gf_alphaG1, gf_betaG1, gf_dtalp1);
 
-  for (int a = 0; a < 3; ++a)
-    apply_upwind(cctkGH, gf_betaG1(a), gf_betaG1, gf_dtbeta1(a));
+  if (!evolveB)
+    for (int a = 0; a < 3; ++a)
+      apply_upwind(cctkGH, gf_betaG1(a), gf_betaG1, gf_dtbeta1(a));
 }
 
 } // namespace Z4c

Z4c/src/adm2.cxx

@@ -246,25 +246,29 @@ extern "C" void Z4c_ADM2(CCTK_ARGUMENTS) {
         // beta^k d_k K_ij, which is added below with apply_upwind.
         gf_dtk1.store(mask, index1, vars.K_rhs);
 
-        // Z4c_ADM reports the advected gauge conditions,
+        // Z4c_ADM reports the complete time derivatives of lapse and shift,
         //
-        //     D   = d/dt alpha  = S   + beta^i d_i alpha
-        //     D^a = d/dt beta^a = S^a + beta^j d_j beta^a
+        //     D = d/dt alpha,   D^a = d/dt beta^a.
         //
-        // so their time derivatives are
+        // Without evolveA, D = S + beta^i d_i alpha with the source
+        // S = -alpha f_mu_L Khat, so
         //
         //     d/dt D   = d/dt S   + (d/dt beta^i) d_i alpha  + beta^i d_i D
         //     d/dt D^a = d/dt S^a + (d/dt beta^j) d_j beta^a + beta^j d_j D^a
         //
-        // z4c_vars' dtalpha_target_rhs and dtbeta_target_rhs cannot be used
-        // for the d/dt S terms: they differentiate the sources using the
-        // pre-advection right hand sides held by z4c_vars. The sources
-        // are recomputed here from the advected derivatives instead. Khat and
-        // Gamt^i are advected with the centred derivatives that z4c_vars
-        // carries rather than the upwinded ones rhs.cxx uses, and Kreiss-
-        // Oliger dissipation is excluded throughout, so these remain the
-        // derivatives of the continuum gauge conditions rather than of the
-        // discrete update.
+        // where d/dt S needs the complete d/dt Khat = Kh_rhs + beta^i d_i Khat.
+        // Khat and Gamt^i are advected with the centred derivatives that
+        // z4c_vars carries rather than the upwinded ones rhs.cxx uses, and
+        // Kreiss-Oliger dissipation is excluded throughout, so these remain
+        // the derivatives of the continuum gauge conditions rather than of the
+        // discrete update. The beta^i d_i D term is added below with
+        // apply_upwind.
+        //
+        // With evolveA, D = A and d/dt D = A_rhs + beta^i d_i A: z4c_vars'
+        // A_rhs is the complete d/dt of the advective lapse condition,
+        // (d/dt beta^i) d_i alpha included, except for the advection of A
+        // itself, which is added below with apply_upwind. Likewise for the
+        // shift.
         const vreal D = gf_dtalp1(mask, index1);
         const vec<vreal, 3> Da = gf_dtbeta1(mask, index1);
 
@@ -278,27 +282,23 @@ extern "C" void Z4c_ADM2(CCTK_ARGUMENTS) {
                  });
         });
 
-        // With evolveA / evolveB the source is the evolved A (B^i), whose
-        // full time derivative is A_rhs plus the advection term beta^i d_i A
-        // that rhs.cxx adds; that advection term is added below with
-        // apply_upwind.
-        const vreal dS =
+        const vreal dt2alp =
             evolveA ? vars.A_rhs
-                    : -f_mu_L * (D * vars.Kh + (1 + vars.alphaG) * dtKh);
-
-        const vec<vreal, 3> dSa([&](int a) ARITH_INLINE {
-          return evolveB ? vars.B_rhs(a) : f_mu_S * dtGamt(a) - eta * Da(a);
+                    : -f_mu_L * (D * vars.Kh + (1 + vars.alphaG) * dtKh) +
+                          sum<3>([&](int i) ARITH_INLINE {
+                            return Da(i) * vars.dalphaG(i);
+                          });
+
+        const vec<vreal, 3> dt2beta([&](int a) ARITH_INLINE {
+          return evolveB ? vars.B_rhs(a)
+                         : f_mu_S * dtGamt(a) - eta * Da(a) +
+                               sum<3>([&](int j) ARITH_INLINE {
+                                 return Da(j) * vars.dbetaG(a)(j);
+                               });
         });
 
-        gf_dt2alp1.store(mask, index1, dS + sum<3>([&](int i) ARITH_INLINE {
-                                         return Da(i) * vars.dalphaG(i);
-                                       }));
-
-        gf_dt2beta1.store(mask, index1, vec<vreal, 3>([&](int a) ARITH_INLINE {
-                            return dSa(a) + sum<3>([&](int j) ARITH_INLINE {
-                                     return Da(j) * vars.dbetaG(a)(j);
-                                   });
-                          }));
+        gf_dt2alp1.store(mask, index1, dt2alp);
+        gf_dt2beta1.store(mask, index1, dt2beta);
       });
 #ifdef __CUDACC__
   nvtxRangeEnd(range);
@@ -313,15 +313,18 @@ extern "C" void Z4c_ADM2(CCTK_ARGUMENTS) {
     for (int b = a; b < 3; ++b)
       apply_upwind(cctkGH, gf_k1(a, b), gf_betaG1, gf_dtk1(a, b));
 
-  // beta^i d_i D and beta^j d_j D^a
-  apply_upwind(cctkGH, gf_dtalp1, gf_betaG1, gf_dt2alp1);
+  // beta^i d_i D and beta^j d_j D^a, completing d/dt D and d/dt D^a when A,
+  // B^i are not evolved
+  if (!evolveA)
+    apply_upwind(cctkGH, gf_dtalp1, gf_betaG1, gf_dt2alp1);
 
-  for (int a = 0; a < 3; ++a)
-    apply_upwind(cctkGH, gf_dtbeta1(a), gf_betaG1, gf_dt2beta1(a));
+  if (!evolveB)
+    for (int a = 0; a < 3; ++a)
+      apply_upwind(cctkGH, gf_dtbeta1(a), gf_betaG1, gf_dt2beta1(a));
 
-  // With evolveA / evolveB the evolved A and B^i are themselves advected in
-  // rhs.cxx, so d/dt A = A_rhs + beta^i d_i A enters d/dt D; the main loop
-  // stored only A_rhs. Mirrors the evolveA / evolveB guards in rhs.cxx.
+  // With evolveA / evolveB, D = A and D^a = B^a, which are advected in
+  // rhs.cxx, so d/dt D = A_rhs + beta^i d_i A; the main loop stored only
+  // A_rhs. Mirrors the evolveA / evolveB guards in rhs.cxx.
   if (evolveA)
     apply_upwind(cctkGH, gf_A1, gf_betaG1, gf_dt2alp1);