From 0e6f632df2e9e0a28e0c938bccfae19f4ad5efaa Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 16 Aug 2019 10:50:29 +0200 Subject: Compute Kdot ourselves, do not rely on ML_BSSN to do it. --- src/qms.c | 57 ++++++++++++++++++++++++--------------------------------- 1 file 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); -- cgit v1.2.3