diff options
Diffstat (limited to 'ML_BSSN/src/ML_BSSN_Dissipation.cc')
| -rw-r--r-- | ML_BSSN/src/ML_BSSN_Dissipation.cc | 1110 |
1 files changed, 1009 insertions, 101 deletions
diff --git a/ML_BSSN/src/ML_BSSN_Dissipation.cc b/ML_BSSN/src/ML_BSSN_Dissipation.cc index 524b6d8..34683db 100644 --- a/ML_BSSN/src/ML_BSSN_Dissipation.cc +++ b/ML_BSSN/src/ML_BSSN_Dissipation.cc @@ -79,7 +79,24 @@ static void ML_BSSN_Dissipation_Body(cGH const * restrict const cctkGH, int cons const char *groups[] = {"ML_BSSN::ML_curv","ML_BSSN::ML_curvrhs","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtlapserhs","ML_BSSN::ML_dtshift","ML_BSSN::ML_dtshiftrhs","ML_BSSN::ML_Gamma","ML_BSSN::ML_Gammarhs","ML_BSSN::ML_lapse","ML_BSSN::ML_lapserhs","ML_BSSN::ML_log_confac","ML_BSSN::ML_log_confacrhs","ML_BSSN::ML_metric","ML_BSSN::ML_metricrhs","ML_BSSN::ML_shift","ML_BSSN::ML_shiftrhs","ML_BSSN::ML_trace_curv","ML_BSSN::ML_trace_curvrhs"}; GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_Dissipation", 18, groups); - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_Dissipation", 3, 3, 3); + switch(fdOrder) + { + case 2: + GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_Dissipation", 2, 2, 2); + break; + + case 4: + GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_Dissipation", 3, 3, 3); + break; + + case 6: + GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_Dissipation", 4, 4, 4); + break; + + case 8: + GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_Dissipation", 5, 5, 5); + break; + } /* Include user-supplied include files */ @@ -107,24 +124,129 @@ static void ML_BSSN_Dissipation_Body(cGH const * restrict const cctkGH, int cons CCTK_REAL_VEC const hdzi = kmul(ToReal(0.5), dzi); /* Initialize predefined quantities */ + CCTK_REAL_VEC const p1o1024dx = kmul(INV(dx),ToReal(0.0009765625)); + CCTK_REAL_VEC const p1o1024dy = kmul(INV(dy),ToReal(0.0009765625)); + CCTK_REAL_VEC const p1o1024dz = kmul(INV(dz),ToReal(0.0009765625)); + CCTK_REAL_VEC const p1o120dx = kmul(INV(dx),ToReal(0.00833333333333333333333333333333)); + CCTK_REAL_VEC const p1o120dy = kmul(INV(dy),ToReal(0.00833333333333333333333333333333)); + CCTK_REAL_VEC const p1o120dz = kmul(INV(dz),ToReal(0.00833333333333333333333333333333)); CCTK_REAL_VEC const p1o12dx = kmul(INV(dx),ToReal(0.0833333333333333333333333333333)); CCTK_REAL_VEC const p1o12dy = kmul(INV(dy),ToReal(0.0833333333333333333333333333333)); CCTK_REAL_VEC const p1o12dz = kmul(INV(dz),ToReal(0.0833333333333333333333333333333)); CCTK_REAL_VEC const p1o144dxdy = kmul(INV(dx),kmul(INV(dy),ToReal(0.00694444444444444444444444444444))); CCTK_REAL_VEC const p1o144dxdz = kmul(INV(dx),kmul(INV(dz),ToReal(0.00694444444444444444444444444444))); CCTK_REAL_VEC const p1o144dydz = kmul(INV(dy),kmul(INV(dz),ToReal(0.00694444444444444444444444444444))); + CCTK_REAL_VEC const p1o1680dx = kmul(INV(dx),ToReal(0.000595238095238095238095238095238)); + CCTK_REAL_VEC const p1o1680dy = kmul(INV(dy),ToReal(0.000595238095238095238095238095238)); + CCTK_REAL_VEC const p1o1680dz = kmul(INV(dz),ToReal(0.000595238095238095238095238095238)); + CCTK_REAL_VEC const p1o16dx = kmul(INV(dx),ToReal(0.0625)); + CCTK_REAL_VEC const p1o16dy = kmul(INV(dy),ToReal(0.0625)); + CCTK_REAL_VEC const p1o16dz = kmul(INV(dz),ToReal(0.0625)); + CCTK_REAL_VEC const p1o180dx2 = kmul(INV(SQR(dx)),ToReal(0.00555555555555555555555555555556)); + CCTK_REAL_VEC const p1o180dy2 = kmul(INV(SQR(dy)),ToReal(0.00555555555555555555555555555556)); + CCTK_REAL_VEC const p1o180dz2 = kmul(INV(SQR(dz)),ToReal(0.00555555555555555555555555555556)); CCTK_REAL_VEC const p1o24dx = kmul(INV(dx),ToReal(0.0416666666666666666666666666667)); CCTK_REAL_VEC const p1o24dy = kmul(INV(dy),ToReal(0.0416666666666666666666666666667)); CCTK_REAL_VEC const p1o24dz = kmul(INV(dz),ToReal(0.0416666666666666666666666666667)); + CCTK_REAL_VEC const p1o256dx = kmul(INV(dx),ToReal(0.00390625)); + CCTK_REAL_VEC const p1o256dy = kmul(INV(dy),ToReal(0.00390625)); + CCTK_REAL_VEC const p1o256dz = kmul(INV(dz),ToReal(0.00390625)); + CCTK_REAL_VEC const p1o2dx = kmul(INV(dx),ToReal(0.5)); + CCTK_REAL_VEC const p1o2dy = kmul(INV(dy),ToReal(0.5)); + CCTK_REAL_VEC const p1o2dz = kmul(INV(dz),ToReal(0.5)); + CCTK_REAL_VEC const p1o3600dxdy = kmul(INV(dx),kmul(INV(dy),ToReal(0.000277777777777777777777777777778))); + CCTK_REAL_VEC const p1o3600dxdz = kmul(INV(dx),kmul(INV(dz),ToReal(0.000277777777777777777777777777778))); + CCTK_REAL_VEC const p1o3600dydz = kmul(INV(dy),kmul(INV(dz),ToReal(0.000277777777777777777777777777778))); + CCTK_REAL_VEC const p1o4dx = kmul(INV(dx),ToReal(0.25)); + CCTK_REAL_VEC const p1o4dxdy = kmul(INV(dx),kmul(INV(dy),ToReal(0.25))); + CCTK_REAL_VEC const p1o4dxdz = kmul(INV(dx),kmul(INV(dz),ToReal(0.25))); + CCTK_REAL_VEC const p1o4dy = kmul(INV(dy),ToReal(0.25)); + CCTK_REAL_VEC const p1o4dydz = kmul(INV(dy),kmul(INV(dz),ToReal(0.25))); + CCTK_REAL_VEC const p1o4dz = kmul(INV(dz),ToReal(0.25)); + CCTK_REAL_VEC const p1o5040dx2 = kmul(INV(SQR(dx)),ToReal(0.000198412698412698412698412698413)); + CCTK_REAL_VEC const p1o5040dy2 = kmul(INV(SQR(dy)),ToReal(0.000198412698412698412698412698413)); + CCTK_REAL_VEC const p1o5040dz2 = kmul(INV(SQR(dz)),ToReal(0.000198412698412698412698412698413)); + CCTK_REAL_VEC const p1o560dx = kmul(INV(dx),ToReal(0.00178571428571428571428571428571)); + CCTK_REAL_VEC const p1o560dy = kmul(INV(dy),ToReal(0.00178571428571428571428571428571)); + CCTK_REAL_VEC const p1o560dz = kmul(INV(dz),ToReal(0.00178571428571428571428571428571)); + CCTK_REAL_VEC const p1o60dx = kmul(INV(dx),ToReal(0.0166666666666666666666666666667)); + CCTK_REAL_VEC const p1o60dy = kmul(INV(dy),ToReal(0.0166666666666666666666666666667)); + CCTK_REAL_VEC const p1o60dz = kmul(INV(dz),ToReal(0.0166666666666666666666666666667)); CCTK_REAL_VEC const p1o64dx = kmul(INV(dx),ToReal(0.015625)); CCTK_REAL_VEC const p1o64dy = kmul(INV(dy),ToReal(0.015625)); CCTK_REAL_VEC const p1o64dz = kmul(INV(dz),ToReal(0.015625)); + CCTK_REAL_VEC const p1o705600dxdy = kmul(INV(dx),kmul(INV(dy),ToReal(1.41723356009070294784580498866e-6))); + CCTK_REAL_VEC const p1o705600dxdz = kmul(INV(dx),kmul(INV(dz),ToReal(1.41723356009070294784580498866e-6))); + CCTK_REAL_VEC const p1o705600dydz = kmul(INV(dy),kmul(INV(dz),ToReal(1.41723356009070294784580498866e-6))); + CCTK_REAL_VEC const p1o840dx = kmul(INV(dx),ToReal(0.00119047619047619047619047619048)); + CCTK_REAL_VEC const p1o840dy = kmul(INV(dy),ToReal(0.00119047619047619047619047619048)); + CCTK_REAL_VEC const p1o840dz = kmul(INV(dz),ToReal(0.00119047619047619047619047619048)); CCTK_REAL_VEC const p1odx = INV(dx); + CCTK_REAL_VEC const p1odx2 = INV(SQR(dx)); CCTK_REAL_VEC const p1ody = INV(dy); + CCTK_REAL_VEC const p1ody2 = INV(SQR(dy)); CCTK_REAL_VEC const p1odz = INV(dz); + CCTK_REAL_VEC const p1odz2 = INV(SQR(dz)); + CCTK_REAL_VEC const pm1o120dx = kmul(INV(dx),ToReal(-0.00833333333333333333333333333333)); + CCTK_REAL_VEC const pm1o120dy = kmul(INV(dy),ToReal(-0.00833333333333333333333333333333)); + CCTK_REAL_VEC const pm1o120dz = kmul(INV(dz),ToReal(-0.00833333333333333333333333333333)); CCTK_REAL_VEC const pm1o12dx2 = kmul(INV(SQR(dx)),ToReal(-0.0833333333333333333333333333333)); CCTK_REAL_VEC const pm1o12dy2 = kmul(INV(SQR(dy)),ToReal(-0.0833333333333333333333333333333)); CCTK_REAL_VEC const pm1o12dz2 = kmul(INV(SQR(dz)),ToReal(-0.0833333333333333333333333333333)); + CCTK_REAL_VEC const pm1o2dx = kmul(INV(dx),ToReal(-0.5)); + CCTK_REAL_VEC const pm1o2dy = kmul(INV(dy),ToReal(-0.5)); + CCTK_REAL_VEC const pm1o2dz = kmul(INV(dz),ToReal(-0.5)); + CCTK_REAL_VEC const pm1o4dx = kmul(INV(dx),ToReal(-0.25)); + CCTK_REAL_VEC const pm1o4dy = kmul(INV(dy),ToReal(-0.25)); + CCTK_REAL_VEC const pm1o4dz = kmul(INV(dz),ToReal(-0.25)); + CCTK_REAL_VEC const pm1o60dx = kmul(INV(dx),ToReal(-0.0166666666666666666666666666667)); + CCTK_REAL_VEC const pm1o60dy = kmul(INV(dy),ToReal(-0.0166666666666666666666666666667)); + CCTK_REAL_VEC const pm1o60dz = kmul(INV(dz),ToReal(-0.0166666666666666666666666666667)); + + /* Jacobian variable pointers */ + bool const use_jacobian = (!CCTK_IsFunctionAliased("MultiPatch_GetMap") || MultiPatch_GetMap(cctkGH) != jacobian_identity_map) + && strlen(jacobian_group) > 0; + if (use_jacobian && strlen(jacobian_derivative_group) == 0) + { + CCTK_WARN (1, "GenericFD::jacobian_group and GenericFD::jacobian_derivative_group must both be set to valid group names"); + } + + CCTK_REAL const *restrict jacobian_ptrs[9]; + if (use_jacobian) GenericFD_GroupDataPointers(cctkGH, jacobian_group, + 9, jacobian_ptrs); + + CCTK_REAL const *restrict const J11 = use_jacobian ? jacobian_ptrs[0] : 0; + CCTK_REAL const *restrict const J12 = use_jacobian ? jacobian_ptrs[1] : 0; + CCTK_REAL const *restrict const J13 = use_jacobian ? jacobian_ptrs[2] : 0; + CCTK_REAL const *restrict const J21 = use_jacobian ? jacobian_ptrs[3] : 0; + CCTK_REAL const *restrict const J22 = use_jacobian ? jacobian_ptrs[4] : 0; + CCTK_REAL const *restrict const J23 = use_jacobian ? jacobian_ptrs[5] : 0; + CCTK_REAL const *restrict const J31 = use_jacobian ? jacobian_ptrs[6] : 0; + CCTK_REAL const *restrict const J32 = use_jacobian ? jacobian_ptrs[7] : 0; + CCTK_REAL const *restrict const J33 = use_jacobian ? jacobian_ptrs[8] : 0; + + CCTK_REAL const *restrict jacobian_derivative_ptrs[18]; + if (use_jacobian) GenericFD_GroupDataPointers(cctkGH, jacobian_derivative_group, + 18, jacobian_derivative_ptrs); + + CCTK_REAL const *restrict const dJ111 = use_jacobian ? jacobian_derivative_ptrs[0] : 0; + CCTK_REAL const *restrict const dJ112 = use_jacobian ? jacobian_derivative_ptrs[1] : 0; + CCTK_REAL const *restrict const dJ113 = use_jacobian ? jacobian_derivative_ptrs[2] : 0; + CCTK_REAL const *restrict const dJ122 = use_jacobian ? jacobian_derivative_ptrs[3] : 0; + CCTK_REAL const *restrict const dJ123 = use_jacobian ? jacobian_derivative_ptrs[4] : 0; + CCTK_REAL const *restrict const dJ133 = use_jacobian ? jacobian_derivative_ptrs[5] : 0; + CCTK_REAL const *restrict const dJ211 = use_jacobian ? jacobian_derivative_ptrs[6] : 0; + CCTK_REAL const *restrict const dJ212 = use_jacobian ? jacobian_derivative_ptrs[7] : 0; + CCTK_REAL const *restrict const dJ213 = use_jacobian ? jacobian_derivative_ptrs[8] : 0; + CCTK_REAL const *restrict const dJ222 = use_jacobian ? jacobian_derivative_ptrs[9] : 0; + CCTK_REAL const *restrict const dJ223 = use_jacobian ? jacobian_derivative_ptrs[10] : 0; + CCTK_REAL const *restrict const dJ233 = use_jacobian ? jacobian_derivative_ptrs[11] : 0; + CCTK_REAL const *restrict const dJ311 = use_jacobian ? jacobian_derivative_ptrs[12] : 0; + CCTK_REAL const *restrict const dJ312 = use_jacobian ? jacobian_derivative_ptrs[13] : 0; + CCTK_REAL const *restrict const dJ313 = use_jacobian ? jacobian_derivative_ptrs[14] : 0; + CCTK_REAL const *restrict const dJ322 = use_jacobian ? jacobian_derivative_ptrs[15] : 0; + CCTK_REAL const *restrict const dJ323 = use_jacobian ? jacobian_derivative_ptrs[16] : 0; + CCTK_REAL const *restrict const dJ333 = use_jacobian ? jacobian_derivative_ptrs[17] : 0; /* Loop over the grid points */ #pragma omp parallel @@ -189,86 +311,872 @@ static void ML_BSSN_Dissipation_Body(cGH const * restrict const cctkGH, int cons CCTK_REAL_VEC Xt3rhsL = vec_load(Xt3rhs[index]); + CCTK_REAL_VEC J11L, J12L, J13L, J21L, J22L, J23L, J31L, J32L, J33L; + + if (use_jacobian) + { + J11L = vec_load(J11[index]); + J12L = vec_load(J12[index]); + J13L = vec_load(J13[index]); + J21L = vec_load(J21[index]); + J22L = vec_load(J22[index]); + J23L = vec_load(J23[index]); + J31L = vec_load(J31[index]); + J32L = vec_load(J32[index]); + J33L = vec_load(J33[index]); + } + /* Include user supplied include files */ /* Precompute derivatives */ - CCTK_REAL_VEC const PDdissipationNth1A = PDdissipationNth1(&A[index]); - CCTK_REAL_VEC const PDdissipationNth2A = PDdissipationNth2(&A[index]); - CCTK_REAL_VEC const PDdissipationNth3A = PDdissipationNth3(&A[index]); - CCTK_REAL_VEC const PDdissipationNth1alpha = PDdissipationNth1(&alpha[index]); - CCTK_REAL_VEC const PDdissipationNth2alpha = PDdissipationNth2(&alpha[index]); - CCTK_REAL_VEC const PDdissipationNth3alpha = PDdissipationNth3(&alpha[index]); - CCTK_REAL_VEC const PDdissipationNth1At11 = PDdissipationNth1(&At11[index]); - CCTK_REAL_VEC const PDdissipationNth2At11 = PDdissipationNth2(&At11[index]); - CCTK_REAL_VEC const PDdissipationNth3At11 = PDdissipationNth3(&At11[index]); - CCTK_REAL_VEC const PDdissipationNth1At12 = PDdissipationNth1(&At12[index]); - CCTK_REAL_VEC const PDdissipationNth2At12 = PDdissipationNth2(&At12[index]); - CCTK_REAL_VEC const PDdissipationNth3At12 = PDdissipationNth3(&At12[index]); - CCTK_REAL_VEC const PDdissipationNth1At13 = PDdissipationNth1(&At13[index]); - CCTK_REAL_VEC const PDdissipationNth2At13 = PDdissipationNth2(&At13[index]); - CCTK_REAL_VEC const PDdissipationNth3At13 = PDdissipationNth3(&At13[index]); - CCTK_REAL_VEC const PDdissipationNth1At22 = PDdissipationNth1(&At22[index]); - CCTK_REAL_VEC const PDdissipationNth2At22 = PDdissipationNth2(&At22[index]); - CCTK_REAL_VEC const PDdissipationNth3At22 = PDdissipationNth3(&At22[index]); - CCTK_REAL_VEC const PDdissipationNth1At23 = PDdissipationNth1(&At23[index]); - CCTK_REAL_VEC const PDdissipationNth2At23 = PDdissipationNth2(&At23[index]); - CCTK_REAL_VEC const PDdissipationNth3At23 = PDdissipationNth3(&At23[index]); - CCTK_REAL_VEC const PDdissipationNth1At33 = PDdissipationNth1(&At33[index]); - CCTK_REAL_VEC const PDdissipationNth2At33 = PDdissipationNth2(&At33[index]); - CCTK_REAL_VEC const PDdissipationNth3At33 = PDdissipationNth3(&At33[index]); - CCTK_REAL_VEC const PDdissipationNth1B1 = PDdissipationNth1(&B1[index]); - CCTK_REAL_VEC const PDdissipationNth2B1 = PDdissipationNth2(&B1[index]); - CCTK_REAL_VEC const PDdissipationNth3B1 = PDdissipationNth3(&B1[index]); - CCTK_REAL_VEC const PDdissipationNth1B2 = PDdissipationNth1(&B2[index]); - CCTK_REAL_VEC const PDdissipationNth2B2 = PDdissipationNth2(&B2[index]); - CCTK_REAL_VEC const PDdissipationNth3B2 = PDdissipationNth3(&B2[index]); - CCTK_REAL_VEC const PDdissipationNth1B3 = PDdissipationNth1(&B3[index]); - CCTK_REAL_VEC const PDdissipationNth2B3 = PDdissipationNth2(&B3[index]); - CCTK_REAL_VEC const PDdissipationNth3B3 = PDdissipationNth3(&B3[index]); - CCTK_REAL_VEC const PDdissipationNth1beta1 = PDdissipationNth1(&beta1[index]); - CCTK_REAL_VEC const PDdissipationNth2beta1 = PDdissipationNth2(&beta1[index]); - CCTK_REAL_VEC const PDdissipationNth3beta1 = PDdissipationNth3(&beta1[index]); - CCTK_REAL_VEC const PDdissipationNth1beta2 = PDdissipationNth1(&beta2[index]); - CCTK_REAL_VEC const PDdissipationNth2beta2 = PDdissipationNth2(&beta2[index]); - CCTK_REAL_VEC const PDdissipationNth3beta2 = PDdissipationNth3(&beta2[index]); - CCTK_REAL_VEC const PDdissipationNth1beta3 = PDdissipationNth1(&beta3[index]); - CCTK_REAL_VEC const PDdissipationNth2beta3 = PDdissipationNth2(&beta3[index]); - CCTK_REAL_VEC const PDdissipationNth3beta3 = PDdissipationNth3(&beta3[index]); - CCTK_REAL_VEC const PDdissipationNth1gt11 = PDdissipationNth1(>11[index]); - CCTK_REAL_VEC const PDdissipationNth2gt11 = PDdissipationNth2(>11[index]); - CCTK_REAL_VEC const PDdissipationNth3gt11 = PDdissipationNth3(>11[index]); - CCTK_REAL_VEC const PDdissipationNth1gt12 = PDdissipationNth1(>12[index]); - CCTK_REAL_VEC const PDdissipationNth2gt12 = PDdissipationNth2(>12[index]); - CCTK_REAL_VEC const PDdissipationNth3gt12 = PDdissipationNth3(>12[index]); - CCTK_REAL_VEC const PDdissipationNth1gt13 = PDdissipationNth1(>13[index]); - CCTK_REAL_VEC const PDdissipationNth2gt13 = PDdissipationNth2(>13[index]); - CCTK_REAL_VEC const PDdissipationNth3gt13 = PDdissipationNth3(>13[index]); - CCTK_REAL_VEC const PDdissipationNth1gt22 = PDdissipationNth1(>22[index]); - CCTK_REAL_VEC const PDdissipationNth2gt22 = PDdissipationNth2(>22[index]); - CCTK_REAL_VEC const PDdissipationNth3gt22 = PDdissipationNth3(>22[index]); - CCTK_REAL_VEC const PDdissipationNth1gt23 = PDdissipationNth1(>23[index]); - CCTK_REAL_VEC const PDdissipationNth2gt23 = PDdissipationNth2(>23[index]); - CCTK_REAL_VEC const PDdissipationNth3gt23 = PDdissipationNth3(>23[index]); - CCTK_REAL_VEC const PDdissipationNth1gt33 = PDdissipationNth1(>33[index]); - CCTK_REAL_VEC const PDdissipationNth2gt33 = PDdissipationNth2(>33[index]); - CCTK_REAL_VEC const PDdissipationNth3gt33 = PDdissipationNth3(>33[index]); - CCTK_REAL_VEC const PDdissipationNth1phi = PDdissipationNth1(&phi[index]); - CCTK_REAL_VEC const PDdissipationNth2phi = PDdissipationNth2(&phi[index]); - CCTK_REAL_VEC const PDdissipationNth3phi = PDdissipationNth3(&phi[index]); - CCTK_REAL_VEC const PDdissipationNth1trK = PDdissipationNth1(&trK[index]); - CCTK_REAL_VEC const PDdissipationNth2trK = PDdissipationNth2(&trK[index]); - CCTK_REAL_VEC const PDdissipationNth3trK = PDdissipationNth3(&trK[index]); - CCTK_REAL_VEC const PDdissipationNth1Xt1 = PDdissipationNth1(&Xt1[index]); - CCTK_REAL_VEC const PDdissipationNth2Xt1 = PDdissipationNth2(&Xt1[index]); - CCTK_REAL_VEC const PDdissipationNth3Xt1 = PDdissipationNth3(&Xt1[index]); - CCTK_REAL_VEC const PDdissipationNth1Xt2 = PDdissipationNth1(&Xt2[index]); - CCTK_REAL_VEC const PDdissipationNth2Xt2 = PDdissipationNth2(&Xt2[index]); - CCTK_REAL_VEC const PDdissipationNth3Xt2 = PDdissipationNth3(&Xt2[index]); - CCTK_REAL_VEC const PDdissipationNth1Xt3 = PDdissipationNth1(&Xt3[index]); - CCTK_REAL_VEC const PDdissipationNth2Xt3 = PDdissipationNth2(&Xt3[index]); - CCTK_REAL_VEC const PDdissipationNth3Xt3 = PDdissipationNth3(&Xt3[index]); + CCTK_REAL_VEC PDdissipationNth1A; + CCTK_REAL_VEC PDdissipationNth2A; + CCTK_REAL_VEC PDdissipationNth3A; + CCTK_REAL_VEC PDdissipationNth1alpha; + CCTK_REAL_VEC PDdissipationNth2alpha; + CCTK_REAL_VEC PDdissipationNth3alpha; + CCTK_REAL_VEC PDdissipationNth1At11; + CCTK_REAL_VEC PDdissipationNth2At11; + CCTK_REAL_VEC PDdissipationNth3At11; + CCTK_REAL_VEC PDdissipationNth1At12; + CCTK_REAL_VEC PDdissipationNth2At12; + CCTK_REAL_VEC PDdissipationNth3At12; + CCTK_REAL_VEC PDdissipationNth1At13; + CCTK_REAL_VEC PDdissipationNth2At13; + CCTK_REAL_VEC PDdissipationNth3At13; + CCTK_REAL_VEC PDdissipationNth1At22; + CCTK_REAL_VEC PDdissipationNth2At22; + CCTK_REAL_VEC PDdissipationNth3At22; + CCTK_REAL_VEC PDdissipationNth1At23; + CCTK_REAL_VEC PDdissipationNth2At23; + CCTK_REAL_VEC PDdissipationNth3At23; + CCTK_REAL_VEC PDdissipationNth1At33; + CCTK_REAL_VEC PDdissipationNth2At33; + CCTK_REAL_VEC PDdissipationNth3At33; + CCTK_REAL_VEC PDdissipationNth1B1; + CCTK_REAL_VEC PDdissipationNth2B1; + CCTK_REAL_VEC PDdissipationNth3B1; + CCTK_REAL_VEC PDdissipationNth1B2; + CCTK_REAL_VEC PDdissipationNth2B2; + CCTK_REAL_VEC PDdissipationNth3B2; + CCTK_REAL_VEC PDdissipationNth1B3; + CCTK_REAL_VEC PDdissipationNth2B3; + CCTK_REAL_VEC PDdissipationNth3B3; + CCTK_REAL_VEC PDdissipationNth1beta1; + CCTK_REAL_VEC PDdissipationNth2beta1; + CCTK_REAL_VEC PDdissipationNth3beta1; + CCTK_REAL_VEC PDdissipationNth1beta2; + CCTK_REAL_VEC PDdissipationNth2beta2; + CCTK_REAL_VEC PDdissipationNth3beta2; + CCTK_REAL_VEC PDdissipationNth1beta3; + CCTK_REAL_VEC PDdissipationNth2beta3; + CCTK_REAL_VEC PDdissipationNth3beta3; + CCTK_REAL_VEC PDdissipationNth1gt11; + CCTK_REAL_VEC PDdissipationNth2gt11; + CCTK_REAL_VEC PDdissipationNth3gt11; + CCTK_REAL_VEC PDdissipationNth1gt12; + CCTK_REAL_VEC PDdissipationNth2gt12; + CCTK_REAL_VEC PDdissipationNth3gt12; + CCTK_REAL_VEC PDdissipationNth1gt13; + CCTK_REAL_VEC PDdissipationNth2gt13; + CCTK_REAL_VEC PDdissipationNth3gt13; + CCTK_REAL_VEC PDdissipationNth1gt22; + CCTK_REAL_VEC PDdissipationNth2gt22; + CCTK_REAL_VEC PDdissipationNth3gt22; + CCTK_REAL_VEC PDdissipationNth1gt23; + CCTK_REAL_VEC PDdissipationNth2gt23; + CCTK_REAL_VEC PDdissipationNth3gt23; + CCTK_REAL_VEC PDdissipationNth1gt33; + CCTK_REAL_VEC PDdissipationNth2gt33; + CCTK_REAL_VEC PDdissipationNth3gt33; + CCTK_REAL_VEC PDdissipationNth1phi; + CCTK_REAL_VEC PDdissipationNth2phi; + CCTK_REAL_VEC PDdissipationNth3phi; + CCTK_REAL_VEC PDdissipationNth1trK; + CCTK_REAL_VEC PDdissipationNth2trK; + CCTK_REAL_VEC PDdissipationNth3trK; + CCTK_REAL_VEC PDdissipationNth1Xt1; + CCTK_REAL_VEC PDdissipationNth2Xt1; + CCTK_REAL_VEC PDdissipationNth3Xt1; + CCTK_REAL_VEC PDdissipationNth1Xt2; + CCTK_REAL_VEC PDdissipationNth2Xt2; + CCTK_REAL_VEC PDdissipationNth3Xt2; + CCTK_REAL_VEC PDdissipationNth1Xt3; + CCTK_REAL_VEC PDdissipationNth2Xt3; + CCTK_REAL_VEC PDdissipationNth3Xt3; + + switch(fdOrder) + { + case 2: + PDdissipationNth1A = PDdissipationNthfdOrder21(&A[index]); + PDdissipationNth2A = PDdissipationNthfdOrder22(&A[index]); + PDdissipationNth3A = PDdissipationNthfdOrder23(&A[index]); + PDdissipationNth1alpha = PDdissipationNthfdOrder21(&alpha[index]); + PDdissipationNth2alpha = PDdissipationNthfdOrder22(&alpha[index]); + PDdissipationNth3alpha = PDdissipationNthfdOrder23(&alpha[index]); + PDdissipationNth1At11 = PDdissipationNthfdOrder21(&At11[index]); + PDdissipationNth2At11 = PDdissipationNthfdOrder22(&At11[index]); + PDdissipationNth3At11 = PDdissipationNthfdOrder23(&At11[index]); + PDdissipationNth1At12 = PDdissipationNthfdOrder21(&At12[index]); + PDdissipationNth2At12 = PDdissipationNthfdOrder22(&At12[index]); + PDdissipationNth3At12 = PDdissipationNthfdOrder23(&At12[index]); + PDdissipationNth1At13 = PDdissipationNthfdOrder21(&At13[index]); + PDdissipationNth2At13 = PDdissipationNthfdOrder22(&At13[index]); + PDdissipationNth3At13 = PDdissipationNthfdOrder23(&At13[index]); + PDdissipationNth1At22 = PDdissipationNthfdOrder21(&At22[index]); + PDdissipationNth2At22 = PDdissipationNthfdOrder22(&At22[index]); + PDdissipationNth3At22 = PDdissipationNthfdOrder23(&At22[index]); + PDdissipationNth1At23 = PDdissipationNthfdOrder21(&At23[index]); + PDdissipationNth2At23 = PDdissipationNthfdOrder22(&At23[index]); + PDdissipationNth3At23 = PDdissipationNthfdOrder23(&At23[index]); + PDdissipationNth1At33 = PDdissipationNthfdOrder21(&At33[index]); + PDdissipationNth2At33 = PDdissipationNthfdOrder22(&At33[index]); + PDdissipationNth3At33 = PDdissipationNthfdOrder23(&At33[index]); + PDdissipationNth1B1 = PDdissipationNthfdOrder21(&B1[index]); + PDdissipationNth2B1 = PDdissipationNthfdOrder22(&B1[index]); + PDdissipationNth3B1 = PDdissipationNthfdOrder23(&B1[index]); + PDdissipationNth1B2 = PDdissipationNthfdOrder21(&B2[index]); + PDdissipationNth2B2 = PDdissipationNthfdOrder22(&B2[index]); + PDdissipationNth3B2 = PDdissipationNthfdOrder23(&B2[index]); + PDdissipationNth1B3 = PDdissipationNthfdOrder21(&B3[index]); + PDdissipationNth2B3 = PDdissipationNthfdOrder22(&B3[index]); + PDdissipationNth3B3 = PDdissipationNthfdOrder23(&B3[index]); + PDdissipationNth1beta1 = PDdissipationNthfdOrder21(&beta1[index]); + PDdissipationNth2beta1 = PDdissipationNthfdOrder22(&beta1[index]); + PDdissipationNth3beta1 = PDdissipationNthfdOrder23(&beta1[index]); + PDdissipationNth1beta2 = PDdissipationNthfdOrder21(&beta2[index]); + PDdissipationNth2beta2 = PDdissipationNthfdOrder22(&beta2[index]); + PDdissipationNth3beta2 = PDdissipationNthfdOrder23(&beta2[index]); + PDdissipationNth1beta3 = PDdissipationNthfdOrder21(&beta3[index]); + PDdissipationNth2beta3 = PDdissipationNthfdOrder22(&beta3[index]); + PDdissipationNth3beta3 = PDdissipationNthfdOrder23(&beta3[index]); + PDdissipationNth1gt11 = PDdissipationNthfdOrder21(>11[index]); + PDdissipationNth2gt11 = PDdissipationNthfdOrder22(>11[index]); + PDdissipationNth3gt11 = PDdissipationNthfdOrder23(>11[index]); + PDdissipationNth1gt12 = PDdissipationNthfdOrder21(>12[index]); + PDdissipationNth2gt12 = PDdissipationNthfdOrder22(>12[index]); + PDdissipationNth3gt12 = PDdissipationNthfdOrder23(>12[index]); + PDdissipationNth1gt13 = PDdissipationNthfdOrder21(>13[index]); + PDdissipationNth2gt13 = PDdissipationNthfdOrder22(>13[index]); + PDdissipationNth3gt13 = PDdissipationNthfdOrder23(>13[index]); + PDdissipationNth1gt22 = PDdissipationNthfdOrder21(>22[index]); + PDdissipationNth2gt22 = PDdissipationNthfdOrder22(>22[index]); + PDdissipationNth3gt22 = PDdissipationNthfdOrder23(>22[index]); + PDdissipationNth1gt23 = PDdissipationNthfdOrder21(>23[index]); + PDdissipationNth2gt23 = PDdissipationNthfdOrder22(>23[index]); + PDdissipationNth3gt23 = PDdissipationNthfdOrder23(>23[index]); + PDdissipationNth1gt33 = PDdissipationNthfdOrder21(>33[index]); + PDdissipationNth2gt33 = PDdissipationNthfdOrder22(>33[index]); + PDdissipationNth3gt33 = PDdissipationNthfdOrder23(>33[index]); + PDdissipationNth1phi = PDdissipationNthfdOrder21(&phi[index]); + PDdissipationNth2phi = PDdissipationNthfdOrder22(&phi[index]); + PDdissipationNth3phi = PDdissipationNthfdOrder23(&phi[index]); + PDdissipationNth1trK = PDdissipationNthfdOrder21(&trK[index]); + PDdissipationNth2trK = PDdissipationNthfdOrder22(&trK[index]); + PDdissipationNth3trK = PDdissipationNthfdOrder23(&trK[index]); + PDdissipationNth1Xt1 = PDdissipationNthfdOrder21(&Xt1[index]); + PDdissipationNth2Xt1 = PDdissipationNthfdOrder22(&Xt1[index]); + PDdissipationNth3Xt1 = PDdissipationNthfdOrder23(&Xt1[index]); + PDdissipationNth1Xt2 = PDdissipationNthfdOrder21(&Xt2[index]); + PDdissipationNth2Xt2 = PDdissipationNthfdOrder22(&Xt2[index]); + PDdissipationNth3Xt2 = PDdissipationNthfdOrder23(&Xt2[index]); + PDdissipationNth1Xt3 = PDdissipationNthfdOrder21(&Xt3[index]); + PDdissipationNth2Xt3 = PDdissipationNthfdOrder22(&Xt3[index]); + PDdissipationNth3Xt3 = PDdissipationNthfdOrder23(&Xt3[index]); + break; + + case 4: + PDdissipationNth1A = PDdissipationNthfdOrder41(&A[index]); + PDdissipationNth2A = PDdissipationNthfdOrder42(&A[index]); + PDdissipationNth3A = PDdissipationNthfdOrder43(&A[index]); + PDdissipationNth1alpha = PDdissipationNthfdOrder41(&alpha[index]); + PDdissipationNth2alpha = PDdissipationNthfdOrder42(&alpha[index]); + PDdissipationNth3alpha = PDdissipationNthfdOrder43(&alpha[index]); + PDdissipationNth1At11 = PDdissipationNthfdOrder41(&At11[index]); + PDdissipationNth2At11 = PDdissipationNthfdOrder42(&At11[index]); + PDdissipationNth3At11 = PDdissipationNthfdOrder43(&At11[index]); + PDdissipationNth1At12 = PDdissipationNthfdOrder41(&At12[index]); + PDdissipationNth2At12 = PDdissipationNthfdOrder42(&At12[index]); + PDdissipationNth3At12 = PDdissipationNthfdOrder43(&At12[index]); + PDdissipationNth1At13 = PDdissipationNthfdOrder41(&At13[index]); + PDdissipationNth2At13 = PDdissipationNthfdOrder42(&At13[index]); + PDdissipationNth3At13 = PDdissipationNthfdOrder43(&At13[index]); + PDdissipationNth1At22 = PDdissipationNthfdOrder41(&At22[index]); + PDdissipationNth2At22 = PDdissipationNthfdOrder42(&At22[index]); + PDdissipationNth3At22 = PDdissipationNthfdOrder43(&At22[index]); + PDdissipationNth1At23 = PDdissipationNthfdOrder41(&At23[index]); + PDdissipationNth2At23 = PDdissipationNthfdOrder42(&At23[index]); + PDdissipationNth3At23 = PDdissipationNthfdOrder43(&At23[index]); + PDdissipationNth1At33 = PDdissipationNthfdOrder41(&At33[index]); + PDdissipationNth2At33 = PDdissipationNthfdOrder42(&At33[index]); + PDdissipationNth3At33 = PDdissipationNthfdOrder43(&At33[index]); + PDdissipationNth1B1 = PDdissipationNthfdOrder41(&B1[index]); + PDdissipationNth2B1 = PDdissipationNthfdOrder42(&B1[index]); + PDdissipationNth3B1 = PDdissipationNthfdOrder43(&B1[index]); + PDdissipationNth1B2 = PDdissipationNthfdOrder41(&B2[index]); + PDdissipationNth2B2 = PDdissipationNthfdOrder42(&B2[index]); + PDdissipationNth3B2 = PDdissipationNthfdOrder43(&B2[index]); + PDdissipationNth1B3 = PDdissipationNthfdOrder41(&B3[index]); + PDdissipationNth2B3 = PDdissipationNthfdOrder42(&B3[index]); + PDdissipationNth3B3 = PDdissipationNthfdOrder43(&B3[index]); + PDdissipationNth1beta1 = PDdissipationNthfdOrder41(&beta1[index]); + PDdissipationNth2beta1 = PDdissipationNthfdOrder42(&beta1[index]); + PDdissipationNth3beta1 = PDdissipationNthfdOrder43(&beta1[index]); + PDdissipationNth1beta2 = PDdissipationNthfdOrder41(&beta2[index]); + PDdissipationNth2beta2 = PDdissipationNthfdOrder42(&beta2[index]); + PDdissipationNth3beta2 = PDdissipationNthfdOrder43(&beta2[index]); + PDdissipationNth1beta3 = PDdissipationNthfdOrder41(&beta3[index]); + PDdissipationNth2beta3 = PDdissipationNthfdOrder42(&beta3[index]); + PDdissipationNth3beta3 = PDdissipationNthfdOrder43(&beta3[index]); + PDdissipationNth1gt11 = PDdissipationNthfdOrder41(>11[index]); + PDdissipationNth2gt11 = PDdissipationNthfdOrder42(>11[index]); + PDdissipationNth3gt11 = PDdissipationNthfdOrder43(>11[index]); + PDdissipationNth1gt12 = PDdissipationNthfdOrder41(>12[index]); + PDdissipationNth2gt12 = PDdissipationNthfdOrder42(>12[index]); + PDdissipationNth3gt12 = PDdissipationNthfdOrder43(>12[index]); + PDdissipationNth1gt13 = PDdissipationNthfdOrder41(>13[index]); + PDdissipationNth2gt13 = PDdissipationNthfdOrder42(>13[index]); + PDdissipationNth3gt13 = PDdissipationNthfdOrder43(>13[index]); + PDdissipationNth1gt22 = PDdissipationNthfdOrder41(>22[index]); + PDdissipationNth2gt22 = PDdissipationNthfdOrder42(>22[index]); + PDdissipationNth3gt22 = PDdissipationNthfdOrder43(>22[index]); + PDdissipationNth1gt23 = PDdissipationNthfdOrder41(>23[index]); + PDdissipationNth2gt23 = PDdissipationNthfdOrder42(>23[index]); + PDdissipationNth3gt23 = PDdissipationNthfdOrder43(>23[index]); + PDdissipationNth1gt33 = PDdissipationNthfdOrder41(>33[index]); + PDdissipationNth2gt33 = PDdissipationNthfdOrder42(>33[index]); + PDdissipationNth3gt33 = PDdissipationNthfdOrder43(>33[index]); + PDdissipationNth1phi = PDdissipationNthfdOrder41(&phi[index]); + PDdissipationNth2phi = PDdissipationNthfdOrder42(&phi[index]); + PDdissipationNth3phi = PDdissipationNthfdOrder43(&phi[index]); + PDdissipationNth1trK = PDdissipationNthfdOrder41(&trK[index]); + PDdissipationNth2trK = PDdissipationNthfdOrder42(&trK[index]); + PDdissipationNth3trK = PDdissipationNthfdOrder43(&trK[index]); + PDdissipationNth1Xt1 = PDdissipationNthfdOrder41(&Xt1[index]); + PDdissipationNth2Xt1 = PDdissipationNthfdOrder42(&Xt1[index]); + PDdissipationNth3Xt1 = PDdissipationNthfdOrder43(&Xt1[index]); + PDdissipationNth1Xt2 = PDdissipationNthfdOrder41(&Xt2[index]); + PDdissipationNth2Xt2 = PDdissipationNthfdOrder42(&Xt2[index]); + PDdissipationNth3Xt2 = PDdissipationNthfdOrder43(&Xt2[index]); + PDdissipationNth1Xt3 = PDdissipationNthfdOrder41(&Xt3[index]); + PDdissipationNth2Xt3 = PDdissipationNthfdOrder42(&Xt3[index]); + PDdissipationNth3Xt3 = PDdissipationNthfdOrder43(&Xt3[index]); + break; + + case 6: + PDdissipationNth1A = PDdissipationNthfdOrder61(&A[index]); + PDdissipationNth2A = PDdissipationNthfdOrder62(&A[index]); + PDdissipationNth3A = PDdissipationNthfdOrder63(&A[index]); + PDdissipationNth1alpha = PDdissipationNthfdOrder61(&alpha[index]); + PDdissipationNth2alpha = PDdissipationNthfdOrder62(&alpha[index]); + PDdissipationNth3alpha = PDdissipationNthfdOrder63(&alpha[index]); + PDdissipationNth1At11 = PDdissipationNthfdOrder61(&At11[index]); + PDdissipationNth2At11 = PDdissipationNthfdOrder62(&At11[index]); + PDdissipationNth3At11 = PDdissipationNthfdOrder63(&At11[index]); + PDdissipationNth1At12 = PDdissipationNthfdOrder61(&At12[index]); + PDdissipationNth2At12 = PDdissipationNthfdOrder62(&At12[index]); + PDdissipationNth3At12 = PDdissipationNthfdOrder63(&At12[index]); + PDdissipationNth1At13 = PDdissipationNthfdOrder61(&At13[index]); + PDdissipationNth2At13 = PDdissipationNthfdOrder62(&At13[index]); + PDdissipationNth3At13 = PDdissipationNthfdOrder63(&At13[index]); + PDdissipationNth1At22 = PDdissipationNthfdOrder61(&At22[index]); + PDdissipationNth2At22 = PDdissipationNthfdOrder62(&At22[index]); + PDdissipationNth3At22 = PDdissipationNthfdOrder63(&At22[index]); + PDdissipationNth1At23 = PDdissipationNthfdOrder61(&At23[index]); + PDdissipationNth2At23 = PDdissipationNthfdOrder62(&At23[index]); + PDdissipationNth3At23 = PDdissipationNthfdOrder63(&At23[index]); + PDdissipationNth1At33 = PDdissipationNthfdOrder61(&At33[index]); + PDdissipationNth2At33 = PDdissipationNthfdOrder62(&At33[index]); + PDdissipationNth3At33 = PDdissipationNthfdOrder63(&At33[index]); + PDdissipationNth1B1 = PDdissipationNthfdOrder61(&B1[index]); + PDdissipationNth2B1 = PDdissipationNthfdOrder62(&B1[index]); + PDdissipationNth3B1 = PDdissipationNthfdOrder63(&B1[index]); + PDdissipationNth1B2 = PDdissipationNthfdOrder61(&B2[index]); + PDdissipationNth2B2 = PDdissipationNthfdOrder62(&B2[index]); + PDdissipationNth3B2 = PDdissipationNthfdOrder63(&B2[index]); + PDdissipationNth1B3 = PDdissipationNthfdOrder61(&B3[index]); + PDdissipationNth2B3 = PDdissipationNthfdOrder62(&B3[index]); + PDdissipationNth3B3 = PDdissipationNthfdOrder63(&B3[index]); + PDdissipationNth1beta1 = PDdissipationNthfdOrder61(&beta1[index]); + PDdissipationNth2beta1 = PDdissipationNthfdOrder62(&beta1[index]); + PDdissipationNth3beta1 = PDdissipationNthfdOrder63(&beta1[index]); + PDdissipationNth1beta2 = PDdissipationNthfdOrder61(&beta2[index]); + PDdissipationNth2beta2 = PDdissipationNthfdOrder62(&beta2[index]); + PDdissipationNth3beta2 = PDdissipationNthfdOrder63(&beta2[index]); + PDdissipationNth1beta3 = PDdissipationNthfdOrder61(&beta3[index]); + PDdissipationNth2beta3 = PDdissipationNthfdOrder62(&beta3[index]); + PDdissipationNth3beta3 = PDdissipationNthfdOrder63(&beta3[index]); + PDdissipationNth1gt11 = PDdissipationNthfdOrder61(>11[index]); + PDdissipationNth2gt11 = PDdissipationNthfdOrder62(>11[index]); + PDdissipationNth3gt11 = PDdissipationNthfdOrder63(>11[index]); + PDdissipationNth1gt12 = PDdissipationNthfdOrder61(>12[index]); + PDdissipationNth2gt12 = PDdissipationNthfdOrder62(>12[index]); + PDdissipationNth3gt12 = PDdissipationNthfdOrder63(>12[index]); + PDdissipationNth1gt13 = PDdissipationNthfdOrder61(>13[index]); + PDdissipationNth2gt13 = PDdissipationNthfdOrder62(>13[index]); + PDdissipationNth3gt13 = PDdissipationNthfdOrder63(>13[index]); + PDdissipationNth1gt22 = PDdissipationNthfdOrder61(>22[index]); + PDdissipationNth2gt22 = PDdissipationNthfdOrder62(>22[index]); + PDdissipationNth3gt22 = PDdissipationNthfdOrder63(>22[index]); + PDdissipationNth1gt23 = PDdissipationNthfdOrder61(>23[index]); + PDdissipationNth2gt23 = PDdissipationNthfdOrder62(>23[index]); + PDdissipationNth3gt23 = PDdissipationNthfdOrder63(>23[index]); + PDdissipationNth1gt33 = PDdissipationNthfdOrder61(>33[index]); + PDdissipationNth2gt33 = PDdissipationNthfdOrder62(>33[index]); + PDdissipationNth3gt33 = PDdissipationNthfdOrder63(>33[index]); + PDdissipationNth1phi = PDdissipationNthfdOrder61(&phi[index]); + PDdissipationNth2phi = PDdissipationNthfdOrder62(&phi[index]); + PDdissipationNth3phi = PDdissipationNthfdOrder63(&phi[index]); + PDdissipationNth1trK = PDdissipationNthfdOrder61(&trK[index]); + PDdissipationNth2trK = PDdissipationNthfdOrder62(&trK[index]); + PDdissipationNth3trK = PDdissipationNthfdOrder63(&trK[index]); + PDdissipationNth1Xt1 = PDdissipationNthfdOrder61(&Xt1[index]); + PDdissipationNth2Xt1 = PDdissipationNthfdOrder62(&Xt1[index]); + PDdissipationNth3Xt1 = PDdissipationNthfdOrder63(&Xt1[index]); + PDdissipationNth1Xt2 = PDdissipationNthfdOrder61(&Xt2[index]); + PDdissipationNth2Xt2 = PDdissipationNthfdOrder62(&Xt2[index]); + PDdissipationNth3Xt2 = PDdissipationNthfdOrder63(&Xt2[index]); + PDdissipationNth1Xt3 = PDdissipationNthfdOrder61(&Xt3[index]); + PDdissipationNth2Xt3 = PDdissipationNthfdOrder62(&Xt3[index]); + PDdissipationNth3Xt3 = PDdissipationNthfdOrder63(&Xt3[index]); + break; + + case 8: + PDdissipationNth1A = PDdissipationNthfdOrder81(&A[index]); + PDdissipationNth2A = PDdissipationNthfdOrder82(&A[index]); + PDdissipationNth3A = PDdissipationNthfdOrder83(&A[index]); + PDdissipationNth1alpha = PDdissipationNthfdOrder81(&alpha[index]); + PDdissipationNth2alpha = PDdissipationNthfdOrder82(&alpha[index]); + PDdissipationNth3alpha = PDdissipationNthfdOrder83(&alpha[index]); + PDdissipationNth1At11 = PDdissipationNthfdOrder81(&At11[index]); + PDdissipationNth2At11 = PDdissipationNthfdOrder82(&At11[index]); + PDdissipationNth3At11 = PDdissipationNthfdOrder83(&At11[index]); + PDdissipationNth1At12 = PDdissipationNthfdOrder81(&At12[index]); + PDdissipationNth2At12 = PDdissipationNthfdOrder82(&At12[index]); + PDdissipationNth3At12 = PDdissipationNthfdOrder83(&At12[index]); + PDdissipationNth1At13 = PDdissipationNthfdOrder81(&At13[index]); + PDdissipationNth2At13 = PDdissipationNthfdOrder82(&At13[index]); + PDdissipationNth3At13 = PDdissipationNthfdOrder83(&At13[index]); + PDdissipationNth1At22 = PDdissipationNthfdOrder81(&At22[index]); + PDdissipationNth2At22 = PDdissipationNthfdOrder82(&At22[index]); + PDdissipationNth3At22 = PDdissipationNthfdOrder83(&At22[index]); + PDdissipationNth1At23 = PDdissipationNthfdOrder81(&At23[index]); + PDdissipationNth2At23 = PDdissipationNthfdOrder82(&At23[index]); + PDdissipationNth3At23 = PDdissipationNthfdOrder83(&At23[index]); + PDdissipationNth1At33 = PDdissipationNthfdOrder81(&At33[index]); + PDdissipationNth2At33 = PDdissipationNthfdOrder82(&At33[index]); + PDdissipationNth3At33 = PDdissipationNthfdOrder83(&At33[index]); + PDdissipationNth1B1 = PDdissipationNthfdOrder81(&B1[index]); + PDdissipationNth2B1 = PDdissipationNthfdOrder82(&B1[index]); + PDdissipationNth3B1 = PDdissipationNthfdOrder83(&B1[index]); + PDdissipationNth1B2 = PDdissipationNthfdOrder81(&B2[index]); + PDdissipationNth2B2 = PDdissipationNthfdOrder82(&B2[index]); + PDdissipationNth3B2 = PDdissipationNthfdOrder83(&B2[index]); + PDdissipationNth1B3 = PDdissipationNthfdOrder81(&B3[index]); + PDdissipationNth2B3 = PDdissipationNthfdOrder82(&B3[index]); + PDdissipationNth3B3 = PDdissipationNthfdOrder83(&B3[index]); + PDdissipationNth1beta1 = PDdissipationNthfdOrder81(&beta1[index]); + PDdissipationNth2beta1 = PDdissipationNthfdOrder82(&beta1[index]); + PDdissipationNth3beta1 = PDdissipationNthfdOrder83(&beta1[index]); + PDdissipationNth1beta2 = PDdissipationNthfdOrder81(&beta2[index]); + PDdissipationNth2beta2 = PDdissipationNthfdOrder82(&beta2[index]); + PDdissipationNth3beta2 = PDdissipationNthfdOrder83(&beta2[index]); + PDdissipationNth1beta3 = PDdissipationNthfdOrder81(&beta3[index]); + PDdissipationNth2beta3 = PDdissipationNthfdOrder82(&beta3[index]); + PDdissipationNth3beta3 = PDdissipationNthfdOrder83(&beta3[index]); + PDdissipationNth1gt11 = PDdissipationNthfdOrder81(>11[index]); + PDdissipationNth2gt11 = PDdissipationNthfdOrder82(>11[index]); + PDdissipationNth3gt11 = PDdissipationNthfdOrder83(>11[index]); + PDdissipationNth1gt12 = PDdissipationNthfdOrder81(>12[index]); + PDdissipationNth2gt12 = PDdissipationNthfdOrder82(>12[index]); + PDdissipationNth3gt12 = PDdissipationNthfdOrder83(>12[index]); + PDdissipationNth1gt13 = PDdissipationNthfdOrder81(>13[index]); + PDdissipationNth2gt13 = PDdissipationNthfdOrder82(>13[index]); + PDdissipationNth3gt13 = PDdissipationNthfdOrder83(>13[index]); + PDdissipationNth1gt22 = PDdissipationNthfdOrder81(>22[index]); + PDdissipationNth2gt22 = PDdissipationNthfdOrder82(>22[index]); + PDdissipationNth3gt22 = PDdissipationNthfdOrder83(>22[index]); + PDdissipationNth1gt23 = PDdissipationNthfdOrder81(>23[index]); + PDdissipationNth2gt23 = PDdissipationNthfdOrder82(>23[index]); + PDdissipationNth3gt23 = PDdissipationNthfdOrder83(>23[index]); + PDdissipationNth1gt33 = PDdissipationNthfdOrder81(>33[index]); + PDdissipationNth2gt33 = PDdissipationNthfdOrder82(>33[index]); + PDdissipationNth3gt33 = PDdissipationNthfdOrder83(>33[index]); + PDdissipationNth1phi = PDdissipationNthfdOrder81(&phi[index]); + PDdissipationNth2phi = PDdissipationNthfdOrder82(&phi[index]); + PDdissipationNth3phi = PDdissipationNthfdOrder83(&phi[index]); + PDdissipationNth1trK = PDdissipationNthfdOrder81(&trK[index]); + PDdissipationNth2trK = PDdissipationNthfdOrder82(&trK[index]); + PDdissipationNth3trK = PDdissipationNthfdOrder83(&trK[index]); + PDdissipationNth1Xt1 = PDdissipationNthfdOrder81(&Xt1[index]); + PDdissipationNth2Xt1 = PDdissipationNthfdOrder82(&Xt1[index]); + PDdissipationNth3Xt1 = PDdissipationNthfdOrder83(&Xt1[index]); + PDdissipationNth1Xt2 = PDdissipationNthfdOrder81(&Xt2[index]); + PDdissipationNth2Xt2 = PDdissipationNthfdOrder82(&Xt2[index]); + PDdissipationNth3Xt2 = PDdissipationNthfdOrder83(&Xt2[index]); + PDdissipationNth1Xt3 = PDdissipationNthfdOrder81(&Xt3[index]); + PDdissipationNth2Xt3 = PDdissipationNthfdOrder82(&Xt3[index]); + PDdissipationNth3Xt3 = PDdissipationNthfdOrder83(&Xt3[index]); + break; + } /* Calculate temporaries and grid functions */ + CCTK_REAL_VEC JacPDdissipationNth1A; + CCTK_REAL_VEC JacPDdissipationNth1alpha; + CCTK_REAL_VEC JacPDdissipationNth1At11; + CCTK_REAL_VEC JacPDdissipationNth1At12; + CCTK_REAL_VEC JacPDdissipationNth1At13; + CCTK_REAL_VEC JacPDdissipationNth1At22; + CCTK_REAL_VEC JacPDdissipationNth1At23; + CCTK_REAL_VEC JacPDdissipationNth1At33; + CCTK_REAL_VEC JacPDdissipationNth1B1; + CCTK_REAL_VEC JacPDdissipationNth1B2; + CCTK_REAL_VEC JacPDdissipationNth1B3; + CCTK_REAL_VEC JacPDdissipationNth1beta1; + CCTK_REAL_VEC JacPDdissipationNth1beta2; + CCTK_REAL_VEC JacPDdissipationNth1beta3; + CCTK_REAL_VEC JacPDdissipationNth1gt11; + CCTK_REAL_VEC JacPDdissipationNth1gt12; + CCTK_REAL_VEC JacPDdissipationNth1gt13; + CCTK_REAL_VEC JacPDdissipationNth1gt22; + CCTK_REAL_VEC JacPDdissipationNth1gt23; + CCTK_REAL_VEC JacPDdissipationNth1gt33; + CCTK_REAL_VEC JacPDdissipationNth1phi; + CCTK_REAL_VEC JacPDdissipationNth1trK; + CCTK_REAL_VEC JacPDdissipationNth1Xt1; + CCTK_REAL_VEC JacPDdissipationNth1Xt2; + CCTK_REAL_VEC JacPDdissipationNth1Xt3; + CCTK_REAL_VEC JacPDdissipationNth2A; + CCTK_REAL_VEC JacPDdissipationNth2alpha; + CCTK_REAL_VEC JacPDdissipationNth2At11; + CCTK_REAL_VEC JacPDdissipationNth2At12; + CCTK_REAL_VEC JacPDdissipationNth2At13; + CCTK_REAL_VEC JacPDdissipationNth2At22; + CCTK_REAL_VEC JacPDdissipationNth2At23; + CCTK_REAL_VEC JacPDdissipationNth2At33; + CCTK_REAL_VEC JacPDdissipationNth2B1; + CCTK_REAL_VEC JacPDdissipationNth2B2; + CCTK_REAL_VEC JacPDdissipationNth2B3; + CCTK_REAL_VEC JacPDdissipationNth2beta1; + CCTK_REAL_VEC JacPDdissipationNth2beta2; + CCTK_REAL_VEC JacPDdissipationNth2beta3; + CCTK_REAL_VEC JacPDdissipationNth2gt11; + CCTK_REAL_VEC JacPDdissipationNth2gt12; + CCTK_REAL_VEC JacPDdissipationNth2gt13; + CCTK_REAL_VEC JacPDdissipationNth2gt22; + CCTK_REAL_VEC JacPDdissipationNth2gt23; + CCTK_REAL_VEC JacPDdissipationNth2gt33; + CCTK_REAL_VEC JacPDdissipationNth2phi; + CCTK_REAL_VEC JacPDdissipationNth2trK; + CCTK_REAL_VEC JacPDdissipationNth2Xt1; + CCTK_REAL_VEC JacPDdissipationNth2Xt2; + CCTK_REAL_VEC JacPDdissipationNth2Xt3; + CCTK_REAL_VEC JacPDdissipationNth3A; + CCTK_REAL_VEC JacPDdissipationNth3alpha; + CCTK_REAL_VEC JacPDdissipationNth3At11; + CCTK_REAL_VEC JacPDdissipationNth3At12; + CCTK_REAL_VEC JacPDdissipationNth3At13; + CCTK_REAL_VEC JacPDdissipationNth3At22; + CCTK_REAL_VEC JacPDdissipationNth3At23; + CCTK_REAL_VEC JacPDdissipationNth3At33; + CCTK_REAL_VEC JacPDdissipationNth3B1; + CCTK_REAL_VEC JacPDdissipationNth3B2; + CCTK_REAL_VEC JacPDdissipationNth3B3; + CCTK_REAL_VEC JacPDdissipationNth3beta1; + CCTK_REAL_VEC JacPDdissipationNth3beta2; + CCTK_REAL_VEC JacPDdissipationNth3beta3; + CCTK_REAL_VEC JacPDdissipationNth3gt11; + CCTK_REAL_VEC JacPDdissipationNth3gt12; + CCTK_REAL_VEC JacPDdissipationNth3gt13; + CCTK_REAL_VEC JacPDdissipationNth3gt22; + CCTK_REAL_VEC JacPDdissipationNth3gt23; + CCTK_REAL_VEC JacPDdissipationNth3gt33; + CCTK_REAL_VEC JacPDdissipationNth3phi; + CCTK_REAL_VEC JacPDdissipationNth3trK; + CCTK_REAL_VEC JacPDdissipationNth3Xt1; + CCTK_REAL_VEC JacPDdissipationNth3Xt2; + CCTK_REAL_VEC JacPDdissipationNth3Xt3; + + if (use_jacobian) + { + JacPDdissipationNth1A = + kmadd(J11L,PDdissipationNth1A,kmadd(J21L,PDdissipationNth2A,kmul(J31L,PDdissipationNth3A))); + + JacPDdissipationNth1alpha = + kmadd(J11L,PDdissipationNth1alpha,kmadd(J21L,PDdissipationNth2alpha,kmul(J31L,PDdissipationNth3alpha))); + + JacPDdissipationNth1At11 = + kmadd(J11L,PDdissipationNth1At11,kmadd(J21L,PDdissipationNth2At11,kmul(J31L,PDdissipationNth3At11))); + + JacPDdissipationNth1At12 = + kmadd(J11L,PDdissipationNth1At12,kmadd(J21L,PDdissipationNth2At12,kmul(J31L,PDdissipationNth3At12))); + + JacPDdissipationNth1At13 = + kmadd(J11L,PDdissipationNth1At13,kmadd(J21L,PDdissipationNth2At13,kmul(J31L,PDdissipationNth3At13))); + + JacPDdissipationNth1At22 = + kmadd(J11L,PDdissipationNth1At22,kmadd(J21L,PDdissipationNth2At22,kmul(J31L,PDdissipationNth3At22))); + + JacPDdissipationNth1At23 = + kmadd(J11L,PDdissipationNth1At23,kmadd(J21L,PDdissipationNth2At23,kmul(J31L,PDdissipationNth3At23))); + + JacPDdissipationNth1At33 = + kmadd(J11L,PDdissipationNth1At33,kmadd(J21L,PDdissipationNth2At33,kmul(J31L,PDdissipationNth3At33))); + + JacPDdissipationNth1B1 = + kmadd(J11L,PDdissipationNth1B1,kmadd(J21L,PDdissipationNth2B1,kmul(J31L,PDdissipationNth3B1))); + + JacPDdissipationNth1B2 = + kmadd(J11L,PDdissipationNth1B2,kmadd(J21L,PDdissipationNth2B2,kmul(J31L,PDdissipationNth3B2))); + + JacPDdissipationNth1B3 = + kmadd(J11L,PDdissipationNth1B3,kmadd(J21L,PDdissipationNth2B3,kmul(J31L,PDdissipationNth3B3))); + + JacPDdissipationNth1beta1 = + kmadd(J11L,PDdissipationNth1beta1,kmadd(J21L,PDdissipationNth2beta1,kmul(J31L,PDdissipationNth3beta1))); + + JacPDdissipationNth1beta2 = + kmadd(J11L,PDdissipationNth1beta2,kmadd(J21L,PDdissipationNth2beta2,kmul(J31L,PDdissipationNth3beta2))); + + JacPDdissipationNth1beta3 = + kmadd(J11L,PDdissipationNth1beta3,kmadd(J21L,PDdissipationNth2beta3,kmul(J31L,PDdissipationNth3beta3))); + + JacPDdissipationNth1gt11 = + kmadd(J11L,PDdissipationNth1gt11,kmadd(J21L,PDdissipationNth2gt11,kmul(J31L,PDdissipationNth3gt11))); + + JacPDdissipationNth1gt12 = + kmadd(J11L,PDdissipationNth1gt12,kmadd(J21L,PDdissipationNth2gt12,kmul(J31L,PDdissipationNth3gt12))); + + JacPDdissipationNth1gt13 = + kmadd(J11L,PDdissipationNth1gt13,kmadd(J21L,PDdissipationNth2gt13,kmul(J31L,PDdissipationNth3gt13))); + + JacPDdissipationNth1gt22 = + kmadd(J11L,PDdissipationNth1gt22,kmadd(J21L,PDdissipationNth2gt22,kmul(J31L,PDdissipationNth3gt22))); + + JacPDdissipationNth1gt23 = + kmadd(J11L,PDdissipationNth1gt23,kmadd(J21L,PDdissipationNth2gt23,kmul(J31L,PDdissipationNth3gt23))); + + JacPDdissipationNth1gt33 = + kmadd(J11L,PDdissipationNth1gt33,kmadd(J21L,PDdissipationNth2gt33,kmul(J31L,PDdissipationNth3gt33))); + + JacPDdissipationNth1phi = + kmadd(J11L,PDdissipationNth1phi,kmadd(J21L,PDdissipationNth2phi,kmul(J31L,PDdissipationNth3phi))); + + JacPDdissipationNth1trK = + kmadd(J11L,PDdissipationNth1trK,kmadd(J21L,PDdissipationNth2trK,kmul(J31L,PDdissipationNth3trK))); + + JacPDdissipationNth1Xt1 = + kmadd(J11L,PDdissipationNth1Xt1,kmadd(J21L,PDdissipationNth2Xt1,kmul(J31L,PDdissipationNth3Xt1))); + + JacPDdissipationNth1Xt2 = + kmadd(J11L,PDdissipationNth1Xt2,kmadd(J21L,PDdissipationNth2Xt2,kmul(J31L,PDdissipationNth3Xt2))); + + JacPDdissipationNth1Xt3 = + kmadd(J11L,PDdissipationNth1Xt3,kmadd(J21L,PDdissipationNth2Xt3,kmul(J31L,PDdissipationNth3Xt3))); + + JacPDdissipationNth2A = + kmadd(J12L,PDdissipationNth1A,kmadd(J22L,PDdissipationNth2A,kmul(J32L,PDdissipationNth3A))); + + JacPDdissipationNth2alpha = + kmadd(J12L,PDdissipationNth1alpha,kmadd(J22L,PDdissipationNth2alpha,kmul(J32L,PDdissipationNth3alpha))); + + JacPDdissipationNth2At11 = + kmadd(J12L,PDdissipationNth1At11,kmadd(J22L,PDdissipationNth2At11,kmul(J32L,PDdissipationNth3At11))); + + JacPDdissipationNth2At12 = + kmadd(J12L,PDdissipationNth1At12,kmadd(J22L,PDdissipationNth2At12,kmul(J32L,PDdissipationNth3At12))); + + JacPDdissipationNth2At13 = + kmadd(J12L,PDdissipationNth1At13,kmadd(J22L,PDdissipationNth2At13,kmul(J32L,PDdissipationNth3At13))); + + JacPDdissipationNth2At22 = + kmadd(J12L,PDdissipationNth1At22,kmadd(J22L,PDdissipationNth2At22,kmul(J32L,PDdissipationNth3At22))); + + JacPDdissipationNth2At23 = + kmadd(J12L,PDdissipationNth1At23,kmadd(J22L,PDdissipationNth2At23,kmul(J32L,PDdissipationNth3At23))); + + JacPDdissipationNth2At33 = + kmadd(J12L,PDdissipationNth1At33,kmadd(J22L,PDdissipationNth2At33,kmul(J32L,PDdissipationNth3At33))); + + JacPDdissipationNth2B1 = + kmadd(J12L,PDdissipationNth1B1,kmadd(J22L,PDdissipationNth2B1,kmul(J32L,PDdissipationNth3B1))); + + JacPDdissipationNth2B2 = + kmadd(J12L,PDdissipationNth1B2,kmadd(J22L,PDdissipationNth2B2,kmul(J32L,PDdissipationNth3B2))); + + JacPDdissipationNth2B3 = + kmadd(J12L,PDdissipationNth1B3,kmadd(J22L,PDdissipationNth2B3,kmul(J32L,PDdissipationNth3B3))); + + JacPDdissipationNth2beta1 = + kmadd(J12L,PDdissipationNth1beta1,kmadd(J22L,PDdissipationNth2beta1,kmul(J32L,PDdissipationNth3beta1))); + + JacPDdissipationNth2beta2 = + kmadd(J12L,PDdissipationNth1beta2,kmadd(J22L,PDdissipationNth2beta2,kmul(J32L,PDdissipationNth3beta2))); + + JacPDdissipationNth2beta3 = + kmadd(J12L,PDdissipationNth1beta3,kmadd(J22L,PDdissipationNth2beta3,kmul(J32L,PDdissipationNth3beta3))); + + JacPDdissipationNth2gt11 = + kmadd(J12L,PDdissipationNth1gt11,kmadd(J22L,PDdissipationNth2gt11,kmul(J32L,PDdissipationNth3gt11))); + + JacPDdissipationNth2gt12 = + kmadd(J12L,PDdissipationNth1gt12,kmadd(J22L,PDdissipationNth2gt12,kmul(J32L,PDdissipationNth3gt12))); + + JacPDdissipationNth2gt13 = + kmadd(J12L,PDdissipationNth1gt13,kmadd(J22L,PDdissipationNth2gt13,kmul(J32L,PDdissipationNth3gt13))); + + JacPDdissipationNth2gt22 = + kmadd(J12L,PDdissipationNth1gt22,kmadd(J22L,PDdissipationNth2gt22,kmul(J32L,PDdissipationNth3gt22))); + + JacPDdissipationNth2gt23 = + kmadd(J12L,PDdissipationNth1gt23,kmadd(J22L,PDdissipationNth2gt23,kmul(J32L,PDdissipationNth3gt23))); + + JacPDdissipationNth2gt33 = + kmadd(J12L,PDdissipationNth1gt33,kmadd(J22L,PDdissipationNth2gt33,kmul(J32L,PDdissipationNth3gt33))); + + JacPDdissipationNth2phi = + kmadd(J12L,PDdissipationNth1phi,kmadd(J22L,PDdissipationNth2phi,kmul(J32L,PDdissipationNth3phi))); + + JacPDdissipationNth2trK = + kmadd(J12L,PDdissipationNth1trK,kmadd(J22L,PDdissipationNth2trK,kmul(J32L,PDdissipationNth3trK))); + + JacPDdissipationNth2Xt1 = + kmadd(J12L,PDdissipationNth1Xt1,kmadd(J22L,PDdissipationNth2Xt1,kmul(J32L,PDdissipationNth3Xt1))); + + JacPDdissipationNth2Xt2 = + kmadd(J12L,PDdissipationNth1Xt2,kmadd(J22L,PDdissipationNth2Xt2,kmul(J32L,PDdissipationNth3Xt2))); + + JacPDdissipationNth2Xt3 = + kmadd(J12L,PDdissipationNth1Xt3,kmadd(J22L,PDdissipationNth2Xt3,kmul(J32L,PDdissipationNth3Xt3))); + + JacPDdissipationNth3A = + kmadd(J13L,PDdissipationNth1A,kmadd(J23L,PDdissipationNth2A,kmul(J33L,PDdissipationNth3A))); + + JacPDdissipationNth3alpha = + kmadd(J13L,PDdissipationNth1alpha,kmadd(J23L,PDdissipationNth2alpha,kmul(J33L,PDdissipationNth3alpha))); + + JacPDdissipationNth3At11 = + kmadd(J13L,PDdissipationNth1At11,kmadd(J23L,PDdissipationNth2At11,kmul(J33L,PDdissipationNth3At11))); + + JacPDdissipationNth3At12 = + kmadd(J13L,PDdissipationNth1At12,kmadd(J23L,PDdissipationNth2At12,kmul(J33L,PDdissipationNth3At12))); + + JacPDdissipationNth3At13 = + kmadd(J13L,PDdissipationNth1At13,kmadd(J23L,PDdissipationNth2At13,kmul(J33L,PDdissipationNth3At13))); + + JacPDdissipationNth3At22 = + kmadd(J13L,PDdissipationNth1At22,kmadd(J23L,PDdissipationNth2At22,kmul(J33L,PDdissipationNth3At22))); + + JacPDdissipationNth3At23 = + kmadd(J13L,PDdissipationNth1At23,kmadd(J23L,PDdissipationNth2At23,kmul(J33L,PDdissipationNth3At23))); + + JacPDdissipationNth3At33 = + kmadd(J13L,PDdissipationNth1At33,kmadd(J23L,PDdissipationNth2At33,kmul(J33L,PDdissipationNth3At33))); + + JacPDdissipationNth3B1 = + kmadd(J13L,PDdissipationNth1B1,kmadd(J23L,PDdissipationNth2B1,kmul(J33L,PDdissipationNth3B1))); + + JacPDdissipationNth3B2 = + kmadd(J13L,PDdissipationNth1B2,kmadd(J23L,PDdissipationNth2B2,kmul(J33L,PDdissipationNth3B2))); + + JacPDdissipationNth3B3 = + kmadd(J13L,PDdissipationNth1B3,kmadd(J23L,PDdissipationNth2B3,kmul(J33L,PDdissipationNth3B3))); + + JacPDdissipationNth3beta1 = + kmadd(J13L,PDdissipationNth1beta1,kmadd(J23L,PDdissipationNth2beta1,kmul(J33L,PDdissipationNth3beta1))); + + JacPDdissipationNth3beta2 = + kmadd(J13L,PDdissipationNth1beta2,kmadd(J23L,PDdissipationNth2beta2,kmul(J33L,PDdissipationNth3beta2))); + + JacPDdissipationNth3beta3 = + kmadd(J13L,PDdissipationNth1beta3,kmadd(J23L,PDdissipationNth2beta3,kmul(J33L,PDdissipationNth3beta3))); + + JacPDdissipationNth3gt11 = + kmadd(J13L,PDdissipationNth1gt11,kmadd(J23L,PDdissipationNth2gt11,kmul(J33L,PDdissipationNth3gt11))); + + JacPDdissipationNth3gt12 = + kmadd(J13L,PDdissipationNth1gt12,kmadd(J23L,PDdissipationNth2gt12,kmul(J33L,PDdissipationNth3gt12))); + + JacPDdissipationNth3gt13 = + kmadd(J13L,PDdissipationNth1gt13,kmadd(J23L,PDdissipationNth2gt13,kmul(J33L,PDdissipationNth3gt13))); + + JacPDdissipationNth3gt22 = + kmadd(J13L,PDdissipationNth1gt22,kmadd(J23L,PDdissipationNth2gt22,kmul(J33L,PDdissipationNth3gt22))); + + JacPDdissipationNth3gt23 = + kmadd(J13L,PDdissipationNth1gt23,kmadd(J23L,PDdissipationNth2gt23,kmul(J33L,PDdissipationNth3gt23))); + + JacPDdissipationNth3gt33 = + kmadd(J13L,PDdissipationNth1gt33,kmadd(J23L,PDdissipationNth2gt33,kmul(J33L,PDdissipationNth3gt33))); + + JacPDdissipationNth3phi = + kmadd(J13L,PDdissipationNth1phi,kmadd(J23L,PDdissipationNth2phi,kmul(J33L,PDdissipationNth3phi))); + + JacPDdissipationNth3trK = + kmadd(J13L,PDdissipationNth1trK,kmadd(J23L,PDdissipationNth2trK,kmul(J33L,PDdissipationNth3trK))); + + JacPDdissipationNth3Xt1 = + kmadd(J13L,PDdissipationNth1Xt1,kmadd(J23L,PDdissipationNth2Xt1,kmul(J33L,PDdissipationNth3Xt1))); + + JacPDdissipationNth3Xt2 = + kmadd(J13L,PDdissipationNth1Xt2,kmadd(J23L,PDdissipationNth2Xt2,kmul(J33L,PDdissipationNth3Xt2))); + + JacPDdissipationNth3Xt3 = + kmadd(J13L,PDdissipationNth1Xt3,kmadd(J23L,PDdissipationNth2Xt3,kmul(J33L,PDdissipationNth3Xt3))); + } + else + { + JacPDdissipationNth1A = PDdissipationNth1A; + + JacPDdissipationNth1alpha = PDdissipationNth1alpha; + + JacPDdissipationNth1At11 = PDdissipationNth1At11; + + JacPDdissipationNth1At12 = PDdissipationNth1At12; + + JacPDdissipationNth1At13 = PDdissipationNth1At13; + + JacPDdissipationNth1At22 = PDdissipationNth1At22; + + JacPDdissipationNth1At23 = PDdissipationNth1At23; + + JacPDdissipationNth1At33 = PDdissipationNth1At33; + + JacPDdissipationNth1B1 = PDdissipationNth1B1; + + JacPDdissipationNth1B2 = PDdissipationNth1B2; + + JacPDdissipationNth1B3 = PDdissipationNth1B3; + + JacPDdissipationNth1beta1 = PDdissipationNth1beta1; + + JacPDdissipationNth1beta2 = PDdissipationNth1beta2; + + JacPDdissipationNth1beta3 = PDdissipationNth1beta3; + + JacPDdissipationNth1gt11 = PDdissipationNth1gt11; + + JacPDdissipationNth1gt12 = PDdissipationNth1gt12; + + JacPDdissipationNth1gt13 = PDdissipationNth1gt13; + + JacPDdissipationNth1gt22 = PDdissipationNth1gt22; + + JacPDdissipationNth1gt23 = PDdissipationNth1gt23; + + JacPDdissipationNth1gt33 = PDdissipationNth1gt33; + + JacPDdissipationNth1phi = PDdissipationNth1phi; + + JacPDdissipationNth1trK = PDdissipationNth1trK; + + JacPDdissipationNth1Xt1 = PDdissipationNth1Xt1; + + JacPDdissipationNth1Xt2 = PDdissipationNth1Xt2; + + JacPDdissipationNth1Xt3 = PDdissipationNth1Xt3; + + JacPDdissipationNth2A = PDdissipationNth2A; + + JacPDdissipationNth2alpha = PDdissipationNth2alpha; + + JacPDdissipationNth2At11 = PDdissipationNth2At11; + + JacPDdissipationNth2At12 = PDdissipationNth2At12; + + JacPDdissipationNth2At13 = PDdissipationNth2At13; + + JacPDdissipationNth2At22 = PDdissipationNth2At22; + + JacPDdissipationNth2At23 = PDdissipationNth2At23; + + JacPDdissipationNth2At33 = PDdissipationNth2At33; + + JacPDdissipationNth2B1 = PDdissipationNth2B1; + + JacPDdissipationNth2B2 = PDdissipationNth2B2; + + JacPDdissipationNth2B3 = PDdissipationNth2B3; + + JacPDdissipationNth2beta1 = PDdissipationNth2beta1; + + JacPDdissipationNth2beta2 = PDdissipationNth2beta2; + + JacPDdissipationNth2beta3 = PDdissipationNth2beta3; + + JacPDdissipationNth2gt11 = PDdissipationNth2gt11; + + JacPDdissipationNth2gt12 = PDdissipationNth2gt12; + + JacPDdissipationNth2gt13 = PDdissipationNth2gt13; + + JacPDdissipationNth2gt22 = PDdissipationNth2gt22; + + JacPDdissipationNth2gt23 = PDdissipationNth2gt23; + + JacPDdissipationNth2gt33 = PDdissipationNth2gt33; + + JacPDdissipationNth2phi = PDdissipationNth2phi; + + JacPDdissipationNth2trK = PDdissipationNth2trK; + + JacPDdissipationNth2Xt1 = PDdissipationNth2Xt1; + + JacPDdissipationNth2Xt2 = PDdissipationNth2Xt2; + + JacPDdissipationNth2Xt3 = PDdissipationNth2Xt3; + + JacPDdissipationNth3A = PDdissipationNth3A; + + JacPDdissipationNth3alpha = PDdissipationNth3alpha; + + JacPDdissipationNth3At11 = PDdissipationNth3At11; + + JacPDdissipationNth3At12 = PDdissipationNth3At12; + + JacPDdissipationNth3At13 = PDdissipationNth3At13; + + JacPDdissipationNth3At22 = PDdissipationNth3At22; + + JacPDdissipationNth3At23 = PDdissipationNth3At23; + + JacPDdissipationNth3At33 = PDdissipationNth3At33; + + JacPDdissipationNth3B1 = PDdissipationNth3B1; + + JacPDdissipationNth3B2 = PDdissipationNth3B2; + + JacPDdissipationNth3B3 = PDdissipationNth3B3; + + JacPDdissipationNth3beta1 = PDdissipationNth3beta1; + + JacPDdissipationNth3beta2 = PDdissipationNth3beta2; + + JacPDdissipationNth3beta3 = PDdissipationNth3beta3; + + JacPDdissipationNth3gt11 = PDdissipationNth3gt11; + + JacPDdissipationNth3gt12 = PDdissipationNth3gt12; + + JacPDdissipationNth3gt13 = PDdissipationNth3gt13; + + JacPDdissipationNth3gt22 = PDdissipationNth3gt22; + + JacPDdissipationNth3gt23 = PDdissipationNth3gt23; + + JacPDdissipationNth3gt33 = PDdissipationNth3gt33; + + JacPDdissipationNth3phi = PDdissipationNth3phi; + + JacPDdissipationNth3trK = PDdissipationNth3trK; + + JacPDdissipationNth3Xt1 = PDdissipationNth3Xt1; + + JacPDdissipationNth3Xt2 = PDdissipationNth3Xt2; + + JacPDdissipationNth3Xt3 = PDdissipationNth3Xt3; + } + CCTK_REAL_VEC epsdiss1 = ToReal(EpsDiss); CCTK_REAL_VEC epsdiss2 = ToReal(EpsDiss); @@ -276,79 +1184,79 @@ static void ML_BSSN_Dissipation_Body(cGH const * restrict const cctkGH, int cons CCTK_REAL_VEC epsdiss3 = ToReal(EpsDiss); phirhsL = - kmadd(epsdiss1,PDdissipationNth1phi,kmadd(epsdiss2,PDdissipationNth2phi,kmadd(epsdiss3,PDdissipationNth3phi,phirhsL))); + kmadd(epsdiss1,JacPDdissipationNth1phi,kmadd(epsdiss2,JacPDdissipationNth2phi,kmadd(epsdiss3,JacPDdissipationNth3phi,phirhsL))); gt11rhsL = - kadd(gt11rhsL,kmadd(epsdiss1,PDdissipationNth1gt11,kmadd(epsdiss2,PDdissipationNth2gt11,kmul(epsdiss3,PDdissipationNth3gt11)))); + kadd(gt11rhsL,kmadd(epsdiss1,JacPDdissipationNth1gt11,kmadd(epsdiss2,JacPDdissipationNth2gt11,kmul(epsdiss3,JacPDdissipationNth3gt11)))); gt12rhsL = - kadd(gt12rhsL,kmadd(epsdiss1,PDdissipationNth1gt12,kmadd(epsdiss2,PDdissipationNth2gt12,kmul(epsdiss3,PDdissipationNth3gt12)))); + kadd(gt12rhsL,kmadd(epsdiss1,JacPDdissipationNth1gt12,kmadd(epsdiss2,JacPDdissipationNth2gt12,kmul(epsdiss3,JacPDdissipationNth3gt12)))); gt13rhsL = - kadd(gt13rhsL,kmadd(epsdiss1,PDdissipationNth1gt13,kmadd(epsdiss2,PDdissipationNth2gt13,kmul(epsdiss3,PDdissipationNth3gt13)))); + kadd(gt13rhsL,kmadd(epsdiss1,JacPDdissipationNth1gt13,kmadd(epsdiss2,JacPDdissipationNth2gt13,kmul(epsdiss3,JacPDdissipationNth3gt13)))); gt22rhsL = - kadd(gt22rhsL,kmadd(epsdiss1,PDdissipationNth1gt22,kmadd(epsdiss2,PDdissipationNth2gt22,kmul(epsdiss3,PDdissipationNth3gt22)))); + kadd(gt22rhsL,kmadd(epsdiss1,JacPDdissipationNth1gt22,kmadd(epsdiss2,JacPDdissipationNth2gt22,kmul(epsdiss3,JacPDdissipationNth3gt22)))); gt23rhsL = - kadd(gt23rhsL,kmadd(epsdiss1,PDdissipationNth1gt23,kmadd(epsdiss2,PDdissipationNth2gt23,kmul(epsdiss3,PDdissipationNth3gt23)))); + kadd(gt23rhsL,kmadd(epsdiss1,JacPDdissipationNth1gt23,kmadd(epsdiss2,JacPDdissipationNth2gt23,kmul(epsdiss3,JacPDdissipationNth3gt23)))); gt33rhsL = - kadd(gt33rhsL,kmadd(epsdiss1,PDdissipationNth1gt33,kmadd(epsdiss2,PDdissipationNth2gt33,kmul(epsdiss3,PDdissipationNth3gt33)))); + kadd(gt33rhsL,kmadd(epsdiss1,JacPDdissipationNth1gt33,kmadd(epsdiss2,JacPDdissipationNth2gt33,kmul(epsdiss3,JacPDdissipationNth3gt33)))); Xt1rhsL = - kmadd(epsdiss1,PDdissipationNth1Xt1,kmadd(epsdiss2,PDdissipationNth2Xt1,kmadd(epsdiss3,PDdissipationNth3Xt1,Xt1rhsL))); + kmadd(epsdiss1,JacPDdissipationNth1Xt1,kmadd(epsdiss2,JacPDdissipationNth2Xt1,kmadd(epsdiss3,JacPDdissipationNth3Xt1,Xt1rhsL))); Xt2rhsL = - kmadd(epsdiss1,PDdissipationNth1Xt2,kmadd(epsdiss2,PDdissipationNth2Xt2,kmadd(epsdiss3,PDdissipationNth3Xt2,Xt2rhsL))); + kmadd(epsdiss1,JacPDdissipationNth1Xt2,kmadd(epsdiss2,JacPDdissipationNth2Xt2,kmadd(epsdiss3,JacPDdissipationNth3Xt2,Xt2rhsL))); Xt3rhsL = - kmadd(epsdiss1,PDdissipationNth1Xt3,kmadd(epsdiss2,PDdissipationNth2Xt3,kmadd(epsdiss3,PDdissipationNth3Xt3,Xt3rhsL))); + kmadd(epsdiss1,JacPDdissipationNth1Xt3,kmadd(epsdiss2,JacPDdissipationNth2Xt3,kmadd(epsdiss3,JacPDdissipationNth3Xt3,Xt3rhsL))); trKrhsL = - kmadd(epsdiss1,PDdissipationNth1trK,kmadd(epsdiss2,PDdissipationNth2trK,kmadd(epsdiss3,PDdissipationNth3trK,trKrhsL))); + kmadd(epsdiss1,JacPDdissipationNth1trK,kmadd(epsdiss2,JacPDdissipationNth2trK,kmadd(epsdiss3,JacPDdissipationNth3trK,trKrhsL))); At11rhsL = - kadd(At11rhsL,kmadd(epsdiss1,PDdissipationNth1At11,kmadd(epsdiss2,PDdissipationNth2At11,kmul(epsdiss3,PDdissipationNth3At11)))); + kadd(At11rhsL,kmadd(epsdiss1,JacPDdissipationNth1At11,kmadd(epsdiss2,JacPDdissipationNth2At11,kmul(epsdiss3,JacPDdissipationNth3At11)))); At12rhsL = - kadd(At12rhsL,kmadd(epsdiss1,PDdissipationNth1At12,kmadd(epsdiss2,PDdissipationNth2At12,kmul(epsdiss3,PDdissipationNth3At12)))); + kadd(At12rhsL,kmadd(epsdiss1,JacPDdissipationNth1At12,kmadd(epsdiss2,JacPDdissipationNth2At12,kmul(epsdiss3,JacPDdissipationNth3At12)))); At13rhsL = - kadd(At13rhsL,kmadd(epsdiss1,PDdissipationNth1At13,kmadd(epsdiss2,PDdissipationNth2At13,kmul(epsdiss3,PDdissipationNth3At13)))); + kadd(At13rhsL,kmadd(epsdiss1,JacPDdissipationNth1At13,kmadd(epsdiss2,JacPDdissipationNth2At13,kmul(epsdiss3,JacPDdissipationNth3At13)))); At22rhsL = - kadd(At22rhsL,kmadd(epsdiss1,PDdissipationNth1At22,kmadd(epsdiss2,PDdissipationNth2At22,kmul(epsdiss3,PDdissipationNth3At22)))); + kadd(At22rhsL,kmadd(epsdiss1,JacPDdissipationNth1At22,kmadd(epsdiss2,JacPDdissipationNth2At22,kmul(epsdiss3,JacPDdissipationNth3At22)))); At23rhsL = - kadd(At23rhsL,kmadd(epsdiss1,PDdissipationNth1At23,kmadd(epsdiss2,PDdissipationNth2At23,kmul(epsdiss3,PDdissipationNth3At23)))); + kadd(At23rhsL,kmadd(epsdiss1,JacPDdissipationNth1At23,kmadd(epsdiss2,JacPDdissipationNth2At23,kmul(epsdiss3,JacPDdissipationNth3At23)))); At33rhsL = - kadd(At33rhsL,kmadd(epsdiss1,PDdissipationNth1At33,kmadd(epsdiss2,PDdissipationNth2At33,kmul(epsdiss3,PDdissipationNth3At33)))); + kadd(At33rhsL,kmadd(epsdiss1,JacPDdissipationNth1At33,kmadd(epsdiss2,JacPDdissipationNth2At33,kmul(epsdiss3,JacPDdissipationNth3At33)))); alpharhsL = - kadd(alpharhsL,kmadd(epsdiss1,PDdissipationNth1alpha,kmadd(epsdiss2,PDdissipationNth2alpha,kmul(epsdiss3,PDdissipationNth3alpha)))); + kadd(alpharhsL,kmadd(epsdiss1,JacPDdissipationNth1alpha,kmadd(epsdiss2,JacPDdissipationNth2alpha,kmul(epsdiss3,JacPDdissipationNth3alpha)))); ArhsL = - kadd(ArhsL,kmadd(epsdiss1,PDdissipationNth1A,kmadd(epsdiss2,PDdissipationNth2A,kmul(epsdiss3,PDdissipationNth3A)))); + kadd(ArhsL,kmadd(epsdiss1,JacPDdissipationNth1A,kmadd(epsdiss2,JacPDdissipationNth2A,kmul(epsdiss3,JacPDdissipationNth3A)))); beta1rhsL = - kadd(beta1rhsL,kmadd(epsdiss1,PDdissipationNth1beta1,kmadd(epsdiss2,PDdissipationNth2beta1,kmul(epsdiss3,PDdissipationNth3beta1)))); + kadd(beta1rhsL,kmadd(epsdiss1,JacPDdissipationNth1beta1,kmadd(epsdiss2,JacPDdissipationNth2beta1,kmul(epsdiss3,JacPDdissipationNth3beta1)))); beta2rhsL = - kadd(beta2rhsL,kmadd(epsdiss1,PDdissipationNth1beta2,kmadd(epsdiss2,PDdissipationNth2beta2,kmul(epsdiss3,PDdissipationNth3beta2)))); + kadd(beta2rhsL,kmadd(epsdiss1,JacPDdissipationNth1beta2,kmadd(epsdiss2,JacPDdissipationNth2beta2,kmul(epsdiss3,JacPDdissipationNth3beta2)))); beta3rhsL = - kadd(beta3rhsL,kmadd(epsdiss1,PDdissipationNth1beta3,kmadd(epsdiss2,PDdissipationNth2beta3,kmul(epsdiss3,PDdissipationNth3beta3)))); + kadd(beta3rhsL,kmadd(epsdiss1,JacPDdissipationNth1beta3,kmadd(epsdiss2,JacPDdissipationNth2beta3,kmul(epsdiss3,JacPDdissipationNth3beta3)))); B1rhsL = - kadd(B1rhsL,kmadd(epsdiss1,PDdissipationNth1B1,kmadd(epsdiss2,PDdissipationNth2B1,kmul(epsdiss3,PDdissipationNth3B1)))); + kadd(B1rhsL,kmadd(epsdiss1,JacPDdissipationNth1B1,kmadd(epsdiss2,JacPDdissipationNth2B1,kmul(epsdiss3,JacPDdissipationNth3B1)))); B2rhsL = - kadd(B2rhsL,kmadd(epsdiss1,PDdissipationNth1B2,kmadd(epsdiss2,PDdissipationNth2B2,kmul(epsdiss3,PDdissipationNth3B2)))); + kadd(B2rhsL,kmadd(epsdiss1,JacPDdissipationNth1B2,kmadd(epsdiss2,JacPDdissipationNth2B2,kmul(epsdiss3,JacPDdissipationNth3B2)))); B3rhsL = - kadd(B3rhsL,kmadd(epsdiss1,PDdissipationNth1B3,kmadd(epsdiss2,PDdissipationNth2B3,kmul(epsdiss3,PDdissipationNth3B3)))); + kadd(B3rhsL,kmadd(epsdiss1,JacPDdissipationNth1B3,kmadd(epsdiss2,JacPDdissipationNth2B3,kmul(epsdiss3,JacPDdissipationNth3B3)))); /* If necessary, store only partial vectors after the first iteration */ |