solve: drop unnecessary variables

1 files changed, 21 insertions, 86 deletions

diff --git a/src/qms_solve.c b/src/qms_solve.c
index c84c1d8..0d74300 100644
--- a/src/qms_solve.c
+++ b/src/qms_solve.c
@@ -57,9 +57,6 @@ enum MetricVars {
     ATXZ,
     ATYZ,
     K,
-    XTX,
-    XTY,
-    XTZ,
     BETAX,
     BETAY,
     BETAZ,
@@ -70,9 +67,6 @@ enum MetricVars {
     KDOT_XY,
     KDOT_XZ,
     KDOT_YZ,
-    XTDOT_X,
-    XTDOT_Y,
-    XTDOT_Z,
     NB_METRIC_VARS,
 };
 
@@ -106,7 +100,6 @@ enum InterpMetricVars {
     I_GTXZ_DZZ,
     I_PHI,
     I_PHI_DX,
-    I_PHI_DY,
     I_PHI_DZ,
     I_PHI_DXX,
     I_PHI_DZZ,
@@ -125,12 +118,9 @@ enum InterpMetricVars {
     I_ATYY_DZ,
     I_ATZZ_DZ,
     I_ATXZ_DZ,
-    I_K,
+    I_TRK,
     I_TRK_DX,
     I_TRK_DZ,
-    I_XTX,
-    I_XTY,
-    I_XTZ,
     I_BETAX,
     I_BETAY,
     I_BETAZ,
@@ -146,23 +136,16 @@ enum InterpMetricVars {
     I_BETAZ_DZZ,
     I_ALPHA,
     I_ALPHA_DX,
-    I_ALPHA_DY,
     I_ALPHA_DZ,
     I_ALPHA_DXX,
-    I_ALPHA_DYY,
     I_ALPHA_DZZ,
-    I_ALPHA_DXY,
     I_ALPHA_DXZ,
-    I_ALPHA_DYZ,
     I_KDOT_XX,
     I_KDOT_YY,
     I_KDOT_ZZ,
     I_KDOT_XY,
     I_KDOT_XZ,
     I_KDOT_YZ,
-    I_XTDOT_X,
-    I_XTDOT_Y,
-    I_XTDOT_Z,
     NB_INTERP_VARS,
 };
 
@@ -221,9 +204,6 @@ static const char *metric_vars[] = {
     [ATYZ]  = "ML_CCZ4::At23",
     [PHI]   = "ML_CCZ4::phi",
     [K]     = "ML_CCZ4::trK",
-    [XTX]   = "ML_CCZ4::Xt1",
-    [XTY]   = "ML_CCZ4::Xt2",
-    [XTZ]   = "ML_CCZ4::Xt3",
     [BETAX] = "ML_CCZ4::beta1",
     [BETAY] = "ML_CCZ4::beta2",
     [BETAZ] = "ML_CCZ4::beta3",
@@ -234,9 +214,6 @@ static const char *metric_vars[] = {
     [KDOT_XY] = "ML_CCZ4::Kdot12",
     [KDOT_XZ] = "ML_CCZ4::Kdot13",
     [KDOT_YZ] = "ML_CCZ4::Kdot23",
-    [XTDOT_X] = "ML_CCZ4::Xtdot1",
-    [XTDOT_Y] = "ML_CCZ4::Xtdot2",
-    [XTDOT_Z] = "ML_CCZ4::Xtdot3",
 #else
     [GTXX]  = "ML_BSSN::gt11",
     [GTYY]  = "ML_BSSN::gt22",
@@ -252,9 +229,6 @@ static const char *metric_vars[] = {
     [ATYZ]  = "ML_BSSN::At23",
     [PHI]   = "ML_BSSN::phi",
     [K]     = "ML_BSSN::trK",
-    [XTX]   = "ML_BSSN::Xt1",
-    [XTY]   = "ML_BSSN::Xt2",
-    [XTZ]   = "ML_BSSN::Xt3",
     [BETAX] = "ML_BSSN::beta1",
     [BETAY] = "ML_BSSN::beta2",
     [BETAZ] = "ML_BSSN::beta3",
@@ -266,9 +240,6 @@ static const char *metric_vars[] = {
     [KDOT_XY] = "ML_BSSN::Kdot12",
     [KDOT_XZ] = "ML_BSSN::Kdot13",
     [KDOT_YZ] = "ML_BSSN::Kdot23",
-    [XTDOT_X] = "ML_BSSN::Xtdot1",
-    [XTDOT_Y] = "ML_BSSN::Xtdot2",
-    [XTDOT_Z] = "ML_BSSN::Xtdot3",
 #endif
 };
 
@@ -302,7 +273,6 @@ static const CCTK_INT interp_operation_indices[] = {
     [I_GTXZ_DZZ] = GTXZ,
     [I_PHI]      = PHI,
     [I_PHI_DX]   = PHI,
-    [I_PHI_DY]   = PHI,
     [I_PHI_DZ]   = PHI,
     [I_PHI_DXX]  = PHI,
     [I_PHI_DZZ]  = PHI,
@@ -321,12 +291,9 @@ static const CCTK_INT interp_operation_indices[] = {
     [I_ATYY_DZ]  = ATYY,
     [I_ATZZ_DZ]  = ATZZ,
     [I_ATXZ_DZ]  = ATXZ,
-    [I_K]        = K,
+    [I_TRK]      = K,
     [I_TRK_DX]   = K,
     [I_TRK_DZ]   = K,
-    [I_XTX]      = XTX,
-    [I_XTY]      = XTY,
-    [I_XTZ]      = XTZ,
     [I_BETAX]     = BETAX,
     [I_BETAY]     = BETAY,
     [I_BETAZ]     = BETAZ,
@@ -342,23 +309,16 @@ static const CCTK_INT interp_operation_indices[] = {
     [I_BETAZ_DZZ] = BETAZ,
     [I_ALPHA]    = ALPHA,
     [I_ALPHA_DX] = ALPHA,
-    [I_ALPHA_DY] = ALPHA,
     [I_ALPHA_DZ] = ALPHA,
     [I_ALPHA_DXX] = ALPHA,
-    [I_ALPHA_DYY] = ALPHA,
     [I_ALPHA_DZZ] = ALPHA,
-    [I_ALPHA_DXY] = ALPHA,
     [I_ALPHA_DXZ] = ALPHA,
-    [I_ALPHA_DYZ] = ALPHA,
     [I_KDOT_XX]   = KDOT_XX,
     [I_KDOT_YY]   = KDOT_YY,
     [I_KDOT_ZZ]   = KDOT_ZZ,
     [I_KDOT_XY]   = KDOT_XY,
     [I_KDOT_XZ]   = KDOT_XZ,
     [I_KDOT_YZ]   = KDOT_YZ,
-    [I_XTDOT_X]   = XTDOT_X,
-    [I_XTDOT_Y]   = XTDOT_Y,
-    [I_XTDOT_Z]   = XTDOT_Z,
 };
 
 /* the operation (plain value or x/y/z-derivative) to apply during interpolation */
@@ -391,7 +351,6 @@ static const CCTK_INT interp_operation_codes[] = {
     [I_GTXZ_DZZ] = 33,
     [I_PHI]      = 0,
     [I_PHI_DX]   = 1,
-    [I_PHI_DY]   = 2,
     [I_PHI_DZ]   = 3,
     [I_PHI_DXX]  = 11,
     [I_PHI_DZZ]  = 33,
@@ -410,12 +369,9 @@ static const CCTK_INT interp_operation_codes[] = {
     [I_ATYY_DZ]  = 3,
     [I_ATZZ_DZ]  = 3,
     [I_ATXZ_DZ]  = 3,
-    [I_K]        = 0,
+    [I_TRK]      = 0,
     [I_TRK_DX]   = 1,
     [I_TRK_DZ]   = 3,
-    [I_XTX]      = 0,
-    [I_XTY]      = 0,
-    [I_XTZ]      = 0,
     [I_BETAX]     = 0,
     [I_BETAY]     = 0,
     [I_BETAZ]     = 0,
@@ -431,23 +387,16 @@ static const CCTK_INT interp_operation_codes[] = {
     [I_BETAZ_DZZ] = 33,
     [I_ALPHA]    = 0,
     [I_ALPHA_DX] = 1,
-    [I_ALPHA_DY] = 2,
     [I_ALPHA_DZ] = 3,
     [I_ALPHA_DXX] = 11,
-    [I_ALPHA_DYY] = 22,
     [I_ALPHA_DZZ] = 33,
-    [I_ALPHA_DXY] = 12,
     [I_ALPHA_DXZ] = 13,
-    [I_ALPHA_DYZ] = 23,
     [I_KDOT_XX]   = 0,
     [I_KDOT_YY]   = 0,
     [I_KDOT_ZZ]   = 0,
     [I_KDOT_XY]   = 0,
     [I_KDOT_XZ]   = 0,
     [I_KDOT_YZ]   = 0,
-    [I_XTDOT_X]   = 0,
-    [I_XTDOT_Y]   = 0,
-    [I_XTDOT_Z]   = 0,
 };
 
 /* interpolate the cactus gridfunctions onto the pseudospectral grid */
@@ -480,7 +429,7 @@ static int calc_eq_coeffs(QMSSolver *ctx)
         double K[3][3], Km[3][3], Ku[3][3], dK[3][3][3];
         double gtu[3][3], g[3][3], gu[3][3];
         double dg[3][3][3], dgu[3][3][3], d2g[3][3][3][3], gudot[3][3];
-        double G[3][3][3], dG[3][3][3][3], Gdot[3][3][3];
+        double G[3][3][3], dG[3][3][3][3], Gdot[3][3][3], X[3];
         double Ric[3][3], Ricm[3][3];
         double dgdot[3][3][3];
         double k2, Kij_Dij_alpha, k_kdot, k3;
@@ -660,7 +609,7 @@ static int calc_eq_coeffs(QMSSolver *ctx)
             { phi_dxz,    0.0,  phi_dzz },
         };
 
-        const double trK     = s->interp_values[I_K][i];
+        const double trK     = s->interp_values[I_TRK][i];
 
         const double dx_trK = s->interp_values[I_TRK_DX][i];
         const double dz_trK = s->interp_values[I_TRK_DZ][i];
@@ -733,12 +682,6 @@ static int calc_eq_coeffs(QMSSolver *ctx)
             },
         };
 
-        const double Xtx  = s->interp_values[I_XTX][i];
-        const double Xtz  = s->interp_values[I_XTZ][i];
-
-        const double Xtdot_x  = s->interp_values[I_XTDOT_X][i];
-        const double Xtdot_z  = s->interp_values[I_XTDOT_Z][i];
-
         const double det = gtxx * gtyy * gtzz + 2 * gtxy * gtyz * gtxz - gtzz * SQR(gtxy) - SQR(gtxz) * gtyy - gtxx * SQR(gtyz);
 
         // \tilde{γ}^{ij}
@@ -822,6 +765,14 @@ static int calc_eq_coeffs(QMSSolver *ctx)
                         dG[j][k][l][m] = 0.5 * val;
                     }
 
+        for (int j = 0; j < 3; j++) {
+            double val = 0.0;
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l < 3; l++)
+                    val += gu[k][l] * G[j][k][l];
+            X[j] = val;
+        }
+
         // ∂_j β_k
         for (int j = 0; j < 3; j++)
             for (int k = 0; k < 3; k++) {
@@ -983,31 +934,15 @@ static int calc_eq_coeffs(QMSSolver *ctx)
             for (int k = 0; k < 3; k++)
                 beta_term += -beta[k] * Ricm[j][k] * dalpha[j];
 
-        {
-            const double gtuxx         = gtu[0][0];
-            const double gtuyy         = gtu[1][1];
-            const double gtuzz         = gtu[2][2];
-            const double gtuxz         = gtu[0][2];
-
-            const double phi_dx        = s->interp_values[I_PHI_DX][i];
-            const double phi_dz        = s->interp_values[I_PHI_DZ][i];
-
-            const double Xtx           = s->interp_values[I_XTX][i];
-            const double Xtz           = s->interp_values[I_XTZ][i];
+        s->eq_coeffs[PSSOLVE_DIFF_ORDER_20][i] = gu[0][0] + ((x <= EPS) ? gu[1][1] : 0.0);
+        s->eq_coeffs[PSSOLVE_DIFF_ORDER_02][i] = gu[2][2];
+        s->eq_coeffs[PSSOLVE_DIFF_ORDER_11][i] = 2.0 * gu[0][2];
+        s->eq_coeffs[PSSOLVE_DIFF_ORDER_10][i] = -X[0] + ((x > EPS) ? gu[1][1] / x : 0.0);
+        s->eq_coeffs[PSSOLVE_DIFF_ORDER_01][i] = -X[2];
+        s->eq_coeffs[PSSOLVE_DIFF_ORDER_00][i] = -k2;
 
-            const double Xx = SQR(phi) * (Xtx + (phi_dx * gtuxx + phi_dz * gtuxz) / phi);
-            const double Xz = SQR(phi) * (Xtz + (phi_dx * gtuxz + phi_dz * gtuzz) / phi);
-
-            s->eq_coeffs[PSSOLVE_DIFF_ORDER_20][i] = SQR(phi) * (gtuxx + ((x <= EPS) ? gtuyy : 0.0));
-            s->eq_coeffs[PSSOLVE_DIFF_ORDER_02][i] = SQR(phi) * gtuzz;
-            s->eq_coeffs[PSSOLVE_DIFF_ORDER_11][i] = SQR(phi) * gtuxz * 2;
-            s->eq_coeffs[PSSOLVE_DIFF_ORDER_10][i] = -Xx + ((x > EPS) ? SQR(phi) * gtuyy / x : 0.0);
-            s->eq_coeffs[PSSOLVE_DIFF_ORDER_01][i] = -Xz;
-            s->eq_coeffs[PSSOLVE_DIFF_ORDER_00][i] = -k2;
-
-            s->rhs[i] = -2 * alpha * Kij_Dij_alpha + Gammadot_term +
-                        2 * alpha * (k_kdot + 2 * alpha * k3) + beta_term;
-        }
+        s->rhs[i] = -2 * alpha * Kij_Dij_alpha + Gammadot_term +
+                    2 * alpha * (k_kdot + 2 * alpha * k3) + beta_term;
     }
 
     return 0;