From 2dd18bf8dfa32531b20cea95be8d32b3d8177f01 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 27 Aug 2022 16:22:06 +0200 Subject: fill_eq_coeffs(): refactor loading variables/calculating derivatives Prepare the code for templatization/vectorization. --- src/qms.c | 341 +++++++++++++++++++++++++++++++++++++------------------------- 1 file changed, 202 insertions(+), 139 deletions(-) diff --git a/src/qms.c b/src/qms.c index ab803c2..30e1c73 100644 --- a/src/qms.c +++ b/src/qms.c @@ -494,34 +494,54 @@ static void print_stats(QMSMGContext *ms) ms->log_level = orig_log_level; } -#define FD4(arr, stride, h_inv) \ - ((-1.0 * arr[2 * stride] + 8.0 * arr[1 * stride] - 8.0 * arr[-1 * stride] + 1.0 * arr[-2 * stride]) * h_inv / 12.0) -#define FD24(arr, stride, h_inv) \ - ((-1.0 * arr[2 * stride] + 16.0 * arr[1 * stride] - 30.0 * arr[0] + 16.0 * arr[-1 * stride] - 1.0 * arr[-2 * stride]) * SQR(h_inv) / 12.0) +#define LOAD(arr, idx) (arr[idx]) -static double diff4_dx(const double *arr, ptrdiff_t stride, double h_inv) +static inline double fd1_4(const double *arr, const ptrdiff_t stride, const double h_inv) { - return FD4(arr, 1, h_inv); -} -static double diff4_dz(const double *arr, ptrdiff_t stride, double h_inv) -{ - return FD4(arr, stride, h_inv); + const double ap2 = LOAD(arr, 2 * stride); + const double ap1 = LOAD(arr, 1 * stride); + const double am1 = LOAD(arr, -1 * stride); + const double am2 = LOAD(arr, -2 * stride); + return (-ap2 + 8.0 * ap1 - 8.0 * am1 + am2) * h_inv / 12.0; } -static double diff4_dxx(const double *arr, ptrdiff_t stride, double h_inv) -{ - return FD24(arr, 1, h_inv); -} -static double diff4_dzz(const double *arr, ptrdiff_t stride, double h_inv) +static inline double fd2_4(const double *arr, const ptrdiff_t stride, const double h_inv) { - return FD24(arr, stride, h_inv); + const double ap2 = LOAD(arr, 2 * stride); + const double ap1 = LOAD(arr, 1 * stride); + const double a0 = LOAD(arr, 0 * stride); + const double am1 = LOAD(arr, -1 * stride); + const double am2 = LOAD(arr, -2 * stride); + return (-ap2 + 16.0 * ap1 - 30.0 * a0 + 16.0 * am1 - am2) * SQR(h_inv) / 12.0; } -static double diff4_dxz(const double *arr, ptrdiff_t stride, double dx_inv, double dz_inv) + +static inline double fd11_4(const double *arr, ptrdiff_t stride_x, ptrdiff_t stride_z, + double dx_inv, double dz_inv) { - return (-1.0 * (-1.0 * arr[ 2 * stride + 2] + 8.0 * arr[ 2 * stride + 1] - 8.0 * arr[ 2 * stride - 1] + 1.0 * arr[ 2 * stride - 2]) - +8.0 * (-1.0 * arr[ 1 * stride + 2] + 8.0 * arr[ 1 * stride + 1] - 8.0 * arr[ 1 * stride - 1] + 1.0 * arr[ 1 * stride - 2]) - -8.0 * (-1.0 * arr[-1 * stride + 2] + 8.0 * arr[-1 * stride + 1] - 8.0 * arr[-1 * stride - 1] + 1.0 * arr[-1 * stride - 2]) - +1.0 * (-1.0 * arr[-2 * stride + 2] + 8.0 * arr[-2 * stride + 1] - 8.0 * arr[-2 * stride - 1] + 1.0 * arr[-2 * stride - 2])) * dx_inv * dz_inv / 144.0; + const double ap2p2 = LOAD(arr, 2 * stride_x + 2 * stride_z); + const double ap2p1 = LOAD(arr, 2 * stride_x + 1 * stride_z); + const double ap2m1 = LOAD(arr, 2 * stride_x - 1 * stride_z); + const double ap2m2 = LOAD(arr, 2 * stride_x - 2 * stride_z); + + const double ap1p2 = LOAD(arr, 1 * stride_x + 2 * stride_z); + const double ap1p1 = LOAD(arr, 1 * stride_x + 1 * stride_z); + const double ap1m1 = LOAD(arr, 1 * stride_x - 1 * stride_z); + const double ap1m2 = LOAD(arr, 1 * stride_x - 2 * stride_z); + + const double am1p2 = LOAD(arr, -1 * stride_x + 2 * stride_z); + const double am1p1 = LOAD(arr, -1 * stride_x + 1 * stride_z); + const double am1m1 = LOAD(arr, -1 * stride_x - 1 * stride_z); + const double am1m2 = LOAD(arr, -1 * stride_x - 2 * stride_z); + + const double am2p2 = LOAD(arr, -2 * stride_x + 2 * stride_z); + const double am2p1 = LOAD(arr, -2 * stride_x + 1 * stride_z); + const double am2m1 = LOAD(arr, -2 * stride_x - 1 * stride_z); + const double am2m2 = LOAD(arr, -2 * stride_x - 2 * stride_z); + + return (-1.0 * (-1.0 * ap2p2 + 8.0 * ap1p2 - 8.0 * am1p2 + 1.0 * am2p2) + +8.0 * (-1.0 * ap2p1 + 8.0 * ap1p1 - 8.0 * am1p1 + 1.0 * am2p1) + -8.0 * (-1.0 * ap2m1 + 8.0 * ap1m1 - 8.0 * am1m1 + 1.0 * am2m1) + +1.0 * (-1.0 * ap2m2 + 8.0 * ap1m2 - 8.0 * am1m2 + 1.0 * am2m2)) * dx_inv * dz_inv / 144.0; } #define FD8(arr, stride, h_inv) \ @@ -564,19 +584,55 @@ static double diff8_dxz(const double *arr, ptrdiff_t stride, double dx_inv, doub } #if 0 -#define diff_dx diff8_dx -#define diff_dz diff8_dz -#define diff_dxx diff8_dxx -#define diff_dzz diff8_dzz -#define diff_dxz diff8_dxz +#define diff1 fd1_8 +#define diff2 fd2_8 +#define diff11 fd11_8 #else -#define diff_dx diff4_dx -#define diff_dz diff4_dz -#define diff_dxx diff4_dxx -#define diff_dzz diff4_dzz -#define diff_dxz diff4_dxz +#define diff1 fd1_4 +#define diff2 fd2_4 +#define diff11 fd11_4 #endif +#define diff_dx(arr, stride, dx_inv, dz_inv) (diff1 (arr, 1, dx_inv)) +#define diff_dz(arr, stride, dx_inv, dz_inv) (diff1 (arr, stride, dz_inv)) +#define diff_dxx(arr, stride, dx_inv, dz_inv) (diff2 (arr, 1, dx_inv)) +#define diff_dzz(arr, stride, dx_inv, dz_inv) (diff2 (arr, stride, dz_inv)) +#define diff_dxz(arr, stride, dx_inv, dz_inv) (diff11(arr, 1, stride, dx_inv, dz_inv)) + +enum RHSVar { + RHS_VAR_GTXX, + RHS_VAR_GTYY, + RHS_VAR_GTZZ, + RHS_VAR_GTXY, + RHS_VAR_GTXZ, + RHS_VAR_GTYZ, + + RHS_VAR_ATXX, + RHS_VAR_ATYY, + RHS_VAR_ATZZ, + RHS_VAR_ATXY, + RHS_VAR_ATXZ, + RHS_VAR_ATYZ, + + RHS_VAR_BETAX, + RHS_VAR_BETAY, + RHS_VAR_BETAZ, + + RHS_VAR_XTX, + RHS_VAR_XTY, + RHS_VAR_XTZ, + + RHS_VAR_TRK, + RHS_VAR_PHI, + RHS_VAR_ALPHA, + + RHS_VAR_X, + RHS_VAR_Y, + RHS_VAR_Z, + + RHS_VAR_NB, +}; + static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solver) { cGH *gh = ctx->gh; @@ -584,35 +640,36 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve const ptrdiff_t stride_z = CCTK_GFINDEX3D(gh, 0, 0, 1); - double *a_x = CCTK_VarDataPtr(gh, 0, "grid::x"); - double *a_z = CCTK_VarDataPtr(gh, 0, "grid::z"); - - const double dx = a_x[1] - a_x[0]; - const double dz = a_z[stride_z] - a_z[0]; - const double dx_inv = 1.0 / dx; - const double dz_inv = 1.0 / dz; - - double *a_gtxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt11"); - double *a_gtyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt22"); - double *a_gtzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt33"); - double *a_gtxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt12"); - double *a_gtxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt13"); - double *a_gtyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt23"); - double *a_Atxx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At11"); - double *a_Atyy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At22"); - double *a_Atzz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At33"); - double *a_Atxy = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At12"); - double *a_Atxz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At13"); - double *a_Atyz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At23"); - double *a_trK = CCTK_VarDataPtr(gh, 0, "ML_BSSN::trK"); - double *a_phi = CCTK_VarDataPtr(gh, 0, "ML_BSSN::phi"); - double *a_Xtx = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt1"); - double *a_Xtz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt3"); - - double *a_alpha = CCTK_VarDataPtr(gh, 0, "ML_BSSN::alpha"); - - double *a_betax = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta1"); - double *a_betaz = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta3"); + const double * const vars[] = { + [RHS_VAR_GTXX] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt11"), + [RHS_VAR_GTYY] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt22"), + [RHS_VAR_GTZZ] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt33"), + [RHS_VAR_GTXY] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt12"), + [RHS_VAR_GTXZ] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt13"), + [RHS_VAR_GTYZ] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::gt23"), + + [RHS_VAR_ATXX] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At11"), + [RHS_VAR_ATYY] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At22"), + [RHS_VAR_ATZZ] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At33"), + [RHS_VAR_ATXY] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At12"), + [RHS_VAR_ATXZ] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At13"), + [RHS_VAR_ATYZ] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::At23"), + + [RHS_VAR_TRK] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::trK"), + [RHS_VAR_PHI] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::phi"), + [RHS_VAR_ALPHA] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::alpha"), + + [RHS_VAR_XTX] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt1"), + [RHS_VAR_XTZ] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::Xt3"), + + [RHS_VAR_BETAX] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta1"), + [RHS_VAR_BETAZ] = CCTK_VarDataPtr(gh, 0, "ML_BSSN::beta3"), + [RHS_VAR_X] = CCTK_VarDataPtr(gh, 0, "grid::x"), + [RHS_VAR_Z] = CCTK_VarDataPtr(gh, 0, "grid::z"), + }; + + const double dx_inv = 1.0 / (vars[RHS_VAR_X][1] - vars[RHS_VAR_X][0]); + const double dz_inv = 1.0 / (vars[RHS_VAR_Z][stride_z] - vars[RHS_VAR_Z][0]); solver->diff_coeffs[MG2D_DIFF_COEFF_20]->boundaries[MG2D_BOUNDARY_0L].flags = MG2D_DC_FLAG_DISCONT; solver->diff_coeffs[MG2D_DIFF_COEFF_10]->boundaries[MG2D_BOUNDARY_0L].flags = MG2D_DC_FLAG_POLE; @@ -624,7 +681,7 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve const int idx_dc = idx_z * solver->diff_coeffs[0]->stride + idx_x; const int idx_rhs = idx_z * solver->rhs_stride + idx_x; - const double x = a_x[idx_src]; + const double x = vars[RHS_VAR_X][idx_src]; const int on_axis = fabs(x) < 1e-13; double K[3][3], Km[3][3], Ku[3][3], dK[3][3][3]; @@ -640,26 +697,33 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve double betal[3], dbetal[3][3], d2betal[3][3][3]; - const double gtxx = a_gtxx[idx_src]; - const double gtyy = a_gtyy[idx_src]; - const double gtzz = a_gtzz[idx_src]; - const double gtxy = a_gtxy[idx_src]; - const double gtxz = a_gtxz[idx_src]; - const double gtyz = a_gtyz[idx_src]; +#define LOAD_VAR(var) ( vars[RHS_VAR_ ## var][idx_src]) +#define DX(var) (diff1 (&vars[RHS_VAR_ ## var][idx_src], 1, dx_inv )) +#define DZ(var) (diff1 (&vars[RHS_VAR_ ## var][idx_src], stride_z, dz_inv)) +#define DXX(var) (diff2 (&vars[RHS_VAR_ ## var][idx_src], 1, dx_inv)) +#define DZZ(var) (diff2 (&vars[RHS_VAR_ ## var][idx_src], stride_z, dz_inv)) +#define DXZ(var) (diff11(&vars[RHS_VAR_ ## var][idx_src], 1, stride_z, dx_inv, dz_inv)) + + const double gtxx = LOAD_VAR(GTXX); + const double gtyy = LOAD_VAR(GTYY); + const double gtzz = LOAD_VAR(GTZZ); + const double gtxy = LOAD_VAR(GTXY); + const double gtxz = LOAD_VAR(GTXZ); + const double gtyz = LOAD_VAR(GTYZ); const double gt[3][3] = {{ gtxx, gtxy, gtxz }, { gtxy, gtyy, gtyz }, { gtxz, gtyz, gtzz }}; - const double dx_gt11 = diff_dx(&a_gtxx[idx_src], 1, dx_inv); - const double dx_gt22 = diff_dx(&a_gtyy[idx_src], 1, dx_inv); - const double dx_gt33 = diff_dx(&a_gtzz[idx_src], 1, dx_inv); - const double dx_gt13 = diff_dx(&a_gtxz[idx_src], 1, dx_inv); + const double dx_gt11 = DX(GTXX); + const double dx_gt22 = DX(GTYY); + const double dx_gt33 = DX(GTZZ); + const double dx_gt13 = DX(GTXZ); - const double dz_gt11 = diff_dz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_gt22 = diff_dz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_gt33 = diff_dz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_gt13 = diff_dz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv); + const double dz_gt11 = DZ(GTXX); + const double dz_gt22 = DZ(GTYY); + const double dz_gt33 = DZ(GTZZ); + const double dz_gt13 = DZ(GTXZ); const double dgt[3][3][3] = { { @@ -679,20 +743,20 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve }, }; - const double dxx_gt11 = diff_dxx(&a_gtxx[idx_src], 1, dx_inv); - const double dxx_gt22 = diff_dxx(&a_gtyy[idx_src], 1, dx_inv); - const double dxx_gt33 = diff_dxx(&a_gtzz[idx_src], 1, dx_inv); - const double dxx_gt13 = diff_dxx(&a_gtxz[idx_src], 1, dx_inv); + const double dxx_gt11 = DXX(GTXX); + const double dxx_gt22 = DXX(GTYY); + const double dxx_gt33 = DXX(GTZZ); + const double dxx_gt13 = DXX(GTXZ); - const double dxz_gt11 = diff_dxz(&a_gtxx[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dxz_gt22 = diff_dxz(&a_gtyy[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dxz_gt33 = diff_dxz(&a_gtzz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dxz_gt13 = diff_dxz(&a_gtxz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); + const double dxz_gt11 = DXZ(GTXX); + const double dxz_gt22 = DXZ(GTYY); + const double dxz_gt33 = DXZ(GTZZ); + const double dxz_gt13 = DXZ(GTXZ); - const double dzz_gt11 = diff_dzz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv); - const double dzz_gt22 = diff_dzz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv); - const double dzz_gt33 = diff_dzz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv); - const double dzz_gt13 = diff_dzz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv); + const double dzz_gt11 = DZZ(GTXX); + const double dzz_gt22 = DZZ(GTYY); + const double dzz_gt33 = DZZ(GTZZ); + const double dzz_gt13 = DZZ(GTXZ); const double d2gt[3][3][3][3] = { { @@ -756,27 +820,27 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve }, }; - const double Atxx = a_Atxx[idx_src]; - const double Atyy = a_Atyy[idx_src]; - const double Atzz = a_Atzz[idx_src]; - const double Atxy = a_Atxy[idx_src]; - const double Atxz = a_Atxz[idx_src]; - const double Atyz = a_Atyz[idx_src]; + const double Atxx = LOAD_VAR(ATXX); + const double Atyy = LOAD_VAR(ATYY); + const double Atzz = LOAD_VAR(ATZZ); + const double Atxy = LOAD_VAR(ATXY); + const double Atxz = LOAD_VAR(ATXZ); + const double Atyz = LOAD_VAR(ATYZ); const double At[3][3] = { { Atxx, Atxy, Atxz }, { Atxy, Atyy, Atyz }, { Atxz, Atyz, Atzz }}; - const double dx_At11 = diff_dx(&a_Atxx[idx_src], 1, dx_inv); - const double dx_At22 = diff_dx(&a_Atyy[idx_src], 1, dx_inv); - const double dx_At33 = diff_dx(&a_Atzz[idx_src], 1, dx_inv); - const double dx_At13 = diff_dx(&a_Atxz[idx_src], 1, dx_inv); + const double dx_At11 = DX(ATXX); + const double dx_At22 = DX(ATYY); + const double dx_At33 = DX(ATZZ); + const double dx_At13 = DX(ATXZ); - const double dz_At11 = diff_dz(&a_Atxx[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_At22 = diff_dz(&a_Atyy[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_At33 = diff_dz(&a_Atzz[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_At13 = diff_dz(&a_Atxz[idx_src], gh->cctk_lsh[0], dz_inv); + const double dz_At11 = DZ(ATXX); + const double dz_At22 = DZ(ATYY); + const double dz_At33 = DZ(ATZZ); + const double dz_At13 = DZ(ATXZ); const double dAt[3][3][3] = { { @@ -796,69 +860,68 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve }, }; - const double phi = a_phi[idx_src]; - - const double dx_phi = diff_dx(&a_phi[idx_src], 1, dx_inv); - const double dz_phi = diff_dz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv); - - const double dphi[3] = { dx_phi, 0.0, dz_phi }; + const double phi = LOAD_VAR(PHI); + const double dphi[3] = { DX(PHI), 0.0, DZ(PHI) }; - const double dxx_phi = diff_dxx(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv); - const double dzz_phi = diff_dzz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv); - const double dxz_phi = diff_dxz(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); + const double dxx_phi = DXX(PHI); + const double dzz_phi = DZZ(PHI); + const double dxz_phi = DXZ(PHI); const double d2phi[3][3] = { - { dxx_phi, 0.0, dxz_phi }, - { 0.0, on_axis ? dxx_phi : dx_phi / x, 0.0 }, - { dxz_phi, 0.0, dzz_phi }, + { dxx_phi, 0.0, dxz_phi }, + { 0.0, on_axis ? dxx_phi : dphi[0] / x, 0.0 }, + { dxz_phi, 0.0, dzz_phi }, }; - const double trK = a_trK[idx_src]; + const double trK = LOAD_VAR(TRK); + const double dtrK[3] = { DX(TRK), 0.0, DZ(TRK) }; - const double dx_trk = diff_dx(&a_trK[idx_src], 1, dx_inv); - const double dz_trk = diff_dz(&a_trK[idx_src], gh->cctk_lsh[0], dz_inv); + const double Xtx = LOAD_VAR(XTX); + const double Xtz = LOAD_VAR(XTZ); - const double dtrK[3] = { dx_trk, 0.0, dz_trk }; - - const double Xtx = a_Xtx[idx_src]; - const double Xtz = a_Xtz[idx_src]; - - const double alpha = a_alpha[idx_src]; - const double dx_alpha = diff_dx(&a_alpha[idx_src], 1, dx_inv); - const double dz_alpha = diff_dz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv); + const double alpha = LOAD_VAR(ALPHA); + const double dx_alpha = DX(ALPHA); + const double dz_alpha = DZ(ALPHA); const double dalpha[3] = { dx_alpha, 0.0, dz_alpha }; - const double dxx_alpha = diff_dxx(&a_alpha[idx_src], 1, dx_inv); - const double dzz_alpha = diff_dzz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv); - const double dxz_alpha = diff_dxz(&a_alpha[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); + const double dxx_alpha = DXX(ALPHA); + const double dzz_alpha = DZZ(ALPHA); + const double dxz_alpha = DXZ(ALPHA); const double d2alpha[3][3] = {{ dxx_alpha, 0, dxz_alpha }, { 0, on_axis ? dxx_alpha : dx_alpha / x, 0 }, { dxz_alpha, 0, dzz_alpha }}; - const double betax = a_betax[idx_src]; - const double betaz = a_betaz[idx_src]; + const double betax = LOAD_VAR(BETAX); + const double betaz = LOAD_VAR(BETAX); const double beta[3] = { betax, 0.0, betaz }; - const double dx_betax = diff_dx(&a_betax[idx_src], 1, dx_inv); - const double dz_betax = diff_dz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv); + const double dx_betax = DX(BETAX); + const double dz_betax = DZ(BETAX); - const double dx_betaz = diff_dx(&a_betaz[idx_src], 1, dx_inv); - const double dz_betaz = diff_dz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv); + const double dx_betaz = DX(BETAZ); + const double dz_betaz = DZ(BETAZ); const double dbeta[3][3] = {{ dx_betax, 0.0, dx_betaz }, { 0.0, on_axis ? dx_betax : betax / x, 0.0 }, { dz_betax, 0.0, dz_betaz }}; - const double dxx_betax = diff_dxx(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv); - const double dxz_betax = diff_dxz(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dzz_betax = diff_dzz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv); + const double dxx_betax = DXX(BETAX); + const double dxz_betax = DXZ(BETAX); + const double dzz_betax = DZZ(BETAX); + + const double dxx_betaz = DXX(BETAZ); + const double dxz_betaz = DXZ(BETAZ); + const double dzz_betaz = DZZ(BETAZ); - const double dxx_betaz = diff_dxx(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv); - const double dxz_betaz = diff_dxz(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dzz_betaz = diff_dzz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv); +#undef LOAD_VAR +#undef DX +#undef DZ +#undef DXX +#undef DZZ +#undef DXZ const double d2beta[3][3][3] = { { -- cgit v1.2.3