AMSS-NCKU/AMSS_NCKU_source/fderivs_c.C

#include "tool.h"

/*
 * C 版 fderivs
 *
 * Fortran:
 * subroutine fderivs(ex,f,fx,fy,fz,X,Y,Z,SYM1,SYM2,SYM3,symmetry,onoff)
 *
 * 约定：
 *   f, fx, fy, fz: ex1*ex2*ex3，按 idx_ex 布局
 *   X: ex1, Y: ex2, Z: ex3
 */
void fderivs(const int ex[3],
             const double *f,
             double *fx, double *fy, double *fz,
             const double *X, const double *Y, const double *Z,
             double SYM1, double SYM2, double SYM3,
             int Symmetry, int onoff)
{
    (void)onoff; // Fortran 里没用到

    const double ZEO = 0.0, ONE = 1.0;
    const double TWO = 2.0, EIT = 8.0;
    const double F12 = 12.0;

    const int NO_SYMM = 0, EQ_SYMM = 1; // OCTANT=2 在本子程序里不直接用

    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];

    // dX = X(2)-X(1) -> C: X[1]-X[0]
    const double dX = X[1] - X[0];
    const double dY = Y[1] - Y[0];
    const double dZ = Z[1] - Z[0];

    // Fortran 1-based bounds
    const int imaxF = ex1;
    const int jmaxF = ex2;
    const int kmaxF = ex3;

    int iminF = 1, jminF = 1, kminF = 1;
    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -1;
    if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -1;
    if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -1;

    // SoA(1:3) = SYM1,SYM2,SYM3
    const double SoA[3] = { SYM1, SYM2, SYM3 };

    // fh: (ex1+2)*(ex2+2)*(ex3+2) because ord=2
    const size_t nx = (size_t)ex1 + 2;
    const size_t ny = (size_t)ex2 + 2;
    const size_t nz = (size_t)ex3 + 2;
    const size_t fh_size = nx * ny * nz;
    static double *fh = NULL;
    static size_t cap = 0;

    if (fh_size > cap) {
        free(fh);
        fh = (double*)aligned_alloc(64, fh_size * sizeof(double));
        cap = fh_size;
    }
    // double *fh = (double*)malloc(fh_size * sizeof(double));
    if (!fh) return;

    // call symmetry_bd(2,ex,f,fh,SoA)
    symmetry_bd(2, ex, f, fh, SoA);

    const double d12dx = ONE / F12 / dX;
    const double d12dy = ONE / F12 / dY;
    const double d12dz = ONE / F12 / dZ;

    const double d2dx  = ONE / TWO / dX;
    const double d2dy  = ONE / TWO / dY;
    const double d2dz  = ONE / TWO / dZ;

    // fx = fy = fz = 0
    const size_t all = (size_t)ex1 * (size_t)ex2 * (size_t)ex3;
    for (size_t p = 0; p < all; ++p) {
        fx[p] = ZEO;
        fy[p] = ZEO;
        fz[p] = ZEO;
    }

    /*
     * 两段式：
     * 1) 先在二阶可用区域计算二阶模板
     * 2) 再在高阶可用区域覆盖为四阶模板
     *
     * 与原 if/elseif 逻辑等价，但减少逐点分支判断。
     */
    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;
    const int j2_hi = ex2 - 2;
    const int k2_hi = ex3 - 2;

    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;
    const int j4_hi = ex2 - 3;
    const int k4_hi = ex3 - 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) {
                    const int iF = i0 + 1;
                    const size_t p = idx_ex(i0, j0, k0, ex);

                    fx[p] = d2dx * (
                        -fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] +
                         fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)]
                    );

                    fy[p] = d2dy * (
                        -fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] +
                         fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)]
                    );

                    fz[p] = d2dz * (
                        -fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] +
                         fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)]
                    );
                }
            }
        }
    }

    if (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi) {
        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) {
                    const int iF = i0 + 1;
                    const size_t p = idx_ex(i0, j0, k0, ex);

                    fx[p] = d12dx * (
                        fh[idx_fh_F_ord2(iF - 2, jF,     kF,     ex)] -
                        EIT * fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] +
                        EIT * fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)] -
                        fh[idx_fh_F_ord2(iF + 2, jF,     kF,     ex)]
                    );

                    fy[p] = d12dy * (
                        fh[idx_fh_F_ord2(iF,     jF - 2, kF,     ex)] -
                        EIT * fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] +
                        EIT * fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)] -
                        fh[idx_fh_F_ord2(iF,     jF + 2, kF,     ex)]
                    );

                    fz[p] = d12dz * (
                        fh[idx_fh_F_ord2(iF,     jF,     kF - 2, ex)] -
                        EIT * fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] +
                        EIT * fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)] -
                        fh[idx_fh_F_ord2(iF,     jF,     kF + 2, ex)]
                    );
                }
            }
        }
    }

    // free(fh);
}