From 1fbc213888666a2ada3bb70fa4484bcfade98a16 Mon Sep 17 00:00:00 2001
From: CGH0S7 <776459475@qq.com>
Date: Sat, 28 Feb 2026 15:50:52 +0800
Subject: [PATCH] bssn_rhs_c: remove gxx/gyy/gzz temporaries in favor of
 dxx/dyy/dzz+1

---
 AMSS_NCKU_source/bssn_rhs_c.C | 132 +++++++++++++++++-----------------
 1 file changed, 64 insertions(+), 68 deletions(-)

diff --git a/AMSS_NCKU_source/bssn_rhs_c.C b/AMSS_NCKU_source/bssn_rhs_c.C
index 9dabd74..ba0f615 100644
--- a/AMSS_NCKU_source/bssn_rhs_c.C
+++ b/AMSS_NCKU_source/bssn_rhs_c.C
@@ -39,7 +39,6 @@ int f_compute_rhs_bssn(int *ex, double &T,
     // printf("nx=%d ny=%d nz=%d all=%d\n", nx, ny, nz, all);
 
     // temp variable
-    double gxx[all],gyy[all],gzz[all];
     double chix[all],chiy[all],chiz[all];
     double gxxx[all],gxyx[all],gxzx[all],gyyx[all],gyzx[all],gzzx[all];
     double gxxy[all],gxyy[all],gxzy[all],gyyy[all],gyzy[all],gzzy[all];
@@ -107,9 +106,6 @@ int f_compute_rhs_bssn(int *ex, double &T,
         for(int i=0;i<all;i+=1){
             alpn1[i] = Lap[i] + 1.0;
             chin1[i] = chi[i] + 1.0;
-            gxx[i] = dxx[i] + 1.0;
-            gyy[i] = dyy[i] + 1.0;
-            gzz[i] = dzz[i] + 1.0;
         }
         // 9ms //
         fderivs(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev);
@@ -129,32 +125,32 @@ int f_compute_rhs_bssn(int *ex, double &T,
         for(int i=0;i<all;i+=1){
             div_beta[i] = betaxx[i] + betayy[i] + betazz[i];
             chi_rhs[i] = F2o3 * chin1[i] * (alpn1[i] * trK[i] - div_beta[i]); 
-            gxx_rhs[i] = -TWO * alpn1[i] * Axx[i] - F2o3 * gxx[i] * div_beta[i] + 
-                    TWO * (gxx[i] * betaxx[i] + gxy[i] * betayx[i] + gxz[i] * betazx[i]);
-            gyy_rhs[i] = -TWO * alpn1[i] * Ayy[i] - F2o3 * gyy[i] * div_beta[i] + 
-                    TWO * (gxy[i] * betaxy[i] + gyy[i] * betayy[i] + gyz[i] * betazy[i]);
-            gzz_rhs[i] = -TWO * alpn1[i] * Azz[i] - F2o3 * gzz[i] * div_beta[i] + 
-                        TWO * (gxz[i] * betaxz[i] + gyz[i] * betayz[i] + gzz[i] * betazz[i]);
+            gxx_rhs[i] = -TWO * alpn1[i] * Axx[i] - F2o3 * (dxx[i] + ONE) * div_beta[i] + 
+                    TWO * ((dxx[i] + ONE) * betaxx[i] + gxy[i] * betayx[i] + gxz[i] * betazx[i]);
+            gyy_rhs[i] = -TWO * alpn1[i] * Ayy[i] - F2o3 * (dyy[i] + ONE) * div_beta[i] + 
+                    TWO * (gxy[i] * betaxy[i] + (dyy[i] + ONE) * betayy[i] + gyz[i] * betazy[i]);
+            gzz_rhs[i] = -TWO * alpn1[i] * Azz[i] - F2o3 * (dzz[i] + ONE) * div_beta[i] + 
+                        TWO * (gxz[i] * betaxz[i] + gyz[i] * betayz[i] + (dzz[i] + ONE) * betazz[i]);
             gxy_rhs[i] = -TWO * alpn1[i] * Axy[i] + F1o3 * gxy[i] * div_beta[i] + 
-                        gxx[i] * betaxy[i] + gxz[i] * betazy[i] + gyy[i] * betayx[i] 
+                        (dxx[i] + ONE) * betaxy[i] + gxz[i] * betazy[i] + (dyy[i] + ONE) * betayx[i] 
                         + gyz[i] * betazx[i] - gxy[i] * betazz[i];        
             gyz_rhs[i] = -TWO * alpn1[i] * Ayz[i] + F1o3 * gyz[i] * div_beta[i] + 
-                        gxy[i] * betaxz[i] + gyy[i] * betayz[i] + gxz[i] * betaxy[i] 
-                        + gzz[i] * betazy[i] - gyz[i] * betaxx[i];
+                        gxy[i] * betaxz[i] + (dyy[i] + ONE) * betayz[i] + gxz[i] * betaxy[i] 
+                        + (dzz[i] + ONE) * betazy[i] - gyz[i] * betaxx[i];
             gxz_rhs[i] = -TWO * alpn1[i] * Axz[i] + F1o3 * gxz[i] * div_beta[i] + 
-                        gxx[i] * betaxz[i] + gxy[i] * betayz[i] + gyz[i] * betayx[i] 
-                        + gzz[i] * betazx[i] - gxz[i] * betayy[i];                  
+                        (dxx[i] + ONE) * betaxz[i] + gxy[i] * betayz[i] + gyz[i] * betayx[i] 
+                        + (dzz[i] + ONE) * betazx[i] - gxz[i] * betayy[i];                  
         }
         // Fused: inverse metric + Gamma constraint + Christoffel (3 loops -> 1)
         for(int i=0;i<all;i+=1){
-            double det = gxx[i] * gyy[i] * gzz[i] + gxy[i] * gyz[i] * gxz[i] + gxz[i] * gxy[i] * gyz[i] -
-                        gxz[i] * gyy[i] * gxz[i] - gxy[i] * gxy[i] * gzz[i] - gxx[i] * gyz[i] * gyz[i];
-            double lg_xx = (gyy[i] * gzz[i] - gyz[i] * gyz[i]) / det;
-            double lg_xy = -(gxy[i] * gzz[i] - gyz[i] * gxz[i]) / det;
-            double lg_xz = (gxy[i] * gyz[i] - gyy[i] * gxz[i]) / det;
-            double lg_yy = (gxx[i] * gzz[i] - gxz[i] * gxz[i]) / det;
-            double lg_yz = -(gxx[i] * gyz[i] - gxy[i] * gxz[i]) / det;
-            double lg_zz = (gxx[i] * gyy[i] - gxy[i] * gxy[i]) / det;
+            double det = (dxx[i] + ONE) * (dyy[i] + ONE) * (dzz[i] + ONE) + gxy[i] * gyz[i] * gxz[i] + gxz[i] * gxy[i] * gyz[i] -
+                        gxz[i] * (dyy[i] + ONE) * gxz[i] - gxy[i] * gxy[i] * (dzz[i] + ONE) - (dxx[i] + ONE) * gyz[i] * gyz[i];
+            double lg_xx = ((dyy[i] + ONE) * (dzz[i] + ONE) - gyz[i] * gyz[i]) / det;
+            double lg_xy = -(gxy[i] * (dzz[i] + ONE) - gyz[i] * gxz[i]) / det;
+            double lg_xz = (gxy[i] * gyz[i] - (dyy[i] + ONE) * gxz[i]) / det;
+            double lg_yy = ((dxx[i] + ONE) * (dzz[i] + ONE) - gxz[i] * gxz[i]) / det;
+            double lg_yz = -((dxx[i] + ONE) * gyz[i] - gxy[i] * gxz[i]) / det;
+            double lg_zz = ((dxx[i] + ONE) * (dyy[i] + ONE) - gxy[i] * gxy[i]) / det;
             gupxx[i] = lg_xx; gupxy[i] = lg_xy; gupxz[i] = lg_xz;
             gupyy[i] = lg_yy; gupyz[i] = lg_yz; gupzz[i] = lg_zz;
 
@@ -368,26 +364,26 @@ int f_compute_rhs_bssn(int *ex, double &T,
         }
         // 4.4ms //
         for (int i=0;i<all;i+=1) {
-            gxxx[i] = gxx[i]*Gamxxx[i] + gxy[i]*Gamyxx[i] + gxz[i]*Gamzxx[i];
-            gxyx[i] = gxx[i]*Gamxxy[i] + gxy[i]*Gamyxy[i] + gxz[i]*Gamzxy[i];
-            gxzx[i] = gxx[i]*Gamxxz[i] + gxy[i]*Gamyxz[i] + gxz[i]*Gamzxz[i];
-            gyyx[i] = gxx[i]*Gamxyy[i] + gxy[i]*Gamyyy[i] + gxz[i]*Gamzyy[i];
-            gyzx[i] = gxx[i]*Gamxyz[i] + gxy[i]*Gamyyz[i] + gxz[i]*Gamzyz[i];
-            gzzx[i] = gxx[i]*Gamxzz[i] + gxy[i]*Gamyzz[i] + gxz[i]*Gamzzz[i];
+            gxxx[i] = (dxx[i] + ONE)*Gamxxx[i] + gxy[i]*Gamyxx[i] + gxz[i]*Gamzxx[i];
+            gxyx[i] = (dxx[i] + ONE)*Gamxxy[i] + gxy[i]*Gamyxy[i] + gxz[i]*Gamzxy[i];
+            gxzx[i] = (dxx[i] + ONE)*Gamxxz[i] + gxy[i]*Gamyxz[i] + gxz[i]*Gamzxz[i];
+            gyyx[i] = (dxx[i] + ONE)*Gamxyy[i] + gxy[i]*Gamyyy[i] + gxz[i]*Gamzyy[i];
+            gyzx[i] = (dxx[i] + ONE)*Gamxyz[i] + gxy[i]*Gamyyz[i] + gxz[i]*Gamzyz[i];
+            gzzx[i] = (dxx[i] + ONE)*Gamxzz[i] + gxy[i]*Gamyzz[i] + gxz[i]*Gamzzz[i];
 
-            gxxy[i] = gxy[i]*Gamxxx[i] + gyy[i]*Gamyxx[i] + gyz[i]*Gamzxx[i];
-            gxyy[i] = gxy[i]*Gamxxy[i] + gyy[i]*Gamyxy[i] + gyz[i]*Gamzxy[i];
-            gxzy[i] = gxy[i]*Gamxxz[i] + gyy[i]*Gamyxz[i] + gyz[i]*Gamzxz[i];
-            gyyy[i] = gxy[i]*Gamxyy[i] + gyy[i]*Gamyyy[i] + gyz[i]*Gamzyy[i];
-            gyzy[i] = gxy[i]*Gamxyz[i] + gyy[i]*Gamyyz[i] + gyz[i]*Gamzyz[i];
-            gzzy[i] = gxy[i]*Gamxzz[i] + gyy[i]*Gamyzz[i] + gyz[i]*Gamzzz[i];
+            gxxy[i] = gxy[i]*Gamxxx[i] + (dyy[i] + ONE)*Gamyxx[i] + gyz[i]*Gamzxx[i];
+            gxyy[i] = gxy[i]*Gamxxy[i] + (dyy[i] + ONE)*Gamyxy[i] + gyz[i]*Gamzxy[i];
+            gxzy[i] = gxy[i]*Gamxxz[i] + (dyy[i] + ONE)*Gamyxz[i] + gyz[i]*Gamzxz[i];
+            gyyy[i] = gxy[i]*Gamxyy[i] + (dyy[i] + ONE)*Gamyyy[i] + gyz[i]*Gamzyy[i];
+            gyzy[i] = gxy[i]*Gamxyz[i] + (dyy[i] + ONE)*Gamyyz[i] + gyz[i]*Gamzyz[i];
+            gzzy[i] = gxy[i]*Gamxzz[i] + (dyy[i] + ONE)*Gamyzz[i] + gyz[i]*Gamzzz[i];
 
-            gxxz[i] = gxz[i]*Gamxxx[i] + gyz[i]*Gamyxx[i] + gzz[i]*Gamzxx[i];
-            gxyz[i] = gxz[i]*Gamxxy[i] + gyz[i]*Gamyxy[i] + gzz[i]*Gamzxy[i];
-            gxzz[i] = gxz[i]*Gamxxz[i] + gyz[i]*Gamyxz[i] + gzz[i]*Gamzxz[i];
-            gyyz[i] = gxz[i]*Gamxyy[i] + gyz[i]*Gamyyy[i] + gzz[i]*Gamzyy[i];
-            gyzz[i] = gxz[i]*Gamxyz[i] + gyz[i]*Gamyyz[i] + gzz[i]*Gamzyz[i];
-            gzzz[i] = gxz[i]*Gamxzz[i] + gyz[i]*Gamyzz[i] + gzz[i]*Gamzzz[i];
+            gxxz[i] = gxz[i]*Gamxxx[i] + gyz[i]*Gamyxx[i] + (dzz[i] + ONE)*Gamzxx[i];
+            gxyz[i] = gxz[i]*Gamxxy[i] + gyz[i]*Gamyxy[i] + (dzz[i] + ONE)*Gamzxy[i];
+            gxzz[i] = gxz[i]*Gamxxz[i] + gyz[i]*Gamyxz[i] + (dzz[i] + ONE)*Gamzxz[i];
+            gyyz[i] = gxz[i]*Gamxyy[i] + gyz[i]*Gamyyy[i] + (dzz[i] + ONE)*Gamzyy[i];
+            gyzz[i] = gxz[i]*Gamxyz[i] + gyz[i]*Gamyyz[i] + (dzz[i] + ONE)*Gamzyz[i];
+            gzzz[i] = gxz[i]*Gamxzz[i] + gyz[i]*Gamyzz[i] + (dzz[i] + ONE)*Gamzzz[i];
         }
         // 22.2ms //
             fdderivs(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev);
@@ -437,7 +433,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
         for (int i = 0; i < all; i += 1) {
             Rxx[i] =
                 -HALF * Rxx[i]
-                + gxx[i] * Gamxx[i] + gxy[i] * Gamyx[i] + gxz[i] * Gamzx[i]
+                + (dxx[i] + ONE) * Gamxx[i] + gxy[i] * Gamyx[i] + gxz[i] * Gamzx[i]
                 + Gamxa[i] * gxxx[i] + Gamya[i] * gxyx[i] + Gamza[i] * gxzx[i]
                 + gupxx[i] * (
                     TWO * (Gamxxx[i] * gxxx[i] + Gamyxx[i] * gxyx[i] + Gamzxx[i] * gxzx[i]) +
@@ -472,7 +468,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
 
             Ryy[i] =
                 -HALF * Ryy[i]
-                + gxy[i] * Gamxy[i] + gyy[i] * Gamyy[i] + gyz[i] * Gamzy[i]
+                + gxy[i] * Gamxy[i] + (dyy[i] + ONE) * Gamyy[i] + gyz[i] * Gamzy[i]
                 + Gamxa[i] * gxyy[i] + Gamya[i] * gyyy[i] + Gamza[i] * gyzy[i]
                 + gupxx[i] * (
                     TWO * (Gamxxy[i] * gxxy[i] + Gamyxy[i] * gxyy[i] + Gamzxy[i] * gxzy[i]) +
@@ -507,7 +503,7 @@ int f_compute_rhs_bssn(int *ex, double &T,
 
             Rzz[i] =
                 -HALF * Rzz[i]
-                + gxz[i] * Gamxz[i] + gyz[i] * Gamyz[i] + gzz[i] * Gamzz[i]
+                + gxz[i] * Gamxz[i] + gyz[i] * Gamyz[i] + (dzz[i] + ONE) * Gamzz[i]
                 + Gamxa[i] * gxzz[i] + Gamya[i] * gyzz[i] + Gamza[i] * gzzz[i]
                 + gupxx[i] * (
                     TWO * (Gamxxz[i] * gxxz[i] + Gamyxz[i] * gxyz[i] + Gamzxz[i] * gxzz[i]) +
@@ -543,8 +539,8 @@ int f_compute_rhs_bssn(int *ex, double &T,
             Rxy[i] =
                 HALF * (
                     -Rxy[i]
-                    + gxx[i] * Gamxy[i] + gxy[i] * Gamyy[i] + gxz[i] * Gamzy[i]
-                    + gxy[i] * Gamxx[i] + gyy[i] * Gamyx[i] + gyz[i] * Gamzx[i]
+                    + (dxx[i] + ONE) * Gamxy[i] + gxy[i] * Gamyy[i] + gxz[i] * Gamzy[i]
+                    + gxy[i] * Gamxx[i] + (dyy[i] + ONE) * Gamyx[i] + gyz[i] * Gamzx[i]
                     + Gamxa[i] * gxyx[i] + Gamya[i] * gyyx[i] + Gamza[i] * gyzx[i]
                     + Gamxa[i] * gxxy[i] + Gamya[i] * gxyy[i] + Gamza[i] * gxzy[i]
                 )
@@ -591,8 +587,8 @@ int f_compute_rhs_bssn(int *ex, double &T,
             Rxz[i] =
                 HALF * (
                     -Rxz[i]
-                    + gxx[i] * Gamxz[i] + gxy[i] * Gamyz[i] + gxz[i] * Gamzz[i]
-                    + gxz[i] * Gamxx[i] + gyz[i] * Gamyx[i] + gzz[i] * Gamzx[i]
+                    + (dxx[i] + ONE) * Gamxz[i] + gxy[i] * Gamyz[i] + gxz[i] * Gamzz[i]
+                    + gxz[i] * Gamxx[i] + gyz[i] * Gamyx[i] + (dzz[i] + ONE) * Gamzx[i]
                     + Gamxa[i] * gxzx[i] + Gamya[i] * gyzx[i] + Gamza[i] * gzzx[i]
                     + Gamxa[i] * gxxz[i] + Gamya[i] * gxyz[i] + Gamza[i] * gxzz[i]
                 )
@@ -639,8 +635,8 @@ int f_compute_rhs_bssn(int *ex, double &T,
             Ryz[i] =
                 HALF * (
                     -Ryz[i]
-                    + gxy[i] * Gamxz[i] + gyy[i] * Gamyz[i] + gyz[i] * Gamzz[i]
-                    + gxz[i] * Gamxy[i] + gyz[i] * Gamyy[i] + gzz[i] * Gamzy[i]
+                    + gxy[i] * Gamxz[i] + (dyy[i] + ONE) * Gamyz[i] + gyz[i] * Gamzz[i]
+                    + gxz[i] * Gamxy[i] + gyz[i] * Gamyy[i] + (dzz[i] + ONE) * Gamzy[i]
                     + Gamxa[i] * gxzy[i] + Gamya[i] * gyzy[i] + Gamza[i] * gzzy[i]
                     + Gamxa[i] * gxyz[i] + Gamya[i] * gyyz[i] + Gamza[i] * gyzz[i]
                 )
@@ -703,9 +699,9 @@ int f_compute_rhs_bssn(int *ex, double &T,
                 + 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) + gxx[i] * f[i] ) / (chin1[i] * TWO);
-            Ryy[i] = Ryy[i] + ( fyy[i] - (chiy[i] * chiy[i]) / (chin1[i] * TWO) + gyy[i] * f[i] ) / (chin1[i] * TWO);
-            Rzz[i] = Rzz[i] + ( fzz[i] - (chiz[i] * chiz[i]) / (chin1[i] * TWO) + gzz[i] * f[i] ) / (chin1[i] * TWO);
+            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);
 
             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);
@@ -724,17 +720,17 @@ int f_compute_rhs_bssn(int *ex, double &T,
             gxxz[i] = (gupxz[i] * chix[i] + gupyz[i] * chiy[i] + gupzz[i] * chiz[i]) / chin1[i];
 
             /* Christoffel 修正项 */
-            Gamxxx[i] = Gamxxx[i] - ( ((chix[i] + chix[i]) / chin1[i]) - gxx[i] * gxxx[i] ) * HALF;
-            Gamyxx[i] = Gamyxx[i] - ( 0.0           - gxx[i] * gxxy[i] ) * HALF;  /* 原式只有 -gxx*gxxy */
-            Gamzxx[i] = Gamzxx[i] - ( 0.0           - gxx[i] * gxxz[i] ) * HALF;
+            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;
 
-            Gamxyy[i] = Gamxyy[i] - ( 0.0           - gyy[i] * gxxx[i] ) * HALF;
-            Gamyyy[i] = Gamyyy[i] - ( ((chiy[i] + chiy[i]) / chin1[i]) - gyy[i] * gxxy[i] ) * HALF;
-            Gamzyy[i] = Gamzyy[i] - ( 0.0           - gyy[i] * gxxz[i] ) * 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;
 
-            Gamxzz[i] = Gamxzz[i] - ( 0.0           - gzz[i] * gxxx[i] ) * HALF;
-            Gamyzz[i] = Gamyzz[i] - ( 0.0           - gzz[i] * gxxy[i] ) * HALF;
-            Gamzzz[i] = Gamzzz[i] - ( ((chiz[i] + chiz[i]) / chin1[i]) - gzz[i] * gxxz[i] ) * 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;
 
             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;
@@ -820,9 +816,9 @@ int f_compute_rhs_bssn(int *ex, double &T,
             double l_fzz = alpn1[i] * (Rzz[i] - EIGHT * PI * Szz[i]) - fzz[i];
 
             /* Aij_rhs = fij - gij * f */
-            Axx_rhs[i] = l_fxx - gxx[i] * f[i];
-            Ayy_rhs[i] = l_fyy - gyy[i] * f[i];
-            Azz_rhs[i] = l_fzz - gzz[i] * f[i];
+            Axx_rhs[i] = l_fxx - (dxx[i] + ONE) * f[i];
+            Ayy_rhs[i] = l_fyy - (dyy[i] + ONE) * f[i];
+            Azz_rhs[i] = l_fzz - (dzz[i] + ONE) * f[i];
             Axy_rhs[i] = l_fxy - gxy[i] * f[i];
             Axz_rhs[i] = l_fxz - gxz[i] * f[i];
             Ayz_rhs[i] = l_fyz - gyz[i] * f[i];
@@ -1060,17 +1056,17 @@ int f_compute_rhs_bssn(int *ex, double &T,
             #endif
         }
         // advection + KO dissipation with shared symmetry buffer
-        lopsided_kodis(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps);
+        lopsided_kodis(ex,X,Y,Z,dxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps);
         lopsided_kodis(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA,eps);
         lopsided_kodis(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,eps);
         lopsided_kodis(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS,eps);
         lopsided_kodis(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,eps);
         lopsided_kodis(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS,eps);
-        lopsided_kodis(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps);
+        lopsided_kodis(ex,X,Y,Z,dyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps);
         lopsided_kodis(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS,eps);
         lopsided_kodis(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,eps);
         lopsided_kodis(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA,eps);
-        lopsided_kodis(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps);
+        lopsided_kodis(ex,X,Y,Z,dzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps);
         lopsided_kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS,eps);
         lopsided_kodis(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,eps);
         lopsided_kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS,eps);