249 lines
12 KiB
C
249 lines
12 KiB
C
#include "tool.h"
|
|
|
|
/*
|
|
* Combined advection (lopsided) + KO dissipation (kodis).
|
|
* Uses one shared symmetry_bd buffer per call.
|
|
*/
|
|
void lopsided_kodis(const int ex[3],
|
|
const double *X, const double *Y, const double *Z,
|
|
const double *f, double *f_rhs,
|
|
const double *Sfx, const double *Sfy, const double *Sfz,
|
|
int Symmetry, const double SoA[3], double eps)
|
|
{
|
|
const double ZEO = 0.0, ONE = 1.0, F3 = 3.0;
|
|
const double F6 = 6.0, F18 = 18.0;
|
|
const double F12 = 12.0, F10 = 10.0, EIT = 8.0;
|
|
const double SIX = 6.0, FIT = 15.0, TWT = 20.0;
|
|
const double cof = 64.0; // 2^6
|
|
|
|
const int NO_SYMM = 0, EQ_SYMM = 1;
|
|
|
|
const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
|
|
|
|
const double dX = X[1] - X[0];
|
|
const double dY = Y[1] - Y[0];
|
|
const double dZ = Z[1] - Z[0];
|
|
|
|
const double d12dx = ONE / F12 / dX;
|
|
const double d12dy = ONE / F12 / dY;
|
|
const double d12dz = ONE / F12 / dZ;
|
|
|
|
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 = -2;
|
|
if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -2;
|
|
if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -2;
|
|
|
|
// fh for Fortran-style domain (-2:ex1,-2:ex2,-2:ex3)
|
|
const size_t nx = (size_t)ex1 + 3;
|
|
const size_t ny = (size_t)ex2 + 3;
|
|
const size_t nz = (size_t)ex3 + 3;
|
|
const size_t fh_size = nx * ny * nz;
|
|
|
|
double *fh = (double*)malloc(fh_size * sizeof(double));
|
|
if (!fh) return;
|
|
|
|
symmetry_bd(3, ex, f, fh, SoA);
|
|
|
|
// Advection (same stencil logic as lopsided_c.C)
|
|
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 - 4) {
|
|
f_rhs[p] += sfx * d12dx *
|
|
(-F3 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
|
|
-F10 * fh[idx_fh_F(iF , jF, kF, ex)]
|
|
+F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
|
|
-F6 * fh[idx_fh_F(iF + 2, jF, kF, ex)]
|
|
+ fh[idx_fh_F(iF + 3, jF, kF, ex)]);
|
|
} else if (i0 <= ex1 - 3) {
|
|
f_rhs[p] += sfx * d12dx *
|
|
( fh[idx_fh_F(iF - 2, jF, kF, ex)]
|
|
-EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
|
|
+EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
|
|
- fh[idx_fh_F(iF + 2, jF, kF, ex)]);
|
|
} else if (i0 <= ex1 - 2) {
|
|
f_rhs[p] -= sfx * d12dx *
|
|
(-F3 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
|
|
-F10 * fh[idx_fh_F(iF , jF, kF, ex)]
|
|
+F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
|
|
-F6 * fh[idx_fh_F(iF - 2, jF, kF, ex)]
|
|
+ fh[idx_fh_F(iF - 3, jF, kF, ex)]);
|
|
}
|
|
} else if (sfx < ZEO) {
|
|
if ((i0 - 2) >= iminF) {
|
|
f_rhs[p] -= sfx * d12dx *
|
|
(-F3 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
|
|
-F10 * fh[idx_fh_F(iF , jF, kF, ex)]
|
|
+F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
|
|
-F6 * fh[idx_fh_F(iF - 2, jF, kF, ex)]
|
|
+ fh[idx_fh_F(iF - 3, jF, kF, ex)]);
|
|
} else if ((i0 - 1) >= iminF) {
|
|
f_rhs[p] += sfx * d12dx *
|
|
( fh[idx_fh_F(iF - 2, jF, kF, ex)]
|
|
-EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
|
|
+EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
|
|
- fh[idx_fh_F(iF + 2, jF, kF, ex)]);
|
|
} else if (i0 >= iminF) {
|
|
f_rhs[p] += sfx * d12dx *
|
|
(-F3 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
|
|
-F10 * fh[idx_fh_F(iF , jF, kF, ex)]
|
|
+F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
|
|
-F6 * fh[idx_fh_F(iF + 2, jF, kF, ex)]
|
|
+ fh[idx_fh_F(iF + 3, jF, kF, ex)]);
|
|
}
|
|
}
|
|
|
|
const double sfy = Sfy[p];
|
|
if (sfy > ZEO) {
|
|
if (j0 <= ex2 - 4) {
|
|
f_rhs[p] += sfy * d12dy *
|
|
(-F3 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
|
|
-F10 * fh[idx_fh_F(iF, jF , kF, ex)]
|
|
+F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
|
|
-F6 * fh[idx_fh_F(iF, jF + 2, kF, ex)]
|
|
+ fh[idx_fh_F(iF, jF + 3, kF, ex)]);
|
|
} else if (j0 <= ex2 - 3) {
|
|
f_rhs[p] += sfy * d12dy *
|
|
( fh[idx_fh_F(iF, jF - 2, kF, ex)]
|
|
-EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
|
|
+EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
|
|
- fh[idx_fh_F(iF, jF + 2, kF, ex)]);
|
|
} else if (j0 <= ex2 - 2) {
|
|
f_rhs[p] -= sfy * d12dy *
|
|
(-F3 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
|
|
-F10 * fh[idx_fh_F(iF, jF , kF, ex)]
|
|
+F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
|
|
-F6 * fh[idx_fh_F(iF, jF - 2, kF, ex)]
|
|
+ fh[idx_fh_F(iF, jF - 3, kF, ex)]);
|
|
}
|
|
} else if (sfy < ZEO) {
|
|
if ((j0 - 2) >= jminF) {
|
|
f_rhs[p] -= sfy * d12dy *
|
|
(-F3 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
|
|
-F10 * fh[idx_fh_F(iF, jF , kF, ex)]
|
|
+F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
|
|
-F6 * fh[idx_fh_F(iF, jF - 2, kF, ex)]
|
|
+ fh[idx_fh_F(iF, jF - 3, kF, ex)]);
|
|
} else if ((j0 - 1) >= jminF) {
|
|
f_rhs[p] += sfy * d12dy *
|
|
( fh[idx_fh_F(iF, jF - 2, kF, ex)]
|
|
-EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
|
|
+EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
|
|
- fh[idx_fh_F(iF, jF + 2, kF, ex)]);
|
|
} else if (j0 >= jminF) {
|
|
f_rhs[p] += sfy * d12dy *
|
|
(-F3 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
|
|
-F10 * fh[idx_fh_F(iF, jF , kF, ex)]
|
|
+F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
|
|
-F6 * fh[idx_fh_F(iF, jF + 2, kF, ex)]
|
|
+ fh[idx_fh_F(iF, jF + 3, kF, ex)]);
|
|
}
|
|
}
|
|
|
|
const double sfz = Sfz[p];
|
|
if (sfz > ZEO) {
|
|
if (k0 <= ex3 - 4) {
|
|
f_rhs[p] += sfz * d12dz *
|
|
(-F3 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
|
|
-F10 * fh[idx_fh_F(iF, jF, kF , ex)]
|
|
+F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
|
|
-F6 * fh[idx_fh_F(iF, jF, kF + 2, ex)]
|
|
+ fh[idx_fh_F(iF, jF, kF + 3, ex)]);
|
|
} else if (k0 <= ex3 - 3) {
|
|
f_rhs[p] += sfz * d12dz *
|
|
( fh[idx_fh_F(iF, jF, kF - 2, ex)]
|
|
-EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
|
|
+EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
|
|
- fh[idx_fh_F(iF, jF, kF + 2, ex)]);
|
|
} else if (k0 <= ex3 - 2) {
|
|
f_rhs[p] -= sfz * d12dz *
|
|
(-F3 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
|
|
-F10 * fh[idx_fh_F(iF, jF, kF , ex)]
|
|
+F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
|
|
-F6 * fh[idx_fh_F(iF, jF, kF - 2, ex)]
|
|
+ fh[idx_fh_F(iF, jF, kF - 3, ex)]);
|
|
}
|
|
} else if (sfz < ZEO) {
|
|
if ((k0 - 2) >= kminF) {
|
|
f_rhs[p] -= sfz * d12dz *
|
|
(-F3 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
|
|
-F10 * fh[idx_fh_F(iF, jF, kF , ex)]
|
|
+F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
|
|
-F6 * fh[idx_fh_F(iF, jF, kF - 2, ex)]
|
|
+ fh[idx_fh_F(iF, jF, kF - 3, ex)]);
|
|
} else if ((k0 - 1) >= kminF) {
|
|
f_rhs[p] += sfz * d12dz *
|
|
( fh[idx_fh_F(iF, jF, kF - 2, ex)]
|
|
-EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
|
|
+EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
|
|
- fh[idx_fh_F(iF, jF, kF + 2, ex)]);
|
|
} else if (k0 >= kminF) {
|
|
f_rhs[p] += sfz * d12dz *
|
|
(-F3 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
|
|
-F10 * fh[idx_fh_F(iF, jF, kF , ex)]
|
|
+F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
|
|
-F6 * fh[idx_fh_F(iF, jF, kF + 2, ex)]
|
|
+ fh[idx_fh_F(iF, jF, kF + 3, ex)]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// KO dissipation (same domain restriction as kodiss_c.C)
|
|
if (eps > ZEO) {
|
|
const int i0_lo = (iminF + 2 > 0) ? iminF + 2 : 0;
|
|
const int j0_lo = (jminF + 2 > 0) ? jminF + 2 : 0;
|
|
const int k0_lo = (kminF + 2 > 0) ? kminF + 2 : 0;
|
|
const int i0_hi = imaxF - 4; // inclusive
|
|
const int j0_hi = jmaxF - 4;
|
|
const int k0_hi = kmaxF - 4;
|
|
|
|
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_term =
|
|
((fh[idx_fh_F(iF - 3, jF, kF, ex)] + fh[idx_fh_F(iF + 3, jF, kF, ex)]) -
|
|
SIX * (fh[idx_fh_F(iF - 2, jF, kF, ex)] + fh[idx_fh_F(iF + 2, jF, kF, ex)]) +
|
|
FIT * (fh[idx_fh_F(iF - 1, jF, kF, ex)] + fh[idx_fh_F(iF + 1, jF, kF, ex)]) -
|
|
TWT * fh[idx_fh_F(iF, jF, kF, ex)]) / dX;
|
|
|
|
const double Dy_term =
|
|
((fh[idx_fh_F(iF, jF - 3, kF, ex)] + fh[idx_fh_F(iF, jF + 3, kF, ex)]) -
|
|
SIX * (fh[idx_fh_F(iF, jF - 2, kF, ex)] + fh[idx_fh_F(iF, jF + 2, kF, ex)]) +
|
|
FIT * (fh[idx_fh_F(iF, jF - 1, kF, ex)] + fh[idx_fh_F(iF, jF + 1, kF, ex)]) -
|
|
TWT * fh[idx_fh_F(iF, jF, kF, ex)]) / dY;
|
|
|
|
const double Dz_term =
|
|
((fh[idx_fh_F(iF, jF, kF - 3, ex)] + fh[idx_fh_F(iF, jF, kF + 3, ex)]) -
|
|
SIX * (fh[idx_fh_F(iF, jF, kF - 2, ex)] + fh[idx_fh_F(iF, jF, kF + 2, ex)]) +
|
|
FIT * (fh[idx_fh_F(iF, jF, kF - 1, ex)] + fh[idx_fh_F(iF, jF, kF + 1, ex)]) -
|
|
TWT * fh[idx_fh_F(iF, jF, kF, ex)]) / dZ;
|
|
|
|
f_rhs[p] += (eps / cof) * (Dx_term + Dy_term + Dz_term);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
free(fh);
|
|
}
|