diff --git a/AMSS_NCKU_source/bssn_rhs_cuda.cu b/AMSS_NCKU_source/bssn_rhs_cuda.cu index 2e18a48..2ee8a8b 100644 --- a/AMSS_NCKU_source/bssn_rhs_cuda.cu +++ b/AMSS_NCKU_source/bssn_rhs_cuda.cu @@ -401,6 +401,8 @@ __device__ __forceinline__ double fetch_sym_ord3_direct(const double *src, + (skF - 1) * d_gp.ex[0] * d_gp.ex[1]]; } +#include "fd_cuda_helpers.cuh" + /* ------------------------------------------------------------------ */ /* GPU buffer management */ /* ------------------------------------------------------------------ */ @@ -1729,45 +1731,10 @@ void kern_fderivs_batched(FDerivTables tables, int field_count) const int jF = j0 + 1; const int kF = k0 + 1; - if ((iF + 2) <= imaxF && (iF - 2) >= iminF && - (jF + 2) <= jmaxF && (jF - 2) >= jminF && - (kF + 2) <= kmaxF && (kF - 2) >= kminF) - { - fx[tid] = d_gp.d12dx * ( - fetch_sym_ord2_direct(src, iF - 2, jF, kF, SoA0, SoA1, SoA2) - - 8.0 * fetch_sym_ord2_direct(src, iF - 1, jF, kF, SoA0, SoA1, SoA2) - + 8.0 * fetch_sym_ord2_direct(src, iF + 1, jF, kF, SoA0, SoA1, SoA2) - - fetch_sym_ord2_direct(src, iF + 2, jF, kF, SoA0, SoA1, SoA2)); - fy[tid] = d_gp.d12dy * ( - fetch_sym_ord2_direct(src, iF, jF - 2, kF, SoA0, SoA1, SoA2) - - 8.0 * fetch_sym_ord2_direct(src, iF, jF - 1, kF, SoA0, SoA1, SoA2) - + 8.0 * fetch_sym_ord2_direct(src, iF, jF + 1, kF, SoA0, SoA1, SoA2) - - fetch_sym_ord2_direct(src, iF, jF + 2, kF, SoA0, SoA1, SoA2)); - fz[tid] = d_gp.d12dz * ( - fetch_sym_ord2_direct(src, iF, jF, kF - 2, SoA0, SoA1, SoA2) - - 8.0 * fetch_sym_ord2_direct(src, iF, jF, kF - 1, SoA0, SoA1, SoA2) - + 8.0 * fetch_sym_ord2_direct(src, iF, jF, kF + 1, SoA0, SoA1, SoA2) - - fetch_sym_ord2_direct(src, iF, jF, kF + 2, SoA0, SoA1, SoA2)); - } - else if ((iF + 1) <= imaxF && (iF - 1) >= iminF && - (jF + 1) <= jmaxF && (jF - 1) >= jminF && - (kF + 1) <= kmaxF && (kF - 1) >= kminF) - { - fx[tid] = d_gp.d2dx * ( - -fetch_sym_ord2_direct(src, iF - 1, jF, kF, SoA0, SoA1, SoA2) - +fetch_sym_ord2_direct(src, iF + 1, jF, kF, SoA0, SoA1, SoA2)); - fy[tid] = d_gp.d2dy * ( - -fetch_sym_ord2_direct(src, iF, jF - 1, kF, SoA0, SoA1, SoA2) - +fetch_sym_ord2_direct(src, iF, jF + 1, kF, SoA0, SoA1, SoA2)); - fz[tid] = d_gp.d2dz * ( - -fetch_sym_ord2_direct(src, iF, jF, kF - 1, SoA0, SoA1, SoA2) - +fetch_sym_ord2_direct(src, iF, jF, kF + 1, SoA0, SoA1, SoA2)); - } - else { - fx[tid] = 0.0; - fy[tid] = 0.0; - fz[tid] = 0.0; - } + fd_compute_first3(src, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + SoA0, SoA1, SoA2, + fx[tid], fy[tid], fz[tid]); } __global__ __launch_bounds__(128, 4) @@ -1807,6 +1774,12 @@ void kern_fdderivs_batched(FDDerivTables tables, int field_count) const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + fd_compute_second6(src, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + SoA0, SoA1, SoA2, + fxx[tid], fxy[tid], fxz[tid], fyy[tid], fyz[tid], fzz[tid]); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -1934,12 +1907,43 @@ void kern_fdderivs_batched(FDDerivTables tables, int field_count) fxx[tid] = 0.0; fxy[tid] = 0.0; fxz[tid] = 0.0; fyy[tid] = 0.0; fyz[tid] = 0.0; fzz[tid] = 0.0; } +#endif } +static void gpu_fderivs_batch(int field_count, + double *const *src_fields, + double *const *fx_fields, + double *const *fy_fields, + double *const *fz_fields, + const int *soa_signs, + int all); +static void gpu_fdderivs_batch(int field_count, + double *const *src_fields, + double *const *fxx_fields, + double *const *fxy_fields, + double *const *fxz_fields, + double *const *fyy_fields, + double *const *fyz_fields, + double *const *fzz_fields, + const int *soa_signs, + int all); +static void gpu_lopsided_kodis_single_batch(double *d_f_adv, double *d_f_ko, double *d_f_rhs, + double *d_Sfx, double *d_Sfy, double *d_Sfz, + double SoA0, double SoA1, double SoA2, + double eps_val, int all); + /* symmetry_bd on GPU for ord=2, then launch fderivs kernel */ static void gpu_fderivs(double *d_f, double *d_fx, double *d_fy, double *d_fz, double SoA0, double SoA1, double SoA2, int all) { +#if ghost_width != 3 + double *src_fields[1] = {d_f}; + double *fx_fields[1] = {d_fx}; + double *fy_fields[1] = {d_fy}; + double *fz_fields[1] = {d_fz}; + const int soa_signs[3] = {(int)SoA0, (int)SoA1, (int)SoA2}; + gpu_fderivs_batch(1, src_fields, fx_fields, fy_fields, fz_fields, soa_signs, all); +#else double *fh = g_buf.d_fh2; const size_t nx = (size_t)g_buf.prev_nx; const size_t ny = (size_t)g_buf.prev_ny; @@ -1948,6 +1952,7 @@ static void gpu_fderivs(double *d_f, double *d_fx, double *d_fy, double *d_fz, kern_symbd_pack_ord2<<>>(d_f, fh, SoA0, SoA1, SoA2); kern_fderivs<<>>(fh, d_fx, d_fy, d_fz); +#endif } /* symmetry_bd on GPU for ord=2, then launch fdderivs kernel */ @@ -1956,6 +1961,18 @@ static void gpu_fdderivs(double *d_f, double *d_fyy, double *d_fyz, double *d_fzz, double SoA0, double SoA1, double SoA2, int all) { +#if ghost_width != 3 + double *src_fields[1] = {d_f}; + double *fxx_fields[1] = {d_fxx}; + double *fxy_fields[1] = {d_fxy}; + double *fxz_fields[1] = {d_fxz}; + double *fyy_fields[1] = {d_fyy}; + double *fyz_fields[1] = {d_fyz}; + double *fzz_fields[1] = {d_fzz}; + const int soa_signs[3] = {(int)SoA0, (int)SoA1, (int)SoA2}; + gpu_fdderivs_batch(1, src_fields, fxx_fields, fxy_fields, fxz_fields, + fyy_fields, fyz_fields, fzz_fields, soa_signs, all); +#else double *fh = g_buf.d_fh2; const size_t nx = (size_t)g_buf.prev_nx; const size_t ny = (size_t)g_buf.prev_ny; @@ -1964,6 +1981,7 @@ static void gpu_fdderivs(double *d_f, kern_symbd_pack_ord2<<>>(d_f, fh, SoA0, SoA1, SoA2); kern_fdderivs<<>>(fh, d_fxx, d_fxy, d_fxz, d_fyy, d_fyz, d_fzz); +#endif } static void gpu_fderivs_batch(int field_count, @@ -2053,6 +2071,12 @@ void kern_phase10_ricci_batched(const double * __restrict__ gupxx, const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + fd_compute_second6(src, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + SoA0, SoA1, SoA2, + fxx, fxy, fxz, fyy, fyz, fzz); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -2176,6 +2200,7 @@ void kern_phase10_ricci_batched(const double * __restrict__ gupxx, - fetch_sym_ord2_direct(src, iF, jF - 1, kF + 1, SoA0, SoA1, SoA2) + fetch_sym_ord2_direct(src, iF, jF + 1, kF + 1, SoA0, SoA1, SoA2)); } +#endif } dst[tid] = gupxx[tid] * fxx + gupyy[tid] * fyy + gupzz[tid] * fzz @@ -2240,6 +2265,16 @@ void kern_phase14_lap_chi_derivs(const double * __restrict__ Lap, const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + fd_compute_second6(Lap, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + 1, 1, 1, + fxx[tid], fxy[tid], fxz[tid], fyy[tid], fyz[tid], fzz[tid]); + fd_compute_first3(chi, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + 1, 1, 1, + chix_out[tid], chiy_out[tid], chiz_out[tid]); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -2393,6 +2428,7 @@ void kern_phase14_lap_chi_derivs(const double * __restrict__ Lap, fyy[tid] = 0.0; fyz[tid] = 0.0; fzz[tid] = 0.0; chix_out[tid] = 0.0; chiy_out[tid] = 0.0; chiz_out[tid] = 0.0; } +#endif } /* Combined ord=3 advection + KO dissipation. @@ -2404,6 +2440,11 @@ static void gpu_lopsided_kodis(double *d_f_adv, double *d_f_ko, double *d_f_rhs, double SoA0, double SoA1, double SoA2, double eps_val, int all) { +#if ghost_width != 3 + gpu_lopsided_kodis_single_batch(d_f_adv, d_f_ko, d_f_rhs, + d_Sfx, d_Sfy, d_Sfz, + SoA0, SoA1, SoA2, eps_val, all); +#else double *fh = g_buf.d_fh3; const size_t nx = (size_t)g_buf.prev_nx; const size_t ny = (size_t)g_buf.prev_ny; @@ -2419,6 +2460,7 @@ static void gpu_lopsided_kodis(double *d_f_adv, double *d_f_ko, double *d_f_rhs, } kern_kodis<<>>(fh, d_f_rhs, eps_val); } +#endif } __global__ __launch_bounds__(128, 4) @@ -2449,6 +2491,22 @@ void kern_lopsided_kodis_batched(const double * __restrict__ Sfx, const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + if (i0 <= nx - 2 && j0 <= ny - 2 && k0 <= nz - 2) { + const double val = + fd_lopsided_axis(adv_src, iF, jF, kF, 0, Sfx[tid], iF, iminF, imaxF, + d_gp.dX, SoA0, SoA1, SoA2) + + fd_lopsided_axis(adv_src, iF, jF, kF, 1, Sfy[tid], jF, jminF, jmaxF, + d_gp.dY, SoA0, SoA1, SoA2) + + fd_lopsided_axis(adv_src, iF, jF, kF, 2, Sfz[tid], kF, kminF, kmaxF, + d_gp.dZ, SoA0, SoA1, SoA2); + rhs[tid] += val; + } + + rhs[tid] += fd_ko_term(ko_src, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + eps_val, SoA0, SoA1, SoA2); +#else if (i0 <= nx - 2 && j0 <= ny - 2 && k0 <= nz - 2) { double val = 0.0; @@ -2631,6 +2689,25 @@ void kern_lopsided_kodis_batched(const double * __restrict__ Sfx, rhs[tid] += (eps_val / cof) * (Dx / d_gp.dX + Dy / d_gp.dY + Dz / d_gp.dZ); } +#endif +} + +static void gpu_lopsided_kodis_single_batch(double *d_f_adv, double *d_f_ko, double *d_f_rhs, + double *d_Sfx, double *d_Sfy, double *d_Sfz, + double SoA0, double SoA1, double SoA2, + double eps_val, int all) +{ + LopsidedKodisTables tables = {}; + tables.adv_fields[0] = d_f_adv; + tables.ko_fields[0] = d_f_ko; + tables.rhs_fields[0] = d_f_rhs; + tables.soa_signs[0] = (int)SoA0; + tables.soa_signs[1] = (int)SoA1; + tables.soa_signs[2] = (int)SoA2; + + dim3 launch_grid((unsigned int)grid((size_t)all), 1u); + kern_lopsided_kodis_batched<<>>( + d_Sfx, d_Sfy, d_Sfz, tables, eps_val); } static void gpu_lopsided_kodis_state_batch(double eps_val, int all, bool include_escalar = false) @@ -4624,6 +4701,12 @@ void kern_phase12_13_chi_correction_fused( const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + fd_compute_second6(chi, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + 1, 1, 1, + cxx, cxy, cxz, cyy, cyz, czz); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -4747,6 +4830,7 @@ void kern_phase12_13_chi_correction_fused( - fetch_sym_ord2_direct(chi, iF, jF - 1, kF + 1, 1, 1, 1) + fetch_sym_ord2_direct(chi, iF, jF + 1, kF + 1, 1, 1, 1)); } +#endif } const double cx = chix[tid]; @@ -4917,6 +5001,12 @@ void kern_phase15_trK_Aij_gauge( double fyy_v = 0.0, fyz_v = 0.0, fzz_v = 0.0; if (!(i0 > nx - 2 || j0 > ny - 2 || k0 > nz - 2)) { +#if ghost_width != 3 + fd_compute_second6(alpn1, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + 1, 1, 1, + fxx_v, fxy_v, fxz_v, fyy_v, fyz_v, fzz_v); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -5040,6 +5130,7 @@ void kern_phase15_trK_Aij_gauge( - fetch_sym_ord2_direct(alpn1, iF, jF - 1, kF + 1, 1, 1, 1) + fetch_sym_ord2_direct(alpn1, iF, jF + 1, kF + 1, 1, 1, 1)); } +#endif } /* raised chi/chi */ @@ -5443,15 +5534,15 @@ static void setup_grid_params(int *ex, gp.imaxF = nx; gp.jmaxF = ny; gp.kmaxF = nz; - if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) gp.kminF = -1; - if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) gp.iminF = -1; - if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) gp.jminF = -1; + if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) gp.kminF = 2 - ghost_width; + if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) gp.iminF = 2 - ghost_width; + if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) gp.jminF = 2 - ghost_width; gp.iminF3 = 1; gp.jminF3 = 1; gp.kminF3 = 1; - if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) gp.kminF3 = -2; - if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) gp.iminF3 = -2; - if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) gp.jminF3 = -2; + if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) gp.kminF3 = 1 - ghost_width; + if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) gp.iminF3 = 1 - ghost_width; + if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) gp.jminF3 = 1 - ghost_width; gp.Symmetry = Symmetry; gp.eps = eps; gp.co = co; diff --git a/AMSS_NCKU_source/fd_cuda_helpers.cuh b/AMSS_NCKU_source/fd_cuda_helpers.cuh new file mode 100644 index 0000000..b43129f --- /dev/null +++ b/AMSS_NCKU_source/fd_cuda_helpers.cuh @@ -0,0 +1,412 @@ +#ifndef AMSS_NCKU_FD_CUDA_HELPERS_CUH +#define AMSS_NCKU_FD_CUDA_HELPERS_CUH + +#ifndef ghost_width +#error "ghost_width must be defined before including fd_cuda_helpers.cuh" +#endif + +#if ghost_width < 2 || ghost_width > 5 +#error "CUDA finite-difference helpers support ghost_width 2..5" +#endif + +#define AMSS_FD_CENTER_RADIUS (ghost_width - 1) +#define AMSS_FD_LK_RADIUS (ghost_width) + +__device__ __forceinline__ int fd_axis_radius(int qF, int qminF, int qmaxF) +{ +#if AMSS_FD_CENTER_RADIUS >= 4 + if (qF - 4 >= qminF && qF + 4 <= qmaxF) return 4; +#endif +#if AMSS_FD_CENTER_RADIUS >= 3 + if (qF - 3 >= qminF && qF + 3 <= qmaxF) return 3; +#endif +#if AMSS_FD_CENTER_RADIUS >= 2 + if (qF - 2 >= qminF && qF + 2 <= qmaxF) return 2; +#endif + if (qF - 1 >= qminF && qF + 1 <= qmaxF) return 1; + return 0; +} + +__device__ __forceinline__ int fd_common_radius(int iF, int jF, int kF, + int iminF, int jminF, int kminF, + int imaxF, int jmaxF, int kmaxF) +{ + int r = fd_axis_radius(iF, iminF, imaxF); + const int ry = fd_axis_radius(jF, jminF, jmaxF); + const int rz = fd_axis_radius(kF, kminF, kmaxF); + if (ry < r) r = ry; + if (rz < r) r = rz; + return r; +} + +__device__ __forceinline__ double fd_first_coef(int r, int off) +{ + switch (r) { + case 1: + if (off == -1) return -1.0; + if (off == 1) return 1.0; + return 0.0; + case 2: + if (off == -2) return 1.0; + if (off == -1) return -8.0; + if (off == 1) return 8.0; + if (off == 2) return -1.0; + return 0.0; + case 3: + if (off == -3) return -1.0; + if (off == -2) return 9.0; + if (off == -1) return -45.0; + if (off == 1) return 45.0; + if (off == 2) return -9.0; + if (off == 3) return 1.0; + return 0.0; + case 4: + if (off == -4) return 3.0; + if (off == -3) return -32.0; + if (off == -2) return 168.0; + if (off == -1) return -672.0; + if (off == 1) return 672.0; + if (off == 2) return -168.0; + if (off == 3) return 32.0; + if (off == 4) return -3.0; + return 0.0; + default: + return 0.0; + } +} + +__device__ __forceinline__ double fd_second_coef(int r, int off) +{ + switch (r) { + case 1: + if (off == -1) return 1.0; + if (off == 0) return -2.0; + if (off == 1) return 1.0; + return 0.0; + case 2: + if (off == -2) return -1.0; + if (off == -1) return 16.0; + if (off == 0) return -30.0; + if (off == 1) return 16.0; + if (off == 2) return -1.0; + return 0.0; + case 3: + if (off == -3) return 2.0; + if (off == -2) return -27.0; + if (off == -1) return 270.0; + if (off == 0) return -490.0; + if (off == 1) return 270.0; + if (off == 2) return -27.0; + if (off == 3) return 2.0; + return 0.0; + case 4: + if (off == -4) return -9.0; + if (off == -3) return 128.0; + if (off == -2) return -1008.0; + if (off == -1) return 8064.0; + if (off == 0) return -14350.0; + if (off == 1) return 8064.0; + if (off == 2) return -1008.0; + if (off == 3) return 128.0; + if (off == 4) return -9.0; + return 0.0; + default: + return 0.0; + } +} + +__device__ __forceinline__ double fd_first_denom(int r) +{ + return (r == 4) ? 840.0 : ((r == 3) ? 60.0 : ((r == 2) ? 12.0 : 2.0)); +} + +__device__ __forceinline__ double fd_second_denom(int r) +{ + return (r == 4) ? 5040.0 : ((r == 3) ? 180.0 : ((r == 2) ? 12.0 : 1.0)); +} + +__device__ __forceinline__ double fd_fetch_axis(const double *src, + int iF, int jF, int kF, + int axis, int off, + int SoA0, int SoA1, int SoA2) +{ + if (axis == 0) iF += off; + else if (axis == 1) jF += off; + else kF += off; + return fetch_sym_ord2_direct(src, iF, jF, kF, SoA0, SoA1, SoA2); +} + +__device__ __forceinline__ double fd_fetch_axis2(const double *src, + int iF, int jF, int kF, + int axis_a, int off_a, + int axis_b, int off_b, + int SoA0, int SoA1, int SoA2) +{ + if (axis_a == 0) iF += off_a; + else if (axis_a == 1) jF += off_a; + else kF += off_a; + if (axis_b == 0) iF += off_b; + else if (axis_b == 1) jF += off_b; + else kF += off_b; + return fetch_sym_ord2_direct(src, iF, jF, kF, SoA0, SoA1, SoA2); +} + +__device__ __forceinline__ double fd_first_axis_radius(const double *src, + int iF, int jF, int kF, + int axis, int r, double h, + int SoA0, int SoA1, int SoA2) +{ + if (r <= 0) return 0.0; + double s = 0.0; +#pragma unroll + for (int off = -4; off <= 4; ++off) { + const double c = fd_first_coef(r, off); + if (c != 0.0) { + s += c * fd_fetch_axis(src, iF, jF, kF, axis, off, SoA0, SoA1, SoA2); + } + } + return s / (fd_first_denom(r) * h); +} + +__device__ __forceinline__ double fd_second_axis_radius(const double *src, + int iF, int jF, int kF, + int axis, int r, double h, + int SoA0, int SoA1, int SoA2) +{ + if (r <= 0) return 0.0; + double s = 0.0; +#pragma unroll + for (int off = -4; off <= 4; ++off) { + const double c = fd_second_coef(r, off); + if (c != 0.0) { + s += c * fd_fetch_axis(src, iF, jF, kF, axis, off, SoA0, SoA1, SoA2); + } + } + return s / (fd_second_denom(r) * h * h); +} + +__device__ __forceinline__ double fd_mixed_axis_radius(const double *src, + int iF, int jF, int kF, + int axis_a, int r_a, double h_a, + int axis_b, int r_b, double h_b, + int SoA0, int SoA1, int SoA2) +{ + if (r_a <= 0 || r_b <= 0) return 0.0; + double s = 0.0; +#pragma unroll + for (int off_a = -4; off_a <= 4; ++off_a) { + const double ca = fd_first_coef(r_a, off_a); + if (ca == 0.0) continue; +#pragma unroll + for (int off_b = -4; off_b <= 4; ++off_b) { + const double cb = fd_first_coef(r_b, off_b); + if (cb != 0.0) { + s += ca * cb * fd_fetch_axis2(src, iF, jF, kF, axis_a, off_a, + axis_b, off_b, SoA0, SoA1, SoA2); + } + } + } + return s / (fd_first_denom(r_a) * fd_first_denom(r_b) * h_a * h_b); +} + +__device__ __forceinline__ void fd_compute_first3(const double *src, + int iF, int jF, int kF, + int iminF, int jminF, int kminF, + int imaxF, int jmaxF, int kmaxF, + int SoA0, int SoA1, int SoA2, + double &fx, double &fy, double &fz) +{ +#if ghost_width == 3 + const int r = fd_common_radius(iF, jF, kF, iminF, jminF, kminF, imaxF, jmaxF, kmaxF); + fx = fd_first_axis_radius(src, iF, jF, kF, 0, r, d_gp.dX, SoA0, SoA1, SoA2); + fy = fd_first_axis_radius(src, iF, jF, kF, 1, r, d_gp.dY, SoA0, SoA1, SoA2); + fz = fd_first_axis_radius(src, iF, jF, kF, 2, r, d_gp.dZ, SoA0, SoA1, SoA2); +#else + fx = fd_first_axis_radius(src, iF, jF, kF, 0, fd_axis_radius(iF, iminF, imaxF), + d_gp.dX, SoA0, SoA1, SoA2); + fy = fd_first_axis_radius(src, iF, jF, kF, 1, fd_axis_radius(jF, jminF, jmaxF), + d_gp.dY, SoA0, SoA1, SoA2); + fz = fd_first_axis_radius(src, iF, jF, kF, 2, fd_axis_radius(kF, kminF, kmaxF), + d_gp.dZ, SoA0, SoA1, SoA2); +#endif +} + +__device__ __forceinline__ void fd_compute_second6(const double *src, + int iF, int jF, int kF, + int iminF, int jminF, int kminF, + int imaxF, int jmaxF, int kmaxF, + int SoA0, int SoA1, int SoA2, + double &fxx, double &fxy, double &fxz, + double &fyy, double &fyz, double &fzz) +{ +#if ghost_width == 3 + const int r = fd_common_radius(iF, jF, kF, iminF, jminF, kminF, imaxF, jmaxF, kmaxF); + const int rx = r, ry = r, rz = r; +#else + const int rx = fd_axis_radius(iF, iminF, imaxF); + const int ry = fd_axis_radius(jF, jminF, jmaxF); + const int rz = fd_axis_radius(kF, kminF, kmaxF); +#endif + fxx = fd_second_axis_radius(src, iF, jF, kF, 0, rx, d_gp.dX, SoA0, SoA1, SoA2); + fyy = fd_second_axis_radius(src, iF, jF, kF, 1, ry, d_gp.dY, SoA0, SoA1, SoA2); + fzz = fd_second_axis_radius(src, iF, jF, kF, 2, rz, d_gp.dZ, SoA0, SoA1, SoA2); + fxy = fd_mixed_axis_radius(src, iF, jF, kF, 0, rx, d_gp.dX, 1, ry, d_gp.dY, SoA0, SoA1, SoA2); + fxz = fd_mixed_axis_radius(src, iF, jF, kF, 0, rx, d_gp.dX, 2, rz, d_gp.dZ, SoA0, SoA1, SoA2); + fyz = fd_mixed_axis_radius(src, iF, jF, kF, 1, ry, d_gp.dY, 2, rz, d_gp.dZ, SoA0, SoA1, SoA2); +} + +__device__ __forceinline__ bool fd_lop_fits(int qF, int qminF, int qmaxF, + int dir, int lo, int hi) +{ + for (int off = lo; off <= hi; ++off) { + const int q = qF + dir * off; + if (q < qminF || q > qmaxF) return false; + } + return true; +} + +__device__ __forceinline__ double fd_lop_fetch_sum(const double *src, + int iF, int jF, int kF, + int axis, int dir, + const double *coef, + int lo, int hi, + int SoA0, int SoA1, int SoA2) +{ + double s = 0.0; + for (int off = lo; off <= hi; ++off) { + const double c = coef[off - lo]; + if (c != 0.0) { + s += c * fd_fetch_axis(src, iF, jF, kF, axis, dir * off, SoA0, SoA1, SoA2); + } + } + return s; +} + +__device__ __forceinline__ double fd_lopsided_axis(const double *src, + int iF, int jF, int kF, + int axis, double speed, + int qF, int qminF, int qmaxF, + double h, + int SoA0, int SoA1, int SoA2) +{ + if (speed == 0.0) return 0.0; + const int dir = (speed > 0.0) ? 1 : -1; + const double mag = (speed > 0.0) ? speed : -speed; + +#if ghost_width == 2 + if (fd_lop_fits(qF, qminF, qmaxF, dir, 0, 2)) { + const double c[] = {-3.0, 4.0, -1.0}; + return mag * fd_lop_fetch_sum(src, iF, jF, kF, axis, dir, c, 0, 2, SoA0, SoA1, SoA2) / (2.0 * h); + } + if (fd_lop_fits(qF, qminF, qmaxF, dir, 0, 1)) { + const double c[] = {-1.0, 1.0}; + return mag * fd_lop_fetch_sum(src, iF, jF, kF, axis, dir, c, 0, 1, SoA0, SoA1, SoA2) / (2.0 * h); + } + return 0.0; +#elif ghost_width == 3 + if (fd_lop_fits(qF, qminF, qmaxF, dir, -1, 3)) { + const double c[] = {-3.0, -10.0, 18.0, -6.0, 1.0}; + return mag * fd_lop_fetch_sum(src, iF, jF, kF, axis, dir, c, -1, 3, SoA0, SoA1, SoA2) / (12.0 * h); + } + const int r = fd_axis_radius(qF, qminF, qmaxF); + return speed * fd_first_axis_radius(src, iF, jF, kF, axis, r, h, SoA0, SoA1, SoA2); +#elif ghost_width == 4 + if (fd_lop_fits(qF, qminF, qmaxF, dir, -2, 4)) { + const double c[] = {2.0, -24.0, -35.0, 80.0, -30.0, 8.0, -1.0}; + return mag * fd_lop_fetch_sum(src, iF, jF, kF, axis, dir, c, -2, 4, SoA0, SoA1, SoA2) / (60.0 * h); + } + if (fd_lop_fits(qF, qminF, qmaxF, dir, -1, 5)) { + const double c[] = {-10.0, -77.0, 150.0, -100.0, 50.0, -15.0, 2.0}; + return mag * fd_lop_fetch_sum(src, iF, jF, kF, axis, dir, c, -1, 5, SoA0, SoA1, SoA2) / (60.0 * h); + } + const int r = fd_axis_radius(qF, qminF, qmaxF); + return speed * fd_first_axis_radius(src, iF, jF, kF, axis, r, h, SoA0, SoA1, SoA2); +#else + if (fd_lop_fits(qF, qminF, qmaxF, dir, -3, 5)) { + const double c[] = {-5.0, 60.0, -420.0, -378.0, 1050.0, -420.0, 140.0, -30.0, 3.0}; + return mag * fd_lop_fetch_sum(src, iF, jF, kF, axis, dir, c, -3, 5, SoA0, SoA1, SoA2) / (840.0 * h); + } + const int r = fd_axis_radius(qF, qminF, qmaxF); + return speed * fd_first_axis_radius(src, iF, jF, kF, axis, r, h, SoA0, SoA1, SoA2); +#endif +} + +__device__ __forceinline__ double fd_ko_coef(int r, int off) +{ + const int a = off < 0 ? -off : off; + if (r == 2) { + if (a == 0) return 6.0; + if (a == 1) return -4.0; + if (a == 2) return 1.0; + } else if (r == 3) { + if (a == 0) return -20.0; + if (a == 1) return 15.0; + if (a == 2) return -6.0; + if (a == 3) return 1.0; + } else if (r == 4) { + if (a == 0) return 70.0; + if (a == 1) return -56.0; + if (a == 2) return 28.0; + if (a == 3) return -8.0; + if (a == 4) return 1.0; + } else if (r == 5) { + if (a == 0) return -252.0; + if (a == 1) return 210.0; + if (a == 2) return -120.0; + if (a == 3) return 45.0; + if (a == 4) return -10.0; + if (a == 5) return 1.0; + } + return 0.0; +} + +__device__ __forceinline__ double fd_ko_axis(const double *src, + int iF, int jF, int kF, + int axis, int r, + int SoA0, int SoA1, int SoA2) +{ + double s = 0.0; +#pragma unroll + for (int off = -5; off <= 5; ++off) { + if (off < -r || off > r) continue; + const double c = fd_ko_coef(r, off); + if (c != 0.0) { + s += c * fd_fetch_axis(src, iF, jF, kF, axis, off, SoA0, SoA1, SoA2); + } + } + return s; +} + +__device__ __forceinline__ double fd_ko_term(const double *src, + int iF, int jF, int kF, + int iminF, int jminF, int kminF, + int imaxF, int jmaxF, int kmaxF, + double eps_val, + int SoA0, int SoA1, int SoA2) +{ + const int r = AMSS_FD_LK_RADIUS; + if (eps_val <= 0.0) return 0.0; +#if ghost_width >= 4 + if (iF - r <= iminF || iF + r >= imaxF || + jF - r <= jminF || jF + r >= jmaxF || + kF - r <= kminF || kF + r >= kmaxF) { + return 0.0; + } +#else + if (iF - r < iminF || iF + r > imaxF || + jF - r < jminF || jF + r > jmaxF || + kF - r < kminF || kF + r > kmaxF) { + return 0.0; + } +#endif + double cof = 1.0; +#pragma unroll + for (int n = 0; n < 2 * r; ++n) cof *= 2.0; + const double sign = (r & 1) ? 1.0 : -1.0; + const double dx = fd_ko_axis(src, iF, jF, kF, 0, r, SoA0, SoA1, SoA2); + const double dy = fd_ko_axis(src, iF, jF, kF, 1, r, SoA0, SoA1, SoA2); + const double dz = fd_ko_axis(src, iF, jF, kF, 2, r, SoA0, SoA1, SoA2); + return sign * eps_val * (dx / d_gp.dX + dy / d_gp.dY + dz / d_gp.dZ) / cof; +} + +#endif diff --git a/AMSS_NCKU_source/makefile b/AMSS_NCKU_source/makefile index 93cf81c..89df7c9 100644 --- a/AMSS_NCKU_source/makefile +++ b/AMSS_NCKU_source/makefile @@ -63,11 +63,11 @@ endif $(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH) # CUDA rewrite of BSSN RHS (drop-in replacement for bssn_rhs_c + stencil helpers) -bssn_rhs_cuda.o: bssn_rhs_cuda.cu bssn_rhs.h macrodef.h +bssn_rhs_cuda.o: bssn_rhs_cuda.cu bssn_rhs.h macrodef.h fd_cuda_helpers.cuh $(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH) # CUDA rewrite of Z4C Cartesian RHS -z4c_rhs_cuda.o: z4c_rhs_cuda.cu z4c_rhs_cuda.h bssn_rhs.h macrodef.h ricci_gamma.h +z4c_rhs_cuda.o: z4c_rhs_cuda.cu z4c_rhs_cuda.h bssn_rhs.h macrodef.h ricci_gamma.h fd_cuda_helpers.cuh $(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH) # C rewrite of BSSN RHS kernel and helpers diff --git a/AMSS_NCKU_source/z4c_rhs_cuda.cu b/AMSS_NCKU_source/z4c_rhs_cuda.cu index 8addc3e..01b4dd8 100644 --- a/AMSS_NCKU_source/z4c_rhs_cuda.cu +++ b/AMSS_NCKU_source/z4c_rhs_cuda.cu @@ -266,6 +266,8 @@ __device__ __forceinline__ double fetch_sym_ord3_direct(const double *src, + (skF - 1) * d_gp.ex[0] * d_gp.ex[1]]; } +#include "fd_cuda_helpers.cuh" + /* ------------------------------------------------------------------ */ /* GPU buffer management */ /* ------------------------------------------------------------------ */ @@ -1419,45 +1421,10 @@ void kern_fderivs_batched(FDerivTables tables, int field_count) const int jF = j0 + 1; const int kF = k0 + 1; - if ((iF + 2) <= imaxF && (iF - 2) >= iminF && - (jF + 2) <= jmaxF && (jF - 2) >= jminF && - (kF + 2) <= kmaxF && (kF - 2) >= kminF) - { - fx[tid] = d_gp.d12dx * ( - fetch_sym_ord2_direct(src, iF - 2, jF, kF, SoA0, SoA1, SoA2) - - 8.0 * fetch_sym_ord2_direct(src, iF - 1, jF, kF, SoA0, SoA1, SoA2) - + 8.0 * fetch_sym_ord2_direct(src, iF + 1, jF, kF, SoA0, SoA1, SoA2) - - fetch_sym_ord2_direct(src, iF + 2, jF, kF, SoA0, SoA1, SoA2)); - fy[tid] = d_gp.d12dy * ( - fetch_sym_ord2_direct(src, iF, jF - 2, kF, SoA0, SoA1, SoA2) - - 8.0 * fetch_sym_ord2_direct(src, iF, jF - 1, kF, SoA0, SoA1, SoA2) - + 8.0 * fetch_sym_ord2_direct(src, iF, jF + 1, kF, SoA0, SoA1, SoA2) - - fetch_sym_ord2_direct(src, iF, jF + 2, kF, SoA0, SoA1, SoA2)); - fz[tid] = d_gp.d12dz * ( - fetch_sym_ord2_direct(src, iF, jF, kF - 2, SoA0, SoA1, SoA2) - - 8.0 * fetch_sym_ord2_direct(src, iF, jF, kF - 1, SoA0, SoA1, SoA2) - + 8.0 * fetch_sym_ord2_direct(src, iF, jF, kF + 1, SoA0, SoA1, SoA2) - - fetch_sym_ord2_direct(src, iF, jF, kF + 2, SoA0, SoA1, SoA2)); - } - else if ((iF + 1) <= imaxF && (iF - 1) >= iminF && - (jF + 1) <= jmaxF && (jF - 1) >= jminF && - (kF + 1) <= kmaxF && (kF - 1) >= kminF) - { - fx[tid] = d_gp.d2dx * ( - -fetch_sym_ord2_direct(src, iF - 1, jF, kF, SoA0, SoA1, SoA2) - +fetch_sym_ord2_direct(src, iF + 1, jF, kF, SoA0, SoA1, SoA2)); - fy[tid] = d_gp.d2dy * ( - -fetch_sym_ord2_direct(src, iF, jF - 1, kF, SoA0, SoA1, SoA2) - +fetch_sym_ord2_direct(src, iF, jF + 1, kF, SoA0, SoA1, SoA2)); - fz[tid] = d_gp.d2dz * ( - -fetch_sym_ord2_direct(src, iF, jF, kF - 1, SoA0, SoA1, SoA2) - +fetch_sym_ord2_direct(src, iF, jF, kF + 1, SoA0, SoA1, SoA2)); - } - else { - fx[tid] = 0.0; - fy[tid] = 0.0; - fz[tid] = 0.0; - } + fd_compute_first3(src, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + SoA0, SoA1, SoA2, + fx[tid], fy[tid], fz[tid]); } __global__ __launch_bounds__(128, 4) @@ -1497,6 +1464,12 @@ void kern_fdderivs_batched(FDDerivTables tables, int field_count) const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + fd_compute_second6(src, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + SoA0, SoA1, SoA2, + fxx[tid], fxy[tid], fxz[tid], fyy[tid], fyz[tid], fzz[tid]); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -1624,12 +1597,43 @@ void kern_fdderivs_batched(FDDerivTables tables, int field_count) fxx[tid] = 0.0; fxy[tid] = 0.0; fxz[tid] = 0.0; fyy[tid] = 0.0; fyz[tid] = 0.0; fzz[tid] = 0.0; } +#endif } +static void gpu_fderivs_batch(int field_count, + double *const *src_fields, + double *const *fx_fields, + double *const *fy_fields, + double *const *fz_fields, + const int *soa_signs, + int all); +static void gpu_fdderivs_batch(int field_count, + double *const *src_fields, + double *const *fxx_fields, + double *const *fxy_fields, + double *const *fxz_fields, + double *const *fyy_fields, + double *const *fyz_fields, + double *const *fzz_fields, + const int *soa_signs, + int all); +static void gpu_lopsided_kodis_single_batch(double *d_f_adv, double *d_f_ko, double *d_f_rhs, + double *d_Sfx, double *d_Sfy, double *d_Sfz, + double SoA0, double SoA1, double SoA2, + double eps_val, int all); + /* symmetry_bd on GPU for ord=2, then launch fderivs kernel */ static void gpu_fderivs(double *d_f, double *d_fx, double *d_fy, double *d_fz, double SoA0, double SoA1, double SoA2, int all) { +#if ghost_width != 3 + double *src_fields[1] = {d_f}; + double *fx_fields[1] = {d_fx}; + double *fy_fields[1] = {d_fy}; + double *fz_fields[1] = {d_fz}; + const int soa_signs[3] = {(int)SoA0, (int)SoA1, (int)SoA2}; + gpu_fderivs_batch(1, src_fields, fx_fields, fy_fields, fz_fields, soa_signs, all); +#else double *fh = g_buf.d_fh2; const size_t nx = (size_t)g_buf.prev_nx; const size_t ny = (size_t)g_buf.prev_ny; @@ -1638,6 +1642,7 @@ static void gpu_fderivs(double *d_f, double *d_fx, double *d_fy, double *d_fz, kern_symbd_pack_ord2<<>>(d_f, fh, SoA0, SoA1, SoA2); kern_fderivs<<>>(fh, d_fx, d_fy, d_fz); +#endif } /* symmetry_bd on GPU for ord=2, then launch fdderivs kernel */ @@ -1646,6 +1651,18 @@ static void gpu_fdderivs(double *d_f, double *d_fyy, double *d_fyz, double *d_fzz, double SoA0, double SoA1, double SoA2, int all) { +#if ghost_width != 3 + double *src_fields[1] = {d_f}; + double *fxx_fields[1] = {d_fxx}; + double *fxy_fields[1] = {d_fxy}; + double *fxz_fields[1] = {d_fxz}; + double *fyy_fields[1] = {d_fyy}; + double *fyz_fields[1] = {d_fyz}; + double *fzz_fields[1] = {d_fzz}; + const int soa_signs[3] = {(int)SoA0, (int)SoA1, (int)SoA2}; + gpu_fdderivs_batch(1, src_fields, fxx_fields, fxy_fields, fxz_fields, + fyy_fields, fyz_fields, fzz_fields, soa_signs, all); +#else double *fh = g_buf.d_fh2; const size_t nx = (size_t)g_buf.prev_nx; const size_t ny = (size_t)g_buf.prev_ny; @@ -1654,6 +1671,7 @@ static void gpu_fdderivs(double *d_f, kern_symbd_pack_ord2<<>>(d_f, fh, SoA0, SoA1, SoA2); kern_fdderivs<<>>(fh, d_fxx, d_fxy, d_fxz, d_fyy, d_fyz, d_fzz); +#endif } static void gpu_fderivs_batch(int field_count, @@ -1743,6 +1761,12 @@ void kern_phase10_ricci_batched(const double * __restrict__ gupxx, const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + fd_compute_second6(src, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + SoA0, SoA1, SoA2, + fxx, fxy, fxz, fyy, fyz, fzz); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -1866,6 +1890,7 @@ void kern_phase10_ricci_batched(const double * __restrict__ gupxx, - fetch_sym_ord2_direct(src, iF, jF - 1, kF + 1, SoA0, SoA1, SoA2) + fetch_sym_ord2_direct(src, iF, jF + 1, kF + 1, SoA0, SoA1, SoA2)); } +#endif } dst[tid] = gupxx[tid] * fxx + gupyy[tid] * fyy + gupzz[tid] * fzz @@ -1930,6 +1955,16 @@ void kern_phase14_lap_chi_derivs(const double * __restrict__ Lap, const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + fd_compute_second6(Lap, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + 1, 1, 1, + fxx[tid], fxy[tid], fxz[tid], fyy[tid], fyz[tid], fzz[tid]); + fd_compute_first3(chi, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + 1, 1, 1, + chix_out[tid], chiy_out[tid], chiz_out[tid]); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -2083,6 +2118,7 @@ void kern_phase14_lap_chi_derivs(const double * __restrict__ Lap, fyy[tid] = 0.0; fyz[tid] = 0.0; fzz[tid] = 0.0; chix_out[tid] = 0.0; chiy_out[tid] = 0.0; chiz_out[tid] = 0.0; } +#endif } /* Combined ord=3 advection + KO dissipation. @@ -2094,6 +2130,11 @@ static void gpu_lopsided_kodis(double *d_f_adv, double *d_f_ko, double *d_f_rhs, double SoA0, double SoA1, double SoA2, double eps_val, int all) { +#if ghost_width != 3 + gpu_lopsided_kodis_single_batch(d_f_adv, d_f_ko, d_f_rhs, + d_Sfx, d_Sfy, d_Sfz, + SoA0, SoA1, SoA2, eps_val, all); +#else double *fh = g_buf.d_fh3; const size_t nx = (size_t)g_buf.prev_nx; const size_t ny = (size_t)g_buf.prev_ny; @@ -2109,6 +2150,7 @@ static void gpu_lopsided_kodis(double *d_f_adv, double *d_f_ko, double *d_f_rhs, } kern_kodis<<>>(fh, d_f_rhs, eps_val); } +#endif } __global__ __launch_bounds__(128, 4) @@ -2141,6 +2183,24 @@ void kern_lopsided_kodis_batched(const double * __restrict__ Sfx, const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + if (do_lopsided && i0 <= nx - 2 && j0 <= ny - 2 && k0 <= nz - 2) { + const double val = + fd_lopsided_axis(adv_src, iF, jF, kF, 0, Sfx[tid], iF, iminF, imaxF, + d_gp.dX, SoA0, SoA1, SoA2) + + fd_lopsided_axis(adv_src, iF, jF, kF, 1, Sfy[tid], jF, jminF, jmaxF, + d_gp.dY, SoA0, SoA1, SoA2) + + fd_lopsided_axis(adv_src, iF, jF, kF, 2, Sfz[tid], kF, kminF, kmaxF, + d_gp.dZ, SoA0, SoA1, SoA2); + rhs[tid] += val; + } + + if (do_kodis) { + rhs[tid] += fd_ko_term(ko_src, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + eps_val, SoA0, SoA1, SoA2); + } +#else if (do_lopsided && i0 <= nx - 2 && j0 <= ny - 2 && k0 <= nz - 2) { double val = 0.0; @@ -2323,6 +2383,25 @@ void kern_lopsided_kodis_batched(const double * __restrict__ Sfx, rhs[tid] += (eps_val / cof) * (Dx / d_gp.dX + Dy / d_gp.dY + Dz / d_gp.dZ); } +#endif +} + +static void gpu_lopsided_kodis_single_batch(double *d_f_adv, double *d_f_ko, double *d_f_rhs, + double *d_Sfx, double *d_Sfy, double *d_Sfz, + double SoA0, double SoA1, double SoA2, + double eps_val, int all) +{ + LopsidedKodisTables tables = {}; + tables.adv_fields[0] = d_f_adv; + tables.ko_fields[0] = d_f_ko; + tables.rhs_fields[0] = d_f_rhs; + tables.soa_signs[0] = (int)SoA0; + tables.soa_signs[1] = (int)SoA1; + tables.soa_signs[2] = (int)SoA2; + + dim3 launch_grid((unsigned int)grid((size_t)all), 1u); + kern_lopsided_kodis_batched<<>>( + d_Sfx, d_Sfy, d_Sfz, tables, eps_val, 1, eps_val > 0.0 ? 1 : 0); } static void gpu_lopsided_kodis_state_batch(double eps_val, int all) @@ -3873,6 +3952,12 @@ void kern_phase12_13_chi_correction_fused( const int jF = j0 + 1; const int kF = k0 + 1; +#if ghost_width != 3 + fd_compute_second6(chi, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + 1, 1, 1, + cxx, cxy, cxz, cyy, cyz, czz); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -3996,6 +4081,7 @@ void kern_phase12_13_chi_correction_fused( - fetch_sym_ord2_direct(chi, iF, jF - 1, kF + 1, 1, 1, 1) + fetch_sym_ord2_direct(chi, iF, jF + 1, kF + 1, 1, 1, 1)); } +#endif } const double cx = chix[tid]; @@ -4166,6 +4252,12 @@ void kern_phase15_trK_Aij_gauge( double fyy_v = 0.0, fyz_v = 0.0, fzz_v = 0.0; if (!(i0 > nx - 2 || j0 > ny - 2 || k0 > nz - 2)) { +#if ghost_width != 3 + fd_compute_second6(alpn1, iF, jF, kF, + iminF, jminF, kminF, imaxF, jmaxF, kmaxF, + 1, 1, 1, + fxx_v, fxy_v, fxz_v, fyy_v, fyz_v, fzz_v); +#else if ((iF + 2) <= imaxF && (iF - 2) >= iminF && (jF + 2) <= jmaxF && (jF - 2) >= jminF && (kF + 2) <= kmaxF && (kF - 2) >= kminF) @@ -4289,6 +4381,7 @@ void kern_phase15_trK_Aij_gauge( - fetch_sym_ord2_direct(alpn1, iF, jF - 1, kF + 1, 1, 1, 1) + fetch_sym_ord2_direct(alpn1, iF, jF + 1, kF + 1, 1, 1, 1)); } +#endif } /* raised chi/chi */ @@ -4626,15 +4719,15 @@ static void setup_grid_params(int *ex, gp.imaxF = nx; gp.jmaxF = ny; gp.kmaxF = nz; - if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) gp.kminF = -1; - if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) gp.iminF = -1; - if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) gp.jminF = -1; + if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) gp.kminF = 2 - ghost_width; + if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) gp.iminF = 2 - ghost_width; + if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) gp.jminF = 2 - ghost_width; gp.iminF3 = 1; gp.jminF3 = 1; gp.kminF3 = 1; - if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) gp.kminF3 = -2; - if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) gp.iminF3 = -2; - if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) gp.jminF3 = -2; + if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) gp.kminF3 = 1 - ghost_width; + if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) gp.iminF3 = 1 - ghost_width; + if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) gp.jminF3 = 1 - ghost_width; gp.Symmetry = Symmetry; gp.eps = eps; gp.co = co;