64-BitBrainstorm_2026/AMSS-NCKU
8
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source 
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>
This commit is contained in:
CGH0S7
2026-05-14 13:00:07 +08:00
parent ac92b75dc1
commit 320a3694c2
1 changed files with 179 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files

182
AMSS_NCKU_source/z4c_rhs_c.C
View File

@@ -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);