Compute Kdot ourselves, do not rely on ML_BSSN to do it.

1 files changed, 24 insertions, 33 deletions

diff --git a/src/qms.c b/src/qms.c
index 7fbde0b..cfad019 100644
--- a/src/qms.c
+++ b/src/qms.c
@@ -371,12 +371,6 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
     double  *a_phi   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::phi");
     double  *a_Xtx   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt1");
     double  *a_Xtz   = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt3");
-    double *a_kdotxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot11");
-    double *a_kdotyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot22");
-    double *a_kdotzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot33");
-    double *a_kdotxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot12");
-    double *a_kdotxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot13");
-    double *a_kdotyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Kdot23");
 
     double  *a_alpha = CCTK_VarDataPtr(gh, 0, "ML_BSSN::alpha");
 
@@ -402,6 +396,8 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
             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 D2alpha[3][3];
+            double Kdot[3][3];
             double k2, Kij_Dij_alpha, k_kdot, k3;
             double Gammadot_term, beta_term;
 
@@ -587,17 +583,6 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
 
             const double dtrK[3] = { dx_trk, 0.0, dz_trk };
 
-            const double kdot_xx = a_kdotxx[idx_src];
-            const double kdot_yy = a_kdotyy[idx_src];
-            const double kdot_zz = a_kdotzz[idx_src];
-            const double kdot_xy = a_kdotxy[idx_src];
-            const double kdot_xz = a_kdotxz[idx_src];
-            const double kdot_yz = a_kdotyz[idx_src];
-
-            const double kdot[3][3] = {{ kdot_xx, kdot_xy, kdot_xz },
-                                       { kdot_xy, kdot_yy, kdot_yz },
-                                       { kdot_xz, kdot_yz, kdot_zz }};
-
             const double  Xtx  = a_Xtx[idx_src];
             const double  Xtz  = a_Xtz[idx_src];
 
@@ -838,20 +823,20 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
                 Gammadot_term += val * dalpha[j];
             }
 
-            Kij_Dij_alpha = 0.0;
             for (int j = 0; j < 3; j++)
                 for (int k = 0; k < 3; k++) {
                     double val = 0.0;
                     for (int l = 0; l < 3; l++)
                         val += G[l][j][k] * dalpha[l];
 
-                    Kij_Dij_alpha += Ku[j][k] * (d2alpha[j][k] - val);
+                    D2alpha[j][k] = d2alpha[j][k] - val;
                 }
 
-            k_kdot = 0.0;
+            Kij_Dij_alpha = 0.0;
             for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++)
-                    k_kdot += kdot[j][k] * Ku[j][k];
+                for (int k = 0; k < 3; k++) {
+                    Kij_Dij_alpha += Ku[j][k] * D2alpha[j][k];
+                }
 
             k3 = 0.0;
             for (int j = 0; j < 3; j++)
@@ -889,6 +874,23 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve
                     Ricm[j][k] = val;
                 }
 
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    double val = 0.0;
+                    for (int l = 0; l < 3; l++)
+                        val += K[j][l] * Km[l][k];
+
+                    Kdot[j][k] = -D2alpha[j][k] + alpha * (Ric[j][k] + trK * K[j][k] - 2 * val);
+                }
+
+            k_kdot = 0.0;
+            for (int j = 0; j < 3; j++)
+                for (int k = 0; k < 3; k++) {
+                    //k_kdot += kdot[j][k] * Ku[j][k];
+                    k_kdot += Kdot[j][k] * Ku[j][k];
+                }
+
+
             beta_term = 0.0;
             for (int j = 0; j < 3; j++)
                 for (int k = 0; k < 3; k++)
@@ -1264,17 +1266,6 @@ void qms_mg_inithist(CCTK_ARGUMENTS)
         W_val1_time[i]  = DBL_MAX;
     }
 
-    memset(Kdot11, 0, sizeof(double) * grid_size);
-    memset(Kdot22, 0, sizeof(double) * grid_size);
-    memset(Kdot33, 0, sizeof(double) * grid_size);
-    memset(Kdot12, 0, sizeof(double) * grid_size);
-    memset(Kdot13, 0, sizeof(double) * grid_size);
-    memset(Kdot23, 0, sizeof(double) * grid_size);
-
-    memset(Xtdot1, 0, sizeof(double) * grid_size);
-    memset(Xtdot2, 0, sizeof(double) * grid_size);
-    memset(Xtdot3, 0, sizeof(double) * grid_size);
-
     memset(W_pred0, 0, sizeof(double) * grid_size);
     memset(W_pred1, 0, sizeof(double) * grid_size);
     memset(W_val0, 0, sizeof(double) * grid_size);