From 28a757ee226406ad91156fa5620c3c047fcd5797 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Thu, 15 Dec 2011 21:45:04 -0500 Subject: Regenerate code --- ML_BSSN/param.ccl | 6 + ML_BSSN/src/ML_BSSN_RHS1.cc | 38 ++- .../src/ML_BSSN_convertToADMBaseDtLapseShift.cc | 304 ++++++++++++++++++++- ...ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc | 35 ++- .../ML_BSSN_convertToADMBaseFakeDtLapseShift.cc | 35 ++- 5 files changed, 380 insertions(+), 38 deletions(-) (limited to 'ML_BSSN') diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl index e73534d..acf1f14 100644 --- a/ML_BSSN/param.ccl +++ b/ML_BSSN/param.ccl @@ -149,6 +149,12 @@ CCTK_INT fdOrder "fdOrder" *:* :: "" } 4 +restricted: +CCTK_INT harmonicShift "Whether to use the harmonic shift" +{ + *:* :: "" +} 0 + private: KEYWORD my_initial_data "my_initial_data" { diff --git a/ML_BSSN/src/ML_BSSN_RHS1.cc b/ML_BSSN/src/ML_BSSN_RHS1.cc index 9e64eab..5b5f80a 100644 --- a/ML_BSSN/src/ML_BSSN_RHS1.cc +++ b/ML_BSSN/src/ML_BSSN_RHS1.cc @@ -1338,17 +1338,35 @@ static void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, CCTK_REAL_VEC theta = kfmin(ToReal(1),kexp(knmsub(rL,INV(ToReal(SpatialShiftGammaCoeffRadius)),ToReal(1)))); - CCTK_REAL_VEC beta1rhsL = - kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + CCTK_REAL_VEC beta1rhsL; + CCTK_REAL_VEC beta2rhsL; + CCTK_REAL_VEC beta3rhsL; - CCTK_REAL_VEC beta2rhsL = - kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); - - CCTK_REAL_VEC beta3rhsL = - kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + if (harmonicShift) + { + beta1rhsL = + kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu11,JacPDstandardNth1alpha,kmadd(gtu12,JacPDstandardNth2alpha,kmul(gtu13,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmadd(JacPDstandardNth1gt11,SQR(gtu11),kmul(JacPDstandardNth1gt22,kmul(SQR(gtu12),ToReal(2)))),kmadd(gtu13,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,kmsub(kmadd(gtu13,JacPDstandardNth1gt33,kmadd(gtu22,JacPDstandardNth2gt23,kmul(gtu23,JacPDstandardNth2gt33))),ToReal(2),kmul(gtu22,JacPDstandardNth3gt22))))),kmadd(gtu11,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu12,JacPDstandardNth2gt11,kmadd(gtu13,JacPDstandardNth3gt11,kmadd(gtu23,kmul(JacPDstandardNth1gt23,ToReal(-2)),knmsub(gtu22,JacPDstandardNth1gt22,kmadd(kmadd(gtu12,JacPDstandardNth1gt12,kmadd(gtu13,JacPDstandardNth1gt13,kmul(gtu22,JacPDstandardNth2gt12))),ToReal(2),kmadd(gtu23,kmul(JacPDstandardNth2gt13,ToReal(2)),kmadd(gtu23,kmul(JacPDstandardNth3gt12,ToReal(2)),kmul(gtu33,kmsub(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33))))))))))),kmul(gtu12,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu22,JacPDstandardNth2gt22,kmadd(gtu23,kmul(JacPDstandardNth3gt22,ToReal(2)),kmadd(gtu33,kmsub(JacPDstandardNth3gt23,ToReal(2),JacPDstandardNth2gt33),kmul(gtu13,kmul(JacPDstandardNth1gt23,ToReal(4))))))))))))))))); + + beta2rhsL = + kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu12,JacPDstandardNth1alpha,kmadd(gtu22,JacPDstandardNth2alpha,kmul(gtu23,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmadd(JacPDstandardNth2gt22,SQR(gtu22),kmul(JacPDstandardNth2gt11,kmul(SQR(gtu12),ToReal(2)))),kmadd(gtu23,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,kmsub(kmadd(gtu11,JacPDstandardNth1gt13,kmadd(gtu13,JacPDstandardNth1gt33,kmul(gtu23,JacPDstandardNth2gt33))),ToReal(2),kmul(gtu11,JacPDstandardNth3gt11))))),kmadd(gtu22,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu23,JacPDstandardNth3gt22,kmadd(kmadd(gtu23,JacPDstandardNth2gt23,kmul(gtu13,kadd(JacPDstandardNth1gt23,ksub(JacPDstandardNth3gt12,JacPDstandardNth2gt13)))),ToReal(2),kmadd(gtu11,kmsub(JacPDstandardNth1gt12,ToReal(2),JacPDstandardNth2gt11),kmul(gtu33,kmsub(JacPDstandardNth3gt23,ToReal(2),JacPDstandardNth2gt33))))))),kmul(gtu12,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu13,kmul(JacPDstandardNth3gt11,ToReal(2)),kmadd(gtu22,kmadd(JacPDstandardNth2gt12,ToReal(2),JacPDstandardNth1gt22),kmadd(gtu33,kmsub(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33),kmul(gtu23,kmul(JacPDstandardNth2gt13,ToReal(4)))))))))))))))))); + + beta3rhsL = + kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu13,JacPDstandardNth1alpha,kmadd(gtu23,JacPDstandardNth2alpha,kmul(gtu33,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmul(kmadd(JacPDstandardNth3gt11,SQR(gtu13),kmul(JacPDstandardNth3gt22,SQR(gtu23))),ToReal(2)),kmadd(gtu23,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu22,JacPDstandardNth2gt22,kmadd(gtu33,JacPDstandardNth2gt33,kmsub(kmadd(gtu11,JacPDstandardNth1gt12,kmadd(gtu12,JacPDstandardNth1gt22,kmul(gtu33,JacPDstandardNth3gt23))),ToReal(2),kmul(gtu11,JacPDstandardNth2gt11)))))),kmadd(gtu33,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,knmsub(gtu22,JacPDstandardNth3gt22,kmadd(kmadd(gtu22,JacPDstandardNth2gt23,kmul(gtu12,kadd(JacPDstandardNth1gt23,ksub(JacPDstandardNth2gt13,JacPDstandardNth3gt12)))),ToReal(2),kmul(gtu11,kmsub(JacPDstandardNth1gt13,ToReal(2),JacPDstandardNth3gt11))))))),kmul(gtu13,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu12,kmul(JacPDstandardNth2gt11,ToReal(2)),kmadd(gtu22,kmsub(JacPDstandardNth2gt12,ToReal(2),JacPDstandardNth1gt22),kmadd(gtu33,kmadd(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33),kmul(gtu23,kmul(JacPDstandardNth3gt12,ToReal(4)))))))))))))))))); + } + else + { + beta1rhsL = + kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + + beta2rhsL = + kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + + beta3rhsL = + kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + } CCTK_REAL_VEC B1rhsL = kmul(knmsub(B1L,kmul(eta,ToReal(BetaDriver)),dotXt1),ToReal(ShiftBCoeff)); diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.cc b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.cc index 6a76e6c..2abc29e 100644 --- a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.cc +++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.cc @@ -223,6 +223,13 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const CCTK_REAL_VEC beta1L = vec_load(beta1[index]); CCTK_REAL_VEC beta2L = vec_load(beta2[index]); CCTK_REAL_VEC beta3L = vec_load(beta3[index]); + CCTK_REAL_VEC gt11L = vec_load(gt11[index]); + CCTK_REAL_VEC gt12L = vec_load(gt12[index]); + CCTK_REAL_VEC gt13L = vec_load(gt13[index]); + CCTK_REAL_VEC gt22L = vec_load(gt22[index]); + CCTK_REAL_VEC gt23L = vec_load(gt23[index]); + CCTK_REAL_VEC gt33L = vec_load(gt33[index]); + CCTK_REAL_VEC phiL = vec_load(phi[index]); CCTK_REAL_VEC rL = vec_load(r[index]); CCTK_REAL_VEC trKL = vec_load(trK[index]); CCTK_REAL_VEC Xt1L = vec_load(Xt1[index]); @@ -248,6 +255,9 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const /* Include user supplied include files */ /* Precompute derivatives */ + CCTK_REAL_VEC PDstandardNth1alpha; + CCTK_REAL_VEC PDstandardNth2alpha; + CCTK_REAL_VEC PDstandardNth3alpha; CCTK_REAL_VEC PDupwindNthAnti1alpha; CCTK_REAL_VEC PDupwindNthSymm1alpha; CCTK_REAL_VEC PDupwindNthAnti2alpha; @@ -272,10 +282,34 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const CCTK_REAL_VEC PDupwindNthSymm2beta3; CCTK_REAL_VEC PDupwindNthAnti3beta3; CCTK_REAL_VEC PDupwindNthSymm3beta3; + CCTK_REAL_VEC PDstandardNth1gt11; + CCTK_REAL_VEC PDstandardNth2gt11; + CCTK_REAL_VEC PDstandardNth3gt11; + CCTK_REAL_VEC PDstandardNth1gt12; + CCTK_REAL_VEC PDstandardNth2gt12; + CCTK_REAL_VEC PDstandardNth3gt12; + CCTK_REAL_VEC PDstandardNth1gt13; + CCTK_REAL_VEC PDstandardNth2gt13; + CCTK_REAL_VEC PDstandardNth3gt13; + CCTK_REAL_VEC PDstandardNth1gt22; + CCTK_REAL_VEC PDstandardNth2gt22; + CCTK_REAL_VEC PDstandardNth3gt22; + CCTK_REAL_VEC PDstandardNth1gt23; + CCTK_REAL_VEC PDstandardNth2gt23; + CCTK_REAL_VEC PDstandardNth3gt23; + CCTK_REAL_VEC PDstandardNth1gt33; + CCTK_REAL_VEC PDstandardNth2gt33; + CCTK_REAL_VEC PDstandardNth3gt33; + CCTK_REAL_VEC PDstandardNth1phi; + CCTK_REAL_VEC PDstandardNth2phi; + CCTK_REAL_VEC PDstandardNth3phi; switch(fdOrder) { case 2: + PDstandardNth1alpha = PDstandardNthfdOrder21(&alpha[index]); + PDstandardNth2alpha = PDstandardNthfdOrder22(&alpha[index]); + PDstandardNth3alpha = PDstandardNthfdOrder23(&alpha[index]); PDupwindNthAnti1alpha = PDupwindNthAntifdOrder21(&alpha[index]); PDupwindNthSymm1alpha = PDupwindNthSymmfdOrder21(&alpha[index]); PDupwindNthAnti2alpha = PDupwindNthAntifdOrder22(&alpha[index]); @@ -300,9 +334,33 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const PDupwindNthSymm2beta3 = PDupwindNthSymmfdOrder22(&beta3[index]); PDupwindNthAnti3beta3 = PDupwindNthAntifdOrder23(&beta3[index]); PDupwindNthSymm3beta3 = PDupwindNthSymmfdOrder23(&beta3[index]); + PDstandardNth1gt11 = PDstandardNthfdOrder21(>11[index]); + PDstandardNth2gt11 = PDstandardNthfdOrder22(>11[index]); + PDstandardNth3gt11 = PDstandardNthfdOrder23(>11[index]); + PDstandardNth1gt12 = PDstandardNthfdOrder21(>12[index]); + PDstandardNth2gt12 = PDstandardNthfdOrder22(>12[index]); + PDstandardNth3gt12 = PDstandardNthfdOrder23(>12[index]); + PDstandardNth1gt13 = PDstandardNthfdOrder21(>13[index]); + PDstandardNth2gt13 = PDstandardNthfdOrder22(>13[index]); + PDstandardNth3gt13 = PDstandardNthfdOrder23(>13[index]); + PDstandardNth1gt22 = PDstandardNthfdOrder21(>22[index]); + PDstandardNth2gt22 = PDstandardNthfdOrder22(>22[index]); + PDstandardNth3gt22 = PDstandardNthfdOrder23(>22[index]); + PDstandardNth1gt23 = PDstandardNthfdOrder21(>23[index]); + PDstandardNth2gt23 = PDstandardNthfdOrder22(>23[index]); + PDstandardNth3gt23 = PDstandardNthfdOrder23(>23[index]); + PDstandardNth1gt33 = PDstandardNthfdOrder21(>33[index]); + PDstandardNth2gt33 = PDstandardNthfdOrder22(>33[index]); + PDstandardNth3gt33 = PDstandardNthfdOrder23(>33[index]); + PDstandardNth1phi = PDstandardNthfdOrder21(&phi[index]); + PDstandardNth2phi = PDstandardNthfdOrder22(&phi[index]); + PDstandardNth3phi = PDstandardNthfdOrder23(&phi[index]); break; case 4: + PDstandardNth1alpha = PDstandardNthfdOrder41(&alpha[index]); + PDstandardNth2alpha = PDstandardNthfdOrder42(&alpha[index]); + PDstandardNth3alpha = PDstandardNthfdOrder43(&alpha[index]); PDupwindNthAnti1alpha = PDupwindNthAntifdOrder41(&alpha[index]); PDupwindNthSymm1alpha = PDupwindNthSymmfdOrder41(&alpha[index]); PDupwindNthAnti2alpha = PDupwindNthAntifdOrder42(&alpha[index]); @@ -327,9 +385,33 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const PDupwindNthSymm2beta3 = PDupwindNthSymmfdOrder42(&beta3[index]); PDupwindNthAnti3beta3 = PDupwindNthAntifdOrder43(&beta3[index]); PDupwindNthSymm3beta3 = PDupwindNthSymmfdOrder43(&beta3[index]); + PDstandardNth1gt11 = PDstandardNthfdOrder41(>11[index]); + PDstandardNth2gt11 = PDstandardNthfdOrder42(>11[index]); + PDstandardNth3gt11 = PDstandardNthfdOrder43(>11[index]); + PDstandardNth1gt12 = PDstandardNthfdOrder41(>12[index]); + PDstandardNth2gt12 = PDstandardNthfdOrder42(>12[index]); + PDstandardNth3gt12 = PDstandardNthfdOrder43(>12[index]); + PDstandardNth1gt13 = PDstandardNthfdOrder41(>13[index]); + PDstandardNth2gt13 = PDstandardNthfdOrder42(>13[index]); + PDstandardNth3gt13 = PDstandardNthfdOrder43(>13[index]); + PDstandardNth1gt22 = PDstandardNthfdOrder41(>22[index]); + PDstandardNth2gt22 = PDstandardNthfdOrder42(>22[index]); + PDstandardNth3gt22 = PDstandardNthfdOrder43(>22[index]); + PDstandardNth1gt23 = PDstandardNthfdOrder41(>23[index]); + PDstandardNth2gt23 = PDstandardNthfdOrder42(>23[index]); + PDstandardNth3gt23 = PDstandardNthfdOrder43(>23[index]); + PDstandardNth1gt33 = PDstandardNthfdOrder41(>33[index]); + PDstandardNth2gt33 = PDstandardNthfdOrder42(>33[index]); + PDstandardNth3gt33 = PDstandardNthfdOrder43(>33[index]); + PDstandardNth1phi = PDstandardNthfdOrder41(&phi[index]); + PDstandardNth2phi = PDstandardNthfdOrder42(&phi[index]); + PDstandardNth3phi = PDstandardNthfdOrder43(&phi[index]); break; case 6: + PDstandardNth1alpha = PDstandardNthfdOrder61(&alpha[index]); + PDstandardNth2alpha = PDstandardNthfdOrder62(&alpha[index]); + PDstandardNth3alpha = PDstandardNthfdOrder63(&alpha[index]); PDupwindNthAnti1alpha = PDupwindNthAntifdOrder61(&alpha[index]); PDupwindNthSymm1alpha = PDupwindNthSymmfdOrder61(&alpha[index]); PDupwindNthAnti2alpha = PDupwindNthAntifdOrder62(&alpha[index]); @@ -354,9 +436,33 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const PDupwindNthSymm2beta3 = PDupwindNthSymmfdOrder62(&beta3[index]); PDupwindNthAnti3beta3 = PDupwindNthAntifdOrder63(&beta3[index]); PDupwindNthSymm3beta3 = PDupwindNthSymmfdOrder63(&beta3[index]); + PDstandardNth1gt11 = PDstandardNthfdOrder61(>11[index]); + PDstandardNth2gt11 = PDstandardNthfdOrder62(>11[index]); + PDstandardNth3gt11 = PDstandardNthfdOrder63(>11[index]); + PDstandardNth1gt12 = PDstandardNthfdOrder61(>12[index]); + PDstandardNth2gt12 = PDstandardNthfdOrder62(>12[index]); + PDstandardNth3gt12 = PDstandardNthfdOrder63(>12[index]); + PDstandardNth1gt13 = PDstandardNthfdOrder61(>13[index]); + PDstandardNth2gt13 = PDstandardNthfdOrder62(>13[index]); + PDstandardNth3gt13 = PDstandardNthfdOrder63(>13[index]); + PDstandardNth1gt22 = PDstandardNthfdOrder61(>22[index]); + PDstandardNth2gt22 = PDstandardNthfdOrder62(>22[index]); + PDstandardNth3gt22 = PDstandardNthfdOrder63(>22[index]); + PDstandardNth1gt23 = PDstandardNthfdOrder61(>23[index]); + PDstandardNth2gt23 = PDstandardNthfdOrder62(>23[index]); + PDstandardNth3gt23 = PDstandardNthfdOrder63(>23[index]); + PDstandardNth1gt33 = PDstandardNthfdOrder61(>33[index]); + PDstandardNth2gt33 = PDstandardNthfdOrder62(>33[index]); + PDstandardNth3gt33 = PDstandardNthfdOrder63(>33[index]); + PDstandardNth1phi = PDstandardNthfdOrder61(&phi[index]); + PDstandardNth2phi = PDstandardNthfdOrder62(&phi[index]); + PDstandardNth3phi = PDstandardNthfdOrder63(&phi[index]); break; case 8: + PDstandardNth1alpha = PDstandardNthfdOrder81(&alpha[index]); + PDstandardNth2alpha = PDstandardNthfdOrder82(&alpha[index]); + PDstandardNth3alpha = PDstandardNthfdOrder83(&alpha[index]); PDupwindNthAnti1alpha = PDupwindNthAntifdOrder81(&alpha[index]); PDupwindNthSymm1alpha = PDupwindNthSymmfdOrder81(&alpha[index]); PDupwindNthAnti2alpha = PDupwindNthAntifdOrder82(&alpha[index]); @@ -381,6 +487,27 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const PDupwindNthSymm2beta3 = PDupwindNthSymmfdOrder82(&beta3[index]); PDupwindNthAnti3beta3 = PDupwindNthAntifdOrder83(&beta3[index]); PDupwindNthSymm3beta3 = PDupwindNthSymmfdOrder83(&beta3[index]); + PDstandardNth1gt11 = PDstandardNthfdOrder81(>11[index]); + PDstandardNth2gt11 = PDstandardNthfdOrder82(>11[index]); + PDstandardNth3gt11 = PDstandardNthfdOrder83(>11[index]); + PDstandardNth1gt12 = PDstandardNthfdOrder81(>12[index]); + PDstandardNth2gt12 = PDstandardNthfdOrder82(>12[index]); + PDstandardNth3gt12 = PDstandardNthfdOrder83(>12[index]); + PDstandardNth1gt13 = PDstandardNthfdOrder81(>13[index]); + PDstandardNth2gt13 = PDstandardNthfdOrder82(>13[index]); + PDstandardNth3gt13 = PDstandardNthfdOrder83(>13[index]); + PDstandardNth1gt22 = PDstandardNthfdOrder81(>22[index]); + PDstandardNth2gt22 = PDstandardNthfdOrder82(>22[index]); + PDstandardNth3gt22 = PDstandardNthfdOrder83(>22[index]); + PDstandardNth1gt23 = PDstandardNthfdOrder81(>23[index]); + PDstandardNth2gt23 = PDstandardNthfdOrder82(>23[index]); + PDstandardNth3gt23 = PDstandardNthfdOrder83(>23[index]); + PDstandardNth1gt33 = PDstandardNthfdOrder81(>33[index]); + PDstandardNth2gt33 = PDstandardNthfdOrder82(>33[index]); + PDstandardNth3gt33 = PDstandardNthfdOrder83(>33[index]); + PDstandardNth1phi = PDstandardNthfdOrder81(&phi[index]); + PDstandardNth2phi = PDstandardNthfdOrder82(&phi[index]); + PDstandardNth3phi = PDstandardNthfdOrder83(&phi[index]); break; } @@ -391,6 +518,30 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const ptrdiff_t dir3 = Sign(beta3L); + CCTK_REAL_VEC JacPDstandardNth1alpha; + CCTK_REAL_VEC JacPDstandardNth1gt11; + CCTK_REAL_VEC JacPDstandardNth1gt12; + CCTK_REAL_VEC JacPDstandardNth1gt13; + CCTK_REAL_VEC JacPDstandardNth1gt22; + CCTK_REAL_VEC JacPDstandardNth1gt23; + CCTK_REAL_VEC JacPDstandardNth1gt33; + CCTK_REAL_VEC JacPDstandardNth1phi; + CCTK_REAL_VEC JacPDstandardNth2alpha; + CCTK_REAL_VEC JacPDstandardNth2gt11; + CCTK_REAL_VEC JacPDstandardNth2gt12; + CCTK_REAL_VEC JacPDstandardNth2gt13; + CCTK_REAL_VEC JacPDstandardNth2gt22; + CCTK_REAL_VEC JacPDstandardNth2gt23; + CCTK_REAL_VEC JacPDstandardNth2gt33; + CCTK_REAL_VEC JacPDstandardNth2phi; + CCTK_REAL_VEC JacPDstandardNth3alpha; + CCTK_REAL_VEC JacPDstandardNth3gt11; + CCTK_REAL_VEC JacPDstandardNth3gt12; + CCTK_REAL_VEC JacPDstandardNth3gt13; + CCTK_REAL_VEC JacPDstandardNth3gt22; + CCTK_REAL_VEC JacPDstandardNth3gt23; + CCTK_REAL_VEC JacPDstandardNth3gt33; + CCTK_REAL_VEC JacPDstandardNth3phi; CCTK_REAL_VEC JacPDupwindNthAnti1alpha; CCTK_REAL_VEC JacPDupwindNthAnti1beta1; CCTK_REAL_VEC JacPDupwindNthAnti1beta2; @@ -418,6 +569,78 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const if (use_jacobian) { + JacPDstandardNth1alpha = + kmadd(J11L,PDstandardNth1alpha,kmadd(J21L,PDstandardNth2alpha,kmul(J31L,PDstandardNth3alpha))); + + JacPDstandardNth1gt11 = + kmadd(J11L,PDstandardNth1gt11,kmadd(J21L,PDstandardNth2gt11,kmul(J31L,PDstandardNth3gt11))); + + JacPDstandardNth1gt12 = + kmadd(J11L,PDstandardNth1gt12,kmadd(J21L,PDstandardNth2gt12,kmul(J31L,PDstandardNth3gt12))); + + JacPDstandardNth1gt13 = + kmadd(J11L,PDstandardNth1gt13,kmadd(J21L,PDstandardNth2gt13,kmul(J31L,PDstandardNth3gt13))); + + JacPDstandardNth1gt22 = + kmadd(J11L,PDstandardNth1gt22,kmadd(J21L,PDstandardNth2gt22,kmul(J31L,PDstandardNth3gt22))); + + JacPDstandardNth1gt23 = + kmadd(J11L,PDstandardNth1gt23,kmadd(J21L,PDstandardNth2gt23,kmul(J31L,PDstandardNth3gt23))); + + JacPDstandardNth1gt33 = + kmadd(J11L,PDstandardNth1gt33,kmadd(J21L,PDstandardNth2gt33,kmul(J31L,PDstandardNth3gt33))); + + JacPDstandardNth1phi = + kmadd(J11L,PDstandardNth1phi,kmadd(J21L,PDstandardNth2phi,kmul(J31L,PDstandardNth3phi))); + + JacPDstandardNth2alpha = + kmadd(J12L,PDstandardNth1alpha,kmadd(J22L,PDstandardNth2alpha,kmul(J32L,PDstandardNth3alpha))); + + JacPDstandardNth2gt11 = + kmadd(J12L,PDstandardNth1gt11,kmadd(J22L,PDstandardNth2gt11,kmul(J32L,PDstandardNth3gt11))); + + JacPDstandardNth2gt12 = + kmadd(J12L,PDstandardNth1gt12,kmadd(J22L,PDstandardNth2gt12,kmul(J32L,PDstandardNth3gt12))); + + JacPDstandardNth2gt13 = + kmadd(J12L,PDstandardNth1gt13,kmadd(J22L,PDstandardNth2gt13,kmul(J32L,PDstandardNth3gt13))); + + JacPDstandardNth2gt22 = + kmadd(J12L,PDstandardNth1gt22,kmadd(J22L,PDstandardNth2gt22,kmul(J32L,PDstandardNth3gt22))); + + JacPDstandardNth2gt23 = + kmadd(J12L,PDstandardNth1gt23,kmadd(J22L,PDstandardNth2gt23,kmul(J32L,PDstandardNth3gt23))); + + JacPDstandardNth2gt33 = + kmadd(J12L,PDstandardNth1gt33,kmadd(J22L,PDstandardNth2gt33,kmul(J32L,PDstandardNth3gt33))); + + JacPDstandardNth2phi = + kmadd(J12L,PDstandardNth1phi,kmadd(J22L,PDstandardNth2phi,kmul(J32L,PDstandardNth3phi))); + + JacPDstandardNth3alpha = + kmadd(J13L,PDstandardNth1alpha,kmadd(J23L,PDstandardNth2alpha,kmul(J33L,PDstandardNth3alpha))); + + JacPDstandardNth3gt11 = + kmadd(J13L,PDstandardNth1gt11,kmadd(J23L,PDstandardNth2gt11,kmul(J33L,PDstandardNth3gt11))); + + JacPDstandardNth3gt12 = + kmadd(J13L,PDstandardNth1gt12,kmadd(J23L,PDstandardNth2gt12,kmul(J33L,PDstandardNth3gt12))); + + JacPDstandardNth3gt13 = + kmadd(J13L,PDstandardNth1gt13,kmadd(J23L,PDstandardNth2gt13,kmul(J33L,PDstandardNth3gt13))); + + JacPDstandardNth3gt22 = + kmadd(J13L,PDstandardNth1gt22,kmadd(J23L,PDstandardNth2gt22,kmul(J33L,PDstandardNth3gt22))); + + JacPDstandardNth3gt23 = + kmadd(J13L,PDstandardNth1gt23,kmadd(J23L,PDstandardNth2gt23,kmul(J33L,PDstandardNth3gt23))); + + JacPDstandardNth3gt33 = + kmadd(J13L,PDstandardNth1gt33,kmadd(J23L,PDstandardNth2gt33,kmul(J33L,PDstandardNth3gt33))); + + JacPDstandardNth3phi = + kmadd(J13L,PDstandardNth1phi,kmadd(J23L,PDstandardNth2phi,kmul(J33L,PDstandardNth3phi))); + JacPDupwindNthAnti1alpha = kmadd(J11L,PDupwindNthAnti1alpha,kmadd(J21L,PDupwindNthAnti2alpha,kmul(J31L,PDupwindNthAnti3alpha))); @@ -492,6 +715,54 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const } else { + JacPDstandardNth1alpha = PDstandardNth1alpha; + + JacPDstandardNth1gt11 = PDstandardNth1gt11; + + JacPDstandardNth1gt12 = PDstandardNth1gt12; + + JacPDstandardNth1gt13 = PDstandardNth1gt13; + + JacPDstandardNth1gt22 = PDstandardNth1gt22; + + JacPDstandardNth1gt23 = PDstandardNth1gt23; + + JacPDstandardNth1gt33 = PDstandardNth1gt33; + + JacPDstandardNth1phi = PDstandardNth1phi; + + JacPDstandardNth2alpha = PDstandardNth2alpha; + + JacPDstandardNth2gt11 = PDstandardNth2gt11; + + JacPDstandardNth2gt12 = PDstandardNth2gt12; + + JacPDstandardNth2gt13 = PDstandardNth2gt13; + + JacPDstandardNth2gt22 = PDstandardNth2gt22; + + JacPDstandardNth2gt23 = PDstandardNth2gt23; + + JacPDstandardNth2gt33 = PDstandardNth2gt33; + + JacPDstandardNth2phi = PDstandardNth2phi; + + JacPDstandardNth3alpha = PDstandardNth3alpha; + + JacPDstandardNth3gt11 = PDstandardNth3gt11; + + JacPDstandardNth3gt12 = PDstandardNth3gt12; + + JacPDstandardNth3gt13 = PDstandardNth3gt13; + + JacPDstandardNth3gt22 = PDstandardNth3gt22; + + JacPDstandardNth3gt23 = PDstandardNth3gt23; + + JacPDstandardNth3gt33 = PDstandardNth3gt33; + + JacPDstandardNth3phi = PDstandardNth3phi; + JacPDupwindNthAnti1alpha = PDupwindNthAnti1alpha; JacPDupwindNthAnti1beta1 = PDupwindNthAnti1beta1; @@ -541,6 +812,23 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const JacPDupwindNthSymm3beta3 = PDupwindNthSymm3beta3; } + CCTK_REAL_VEC detgt = ToReal(1); + + CCTK_REAL_VEC gtu11 = kmul(INV(detgt),kmsub(gt22L,gt33L,SQR(gt23L))); + + CCTK_REAL_VEC gtu12 = + kmul(INV(detgt),kmsub(gt13L,gt23L,kmul(gt12L,gt33L))); + + CCTK_REAL_VEC gtu13 = + kmul(INV(detgt),kmsub(gt12L,gt23L,kmul(gt13L,gt22L))); + + CCTK_REAL_VEC gtu22 = kmul(INV(detgt),kmsub(gt11L,gt33L,SQR(gt13L))); + + CCTK_REAL_VEC gtu23 = + kmul(INV(detgt),kmsub(gt12L,gt13L,kmul(gt11L,gt23L))); + + CCTK_REAL_VEC gtu33 = kmul(INV(detgt),kmsub(gt11L,gt22L,SQR(gt12L))); + CCTK_REAL_VEC eta = kfmin(ToReal(1),kmul(INV(rL),ToReal(SpatialBetaDriverRadius))); @@ -551,19 +839,19 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const kmsub(kmadd(beta1L,JacPDupwindNthAnti1alpha,kmadd(beta2L,JacPDupwindNthAnti2alpha,kmadd(beta3L,JacPDupwindNthAnti3alpha,kmadd(JacPDupwindNthSymm1alpha,kfabs(beta1L),kmadd(JacPDupwindNthSymm2alpha,kfabs(beta2L),kmul(JacPDupwindNthSymm3alpha,kfabs(beta3L))))))),ToReal(LapseAdvectionCoeff),kmul(kpow(alphaL,harmonicN),kmul(ToReal(harmonicF),kmadd(ksub(AL,trKL),ToReal(LapseACoeff),trKL)))); CCTK_REAL_VEC dtbetaxL = - kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta1,kmadd(beta2L,JacPDupwindNthAnti2beta1,kmadd(beta3L,JacPDupwindNthAnti3beta1,kmadd(JacPDupwindNthSymm1beta1,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta1,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta1,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta1,kmadd(beta2L,JacPDupwindNthAnti2beta1,kmadd(beta3L,JacPDupwindNthAnti3beta1,kmadd(JacPDupwindNthSymm1beta1,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta1,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta1,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu11,JacPDstandardNth1alpha,kmadd(gtu12,JacPDstandardNth2alpha,kmul(gtu13,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmadd(JacPDstandardNth1gt11,SQR(gtu11),kmul(JacPDstandardNth1gt22,kmul(SQR(gtu12),ToReal(2)))),kmadd(gtu13,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,kmsub(kmadd(gtu13,JacPDstandardNth1gt33,kmadd(gtu22,JacPDstandardNth2gt23,kmul(gtu23,JacPDstandardNth2gt33))),ToReal(2),kmul(gtu22,JacPDstandardNth3gt22))))),kmadd(gtu11,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu12,JacPDstandardNth2gt11,kmadd(gtu13,JacPDstandardNth3gt11,kmadd(gtu23,kmul(JacPDstandardNth1gt23,ToReal(-2)),knmsub(gtu22,JacPDstandardNth1gt22,kmadd(kmadd(gtu12,JacPDstandardNth1gt12,kmadd(gtu13,JacPDstandardNth1gt13,kmul(gtu22,JacPDstandardNth2gt12))),ToReal(2),kmadd(gtu23,kmul(JacPDstandardNth2gt13,ToReal(2)),kmadd(gtu23,kmul(JacPDstandardNth3gt12,ToReal(2)),kmul(gtu33,kmsub(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33))))))))))),kmul(gtu12,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu22,JacPDstandardNth2gt22,kmadd(gtu23,kmul(JacPDstandardNth3gt22,ToReal(2)),kmadd(gtu33,kmsub(JacPDstandardNth3gt23,ToReal(2),JacPDstandardNth2gt33),kmul(gtu13,kmul(JacPDstandardNth1gt23,ToReal(4))))))))))))))))),kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)))); + ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))))); CCTK_REAL_VEC dtbetayL = - kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta2,kmadd(beta2L,JacPDupwindNthAnti2beta2,kmadd(beta3L,JacPDupwindNthAnti3beta2,kmadd(JacPDupwindNthSymm1beta2,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta2,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta2,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta2,kmadd(beta2L,JacPDupwindNthAnti2beta2,kmadd(beta3L,JacPDupwindNthAnti3beta2,kmadd(JacPDupwindNthSymm1beta2,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta2,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta2,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu12,JacPDstandardNth1alpha,kmadd(gtu22,JacPDstandardNth2alpha,kmul(gtu23,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmadd(JacPDstandardNth2gt22,SQR(gtu22),kmul(JacPDstandardNth2gt11,kmul(SQR(gtu12),ToReal(2)))),kmadd(gtu23,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,kmsub(kmadd(gtu11,JacPDstandardNth1gt13,kmadd(gtu13,JacPDstandardNth1gt33,kmul(gtu23,JacPDstandardNth2gt33))),ToReal(2),kmul(gtu11,JacPDstandardNth3gt11))))),kmadd(gtu22,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu23,JacPDstandardNth3gt22,kmadd(kmadd(gtu23,JacPDstandardNth2gt23,kmul(gtu13,kadd(JacPDstandardNth1gt23,ksub(JacPDstandardNth3gt12,JacPDstandardNth2gt13)))),ToReal(2),kmadd(gtu11,kmsub(JacPDstandardNth1gt12,ToReal(2),JacPDstandardNth2gt11),kmul(gtu33,kmsub(JacPDstandardNth3gt23,ToReal(2),JacPDstandardNth2gt33))))))),kmul(gtu12,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu13,kmul(JacPDstandardNth3gt11,ToReal(2)),kmadd(gtu22,kmadd(JacPDstandardNth2gt12,ToReal(2),JacPDstandardNth1gt22),kmadd(gtu33,kmsub(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33),kmul(gtu23,kmul(JacPDstandardNth2gt13,ToReal(4)))))))))))))))))),kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)))); + ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))))); CCTK_REAL_VEC dtbetazL = - kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta3,kmadd(beta2L,JacPDupwindNthAnti2beta3,kmadd(beta3L,JacPDupwindNthAnti3beta3,kmadd(JacPDupwindNthSymm1beta3,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta3,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta3,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta3,kmadd(beta2L,JacPDupwindNthAnti2beta3,kmadd(beta3L,JacPDupwindNthAnti3beta3,kmadd(JacPDupwindNthSymm1beta3,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta3,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta3,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu13,JacPDstandardNth1alpha,kmadd(gtu23,JacPDstandardNth2alpha,kmul(gtu33,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmul(kmadd(JacPDstandardNth3gt11,SQR(gtu13),kmul(JacPDstandardNth3gt22,SQR(gtu23))),ToReal(2)),kmadd(gtu23,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu22,JacPDstandardNth2gt22,kmadd(gtu33,JacPDstandardNth2gt33,kmsub(kmadd(gtu11,JacPDstandardNth1gt12,kmadd(gtu12,JacPDstandardNth1gt22,kmul(gtu33,JacPDstandardNth3gt23))),ToReal(2),kmul(gtu11,JacPDstandardNth2gt11)))))),kmadd(gtu33,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,knmsub(gtu22,JacPDstandardNth3gt22,kmadd(kmadd(gtu22,JacPDstandardNth2gt23,kmul(gtu12,kadd(JacPDstandardNth1gt23,ksub(JacPDstandardNth2gt13,JacPDstandardNth3gt12)))),ToReal(2),kmul(gtu11,kmsub(JacPDstandardNth1gt13,ToReal(2),JacPDstandardNth3gt11))))))),kmul(gtu13,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu12,kmul(JacPDstandardNth2gt11,ToReal(2)),kmadd(gtu22,kmsub(JacPDstandardNth2gt12,ToReal(2),JacPDstandardNth1gt22),kmadd(gtu33,kmadd(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33),kmul(gtu23,kmul(JacPDstandardNth3gt12,ToReal(4)))))))))))))))))),kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)))); + ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))))); /* If necessary, store only partial vectors after the first iteration */ @@ -625,8 +913,8 @@ extern "C" void ML_BSSN_convertToADMBaseDtLapseShift(CCTK_ARGUMENTS) return; } - const char *groups[] = {"ADMBase::dtlapse","ADMBase::dtshift","grid::coordinates","Grid::coordinates","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"}; - GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertToADMBaseDtLapseShift", 10, groups); + const char *groups[] = {"ADMBase::dtlapse","ADMBase::dtshift","grid::coordinates","Grid::coordinates","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"}; + GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertToADMBaseDtLapseShift", 12, groups); switch(fdOrder) { diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc index db5de18..a957b55 100644 --- a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc +++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc @@ -260,17 +260,32 @@ static void ML_BSSN_convertToADMBaseDtLapseShiftBoundary_Body(cGH const * restri CCTK_REAL_VEC dtalpL = kneg(kmul(kpow(alphaL,harmonicN),kmul(ToReal(harmonicF),kmadd(ksub(AL,trKL),ToReal(LapseACoeff),trKL)))); - CCTK_REAL_VEC dtbetaxL = - kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + CCTK_REAL_VEC dtbetaxL; + CCTK_REAL_VEC dtbetayL; + CCTK_REAL_VEC dtbetazL; - CCTK_REAL_VEC dtbetayL = - kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); - - CCTK_REAL_VEC dtbetazL = - kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + if (harmonicShift) + { + dtbetaxL = ToReal(0); + + dtbetayL = ToReal(0); + + dtbetazL = ToReal(0); + } + else + { + dtbetaxL = + kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + + dtbetayL = + kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + + dtbetazL = + kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + } /* If necessary, store only partial vectors after the first iteration */ diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.cc b/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.cc index a93eb01..0d8065d 100644 --- a/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.cc +++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.cc @@ -245,17 +245,32 @@ static void ML_BSSN_convertToADMBaseFakeDtLapseShift_Body(cGH const * restrict c CCTK_REAL_VEC dtalpL = kneg(kmul(kpow(alphaL,harmonicN),kmul(ToReal(harmonicF),kmadd(ksub(AL,trKL),ToReal(LapseACoeff),trKL)))); - CCTK_REAL_VEC dtbetaxL = - kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + CCTK_REAL_VEC dtbetaxL; + CCTK_REAL_VEC dtbetayL; + CCTK_REAL_VEC dtbetazL; - CCTK_REAL_VEC dtbetayL = - kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); - - CCTK_REAL_VEC dtbetazL = - kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + - ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + if (harmonicShift) + { + dtbetaxL = ToReal(0); + + dtbetayL = ToReal(0); + + dtbetazL = ToReal(0); + } + else + { + dtbetaxL = + kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + + dtbetayL = + kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + + dtbetazL = + kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + + ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))); + } /* If necessary, store only partial vectors after the first iteration */ -- cgit v1.2.3