413 lines
16 KiB
Plaintext
413 lines
16 KiB
Plaintext
#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
|