64-BitBrainstorm_2026/AMSS-NCKU
8
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fix eighth-order C derivative and lopsided stencils

Browse Source
This commit is contained in:
CGH0S7
2026-05-14 20:32:53 +08:00
parent 3b8774c1b1
commit 2bbde059db
5 changed files with 140 additions and 222 deletions
Download Patch File Download Diff File Expand all files Collapse all files

118
AMSS_NCKU_source/fdderivs_c.C
View File

@@ -496,7 +496,7 @@ void fdderivs(const int ex[3],
for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
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;
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 + 1;
const size_t p = idx_ex(i0, j0, k0, ex);
@@ -628,7 +628,7 @@ void fdderivs(const int ex[3],
#elif (ghost_width == 5)
/* ---- 8th-order ----------------------------------------------------- */
{
const int ord = 4;
const int ord = 5;
int iminF = 1, jminF = 1, kminF = 1;
if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -3;
@@ -716,18 +716,18 @@ void fdderivs(const int ex[3],
for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
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;
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 + 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)]);
fyy[p] = Sdydy * (fh[idx_fh_F_ord4(iF,jF-1,kF,ex)] - TWO*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + fh[idx_fh_F_ord4(iF,jF+1,kF,ex)]);
fzz[p] = Sdzdz * (fh[idx_fh_F_ord4(iF,jF,kF-1,ex)] - TWO*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + fh[idx_fh_F_ord4(iF,jF,kF+1,ex)]);
fxx[p] = Sdxdx * (fh[idx_fh_F_ord5(iF-1,jF,kF,ex)] - TWO*fh[idx_fh_F_ord5(iF,jF,kF,ex)] + fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]);
fyy[p] = Sdydy * (fh[idx_fh_F_ord5(iF,jF-1,kF,ex)] - TWO*fh[idx_fh_F_ord5(iF,jF,kF,ex)] + fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]);
fzz[p] = Sdzdz * (fh[idx_fh_F_ord5(iF,jF,kF-1,ex)] - TWO*fh[idx_fh_F_ord5(iF,jF,kF,ex)] + fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]);
fxy[p] = Sdxdy * (fh[idx_fh_F_ord4(iF-1,jF-1,kF,ex)] - fh[idx_fh_F_ord4(iF+1,jF-1,kF,ex)] - fh[idx_fh_F_ord4(iF-1,jF+1,kF,ex)] + fh[idx_fh_F_ord4(iF+1,jF+1,kF,ex)]);
fxz[p] = Sdxdz * (fh[idx_fh_F_ord4(iF-1,jF,kF-1,ex)] - fh[idx_fh_F_ord4(iF+1,jF,kF-1,ex)] - fh[idx_fh_F_ord4(iF-1,jF,kF+1,ex)] + fh[idx_fh_F_ord4(iF+1,jF,kF+1,ex)]);
fyz[p] = Sdydz * (fh[idx_fh_F_ord4(iF,jF-1,kF-1,ex)] - fh[idx_fh_F_ord4(iF,jF+1,kF-1,ex)] - fh[idx_fh_F_ord4(iF,jF-1,kF+1,ex)] + fh[idx_fh_F_ord4(iF,jF+1,kF+1,ex)]);
fxy[p] = Sdxdy * (fh[idx_fh_F_ord5(iF-1,jF-1,kF,ex)] - fh[idx_fh_F_ord5(iF+1,jF-1,kF,ex)] - fh[idx_fh_F_ord5(iF-1,jF+1,kF,ex)] + fh[idx_fh_F_ord5(iF+1,jF+1,kF,ex)]);
fxz[p] = Sdxdz * (fh[idx_fh_F_ord5(iF-1,jF,kF-1,ex)] - fh[idx_fh_F_ord5(iF+1,jF,kF-1,ex)] - fh[idx_fh_F_ord5(iF-1,jF,kF+1,ex)] + fh[idx_fh_F_ord5(iF+1,jF,kF+1,ex)]);
fyz[p] = Sdydz * (fh[idx_fh_F_ord5(iF,jF-1,kF-1,ex)] - fh[idx_fh_F_ord5(iF,jF+1,kF-1,ex)] - fh[idx_fh_F_ord5(iF,jF-1,kF+1,ex)] + fh[idx_fh_F_ord5(iF,jF+1,kF+1,ex)]);
}
}
}
@@ -740,34 +740,34 @@ void fdderivs(const int ex[3],
for (int j0 = j4_lo; j0 <= j4_hi; ++j0) {
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;
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 + 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)]);
fyy[p] = Fdydy * (-fh[idx_fh_F_ord4(iF,jF-2,kF,ex)] + F16*fh[idx_fh_F_ord4(iF,jF-1,kF,ex)] - F30*fh[idx_fh_F_ord4(iF,jF,kF,ex)] - fh[idx_fh_F_ord4(iF,jF+2,kF,ex)] + F16*fh[idx_fh_F_ord4(iF,jF+1,kF,ex)]);
fzz[p] = Fdzdz * (-fh[idx_fh_F_ord4(iF,jF,kF-2,ex)] + F16*fh[idx_fh_F_ord4(iF,jF,kF-1,ex)] - F30*fh[idx_fh_F_ord4(iF,jF,kF,ex)] - fh[idx_fh_F_ord4(iF,jF,kF+2,ex)] + F16*fh[idx_fh_F_ord4(iF,jF,kF+1,ex)]);
fxx[p] = Fdxdx * (-fh[idx_fh_F_ord5(iF-2,jF,kF,ex)] + F16*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)] - F30*fh[idx_fh_F_ord5(iF,jF,kF,ex)] - fh[idx_fh_F_ord5(iF+2,jF,kF,ex)] + F16*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]);
fyy[p] = Fdydy * (-fh[idx_fh_F_ord5(iF,jF-2,kF,ex)] + F16*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)] - F30*fh[idx_fh_F_ord5(iF,jF,kF,ex)] - fh[idx_fh_F_ord5(iF,jF+2,kF,ex)] + F16*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]);
fzz[p] = Fdzdz * (-fh[idx_fh_F_ord5(iF,jF,kF-2,ex)] + F16*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)] - F30*fh[idx_fh_F_ord5(iF,jF,kF,ex)] - fh[idx_fh_F_ord5(iF,jF,kF+2,ex)] + F16*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]);
{
const double t_jm2 = (fh[idx_fh_F_ord4(iF-2,jF-2,kF,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF-2,kF,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF-2,kF,ex)]-fh[idx_fh_F_ord4(iF+2,jF-2,kF,ex)]);
const double t_jm1 = (fh[idx_fh_F_ord4(iF-2,jF-1,kF,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF-1,kF,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF-1,kF,ex)]-fh[idx_fh_F_ord4(iF+2,jF-1,kF,ex)]);
const double t_jp1 = (fh[idx_fh_F_ord4(iF-2,jF+1,kF,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF+1,kF,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF+1,kF,ex)]-fh[idx_fh_F_ord4(iF+2,jF+1,kF,ex)]);
const double t_jp2 = (fh[idx_fh_F_ord4(iF-2,jF+2,kF,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF+2,kF,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF+2,kF,ex)]-fh[idx_fh_F_ord4(iF+2,jF+2,kF,ex)]);
const double t_jm2 = (fh[idx_fh_F_ord5(iF-2,jF-2,kF,ex)]-F8*fh[idx_fh_F_ord5(iF-1,jF-2,kF,ex)]+F8*fh[idx_fh_F_ord5(iF+1,jF-2,kF,ex)]-fh[idx_fh_F_ord5(iF+2,jF-2,kF,ex)]);
const double t_jm1 = (fh[idx_fh_F_ord5(iF-2,jF-1,kF,ex)]-F8*fh[idx_fh_F_ord5(iF-1,jF-1,kF,ex)]+F8*fh[idx_fh_F_ord5(iF+1,jF-1,kF,ex)]-fh[idx_fh_F_ord5(iF+2,jF-1,kF,ex)]);
const double t_jp1 = (fh[idx_fh_F_ord5(iF-2,jF+1,kF,ex)]-F8*fh[idx_fh_F_ord5(iF-1,jF+1,kF,ex)]+F8*fh[idx_fh_F_ord5(iF+1,jF+1,kF,ex)]-fh[idx_fh_F_ord5(iF+2,jF+1,kF,ex)]);
const double t_jp2 = (fh[idx_fh_F_ord5(iF-2,jF+2,kF,ex)]-F8*fh[idx_fh_F_ord5(iF-1,jF+2,kF,ex)]+F8*fh[idx_fh_F_ord5(iF+1,jF+2,kF,ex)]-fh[idx_fh_F_ord5(iF+2,jF+2,kF,ex)]);
fxy[p] = Fdxdy * (t_jm2 - F8*t_jm1 + F8*t_jp1 - t_jp2);
}
{
const double t_km2 = (fh[idx_fh_F_ord4(iF-2,jF,kF-2,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF,kF-2,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF,kF-2,ex)]-fh[idx_fh_F_ord4(iF+2,jF,kF-2,ex)]);
const double t_km1 = (fh[idx_fh_F_ord4(iF-2,jF,kF-1,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF,kF-1,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF,kF-1,ex)]-fh[idx_fh_F_ord4(iF+2,jF,kF-1,ex)]);
const double t_kp1 = (fh[idx_fh_F_ord4(iF-2,jF,kF+1,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF,kF+1,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF,kF+1,ex)]-fh[idx_fh_F_ord4(iF+2,jF,kF+1,ex)]);
const double t_kp2 = (fh[idx_fh_F_ord4(iF-2,jF,kF+2,ex)]-F8*fh[idx_fh_F_ord4(iF-1,jF,kF+2,ex)]+F8*fh[idx_fh_F_ord4(iF+1,jF,kF+2,ex)]-fh[idx_fh_F_ord4(iF+2,jF,kF+2,ex)]);
const double t_km2 = (fh[idx_fh_F_ord5(iF-2,jF,kF-2,ex)]-F8*fh[idx_fh_F_ord5(iF-1,jF,kF-2,ex)]+F8*fh[idx_fh_F_ord5(iF+1,jF,kF-2,ex)]-fh[idx_fh_F_ord5(iF+2,jF,kF-2,ex)]);
const double t_km1 = (fh[idx_fh_F_ord5(iF-2,jF,kF-1,ex)]-F8*fh[idx_fh_F_ord5(iF-1,jF,kF-1,ex)]+F8*fh[idx_fh_F_ord5(iF+1,jF,kF-1,ex)]-fh[idx_fh_F_ord5(iF+2,jF,kF-1,ex)]);
const double t_kp1 = (fh[idx_fh_F_ord5(iF-2,jF,kF+1,ex)]-F8*fh[idx_fh_F_ord5(iF-1,jF,kF+1,ex)]+F8*fh[idx_fh_F_ord5(iF+1,jF,kF+1,ex)]-fh[idx_fh_F_ord5(iF+2,jF,kF+1,ex)]);
const double t_kp2 = (fh[idx_fh_F_ord5(iF-2,jF,kF+2,ex)]-F8*fh[idx_fh_F_ord5(iF-1,jF,kF+2,ex)]+F8*fh[idx_fh_F_ord5(iF+1,jF,kF+2,ex)]-fh[idx_fh_F_ord5(iF+2,jF,kF+2,ex)]);
fxz[p] = Fdxdz * (t_km2 - F8*t_km1 + F8*t_kp1 - t_kp2);
}
{
const double t_km2 = (fh[idx_fh_F_ord4(iF,jF-2,kF-2,ex)]-F8*fh[idx_fh_F_ord4(iF,jF-1,kF-2,ex)]+F8*fh[idx_fh_F_ord4(iF,jF+1,kF-2,ex)]-fh[idx_fh_F_ord4(iF,jF+2,kF-2,ex)]);
const double t_km1 = (fh[idx_fh_F_ord4(iF,jF-2,kF-1,ex)]-F8*fh[idx_fh_F_ord4(iF,jF-1,kF-1,ex)]+F8*fh[idx_fh_F_ord4(iF,jF+1,kF-1,ex)]-fh[idx_fh_F_ord4(iF,jF+2,kF-1,ex)]);
const double t_kp1 = (fh[idx_fh_F_ord4(iF,jF-2,kF+1,ex)]-F8*fh[idx_fh_F_ord4(iF,jF-1,kF+1,ex)]+F8*fh[idx_fh_F_ord4(iF,jF+1,kF+1,ex)]-fh[idx_fh_F_ord4(iF,jF+2,kF+1,ex)]);
const double t_kp2 = (fh[idx_fh_F_ord4(iF,jF-2,kF+2,ex)]-F8*fh[idx_fh_F_ord4(iF,jF-1,kF+2,ex)]+F8*fh[idx_fh_F_ord4(iF,jF+1,kF+2,ex)]-fh[idx_fh_F_ord4(iF,jF+2,kF+2,ex)]);
const double t_km2 = (fh[idx_fh_F_ord5(iF,jF-2,kF-2,ex)]-F8*fh[idx_fh_F_ord5(iF,jF-1,kF-2,ex)]+F8*fh[idx_fh_F_ord5(iF,jF+1,kF-2,ex)]-fh[idx_fh_F_ord5(iF,jF+2,kF-2,ex)]);
const double t_km1 = (fh[idx_fh_F_ord5(iF,jF-2,kF-1,ex)]-F8*fh[idx_fh_F_ord5(iF,jF-1,kF-1,ex)]+F8*fh[idx_fh_F_ord5(iF,jF+1,kF-1,ex)]-fh[idx_fh_F_ord5(iF,jF+2,kF-1,ex)]);
const double t_kp1 = (fh[idx_fh_F_ord5(iF,jF-2,kF+1,ex)]-F8*fh[idx_fh_F_ord5(iF,jF-1,kF+1,ex)]+F8*fh[idx_fh_F_ord5(iF,jF+1,kF+1,ex)]-fh[idx_fh_F_ord5(iF,jF+2,kF+1,ex)]);
const double t_kp2 = (fh[idx_fh_F_ord5(iF,jF-2,kF+2,ex)]-F8*fh[idx_fh_F_ord5(iF,jF-1,kF+2,ex)]+F8*fh[idx_fh_F_ord5(iF,jF+1,kF+2,ex)]-fh[idx_fh_F_ord5(iF,jF+2,kF+2,ex)]);
fyz[p] = Fdydz * (t_km2 - F8*t_km1 + F8*t_kp1 - t_kp2);
}
}
@@ -787,15 +787,15 @@ void fdderivs(const int ex[3],
const size_t p = idx_ex(i0, j0, k0, ex);
fxx[p] = Xdxdx * (
TWO * fh[idx_fh_F_ord4(iF-3,jF,kF,ex)] - F27*fh[idx_fh_F_ord4(iF-2,jF,kF,ex)] + F270*fh[idx_fh_F_ord4(iF-1,jF,kF,ex)] - F490*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord4(iF+1,jF,kF,ex)] - F27*fh[idx_fh_F_ord4(iF+2,jF,kF,ex)] + TWO*fh[idx_fh_F_ord4(iF+3,jF,kF,ex)]);
TWO * fh[idx_fh_F_ord5(iF-3,jF,kF,ex)] - F27*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)] + F270*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)] - F490*fh[idx_fh_F_ord5(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)] - F27*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)] + TWO*fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]);
fyy[p] = Xdydy * (
TWO * fh[idx_fh_F_ord4(iF,jF-3,kF,ex)] - F27*fh[idx_fh_F_ord4(iF,jF-2,kF,ex)] + F270*fh[idx_fh_F_ord4(iF,jF-1,kF,ex)] - F490*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord4(iF,jF+1,kF,ex)] - F27*fh[idx_fh_F_ord4(iF,jF+2,kF,ex)] + TWO*fh[idx_fh_F_ord4(iF,jF+3,kF,ex)]);
TWO * fh[idx_fh_F_ord5(iF,jF-3,kF,ex)] - F27*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)] + F270*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)] - F490*fh[idx_fh_F_ord5(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)] - F27*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)] + TWO*fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]);
fzz[p] = Xdzdz * (
TWO * fh[idx_fh_F_ord4(iF,jF,kF-3,ex)] - F27*fh[idx_fh_F_ord4(iF,jF,kF-2,ex)] + F270*fh[idx_fh_F_ord4(iF,jF,kF-1,ex)] - F490*fh[idx_fh_F_ord4(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord4(iF,jF,kF+1,ex)] - F27*fh[idx_fh_F_ord4(iF,jF,kF+2,ex)] + TWO*fh[idx_fh_F_ord4(iF,jF,kF+3,ex)]);
TWO * fh[idx_fh_F_ord5(iF,jF,kF-3,ex)] - F27*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)] + F270*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)] - F490*fh[idx_fh_F_ord5(iF,jF,kF,ex)] + F270*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)] - F27*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)] + TWO*fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]);
{
#define XSTEN6_8(JF, KF_DUMMY) \
(-fh[idx_fh_F_ord4(iF-3,JF,KF_DUMMY,ex)] + F9*fh[idx_fh_F_ord4(iF-2,JF,KF_DUMMY,ex)] - F45*fh[idx_fh_F_ord4(iF-1,JF,KF_DUMMY,ex)] + F45*fh[idx_fh_F_ord4(iF+1,JF,KF_DUMMY,ex)] - F9*fh[idx_fh_F_ord4(iF+2,JF,KF_DUMMY,ex)] + fh[idx_fh_F_ord4(iF+3,JF,KF_DUMMY,ex)])
(-fh[idx_fh_F_ord5(iF-3,JF,KF_DUMMY,ex)] + F9*fh[idx_fh_F_ord5(iF-2,JF,KF_DUMMY,ex)] - F45*fh[idx_fh_F_ord5(iF-1,JF,KF_DUMMY,ex)] + F45*fh[idx_fh_F_ord5(iF+1,JF,KF_DUMMY,ex)] - F9*fh[idx_fh_F_ord5(iF+2,JF,KF_DUMMY,ex)] + fh[idx_fh_F_ord5(iF+3,JF,KF_DUMMY,ex)])
fxy[p] = Xdxdy * (
-XSTEN6_8(jF-3,kF) + F9*XSTEN6_8(jF-2,kF) - F45*XSTEN6_8(jF-1,kF) + F45*XSTEN6_8(jF+1,kF) - F9*XSTEN6_8(jF+2,kF) + XSTEN6_8(jF+3,kF));
fxz[p] = Xdxdz * (
@@ -804,7 +804,7 @@ void fdderivs(const int ex[3],
}
{
#define YSTEN6_8(JF, KF_DUMMY) \
(-fh[idx_fh_F_ord4(iF,JF-3,KF_DUMMY,ex)] + F9*fh[idx_fh_F_ord4(iF,JF-2,KF_DUMMY,ex)] - F45*fh[idx_fh_F_ord4(iF,JF-1,KF_DUMMY,ex)] + F45*fh[idx_fh_F_ord4(iF,JF+1,KF_DUMMY,ex)] - F9*fh[idx_fh_F_ord4(iF,JF+2,KF_DUMMY,ex)] + fh[idx_fh_F_ord4(iF,JF+3,KF_DUMMY,ex)])
(-fh[idx_fh_F_ord5(iF,JF-3,KF_DUMMY,ex)] + F9*fh[idx_fh_F_ord5(iF,JF-2,KF_DUMMY,ex)] - F45*fh[idx_fh_F_ord5(iF,JF-1,KF_DUMMY,ex)] + F45*fh[idx_fh_F_ord5(iF,JF+1,KF_DUMMY,ex)] - F9*fh[idx_fh_F_ord5(iF,JF+2,KF_DUMMY,ex)] + fh[idx_fh_F_ord5(iF,JF+3,KF_DUMMY,ex)])
fyz[p] = Xdydz * (
-YSTEN6_8(jF,kF-3) + F9*YSTEN6_8(jF,kF-2) - F45*YSTEN6_8(jF,kF-1) + F45*YSTEN6_8(jF,kF+1) - F9*YSTEN6_8(jF,kF+2) + YSTEN6_8(jF,kF+3));
#undef YSTEN6_8
@@ -826,44 +826,44 @@ void fdderivs(const int ex[3],
/* Diagonal: [-9,+128,-1008,+8064,-14350,+8064,-1008,+128,-9] / (5040*dx^2) */
fxx[p] = Edxdx * (
-(double)9 * fh[idx_fh_F_ord4(iF - 4, jF, kF, ex)] +
F128 * fh[idx_fh_F_ord4(iF - 3, jF, kF, ex)] -
F1008 * fh[idx_fh_F_ord4(iF - 2, jF, kF, ex)] +
F8064 * fh[idx_fh_F_ord4(iF - 1, jF, kF, ex)] -
F14350* fh[idx_fh_F_ord4(iF, jF, kF, ex)] +
F8064 * fh[idx_fh_F_ord4(iF + 1, jF, kF, ex)] -
F1008 * fh[idx_fh_F_ord4(iF + 2, jF, kF, ex)] +
F128 * fh[idx_fh_F_ord4(iF + 3, jF, kF, ex)] -
(double)9 * fh[idx_fh_F_ord4(iF + 4, jF, kF, ex)]);
-(double)9 * fh[idx_fh_F_ord5(iF - 4, jF, kF, ex)] +
F128 * fh[idx_fh_F_ord5(iF - 3, jF, kF, ex)] -
F1008 * fh[idx_fh_F_ord5(iF - 2, jF, kF, ex)] +
F8064 * fh[idx_fh_F_ord5(iF - 1, jF, kF, ex)] -
F14350* fh[idx_fh_F_ord5(iF, jF, kF, ex)] +
F8064 * fh[idx_fh_F_ord5(iF + 1, jF, kF, ex)] -
F1008 * fh[idx_fh_F_ord5(iF + 2, jF, kF, ex)] +
F128 * fh[idx_fh_F_ord5(iF + 3, jF, kF, ex)] -
(double)9 * fh[idx_fh_F_ord5(iF + 4, jF, kF, ex)]);
fyy[p] = Edydy * (
-(double)9 * fh[idx_fh_F_ord4(iF, jF - 4, kF, ex)] +
F128 * fh[idx_fh_F_ord4(iF, jF - 3, kF, ex)] -
F1008 * fh[idx_fh_F_ord4(iF, jF - 2, kF, ex)] +
F8064 * fh[idx_fh_F_ord4(iF, jF - 1, kF, ex)] -
F14350* fh[idx_fh_F_ord4(iF, jF, kF, ex)] +
F8064 * fh[idx_fh_F_ord4(iF, jF + 1, kF, ex)] -
F1008 * fh[idx_fh_F_ord4(iF, jF + 2, kF, ex)] +
F128 * fh[idx_fh_F_ord4(iF, jF + 3, kF, ex)] -
(double)9 * fh[idx_fh_F_ord4(iF, jF + 4, kF, ex)]);
-(double)9 * fh[idx_fh_F_ord5(iF, jF - 4, kF, ex)] +
F128 * fh[idx_fh_F_ord5(iF, jF - 3, kF, ex)] -
F1008 * fh[idx_fh_F_ord5(iF, jF - 2, kF, ex)] +
F8064 * fh[idx_fh_F_ord5(iF, jF - 1, kF, ex)] -
F14350* fh[idx_fh_F_ord5(iF, jF, kF, ex)] +
F8064 * fh[idx_fh_F_ord5(iF, jF + 1, kF, ex)] -
F1008 * fh[idx_fh_F_ord5(iF, jF + 2, kF, ex)] +
F128 * fh[idx_fh_F_ord5(iF, jF + 3, kF, ex)] -
(double)9 * fh[idx_fh_F_ord5(iF, jF + 4, kF, ex)]);
fzz[p] = Edzdz * (
-(double)9 * fh[idx_fh_F_ord4(iF, jF, kF - 4, ex)] +
F128 * fh[idx_fh_F_ord4(iF, jF, kF - 3, ex)] -
F1008 * fh[idx_fh_F_ord4(iF, jF, kF - 2, ex)] +
F8064 * fh[idx_fh_F_ord4(iF, jF, kF - 1, ex)] -
F14350* fh[idx_fh_F_ord4(iF, jF, kF, ex)] +
F8064 * fh[idx_fh_F_ord4(iF, jF, kF + 1, ex)] -
F1008 * fh[idx_fh_F_ord4(iF, jF, kF + 2, ex)] +
F128 * fh[idx_fh_F_ord4(iF, jF, kF + 3, ex)] -
(double)9 * fh[idx_fh_F_ord4(iF, jF, kF + 4, ex)]);
-(double)9 * fh[idx_fh_F_ord5(iF, jF, kF - 4, ex)] +
F128 * fh[idx_fh_F_ord5(iF, jF, kF - 3, ex)] -
F1008 * fh[idx_fh_F_ord5(iF, jF, kF - 2, ex)] +
F8064 * fh[idx_fh_F_ord5(iF, jF, kF - 1, ex)] -
F14350* fh[idx_fh_F_ord5(iF, jF, kF, ex)] +
F8064 * fh[idx_fh_F_ord5(iF, jF, kF + 1, ex)] -
F1008 * fh[idx_fh_F_ord5(iF, jF, kF + 2, ex)] +
F128 * fh[idx_fh_F_ord5(iF, jF, kF + 3, ex)] -
(double)9 * fh[idx_fh_F_ord5(iF, jF, kF + 4, ex)]);
/* Mixed: 9x9 outer product.
x-stencil: +3*f(i-4)-32*f(i-3)+168*f(i-2)-672*f(i-1)+672*f(i+1)-168*f(i+2)+32*f(i+3)-3*f(i+4)
y/z weights: same [+3,-32,+168,-672,+672,-168,+32,-3] / 705600 */
{
#define XSTEN8(JF, KF_DUMMY) \
(+(double)3*fh[idx_fh_F_ord4(iF-4,JF,KF_DUMMY,ex)] - F32*fh[idx_fh_F_ord4(iF-3,JF,KF_DUMMY,ex)] + F168*fh[idx_fh_F_ord4(iF-2,JF,KF_DUMMY,ex)] - F672*fh[idx_fh_F_ord4(iF-1,JF,KF_DUMMY,ex)] + F672*fh[idx_fh_F_ord4(iF+1,JF,KF_DUMMY,ex)] - F168*fh[idx_fh_F_ord4(iF+2,JF,KF_DUMMY,ex)] + F32*fh[idx_fh_F_ord4(iF+3,JF,KF_DUMMY,ex)] - (double)3*fh[idx_fh_F_ord4(iF+4,JF,KF_DUMMY,ex)])
(+(double)3*fh[idx_fh_F_ord5(iF-4,JF,KF_DUMMY,ex)] - F32*fh[idx_fh_F_ord5(iF-3,JF,KF_DUMMY,ex)] + F168*fh[idx_fh_F_ord5(iF-2,JF,KF_DUMMY,ex)] - F672*fh[idx_fh_F_ord5(iF-1,JF,KF_DUMMY,ex)] + F672*fh[idx_fh_F_ord5(iF+1,JF,KF_DUMMY,ex)] - F168*fh[idx_fh_F_ord5(iF+2,JF,KF_DUMMY,ex)] + F32*fh[idx_fh_F_ord5(iF+3,JF,KF_DUMMY,ex)] - (double)3*fh[idx_fh_F_ord5(iF+4,JF,KF_DUMMY,ex)])
fxy[p] = Edxdy * (
+(double)3*XSTEN8(jF-4,kF) - F32*XSTEN8(jF-3,kF) + F168*XSTEN8(jF-2,kF) - F672*XSTEN8(jF-1,kF) + F672*XSTEN8(jF+1,kF) - F168*XSTEN8(jF+2,kF) + F32*XSTEN8(jF+3,kF) - (double)3*XSTEN8(jF+4,kF));
@@ -875,7 +875,7 @@ void fdderivs(const int ex[3],
}
{
#define YSTEN8(JF, KF_DUMMY) \
(+(double)3*fh[idx_fh_F_ord4(iF,JF-4,KF_DUMMY,ex)] - F32*fh[idx_fh_F_ord4(iF,JF-3,KF_DUMMY,ex)] + F168*fh[idx_fh_F_ord4(iF,JF-2,KF_DUMMY,ex)] - F672*fh[idx_fh_F_ord4(iF,JF-1,KF_DUMMY,ex)] + F672*fh[idx_fh_F_ord4(iF,JF+1,KF_DUMMY,ex)] - F168*fh[idx_fh_F_ord4(iF,JF+2,KF_DUMMY,ex)] + F32*fh[idx_fh_F_ord4(iF,JF+3,KF_DUMMY,ex)] - (double)3*fh[idx_fh_F_ord4(iF,JF+4,KF_DUMMY,ex)])
(+(double)3*fh[idx_fh_F_ord5(iF,JF-4,KF_DUMMY,ex)] - F32*fh[idx_fh_F_ord5(iF,JF-3,KF_DUMMY,ex)] + F168*fh[idx_fh_F_ord5(iF,JF-2,KF_DUMMY,ex)] - F672*fh[idx_fh_F_ord5(iF,JF-1,KF_DUMMY,ex)] + F672*fh[idx_fh_F_ord5(iF,JF+1,KF_DUMMY,ex)] - F168*fh[idx_fh_F_ord5(iF,JF+2,KF_DUMMY,ex)] + F32*fh[idx_fh_F_ord5(iF,JF+3,KF_DUMMY,ex)] - (double)3*fh[idx_fh_F_ord5(iF,JF+4,KF_DUMMY,ex)])
fyz[p] = Edydz * (
+(double)3*YSTEN8(jF,kF-4) - F32*YSTEN8(jF,kF-3) + F168*YSTEN8(jF,kF-2) - F672*YSTEN8(jF,kF-1) + F672*YSTEN8(jF,kF+1) - F168*YSTEN8(jF,kF+2) + F32*YSTEN8(jF,kF+3) - (double)3*YSTEN8(jF,kF+4));