Fix C derivative ghost-buffer indexing across FD orders
This commit is contained in:
@@ -93,11 +93,11 @@ void fdderivs(const int ex[3],
|
||||
|
||||
if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
|
||||
for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
|
||||
const int kF = k0;
|
||||
const int kF = k0 + 1;
|
||||
for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
|
||||
const int jF = j0;
|
||||
const int jF = j0 + 1;
|
||||
for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
|
||||
const int iF = i0;
|
||||
const int iF = i0 + 1;
|
||||
const size_t p = idx_ex(i0, j0, k0, ex);
|
||||
|
||||
fxx[p] = Sdxdx * (
|
||||
@@ -471,19 +471,19 @@ void fdderivs(const int ex[3],
|
||||
}
|
||||
|
||||
/* loop bounds for each pass (from widest to narrowest) */
|
||||
const int i2_lo = (iminF + 1 > 0) ? (iminF + 1) : 0;
|
||||
const int j2_lo = (jminF + 1 > 0) ? (jminF + 1) : 0;
|
||||
const int k2_lo = (kminF + 1 > 0) ? (kminF + 1) : 0;
|
||||
const int i2_lo = (iminF > 0) ? iminF : 0;
|
||||
const int j2_lo = (jminF > 0) ? jminF : 0;
|
||||
const int k2_lo = (kminF > 0) ? kminF : 0;
|
||||
const int i2_hi = ex1 - 2, j2_hi = ex2 - 2, k2_hi = ex3 - 2;
|
||||
|
||||
const int i4_lo = (iminF + 2 > 0) ? (iminF + 2) : 0;
|
||||
const int j4_lo = (jminF + 2 > 0) ? (jminF + 2) : 0;
|
||||
const int k4_lo = (kminF + 2 > 0) ? (kminF + 2) : 0;
|
||||
const int i4_lo = (iminF + 1 > 0) ? (iminF + 1) : 0;
|
||||
const int j4_lo = (jminF + 1 > 0) ? (jminF + 1) : 0;
|
||||
const int k4_lo = (kminF + 1 > 0) ? (kminF + 1) : 0;
|
||||
const int i4_hi = ex1 - 3, j4_hi = ex2 - 3, k4_hi = ex3 - 3;
|
||||
|
||||
const int i6_lo = (iminF + 3 > 0) ? (iminF + 3) : 0;
|
||||
const int j6_lo = (jminF + 3 > 0) ? (jminF + 3) : 0;
|
||||
const int k6_lo = (kminF + 3 > 0) ? (kminF + 3) : 0;
|
||||
const int i6_lo = (iminF + 2 > 0) ? (iminF + 2) : 0;
|
||||
const int j6_lo = (jminF + 2 > 0) ? (jminF + 2) : 0;
|
||||
const int k6_lo = (kminF + 2 > 0) ? (kminF + 2) : 0;
|
||||
const int i6_hi = ex1 - 4, j6_hi = ex2 - 4, k6_hi = ex3 - 4;
|
||||
|
||||
const int has4 = (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi);
|
||||
@@ -492,13 +492,13 @@ void fdderivs(const int ex[3],
|
||||
/* 2nd-order: skip 4th+6th overlap */
|
||||
if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
|
||||
for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
|
||||
const int kF = k0 + 2;
|
||||
const int kF = k0 + 1;
|
||||
for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
|
||||
const int jF = j0 + 2;
|
||||
const int jF = j0 + 1;
|
||||
for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
|
||||
_Bool in4 = has4 && i0>=i4_lo && i0<=i4_hi && j0>=j4_lo && j0<=j4_hi && k0>=k4_lo && k0<=k4_hi;
|
||||
if (in4) continue;
|
||||
const int iF = i0 + 2;
|
||||
const int iF = i0 + 1;
|
||||
const size_t p = idx_ex(i0, j0, k0, ex);
|
||||
|
||||
fxx[p] = Sdxdx * (fh[idx_fh_F(iF - 1, jF, kF, ex)] - TWO*fh[idx_fh_F(iF,jF,kF,ex)] + fh[idx_fh_F(iF + 1, jF, kF, ex)]);
|
||||
@@ -516,12 +516,12 @@ void fdderivs(const int ex[3],
|
||||
/* 4th-order: skip 6th overlap */
|
||||
if (has4) {
|
||||
for (int k0 = k4_lo; k0 <= k4_hi; ++k0) {
|
||||
const int kF = k0 + 2;
|
||||
const int kF = k0 + 1;
|
||||
for (int j0 = j4_lo; j0 <= j4_hi; ++j0) {
|
||||
const int jF = j0 + 2;
|
||||
const int jF = j0 + 1;
|
||||
for (int i0 = i4_lo; i0 <= i4_hi; ++i0) {
|
||||
if (has6 && i0>=i6_lo && i0<=i6_hi && j0>=j6_lo && j0<=j6_hi && k0>=k6_lo && k0<=k6_hi) continue;
|
||||
const int iF = i0 + 2;
|
||||
const int iF = i0 + 1;
|
||||
const size_t p = idx_ex(i0, j0, k0, ex);
|
||||
|
||||
fxx[p] = Fdxdx * (-fh[idx_fh_F(iF - 2, jF, kF, ex)] + F16*fh[idx_fh_F(iF-1,jF,kF,ex)] - F30*fh[idx_fh_F(iF,jF,kF,ex)] - fh[idx_fh_F(iF+2,jF,kF,ex)] + F16*fh[idx_fh_F(iF+1,jF,kF,ex)]);
|
||||
@@ -557,11 +557,11 @@ void fdderivs(const int ex[3],
|
||||
/* 6th-order: interior only */
|
||||
if (has6) {
|
||||
for (int k0 = k6_lo; k0 <= k6_hi; ++k0) {
|
||||
const int kF = k0 + 2;
|
||||
const int kF = k0 + 1;
|
||||
for (int j0 = j6_lo; j0 <= j6_hi; ++j0) {
|
||||
const int jF = j0 + 2;
|
||||
const int jF = j0 + 1;
|
||||
for (int i0 = i6_lo; i0 <= i6_hi; ++i0) {
|
||||
const int iF = i0 + 2;
|
||||
const int iF = i0 + 1;
|
||||
const size_t p = idx_ex(i0, j0, k0, ex);
|
||||
|
||||
/* Diagonal: [+2,-27,+270,-490,+270,-27,+2] / (180*dx^2) */
|
||||
@@ -685,24 +685,24 @@ void fdderivs(const int ex[3],
|
||||
}
|
||||
|
||||
/* Loop bounds for each pass */
|
||||
const int i2_lo = (iminF + 2 > 0) ? (iminF + 2) : 0;
|
||||
const int j2_lo = (jminF + 2 > 0) ? (jminF + 2) : 0;
|
||||
const int k2_lo = (kminF + 2 > 0) ? (kminF + 2) : 0;
|
||||
const int i2_lo = (iminF > 0) ? iminF : 0;
|
||||
const int j2_lo = (jminF > 0) ? jminF : 0;
|
||||
const int k2_lo = (kminF > 0) ? kminF : 0;
|
||||
const int i2_hi = ex1 - 2, j2_hi = ex2 - 2, k2_hi = ex3 - 2;
|
||||
|
||||
const int i4_lo = (iminF + 3 > 0) ? (iminF + 3) : 0;
|
||||
const int j4_lo = (jminF + 3 > 0) ? (jminF + 3) : 0;
|
||||
const int k4_lo = (kminF + 3 > 0) ? (kminF + 3) : 0;
|
||||
const int i4_lo = (iminF + 1 > 0) ? (iminF + 1) : 0;
|
||||
const int j4_lo = (jminF + 1 > 0) ? (jminF + 1) : 0;
|
||||
const int k4_lo = (kminF + 1 > 0) ? (kminF + 1) : 0;
|
||||
const int i4_hi = ex1 - 3, j4_hi = ex2 - 3, k4_hi = ex3 - 3;
|
||||
|
||||
const int i6_lo = (iminF + 4 > 0) ? (iminF + 4) : 0;
|
||||
const int j6_lo = (jminF + 4 > 0) ? (jminF + 4) : 0;
|
||||
const int k6_lo = (kminF + 4 > 0) ? (kminF + 4) : 0;
|
||||
const int i6_lo = (iminF + 2 > 0) ? (iminF + 2) : 0;
|
||||
const int j6_lo = (jminF + 2 > 0) ? (jminF + 2) : 0;
|
||||
const int k6_lo = (kminF + 2 > 0) ? (kminF + 2) : 0;
|
||||
const int i6_hi = ex1 - 4, j6_hi = ex2 - 4, k6_hi = ex3 - 4;
|
||||
|
||||
const int i8_lo = (iminF + 5 > 0) ? (iminF + 5) : 0;
|
||||
const int j8_lo = (jminF + 5 > 0) ? (jminF + 5) : 0;
|
||||
const int k8_lo = (kminF + 5 > 0) ? (kminF + 5) : 0;
|
||||
const int i8_lo = (iminF + 3 > 0) ? (iminF + 3) : 0;
|
||||
const int j8_lo = (jminF + 3 > 0) ? (jminF + 3) : 0;
|
||||
const int k8_lo = (kminF + 3 > 0) ? (kminF + 3) : 0;
|
||||
const int i8_hi = ex1 - 5, j8_hi = ex2 - 5, k8_hi = ex3 - 5;
|
||||
|
||||
const int has4 = (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi);
|
||||
@@ -712,13 +712,13 @@ void fdderivs(const int ex[3],
|
||||
/* 2nd-order: skip 4th+6th+8th overlap */
|
||||
if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
|
||||
for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
|
||||
const int kF = k0 + 3;
|
||||
const int kF = k0 + 1;
|
||||
for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
|
||||
const int jF = j0 + 3;
|
||||
const int jF = j0 + 1;
|
||||
for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
|
||||
_Bool in4 = has4 && i0>=i4_lo && i0<=i4_hi && j0>=j4_lo && j0<=j4_hi && k0>=k4_lo && k0<=k4_hi;
|
||||
if (in4) continue;
|
||||
const int iF = i0 + 3;
|
||||
const int iF = i0 + 1;
|
||||
const size_t p = idx_ex(i0, j0, k0, ex);
|
||||
|
||||
fxx[p] = Sdxdx * (fh[idx_fh_F_ord4(iF-1,jF,kF,ex)] - TWO*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + fh[idx_fh_F_ord4(iF+1,jF,kF,ex)]);
|
||||
@@ -736,13 +736,13 @@ void fdderivs(const int ex[3],
|
||||
/* 4th-order: skip 6th+8th overlap */
|
||||
if (has4) {
|
||||
for (int k0 = k4_lo; k0 <= k4_hi; ++k0) {
|
||||
const int kF = k0 + 3;
|
||||
const int kF = k0 + 1;
|
||||
for (int j0 = j4_lo; j0 <= j4_hi; ++j0) {
|
||||
const int jF = j0 + 3;
|
||||
const int jF = j0 + 1;
|
||||
for (int i0 = i4_lo; i0 <= i4_hi; ++i0) {
|
||||
_Bool in6 = has6 && i0>=i6_lo && i0<=i6_hi && j0>=j6_lo && j0<=j6_hi && k0>=k6_lo && k0<=k6_hi;
|
||||
if (in6) continue;
|
||||
const int iF = i0 + 3;
|
||||
const int iF = i0 + 1;
|
||||
const size_t p = idx_ex(i0, j0, k0, ex);
|
||||
|
||||
fxx[p] = Fdxdx * (-fh[idx_fh_F_ord4(iF-2,jF,kF,ex)] + F16*fh[idx_fh_F_ord4(iF-1,jF,kF,ex)] - F30*fh[idx_fh_F_ord4(iF,jF,kF,ex)] - fh[idx_fh_F_ord4(iF+2,jF,kF,ex)] + F16*fh[idx_fh_F_ord4(iF+1,jF,kF,ex)]);
|
||||
@@ -778,12 +778,12 @@ void fdderivs(const int ex[3],
|
||||
/* 6th-order: skip 8th overlap */
|
||||
if (has6) {
|
||||
for (int k0 = k6_lo; k0 <= k6_hi; ++k0) {
|
||||
const int kF = k0 + 3;
|
||||
const int kF = k0 + 1;
|
||||
for (int j0 = j6_lo; j0 <= j6_hi; ++j0) {
|
||||
const int jF = j0 + 3;
|
||||
const int jF = j0 + 1;
|
||||
for (int i0 = i6_lo; i0 <= i6_hi; ++i0) {
|
||||
if (has8 && i0>=i8_lo && i0<=i8_hi && j0>=j8_lo && j0<=j8_hi && k0>=k8_lo && k0<=k8_hi) continue;
|
||||
const int iF = i0 + 3;
|
||||
const int iF = i0 + 1;
|
||||
const size_t p = idx_ex(i0, j0, k0, ex);
|
||||
|
||||
fxx[p] = Xdxdx * (
|
||||
@@ -817,11 +817,11 @@ void fdderivs(const int ex[3],
|
||||
/* 8th-order: interior only */
|
||||
if (has8) {
|
||||
for (int k0 = k8_lo; k0 <= k8_hi; ++k0) {
|
||||
const int kF = k0 + 3;
|
||||
const int kF = k0 + 1;
|
||||
for (int j0 = j8_lo; j0 <= j8_hi; ++j0) {
|
||||
const int jF = j0 + 3;
|
||||
const int jF = j0 + 1;
|
||||
for (int i0 = i8_lo; i0 <= i8_hi; ++i0) {
|
||||
const int iF = i0 + 3;
|
||||
const int iF = i0 + 1;
|
||||
const size_t p = idx_ex(i0, j0, k0, ex);
|
||||
|
||||
/* Diagonal: [-9,+128,-1008,+8064,-14350,+8064,-1008,+128,-9] / (5040*dx^2) */
|
||||
|
||||
Reference in New Issue
Block a user