diff options
Diffstat (limited to 'ML_BSSN/src/ML_BSSN_RHS.cc')
| -rw-r--r-- | ML_BSSN/src/ML_BSSN_RHS.cc | 175 |
1 files changed, 68 insertions, 107 deletions
diff --git a/ML_BSSN/src/ML_BSSN_RHS.cc b/ML_BSSN/src/ML_BSSN_RHS.cc index 28c80ea..098f585 100644 --- a/ML_BSSN/src/ML_BSSN_RHS.cc +++ b/ML_BSSN/src/ML_BSSN_RHS.cc @@ -27,6 +27,11 @@ #define CUB(x) (kmul(x,SQR(x))) #define QAD(x) (SQR(SQR(x))) +extern "C" { +double TMP_fce_GammaDotX(double x, double z,double t,double tau); +double TMP_fce_GammaDotZ(double x, double z,double t,double tau); +} + extern "C" void ML_BSSN_RHS_SelectBCs(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; @@ -53,9 +58,6 @@ extern "C" void ML_BSSN_RHS_SelectBCs(CCTK_ARGUMENTS) ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_metricrhs","flat"); if (ierr < 0) CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_metricrhs."); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_phidot","flat"); - if (ierr < 0) - CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_phidot."); ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_BSSN::ML_shiftrhs","flat"); if (ierr < 0) CCTK_WARN(1, "Failed to register flat BC for ML_BSSN::ML_shiftrhs."); @@ -995,7 +997,6 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx; PDstandardNth12phi = 0; PDstandardNth23phi = 0; - PDstandardNth2trK = 0; PDstandardNth2Xt1 = 0; PDstandardNth2Xt2 = origin ? PDstandardNth1Xt1 : Xt1[index] / xx; PDstandardNth2Xt3 = 0; @@ -2361,10 +2362,6 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co kneg(kmul(IfThen(conformalMethod == 1,kmul(phiL,ToReal(0.333333333333333333333333333333)),ToReal(-0.166666666666666666666666666667)),kadd(JacPDstandardNth1beta1,kadd(JacPDstandardNth2beta2,knmsub(alphaL,trKL,JacPDstandardNth3beta3))))); - CCTK_REAL_VEC phidotL CCTK_ATTRIBUTE_UNUSED = - kneg(kmul(IfThen(conformalMethod == - 1,kmul(phiL,ToReal(0.333333333333333333333333333333)),ToReal(-0.166666666666666666666666666667)),kadd(JacPDstandardNth1beta1,kadd(JacPDstandardNth2beta2,knmsub(alphaL,trKL,JacPDstandardNth3beta3))))); - CCTK_REAL_VEC gt11rhsL CCTK_ATTRIBUTE_UNUSED = kmadd(alphaL,kmul(At11L,ToReal(-2)),kmadd(gt11L,kmul(kadd(JacPDstandardNth1beta1,kadd(JacPDstandardNth2beta2,JacPDstandardNth3beta3)),ToReal(-0.666666666666666666666666666667)),kmul(kmadd(gt11L,JacPDstandardNth1beta1,kmadd(gt12L,JacPDstandardNth1beta2,kmul(gt13L,JacPDstandardNth1beta3))),ToReal(2)))); @@ -2398,12 +2395,6 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co CCTK_REAL_VEC Xt3rhsL CCTK_ATTRIBUTE_UNUSED = dotXt3; - CCTK_REAL_VEC Xtdot1L CCTK_ATTRIBUTE_UNUSED = dotXt1; - - CCTK_REAL_VEC Xtdot2L CCTK_ATTRIBUTE_UNUSED = dotXt2; - - CCTK_REAL_VEC Xtdot3L CCTK_ATTRIBUTE_UNUSED = dotXt3; - term1L = kneg(kmul(em4phi,knmsub(JacPDstandardNth1alpha,Xtn1,knmsub(JacPDstandardNth2alpha,Xtn2,knmsub(JacPDstandardNth3alpha,Xtn3,kmadd(gtu11,kmadd(cdphi1,kmul(JacPDstandardNth1alpha,ToReal(2)),JacPDstandardNth11alpha),kmadd(gtu12,kadd(JacPDstandardNth12alpha,kadd(JacPDstandardNth21alpha,kmadd(cdphi2,kmul(JacPDstandardNth1alpha,ToReal(2)),kmul(cdphi1,kmul(JacPDstandardNth2alpha,ToReal(2)))))),kmadd(gtu22,kmadd(cdphi2,kmul(JacPDstandardNth2alpha,ToReal(2)),JacPDstandardNth22alpha),kmadd(gtu13,kadd(JacPDstandardNth13alpha,kadd(JacPDstandardNth31alpha,kmadd(cdphi3,kmul(JacPDstandardNth1alpha,ToReal(2)),kmul(cdphi1,kmul(JacPDstandardNth3alpha,ToReal(2)))))),kmadd(gtu23,kadd(JacPDstandardNth23alpha,kadd(JacPDstandardNth32alpha,kmadd(cdphi3,kmul(JacPDstandardNth2alpha,ToReal(2)),kmul(cdphi2,kmul(JacPDstandardNth3alpha,ToReal(2)))))),kmul(gtu33,kmadd(cdphi3,kmul(JacPDstandardNth3alpha,ToReal(2)),JacPDstandardNth33alpha)))))))))))); @@ -2457,24 +2448,6 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co CCTK_REAL_VEC At33rhsL CCTK_ATTRIBUTE_UNUSED = kmadd(At33L,kmul(kadd(JacPDstandardNth1beta1,kadd(JacPDstandardNth2beta2,JacPDstandardNth3beta3)),ToReal(-0.666666666666666666666666666667)),kmadd(em4phi,kmadd(g33,kmul(trAts,ToReal(-0.333333333333333333333333333333)),Ats33),kmadd(kmadd(At13L,JacPDstandardNth3beta1,kmadd(At23L,JacPDstandardNth3beta2,kmul(At33L,JacPDstandardNth3beta3))),ToReal(2),kmul(alphaL,kmadd(At33L,trKL,kmadd(kmadd(At13L,Atm13,kmadd(At23L,Atm23,kmul(At33L,Atm33))),ToReal(-2),kmul(em4phi,kmul(ToReal(-8),kmul(kmadd(g33,kmul(trS,ToReal(-0.333333333333333333333333333333)),eTzzL),ToReal(3.14159265358979323846264338328)))))))))); - CCTK_REAL_VEC Kdot11L CCTK_ATTRIBUTE_UNUSED = - kmadd(G111,JacPDstandardNth1alpha,kmadd(G211,JacPDstandardNth2alpha,kmadd(G311,JacPDstandardNth3alpha,kmsub(alphaL,kmadd(trKL,K11,kmadd(kmadd(K11,Km11,kmadd(K12,Km12,kmul(K13,Km13))),ToReal(-2),R11)),JacPDstandardNth11alpha)))); - - CCTK_REAL_VEC Kdot12L CCTK_ATTRIBUTE_UNUSED = - kmadd(G112,JacPDstandardNth1alpha,kmadd(G212,JacPDstandardNth2alpha,kmadd(G312,JacPDstandardNth3alpha,kmsub(alphaL,kmadd(trKL,K12,kmadd(kmadd(K11,Km21,kmadd(K12,Km22,kmul(K13,Km23))),ToReal(-2),R12)),JacPDstandardNth12alpha)))); - - CCTK_REAL_VEC Kdot13L CCTK_ATTRIBUTE_UNUSED = - kmadd(G113,JacPDstandardNth1alpha,kmadd(G213,JacPDstandardNth2alpha,kmadd(G313,JacPDstandardNth3alpha,kmsub(alphaL,kmadd(trKL,K13,kmadd(kmadd(K11,Km31,kmadd(K12,Km32,kmul(K13,Km33))),ToReal(-2),R13)),JacPDstandardNth13alpha)))); - - CCTK_REAL_VEC Kdot22L CCTK_ATTRIBUTE_UNUSED = - kmadd(G122,JacPDstandardNth1alpha,kmadd(G222,JacPDstandardNth2alpha,kmadd(G322,JacPDstandardNth3alpha,kmsub(alphaL,kmadd(trKL,K22,kmadd(kmadd(K12,Km21,kmadd(K22,Km22,kmul(K23,Km23))),ToReal(-2),R22)),JacPDstandardNth22alpha)))); - - CCTK_REAL_VEC Kdot23L CCTK_ATTRIBUTE_UNUSED = - kmadd(G123,JacPDstandardNth1alpha,kmadd(G223,JacPDstandardNth2alpha,kmadd(G323,JacPDstandardNth3alpha,kmsub(alphaL,kmadd(trKL,K23,kmadd(kmadd(K12,Km31,kmadd(K22,Km32,kmul(K23,Km33))),ToReal(-2),R23)),JacPDstandardNth23alpha)))); - - CCTK_REAL_VEC Kdot33L CCTK_ATTRIBUTE_UNUSED = - kmadd(G133,JacPDstandardNth1alpha,kmadd(G233,JacPDstandardNth2alpha,kmadd(G333,JacPDstandardNth3alpha,kmsub(alphaL,kmadd(trKL,K33,kmadd(kmadd(K13,Km31,kmadd(K23,Km32,kmul(K33,Km33))),ToReal(-2),R33)),JacPDstandardNth33alpha)))); - CCTK_REAL_VEC eta CCTK_ATTRIBUTE_UNUSED = kdiv(ToReal(SpatialBetaDriverRadius),kfmax(rL,ToReal(SpatialBetaDriverRadius))); @@ -2495,6 +2468,7 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co CCTK_REAL_VEC beta3rhsL CCTK_ATTRIBUTE_UNUSED; if (harmonicShift) +#if 0 { beta1rhsL = kmul(alphaL,kmul(em4phi,kmul(ToReal(-0.5),kmadd(gtu11,kmadd(alphaL,kmadd(kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu12,kadd(JacPDstandardNth1gt12,JacPDstandardNth2gt11),kmadd(gtu22,JacPDstandardNth2gt12,kmadd(gtu13,kadd(JacPDstandardNth1gt13,JacPDstandardNth3gt11),kmadd(gtu23,kadd(JacPDstandardNth2gt13,JacPDstandardNth3gt12),kmadd(gtu33,JacPDstandardNth3gt13,kmul(JacPDstandardNth1phi,IfThen(conformalMethod @@ -2523,19 +2497,57 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co == 1,kdiv(ToReal(1),phiL),ToReal(-2))))))))),ToReal(-2),Ddetgt3),kmul(JacPDstandardNth3alpha,ToReal(2))))))))); } +#else + { + double gu[3][3] = {{ gu11, gu12, gu13 }, + { gu12, gu22, gu23 }, + { gu13, gu23, gu33 }}; + double G[3][3][3] = {{{ G111, G112, G113 }, + { G112, G122, G123 }, + { G113, G123, G133 }}, + {{ G211, G212, G213 }, + { G212, G222, G223 }, + { G213, G223, G233 }}, + {{ G311, G312, G313 }, + { G312, G322, G323 }, + { G313, G323, G333 }}}; + + double X1 = 0.0, X3 = 0.0; + + double dux_alpha, duz_alpha; + + for (int fn = 0; fn < 3; fn++) + for (int bz = 0; bz < 3; bz++) { + X1 += gu[fn][bz] * G[0][fn][bz]; + X3 += gu[fn][bz] * G[2][fn][bz]; + } + + dux_alpha = gu[0][0] * PDstandardNth1alpha + gu[0][2] * PDstandardNth3alpha; + duz_alpha = gu[2][0] * PDstandardNth1alpha + gu[2][2] * PDstandardNth3alpha; + + beta1rhsL = alphaL * alphaL * X1 - alphaL * dux_alpha - alphaL * BetaDriver * beta1L; + beta2rhsL = 0.0; + beta3rhsL = alphaL * alphaL * X3 - alphaL * duz_alpha - alphaL * BetaDriver * beta3L; + } +#endif else { - beta1rhsL = - kmul(theta,kmul(kadd(Xt1L,kmadd(ksub(B1L,Xt1L),ToReal(ShiftBCoeff),kmul(beta1L,kmul(eta,ToReal(BetaDriver*(-1 - + ShiftBCoeff)))))),ToReal(ShiftGammaCoeff))); - - beta2rhsL = - kmul(theta,kmul(kadd(Xt2L,kmadd(ksub(B2L,Xt2L),ToReal(ShiftBCoeff),kmul(beta2L,kmul(eta,ToReal(BetaDriver*(-1 - + ShiftBCoeff)))))),ToReal(ShiftGammaCoeff))); - - beta3rhsL = - kmul(theta,kmul(kadd(Xt3L,kmadd(ksub(B3L,Xt3L),ToReal(ShiftBCoeff),kmul(beta3L,kmul(eta,ToReal(BetaDriver*(-1 - + ShiftBCoeff)))))),ToReal(ShiftGammaCoeff))); + double xx = x[index]; + double zz = z[index]; + double t = cctkGH->cctk_time; + //double rr = sqrt(xx * xx + zz * zz); + //double phi = atan2(zz, xx); + + //double gauge_r = -20.0 * (1.0 + 0.0 * cos(2 * phi) * rr * rr) * rr / (1.0 + rr * rr * rr * rr); + + beta1rhsL = theta * ShiftGammaCoeff * + (ShiftBCoeff * B1L + (1.0 - ShiftBCoeff) * (Xt1L - eta * BetaDriver * beta1L)); + + beta2rhsL = theta * ShiftGammaCoeff * + (ShiftBCoeff * B2L + (1.0 - ShiftBCoeff) * (Xt2L - eta * BetaDriver * beta2L)); + + beta3rhsL = theta * ShiftGammaCoeff * + (ShiftBCoeff * B3L + (1.0 - ShiftBCoeff) * (Xt3L - eta * BetaDriver * beta3L)); } CCTK_REAL_VEC B1rhsL CCTK_ATTRIBUTE_UNUSED = @@ -2567,13 +2579,6 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co vec_store_nta_partial(gt22rhs[index],gt22rhsL); vec_store_nta_partial(gt23rhs[index],gt23rhsL); vec_store_nta_partial(gt33rhs[index],gt33rhsL); - vec_store_nta_partial(Kdot11[index],Kdot11L); - vec_store_nta_partial(Kdot12[index],Kdot12L); - vec_store_nta_partial(Kdot13[index],Kdot13L); - vec_store_nta_partial(Kdot22[index],Kdot22L); - vec_store_nta_partial(Kdot23[index],Kdot23L); - vec_store_nta_partial(Kdot33[index],Kdot33L); - vec_store_nta_partial(phidot[index],phidotL); vec_store_nta_partial(phirhs[index],phirhsL); vec_store_nta_partial(term1[index],term1L); vec_store_nta_partial(term2[index],term2L); @@ -2582,9 +2587,19 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co vec_store_nta_partial(Xt1rhs[index],Xt1rhsL); vec_store_nta_partial(Xt2rhs[index],Xt2rhsL); vec_store_nta_partial(Xt3rhs[index],Xt3rhsL); - vec_store_nta_partial(Xtdot1[index],Xtdot1L); - vec_store_nta_partial(Xtdot2[index],Xtdot2L); - vec_store_nta_partial(Xtdot3[index],Xtdot3L); + + //if (z[index] == 0.0 && x[index - 5] == 0.0) { + // fprintf(stderr, "%d %g", cctkGH->cctk_levfac[0], cctkGH->cctk_time); + // for (int i = 0; i < 5; i++) { + // double xx = x[index + i - 5]; + // int origin = fabs(xx) < 1e-8; + // fprintf(stderr, " %16.16g ", PDstandardNthfdOrder433(&alpha[index + i - 5])); + // } + // fprintf(stderr, " | "); + // for (int i = 0; i < 5; i++) + // fprintf(stderr, " %g ", term1[index - 5 + i]); + // fprintf(stderr, "\n\n"); + //} } CCTK_ENDLOOP3STR(ML_BSSN_RHS); } @@ -2605,60 +2620,6 @@ extern "C" void ML_BSSN_RHS(CCTK_ARGUMENTS) return; } - const char* const groups[] = { - "grid::coordinates", - "ML_BSSN::ML_curv", - "ML_BSSN::ML_curvrhs", - "ML_BSSN::ML_dtshift", - "ML_BSSN::ML_dtshiftrhs", - "ML_BSSN::ML_Gamma", - "ML_BSSN::ML_Gammarhs", - "ML_BSSN::ML_Kdot", - "ML_BSSN::ML_lapse", - "ML_BSSN::ML_log_confac", - "ML_BSSN::ML_log_confacrhs", - "ML_BSSN::ML_metric", - "ML_BSSN::ML_metricrhs", - "ML_BSSN::ML_phidot", - "ML_BSSN::ML_shift", - "ML_BSSN::ML_shiftrhs", - "ML_BSSN::ML_term1", - "ML_BSSN::ML_term2", - "ML_BSSN::ML_term3", - "ML_BSSN::ML_trace_curv", - "ML_BSSN::ML_trace_curvrhs", - "ML_BSSN::ML_Xtdot"}; - GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_RHS", 22, groups); - - switch (fdOrder) - { - case 2: - { - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_RHS", 1, 1, 1); - break; - } - - case 4: - { - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_RHS", 2, 2, 2); - break; - } - - case 6: - { - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_RHS", 3, 3, 3); - break; - } - - case 8: - { - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_RHS", 4, 4, 4); - break; - } - default: - CCTK_BUILTIN_UNREACHABLE(); - } - GenericFD_LoopOverInterior(cctkGH, ML_BSSN_RHS_Body); if (verbose > 1) |