From 45ae92fd929c16196b6f41d5a71bf31362eba647 Mon Sep 17 00:00:00 2001 From: CGH0S7 <776459475@qq.com> Date: Thu, 14 May 2026 20:32:53 +0800 Subject: [PATCH] Fix eighth-order C derivative and lopsided stencils --- AMSS_NCKU_source/fdderivs_c.C | 118 +++++++++++++-------------- AMSS_NCKU_source/fdderivs_sh_c.C | 6 +- AMSS_NCKU_source/fderivs_c.C | 122 ++++++++++++++-------------- AMSS_NCKU_source/lopsided_c.C | 24 +++--- AMSS_NCKU_source/lopsided_kodis_c.C | 92 ++------------------- 5 files changed, 140 insertions(+), 222 deletions(-) diff --git a/AMSS_NCKU_source/fdderivs_c.C b/AMSS_NCKU_source/fdderivs_c.C index 05db174..4d14f8f 100644 --- a/AMSS_NCKU_source/fdderivs_c.C +++ b/AMSS_NCKU_source/fdderivs_c.C @@ -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)); diff --git a/AMSS_NCKU_source/fdderivs_sh_c.C b/AMSS_NCKU_source/fdderivs_sh_c.C index 955319d..83f517a 100644 --- a/AMSS_NCKU_source/fdderivs_sh_c.C +++ b/AMSS_NCKU_source/fdderivs_sh_c.C @@ -163,7 +163,7 @@ extern "C" void fdderivs_sh_(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+1; 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;if(in4)continue; + 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+1;const size_t p=idx_ex(i0,j0,k0,ex); fxx[p]=Sdxdx*(FH(iF-1,jF,kF)-TWO*FH(iF,jF,kF)+FH(iF+1,jF,kF)); fyy[p]=Sdydy*(FH(iF,jF-1,kF)-TWO*FH(iF,jF,kF)+FH(iF,jF+1,kF)); @@ -248,7 +248,7 @@ extern "C" void fdderivs_sh_(const int ex[3], /* 2nd-order pass */ 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+1; 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;if(in4)continue; + 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+1;const size_t p=idx_ex(i0,j0,k0,ex); fxx[p]=Sdxdx*(FH(iF-1,jF,kF)-TWO*FH(iF,jF,kF)+FH(iF+1,jF,kF)); fyy[p]=Sdydy*(FH(iF,jF-1,kF)-TWO*FH(iF,jF,kF)+FH(iF,jF+1,kF)); @@ -260,7 +260,7 @@ extern "C" void fdderivs_sh_(const int ex[3], /* 4th-order pass */ if(has4){for(int k0=k4_lo;k0<=k4_hi;++k0){const int kF=k0+1; 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;if(in6)continue; + 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+1;const size_t p=idx_ex(i0,j0,k0,ex); fxx[p]=Fdxdx*(-FH(iF-2,jF,kF)+F16*FH(iF-1,jF,kF)-F30*FH(iF,jF,kF)-FH(iF+2,jF,kF)+F16*FH(iF+1,jF,kF)); fyy[p]=Fdydy*(-FH(iF,jF-2,kF)+F16*FH(iF,jF-1,kF)-F30*FH(iF,jF,kF)-FH(iF,jF+2,kF)+F16*FH(iF,jF+1,kF)); diff --git a/AMSS_NCKU_source/fderivs_c.C b/AMSS_NCKU_source/fderivs_c.C index baa2eb6..bc59e37 100644 --- a/AMSS_NCKU_source/fderivs_c.C +++ b/AMSS_NCKU_source/fderivs_c.C @@ -398,7 +398,7 @@ void fderivs(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; @@ -484,14 +484,14 @@ void fderivs(const int ex[3], const size_t p = idx_ex(i0, j0, k0, ex); fx[p] = d2dx * ( - -fh[idx_fh_F_ord4(iF - 1, jF, kF, ex)] + - fh[idx_fh_F_ord4(iF + 1, jF, kF, ex)]); + -fh[idx_fh_F_ord5(iF - 1, jF, kF, ex)] + + fh[idx_fh_F_ord5(iF + 1, jF, kF, ex)]); fy[p] = d2dy * ( - -fh[idx_fh_F_ord4(iF, jF - 1, kF, ex)] + - fh[idx_fh_F_ord4(iF, jF + 1, kF, ex)]); + -fh[idx_fh_F_ord5(iF, jF - 1, kF, ex)] + + fh[idx_fh_F_ord5(iF, jF + 1, kF, ex)]); fz[p] = d2dz * ( - -fh[idx_fh_F_ord4(iF, jF, kF - 1, ex)] + - fh[idx_fh_F_ord4(iF, jF, kF + 1, ex)]); + -fh[idx_fh_F_ord5(iF, jF, kF - 1, ex)] + + fh[idx_fh_F_ord5(iF, jF, kF + 1, ex)]); } } } @@ -507,22 +507,22 @@ void fderivs(const int ex[3], const size_t p = idx_ex(i0, j0, k0, ex); fx[p] = d12dx * ( - fh[idx_fh_F_ord4(iF - 2, jF, kF, ex)] - - EIT * fh[idx_fh_F_ord4(iF - 1, jF, kF, ex)] + - EIT * fh[idx_fh_F_ord4(iF + 1, jF, kF, ex)] - - fh[idx_fh_F_ord4(iF + 2, jF, kF, ex)]); + fh[idx_fh_F_ord5(iF - 2, jF, kF, ex)] - + EIT * fh[idx_fh_F_ord5(iF - 1, jF, kF, ex)] + + EIT * fh[idx_fh_F_ord5(iF + 1, jF, kF, ex)] - + fh[idx_fh_F_ord5(iF + 2, jF, kF, ex)]); fy[p] = d12dy * ( - fh[idx_fh_F_ord4(iF, jF - 2, kF, ex)] - - EIT * fh[idx_fh_F_ord4(iF, jF - 1, kF, ex)] + - EIT * fh[idx_fh_F_ord4(iF, jF + 1, kF, ex)] - - fh[idx_fh_F_ord4(iF, jF + 2, kF, ex)]); + fh[idx_fh_F_ord5(iF, jF - 2, kF, ex)] - + EIT * fh[idx_fh_F_ord5(iF, jF - 1, kF, ex)] + + EIT * fh[idx_fh_F_ord5(iF, jF + 1, kF, ex)] - + fh[idx_fh_F_ord5(iF, jF + 2, kF, ex)]); fz[p] = d12dz * ( - fh[idx_fh_F_ord4(iF, jF, kF - 2, ex)] - - EIT * fh[idx_fh_F_ord4(iF, jF, kF - 1, ex)] + - EIT * fh[idx_fh_F_ord4(iF, jF, kF + 1, ex)] - - fh[idx_fh_F_ord4(iF, jF, kF + 2, ex)]); + fh[idx_fh_F_ord5(iF, jF, kF - 2, ex)] - + EIT * fh[idx_fh_F_ord5(iF, jF, kF - 1, ex)] + + EIT * fh[idx_fh_F_ord5(iF, jF, kF + 1, ex)] - + fh[idx_fh_F_ord5(iF, jF, kF + 2, ex)]); } } } @@ -538,28 +538,28 @@ void fderivs(const int ex[3], const size_t p = idx_ex(i0, j0, k0, ex); fx[p] = d60dx * ( - -fh[idx_fh_F_ord4(iF - 3, jF, kF, ex)] + - F9 * fh[idx_fh_F_ord4(iF - 2, jF, kF, ex)] - - F45 * fh[idx_fh_F_ord4(iF - 1, jF, kF, ex)] + - F45 * fh[idx_fh_F_ord4(iF + 1, jF, kF, ex)] - - F9 * fh[idx_fh_F_ord4(iF + 2, jF, kF, ex)] + - fh[idx_fh_F_ord4(iF + 3, jF, kF, ex)]); + -fh[idx_fh_F_ord5(iF - 3, jF, kF, ex)] + + F9 * fh[idx_fh_F_ord5(iF - 2, jF, kF, ex)] - + F45 * fh[idx_fh_F_ord5(iF - 1, jF, kF, ex)] + + F45 * fh[idx_fh_F_ord5(iF + 1, jF, kF, ex)] - + F9 * fh[idx_fh_F_ord5(iF + 2, jF, kF, ex)] + + fh[idx_fh_F_ord5(iF + 3, jF, kF, ex)]); fy[p] = d60dy * ( - -fh[idx_fh_F_ord4(iF, jF - 3, kF, ex)] + - F9 * fh[idx_fh_F_ord4(iF, jF - 2, kF, ex)] - - F45 * fh[idx_fh_F_ord4(iF, jF - 1, kF, ex)] + - F45 * fh[idx_fh_F_ord4(iF, jF + 1, kF, ex)] - - F9 * fh[idx_fh_F_ord4(iF, jF + 2, kF, ex)] + - fh[idx_fh_F_ord4(iF, jF + 3, kF, ex)]); + -fh[idx_fh_F_ord5(iF, jF - 3, kF, ex)] + + F9 * fh[idx_fh_F_ord5(iF, jF - 2, kF, ex)] - + F45 * fh[idx_fh_F_ord5(iF, jF - 1, kF, ex)] + + F45 * fh[idx_fh_F_ord5(iF, jF + 1, kF, ex)] - + F9 * fh[idx_fh_F_ord5(iF, jF + 2, kF, ex)] + + fh[idx_fh_F_ord5(iF, jF + 3, kF, ex)]); fz[p] = d60dz * ( - -fh[idx_fh_F_ord4(iF, jF, kF - 3, ex)] + - F9 * fh[idx_fh_F_ord4(iF, jF, kF - 2, ex)] - - F45 * fh[idx_fh_F_ord4(iF, jF, kF - 1, ex)] + - F45 * fh[idx_fh_F_ord4(iF, jF, kF + 1, ex)] - - F9 * fh[idx_fh_F_ord4(iF, jF, kF + 2, ex)] + - fh[idx_fh_F_ord4(iF, jF, kF + 3, ex)]); + -fh[idx_fh_F_ord5(iF, jF, kF - 3, ex)] + + F9 * fh[idx_fh_F_ord5(iF, jF, kF - 2, ex)] - + F45 * fh[idx_fh_F_ord5(iF, jF, kF - 1, ex)] + + F45 * fh[idx_fh_F_ord5(iF, jF, kF + 1, ex)] - + F9 * fh[idx_fh_F_ord5(iF, jF, kF + 2, ex)] + + fh[idx_fh_F_ord5(iF, jF, kF + 3, ex)]); } } } @@ -576,34 +576,34 @@ void fderivs(const int ex[3], const size_t p = idx_ex(i0, j0, k0, ex); fx[p] = d840dx * ( - +(double)3 * fh[idx_fh_F_ord4(iF - 4, jF, kF, ex)] - - F32 * fh[idx_fh_F_ord4(iF - 3, jF, kF, ex)] + - F168 * fh[idx_fh_F_ord4(iF - 2, jF, kF, ex)] - - F672 * fh[idx_fh_F_ord4(iF - 1, jF, kF, ex)] + - F672 * fh[idx_fh_F_ord4(iF + 1, jF, kF, ex)] - - F168 * fh[idx_fh_F_ord4(iF + 2, jF, kF, ex)] + - F32 * fh[idx_fh_F_ord4(iF + 3, jF, kF, ex)] - - (double)3 * fh[idx_fh_F_ord4(iF + 4, jF, kF, ex)]); + +(double)3 * fh[idx_fh_F_ord5(iF - 4, jF, kF, ex)] - + F32 * fh[idx_fh_F_ord5(iF - 3, jF, kF, ex)] + + F168 * fh[idx_fh_F_ord5(iF - 2, jF, kF, ex)] - + F672 * fh[idx_fh_F_ord5(iF - 1, jF, kF, ex)] + + F672 * fh[idx_fh_F_ord5(iF + 1, jF, kF, ex)] - + F168 * fh[idx_fh_F_ord5(iF + 2, jF, kF, ex)] + + F32 * fh[idx_fh_F_ord5(iF + 3, jF, kF, ex)] - + (double)3 * fh[idx_fh_F_ord5(iF + 4, jF, kF, ex)]); fy[p] = d840dy * ( - +(double)3 * fh[idx_fh_F_ord4(iF, jF - 4, kF, ex)] - - F32 * fh[idx_fh_F_ord4(iF, jF - 3, kF, ex)] + - F168 * fh[idx_fh_F_ord4(iF, jF - 2, kF, ex)] - - F672 * fh[idx_fh_F_ord4(iF, jF - 1, kF, ex)] + - F672 * fh[idx_fh_F_ord4(iF, jF + 1, kF, ex)] - - F168 * fh[idx_fh_F_ord4(iF, jF + 2, kF, ex)] + - F32 * fh[idx_fh_F_ord4(iF, jF + 3, kF, ex)] - - (double)3 * fh[idx_fh_F_ord4(iF, jF + 4, kF, ex)]); + +(double)3 * fh[idx_fh_F_ord5(iF, jF - 4, kF, ex)] - + F32 * fh[idx_fh_F_ord5(iF, jF - 3, kF, ex)] + + F168 * fh[idx_fh_F_ord5(iF, jF - 2, kF, ex)] - + F672 * fh[idx_fh_F_ord5(iF, jF - 1, kF, ex)] + + F672 * fh[idx_fh_F_ord5(iF, jF + 1, kF, ex)] - + F168 * fh[idx_fh_F_ord5(iF, jF + 2, kF, ex)] + + F32 * fh[idx_fh_F_ord5(iF, jF + 3, kF, ex)] - + (double)3 * fh[idx_fh_F_ord5(iF, jF + 4, kF, ex)]); fz[p] = d840dz * ( - +(double)3 * fh[idx_fh_F_ord4(iF, jF, kF - 4, ex)] - - F32 * fh[idx_fh_F_ord4(iF, jF, kF - 3, ex)] + - F168 * fh[idx_fh_F_ord4(iF, jF, kF - 2, ex)] - - F672 * fh[idx_fh_F_ord4(iF, jF, kF - 1, ex)] + - F672 * fh[idx_fh_F_ord4(iF, jF, kF + 1, ex)] - - F168 * fh[idx_fh_F_ord4(iF, jF, kF + 2, ex)] + - F32 * fh[idx_fh_F_ord4(iF, jF, kF + 3, ex)] - - (double)3 * fh[idx_fh_F_ord4(iF, jF, kF + 4, ex)]); + +(double)3 * fh[idx_fh_F_ord5(iF, jF, kF - 4, ex)] - + F32 * fh[idx_fh_F_ord5(iF, jF, kF - 3, ex)] + + F168 * fh[idx_fh_F_ord5(iF, jF, kF - 2, ex)] - + F672 * fh[idx_fh_F_ord5(iF, jF, kF - 1, ex)] + + F672 * fh[idx_fh_F_ord5(iF, jF, kF + 1, ex)] - + F168 * fh[idx_fh_F_ord5(iF, jF, kF + 2, ex)] + + F32 * fh[idx_fh_F_ord5(iF, jF, kF + 3, ex)] - + (double)3 * fh[idx_fh_F_ord5(iF, jF, kF + 4, ex)]); } } } diff --git a/AMSS_NCKU_source/lopsided_c.C b/AMSS_NCKU_source/lopsided_c.C index 2dd839e..60e5798 100644 --- a/AMSS_NCKU_source/lopsided_c.C +++ b/AMSS_NCKU_source/lopsided_c.C @@ -503,25 +503,25 @@ void lopsided(const int ex[3], f_rhs[p] += sfx * d2dx * ( -fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]); } else if (sfx < ZEO) { - if ((i0-5)>=iminF && i0<=ex1-2) // i-5>=imin && i+3<=imax + if ((i0-4)>=iminF && i0<=ex1-4) f_rhs[p] -= sfx * d840dx * ( -F5*fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]+F60*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)] -F420*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)] +F1050*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]-F420*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)] +F140*fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]-F30*fh[idx_fh_F_ord5(iF-4,jF,kF,ex)] +F3*fh[idx_fh_F_ord5(iF-5,jF,kF,ex)]); - else if ((i0-4)>=iminF && i0<=ex1-2) // 8th centered + else if ((i0-3)>=iminF && i0<=ex1-5) // 8th centered f_rhs[p] -= sfx * d840dx * ( +F3*fh[idx_fh_F_ord5(iF-4,jF,kF,ex)]-F32*fh[idx_fh_F_ord5(iF-3,jF,kF,ex)] +F168*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-F672*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)] +F672*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-F168*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)] +F32*fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]-F3*fh[idx_fh_F_ord5(iF+4,jF,kF,ex)]); - else if ((i0-3)>=iminF && i0<=ex1-2) // 6th centered + else if ((i0-2)>=iminF && i0<=ex1-4) // 6th centered f_rhs[p] -= sfx * d60dx * ( -fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]+F9*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)] -F45*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+F45*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)] -F9*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]); - else if ((i0-2)>=iminF && i0<=ex1-2) // 4th centered + else if ((i0-1)>=iminF && i0<=ex1-3) // 4th centered f_rhs[p] -= sfx * d12dx * ( fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-EIT*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)] +EIT*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]); @@ -543,13 +543,13 @@ void lopsided(const int ex[3], else if (j0<=ex2-2 && j0>=jminF) f_rhs[p] += sfy * d2dy*(-fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]); } else if (sfy < ZEO) { - if ((j0-5)>=jminF && j0<=ex2-2) + if ((j0-4)>=jminF && j0<=ex2-4) f_rhs[p] -= sfy * d840dy*(-F5*fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]+F60*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]-F420*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)]+F1050*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]-F420*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]+F140*fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]-F30*fh[idx_fh_F_ord5(iF,jF-4,kF,ex)]+F3*fh[idx_fh_F_ord5(iF,jF-5,kF,ex)]); - else if ((j0-4)>=jminF && j0<=ex2-2) + else if ((j0-3)>=jminF && j0<=ex2-5) f_rhs[p] -= sfy * d840dy*(+F3*fh[idx_fh_F_ord5(iF,jF-4,kF,ex)]-F32*fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]+F168*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-F672*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+F672*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F168*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]+F32*fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]-F3*fh[idx_fh_F_ord5(iF,jF+4,kF,ex)]); - else if ((j0-3)>=jminF && j0<=ex2-2) + else if ((j0-2)>=jminF && j0<=ex2-4) f_rhs[p] -= sfy * d60dy*(-fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]+F9*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-F45*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+F45*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F9*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]); - else if ((j0-2)>=jminF && j0<=ex2-2) + else if ((j0-1)>=jminF && j0<=ex2-3) f_rhs[p] -= sfy * d12dy*(fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-EIT*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+EIT*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]); else if ((j0-1)>=jminF && j0<=ex2-2) f_rhs[p] -= sfy * d2dy*(-fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]); @@ -568,13 +568,13 @@ void lopsided(const int ex[3], else if (k0<=ex3-2 && k0>=kminF) f_rhs[p] += sfz * d2dz*(-fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]); } else if (sfz < ZEO) { - if ((k0-5)>=kminF && k0<=ex3-2) + if ((k0-4)>=kminF && k0<=ex3-4) f_rhs[p] -= sfz * d840dz*(-F5*fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]+F60*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]-F420*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)]+F1050*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]-F420*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]+F140*fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]-F30*fh[idx_fh_F_ord5(iF,jF,kF-4,ex)]+F3*fh[idx_fh_F_ord5(iF,jF,kF-5,ex)]); - else if ((k0-4)>=kminF && k0<=ex3-2) + else if ((k0-3)>=kminF && k0<=ex3-5) f_rhs[p] -= sfz * d840dz*(+F3*fh[idx_fh_F_ord5(iF,jF,kF-4,ex)]-F32*fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]+F168*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-F672*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+F672*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F168*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]+F32*fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]-F3*fh[idx_fh_F_ord5(iF,jF,kF+4,ex)]); - else if ((k0-3)>=kminF && k0<=ex3-2) + else if ((k0-2)>=kminF && k0<=ex3-4) f_rhs[p] -= sfz * d60dz*(-fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]+F9*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-F45*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+F45*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F9*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]); - else if ((k0-2)>=kminF && k0<=ex3-2) + else if ((k0-1)>=kminF && k0<=ex3-3) f_rhs[p] -= sfz * d12dz*(fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-EIT*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+EIT*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]); else if ((k0-1)>=kminF && k0<=ex3-2) f_rhs[p] -= sfz * d2dz*(-fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]); diff --git a/AMSS_NCKU_source/lopsided_kodis_c.C b/AMSS_NCKU_source/lopsided_kodis_c.C index a8c9f83..57f7bfa 100644 --- a/AMSS_NCKU_source/lopsided_kodis_c.C +++ b/AMSS_NCKU_source/lopsided_kodis_c.C @@ -3,7 +3,9 @@ /* * C 版 lopsided_kodis — combined upwind advection + KO dissipation. - * Uses one shared symmetry_bd buffer (ord = ghost_width for both components). + * Uses one shared symmetry_bd buffer (ord = ghost_width for both components) + * where a stable merged stencil is available. The 8th-order path delegates to + * the separate lopsided + kodis kernels, matching the original Fortran flow. * * FD order selection via ghost_width: * 2 → 2nd-order advection + r=2 KO (cof=16, sign=-) @@ -292,92 +294,8 @@ void lopsided_kodis(const int ex[3], } #elif (ghost_width == 5) { - const int ord = 5; - int iminF = 1, jminF = 1, kminF = 1; - if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -4; - if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -4; - if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -4; - - const size_t nx = (size_t)ex1 + ord; - const size_t ny = (size_t)ex2 + ord; - const size_t nz = (size_t)ex3 + ord; - double *fh = (double*)malloc(nx*ny*nz*sizeof(double)); - if (!fh) return; - symmetry_bd(ord, ex, f, fh, SoA); - - const double d840dx=ONE/F840/dX, d840dy=ONE/F840/dY, d840dz=ONE/F840/dZ; - const double d60dx=ONE/F60/dX, d60dy=ONE/F60/dY, d60dz=ONE/F60/dZ; - const double d12dx=ONE/F12/dX, d12dy=ONE/F12/dY, d12dz=ONE/F12/dZ; - const double d2dx=ONE/TWO/dX, d2dy=ONE/TWO/dY, d2dz=ONE/TWO/dZ; - - /* ---- advection (8th-order lopsided) ---- */ - for (int k0=0;k0<=ex3-2;++k0) { const int kF=k0+1; - for (int j0=0;j0<=ex2-2;++j0) { const int jF=j0+1; - for (int i0=0;i0<=ex1-2;++i0) { const int iF=i0+1; - const size_t p=idx_ex(i0,j0,k0,ex); - const double sfx=Sfx[p]; - if (sfx>ZEO) { - if (i0<=ex1-6&&(i0-2)>=iminF) f_rhs[p]+=sfx*d840dx*(-F5*fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]+F60*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-F420*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)]+F1050*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-F420*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]+F140*fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]-F30*fh[idx_fh_F_ord5(iF+4,jF,kF,ex)]+F3*fh[idx_fh_F_ord5(iF+5,jF,kF,ex)]); - else if (i0<=ex1-5&&(i0-3)>=iminF) f_rhs[p]+=sfx*d840dx*(+F3*fh[idx_fh_F_ord5(iF-4,jF,kF,ex)]-F32*fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]+F168*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-F672*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+F672*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-F168*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]+F32*fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]-F3*fh[idx_fh_F_ord5(iF+4,jF,kF,ex)]); - else if (i0<=ex1-4&&(i0-2)>=iminF) f_rhs[p]+=sfx*d60dx*(-fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]+F9*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-F45*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+F45*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-F9*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]); - else if (i0<=ex1-3&&(i0-1)>=iminF) f_rhs[p]+=sfx*d12dx*(fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-EIT*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+EIT*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]); - else if (i0<=ex1-2&&i0>=iminF) f_rhs[p]+=sfx*d2dx*(-fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]); - } else if (sfx=iminF&&i0<=ex1-2) f_rhs[p]-=sfx*d840dx*(-F5*fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]+F60*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]-F420*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)]+F1050*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]-F420*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]+F140*fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]-F30*fh[idx_fh_F_ord5(iF-4,jF,kF,ex)]+F3*fh[idx_fh_F_ord5(iF-5,jF,kF,ex)]); - else if ((i0-4)>=iminF&&i0<=ex1-2) f_rhs[p]-=sfx*d840dx*(+F3*fh[idx_fh_F_ord5(iF-4,jF,kF,ex)]-F32*fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]+F168*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-F672*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+F672*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-F168*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]+F32*fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]-F3*fh[idx_fh_F_ord5(iF+4,jF,kF,ex)]); - else if ((i0-3)>=iminF&&i0<=ex1-2) f_rhs[p]-=sfx*d60dx*(-fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]+F9*fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-F45*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+F45*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-F9*fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+3,jF,kF,ex)]); - else if ((i0-2)>=iminF&&i0<=ex1-2) f_rhs[p]-=sfx*d12dx*(fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]-EIT*fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+EIT*fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]-fh[idx_fh_F_ord5(iF+2,jF,kF,ex)]); - else if ((i0-1)>=iminF&&i0<=ex1-2) f_rhs[p]-=sfx*d2dx*(-fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+1,jF,kF,ex)]); - } - const double sfy=Sfy[p]; - if (sfy>ZEO) { - if (j0<=ex2-6&&(j0-2)>=jminF) f_rhs[p]+=sfy*d840dy*(-F5*fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]+F60*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-F420*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)]+F1050*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F420*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]+F140*fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]-F30*fh[idx_fh_F_ord5(iF,jF+4,kF,ex)]+F3*fh[idx_fh_F_ord5(iF,jF+5,kF,ex)]); - else if (j0<=ex2-5&&(j0-3)>=jminF) f_rhs[p]+=sfy*d840dy*(+F3*fh[idx_fh_F_ord5(iF,jF-4,kF,ex)]-F32*fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]+F168*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-F672*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+F672*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F168*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]+F32*fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]-F3*fh[idx_fh_F_ord5(iF,jF+4,kF,ex)]); - else if (j0<=ex2-4&&(j0-2)>=jminF) f_rhs[p]+=sfy*d60dy*(-fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]+F9*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-F45*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+F45*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F9*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]); - else if (j0<=ex2-3&&(j0-1)>=jminF) f_rhs[p]+=sfy*d12dy*(fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-EIT*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+EIT*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]); - else if (j0<=ex2-2&&j0>=jminF) f_rhs[p]+=sfy*d2dy*(-fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]); - } else if (sfy=jminF&&j0<=ex2-2) f_rhs[p]-=sfy*d840dy*(-F5*fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]+F60*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]-F420*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)]+F1050*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]-F420*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]+F140*fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]-F30*fh[idx_fh_F_ord5(iF,jF-4,kF,ex)]+F3*fh[idx_fh_F_ord5(iF,jF-5,kF,ex)]); - else if ((j0-4)>=jminF&&j0<=ex2-2) f_rhs[p]-=sfy*d840dy*(+F3*fh[idx_fh_F_ord5(iF,jF-4,kF,ex)]-F32*fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]+F168*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-F672*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+F672*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F168*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]+F32*fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]-F3*fh[idx_fh_F_ord5(iF,jF+4,kF,ex)]); - else if ((j0-3)>=jminF&&j0<=ex2-2) f_rhs[p]-=sfy*d60dy*(-fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]+F9*fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-F45*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+F45*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-F9*fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+3,kF,ex)]); - else if ((j0-2)>=jminF&&j0<=ex2-2) f_rhs[p]-=sfy*d12dy*(fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]-EIT*fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+EIT*fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]-fh[idx_fh_F_ord5(iF,jF+2,kF,ex)]); - else if ((j0-1)>=jminF&&j0<=ex2-2) f_rhs[p]-=sfy*d2dy*(-fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+1,kF,ex)]); - } - const double sfz=Sfz[p]; - if (sfz>ZEO) { - if (k0<=ex3-6&&(k0-2)>=kminF) f_rhs[p]+=sfz*d840dz*(-F5*fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]+F60*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-F420*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)]+F1050*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F420*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]+F140*fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]-F30*fh[idx_fh_F_ord5(iF,jF,kF+4,ex)]+F3*fh[idx_fh_F_ord5(iF,jF,kF+5,ex)]); - else if (k0<=ex3-5&&(k0-3)>=kminF) f_rhs[p]+=sfz*d840dz*(+F3*fh[idx_fh_F_ord5(iF,jF,kF-4,ex)]-F32*fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]+F168*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-F672*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+F672*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F168*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]+F32*fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]-F3*fh[idx_fh_F_ord5(iF,jF,kF+4,ex)]); - else if (k0<=ex3-4&&(k0-2)>=kminF) f_rhs[p]+=sfz*d60dz*(-fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]+F9*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-F45*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+F45*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F9*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]); - else if (k0<=ex3-3&&(k0-1)>=kminF) f_rhs[p]+=sfz*d12dz*(fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-EIT*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+EIT*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]); - else if (k0<=ex3-2&&k0>=kminF) f_rhs[p]+=sfz*d2dz*(-fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]); - } else if (sfz=kminF&&k0<=ex3-2) f_rhs[p]-=sfz*d840dz*(-F5*fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]+F60*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]-F420*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F378*fh[idx_fh_F_ord5(iF,jF,kF,ex)]+F1050*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]-F420*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]+F140*fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]-F30*fh[idx_fh_F_ord5(iF,jF,kF-4,ex)]+F3*fh[idx_fh_F_ord5(iF,jF,kF-5,ex)]); - else if ((k0-4)>=kminF&&k0<=ex3-2) f_rhs[p]-=sfz*d840dz*(+F3*fh[idx_fh_F_ord5(iF,jF,kF-4,ex)]-F32*fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]+F168*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-F672*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+F672*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F168*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]+F32*fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]-F3*fh[idx_fh_F_ord5(iF,jF,kF+4,ex)]); - else if ((k0-3)>=kminF&&k0<=ex3-2) f_rhs[p]-=sfz*d60dz*(-fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]+F9*fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-F45*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+F45*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-F9*fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+3,ex)]); - else if ((k0-2)>=kminF&&k0<=ex3-2) f_rhs[p]-=sfz*d12dz*(fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]-EIT*fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+EIT*fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]-fh[idx_fh_F_ord5(iF,jF,kF+2,ex)]); - else if ((k0-1)>=kminF&&k0<=ex3-2) f_rhs[p]-=sfz*d2dz*(-fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+1,ex)]); - } - }}} - - /* ---- KO dissipation (r=5, cof=1024, sign=+) ---- */ - if (eps > ZEO) { - const double cof = 1024.0; - const double F10k=10.0, F45k=45.0, F120=120.0, F210=210.0, F252=252.0; - const int i0_lo=(iminF+4>0)?iminF+4:0, j0_lo=(jminF+4>0)?jminF+4:0, k0_lo=(kminF+4>0)?kminF+4:0; - const int i0_hi=imaxF-6, j0_hi=jmaxF-6, k0_hi=kmaxF-6; - if (!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)) { - for (int k0=k0_lo;k0<=k0_hi;++k0) { const int kF=k0+1; - for (int j0=j0_lo;j0<=j0_hi;++j0) { const int jF=j0+1; - for (int i0=i0_lo;i0<=i0_hi;++i0) { const int iF=i0+1; - const size_t p=idx_ex(i0,j0,k0,ex); - const double Dx=((fh[idx_fh_F_ord5(iF-5,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+5,jF,kF,ex)])-F10k*(fh[idx_fh_F_ord5(iF-4,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+4,jF,kF,ex)])+F45k*(fh[idx_fh_F_ord5(iF-3,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+3,jF,kF,ex)])-F120*(fh[idx_fh_F_ord5(iF-2,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+2,jF,kF,ex)])+F210*(fh[idx_fh_F_ord5(iF-1,jF,kF,ex)]+fh[idx_fh_F_ord5(iF+1,jF,kF,ex)])-F252*fh[idx_fh_F_ord5(iF,jF,kF,ex)])/dX; - const double Dy=((fh[idx_fh_F_ord5(iF,jF-5,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+5,kF,ex)])-F10k*(fh[idx_fh_F_ord5(iF,jF-4,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+4,kF,ex)])+F45k*(fh[idx_fh_F_ord5(iF,jF-3,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+3,kF,ex)])-F120*(fh[idx_fh_F_ord5(iF,jF-2,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+2,kF,ex)])+F210*(fh[idx_fh_F_ord5(iF,jF-1,kF,ex)]+fh[idx_fh_F_ord5(iF,jF+1,kF,ex)])-F252*fh[idx_fh_F_ord5(iF,jF,kF,ex)])/dY; - const double Dz=((fh[idx_fh_F_ord5(iF,jF,kF-5,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+5,ex)])-F10k*(fh[idx_fh_F_ord5(iF,jF,kF-4,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+4,ex)])+F45k*(fh[idx_fh_F_ord5(iF,jF,kF-3,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+3,ex)])-F120*(fh[idx_fh_F_ord5(iF,jF,kF-2,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+2,ex)])+F210*(fh[idx_fh_F_ord5(iF,jF,kF-1,ex)]+fh[idx_fh_F_ord5(iF,jF,kF+1,ex)])-F252*fh[idx_fh_F_ord5(iF,jF,kF,ex)])/dZ; - f_rhs[p] += (eps/cof)*(Dx+Dy+Dz); - }}} - } - } - free(fh); + lopsided(ex, X, Y, Z, f, f_rhs, Sfx, Sfy, Sfz, Symmetry, SoA); + if (eps > ZEO) kodis(ex, X, Y, Z, f, f_rhs, SoA, Symmetry, eps); return; } #else