Add full GAUGE 2-7 support to Z4C C RHS kernel (z4c_rhs_c.C)

Previously only GAUGE 0 and 1 were supported with a compile error for 2-7.
Now supports all 8 gauge choices matching BSSN Fortran formulas:
- GAUGE 2: variable-eta gamma-driver, chi-sqrt denominator
- GAUGE 3: variable-eta gamma-driver, chi-linear denominator
- GAUGE 4: first-order variable-eta, chi-sqrt denominator
- GAUGE 5: first-order variable-eta, chi-linear denominator
- GAUGE 6: Jason's rational position-dependent damping
- GAUGE 7: Jason's exponential position-dependent damping
Also fixes dtSf advection/dissipation guards for gauges where dtSf is active.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

AMSS_NCKU_source/z4c_rhs_c.C

@@ -327,6 +327,9 @@ static int compute_rhs_z4c_cartesian(
     double Axxx[all], Axxy[all], Axxz[all], Axyx[all], Axyy[all], Axyz[all];
     double Axzx[all], Axzy[all], Axzz[all], Ayyx[all], Ayyy[all], Ayyz[all];
     double Ayzx[all], Ayzy[all], Ayzz[all], Azzx[all], Azzy[all], Azzz[all];
+#if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5)
+    double reta[all];
+#endif
 
     const double SSS[3] = {1.0, 1.0, 1.0};
     const double AAS[3] = {-1.0, -1.0, 1.0};
@@ -476,8 +479,181 @@ static int compute_rhs_z4c_cartesian(
         dtSfy_rhs[idx] = ZEO;
         dtSfz_rhs[idx] = ZEO;
     }
+#elif (GAUGE == 2)
+    /* Variable-eta gamma-driver, chi-sqrt denominator */
+    for (int idx = 0; idx < all; ++idx)
+    {
+        const double chin1i = chin1[idx];
+        const double det = gxx[idx] * gyy[idx] * gzz[idx]
+                         + gxy[idx] * gyz[idx] * gxz[idx] * 2.0
+                         - gxz[idx] * gyy[idx] * gxz[idx]
+                         - gxy[idx] * gxy[idx] * gzz[idx]
+                         - gxx[idx] * gyz[idx] * gyz[idx];
+        const double idet = ONE / det;
+        const double upxx = (gyy[idx] * gzz[idx] - gyz[idx] * gyz[idx]) * idet;
+        const double upxy = -(gxy[idx] * gzz[idx] - gyz[idx] * gxz[idx]) * idet;
+        const double upxz = (gxy[idx] * gyz[idx] - gyy[idx] * gxz[idx]) * idet;
+        const double upyy = (gxx[idx] * gzz[idx] - gxz[idx] * gxz[idx]) * idet;
+        const double upyz = -(gxx[idx] * gyz[idx] - gxy[idx] * gxz[idx]) * idet;
+        const double upzz = (gxx[idx] * gyy[idx] - gxy[idx] * gxy[idx]) * idet;
+        const double grdchi2 =
+            upxx * chix[idx] * chix[idx] + upyy * chiy[idx] * chiy[idx] + upzz * chiz[idx] * chiz[idx]
+          + TWO * (upxy * chix[idx] * chiy[idx] + upxz * chix[idx] * chiz[idx] + upyz * chiy[idx] * chiz[idx]);
+        const double sqchi = sqrt(chin1i);
+        reta[idx] = 1.31 / TWO * sqrt(grdchi2 / chin1i) / ((ONE - sqchi) * (ONE - sqchi));
+        betax_rhs[idx] = FF * dtSfx[idx];
+        betay_rhs[idx] = FF * dtSfy[idx];
+        betaz_rhs[idx] = FF * dtSfz[idx];
+        dtSfx_rhs[idx] = Gamx_rhs[idx] - reta[idx] * dtSfx[idx];
+        dtSfy_rhs[idx] = Gamy_rhs[idx] - reta[idx] * dtSfy[idx];
+        dtSfz_rhs[idx] = Gamz_rhs[idx] - reta[idx] * dtSfz[idx];
+    }
+#elif (GAUGE == 3)
+    /* Variable-eta gamma-driver, chi-linear denominator */
+    for (int idx = 0; idx < all; ++idx)
+    {
+        const double chin1i = chin1[idx];
+        const double det = gxx[idx] * gyy[idx] * gzz[idx]
+                         + gxy[idx] * gyz[idx] * gxz[idx] * 2.0
+                         - gxz[idx] * gyy[idx] * gxz[idx]
+                         - gxy[idx] * gxy[idx] * gzz[idx]
+                         - gxx[idx] * gyz[idx] * gyz[idx];
+        const double idet = ONE / det;
+        const double upxx = (gyy[idx] * gzz[idx] - gyz[idx] * gyz[idx]) * idet;
+        const double upxy = -(gxy[idx] * gzz[idx] - gyz[idx] * gxz[idx]) * idet;
+        const double upxz = (gxy[idx] * gyz[idx] - gyy[idx] * gxz[idx]) * idet;
+        const double upyy = (gxx[idx] * gzz[idx] - gxz[idx] * gxz[idx]) * idet;
+        const double upyz = -(gxx[idx] * gyz[idx] - gxy[idx] * gxz[idx]) * idet;
+        const double upzz = (gxx[idx] * gyy[idx] - gxy[idx] * gxy[idx]) * idet;
+        const double grdchi2 =
+            upxx * chix[idx] * chix[idx] + upyy * chiy[idx] * chiy[idx] + upzz * chiz[idx] * chiz[idx]
+          + TWO * (upxy * chix[idx] * chiy[idx] + upxz * chix[idx] * chiz[idx] + upyz * chiy[idx] * chiz[idx]);
+        reta[idx] = 1.31 / TWO * sqrt(grdchi2 / chin1i) / ((ONE - chin1i) * (ONE - chin1i));
+        betax_rhs[idx] = FF * dtSfx[idx];
+        betay_rhs[idx] = FF * dtSfy[idx];
+        betaz_rhs[idx] = FF * dtSfz[idx];
+        dtSfx_rhs[idx] = Gamx_rhs[idx] - reta[idx] * dtSfx[idx];
+        dtSfy_rhs[idx] = Gamy_rhs[idx] - reta[idx] * dtSfy[idx];
+        dtSfz_rhs[idx] = Gamz_rhs[idx] - reta[idx] * dtSfz[idx];
+    }
+#elif (GAUGE == 4)
+    /* Variable-eta gamma-driver, first-order, chi-sqrt denominator */
+    for (int idx = 0; idx < all; ++idx)
+    {
+        const double chin1i = chin1[idx];
+        const double det = gxx[idx] * gyy[idx] * gzz[idx]
+                         + gxy[idx] * gyz[idx] * gxz[idx] * 2.0
+                         - gxz[idx] * gyy[idx] * gxz[idx]
+                         - gxy[idx] * gxy[idx] * gzz[idx]
+                         - gxx[idx] * gyz[idx] * gyz[idx];
+        const double idet = ONE / det;
+        const double upxx = (gyy[idx] * gzz[idx] - gyz[idx] * gyz[idx]) * idet;
+        const double upxy = -(gxy[idx] * gzz[idx] - gyz[idx] * gxz[idx]) * idet;
+        const double upxz = (gxy[idx] * gyz[idx] - gyy[idx] * gxz[idx]) * idet;
+        const double upyy = (gxx[idx] * gzz[idx] - gxz[idx] * gxz[idx]) * idet;
+        const double upyz = -(gxx[idx] * gyz[idx] - gxy[idx] * gxz[idx]) * idet;
+        const double upzz = (gxx[idx] * gyy[idx] - gxy[idx] * gxy[idx]) * idet;
+        const double grdchi2 =
+            upxx * chix[idx] * chix[idx] + upyy * chiy[idx] * chiy[idx] + upzz * chiz[idx] * chiz[idx]
+          + TWO * (upxy * chix[idx] * chiy[idx] + upxz * chix[idx] * chiz[idx] + upyz * chiy[idx] * chiz[idx]);
+        const double sqchi = sqrt(chin1i);
+        reta[idx] = 1.31 / TWO * sqrt(grdchi2 / chin1i) / ((ONE - sqchi) * (ONE - sqchi));
+        betax_rhs[idx] = Gamx_rhs[idx] - reta[idx] * betax[idx];
+        betay_rhs[idx] = Gamy_rhs[idx] - reta[idx] * betay[idx];
+        betaz_rhs[idx] = Gamz_rhs[idx] - reta[idx] * betaz[idx];
+        dtSfx_rhs[idx] = ZEO;
+        dtSfy_rhs[idx] = ZEO;
+        dtSfz_rhs[idx] = ZEO;
+    }
+#elif (GAUGE == 5)
+    /* Variable-eta gamma-driver, first-order, chi-linear denominator */
+    for (int idx = 0; idx < all; ++idx)
+    {
+        const double chin1i = chin1[idx];
+        const double det = gxx[idx] * gyy[idx] * gzz[idx]
+                         + gxy[idx] * gyz[idx] * gxz[idx] * 2.0
+                         - gxz[idx] * gyy[idx] * gxz[idx]
+                         - gxy[idx] * gxy[idx] * gzz[idx]
+                         - gxx[idx] * gyz[idx] * gyz[idx];
+        const double idet = ONE / det;
+        const double upxx = (gyy[idx] * gzz[idx] - gyz[idx] * gyz[idx]) * idet;
+        const double upxy = -(gxy[idx] * gzz[idx] - gyz[idx] * gxz[idx]) * idet;
+        const double upxz = (gxy[idx] * gyz[idx] - gyy[idx] * gxz[idx]) * idet;
+        const double upyy = (gxx[idx] * gzz[idx] - gxz[idx] * gxz[idx]) * idet;
+        const double upyz = -(gxx[idx] * gyz[idx] - gxy[idx] * gxz[idx]) * idet;
+        const double upzz = (gxx[idx] * gyy[idx] - gxy[idx] * gxy[idx]) * idet;
+        const double grdchi2 =
+            upxx * chix[idx] * chix[idx] + upyy * chiy[idx] * chiy[idx] + upzz * chiz[idx] * chiz[idx]
+          + TWO * (upxy * chix[idx] * chiy[idx] + upxz * chix[idx] * chiz[idx] + upyz * chiy[idx] * chiz[idx]);
+        reta[idx] = 1.31 / TWO * sqrt(grdchi2 / chin1i) / ((ONE - chin1i) * (ONE - chin1i));
+        betax_rhs[idx] = Gamx_rhs[idx] - reta[idx] * betax[idx];
+        betay_rhs[idx] = Gamy_rhs[idx] - reta[idx] * betay[idx];
+        betaz_rhs[idx] = Gamz_rhs[idx] - reta[idx] * betaz[idx];
+        dtSfx_rhs[idx] = ZEO;
+        dtSfy_rhs[idx] = ZEO;
+        dtSfz_rhs[idx] = ZEO;
+    }
+#elif (GAUGE == 6 || GAUGE == 7)
+    {
+        /* Jason's position-dependent damping: rational (6) or exponential (7) */
+        int BHN = 0;
+        double Porg[9] = {0.0};
+        double Mass[3] = {0.0};
+    #ifdef fortran1
+        extern "C" { void getpbh(int &, double *, double *); }
+    #elif defined(fortran2)
+        extern "C" { void GETPBH(int &, double *, double *); }
+    #else
+        extern "C" { void getpbh_(int &, double *, double *); }
+    #endif
+        {
+    #ifdef fortran1
+            getpbh(BHN, Porg, Mass);
+    #elif defined(fortran2)
+            GETPBH(BHN, Porg, Mass);
+    #else
+            getpbh_(BHN, Porg, Mass);
+    #endif
+        }
+        if (BHN == 2)
+        {
+            const double M = Mass[0] + Mass[1];
+            const double A = 2.0 / M;
+            const double w1 = 12.0, w2 = 12.0;
+            const double C1 = 1.0 / Mass[0] - A;
+            const double C2 = 1.0 / Mass[1] - A;
+            const double BH_sep2 = (Porg[3] - Porg[0]) * (Porg[3] - Porg[0])
+                                 + (Porg[4] - Porg[1]) * (Porg[4] - Porg[1])
+                                 + (Porg[5] - Porg[2]) * (Porg[5] - Porg[2]);
+            const double inv_BH_sep2 = 1.0 / BH_sep2;
+            for (int k0 = 0; k0 < nz; ++k0) {
+            for (int j0 = 0; j0 < ny; ++j0) {
+            for (int i0 = 0; i0 < nx; ++i0) {
+                const size_t idx = idx_ex(i0, j0, k0, ex);
+                const double xp = X[i0], yp = Y[j0], zp = Z[k0];
+                const double r1 = ((Porg[0]-xp)*(Porg[0]-xp) + (Porg[1]-yp)*(Porg[1]-yp) + (Porg[2]-zp)*(Porg[2]-zp)) * inv_BH_sep2;
+                const double r2 = ((Porg[3]-xp)*(Porg[3]-xp) + (Porg[4]-yp)*(Porg[4]-yp) + (Porg[5]-zp)*(Porg[5]-zp)) * inv_BH_sep2;
+            #if (GAUGE == 6)
+                const double reta_val = A + C1 / (1.0 + w1 * r1) + C2 / (1.0 + w2 * r2);
+            #else
+                const double reta_val = A + C1 * exp(-w1 * r1) + C2 * exp(-w2 * r2);
+            #endif
+                betax_rhs[idx] = FF * dtSfx[idx];
+                betay_rhs[idx] = FF * dtSfy[idx];
+                betaz_rhs[idx] = FF * dtSfz[idx];
+                dtSfx_rhs[idx] = Gamx_rhs[idx] - reta_val * dtSfx[idx];
+                dtSfy_rhs[idx] = Gamy_rhs[idx] - reta_val * dtSfy[idx];
+                dtSfz_rhs[idx] = Gamz_rhs[idx] - reta_val * dtSfz[idx];
+            }}}
+        }
+        else
+        {
+            fprintf(stderr, "z4c_rhs_c: GAUGE %d requires BHN=2, got BHN=%d\n", (int)GAUGE, BHN);
+            return 1;
+        }
+    }
 #else
-#error "z4c_rhs_c.C currently supports GAUGE == 0 or GAUGE == 1 for Z4C"
+#error "z4c_rhs_c.C: unsupported GAUGE value"
 #endif
 
     lopsided(ex, X, Y, Z, gxx, gxx_rhs, betax, betay, betaz, Symmetry, SSS);
@@ -505,7 +681,7 @@ static int compute_rhs_z4c_cartesian(
     lopsided(ex, X, Y, Z, betax, betax_rhs, betax, betay, betaz, Symmetry, ASS);
     lopsided(ex, X, Y, Z, betay, betay_rhs, betax, betay, betaz, Symmetry, SAS);
     lopsided(ex, X, Y, Z, betaz, betaz_rhs, betax, betay, betaz, Symmetry, SSA);
-#if (GAUGE == 0)
+#if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
     lopsided(ex, X, Y, Z, dtSfx, dtSfx_rhs, betax, betay, betaz, Symmetry, ASS);
     lopsided(ex, X, Y, Z, dtSfy, dtSfy_rhs, betax, betay, betaz, Symmetry, SAS);
     lopsided(ex, X, Y, Z, dtSfz, dtSfz_rhs, betax, betay, betaz, Symmetry, SSA);
@@ -552,7 +728,7 @@ static int compute_rhs_z4c_cartesian(
         kodis(ex, X, Y, Z, betax, betax_rhs, ASS, Symmetry, eps);
         kodis(ex, X, Y, Z, betay, betay_rhs, SAS, Symmetry, eps);
         kodis(ex, X, Y, Z, betaz, betaz_rhs, SSA, Symmetry, eps);
-#if (GAUGE == 0)
+#if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
         kodis(ex, X, Y, Z, dtSfx, dtSfx_rhs, ASS, Symmetry, eps);
         kodis(ex, X, Y, Z, dtSfy, dtSfy_rhs, SAS, Symmetry, eps);
         kodis(ex, X, Y, Z, dtSfz, dtSfz_rhs, SSA, Symmetry, eps);