rustyverse/rustytorch
2
0
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases Wiki Activity

embedded3 R6-2 step 2: the band's surface velocity at the foot (ub) and the rebuild's imposition on the device (RTX_E3_UB_DEVICE=1, default off)

Browse Source 
e3_geom_ub / e3_geom_impose / e3_geom_seam: the flag body's surface velocity
(DeviceSdf::vel, the centreline's velocity per point) at the foot from the
trilinear interpolant of the device corner phi, and at the face centre for the
band's solid faces of the uploaded field (a cut mask has no ghosts); the host
mirror of ub by a packed download. RTX_E3_BAND_CHECK=1 compares the device ub
with Mask::surface_velocity_at and the imposed u, v, w with impose_from, bit for
bit; check_against now compares ub on every band face (zeros included).

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
This commit is contained in:
Omar Sobh
2026-09-25 01:30:29 -05:00
co-authored by Claude Opus 5.5
parent 6ae5312b6f
commit 2c054ccc91
8 changed files with 583 additions and 94 deletions
Download Patch File Download Diff File Expand all files Collapse all files
crates/specialized/rtx-cfd/src/kernels/cuda/e3_geom.cu
+156
View File
@@ -275,3 +275,159 @@ extern "C" __global__ void e3_geom_gather(
long long s = (long long) w * idx[t];
for (int c = 0; c < w; ++c) out[(long long) w * t + c] = src[s + c];
}
/*
* R6-2 step 2: the body's surface velocity on the device (the flag test's
* host closure: the centreline's velocity interpolated at the capsule's
* closest point where the span-cut capsule is not farther than the circle,
* zero on the circle, no z component), and its two uses on a moving cut
* mask:
* e3_geom_ub the imposition band's faces: the velocity at the foot
* (`Mask::surface_velocity_at`: the face centre plus its
* centroid shift when the centroid foot is on, φ and its
* unit gradient from the trilinear interpolant of the
* corner φ, the foot x − s n) into the dense table and a
* packed copy for the host mirror;
* e3_geom_impose the band's solid faces (not open at the instant, off the
* domain sides the host loops skip): the velocity at the
* face centre into the field (`Mask::impose_from`, whose
* ghost lists are empty on a cut mask);
* e3_geom_seam the periodic seam: w at k = nz takes k = 0's where both
* are solid.
*/
__device__ double body_velocity(double x, double y, double z, int c, const GeomSdf& s,
const double* __restrict__ poly, const double* __restrict__ vel)
{
double best = 1.0 / 0.0, vx = 0.0, vy = 0.0;
for (int m = 0; m + 1 < s.npts; ++m) {
double ax = poly[2 * m], ay = poly[2 * m + 1];
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
double ex = bx - ax, ey = by - ay;
double l2 = ex * ex + ey * ey;
double u = ((x - ax) * ex + (y - ay) * ey) / l2;
if (u < 0.0) u = 0.0;
if (u > 1.0) u = 1.0;
double px = ax + u * ex, py = ay + u * ey;
double qx = x - px, qy = y - py;
double d = sqrt(qx * qx + qy * qy);
if (d < best) {
best = d;
double avx = vel[2 * m], avy = vel[2 * m + 1];
double bvx = vel[2 * m + 2], bvy = vel[2 * m + 3];
vx = avx + u * (bvx - avx);
vy = avy + u * (bvy - avy);
}
}
double df = best - s.half;
if (s.flag_cut) df = span_cut(df, z, s);
double ex0 = x - s.cx, ey0 = y - s.cy;
double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc;
if (s.cyl_cut) dc = span_cut(dc, z, s);
if (df <= dc) return c == 0 ? vx : (c == 1 ? vy : 0.0);
return 0.0;
}
__device__ __forceinline__ double ub_node(const GeomGrid& g, const double* __restrict__ phi,
long long i, long long j, long long k)
{
if (i < 0) i = 0;
if (i > g.nx) i = g.nx;
if (j < 0) j = 0;
if (j > g.ny) j = g.ny;
if (k < 0) k = 0;
if (k > g.nz) k = g.nz;
return phi[(k * (g.ny + 1) + j) * (g.nx + 1) + i];
}
__device__ __forceinline__ double ub_lerp(double a, double b, double f) { return a + f * (b - a); }
/* The face centre of face `f` of component `c`. */
__device__ __forceinline__ void face_centre(const GeomGrid& g, int c, long long f, double* x)
{
long long ni = g.nx + (c == 0), nj = g.ny + (c == 1);
long long k = f / (nj * ni), j = (f / ni) % nj, i = f % ni;
x[0] = ((double) i + (c == 0 ? 0.0 : 0.5)) * g.dx;
x[1] = ((double) j + (c == 1 ? 0.0 : 0.5)) * g.dy;
x[2] = ((double) k + (c == 2 ? 0.0 : 0.5)) * g.dz;
}
extern "C" __global__ void e3_geom_ub(
GeomGrid g, GeomSdf s, const double* __restrict__ poly, const double* __restrict__ vel,
const double* __restrict__ phi, int c, const unsigned int* __restrict__ faces, long long nf,
const double* __restrict__ shift, int use_shift,
double* __restrict__ ub, double* __restrict__ packed)
{
long long q = (long long) blockIdx.x * blockDim.x + threadIdx.x;
if (q >= nf) return;
long long f = faces[q];
double xc[3];
face_centre(g, c, f, xc);
double x0 = xc[0], x1 = xc[1], x2 = xc[2];
if (use_shift) {
x0 = x0 + shift[3 * f];
x1 = x1 + shift[3 * f + 1];
x2 = x2 + shift[3 * f + 2];
}
/* the interpolant's value and unit gradient */
double gx = x0 / g.dx, gy = x1 / g.dy, gz = x2 / g.dz;
long long i0 = (long long) floor(gx), j0 = (long long) floor(gy), k0 = (long long) floor(gz);
double fx = gx - (double) i0, fy = gy - (double) j0, fz = gz - (double) k0;
double c000 = ub_node(g, phi, i0, j0, k0), c100 = ub_node(g, phi, i0 + 1, j0, k0);
double c010 = ub_node(g, phi, i0, j0 + 1, k0), c110 = ub_node(g, phi, i0 + 1, j0 + 1, k0);
double c001 = ub_node(g, phi, i0, j0, k0 + 1), c101 = ub_node(g, phi, i0 + 1, j0, k0 + 1);
double c011 = ub_node(g, phi, i0, j0 + 1, k0 + 1), c111 = ub_node(g, phi, i0 + 1, j0 + 1, k0 + 1);
double c00 = ub_lerp(c000, c100, fx);
double c10 = ub_lerp(c010, c110, fx);
double c01 = ub_lerp(c001, c101, fx);
double c11 = ub_lerp(c011, c111, fx);
double c0 = ub_lerp(c00, c10, fy);
double c1 = ub_lerp(c01, c11, fy);
double sd = ub_lerp(c0, c1, fz);
double dx0 = ub_lerp(c100 - c000, c110 - c010, fy);
double dx1 = ub_lerp(c101 - c001, c111 - c011, fy);
double px = ub_lerp(dx0, dx1, fz) / g.dx;
double dy0 = ub_lerp(c010 - c000, c110 - c100, fx);
double dy1 = ub_lerp(c011 - c001, c111 - c101, fx);
double py = ub_lerp(dy0, dy1, fz) / g.dy;
double dz0 = ub_lerp(c001 - c000, c101 - c100, fx);
double dz1 = ub_lerp(c011 - c010, c111 - c110, fx);
double pz = ub_lerp(dz0, dz1, fy) / g.dz;
double norm = sqrt(px * px + py * py + pz * pz);
double n0 = 1.0, n1 = 0.0, n2 = 0.0;
if (norm > 0.0) {
n0 = px / norm;
n1 = py / norm;
n2 = pz / norm;
}
double v = body_velocity(x0 - sd * n0, x1 - sd * n1, x2 - sd * n2, c, s, poly, vel);
ub[f] = v;
packed[q] = v;
}
extern "C" __global__ void e3_geom_impose(
GeomGrid g, GeomSdf s, const double* __restrict__ poly, const double* __restrict__ vel,
int c, const unsigned int* __restrict__ faces, long long nf,
const int* __restrict__ open, double* __restrict__ field)
{
long long q = (long long) blockIdx.x * blockDim.x + threadIdx.x;
if (q >= nf) return;
long long f = faces[q];
if (open[f] != 0) return;
long long ni = g.nx + (c == 0), nj = g.ny + (c == 1);
long long j = (f / ni) % nj, i = f % ni;
/* the host loops: u over i in 1..nx, v over j in 1..ny, every w */
if (c == 0 && (i == 0 || i == g.nx)) return;
if (c == 1 && (j == 0 || j == g.ny)) return;
double xc[3];
face_centre(g, c, f, xc);
field[f] = body_velocity(xc[0], xc[1], xc[2], c, s, poly, vel);
}
extern "C" __global__ void e3_geom_seam(GeomGrid g, const int* __restrict__ open_w, double* __restrict__ w)
{
long long n = (long long) blockIdx.x * blockDim.x + threadIdx.x;
long long nxy = (long long) g.nx * g.ny;
if (n >= nxy) return;
long long fn = (long long) g.nz * nxy + n;
if (open_w[n] == 0 && open_w[fn] == 0) w[fn] = w[n];
}