Omar SobhandClaude Fable 5.1
77ecc54b24
CI / Test (macos-latest) (push) Blocked by required conditions
CI / Test (ubuntu-latest) (push) Blocked by required conditions
CI / Python Bindings (maturin) (macos-latest) (push) Blocked by required conditions
CI / Python Bindings (maturin) (ubuntu-latest) (push) Blocked by required conditions
CI / WASM Build + Size Check (push) Blocked by required conditions
CI / Distributed Training Tests (push) Blocked by required conditions
CI / CI Success (push) Blocked by required conditions
CI / Build (macos-latest) (push) Waiting to run
CI / Build CPU-Only (Explicit) (push) Failing after 4s
Documentation / Build API Documentation (push) Failing after 4s
Documentation / Build User Guide (push) Successful in 5s
CI / Format Check (push) Failing after 12s
CI / Clippy Check (push) Failing after 37s
CI / Build (ubuntu-latest) (push) Failing after 2m48s
Performance Benchmarks / Run Benchmarks (push) Successful in 4m4s
Co-Authored-By: Claude Fable 5.1 <[email protected]>
491 lines
25 KiB
Plaintext
491 lines
25 KiB
Plaintext
/**
|
||
* embedded3 item 9b: the apertured cut-cell wall's maps on the device — the
|
||
* host `step/cut_predictor.rs` and `cutwall.rs` expression for expression
|
||
* (FMA contraction off): the predictor on the open faces of one component
|
||
* (the momentum volume `V_u = α h A`, mass fluxes averaged from the two
|
||
* cells, upwind + TVD, apertured diffusion, the projection's pressure
|
||
* force, the net mass flux times the face's value, the implicit wall
|
||
* shear, the inertia floor), the apertured continuity with the wall flux
|
||
* and the fold of the virtually merged small cells into their masters,
|
||
* the corrections on the open faces reading the owner's p', the pressure
|
||
* update. Static bodies (the tables are built once).
|
||
*
|
||
* Appended after `e3_step.cu` at load (shares E3Params / E3Ptrs and the
|
||
* helpers there).
|
||
*/
|
||
#define CUT_INERTIA_FLOOR 0.1
|
||
#define CUT_DISTANCE_FLOOR 0.05
|
||
#define CUT_DISTANCE_FLOOR_FINE 0.01
|
||
#define CUT_TRANSVERSE_FLOOR 0.2
|
||
|
||
struct E3Cut {
|
||
const double *a_u, *a_v, *a_w; /* apertures per face */
|
||
const double *d_u, *d_v, *d_w; /* φ at the face centres */
|
||
const double *ub_u, *ub_v, *ub_w; /* surface velocity component at the foot per face */
|
||
const double *wall_flux; /* per cell, compatible */
|
||
const int *open_u, *open_v, *open_w;/* unknown faces */
|
||
const int *active; /* cells with an equation */
|
||
const unsigned int *owner; /* the master of a merged cell (itself otherwise) */
|
||
const unsigned int *fold_ptr, *fold_idx; /* CSR: the slaves of every cell */
|
||
double *cell_flux; /* scratch per cell */
|
||
const double *s_u, *s_v, *s_w; /* open-part centroid shifts per face, 3 interleaved (S2-5) */
|
||
const double *vn_u, *vn_v, *vn_w; /* the wall's normal velocity into the fluid per face (A3-i) */
|
||
};
|
||
|
||
/* f(xi) = xi / (1 - exp(-xi)), f(0) = 1 exactly (closure.rs advancing_factor). */
|
||
__device__ __forceinline__ double adv_factor(double xi)
|
||
{
|
||
return fabs(xi) < 1e-8 ? 1.0 + 0.5 * xi : xi / (-expm1(-xi));
|
||
}
|
||
|
||
/* Face index of component c at lattice (i, j, k); −1 outside (z wraps when periodic). */
|
||
__device__ __forceinline__ int cut_face(const E3Params& g, int c, int i, int j, int k)
|
||
{
|
||
int ni = c == 0 ? g.nx + 1 : g.nx;
|
||
int nj = c == 1 ? g.ny + 1 : g.ny;
|
||
int nk = c == 2 ? g.nz + 1 : g.nz;
|
||
if (i < 0 || i >= ni || j < 0 || j >= nj) return -1;
|
||
if (g.periodic_z) { k = ((k % g.nz) + g.nz) % g.nz; }
|
||
else if (k < 0 || k >= nk) return -1;
|
||
if (c == 0) return (k * g.ny + j) * (g.nx + 1) + i;
|
||
if (c == 1) return (k * (g.ny + 1) + j) * g.nx + i;
|
||
return (k * g.ny + j) * g.nx + i;
|
||
}
|
||
|
||
__device__ __forceinline__ int cut_cell(const E3Params& g, int i, int j, int k)
|
||
{
|
||
if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return -1;
|
||
if (g.periodic_z) { k = ((k % g.nz) + g.nz) % g.nz; }
|
||
else if (k < 0 || k >= g.nz) return -1;
|
||
return (k * g.ny + j) * g.nx + i;
|
||
}
|
||
|
||
__device__ __forceinline__ const double* cut_ap(const E3Cut& m, int c) { return c == 0 ? m.a_u : (c == 1 ? m.a_v : m.a_w); }
|
||
__device__ __forceinline__ const double* cut_old(const E3Ptrs& f, int c) { return c == 0 ? f.uo : (c == 1 ? f.vo : f.wo); }
|
||
|
||
/* Aperture of component `cc` at the lattice point, or `dflt` outside. */
|
||
__device__ __forceinline__ double cut_ap_at(const E3Params& g, const E3Cut& m, int cc, int i, int j, int k, double dflt)
|
||
{
|
||
int fidx = cut_face(g, cc, i, j, k);
|
||
return fidx < 0 ? dflt : cut_ap(m, cc)[fidx];
|
||
}
|
||
|
||
/* The control volume's geometry of the face of component c at (i, j, k)
|
||
(cutwall.rs cv_geometry): side apertures, the closing wall vector, the
|
||
floored wall distance. */
|
||
__device__ void cut_cv(const E3Params& g, const E3Cut& m, int c, int i, int j, int k, int fidx,
|
||
double* alpha_out, double* apm, double* app, double* wall, double* distance_out)
|
||
{
|
||
double h[3] = { g.dx, g.dy, g.dz };
|
||
double area[3] = { g.dy * g.dz, g.dx * g.dz, g.dx * g.dy };
|
||
double alpha = cut_ap(m, c)[fidx];
|
||
int cm[3] = { i, j, k }; cm[c] -= 1; /* cell minus */
|
||
int cp[3] = { i, j, k }; /* cell plus */
|
||
for (int d = 0; d < 3; ++d) {
|
||
if (d == c) {
|
||
int qm[3] = { i, j, k }; qm[c] -= 1;
|
||
int qp[3] = { i, j, k }; qp[c] += 1;
|
||
double am = cut_ap_at(g, m, c, qm[0], qm[1], qm[2], alpha);
|
||
double apl = cut_ap_at(g, m, c, qp[0], qp[1], qp[2], alpha);
|
||
apm[d] = 0.5 * (am + alpha);
|
||
app[d] = 0.5 * (alpha + apl);
|
||
} else {
|
||
int q1[3] = { cm[0], cm[1], cm[2] }; q1[d] += 1;
|
||
int q2[3] = { cp[0], cp[1], cp[2] }; q2[d] += 1;
|
||
apm[d] = 0.5 * (cut_ap_at(g, m, d, cm[0], cm[1], cm[2], 1.0) + cut_ap_at(g, m, d, cp[0], cp[1], cp[2], 1.0));
|
||
app[d] = 0.5 * (cut_ap_at(g, m, d, q1[0], q1[1], q1[2], 1.0) + cut_ap_at(g, m, d, q2[0], q2[1], q2[2], 1.0));
|
||
}
|
||
}
|
||
for (int d = 0; d < 3; ++d) wall[d] = -(app[d] - apm[d]) * area[d];
|
||
double h_min = fmin(fmin(g.dx, g.dy), g.dz);
|
||
const double* dist = c == 0 ? m.d_u : (c == 1 ? m.d_v : m.d_w);
|
||
/* bit 4 of wall_order: the oblique wall distance (the in-plane part of the wall normal) */
|
||
double n_t = 1.0;
|
||
if (g.wall_order & 16) {
|
||
double a_w0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
|
||
if (a_w0 > 0.0) { double n_c = wall[c] / a_w0; n_t = sqrt(fmax(1.0 - n_c * n_c, 0.0)); }
|
||
}
|
||
/* bit 8 of wall_order: the fine distance floor (S2-6) */
|
||
double floor_h = (g.wall_order & 256) ? CUT_DISTANCE_FLOOR_FINE : CUT_DISTANCE_FLOOR;
|
||
*alpha_out = alpha;
|
||
*distance_out = fmax(dist[fidx] + 0.5 * h[c] * (1.0 - alpha) * n_t, floor_h * h_min);
|
||
}
|
||
|
||
/* S2-6 (bit 7 of wall_order): the transverse part of the centroid diffusion
|
||
between the face (fidx, its geometry given) and its open neighbour fq at
|
||
q along d — grad(u) . ds_perp with the cut faces' own wall-normal
|
||
gradients (u - U_b)/d_f (cut_predictor.rs `transverse`). */
|
||
__device__ double cut_transverse(const E3Params& g, const E3Cut& m, int c, int d, int fidx, int fq, const int* q,
|
||
double sign, double u0, double ub, double alpha, const double* wall, double distance,
|
||
double uq, const double* sh, const double* ubt)
|
||
{
|
||
double ds[3] = { 0.0, 0.0, 0.0 };
|
||
int any = 0;
|
||
for (int e = 0; e < 3; ++e) {
|
||
if (e == d) continue;
|
||
ds[e] = sign * (sh[3 * fq + e] - sh[3 * fidx + e]);
|
||
if (ds[e] != 0.0) any = 1;
|
||
}
|
||
if (!any) return 0.0;
|
||
double alpha_q, apm_q[3], app_q[3], wall_q[3], distance_q;
|
||
cut_cv(g, m, c, q[0], q[1], q[2], fq, &alpha_q, apm_q, app_q, wall_q, &distance_q);
|
||
/* Explicit, coefficient ~ 1/d_f: the gradient from the faces at least
|
||
CUT_TRANSVERSE_FLOOR h off the wall — a FULL neighbour too, over its own
|
||
distance along the cut partner's normal; only when neither is that far,
|
||
from the cut faces over the floored distance (cut_predictor.rs). */
|
||
double h_min = fmin(fmin(g.dx, g.dy), g.dz);
|
||
double d_min = CUT_TRANSVERSE_FLOOR * h_min;
|
||
double a0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
|
||
double aq = sqrt(wall_q[0] * wall_q[0] + wall_q[1] * wall_q[1] + wall_q[2] * wall_q[2]);
|
||
int cut0 = a0 > 0.0 && alpha < 1.0, cutq = aq > 0.0 && alpha_q < 1.0;
|
||
if (!cut0 && !cutq) return 0.0;
|
||
double n0[3], nq[3];
|
||
for (int e = 0; e < 3; ++e) { n0[e] = cut0 ? wall[e] / a0 : 0.0; nq[e] = cutq ? wall_q[e] / aq : 0.0; }
|
||
if (!cut0) for (int e = 0; e < 3; ++e) n0[e] = nq[e];
|
||
if (!cutq) for (int e = 0; e < 3; ++e) nq[e] = n0[e];
|
||
int far = distance >= d_min || distance_q >= d_min;
|
||
double gr[3] = { 0.0, 0.0, 0.0 };
|
||
double count = 0.0;
|
||
if (!((far && distance < d_min) || (!far && !cut0))) {
|
||
double slope = (u0 - ub) / fmax(distance, d_min);
|
||
for (int e = 0; e < 3; ++e) gr[e] -= slope * n0[e];
|
||
count += 1.0;
|
||
}
|
||
if (!((far && distance_q < d_min) || (!far && !cutq))) {
|
||
double slope = (uq - ubt[fq]) / fmax(distance_q, d_min);
|
||
for (int e = 0; e < 3; ++e) gr[e] -= slope * nq[e];
|
||
count += 1.0;
|
||
}
|
||
if (count == 0.0) return 0.0;
|
||
return (gr[0] * ds[0] + gr[1] * ds[1] + gr[2] * ds[2]) / count;
|
||
}
|
||
|
||
/* The predicted value of the open face of component c at (i, j, k). */
|
||
__device__ double cut_face_update(const E3Params& g, const E3Ptrs& f, const E3Cut& m, int c, int i, int j, int k, int fidx)
|
||
{
|
||
double h[3] = { g.dx, g.dy, g.dz };
|
||
double n[3] = { (double)g.nx, (double)g.ny, (double)g.nz };
|
||
double area[3] = { g.dy * g.dz, g.dx * g.dz, g.dx * g.dy };
|
||
double rho = g.rho, mu = g.nu * g.rho;
|
||
int sides[3][2] = { { g.bx0, g.bx1 }, { g.by0, g.by1 }, { g.bz0, g.bz1 } };
|
||
int p[3] = { i, j, k };
|
||
const double* old_c = cut_old(f, c);
|
||
double u0 = old_c[fidx];
|
||
/* the control volume's geometry (cutwall.rs cv_geometry) */
|
||
double alpha, apm[3], app[3], wall[3], distance;
|
||
cut_cv(g, m, c, i, j, k, fidx, &alpha, apm, app, wall, &distance);
|
||
int cm[3] = { i, j, k }; cm[c] -= 1; /* cell minus */
|
||
int cp[3] = { i, j, k }; /* cell plus */
|
||
const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w);
|
||
double ub = ubt[fidx];
|
||
/* face position */
|
||
double x[3];
|
||
for (int d = 0; d < 3; ++d) x[d] = (p[d] + (d == c ? 0.0 : 0.5)) * h[d];
|
||
|
||
double mass_out = 0.0, conv = 0.0, diff = 0.0, wall_implicit = 0.0, wall_rhs = 0.0;
|
||
for (int d = 0; d < 3; ++d) {
|
||
double a_d = area[d];
|
||
int q[3];
|
||
/* neighbouring faces of this component along d */
|
||
q[0] = i; q[1] = j; q[2] = k; q[d] += 1; int f_up1 = cut_face(g, c, q[0], q[1], q[2]);
|
||
q[0] = i; q[1] = j; q[2] = k; q[d] += 2; int f_up2 = cut_face(g, c, q[0], q[1], q[2]);
|
||
q[0] = i; q[1] = j; q[2] = k; q[d] -= 1; int f_dn1 = cut_face(g, c, q[0], q[1], q[2]);
|
||
q[0] = i; q[1] = j; q[2] = k; q[d] -= 2; int f_dn2 = cut_face(g, c, q[0], q[1], q[2]);
|
||
double up1 = f_up1 >= 0 ? old_c[f_up1] : 0.0;
|
||
double up2 = f_up2 >= 0 ? old_c[f_up2] : 0.0;
|
||
double dn1 = f_dn1 >= 0 ? old_c[f_dn1] : 0.0;
|
||
double dn2 = f_dn2 >= 0 ? old_c[f_dn2] : 0.0;
|
||
double m_plus, m_minus;
|
||
if (d == c) {
|
||
double f_up = (f_up1 >= 0 ? cut_ap(m, c)[f_up1] : alpha) * (f_up1 >= 0 ? up1 : u0);
|
||
double f_dn = (f_dn1 >= 0 ? cut_ap(m, c)[f_dn1] : alpha) * (f_dn1 >= 0 ? dn1 : u0);
|
||
double f0 = alpha * u0;
|
||
m_plus = 0.5 * (f0 + f_up) * a_d;
|
||
m_minus = 0.5 * (f_dn + f0) * a_d;
|
||
} else {
|
||
const double* old_d = cut_old(f, d);
|
||
int q1[3] = { cm[0], cm[1], cm[2] }; q1[d] += 1;
|
||
int q2[3] = { cp[0], cp[1], cp[2] }; q2[d] += 1;
|
||
int fa = cut_face(g, d, q1[0], q1[1], q1[2]);
|
||
int fb = cut_face(g, d, q2[0], q2[1], q2[2]);
|
||
int fc = cut_face(g, d, cm[0], cm[1], cm[2]);
|
||
int fd = cut_face(g, d, cp[0], cp[1], cp[2]);
|
||
double fla = fa >= 0 ? cut_ap(m, d)[fa] * old_d[fa] : 0.0;
|
||
double flb = fb >= 0 ? cut_ap(m, d)[fb] * old_d[fb] : 0.0;
|
||
double flc = fc >= 0 ? cut_ap(m, d)[fc] * old_d[fc] : 0.0;
|
||
double fld = fd >= 0 ? cut_ap(m, d)[fd] * old_d[fd] : 0.0;
|
||
m_plus = 0.5 * (fla + flb) * a_d;
|
||
m_minus = 0.5 * (flc + fld) * a_d;
|
||
}
|
||
mass_out += m_plus - m_minus;
|
||
/* beyond a domain side along d */
|
||
double beyond_p = u0, beyond_m = u0;
|
||
if (sides[d][1] == SIDE_VELOCITY || sides[d][0] == SIDE_VELOCITY) {
|
||
/* the boundary tables: [side][component] at the face's transverse position */
|
||
const double* tp = 0; const double* tm = 0; int idx_p = 0, idx_m = 0;
|
||
if (d == 0) { tp = c == 0 ? f.bx1u : (c == 1 ? f.bx1v : f.bx1w); tm = c == 0 ? f.bx0u : (c == 1 ? f.bx0v : f.bx0w);
|
||
int nj = c == 1 ? g.ny + 1 : g.ny; idx_p = idx_m = k * nj + j; }
|
||
else if (d == 1) { tp = c == 0 ? f.by1u : (c == 1 ? f.by1v : f.by1w); tm = c == 0 ? f.by0u : (c == 1 ? f.by0v : f.by0w);
|
||
int ni = c == 0 ? g.nx + 1 : g.nx; idx_p = idx_m = k * ni + i; }
|
||
else { tp = c == 0 ? f.bz1u : (c == 1 ? f.bz1v : f.bz1w); tm = c == 0 ? f.bz0u : (c == 1 ? f.bz0v : f.bz0w);
|
||
int ni = c == 0 ? g.nx + 1 : g.nx; idx_p = idx_m = j * ni + i; }
|
||
if (sides[d][1] == SIDE_VELOCITY) beyond_p = tp[idx_p];
|
||
if (sides[d][0] == SIDE_VELOCITY) beyond_m = tm[idx_m];
|
||
}
|
||
(void)n; (void)x;
|
||
/* convection through the plus face */
|
||
double u_plus, delta_plus;
|
||
/* bit 10 of wall_order (A2): a prescribed neighbour face exchanges no convective momentum */
|
||
const int* opn = c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w);
|
||
int presc_up = (g.wall_order & 1024) && f_up1 >= 0 && !opn[f_up1];
|
||
int presc_dn = (g.wall_order & 1024) && f_dn1 >= 0 && !opn[f_dn1];
|
||
if (presc_up) { u_plus = u0; delta_plus = 0.0; }
|
||
else if (f_up1 >= 0) {
|
||
if (g.scheme == SCHEME_UPWIND) delta_plus = 0.0;
|
||
else if (m_plus >= 0.0) delta_plus = face_corr3(g.scheme, f_dn1 >= 0, dn1, u0, up1);
|
||
else delta_plus = face_corr3(g.scheme, f_up2 >= 0, up2, up1, u0);
|
||
u_plus = upwind3(m_plus, u0, up1);
|
||
} else { u_plus = upwind3(m_plus, u0, beyond_p); delta_plus = 0.0; }
|
||
double u_minus, delta_minus;
|
||
if (presc_dn) { u_minus = u0; delta_minus = 0.0; }
|
||
else if (f_dn1 >= 0) {
|
||
if (g.scheme == SCHEME_UPWIND) delta_minus = 0.0;
|
||
else if (m_minus >= 0.0) delta_minus = face_corr3(g.scheme, f_dn2 >= 0, dn2, dn1, u0);
|
||
else delta_minus = face_corr3(g.scheme, f_up1 >= 0, up1, u0, dn1);
|
||
u_minus = upwind3(m_minus, dn1, u0);
|
||
} else { u_minus = upwind3(m_minus, beyond_m, u0); delta_minus = 0.0; }
|
||
conv += m_plus * (u_plus + delta_plus) - m_minus * (u_minus + delta_minus);
|
||
/* diffusion */
|
||
double g_minus = apm[d], g_plus = app[d];
|
||
/* bit 5 of wall_order: the exchange with a SOLID neighbour over the axis
|
||
distance to the wall, implicit (S2-5) */
|
||
double delta_x = h[d];
|
||
if (g.wall_order & 32) {
|
||
double a_w0 = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
|
||
if (a_w0 > 0.0) {
|
||
double n_d = fabs(wall[d]) / a_w0;
|
||
if (n_d >= 1e-12) delta_x = fmin(distance / n_d, h[d]);
|
||
}
|
||
}
|
||
int solid_up = (g.wall_order & 32) && f_up1 >= 0 && cut_ap(m, c)[f_up1] == 0.0;
|
||
int solid_dn = (g.wall_order & 32) && f_dn1 >= 0 && cut_ap(m, c)[f_dn1] == 0.0;
|
||
/* bit 6 of wall_order: cross diffusion over the open-part centroid spacing,
|
||
the part beyond 1/h point-implicit (S2-5) */
|
||
int centroid = (g.wall_order & 64) && d != c;
|
||
const double* sh = c == 0 ? m.s_u : (c == 1 ? m.s_v : m.s_w);
|
||
/* bit 7 of wall_order: the transverse centroid correction, cross and own direction (S2-6) */
|
||
if ((g.wall_order & 128) && (g.wall_order & 64)) {
|
||
if (f_up1 >= 0 && !solid_up && cut_ap(m, c)[f_up1] > 0.0) {
|
||
int qn[3] = { i, j, k }; qn[d] += 1;
|
||
double dl = h[d];
|
||
if (d != c) { dl = h[d] + (sh[3 * f_up1 + d] - sh[3 * fidx + d]); dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]); }
|
||
diff -= mu * g_plus * a_d * cut_transverse(g, m, c, d, fidx, f_up1, qn, 1.0, u0, ub, alpha, wall, distance, up1, sh, ubt) / dl;
|
||
}
|
||
if (f_dn1 >= 0 && !solid_dn && cut_ap(m, c)[f_dn1] > 0.0) {
|
||
int qn[3] = { i, j, k }; qn[d] -= 1;
|
||
double dl = h[d];
|
||
if (d != c) { dl = h[d] - (sh[3 * f_dn1 + d] - sh[3 * fidx + d]); dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]); }
|
||
diff += mu * g_minus * a_d * cut_transverse(g, m, c, d, fidx, f_dn1, qn, -1.0, u0, ub, alpha, wall, distance, dn1, sh, ubt) / dl;
|
||
}
|
||
}
|
||
if (solid_up) { double kx = mu * g_plus * a_d / delta_x; wall_implicit += kx; wall_rhs += kx * up1; }
|
||
else if (f_up1 >= 0 && centroid) {
|
||
double dl = h[d] + (sh[3 * f_up1 + d] - sh[3 * fidx + d]);
|
||
dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]);
|
||
double kx = mu * g_plus * a_d * (1.0 / dl - 1.0 / h[d]);
|
||
if (kx > 0.0) { wall_implicit += kx; wall_rhs += kx * up1; diff += mu * g_plus * a_d * (up1 - u0) / h[d]; }
|
||
else diff += (mu * g_plus * a_d / h[d] + kx) * (up1 - u0);
|
||
}
|
||
else if (f_up1 >= 0) diff += mu * g_plus * a_d * (up1 - u0) / h[d];
|
||
else if (sides[d][1] == SIDE_VELOCITY) diff += mu * g_plus * a_d * (beyond_p - u0) / (0.5 * h[d]);
|
||
if (solid_dn) { double kx = mu * g_minus * a_d / delta_x; wall_implicit += kx; wall_rhs += kx * dn1; }
|
||
else if (f_dn1 >= 0 && centroid) {
|
||
double dl = h[d] - (sh[3 * f_dn1 + d] - sh[3 * fidx + d]);
|
||
dl = fmin(fmax(dl, 0.25 * h[d]), 2.0 * h[d]);
|
||
double kx = mu * g_minus * a_d * (1.0 / dl - 1.0 / h[d]);
|
||
if (kx > 0.0) { wall_implicit += kx; wall_rhs += kx * dn1; diff -= mu * g_minus * a_d * (u0 - dn1) / h[d]; }
|
||
else diff -= (mu * g_minus * a_d / h[d] + kx) * (u0 - dn1);
|
||
}
|
||
else if (f_dn1 >= 0) diff -= mu * g_minus * a_d * (u0 - dn1) / h[d];
|
||
else if (sides[d][0] == SIDE_VELOCITY) diff -= mu * g_minus * a_d * (u0 - beyond_m) / (0.5 * h[d]);
|
||
}
|
||
/* The mass fluxes above are volume fluxes: the momentum flux carries rho. */
|
||
conv = rho * (conv - mass_out * u0);
|
||
int cpi = cut_cell(g, cp[0], cp[1], cp[2]);
|
||
int cmi = cut_cell(g, cm[0], cm[1], cm[2]);
|
||
double p_plus = cpi >= 0 ? f.p[cpi] : 0.0;
|
||
double p_minus = cmi >= 0 ? f.p[cmi] : 0.0;
|
||
double pressure = -(p_plus - p_minus) * alpha * area[c];
|
||
double v_u = alpha * h[c] * area[c];
|
||
const double* src = c == 0 ? f.su : (c == 1 ? f.sv : f.sw);
|
||
double source = src[fidx] * v_u;
|
||
double a_w = sqrt(wall[0] * wall[0] + wall[1] * wall[1] + wall[2] * wall[2]);
|
||
/* The wall gradient: one-point (order 1) or quadratic through the next
|
||
open face away from the body along the wall normal's dominant axis
|
||
(order 2; the neighbour's old value explicit). */
|
||
double c1 = 1.0 / distance, shear_explicit = 0.0;
|
||
/* bit 9 of wall_order: the advancing-wall closure on the one-point coefficient (A3-i) */
|
||
if (g.wall_order & 512) {
|
||
const double* vnt = c == 0 ? m.vn_u : (c == 1 ? m.vn_v : m.vn_w);
|
||
c1 = adv_factor(vnt[fidx] * distance / g.nu) / distance;
|
||
}
|
||
if ((g.wall_order & 15) >= 2 && a_w > 0.0) {
|
||
double nw[3] = { wall[0] / a_w, wall[1] / a_w, wall[2] / a_w };
|
||
int d = 0;
|
||
for (int kk = 1; kk < 3; ++kk) if (fabs(nw[kk]) > fabs(nw[d])) d = kk;
|
||
int q[3] = { i, j, k };
|
||
q[d] -= nw[d] > 0.0 ? 1 : -1;
|
||
int fn = cut_face(g, c, q[0], q[1], q[2]);
|
||
if (fn >= 0 && cut_ap(m, c)[fn] > 0.0) {
|
||
double d1 = distance, d2 = d1 + h[d] * fabs(nw[d]);
|
||
c1 = d2 / (d1 * (d2 - d1));
|
||
double c2 = -d1 / (d2 * (d2 - d1));
|
||
shear_explicit = mu * a_w * c2 * (old_c[fn] - ub);
|
||
}
|
||
}
|
||
double shear = mu * a_w * c1;
|
||
double v_eff = fmax(alpha, CUT_INERTIA_FLOOR) * h[c] * area[c];
|
||
double inertia = rho * v_eff / g.dt;
|
||
return (inertia * u0 - conv + diff + pressure + source + shear * ub - shear_explicit + wall_rhs) / (inertia + shear + wall_implicit);
|
||
}
|
||
|
||
/* The predictor on the open interior faces of component c. */
|
||
extern "C" __global__ void e3_cut_predict(E3Params g, E3Ptrs f, E3Cut m, int c)
|
||
{
|
||
int t = blockIdx.x * blockDim.x + threadIdx.x;
|
||
int ni = c == 0 ? g.nx + 1 : g.nx;
|
||
int nj = c == 1 ? g.ny + 1 : g.ny;
|
||
int nk = c == 2 ? g.nz + 1 : g.nz;
|
||
if (t >= ni * nj * nk) return;
|
||
int i = t % ni; int j = (t / ni) % nj; int k = t / (ni * nj);
|
||
if (c == 0 && (i == 0 || i == g.nx)) return;
|
||
if (c == 1 && (j == 0 || j == g.ny)) return;
|
||
if (c == 2) { if (g.periodic_z) { if (k == g.nz) return; } else if (k == 0 || k == g.nz) return; }
|
||
const int* open = c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w);
|
||
if (!open[t]) return;
|
||
double val = cut_face_update(g, f, m, c, i, j, k, t);
|
||
double* out = c == 0 ? f.u : (c == 1 ? f.v : f.w);
|
||
out[t] = val;
|
||
}
|
||
|
||
/* cell_flux = ρ (Σ apertured flux out + wall flux) on the active cells (0 elsewhere). */
|
||
extern "C" __global__ void e3_cut_divergence(E3Params g, E3Ptrs f, E3Cut m, int starred)
|
||
{
|
||
int t = blockIdx.x * blockDim.x + threadIdx.x;
|
||
int total = g.nx * g.ny * g.nz;
|
||
if (t >= total) return;
|
||
if (!m.active[t]) { m.cell_flux[t] = 0.0; return; }
|
||
int i = t % g.nx; int j = (t / g.nx) % g.ny; int k = t / (g.nx * g.ny);
|
||
const double* u = starred ? f.us : f.u;
|
||
const double* v = starred ? f.vs : f.v;
|
||
const double* w = starred ? f.ws : f.w;
|
||
int ue = uf3(g, k, j, i + 1), uw = uf3(g, k, j, i);
|
||
int vn = vf3(g, k, j + 1, i), vs = vf3(g, k, j, i);
|
||
int wt = wf3(g, k + 1, j, i), wb = wf3(g, k, j, i);
|
||
double divergence_flux = g.rho
|
||
* ((m.a_u[ue] * u[ue] - m.a_u[uw] * u[uw]) * (g.dy * g.dz)
|
||
+ (m.a_v[vn] * v[vn] - m.a_v[vs] * v[vs]) * (g.dx * g.dz)
|
||
+ (m.a_w[wt] * w[wt] - m.a_w[wb] * w[wb]) * (g.dx * g.dy));
|
||
divergence_flux += g.rho * m.wall_flux[t];
|
||
m.cell_flux[t] = divergence_flux;
|
||
}
|
||
|
||
/* The merged fluxes: sp = −(own + slaves) on the owning cells when `write_sp`;
|
||
* partial[block] = Σ |merged flux| (the source scale / the mass imbalance). */
|
||
extern "C" __global__ void e3_cut_fold(E3Params g, E3Ptrs f, E3Cut m, int write_sp, double* __restrict__ partial)
|
||
{
|
||
int t = blockIdx.x * blockDim.x + threadIdx.x;
|
||
int total = g.nx * g.ny * g.nz;
|
||
double v = 0.0;
|
||
if (t < total) {
|
||
if (m.active[t] && m.owner[t] == (unsigned int)t) {
|
||
double s = m.cell_flux[t];
|
||
for (unsigned int l = m.fold_ptr[t]; l < m.fold_ptr[t + 1]; ++l) s += m.cell_flux[m.fold_idx[l]];
|
||
if (write_sp) f.sp[t] = -s;
|
||
v = fabs(s);
|
||
} else if (write_sp) {
|
||
f.sp[t] = 0.0;
|
||
}
|
||
}
|
||
__shared__ double sh[256];
|
||
sh[threadIdx.x] = v;
|
||
__syncthreads();
|
||
for (int s = 128; s > 0; s >>= 1) { if (threadIdx.x < s) sh[threadIdx.x] += sh[threadIdx.x + s]; __syncthreads(); }
|
||
if (threadIdx.x == 0) partial[blockIdx.x] = sh[0];
|
||
}
|
||
|
||
/* The corrections on the open interior faces of component c from the owners' p'. */
|
||
extern "C" __global__ void e3_cut_correct(E3Params g, E3Ptrs f, E3Cut m, int c)
|
||
{
|
||
int t = blockIdx.x * blockDim.x + threadIdx.x;
|
||
int ni = c == 0 ? g.nx + 1 : g.nx;
|
||
int nj = c == 1 ? g.ny + 1 : g.ny;
|
||
int nk = c == 2 ? g.nz + 1 : g.nz;
|
||
if (t >= ni * nj * nk) return;
|
||
int i = t % ni; int j = (t / ni) % nj; int k = t / (ni * nj);
|
||
double cc = g.dt / g.rho;
|
||
const int* open = c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w);
|
||
double* out = c == 0 ? f.u : (c == 1 ? f.v : f.w);
|
||
const double* star = c == 0 ? f.us : (c == 1 ? f.vs : f.ws);
|
||
int interior = c == 0 ? (i >= 1 && i < g.nx) : (c == 1 ? (j >= 1 && j < g.ny)
|
||
: (g.periodic_z ? (k < g.nz) : (k >= 1 && k < g.nz)));
|
||
if (interior) {
|
||
if (!open[t]) return;
|
||
int cp[3] = { i, j, k }, cm[3] = { i, j, k }; cm[c] -= 1;
|
||
int a = cut_cell(g, cp[0], cp[1], cp[2]);
|
||
int b = cut_cell(g, cm[0], cm[1], cm[2]);
|
||
double h = c == 0 ? g.dx : (c == 1 ? g.dy : g.dz);
|
||
double dp = (f.pp[m.owner[a]] - f.pp[m.owner[b]]) / h;
|
||
out[t] = star[t] - cc * dp;
|
||
return;
|
||
}
|
||
/* outlet faces against 0 outside (no body reaches an outlet) */
|
||
int s0 = c == 0 ? g.bx0 : (c == 1 ? g.by0 : g.bz0);
|
||
int s1 = c == 0 ? g.bx1 : (c == 1 ? g.by1 : g.bz1);
|
||
int own = c == 0 ? i : (c == 1 ? j : k);
|
||
int nn = c == 0 ? g.nx : (c == 1 ? g.ny : g.nz);
|
||
double h = c == 0 ? g.dx : (c == 1 ? g.dy : g.dz);
|
||
if (own == 0 && s0 == SIDE_OUTLET) {
|
||
int cp[3] = { i, j, k };
|
||
int a = cut_cell(g, cp[0], cp[1], cp[2]);
|
||
out[t] = star[t] - cc * (f.pp[a] - 0.0) / (0.5 * h);
|
||
} else if (own == nn && s1 == SIDE_OUTLET) {
|
||
int cm[3] = { i, j, k }; cm[c] -= 1;
|
||
int b = cut_cell(g, cm[0], cm[1], cm[2]);
|
||
out[t] = star[t] - cc * (0.0 - f.pp[b]) / (0.5 * h);
|
||
}
|
||
}
|
||
|
||
/* p += p'(owner) on the active cells. */
|
||
extern "C" __global__ void e3_cut_add_p(E3Params g, E3Ptrs f, E3Cut m)
|
||
{
|
||
int t = blockIdx.x * blockDim.x + threadIdx.x;
|
||
if (t >= g.nx * g.ny * g.nz) return;
|
||
if (!m.active[t]) return;
|
||
f.p[t] += f.pp[m.owner[t]];
|
||
}
|
||
|
||
/* The prescribed interior faces of component c take the surface velocity
|
||
* (the host `impose` for a cut mask; `open` = the instantaneous kinds). */
|
||
extern "C" __global__ void e3_cut_impose(E3Params g, E3Ptrs f, E3Cut m, int c)
|
||
{
|
||
int t = blockIdx.x * blockDim.x + threadIdx.x;
|
||
int ni = c == 0 ? g.nx + 1 : g.nx;
|
||
int nj = c == 1 ? g.ny + 1 : g.ny;
|
||
int nk = c == 2 ? g.nz + 1 : g.nz;
|
||
if (t >= ni * nj * nk) return;
|
||
int i = t % ni; int j = (t / ni) % nj; int k = t / (ni * nj);
|
||
if (c == 0 && (i == 0 || i == g.nx)) return;
|
||
if (c == 1 && (j == 0 || j == g.ny)) return;
|
||
if (c == 2) { if (g.periodic_z) { if (k == g.nz) return; } else if (k == 0 || k == g.nz) return; }
|
||
const int* open = c == 0 ? m.open_u : (c == 1 ? m.open_v : m.open_w);
|
||
if (open[t]) return;
|
||
/* only within the imposition band (host `impose`: 4 cells) */
|
||
const double* dist = c == 0 ? m.d_u : (c == 1 ? m.d_v : m.d_w);
|
||
double h_min = fmin(fmin(g.dx, g.dy), g.dz);
|
||
if (fabs(dist[t]) > 4.0 * h_min) return;
|
||
const double* ubt = c == 0 ? m.ub_u : (c == 1 ? m.ub_v : m.ub_w);
|
||
double* out = c == 0 ? f.u : (c == 1 ? f.v : f.w);
|
||
out[t] = ubt[t];
|
||
}
|