diff options
Diffstat (limited to 'ML_BSSN_MP_O8/src/ML_BSSN_MP_O8_convertToADMBaseDtLapseShift.cc')
| -rw-r--r-- | ML_BSSN_MP_O8/src/ML_BSSN_MP_O8_convertToADMBaseDtLapseShift.cc | 918 |
1 files changed, 0 insertions, 918 deletions
diff --git a/ML_BSSN_MP_O8/src/ML_BSSN_MP_O8_convertToADMBaseDtLapseShift.cc b/ML_BSSN_MP_O8/src/ML_BSSN_MP_O8_convertToADMBaseDtLapseShift.cc deleted file mode 100644 index ceecbd6..0000000 --- a/ML_BSSN_MP_O8/src/ML_BSSN_MP_O8_convertToADMBaseDtLapseShift.cc +++ /dev/null @@ -1,918 +0,0 @@ -/* File produced by Kranc */ - -#define KRANC_C - -#include <assert.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" -#include "GenericFD.h" -#include "Differencing.h" -#include "cctk_Loop.h" -#include "loopcontrol.h" -#include "vectors.h" - -/* Define macros used in calculations */ -#define INITVALUE (42) -#define QAD(x) (SQR(SQR(x))) -#define INV(x) (kdiv(ToReal(1.0),x)) -#define SQR(x) (kmul(x,x)) -#define CUB(x) (kmul(x,SQR(x))) - -extern "C" void ML_BSSN_MP_O8_convertToADMBaseDtLapseShift_SelectBCs(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - CCTK_INT ierr = 0; - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ADMBase::dtlapse","flat"); - if (ierr < 0) - CCTK_WARN(1, "Failed to register flat BC for ADMBase::dtlapse."); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ADMBase::dtshift","flat"); - if (ierr < 0) - CCTK_WARN(1, "Failed to register flat BC for ADMBase::dtshift."); - return; -} - -static void ML_BSSN_MP_O8_convertToADMBaseDtLapseShift_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - - /* Include user-supplied include files */ - - /* Initialise finite differencing variables */ - ptrdiff_t const di = 1; - ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); - ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); - ptrdiff_t const cdi = sizeof(CCTK_REAL) * di; - ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj; - ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk; - CCTK_REAL_VEC const dx = ToReal(CCTK_DELTA_SPACE(0)); - CCTK_REAL_VEC const dy = ToReal(CCTK_DELTA_SPACE(1)); - CCTK_REAL_VEC const dz = ToReal(CCTK_DELTA_SPACE(2)); - CCTK_REAL_VEC const dt = ToReal(CCTK_DELTA_TIME); - CCTK_REAL_VEC const t = ToReal(cctk_time); - CCTK_REAL_VEC const dxi = INV(dx); - CCTK_REAL_VEC const dyi = INV(dy); - CCTK_REAL_VEC const dzi = INV(dz); - CCTK_REAL_VEC const khalf = ToReal(0.5); - CCTK_REAL_VEC const kthird = ToReal(1.0/3.0); - CCTK_REAL_VEC const ktwothird = ToReal(2.0/3.0); - CCTK_REAL_VEC const kfourthird = ToReal(4.0/3.0); - CCTK_REAL_VEC const keightthird = ToReal(8.0/3.0); - CCTK_REAL_VEC const hdxi = kmul(ToReal(0.5), dxi); - CCTK_REAL_VEC const hdyi = kmul(ToReal(0.5), dyi); - 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(kmul(dx,dy)),ToReal(0.00694444444444444444444444444444)); - CCTK_REAL_VEC const p1o144dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.00694444444444444444444444444444)); - CCTK_REAL_VEC const p1o144dydz = kmul(INV(kmul(dy,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(kmul(dx,dy)),ToReal(0.000277777777777777777777777777778)); - CCTK_REAL_VEC const p1o3600dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.000277777777777777777777777777778)); - CCTK_REAL_VEC const p1o3600dydz = kmul(INV(kmul(dy,dz)),ToReal(0.000277777777777777777777777777778)); - CCTK_REAL_VEC const p1o4dx = kmul(INV(dx),ToReal(0.25)); - CCTK_REAL_VEC const p1o4dxdy = kmul(INV(kmul(dx,dy)),ToReal(0.25)); - CCTK_REAL_VEC const p1o4dxdz = kmul(INV(kmul(dx,dz)),ToReal(0.25)); - CCTK_REAL_VEC const p1o4dy = kmul(INV(dy),ToReal(0.25)); - CCTK_REAL_VEC const p1o4dydz = kmul(INV(kmul(dy,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(kmul(dx,dy)),ToReal(1.41723356009070294784580498866e-6)); - CCTK_REAL_VEC const p1o705600dxdz = kmul(INV(kmul(dx,dz)),ToReal(1.41723356009070294784580498866e-6)); - CCTK_REAL_VEC const p1o705600dydz = kmul(INV(kmul(dy,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; - - /* Assign local copies of arrays functions */ - - - - /* Calculate temporaries and arrays functions */ - - /* Copy local copies back to grid functions */ - - /* Loop over the grid points */ - #pragma omp parallel - LC_LOOP3VEC(ML_BSSN_MP_O8_convertToADMBaseDtLapseShift, - i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2], - cctk_lsh[0],cctk_lsh[1],cctk_lsh[2], - CCTK_REAL_VEC_SIZE) - { - ptrdiff_t const index = di*i + dj*j + dk*k; - - /* Assign local copies of grid functions */ - - CCTK_REAL_VEC AL = vec_load(A[index]); - CCTK_REAL_VEC alphaL = vec_load(alpha[index]); - CCTK_REAL_VEC B1L = vec_load(B1[index]); - CCTK_REAL_VEC B2L = vec_load(B2[index]); - CCTK_REAL_VEC B3L = vec_load(B3[index]); - 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]); - CCTK_REAL_VEC Xt2L = vec_load(Xt2[index]); - CCTK_REAL_VEC Xt3L = vec_load(Xt3[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 PDstandardNth1alpha; - CCTK_REAL_VEC PDstandardNth2alpha; - CCTK_REAL_VEC PDstandardNth3alpha; - CCTK_REAL_VEC PDupwindNthAnti1alpha; - CCTK_REAL_VEC PDupwindNthSymm1alpha; - CCTK_REAL_VEC PDupwindNthAnti2alpha; - CCTK_REAL_VEC PDupwindNthSymm2alpha; - CCTK_REAL_VEC PDupwindNthAnti3alpha; - CCTK_REAL_VEC PDupwindNthSymm3alpha; - CCTK_REAL_VEC PDupwindNthAnti1beta1; - CCTK_REAL_VEC PDupwindNthSymm1beta1; - CCTK_REAL_VEC PDupwindNthAnti2beta1; - CCTK_REAL_VEC PDupwindNthSymm2beta1; - CCTK_REAL_VEC PDupwindNthAnti3beta1; - CCTK_REAL_VEC PDupwindNthSymm3beta1; - CCTK_REAL_VEC PDupwindNthAnti1beta2; - CCTK_REAL_VEC PDupwindNthSymm1beta2; - CCTK_REAL_VEC PDupwindNthAnti2beta2; - CCTK_REAL_VEC PDupwindNthSymm2beta2; - CCTK_REAL_VEC PDupwindNthAnti3beta2; - CCTK_REAL_VEC PDupwindNthSymm3beta2; - CCTK_REAL_VEC PDupwindNthAnti1beta3; - CCTK_REAL_VEC PDupwindNthSymm1beta3; - CCTK_REAL_VEC PDupwindNthAnti2beta3; - 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]); - PDupwindNthSymm2alpha = PDupwindNthSymmfdOrder22(&alpha[index]); - PDupwindNthAnti3alpha = PDupwindNthAntifdOrder23(&alpha[index]); - PDupwindNthSymm3alpha = PDupwindNthSymmfdOrder23(&alpha[index]); - PDupwindNthAnti1beta1 = PDupwindNthAntifdOrder21(&beta1[index]); - PDupwindNthSymm1beta1 = PDupwindNthSymmfdOrder21(&beta1[index]); - PDupwindNthAnti2beta1 = PDupwindNthAntifdOrder22(&beta1[index]); - PDupwindNthSymm2beta1 = PDupwindNthSymmfdOrder22(&beta1[index]); - PDupwindNthAnti3beta1 = PDupwindNthAntifdOrder23(&beta1[index]); - PDupwindNthSymm3beta1 = PDupwindNthSymmfdOrder23(&beta1[index]); - PDupwindNthAnti1beta2 = PDupwindNthAntifdOrder21(&beta2[index]); - PDupwindNthSymm1beta2 = PDupwindNthSymmfdOrder21(&beta2[index]); - PDupwindNthAnti2beta2 = PDupwindNthAntifdOrder22(&beta2[index]); - PDupwindNthSymm2beta2 = PDupwindNthSymmfdOrder22(&beta2[index]); - PDupwindNthAnti3beta2 = PDupwindNthAntifdOrder23(&beta2[index]); - PDupwindNthSymm3beta2 = PDupwindNthSymmfdOrder23(&beta2[index]); - PDupwindNthAnti1beta3 = PDupwindNthAntifdOrder21(&beta3[index]); - PDupwindNthSymm1beta3 = PDupwindNthSymmfdOrder21(&beta3[index]); - PDupwindNthAnti2beta3 = PDupwindNthAntifdOrder22(&beta3[index]); - 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]); - PDupwindNthSymm2alpha = PDupwindNthSymmfdOrder42(&alpha[index]); - PDupwindNthAnti3alpha = PDupwindNthAntifdOrder43(&alpha[index]); - PDupwindNthSymm3alpha = PDupwindNthSymmfdOrder43(&alpha[index]); - PDupwindNthAnti1beta1 = PDupwindNthAntifdOrder41(&beta1[index]); - PDupwindNthSymm1beta1 = PDupwindNthSymmfdOrder41(&beta1[index]); - PDupwindNthAnti2beta1 = PDupwindNthAntifdOrder42(&beta1[index]); - PDupwindNthSymm2beta1 = PDupwindNthSymmfdOrder42(&beta1[index]); - PDupwindNthAnti3beta1 = PDupwindNthAntifdOrder43(&beta1[index]); - PDupwindNthSymm3beta1 = PDupwindNthSymmfdOrder43(&beta1[index]); - PDupwindNthAnti1beta2 = PDupwindNthAntifdOrder41(&beta2[index]); - PDupwindNthSymm1beta2 = PDupwindNthSymmfdOrder41(&beta2[index]); - PDupwindNthAnti2beta2 = PDupwindNthAntifdOrder42(&beta2[index]); - PDupwindNthSymm2beta2 = PDupwindNthSymmfdOrder42(&beta2[index]); - PDupwindNthAnti3beta2 = PDupwindNthAntifdOrder43(&beta2[index]); - PDupwindNthSymm3beta2 = PDupwindNthSymmfdOrder43(&beta2[index]); - PDupwindNthAnti1beta3 = PDupwindNthAntifdOrder41(&beta3[index]); - PDupwindNthSymm1beta3 = PDupwindNthSymmfdOrder41(&beta3[index]); - PDupwindNthAnti2beta3 = PDupwindNthAntifdOrder42(&beta3[index]); - 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]); - PDupwindNthSymm2alpha = PDupwindNthSymmfdOrder62(&alpha[index]); - PDupwindNthAnti3alpha = PDupwindNthAntifdOrder63(&alpha[index]); - PDupwindNthSymm3alpha = PDupwindNthSymmfdOrder63(&alpha[index]); - PDupwindNthAnti1beta1 = PDupwindNthAntifdOrder61(&beta1[index]); - PDupwindNthSymm1beta1 = PDupwindNthSymmfdOrder61(&beta1[index]); - PDupwindNthAnti2beta1 = PDupwindNthAntifdOrder62(&beta1[index]); - PDupwindNthSymm2beta1 = PDupwindNthSymmfdOrder62(&beta1[index]); - PDupwindNthAnti3beta1 = PDupwindNthAntifdOrder63(&beta1[index]); - PDupwindNthSymm3beta1 = PDupwindNthSymmfdOrder63(&beta1[index]); - PDupwindNthAnti1beta2 = PDupwindNthAntifdOrder61(&beta2[index]); - PDupwindNthSymm1beta2 = PDupwindNthSymmfdOrder61(&beta2[index]); - PDupwindNthAnti2beta2 = PDupwindNthAntifdOrder62(&beta2[index]); - PDupwindNthSymm2beta2 = PDupwindNthSymmfdOrder62(&beta2[index]); - PDupwindNthAnti3beta2 = PDupwindNthAntifdOrder63(&beta2[index]); - PDupwindNthSymm3beta2 = PDupwindNthSymmfdOrder63(&beta2[index]); - PDupwindNthAnti1beta3 = PDupwindNthAntifdOrder61(&beta3[index]); - PDupwindNthSymm1beta3 = PDupwindNthSymmfdOrder61(&beta3[index]); - PDupwindNthAnti2beta3 = PDupwindNthAntifdOrder62(&beta3[index]); - 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]); - PDupwindNthSymm2alpha = PDupwindNthSymmfdOrder82(&alpha[index]); - PDupwindNthAnti3alpha = PDupwindNthAntifdOrder83(&alpha[index]); - PDupwindNthSymm3alpha = PDupwindNthSymmfdOrder83(&alpha[index]); - PDupwindNthAnti1beta1 = PDupwindNthAntifdOrder81(&beta1[index]); - PDupwindNthSymm1beta1 = PDupwindNthSymmfdOrder81(&beta1[index]); - PDupwindNthAnti2beta1 = PDupwindNthAntifdOrder82(&beta1[index]); - PDupwindNthSymm2beta1 = PDupwindNthSymmfdOrder82(&beta1[index]); - PDupwindNthAnti3beta1 = PDupwindNthAntifdOrder83(&beta1[index]); - PDupwindNthSymm3beta1 = PDupwindNthSymmfdOrder83(&beta1[index]); - PDupwindNthAnti1beta2 = PDupwindNthAntifdOrder81(&beta2[index]); - PDupwindNthSymm1beta2 = PDupwindNthSymmfdOrder81(&beta2[index]); - PDupwindNthAnti2beta2 = PDupwindNthAntifdOrder82(&beta2[index]); - PDupwindNthSymm2beta2 = PDupwindNthSymmfdOrder82(&beta2[index]); - PDupwindNthAnti3beta2 = PDupwindNthAntifdOrder83(&beta2[index]); - PDupwindNthSymm3beta2 = PDupwindNthSymmfdOrder83(&beta2[index]); - PDupwindNthAnti1beta3 = PDupwindNthAntifdOrder81(&beta3[index]); - PDupwindNthSymm1beta3 = PDupwindNthSymmfdOrder81(&beta3[index]); - PDupwindNthAnti2beta3 = PDupwindNthAntifdOrder82(&beta3[index]); - 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; - } - - /* Calculate temporaries and grid functions */ - ptrdiff_t dir1 = Sign(beta1L); - - ptrdiff_t dir2 = Sign(beta2L); - - 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; - CCTK_REAL_VEC JacPDupwindNthAnti1beta3; - CCTK_REAL_VEC JacPDupwindNthAnti2alpha; - CCTK_REAL_VEC JacPDupwindNthAnti2beta1; - CCTK_REAL_VEC JacPDupwindNthAnti2beta2; - CCTK_REAL_VEC JacPDupwindNthAnti2beta3; - CCTK_REAL_VEC JacPDupwindNthAnti3alpha; - CCTK_REAL_VEC JacPDupwindNthAnti3beta1; - CCTK_REAL_VEC JacPDupwindNthAnti3beta2; - CCTK_REAL_VEC JacPDupwindNthAnti3beta3; - CCTK_REAL_VEC JacPDupwindNthSymm1alpha; - CCTK_REAL_VEC JacPDupwindNthSymm1beta1; - CCTK_REAL_VEC JacPDupwindNthSymm1beta2; - CCTK_REAL_VEC JacPDupwindNthSymm1beta3; - CCTK_REAL_VEC JacPDupwindNthSymm2alpha; - CCTK_REAL_VEC JacPDupwindNthSymm2beta1; - CCTK_REAL_VEC JacPDupwindNthSymm2beta2; - CCTK_REAL_VEC JacPDupwindNthSymm2beta3; - CCTK_REAL_VEC JacPDupwindNthSymm3alpha; - CCTK_REAL_VEC JacPDupwindNthSymm3beta1; - CCTK_REAL_VEC JacPDupwindNthSymm3beta2; - CCTK_REAL_VEC JacPDupwindNthSymm3beta3; - - 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))); - - JacPDupwindNthAnti1beta1 = - kmadd(J11L,PDupwindNthAnti1beta1,kmadd(J21L,PDupwindNthAnti2beta1,kmul(J31L,PDupwindNthAnti3beta1))); - - JacPDupwindNthAnti1beta2 = - kmadd(J11L,PDupwindNthAnti1beta2,kmadd(J21L,PDupwindNthAnti2beta2,kmul(J31L,PDupwindNthAnti3beta2))); - - JacPDupwindNthAnti1beta3 = - kmadd(J11L,PDupwindNthAnti1beta3,kmadd(J21L,PDupwindNthAnti2beta3,kmul(J31L,PDupwindNthAnti3beta3))); - - JacPDupwindNthSymm1alpha = - kmadd(J11L,PDupwindNthSymm1alpha,kmadd(J21L,PDupwindNthSymm2alpha,kmul(J31L,PDupwindNthSymm3alpha))); - - JacPDupwindNthSymm1beta1 = - kmadd(J11L,PDupwindNthSymm1beta1,kmadd(J21L,PDupwindNthSymm2beta1,kmul(J31L,PDupwindNthSymm3beta1))); - - JacPDupwindNthSymm1beta2 = - kmadd(J11L,PDupwindNthSymm1beta2,kmadd(J21L,PDupwindNthSymm2beta2,kmul(J31L,PDupwindNthSymm3beta2))); - - JacPDupwindNthSymm1beta3 = - kmadd(J11L,PDupwindNthSymm1beta3,kmadd(J21L,PDupwindNthSymm2beta3,kmul(J31L,PDupwindNthSymm3beta3))); - - JacPDupwindNthAnti2alpha = - kmadd(J12L,PDupwindNthAnti1alpha,kmadd(J22L,PDupwindNthAnti2alpha,kmul(J32L,PDupwindNthAnti3alpha))); - - JacPDupwindNthAnti2beta1 = - kmadd(J12L,PDupwindNthAnti1beta1,kmadd(J22L,PDupwindNthAnti2beta1,kmul(J32L,PDupwindNthAnti3beta1))); - - JacPDupwindNthAnti2beta2 = - kmadd(J12L,PDupwindNthAnti1beta2,kmadd(J22L,PDupwindNthAnti2beta2,kmul(J32L,PDupwindNthAnti3beta2))); - - JacPDupwindNthAnti2beta3 = - kmadd(J12L,PDupwindNthAnti1beta3,kmadd(J22L,PDupwindNthAnti2beta3,kmul(J32L,PDupwindNthAnti3beta3))); - - JacPDupwindNthSymm2alpha = - kmadd(J12L,PDupwindNthSymm1alpha,kmadd(J22L,PDupwindNthSymm2alpha,kmul(J32L,PDupwindNthSymm3alpha))); - - JacPDupwindNthSymm2beta1 = - kmadd(J12L,PDupwindNthSymm1beta1,kmadd(J22L,PDupwindNthSymm2beta1,kmul(J32L,PDupwindNthSymm3beta1))); - - JacPDupwindNthSymm2beta2 = - kmadd(J12L,PDupwindNthSymm1beta2,kmadd(J22L,PDupwindNthSymm2beta2,kmul(J32L,PDupwindNthSymm3beta2))); - - JacPDupwindNthSymm2beta3 = - kmadd(J12L,PDupwindNthSymm1beta3,kmadd(J22L,PDupwindNthSymm2beta3,kmul(J32L,PDupwindNthSymm3beta3))); - - JacPDupwindNthAnti3alpha = - kmadd(J13L,PDupwindNthAnti1alpha,kmadd(J23L,PDupwindNthAnti2alpha,kmul(J33L,PDupwindNthAnti3alpha))); - - JacPDupwindNthAnti3beta1 = - kmadd(J13L,PDupwindNthAnti1beta1,kmadd(J23L,PDupwindNthAnti2beta1,kmul(J33L,PDupwindNthAnti3beta1))); - - JacPDupwindNthAnti3beta2 = - kmadd(J13L,PDupwindNthAnti1beta2,kmadd(J23L,PDupwindNthAnti2beta2,kmul(J33L,PDupwindNthAnti3beta2))); - - JacPDupwindNthAnti3beta3 = - kmadd(J13L,PDupwindNthAnti1beta3,kmadd(J23L,PDupwindNthAnti2beta3,kmul(J33L,PDupwindNthAnti3beta3))); - - JacPDupwindNthSymm3alpha = - kmadd(J13L,PDupwindNthSymm1alpha,kmadd(J23L,PDupwindNthSymm2alpha,kmul(J33L,PDupwindNthSymm3alpha))); - - JacPDupwindNthSymm3beta1 = - kmadd(J13L,PDupwindNthSymm1beta1,kmadd(J23L,PDupwindNthSymm2beta1,kmul(J33L,PDupwindNthSymm3beta1))); - - JacPDupwindNthSymm3beta2 = - kmadd(J13L,PDupwindNthSymm1beta2,kmadd(J23L,PDupwindNthSymm2beta2,kmul(J33L,PDupwindNthSymm3beta2))); - - JacPDupwindNthSymm3beta3 = - kmadd(J13L,PDupwindNthSymm1beta3,kmadd(J23L,PDupwindNthSymm2beta3,kmul(J33L,PDupwindNthSymm3beta3))); - } - 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; - - JacPDupwindNthAnti1beta2 = PDupwindNthAnti1beta2; - - JacPDupwindNthAnti1beta3 = PDupwindNthAnti1beta3; - - JacPDupwindNthSymm1alpha = PDupwindNthSymm1alpha; - - JacPDupwindNthSymm1beta1 = PDupwindNthSymm1beta1; - - JacPDupwindNthSymm1beta2 = PDupwindNthSymm1beta2; - - JacPDupwindNthSymm1beta3 = PDupwindNthSymm1beta3; - - JacPDupwindNthAnti2alpha = PDupwindNthAnti2alpha; - - JacPDupwindNthAnti2beta1 = PDupwindNthAnti2beta1; - - JacPDupwindNthAnti2beta2 = PDupwindNthAnti2beta2; - - JacPDupwindNthAnti2beta3 = PDupwindNthAnti2beta3; - - JacPDupwindNthSymm2alpha = PDupwindNthSymm2alpha; - - JacPDupwindNthSymm2beta1 = PDupwindNthSymm2beta1; - - JacPDupwindNthSymm2beta2 = PDupwindNthSymm2beta2; - - JacPDupwindNthSymm2beta3 = PDupwindNthSymm2beta3; - - JacPDupwindNthAnti3alpha = PDupwindNthAnti3alpha; - - JacPDupwindNthAnti3beta1 = PDupwindNthAnti3beta1; - - JacPDupwindNthAnti3beta2 = PDupwindNthAnti3beta2; - - JacPDupwindNthAnti3beta3 = PDupwindNthAnti3beta3; - - JacPDupwindNthSymm3alpha = PDupwindNthSymm3alpha; - - JacPDupwindNthSymm3beta1 = PDupwindNthSymm3beta1; - - JacPDupwindNthSymm3beta2 = PDupwindNthSymm3beta2; - - 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 = - kmul(INV(kfmax(rL,ToReal(SpatialBetaDriverRadius))),ToReal(SpatialBetaDriverRadius)); - - CCTK_REAL_VEC theta = - kfmin(ToReal(1),kexp(knmsub(rL,INV(ToReal(SpatialShiftGammaCoeffRadius)),ToReal(1)))); - - CCTK_REAL_VEC dtalpL = - 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),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(phiL,kmul(kmadd(gtu11,JacPDstandardNth1alpha,kmadd(gtu12,JacPDstandardNth2alpha,kmul(gtu13,JacPDstandardNth3alpha))),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(ksub(B1L,Xt1L),ToReal(ShiftBCoeff),kmul(beta1L,kmul(eta,kmul(ToReal(BetaDriver),kadd(ToReal(-1),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),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(phiL,kmul(kmadd(gtu12,JacPDstandardNth1alpha,kmadd(gtu22,JacPDstandardNth2alpha,kmul(gtu23,JacPDstandardNth3alpha))),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(ksub(B2L,Xt2L),ToReal(ShiftBCoeff),kmul(beta2L,kmul(eta,kmul(ToReal(BetaDriver),kadd(ToReal(-1),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),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(phiL,kmul(kmadd(gtu13,JacPDstandardNth1alpha,kmadd(gtu23,JacPDstandardNth2alpha,kmul(gtu33,JacPDstandardNth3alpha))),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(ksub(B3L,Xt3L),ToReal(ShiftBCoeff),kmul(beta3L,kmul(eta,kmul(ToReal(BetaDriver),kadd(ToReal(-1),ToReal(ShiftBCoeff))))))),ToReal(ShiftGammaCoeff))))); - - /* Copy local copies back to grid functions */ - vec_store_partial_prepare(i,lc_imin,lc_imax); - vec_store_nta_partial(dtalp[index],dtalpL); - vec_store_nta_partial(dtbetax[index],dtbetaxL); - vec_store_nta_partial(dtbetay[index],dtbetayL); - vec_store_nta_partial(dtbetaz[index],dtbetazL); - } - LC_ENDLOOP3VEC(ML_BSSN_MP_O8_convertToADMBaseDtLapseShift); -} - -extern "C" void ML_BSSN_MP_O8_convertToADMBaseDtLapseShift(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_MP_O8_convertToADMBaseDtLapseShift_Body"); - } - - if (cctk_iteration % ML_BSSN_MP_O8_convertToADMBaseDtLapseShift_calc_every != ML_BSSN_MP_O8_convertToADMBaseDtLapseShift_calc_offset) - { - return; - } - - const char *const groups[] = { - "ADMBase::dtlapse", - "ADMBase::dtshift", - "grid::coordinates", - "Grid::coordinates", - "ML_BSSN_MP_O8::ML_dtlapse", - "ML_BSSN_MP_O8::ML_dtshift", - "ML_BSSN_MP_O8::ML_Gamma", - "ML_BSSN_MP_O8::ML_lapse", - "ML_BSSN_MP_O8::ML_log_confac", - "ML_BSSN_MP_O8::ML_metric", - "ML_BSSN_MP_O8::ML_shift", - "ML_BSSN_MP_O8::ML_trace_curv"}; - GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_MP_O8_convertToADMBaseDtLapseShift", 12, groups); - - switch(fdOrder) - { - case 2: - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_MP_O8_convertToADMBaseDtLapseShift", 2, 2, 2); - break; - - case 4: - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_MP_O8_convertToADMBaseDtLapseShift", 3, 3, 3); - break; - - case 6: - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_MP_O8_convertToADMBaseDtLapseShift", 4, 4, 4); - break; - - case 8: - GenericFD_EnsureStencilFits(cctkGH, "ML_BSSN_MP_O8_convertToADMBaseDtLapseShift", 5, 5, 5); - break; - } - - GenericFD_LoopOverInterior(cctkGH, ML_BSSN_MP_O8_convertToADMBaseDtLapseShift_Body); - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Leaving ML_BSSN_MP_O8_convertToADMBaseDtLapseShift_Body"); - } -} |