From 4cf4cc528c275dd583ea592ffb3a351f7d4d5f12 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 26 Aug 2022 23:13:39 +0200 Subject: Add 8th order FD implementation. Disabled for now. --- src/qms.c | 153 ++++++++++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 103 insertions(+), 50 deletions(-) diff --git a/src/qms.c b/src/qms.c index d66a614..0634929 100644 --- a/src/qms.c +++ b/src/qms.c @@ -497,6 +497,59 @@ static double diff4_dxz(const double *arr, ptrdiff_t stride, double dx_inv, doub +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; } +#define FD8(arr, stride, h_inv) \ + ((-3.0 * arr[4 * stride] + 32.0 * arr[3 * stride] - 168.0 * arr[2 * stride] + 672.0 * arr[1 * stride] \ + - 672.0 * arr[-1 * stride] + 168.0 * arr[-2 * stride] - 32.0 * arr[-3 * stride] + 3.0 * arr[-4 * stride]) * h_inv / 840.0) +#define FD28(arr, stride, h_inv) \ + ((-9.0 * arr[ 4 * stride] + 128.0 * arr[ 3 * stride] - 1008.0 * arr[ 2 * stride] + 8064.0 * arr[ 1 * stride] \ + -9.0 * arr[-4 * stride] + 128.0 * arr[-3 * stride] - 1008.0 * arr[-2 * stride] + 8064.0 * arr[-1 * stride] \ + - 14350.0 * arr[0]) * SQR(h_inv) / 5040.0) + +static double diff8_dx(const double *arr, ptrdiff_t stride, double h_inv) +{ + return FD8(arr, 1, h_inv); +} +static double diff8_dz(const double *arr, ptrdiff_t stride, double h_inv) +{ + return FD8(arr, stride, h_inv); +} + +static double diff8_dxx(const double *arr, ptrdiff_t stride, double h_inv) +{ + return FD28(arr, 1, h_inv); +} +static double diff8_dzz(const double *arr, ptrdiff_t stride, double h_inv) +{ + return FD28(arr, stride, h_inv); +} +static double diff8_dxz(const double *arr, ptrdiff_t stride, double dx_inv, double dz_inv) +{ + static const double fd_coeffs[] = { 3.0, -32.0, 168.0, -672.0, 0.0, 672.0, -168.0, 32.0, -3.0 }; + static const int stencil = ARRAY_ELEMS(fd_coeffs) / 2; + double val = 0.0; + for (int j = -stencil; j <= stencil; j++) { + double tmp = 0.0; + for (int i = -stencil; i <= stencil; i++) + tmp += fd_coeffs[i + stencil] * arr[j * stride + i]; + val += fd_coeffs[j + stencil] * tmp; + } + return val * dx_inv * dz_inv / SQR(840.0); +} + +#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 +#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 +#endif + static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solver) { cGH *gh = ctx->gh; @@ -571,15 +624,15 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve { gtxy, gtyy, gtyz }, { gtxz, gtyz, gtzz }}; - const double dx_gt11 = diff4_dx(&a_gtxx[idx_src], 1, dx_inv); - const double dx_gt22 = diff4_dx(&a_gtyy[idx_src], 1, dx_inv); - const double dx_gt33 = diff4_dx(&a_gtzz[idx_src], 1, dx_inv); - const double dx_gt13 = diff4_dx(&a_gtxz[idx_src], 1, dx_inv); + 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 dz_gt11 = diff4_dz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_gt22 = diff4_dz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_gt33 = diff4_dz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_gt13 = diff4_dz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv); + 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 dgt[3][3][3] = { { @@ -599,20 +652,20 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve }, }; - const double dxx_gt11 = diff4_dxx(&a_gtxx[idx_src], 1, dx_inv); - const double dxx_gt22 = diff4_dxx(&a_gtyy[idx_src], 1, dx_inv); - const double dxx_gt33 = diff4_dxx(&a_gtzz[idx_src], 1, dx_inv); - const double dxx_gt13 = diff4_dxx(&a_gtxz[idx_src], 1, dx_inv); + 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 dxz_gt11 = diff4_dxz(&a_gtxx[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dxz_gt22 = diff4_dxz(&a_gtyy[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dxz_gt33 = diff4_dxz(&a_gtzz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dxz_gt13 = diff4_dxz(&a_gtxz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); + 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 dzz_gt11 = diff4_dzz(&a_gtxx[idx_src], gh->cctk_lsh[0], dz_inv); - const double dzz_gt22 = diff4_dzz(&a_gtyy[idx_src], gh->cctk_lsh[0], dz_inv); - const double dzz_gt33 = diff4_dzz(&a_gtzz[idx_src], gh->cctk_lsh[0], dz_inv); - const double dzz_gt13 = diff4_dzz(&a_gtxz[idx_src], gh->cctk_lsh[0], dz_inv); + 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 d2gt[3][3][3][3] = { { @@ -688,15 +741,15 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve { Atxy, Atyy, Atyz }, { Atxz, Atyz, Atzz }}; - const double dx_At11 = diff4_dx(&a_Atxx[idx_src], 1, dx_inv); - const double dx_At22 = diff4_dx(&a_Atyy[idx_src], 1, dx_inv); - const double dx_At33 = diff4_dx(&a_Atzz[idx_src], 1, dx_inv); - const double dx_At13 = diff4_dx(&a_Atxz[idx_src], 1, dx_inv); + 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 dz_At11 = diff4_dz(&a_Atxx[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_At22 = diff4_dz(&a_Atyy[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_At33 = diff4_dz(&a_Atzz[idx_src], gh->cctk_lsh[0], dz_inv); - const double dz_At13 = diff4_dz(&a_Atxz[idx_src], gh->cctk_lsh[0], dz_inv); + 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 dAt[3][3][3] = { { @@ -718,14 +771,14 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve const double phi = a_phi[idx_src]; - const double dx_phi = diff4_dx(&a_phi[idx_src], 1, dx_inv); - const double dz_phi = diff4_dz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv); + 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 dxx_phi = diff4_dxx(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv); - const double dzz_phi = diff4_dzz(&a_phi[idx_src], gh->cctk_lsh[0], dz_inv); - const double dxz_phi = diff4_dxz(&a_phi[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); + 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 d2phi[3][3] = { { dxx_phi, 0.0, dxz_phi }, @@ -735,8 +788,8 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve const double trK = a_trK[idx_src]; - const double dx_trk = diff4_dx(&a_trK[idx_src], 1, dx_inv); - const double dz_trk = diff4_dz(&a_trK[idx_src], gh->cctk_lsh[0], dz_inv); + 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 dtrK[3] = { dx_trk, 0.0, dz_trk }; @@ -744,14 +797,14 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve const double Xtz = a_Xtz[idx_src]; const double alpha = a_alpha[idx_src]; - const double dx_alpha = diff4_dx(&a_alpha[idx_src], 1, dx_inv); - const double dz_alpha = diff4_dz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv); + 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 dalpha[3] = { dx_alpha, 0.0, dz_alpha }; - const double dxx_alpha = diff4_dxx(&a_alpha[idx_src], 1, dx_inv); - const double dzz_alpha = diff4_dzz(&a_alpha[idx_src], gh->cctk_lsh[0], dz_inv); - const double dxz_alpha = diff4_dxz(&a_alpha[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); + 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 d2alpha[3][3] = {{ dxx_alpha, 0, dxz_alpha }, { 0, on_axis ? dxx_alpha : dx_alpha / x, 0 }, @@ -762,23 +815,23 @@ static void fill_eq_coeffs(QMSMGContext *ctx, CoordPatch *cp, MG2DContext *solve const double beta[3] = { betax, 0.0, betaz }; - const double dx_betax = diff4_dx(&a_betax[idx_src], 1, dx_inv); - const double dz_betax = diff4_dz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv); + 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_betaz = diff4_dx(&a_betaz[idx_src], 1, dx_inv); - const double dz_betaz = diff4_dz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv); + 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 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 = diff4_dxx(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv); - const double dxz_betax = diff4_dxz(&a_betax[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dzz_betax = diff4_dzz(&a_betax[idx_src], gh->cctk_lsh[0], dz_inv); + 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_betaz = diff4_dxx(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv); - const double dxz_betaz = diff4_dxz(&a_betaz[idx_src], gh->cctk_lsh[0], dx_inv, dz_inv); - const double dzz_betaz = diff4_dzz(&a_betaz[idx_src], gh->cctk_lsh[0], dz_inv); + 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); const double d2beta[3][3][3] = { { -- cgit v1.2.3