From 65e0f95f40651f4a9e6d3ffc8ec39cdd816dce6b Mon Sep 17 00:00:00 2001
From: CGH0S7 <776459475@qq.com>
Date: Mon, 13 Apr 2026 15:14:31 +0800
Subject: [PATCH] Localize chi Ricci intermediates in RHS

---
 AMSS_NCKU_source/bssn_rhs_c.C | 86 ++++++++++++++++++-----------------
 1 file changed, 45 insertions(+), 41 deletions(-)

diff --git a/AMSS_NCKU_source/bssn_rhs_c.C b/AMSS_NCKU_source/bssn_rhs_c.C
index 3eaf764..55dd118 100644
--- a/AMSS_NCKU_source/bssn_rhs_c.C
+++ b/AMSS_NCKU_source/bssn_rhs_c.C
@@ -688,26 +688,30 @@ int f_compute_rhs_bssn(int *ex, double &T,
 
         // 7ms //
         for (int i=0;i<all;i+=1) {
-            fxx[i] = fxx[i] - Gamxxx[i] * chix[i] - Gamyxx[i] * chiy[i] - Gamzxx[i] * chiz[i];
-            fxy[i] = fxy[i] - Gamxxy[i] * chix[i] - Gamyxy[i] * chiy[i] - Gamzxy[i] * chiz[i];
-            fxz[i] = fxz[i] - Gamxxz[i] * chix[i] - Gamyxz[i] * chiy[i] - Gamzxz[i] * chiz[i];
-            fyy[i] = fyy[i] - Gamxyy[i] * chix[i] - Gamyyy[i] * chiy[i] - Gamzyy[i] * chiz[i];
-            fyz[i] = fyz[i] - Gamxyz[i] * chix[i] - Gamyyz[i] * chiy[i] - Gamzyz[i] * chiz[i];
-            fzz[i] = fzz[i] - Gamxzz[i] * chix[i] - Gamyzz[i] * chiy[i] - Gamzzz[i] * chiz[i];
-            f[i] =
-                gupxx[i] * (fxx[i] - (F3o2 / chin1[i]) * chix[i] * chix[i])
-                + gupyy[i] * (fyy[i] - (F3o2 / chin1[i]) * chiy[i] * chiy[i])
-                + gupzz[i] * (fzz[i] - (F3o2 / chin1[i]) * chiz[i] * chiz[i])
-                + TWO * gupxy[i] * (fxy[i] - (F3o2 / chin1[i]) * chix[i] * chiy[i])
-                + TWO * gupxz[i] * (fxz[i] - (F3o2 / chin1[i]) * chix[i] * chiz[i])
-                + TWO * gupyz[i] * (fyz[i] - (F3o2 / chin1[i]) * chiy[i] * chiz[i]);
-            Rxx[i] = Rxx[i] + ( fxx[i] - (chix[i] * chix[i]) / (chin1[i] * TWO) + (dxx[i] + ONE) * f[i] ) / (chin1[i] * TWO);
-            Ryy[i] = Ryy[i] + ( fyy[i] - (chiy[i] * chiy[i]) / (chin1[i] * TWO) + (dyy[i] + ONE) * f[i] ) / (chin1[i] * TWO);
-            Rzz[i] = Rzz[i] + ( fzz[i] - (chiz[i] * chiz[i]) / (chin1[i] * TWO) + (dzz[i] + ONE) * f[i] ) / (chin1[i] * TWO);
+            const double inv_chin1 = ONE / chin1[i];
+            const double half_inv_chin1 = HALF * inv_chin1;
+            const double scaled_inv = F3o2 * inv_chin1;
+            const double cxx = fxx[i] - Gamxxx[i] * chix[i] - Gamyxx[i] * chiy[i] - Gamzxx[i] * chiz[i];
+            const double cxy = fxy[i] - Gamxxy[i] * chix[i] - Gamyxy[i] * chiy[i] - Gamzxy[i] * chiz[i];
+            const double cxz = fxz[i] - Gamxxz[i] * chix[i] - Gamyxz[i] * chiy[i] - Gamzxz[i] * chiz[i];
+            const double cyy = fyy[i] - Gamxyy[i] * chix[i] - Gamyyy[i] * chiy[i] - Gamzyy[i] * chiz[i];
+            const double cyz = fyz[i] - Gamxyz[i] * chix[i] - Gamyyz[i] * chiy[i] - Gamzyz[i] * chiz[i];
+            const double czz = fzz[i] - Gamxzz[i] * chix[i] - Gamyzz[i] * chiy[i] - Gamzzz[i] * chiz[i];
+            const double ricci_chi =
+                gupxx[i] * (cxx - scaled_inv * chix[i] * chix[i])
+                + gupyy[i] * (cyy - scaled_inv * chiy[i] * chiy[i])
+                + gupzz[i] * (czz - scaled_inv * chiz[i] * chiz[i])
+                + TWO * gupxy[i] * (cxy - scaled_inv * chix[i] * chiy[i])
+                + TWO * gupxz[i] * (cxz - scaled_inv * chix[i] * chiz[i])
+                + TWO * gupyz[i] * (cyz - scaled_inv * chiy[i] * chiz[i]);
+            f[i] = ricci_chi;
+            Rxx[i] = Rxx[i] + ( cxx - half_inv_chin1 * chix[i] * chix[i] + (dxx[i] + ONE) * ricci_chi ) * half_inv_chin1;
+            Ryy[i] = Ryy[i] + ( cyy - half_inv_chin1 * chiy[i] * chiy[i] + (dyy[i] + ONE) * ricci_chi ) * half_inv_chin1;
+            Rzz[i] = Rzz[i] + ( czz - half_inv_chin1 * chiz[i] * chiz[i] + (dzz[i] + ONE) * ricci_chi ) * half_inv_chin1;
 
-            Rxy[i] = Rxy[i] + ( fxy[i] - (chix[i] * chiy[i]) / (chin1[i] * TWO) + gxy[i] * f[i] ) / (chin1[i] * TWO);
-            Rxz[i] = Rxz[i] + ( fxz[i] - (chix[i] * chiz[i]) / (chin1[i] * TWO) + gxz[i] * f[i] ) / (chin1[i] * TWO);
-            Ryz[i] = Ryz[i] + ( fyz[i] - (chiy[i] * chiz[i]) / (chin1[i] * TWO) + gyz[i] * f[i] ) / (chin1[i] * TWO);
+            Rxy[i] = Rxy[i] + ( cxy - half_inv_chin1 * chix[i] * chiy[i] + gxy[i] * ricci_chi ) * half_inv_chin1;
+            Rxz[i] = Rxz[i] + ( cxz - half_inv_chin1 * chix[i] * chiz[i] + gxz[i] * ricci_chi ) * half_inv_chin1;
+            Ryz[i] = Ryz[i] + ( cyz - half_inv_chin1 * chiy[i] * chiz[i] + gyz[i] * ricci_chi ) * half_inv_chin1;
         }
 
         // 24ms //
@@ -716,35 +720,35 @@ int f_compute_rhs_bssn(int *ex, double &T,
 
         // 6ms //
         for (int i=0;i<all;i+=1) {
-            /* gxxx,gxxy,gxxz (这里是“升指标后的chi导数/chi”那类量，你沿用原变量名即可) */
-            gxxx[i] = (gupxx[i] * chix[i] + gupxy[i] * chiy[i] + gupxz[i] * chiz[i]) / chin1[i];
-            gxxy[i] = (gupxy[i] * chix[i] + gupyy[i] * chiy[i] + gupyz[i] * chiz[i]) / chin1[i];
-            gxxz[i] = (gupxz[i] * chix[i] + gupyz[i] * chiy[i] + gupzz[i] * chiz[i]) / chin1[i];
+            const double inv_chin1 = ONE / chin1[i];
+            const double gchi_x = (gupxx[i] * chix[i] + gupxy[i] * chiy[i] + gupxz[i] * chiz[i]) * inv_chin1;
+            const double gchi_y = (gupxy[i] * chix[i] + gupyy[i] * chiy[i] + gupyz[i] * chiz[i]) * inv_chin1;
+            const double gchi_z = (gupxz[i] * chix[i] + gupyz[i] * chiy[i] + gupzz[i] * chiz[i]) * inv_chin1;
 
             /* Christoffel 修正项 */
-            Gamxxx[i] = Gamxxx[i] - ( ((chix[i] + chix[i]) / chin1[i]) - (dxx[i] + ONE) * gxxx[i] ) * HALF;
-            Gamyxx[i] = Gamyxx[i] - ( 0.0           - (dxx[i] + ONE) * gxxy[i] ) * HALF;  /* 原式只有 -gxx*gxxy */
-            Gamzxx[i] = Gamzxx[i] - ( 0.0           - (dxx[i] + ONE) * gxxz[i] ) * HALF;
+            Gamxxx[i] = Gamxxx[i] - ( ((chix[i] + chix[i]) * inv_chin1) - (dxx[i] + ONE) * gchi_x ) * HALF;
+            Gamyxx[i] = Gamyxx[i] - ( 0.0           - (dxx[i] + ONE) * gchi_y ) * HALF;  /* 原式只有 -gxx*gxxy */
+            Gamzxx[i] = Gamzxx[i] - ( 0.0           - (dxx[i] + ONE) * gchi_z ) * HALF;
 
-            Gamxyy[i] = Gamxyy[i] - ( 0.0           - (dyy[i] + ONE) * gxxx[i] ) * HALF;
-            Gamyyy[i] = Gamyyy[i] - ( ((chiy[i] + chiy[i]) / chin1[i]) - (dyy[i] + ONE) * gxxy[i] ) * HALF;
-            Gamzyy[i] = Gamzyy[i] - ( 0.0           - (dyy[i] + ONE) * gxxz[i] ) * HALF;
+            Gamxyy[i] = Gamxyy[i] - ( 0.0           - (dyy[i] + ONE) * gchi_x ) * HALF;
+            Gamyyy[i] = Gamyyy[i] - ( ((chiy[i] + chiy[i]) * inv_chin1) - (dyy[i] + ONE) * gchi_y ) * HALF;
+            Gamzyy[i] = Gamzyy[i] - ( 0.0           - (dyy[i] + ONE) * gchi_z ) * HALF;
 
-            Gamxzz[i] = Gamxzz[i] - ( 0.0           - (dzz[i] + ONE) * gxxx[i] ) * HALF;
-            Gamyzz[i] = Gamyzz[i] - ( 0.0           - (dzz[i] + ONE) * gxxy[i] ) * HALF;
-            Gamzzz[i] = Gamzzz[i] - ( ((chiz[i] + chiz[i]) / chin1[i]) - (dzz[i] + ONE) * gxxz[i] ) * HALF;
+            Gamxzz[i] = Gamxzz[i] - ( 0.0           - (dzz[i] + ONE) * gchi_x ) * HALF;
+            Gamyzz[i] = Gamyzz[i] - ( 0.0           - (dzz[i] + ONE) * gchi_y ) * HALF;
+            Gamzzz[i] = Gamzzz[i] - ( ((chiz[i] + chiz[i]) * inv_chin1) - (dzz[i] + ONE) * gchi_z ) * HALF;
 
-            Gamxxy[i] = Gamxxy[i] - ( ( chiy[i] / chin1[i])      - gxy[i] * gxxx[i] ) * HALF;
-            Gamyxy[i] = Gamyxy[i] - ( ( chix[i] / chin1[i])      - gxy[i] * gxxy[i] ) * HALF;
-            Gamzxy[i] = Gamzxy[i] - ( 0.0           - gxy[i] * gxxz[i] ) * HALF;
+            Gamxxy[i] = Gamxxy[i] - ( ( chiy[i] * inv_chin1)      - gxy[i] * gchi_x ) * HALF;
+            Gamyxy[i] = Gamyxy[i] - ( ( chix[i] * inv_chin1)      - gxy[i] * gchi_y ) * HALF;
+            Gamzxy[i] = Gamzxy[i] - ( 0.0           - gxy[i] * gchi_z ) * HALF;
 
-            Gamxxz[i] = Gamxxz[i] - ( ( chiz[i] / chin1[i])      - gxz[i] * gxxx[i] ) * HALF;
-            Gamyxz[i] = Gamyxz[i] - ( 0.0           - gxz[i] * gxxy[i] ) * HALF;
-            Gamzxz[i] = Gamzxz[i] - ( ( chix[i] / chin1[i])      - gxz[i] * gxxz[i] ) * HALF;
+            Gamxxz[i] = Gamxxz[i] - ( ( chiz[i] * inv_chin1)      - gxz[i] * gchi_x ) * HALF;
+            Gamyxz[i] = Gamyxz[i] - ( 0.0           - gxz[i] * gchi_y ) * HALF;
+            Gamzxz[i] = Gamzxz[i] - ( ( chix[i] * inv_chin1)      - gxz[i] * gchi_z ) * HALF;
 
-            Gamxyz[i] = Gamxyz[i] - ( 0.0           - gyz[i] * gxxx[i] ) * HALF;
-            Gamyyz[i] = Gamyyz[i] - ( ( chiz[i] / chin1[i])      - gyz[i] * gxxy[i] ) * HALF;
-            Gamzyz[i] = Gamzyz[i] - ( ( chiy[i] / chin1[i])      - gyz[i] * gxxz[i] ) * HALF;
+            Gamxyz[i] = Gamxyz[i] - ( 0.0           - gyz[i] * gchi_x ) * HALF;
+            Gamyyz[i] = Gamyyz[i] - ( ( chiz[i] * inv_chin1)      - gyz[i] * gchi_y ) * HALF;
+            Gamzyz[i] = Gamzyz[i] - ( ( chiy[i] * inv_chin1)      - gyz[i] * gchi_z ) * HALF;
 
             /* fxx..fyz 修正：减去 Γ * ∂Lap */
             fxx[i] = fxx[i] - Gamxxx[i] * Lapx[i] - Gamyxx[i] * Lapy[i] - Gamzxx[i] * Lapz[i];