247 lines
9.1 KiB
C
247 lines
9.1 KiB
C
#ifndef SHARE_FUNC_H
|
||
#define SHARE_FUNC_H
|
||
|
||
#include <stdlib.h>
|
||
#include <stddef.h>
|
||
#include <math.h>
|
||
#include <stdio.h>
|
||
#include <string.h>
|
||
/* 主网格：0-based -> 1D */
|
||
static inline size_t idx_ex(int i0, int j0, int k0, const int ex[3]) {
|
||
const int ex1 = ex[0], ex2 = ex[1];
|
||
return (size_t)i0 + (size_t)j0 * (size_t)ex1 + (size_t)k0 * (size_t)ex1 * (size_t)ex2;
|
||
}
|
||
|
||
/*
|
||
* fh 对应 Fortran: fh(-1:ex1, -1:ex2, -1:ex3)
|
||
* ord=2 => shift=1
|
||
* iF/jF/kF 为 Fortran 索引（可为 -1,0,1..ex）
|
||
*/
|
||
static inline size_t idx_fh_F_ord2(int iF, int jF, int kF, const int ex[3]) {
|
||
const int shift = 1;
|
||
const int nx = ex[0] + 2; // ex1 + ord
|
||
const int ny = ex[1] + 2;
|
||
|
||
const int ii = iF + shift; // 0..ex1+1
|
||
const int jj = jF + shift; // 0..ex2+1
|
||
const int kk = kF + shift; // 0..ex3+1
|
||
|
||
return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
|
||
}
|
||
|
||
/*
|
||
* fh 对应 Fortran: fh(-2:ex1, -2:ex2, -2:ex3)
|
||
* ord=3 => shift=2
|
||
* iF/jF/kF 是 Fortran 索引（可为负）
|
||
*/
|
||
static inline size_t idx_fh_F(int iF, int jF, int kF, const int ex[3]) {
|
||
const int shift = 2; // ord=3 -> -2..ex
|
||
const int nx = ex[0] + 3; // ex1 + ord
|
||
const int ny = ex[1] + 3;
|
||
|
||
const int ii = iF + shift; // 0..ex1+2
|
||
const int jj = jF + shift; // 0..ex2+2
|
||
const int kk = kF + shift; // 0..ex3+2
|
||
|
||
return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
|
||
}
|
||
|
||
/*
|
||
* func: (1..extc1, 1..extc2, 1..extc3) 1-based in Fortran
|
||
* funcc: (-ord+1..extc1, -ord+1..extc2, -ord+1..extc3) in Fortran
|
||
*
|
||
* C 里我们把：
|
||
* func 视为 0-based: i0=0..extc1-1, j0=0..extc2-1, k0=0..extc3-1
|
||
* funcc 用“平移下标”存为一维数组：
|
||
* iF in [-ord+1..extc1] -> ii = iF + (ord-1) in [0..extc1+ord-1]
|
||
* 总长度 nx = extc1 + ord
|
||
* 同理 ny = extc2 + ord, nz = extc3 + ord
|
||
*/
|
||
|
||
static inline size_t idx_func0(int i0, int j0, int k0, const int extc[3]) {
|
||
const int nx = extc[0], ny = extc[1];
|
||
return (size_t)i0 + (size_t)j0 * (size_t)nx + (size_t)k0 * (size_t)nx * (size_t)ny;
|
||
}
|
||
|
||
static inline size_t idx_funcc_F(int iF, int jF, int kF, int ord, const int extc[3]) {
|
||
const int shift = ord - 1; // iF = -shift .. extc1
|
||
const int nx = extc[0] + ord; // [-shift..extc1] 共 extc1+ord 个
|
||
const int ny = extc[1] + ord;
|
||
|
||
const int ii = iF + shift; // 0..extc1+shift
|
||
const int jj = jF + shift; // 0..extc2+shift
|
||
const int kk = kF + shift; // 0..extc3+shift
|
||
|
||
return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
|
||
}
|
||
|
||
/*
|
||
* 等价于 Fortran:
|
||
* funcc(1:extc1,1:extc2,1:extc3)=func
|
||
* do i=0,ord-1
|
||
* funcc(-i,1:extc2,1:extc3) = funcc(i+1,1:extc2,1:extc3)*SoA(1)
|
||
* enddo
|
||
* do i=0,ord-1
|
||
* funcc(:,-i,1:extc3) = funcc(:,i+1,1:extc3)*SoA(2)
|
||
* enddo
|
||
* do i=0,ord-1
|
||
* funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
|
||
* enddo
|
||
*/
|
||
static inline void symmetry_bd_impl(int ord,
|
||
int shift,
|
||
const int extc[3],
|
||
const double *__restrict func,
|
||
double *__restrict funcc,
|
||
const double SoA[3])
|
||
{
|
||
const int extc1 = extc[0], extc2 = extc[1], extc3 = extc[2];
|
||
const int nx = extc1 + ord;
|
||
const int ny = extc2 + ord;
|
||
|
||
const size_t snx = (size_t)nx;
|
||
const size_t splane = (size_t)nx * (size_t)ny;
|
||
const size_t interior_i = (size_t)shift + 1u; /* iF = 1 */
|
||
const size_t interior_j = ((size_t)shift + 1u) * snx; /* jF = 1 */
|
||
const size_t interior_k = ((size_t)shift + 1u) * splane; /* kF = 1 */
|
||
const size_t interior0 = interior_k + interior_j + interior_i;
|
||
|
||
/* 1) funcc(1:extc1,1:extc2,1:extc3) = func */
|
||
for (int k0 = 0; k0 < extc3; ++k0) {
|
||
const double *src_k = func + (size_t)k0 * (size_t)extc2 * (size_t)extc1;
|
||
const size_t dst_k0 = interior0 + (size_t)k0 * splane;
|
||
for (int j0 = 0; j0 < extc2; ++j0) {
|
||
const double *src = src_k + (size_t)j0 * (size_t)extc1;
|
||
double *dst = funcc + dst_k0 + (size_t)j0 * snx;
|
||
memcpy(dst, src, (size_t)extc1 * sizeof(double));
|
||
}
|
||
}
|
||
|
||
/* 2) funcc(-i,1:extc2,1:extc3) = funcc(i+1,1:extc2,1:extc3)*SoA(1) */
|
||
const double s1 = SoA[0];
|
||
if (s1 == 1.0) {
|
||
for (int ii = 0; ii < ord; ++ii) {
|
||
const size_t dst_i = (size_t)(shift - ii);
|
||
const size_t src_i = (size_t)(shift + ii + 1);
|
||
for (int k0 = 0; k0 < extc3; ++k0) {
|
||
const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
|
||
for (int j0 = 0; j0 < extc2; ++j0) {
|
||
const size_t off = kbase + (size_t)j0 * snx;
|
||
funcc[off + dst_i] = funcc[off + src_i];
|
||
}
|
||
}
|
||
}
|
||
} else if (s1 == -1.0) {
|
||
for (int ii = 0; ii < ord; ++ii) {
|
||
const size_t dst_i = (size_t)(shift - ii);
|
||
const size_t src_i = (size_t)(shift + ii + 1);
|
||
for (int k0 = 0; k0 < extc3; ++k0) {
|
||
const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
|
||
for (int j0 = 0; j0 < extc2; ++j0) {
|
||
const size_t off = kbase + (size_t)j0 * snx;
|
||
funcc[off + dst_i] = -funcc[off + src_i];
|
||
}
|
||
}
|
||
}
|
||
} else {
|
||
for (int ii = 0; ii < ord; ++ii) {
|
||
const size_t dst_i = (size_t)(shift - ii);
|
||
const size_t src_i = (size_t)(shift + ii + 1);
|
||
for (int k0 = 0; k0 < extc3; ++k0) {
|
||
const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
|
||
for (int j0 = 0; j0 < extc2; ++j0) {
|
||
const size_t off = kbase + (size_t)j0 * snx;
|
||
funcc[off + dst_i] = funcc[off + src_i] * s1;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
/* 3) funcc(:,-j,1:extc3) = funcc(:,j+1,1:extc3)*SoA(2) */
|
||
const double s2 = SoA[1];
|
||
if (s2 == 1.0) {
|
||
for (int jj = 0; jj < ord; ++jj) {
|
||
const size_t dst_j = (size_t)(shift - jj) * snx;
|
||
const size_t src_j = (size_t)(shift + jj + 1) * snx;
|
||
for (int k0 = 0; k0 < extc3; ++k0) {
|
||
const size_t kbase = interior_k + (size_t)k0 * splane;
|
||
double *dst = funcc + kbase + dst_j;
|
||
const double *src = funcc + kbase + src_j;
|
||
for (int i = 0; i < nx; ++i) dst[i] = src[i];
|
||
}
|
||
}
|
||
} else if (s2 == -1.0) {
|
||
for (int jj = 0; jj < ord; ++jj) {
|
||
const size_t dst_j = (size_t)(shift - jj) * snx;
|
||
const size_t src_j = (size_t)(shift + jj + 1) * snx;
|
||
for (int k0 = 0; k0 < extc3; ++k0) {
|
||
const size_t kbase = interior_k + (size_t)k0 * splane;
|
||
double *dst = funcc + kbase + dst_j;
|
||
const double *src = funcc + kbase + src_j;
|
||
for (int i = 0; i < nx; ++i) dst[i] = -src[i];
|
||
}
|
||
}
|
||
} else {
|
||
for (int jj = 0; jj < ord; ++jj) {
|
||
const size_t dst_j = (size_t)(shift - jj) * snx;
|
||
const size_t src_j = (size_t)(shift + jj + 1) * snx;
|
||
for (int k0 = 0; k0 < extc3; ++k0) {
|
||
const size_t kbase = interior_k + (size_t)k0 * splane;
|
||
double *dst = funcc + kbase + dst_j;
|
||
const double *src = funcc + kbase + src_j;
|
||
for (int i = 0; i < nx; ++i) dst[i] = src[i] * s2;
|
||
}
|
||
}
|
||
}
|
||
|
||
/* 4) funcc(:,:,-k) = funcc(:,:,k+1)*SoA(3) */
|
||
const double s3 = SoA[2];
|
||
if (s3 == 1.0) {
|
||
for (int kk = 0; kk < ord; ++kk) {
|
||
const size_t dst_k = (size_t)(shift - kk) * splane;
|
||
const size_t src_k = (size_t)(shift + kk + 1) * splane;
|
||
double *dst = funcc + dst_k;
|
||
const double *src = funcc + src_k;
|
||
for (size_t p = 0; p < splane; ++p) dst[p] = src[p];
|
||
}
|
||
} else if (s3 == -1.0) {
|
||
for (int kk = 0; kk < ord; ++kk) {
|
||
const size_t dst_k = (size_t)(shift - kk) * splane;
|
||
const size_t src_k = (size_t)(shift + kk + 1) * splane;
|
||
double *dst = funcc + dst_k;
|
||
const double *src = funcc + src_k;
|
||
for (size_t p = 0; p < splane; ++p) dst[p] = -src[p];
|
||
}
|
||
} else {
|
||
for (int kk = 0; kk < ord; ++kk) {
|
||
const size_t dst_k = (size_t)(shift - kk) * splane;
|
||
const size_t src_k = (size_t)(shift + kk + 1) * splane;
|
||
double *dst = funcc + dst_k;
|
||
const double *src = funcc + src_k;
|
||
for (size_t p = 0; p < splane; ++p) dst[p] = src[p] * s3;
|
||
}
|
||
}
|
||
}
|
||
|
||
static inline void symmetry_bd(int ord,
|
||
const int extc[3],
|
||
const double *func,
|
||
double *funcc,
|
||
const double SoA[3])
|
||
{
|
||
if (ord <= 0) return;
|
||
|
||
/* Fast paths used by current C kernels: ord=2 (derivs), ord=3 (lopsided/KO). */
|
||
if (ord == 2) {
|
||
symmetry_bd_impl(2, 1, extc, func, funcc, SoA);
|
||
return;
|
||
}
|
||
if (ord == 3) {
|
||
symmetry_bd_impl(3, 2, extc, func, funcc, SoA);
|
||
return;
|
||
}
|
||
|
||
symmetry_bd_impl(ord, ord - 1, extc, func, funcc, SoA);
|
||
}
|
||
#endif
|