diff --git a/docs/theory/2d_solver/src/cpu.rs b/docs/theory/2d_solver/src/cpu.rs
index da3b7f0..9017e3a 100644
--- a/docs/theory/2d_solver/src/cpu.rs
+++ b/docs/theory/2d_solver/src/cpu.rs
@@ -8,7 +8,7 @@
use rayon::prelude::*;
-use crate::math::{self, Body, Kbc, Link, LinkKind, WallModel, R, CX, CY, OPP, Q};
+use crate::math::{self, Kbc, Link, LinkKind, Scene, WallModel, MAX_BODY_BUCKETS, R, CX, CY, OPP, Q};
use crate::{Collision, FieldKind, Spec, StepRec};
// ─────────────────────────────────────────────────────────────────────────────
@@ -30,21 +30,30 @@ pub struct Geom {
}
impl Geom {
- /// Собрать геометрию уровня по телу. `solid` определяется знаком SDF: φ ≤ 0 — тело.
- pub fn build(nx: usize, ny: usize, body: &Body) -> Self {
+ /// Собрать геометрию уровня по сцене. `solid` определяется знаком SDF сцены (минимум по
+ /// телам): φ ≤ 0 — тело. Каждому линку проставляется индекс тела, которого он касается,
+ /// чтобы сила считалась по телам отдельно.
+ pub fn build(nx: usize, ny: usize, scene: &Scene) -> Self {
let n = nx * ny;
let mut phi = vec![0.0; n];
let mut solid = vec![false; n];
+ let mut owner = vec![0u8; n];
for y in 0..ny {
for x in 0..nx {
- let p = body.sdf(x as R, y as R);
- phi[y * nx + x] = p;
- solid[y * nx + x] = p <= 0.0;
+ let k = y * nx + x;
+ let p = scene.sdf(x as R, y as R);
+ phi[k] = p;
+ solid[k] = p <= 0.0;
+ if solid[k] {
+ owner[k] =
+ scene.nearest(x as R, y as R).min(MAX_BODY_BUCKETS - 1) as u8;
+ }
}
}
- let links = build_links(nx, ny, &solid, &phi);
+ let links = build_links(nx, ny, &solid, &phi, &owner);
let wall_nodes = build_wall_nodes(&links);
- Geom { nx, ny, solid, links, wall_nodes, body_cx: body.cx, body_cy: body.cy }
+ let (bcx, bcy) = scene.center();
+ Geom { nx, ny, solid, links, wall_nodes, body_cx: bcx, body_cy: bcy }
}
}
@@ -59,7 +68,7 @@ fn nb(nx: usize, ny: usize, x: usize, y: usize, i: usize, sign: i32) -> usize {
/// Для каждого жидкого узла и каждого направления, упирающегося в тело, решаем, какая из трёх
/// формул Bouzidi применима, и считаем долю пересечения q из SDF. Делается один раз.
-fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R]) -> Vec {
+fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R], owner: &[u8]) -> Vec {
let mut links = Vec::new();
for y in 0..ny {
for x in 0..nx {
@@ -88,7 +97,7 @@ fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R]) -> Vec {
ib: OPP[i] as u8,
kind,
q,
- body: true,
+ body: owner[s],
});
}
}
@@ -435,26 +444,26 @@ impl Level {
}
}
- /// Сила и момент на теле по GMEM, суммой по линкам тела.
- fn force(&self) -> (R, R, R) {
- let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0);
+ /// Сила и момент по GMEM, разложенные по телам: `out[b] = (F_x, F_y, T_z)` для тела b.
+ /// Момент считается вокруг центра ПЕРВОГО тела — общей точки отсчёта для всей сцены.
+ fn force(&self) -> [[R; 3]; MAX_BODY_BUCKETS] {
+ let mut out = [[0.0; 3]; MAX_BODY_BUCKETS];
let nx = self.nx;
for l in &self.geom.links {
- if !l.body {
- continue;
- }
let node = l.node as usize;
let i = l.i as usize;
- let (dfx, dfy) = math::gmem_link(i, self.post[node][i], self.f[node][l.ib as usize], 0.0, 0.0);
- fx += dfx;
- fy += dfy;
+ let (dfx, dfy) =
+ math::gmem_link(i, self.post[node][i], self.f[node][l.ib as usize], 0.0, 0.0);
// плечо — до ТОЧКИ ПЕРЕСЕЧЕНИЯ линка со стенкой r_w = r_f + q·c_i, а не до узла:
// узел дал бы ошибку плеча до целой ячейки, а q уже посчитан
let rx = (node % nx) as R - self.geom.body_cx + l.q * CX[i] as R;
let ry = (node / nx) as R - self.geom.body_cy + l.q * CY[i] as R;
- tz += rx * dfy - ry * dfx;
+ let b = (l.body as usize).min(MAX_BODY_BUCKETS - 1);
+ out[b][0] += dfx;
+ out[b][1] += dfy;
+ out[b][2] += rx * dfy - ry * dfx;
}
- (fx, fy, tz)
+ out
}
#[inline]
@@ -463,6 +472,20 @@ impl Level {
}
}
+/// Свести стартовую скорость к нулю на подходе к телу, чтобы на первом шаге не возникло
+/// разрыва. Вдали от тела поле не трогается вовсе.
+fn taper_start(lvl: &mut Level, scene: &Scene, u0: (R, R), width: R) {
+ let nx = lvl.nx;
+ for (n, cell) in lvl.f.iter_mut().enumerate() {
+ let (x, y) = ((n % nx) as R, (n / nx) as R);
+ let k = math::wall_taper(scene.sdf(x, y), width);
+ if k < 1.0 {
+ *cell = math::feq(1.0, u0.0 * k, u0.1 * k);
+ }
+ }
+ lvl.post.copy_from_slice(&lvl.f);
+}
+
// ─────────────────────────────────────────────────────────────────────────────
// Статистика KBC за шаг
// ─────────────────────────────────────────────────────────────────────────────
@@ -730,22 +753,15 @@ pub struct Sim {
impl Sim {
pub fn new(spec: Spec) -> Sim {
let (nx, ny) = (spec.nx, spec.ny);
- let geom0 = Geom::build(nx, ny, &spec.body);
+ let geom0 = Geom::build(nx, ny, &spec.scene);
- // губка: плавный рост вязкости в последних sponge_len столбцах. Канал
- // «скорость-вход + давление-выход» — недодемпфированный акустический резонатор,
- // губка гасит и вихри, и акустику до прихода на выход.
- let beta0: Vec = (0..nx)
- .map(|x| {
- if spec.sponge_len == 0 {
- return spec.beta0;
- }
- let start = nx - 1 - spec.sponge_len;
- let s = math::smoothstep((x as R - start as R) / spec.sponge_len as R);
- let nu = math::nu_of_beta(spec.beta0);
- math::beta_of_nu(nu * (1.0 + (spec.sponge_mult - 1.0) * s))
- })
- .collect();
+ let beta0 = math::beta_profile(
+ nx,
+ spec.beta0,
+ spec.sponge_in,
+ spec.sponge_len,
+ spec.sponge_mult,
+ );
// стартовое поле: либо сразу набегающий поток, либо покой
let u0 = if spec.init_uniform {
@@ -754,23 +770,30 @@ impl Sim {
} else {
(0.0, 0.0)
};
- let l0 = Level::new(nx, ny, beta0, geom0, u0);
+ let mut l0 = Level::new(nx, ny, beta0, geom0, u0);
+ if spec.init_uniform {
+ taper_start(&mut l0, &spec.scene, u0, spec.init_taper);
+ }
let fluid_count = l0.geom.solid.iter().filter(|s| !**s).count() as R;
let (l1, patch) = if spec.refine > 1 {
let (ax, bx, ay, by) = spec.patch.expect("refine > 1 требует патч");
let r = spec.refine;
- let body1 = spec.body.refined(r as R, ax as R, ay as R);
+ let scene1 = spec.scene.refined(r as R, ax as R, ay as R);
let nfx = r * (bx - ax) + 1;
let nfy = r * (by - ay) + 1;
- let geom1 = Geom::build(nfx, nfy, &body1);
+ let geom1 = Geom::build(nfx, nfy, &scene1);
// τ_f = r(τ_c − ½) + ½ ⇒ одинаковая ν на обоих уровнях
let tau0 = 1.0 / (2.0 * spec.beta0);
let tau1 = r as R * (tau0 - 0.5) + 0.5;
let beta1 = 1.0 / (2.0 * tau1);
let r01 = tau1 / (r as R * tau0);
let patch = Patch::new(&spec, &l0.geom, &geom1.solid, r01);
- (Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1, u0)), Some(patch))
+ let mut lvl1 = Level::new(nfx, nfy, vec![beta1; nfx], geom1, u0);
+ if spec.init_uniform {
+ taper_start(&mut lvl1, &scene1, u0, spec.init_taper * r as R);
+ }
+ (Some(lvl1), Some(patch))
} else {
(None, None)
};
@@ -842,7 +865,7 @@ impl Sim {
self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap);
// ── уровень 1: r подшагов с временной интерполяцией рамки ──
- let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0);
+ let mut fb = [[0.0 as R; 3]; MAX_BODY_BUCKETS];
if let (Some(l1), Some(p)) = (self.l1.as_mut(), self.patch.as_ref()) {
p.ghost_values(&self.pre, &mut self.gh_old);
p.ghost_values(&self.l0.f, &mut self.gh_new);
@@ -852,23 +875,24 @@ impl Sim {
l1.apply_wall(wall);
// силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на
// последнем подшаге даёт лишний шум в рядах при том же среднем
- let (a, b, c) = l1.force();
- fx += a;
- fy += b;
- tz += c;
+ let g = l1.force();
+ for k in 0..MAX_BODY_BUCKETS {
+ for c in 0..3 {
+ fb[k][c] += g[k][c];
+ }
+ }
let w = (s + 1) as R / p.r as R;
p.fill(&mut l1.f, &self.gh_old, &self.gh_new, w);
}
let inv = 1.0 / p.r as R;
- fx *= inv;
- fy *= inv;
- tz *= inv;
+ for k in 0..MAX_BODY_BUCKETS {
+ for c in 0..3 {
+ fb[k][c] *= inv;
+ }
+ }
p.restrict_to(&l1.f, &mut self.l0.f);
} else {
- let (a, b, c) = self.l0.force();
- fx = a;
- fy = b;
- tz = c;
+ fb = self.l0.force();
}
// ── диагностика ──
@@ -892,11 +916,18 @@ impl Sim {
.reduce(|| (0.0, 0.0), |a, b| (a.0 + b.0, a.1.max(b.1)));
self.step_index += 1;
+ let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0);
+ for k in 0..MAX_BODY_BUCKETS {
+ fx += fb[k][0];
+ fy += fb[k][1];
+ tz += fb[k][2];
+ }
StepRec {
step: t,
fx,
fy,
tz,
+ body: fb,
uy_probe: probe.2,
rho_mean: rho_sum / self.fluid_count,
max_u,
@@ -911,8 +942,8 @@ impl Sim {
/// Характерный размер тела в единицах того уровня, где снимается сила.
pub fn force_ref_size(&self) -> R {
match &self.patch {
- Some(p) => self.spec.body.d * p.r as R,
- None => self.spec.body.d,
+ Some(p) => self.spec.scene.ref_size() * p.r as R,
+ None => self.spec.scene.ref_size(),
}
}
@@ -958,6 +989,22 @@ impl Sim {
(out, &self.l0.geom.solid)
}
+ /// Срез вдоль осевой линии: (ρ, u_x) по каждому столбцу. Из последовательности таких
+ /// срезов складывается x–t диаграмма, по которой видно, бежит возмущение со скоростью
+ /// звука или конвекции и есть ли стоячие узлы.
+ pub fn sample_centerline(&self) -> (Vec, Vec) {
+ let y = self.l0.ny / 2;
+ let nx = self.l0.nx;
+ let mut rho = Vec::with_capacity(nx);
+ let mut ux = Vec::with_capacity(nx);
+ for x in 0..nx {
+ let (r, a, _) = self.l0.macros_at(y * nx + x);
+ rho.push(r);
+ ux.push(a);
+ }
+ (rho, ux)
+ }
+
/// Есть ли в поле NaN/inf — признак развала счёта.
pub fn is_finite(&self) -> bool {
self.l0.f.par_iter().all(|c| c.iter().all(|v| v.is_finite()))
@@ -1052,8 +1099,8 @@ mod tests {
fn wall_is_subgrid_not_staircase() {
let mut means = Vec::new();
for off in [0.0, 0.25, 0.5] {
- let body = Body::new(math::ShapeKind::Cylinder, 60.0, 40.0 + off, 24.0, 0.0, 0.3);
- let g = Geom::build(120, 80, &body);
+ let body = math::Body::new(math::ShapeKind::Cylinder, 60.0, 40.0 + off, 24.0, 0.0, 0.3);
+ let g = Geom::build(120, 80, &Scene::single(body));
let w = g.wall_stats();
assert!(w.links > 100, "смещение {off}: линков всего {}", w.links);
assert_eq!(w.simple, 0, "смещение {off}: {} ступенчатых линков", w.simple);
@@ -1069,8 +1116,8 @@ mod tests {
/// без пересечений, в том же порядке.
#[test]
fn wall_node_index_partitions_links() {
- let body = Body::new(math::ShapeKind::Naca, 60.0, 40.0, 30.0, 12.0, 0.25);
- let g = Geom::build(120, 80, &body);
+ let body = math::Body::new(math::ShapeKind::Naca, 60.0, 40.0, 30.0, 12.0, 0.25);
+ let g = Geom::build(120, 80, &Scene::single(body));
let mut covered = 0usize;
for w in &g.wall_nodes {
assert_eq!(w.first as usize, covered, "разрыв в индексе граничных узлов");
diff --git a/docs/theory/2d_solver/src/gpu.rs b/docs/theory/2d_solver/src/gpu.rs
index 42824d9..50b96ea 100644
--- a/docs/theory/2d_solver/src/gpu.rs
+++ b/docs/theory/2d_solver/src/gpu.rs
@@ -102,7 +102,8 @@ struct GLink {
ib: u32,
kind: u32,
q: f32,
- _p: [u32; 2],
+ body: u32,
+ _p: u32,
}
#[repr(C)]
@@ -134,6 +135,8 @@ struct Results {
cnt: f32,
degen: f32,
xineg: f32,
+ /// Сила и момент по телам: [b*3 + c].
+ body: [f32; 12],
}
// ─────────────────────────────────────────────────────────────────────────────
@@ -146,7 +149,7 @@ struct LevelParams { n:u32, nx:u32, ny:u32, nlinks:u32, bcx:f32, bcy:f32, probe_
struct Dyn { ux_in:f32, uy_in:f32, rho_out:f32, kbc_model:u32, outlet_extrap:u32, nparts:u32, refine:u32, collision:u32, slot:u32, wall_hrr:u32, dp2:u32, dp3:u32 };
struct Substep { idx:u32, p0:u32, p1:u32, p2:u32 };
struct AmrParams { nghost:u32, nrestrict:u32, ccount:u32, fcount:u32, r01:f32, rfc:f32, w:f32, pad:f32 };
-struct GLink { node:u32, far:u32, i:u32, ib:u32, kind:u32, q:f32, p0:u32, p1:u32 };
+struct GLink { node:u32, far:u32, i:u32, ib:u32, kind:u32, q:f32, body:u32, p1:u32 };
struct GGhost { fine:u32, c00:u32, c10:u32, c01:u32, c11:u32, tx:f32, ty:f32, p0:u32 };
struct Partial { rho:f32, maxu:f32, gsum:f32, gmin:f32, gmax:f32, cnt:f32, degen:f32, xineg:f32 };
@@ -184,7 +187,9 @@ const CS2 : f32 = 0.3333333333;
// рабочие слоты силы по подшагам. История нужна, чтобы не синхронизироваться с устройством
// на каждом шаге: результаты копятся и читаются пачкой (см. HIST в gpu.rs).
const HIST: u32 = 128u;
-const FORCE_BASE: u32 = 1536u; // = HIST*12u
+const SLOT: u32 = 24u; // чисел на шаг: 12 общих + 4 тела по 3
+const MAXB: u32 = 4u; // вёдер силы по телам
+const FORCE_BASE: u32 = 3072u; // = HIST*SLOT
// Компенсированное сложение Кэхена–Ноймайера: возвращает (сумма, накопленная поправка).
// Наивная сумма по 10^5…10^7 значений в f32 съедает ~log2(N) бит; здесь потеря не копится,
@@ -404,9 +409,10 @@ fn k_force(@builtin(local_invocation_id) lid: vec3) {
let t = lid.x;
var cx = array(0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0);
var cy = array(0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0);
- var sx = vec2(0.0, 0.0);
- var sy = vec2(0.0, 0.0);
- var sz = vec2(0.0, 0.0);
+ // накопители на каждое тело держатся в регистрах, редукция потом идёт по телам подряд
+ var pfx = array(0.0, 0.0, 0.0, 0.0);
+ var pfy = array(0.0, 0.0, 0.0, 0.0);
+ var ptz = array(0.0, 0.0, 0.0, 0.0);
var k = t;
loop {
if (k >= P.nlinks) { break; }
@@ -415,31 +421,36 @@ fn k_force(@builtin(local_invocation_id) lid: vec3) {
let fb = f[L.ib*P.n + L.node];
let dfx = cx[L.i]*fp + cx[L.i]*fb;
let dfy = cy[L.i]*fp + cy[L.i]*fb;
- sx = kadd(sx.x, sx.y, dfx);
- sy = kadd(sy.x, sy.y, dfy);
// плечо до точки пересечения линка со стенкой, а не до узла
let rx = f32(L.node % P.nx) - P.bcx + L.q*cx[L.i];
let ry = f32(L.node / P.nx) - P.bcy + L.q*cy[L.i];
- sz = kadd(sz.x, sz.y, rx*dfy - ry*dfx);
+ let b = min(L.body, MAXB - 1u);
+ pfx[b] = pfx[b] + dfx;
+ pfy[b] = pfy[b] + dfy;
+ ptz[b] = ptz[b] + rx*dfy - ry*dfx;
k = k + 256u;
}
- wfx[t] = sx.x + sx.y; wfy[t] = sy.x + sy.y; wtz[t] = sz.x + sz.y;
- workgroupBarrier();
- var s = 128u;
- loop {
- if (s == 0u) { break; }
- if (t < s) {
- wfx[t] = wfx[t] + wfx[t + s];
- wfy[t] = wfy[t] + wfy[t + s];
- wtz[t] = wtz[t] + wtz[t + s];
+ for (var b = 0u; b < MAXB; b = b + 1u) {
+ wfx[t] = pfx[b]; wfy[t] = pfy[b]; wtz[t] = ptz[b];
+ workgroupBarrier();
+ var s = 128u;
+ loop {
+ if (s == 0u) { break; }
+ if (t < s) {
+ wfx[t] = wfx[t] + wfx[t + s];
+ wfy[t] = wfy[t] + wfy[t + s];
+ wtz[t] = wtz[t] + wtz[t + s];
+ }
+ workgroupBarrier();
+ s = s >> 1u;
+ }
+ if (t == 0u) {
+ let o = FORCE_BASE + (S.idx*MAXB + b)*4u;
+ results[o + 0u] = wfx[0];
+ results[o + 1u] = wfy[0];
+ results[o + 2u] = wtz[0];
}
workgroupBarrier();
- s = s >> 1u;
- }
- if (t == 0u) {
- results[FORCE_BASE + S.idx*4u + 0u] = wfx[0];
- results[FORCE_BASE + S.idx*4u + 1u] = wfy[0];
- results[FORCE_BASE + S.idx*4u + 2u] = wtz[0];
}
}
@@ -477,7 +488,7 @@ fn k_stats1(@builtin(global_invocation_id) gid: vec3,
var fv: array;
load9(&fv, nd, P.n);
let m = macros9(fv);
- results[D.slot*12u + 3u] = m.z;
+ results[D.slot*SLOT + 3u] = m.z;
}
wp[t] = p;
workgroupBarrier();
@@ -554,17 +565,26 @@ fn k_stats2(@builtin(local_invocation_id) lid: vec3) {
qgmin = min(qgmin, o.gmin);
qgmax = max(qgmax, o.gmax);
}
- var fx = 0.0; var fy = 0.0; var tz = 0.0;
- for (var s = 0u; s < D.refine; s = s + 1u) {
- fx = fx + results[FORCE_BASE + s*4u + 0u];
- fy = fy + results[FORCE_BASE + s*4u + 1u];
- tz = tz + results[FORCE_BASE + s*4u + 2u];
- }
+ // сила усредняется по подшагам L1 и раскладывается по телам; общая — их сумма
let inv = 1.0 / f32(D.refine);
- let o = D.slot*12u;
- results[o + 0u] = fx*inv;
- results[o + 1u] = fy*inv;
- results[o + 2u] = tz*inv;
+ let o = D.slot*SLOT;
+ var tfx = 0.0; var tfy = 0.0; var ttz = 0.0;
+ for (var b = 0u; b < MAXB; b = b + 1u) {
+ var fx = 0.0; var fy = 0.0; var tz = 0.0;
+ for (var sb = 0u; sb < D.refine; sb = sb + 1u) {
+ let g = FORCE_BASE + (sb*MAXB + b)*4u;
+ fx = fx + results[g + 0u];
+ fy = fy + results[g + 1u];
+ tz = tz + results[g + 2u];
+ }
+ results[o + 12u + b*3u + 0u] = fx*inv;
+ results[o + 12u + b*3u + 1u] = fy*inv;
+ results[o + 12u + b*3u + 2u] = tz*inv;
+ tfx = tfx + fx*inv; tfy = tfy + fy*inv; ttz = ttz + tz*inv;
+ }
+ results[o + 0u] = tfx;
+ results[o + 1u] = tfy;
+ results[o + 2u] = ttz;
results[o + 4u] = qrho.x + qrho.y;
results[o + 5u] = qmaxu;
results[o + 6u] = qgsum.x + qgsum.y;
@@ -784,11 +804,30 @@ impl GpuLevel {
}
}
-fn soa_equilibrium(n: usize, u0: (R, R)) -> Vec {
- let fe = math::feq(1.0, u0.0, u0.1);
+/// Стартовое поле в раскладке SoA. Если задана сцена, скорость сводится к нулю на подходе
+/// к телу (см. `math::wall_taper`) — иначе на первом шаге возникает разрыв и импульс сжатия.
+fn soa_equilibrium(
+ n: usize,
+ nx: usize,
+ u0: (R, R),
+ taper: Option<(&math::Scene, R)>,
+) -> Vec {
let mut v = vec![0.0f32; 9 * n];
- for i in 0..9 {
- for k in 0..n {
+ let plain = math::feq(1.0, u0.0, u0.1);
+ for k in 0..n {
+ let fe = match taper {
+ Some((scene, w)) => {
+ let (x, y) = ((k % nx) as R, (k / nx) as R);
+ let c = math::wall_taper(scene.sdf(x, y), w);
+ if c < 1.0 {
+ math::feq(1.0, u0.0 * c, u0.1 * c)
+ } else {
+ plain
+ }
+ }
+ None => plain,
+ };
+ for i in 0..9 {
v[i * n + k] = fe[i] as f32;
}
}
@@ -807,9 +846,10 @@ fn make_level(
flags: u32,
u0: (R, R),
wall_layout: &wgpu::BindGroupLayout,
+ taper: Option<(&math::Scene, R)>,
) -> GpuLevel {
let n = nx * ny;
- let init = soa_equilibrium(n, u0);
+ let init = soa_equilibrium(n, nx, u0, taper);
let mkf = |label: &str| {
device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
label: Some(label),
@@ -841,7 +881,8 @@ fn make_level(
math::LinkKind::Simple => 2,
},
q: l.q as f32,
- _p: [0; 2],
+ body: l.body as u32,
+ _p: 0,
})
.collect();
let links_buf = storage_init(device, "links", bytemuck::cast_slice(&links));
@@ -997,20 +1038,14 @@ impl Sim {
}
// ── топология берётся из процессорного бэкенда, а не строится заново ──
- let geom0 = cpu::Geom::build(spec.nx, spec.ny, &spec.body);
+ let geom0 = cpu::Geom::build(spec.nx, spec.ny, &spec.scene);
let fluid_count = geom0.solid.iter().filter(|s| !**s).count() as R;
- let beta0: Vec = (0..spec.nx)
- .map(|x| {
- if spec.sponge_len == 0 {
- return spec.beta0 as f32;
- }
- let start = spec.nx - 1 - spec.sponge_len;
- let s = math::smoothstep((x as R - start as R) / spec.sponge_len as R);
- let nu = math::nu_of_beta(spec.beta0);
- math::beta_of_nu(nu * (1.0 + (spec.sponge_mult - 1.0) * s)) as f32
- })
- .collect();
+ let beta0: Vec =
+ math::beta_profile(spec.nx, spec.beta0, spec.sponge_in, spec.sponge_len, spec.sponge_mult)
+ .into_iter()
+ .map(|v| v as f32)
+ .collect();
let module = device.create_shader_module(wgpu::ShaderModuleDescriptor {
label: Some("kbc2d.wgsl"),
@@ -1122,6 +1157,7 @@ impl Sim {
flags0,
u0,
&bgl_wall,
+ if spec.init_uniform { Some((&spec.scene, spec.init_taper)) } else { None },
);
let pre = device.create_buffer(&wgpu::BufferDescriptor {
@@ -1137,10 +1173,10 @@ impl Sim {
if spec.refine > 1 {
let (ax, bx, ay, by) = spec.patch.ok_or("refine > 1 требует патч")?;
let r = spec.refine;
- let body1 = spec.body.refined(r as R, ax as R, ay as R);
+ let scene1 = spec.scene.refined(r as R, ax as R, ay as R);
let nfx = r * (bx - ax) + 1;
let nfy = r * (by - ay) + 1;
- let geom1 = cpu::Geom::build(nfx, nfy, &body1);
+ let geom1 = cpu::Geom::build(nfx, nfy, &scene1);
let tau0 = 1.0 / (2.0 * spec.beta0);
let tau1 = r as R * (tau0 - 0.5) + 0.5;
let beta1 = vec![(1.0 / (2.0 * tau1)) as f32; nfx];
@@ -1151,7 +1187,12 @@ impl Sim {
let flags1 = if probe_on_fine { 2 } else { 0 };
let lvl1 =
make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1,
- u0, &bgl_wall);
+ u0, &bgl_wall,
+ if spec.init_uniform {
+ Some((&scene1, spec.init_taper * r as R))
+ } else {
+ None
+ });
let ghosts: Vec = patch
.ghosts()
@@ -1262,8 +1303,8 @@ impl Sim {
});
let d_ref = match spec.patch {
- Some(_) if spec.refine > 1 => spec.body.d * spec.refine as R,
- _ => spec.body.d,
+ Some(_) if spec.refine > 1 => spec.scene.ref_size() * spec.refine as R,
+ _ => spec.scene.ref_size(),
};
Ok(Sim {
@@ -1516,11 +1557,18 @@ impl Sim {
fn make_rec(&self, t: u64, res: &Results) -> StepRec {
let cnt = if res.cnt > 0.0 { res.cnt as R } else { self.fluid_count };
+ let mut body = [[0.0 as R; 3]; crate::MAX_BODY_BUCKETS];
+ for b in 0..crate::MAX_BODY_BUCKETS.min(4) {
+ for c in 0..3 {
+ body[b][c] = res.body[b * 3 + c] as R;
+ }
+ }
StepRec {
step: t,
fx: res.fx as R,
fy: res.fy as R,
tz: res.tz as R,
+ body,
uy_probe: res.uy_probe as R,
rho_mean: res.rho_sum as R / self.fluid_count,
max_u: res.max_u as R,
@@ -1564,6 +1612,24 @@ impl Sim {
out
}
+ /// Срез вдоль осевой линии: (ρ, u_x) по столбцам. На GPU это полное скачивание поля,
+ /// поэтому x–t диагностика включается редким шагом и только там, где нужна.
+ pub fn sample_centerline(&self) -> (Vec, Vec) {
+ let (nx, ny, n) = (self.l0.nx, self.l0.ny, self.l0.n);
+ let raw = self.download_l0();
+ let y = ny / 2;
+ let mut rho = Vec::with_capacity(nx);
+ let mut ux = Vec::with_capacity(nx);
+ for x in 0..nx {
+ let k = y * nx + x;
+ let c: [R; 9] = std::array::from_fn(|i| raw[i * n + k] as R);
+ let (r, a, _) = math::macros(&c);
+ rho.push(r);
+ ux.push(a);
+ }
+ (rho, ux)
+ }
+
pub fn is_finite(&self) -> bool {
self.download_l0().iter().all(|v| v.is_finite())
}
@@ -1616,7 +1682,7 @@ impl Sim {
}
}
}
- let solid = cpu::Geom::build(nx, ny, &self.spec.body).solid;
+ let solid = cpu::Geom::build(nx, ny, &self.spec.scene).solid;
(out, solid)
}
}
diff --git a/docs/theory/2d_solver/src/main.rs b/docs/theory/2d_solver/src/main.rs
index d083d49..770ef23 100644
--- a/docs/theory/2d_solver/src/main.rs
+++ b/docs/theory/2d_solver/src/main.rs
@@ -19,7 +19,7 @@ use std::time::Instant;
use clap::Parser;
-use math::{Body, ShapeKind, Units, R};
+use math::{Body, Scene, ShapeKind, Units, MAX_BODY_BUCKETS, R};
// ─────────────────────────────────────────────────────────────────────────────
// Контракт оркестратор ↔ бэкенд
@@ -57,9 +57,13 @@ impl FieldKind {
#[derive(Clone, Copy, Debug, Default)]
pub struct StepRec {
pub step: u64,
+ /// Суммарная сила и момент по всей сцене.
pub fx: R,
pub fy: R,
pub tz: R,
+ /// Она же по телам: `body[b] = (F_x, F_y, T_z)`. Тела за пределом набора вёдер
+ /// сваливаются в последнее, поэтому сумма по вёдрам всегда точна.
+ pub body: [[R; 3]; MAX_BODY_BUCKETS],
pub uy_probe: R,
pub rho_mean: R,
pub max_u: R,
@@ -75,7 +79,7 @@ pub struct StepRec {
pub struct Spec {
pub nx: usize,
pub ny: usize,
- pub body: Body,
+ pub scene: Scene,
pub refine: usize,
/// Границы патча измельчения в координатах L0: (ax, bx, ay, by).
pub patch: Option<(usize, usize, usize, usize)>,
@@ -90,7 +94,11 @@ pub struct Spec {
pub outlet_extrapolate: bool,
/// Начальное поле уже несёт набегающий поток (иначе — покой с разгоном входа).
pub init_uniform: bool,
+ /// Ширина сведения стартовой скорости к нулю у тела, клеток.
+ pub init_taper: R,
pub sponge_len: usize,
+ /// Длина губки после входа, столбцов.
+ pub sponge_in: usize,
pub sponge_mult: R,
pub collision: Collision,
/// Что входит в сдвиговую часть s (см. `math::KbcModel`).
@@ -167,6 +175,20 @@ struct Cli {
/// Положение центра тела по y, ячеек (по умолчанию ny/2)
#[arg(long, help_heading = "Тело")]
body_y: Option,
+ /// Четырёхзначное обозначение профиля NACA, например 4412 (кривизна 4%, её максимум на
+ /// 40% хорды, толщина 12%). Действует при --shape naca и перекрывает --thickness.
+ #[arg(long, help_heading = "Тело")]
+ naca: Option,
+ /// Вершины произвольного многоугольника при --shape polygon: "x1,y1;x2,y2;..." в клетках,
+ /// отсчёт от центра тела. Обход любой, SDF точный.
+ #[arg(long, allow_hyphen_values = true, help_heading = "Тело")]
+ poly: Option,
+ /// НЕСКОЛЬКО ТЕЛ вместо одного. Список через точку с запятой, каждое — форма и параметры:
+ /// "cylinder:d=24,x=120,y=120; cylinder:d=24,x=200,y=120". Ключи: d — размер, x/y — центр,
+ /// a — угол атаки (град), t — относительная толщина, naca — код профиля,
+ /// p — вершины многоугольника через |. Перекрывает --shape и прочие одиночные ключи.
+ #[arg(long, allow_hyphen_values = true, help_heading = "Тело")]
+ bodies: Option,
// ── время ──
/// Число шагов симуляции. Взаимоисключающе с --time.
@@ -192,6 +214,17 @@ struct Cli {
#[arg(long, default_value = "uniform", value_parser = ["uniform", "rest"],
help_heading = "Время")]
init: String,
+ /// Ширина полосы у тела, на которой стартовая скорость сводится к нулю, клеток.
+ ///
+ /// ЗАМЕРЕНО, что это НЕ ПОМОГАЕТ, и ключ оставлен только для повторной проверки.
+ /// Однородный старт заливает потоком место, где стоит тело, и на первом шаге у стенки
+ /// рождается возмущение плотности ~2.5·10⁻². Сглаживание монотонно давит именно его
+ /// (при ширине 24 клетки — до 3.1·10⁻³, восьмикратно), но ПИК ЗА ПРОГОН при этом даже
+ /// подрастает: 2.54·10⁻² → 3.22·10⁻². Возмущение просто переносится во времени — поток
+ /// всё равно обязан разогнаться вокруг тела, и энергия этого переходного процесса задана
+ /// физикой, а не гладкостью начального поля. Продольную моду ест губка, а не это.
+ #[arg(long, default_value_t = 0.0, help_heading = "Время")]
+ init_taper: R,
// ── численная схема ──
/// Оператор столкновения
@@ -221,6 +254,10 @@ struct Cli {
/// 0 — выключить. Губка гасит продольную моду, которую вязкость сама не гасит.
#[arg(long, help_heading = "Схема")]
sponge_len: Option,
+ /// Длина губки ПОСЛЕ входа, столбцов (0 — выключена). Вход по скорости отражает продольные
+ /// волны как жёсткий поршень; губка у входа гасит их до отражения.
+ #[arg(long, default_value_t = 0, help_heading = "Схема")]
+ sponge_in: usize,
/// Во сколько раз губка поднимает вязкость
#[arg(long, default_value_t = 30.0, help_heading = "Схема")]
sponge_mult: R,
@@ -273,9 +310,22 @@ struct Cli {
/// Число окон в отчёте о сходимости
#[arg(long, default_value_t = 10, help_heading = "Вывод")]
windows: usize,
+ /// Прореживание временных рядов: хранить каждую N-ю запись. Статистика и спектр от этого
+ /// не страдают, пока N много меньше периода схода вихрей, зато сверхдлинные прогоны
+ /// перестают съедать память (на 5·10⁶ шагов полный ряд — сотни мегабайт).
+ #[arg(long, default_value_t = 1, help_heading = "Вывод")]
+ series_every: u64,
/// Выгрузить временные ряды в CSV
#[arg(long, help_heading = "Вывод")]
csv: Option,
+ /// Файл x–t диаграммы: срез ⟨ρ⟩ и u_x вдоль осевой линии, по строке на срез. По наклону
+ /// полос видно, бежит возмущение со скоростью звука или конвекции, а стоячие узлы
+ /// проявляются сразу. Основной инструмент разбора стартовой акустики.
+ #[arg(long, help_heading = "Вывод")]
+ xt: Option,
+ /// Период записи срезов x–t, шагов
+ #[arg(long, default_value_t = 200, help_heading = "Вывод")]
+ xt_every: u64,
}
// ─────────────────────────────────────────────────────────────────────────────
@@ -292,6 +342,65 @@ fn parse_pair(s: &str, what: &str) -> Result<(R, R), String> {
Ok((a, b))
}
+/// Разбор одного тела из мини-синтаксиса `форма:ключ=значение,...`.
+fn parse_body(spec: &str, def_x: R, def_y: R) -> Result {
+ let (name, rest) = match spec.split_once(':') {
+ Some((a, b)) => (a.trim(), b),
+ None => (spec.trim(), ""),
+ };
+ let kind = ShapeKind::from_str(name).ok_or(format!("неизвестная форма «{name}»"))?;
+ let (mut d, mut x, mut y, mut a, mut t) = (24.0, def_x, def_y, 0.0, 0.3);
+ let (mut camber, mut cpos) = (0.0, 0.0);
+ let mut verts: Vec<[R; 2]> = Vec::new();
+ for kv in rest.split(',').filter(|z| !z.trim().is_empty()) {
+ let (k, v) = kv
+ .split_once('=')
+ .ok_or(format!("в «{}» ожидалось ключ=значение", kv.trim()))?;
+ let (k, v) = (k.trim(), v.trim());
+ let num = |v: &str| v.parse::().map_err(|e| format!("{k}: {e}"));
+ match k {
+ "d" | "size" | "c" => d = num(v)?,
+ "x" => x = num(v)?,
+ "y" => y = num(v)?,
+ "a" | "angle" => a = num(v)?,
+ "t" | "thickness" => t = num(v)?,
+ "naca" => {
+ let (m, pp, tt) =
+ Body::naca_code(v).ok_or(format!("naca: ожидалось четыре цифры, «{v}»"))?;
+ camber = m;
+ cpos = pp;
+ t = tt;
+ }
+ "p" | "poly" => {
+ let n: Vec<&str> = v.split('|').collect();
+ if n.len() < 6 || n.len() % 2 != 0 {
+ return Err("p: нужно чётное число координат, минимум три вершины".into());
+ }
+ for pair in n.chunks(2) {
+ verts.push([num(pair[0])?, num(pair[1])?]);
+ }
+ }
+ _ => return Err(format!("неизвестный ключ «{k}»")),
+ }
+ }
+ if kind == ShapeKind::Polygon && verts.len() < 3 {
+ return Err("для polygon нужен ключ p с вершинами".into());
+ }
+ Ok(Body::full(kind, x, y, d, a, t, camber, cpos, &verts))
+}
+
+fn parse_poly(v: &str) -> Result, String> {
+ let mut out = Vec::new();
+ for pt in v.split(';').filter(|z| !z.trim().is_empty()) {
+ let (a, b) = parse_pair(pt, "--poly")?;
+ out.push([a, b]);
+ }
+ if out.len() < 3 {
+ return Err("--poly: нужно минимум три вершины".into());
+ }
+ Ok(out)
+}
+
fn build_spec(cli: &Cli) -> Result {
if cli.nx < 16 || cli.ny < 16 {
return Err("сетка меньше 16×16 не имеет смысла".into());
@@ -306,12 +415,51 @@ fn build_spec(cli: &Cli) -> Result {
return Err("--refine обязан быть ≥ 1".into());
}
- let kind = ShapeKind::from_str(&cli.shape).ok_or("неизвестная форма")?;
let cx = cli.body_x.unwrap_or(cli.nx as R / 4.0);
let cy = cli.body_y.unwrap_or(cli.ny as R / 2.0);
- let body = Body::new(kind, cx, cy, cli.size, cli.body_angle, cli.thickness);
+ let scene = match &cli.bodies {
+ Some(list) => {
+ let mut v = Vec::new();
+ for one in list.split(';').filter(|z| !z.trim().is_empty()) {
+ v.push(parse_body(one, cx, cy)?);
+ }
+ if v.is_empty() {
+ return Err("--bodies: список пуст".into());
+ }
+ Scene { bodies: v }
+ }
+ None => {
+ let kind = ShapeKind::from_str(&cli.shape).ok_or("неизвестная форма")?;
+ let (mut camber, mut cpos, mut thick) = (0.0, 0.0, cli.thickness);
+ if let Some(code) = &cli.naca {
+ let (m, pp, tt) = Body::naca_code(code)
+ .ok_or("--naca: ожидалось четыре цифры, например 4412")?;
+ camber = m;
+ cpos = pp;
+ thick = tt;
+ }
+ let verts = match &cli.poly {
+ Some(v) => parse_poly(v)?,
+ None => Vec::new(),
+ };
+ if kind == ShapeKind::Polygon && verts.is_empty() {
+ return Err("--shape polygon требует --poly с вершинами".into());
+ }
+ Scene::single(Body::full(
+ kind,
+ cx,
+ cy,
+ cli.size,
+ cli.body_angle,
+ thick,
+ camber,
+ cpos,
+ &verts,
+ ))
+ }
+ };
- let units = Units::new(cli.dx, cli.u_phys, cli.u_lat, cli.re, cli.size);
+ let units = Units::new(cli.dx, cli.u_phys, cli.u_lat, cli.re, scene.ref_size());
let beta0 = math::beta_of_nu(units.nu_lat);
let tau0 = 1.0 / (2.0 * beta0);
if tau0 <= 0.5 {
@@ -364,6 +512,17 @@ fn build_spec(cli: &Cli) -> Result {
// как открытый конец. Собственная мода затухает как ν(π/2Nx)², то есть практически не
// затухает, и любой стартовый удар остаётся в домене навсегда. Губка — единственное, что
// её реально ест, поэтому по умолчанию она включена и подобрана под длину домена.
+ // губка входа не должна доставать до патча измельчения
+ if cli.sponge_in > 0 {
+ if let Some((ax, _, _, _)) = patch {
+ if cli.sponge_in + 2 >= ax {
+ return Err(format!(
+ "губка входа ({} столбцов) достаёт до патча (ax={ax})",
+ cli.sponge_in
+ ));
+ }
+ }
+ }
let sponge_len = match cli.sponge_len {
Some(v) => {
if v > 0 {
@@ -421,7 +580,7 @@ fn build_spec(cli: &Cli) -> Result {
Ok(Spec {
nx: cli.nx,
ny: cli.ny,
- body,
+ scene,
refine: cli.refine,
patch,
units,
@@ -433,7 +592,9 @@ fn build_spec(cli: &Cli) -> Result {
pert_dur: cli.pert_dur,
outlet_extrapolate: cli.outlet == "extrapolate",
init_uniform: cli.init == "uniform",
+ init_taper: cli.init_taper.max(0.0),
sponge_len,
+ sponge_in: cli.sponge_in,
sponge_mult: cli.sponge_mult,
collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc },
kbc_model: math::KbcModel::from_str(&cli.kbc_model).ok_or("неизвестная модель KBC")?,
@@ -482,6 +643,13 @@ impl Backend {
Backend::Gpu(s) => s.sample_field(k),
}
}
+ fn sample_centerline(&mut self) -> (Vec, Vec) {
+ match self {
+ Backend::Cpu(s) => s.sample_centerline(),
+ #[cfg(feature = "gpu")]
+ Backend::Gpu(s) => s.sample_centerline(),
+ }
+ }
fn is_finite(&mut self) -> bool {
match self {
Backend::Cpu(s) => s.is_finite(),
@@ -598,6 +766,19 @@ fn run(cli: Cli) -> Result<(), String> {
None => None,
};
+ // ── x–t диагностика ──
+ let mut xt_file = match &cli.xt {
+ Some(path) => {
+ let mut f = std::io::BufWriter::new(
+ std::fs::File::create(path).map_err(|e| format!("не создать {path}: {e}"))?,
+ );
+ writeln!(f, "# срез вдоль осевой линии, столбцы: step,t_phys_s,field,v0..v{}", spec.nx - 1)
+ .map_err(|e| e.to_string())?;
+ Some(f)
+ }
+ None => None,
+ };
+
// ── цикл ──
let d_ref = back.force_ref_size();
let mut recs: Vec = Vec::with_capacity(spec.steps as usize + 1);
@@ -611,16 +792,24 @@ fn run(cli: Cli) -> Result<(), String> {
})
.unwrap_or(0)) as f64;
+ let series_every = cli.series_every.max(1);
+ let mut scratch: Vec = Vec::with_capacity(256);
+ let mut rec = StepRec::default();
for t in 0..=spec.steps {
- back.advance(&mut recs);
+ back.advance(&mut scratch);
let need_frame = writer.is_some() && t % plan.stride == 0;
let need_report = cli.report_every > 0 && !quiet && t % cli.report_every == 0;
// и кадр, и живая строка требуют актуальных данных — синхронизируемся только здесь
if need_frame || need_report {
- back.flush(&mut recs);
+ back.flush(&mut scratch);
+ }
+ for r in scratch.drain(..) {
+ if r.step % series_every == 0 {
+ recs.push(r);
+ }
+ rec = r;
}
- let rec = *recs.last().unwrap_or(&StepRec::default());
if need_frame {
if let Some(w) = writer.as_mut() {
@@ -634,18 +823,36 @@ fn run(cli: Cli) -> Result<(), String> {
if need_report {
live_line(&spec, &rec, d_ref, t_start.elapsed().as_secs_f64(), nodes_per_step, full);
}
+ if let Some(f) = xt_file.as_mut() {
+ if t % cli.xt_every.max(1) == 0 {
+ let (rho, ux) = back.sample_centerline();
+ let tp = spec.units.time_of_step(t);
+ for (name, v) in [("rho", &rho), ("ux", &ux)] {
+ write!(f, "{t},{tp:.9e},{name}").map_err(|e| e.to_string())?;
+ for z in v {
+ write!(f, ",{z:.7e}").map_err(|e| e.to_string())?;
+ }
+ writeln!(f).map_err(|e| e.to_string())?;
+ }
+ }
+ }
// Развал счёта ловится по УЖЕ ПОСЧИТАННЫМ величинам, а не полным проходом по полю:
// на сетке 4096×2048 такой проход тянет с устройства сотни мегабайт, и делать его
// регулярно нельзя. Полная проверка — один раз, для подтверждения.
if !rec.max_u.is_finite() || !rec.rho_mean.is_finite() {
- back.flush(&mut recs);
+ back.flush(&mut scratch);
eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ на шаге {t}: в поле появились NaN/inf.");
blew_up = true;
break;
}
}
- back.flush(&mut recs);
+ back.flush(&mut scratch);
+ for r in scratch.drain(..) {
+ if r.step % series_every == 0 {
+ recs.push(r);
+ }
+ }
if !blew_up && !back.is_finite() {
eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ: итоговое поле содержит NaN/inf.");
blew_up = true;
@@ -690,7 +897,7 @@ fn default_range(k: FieldKind, u_lat: R, d_lat: R) -> gif::Range {
fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) {
let u = &spec.units;
let tau = 1.0 / (2.0 * spec.beta0);
- let blockage = spec.body.frontal_extent() / spec.ny as R;
+ let blockage = spec.scene.frontal_extent() / spec.ny as R;
println!("╔══════════════════════════════════════════════════════════════════════════╗");
println!("║ KBC-2D — обтекание тела в канале, D2Q9 + энтропийное столкновение KBC-D ║");
@@ -701,7 +908,12 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) {
println!(" размер ячейки {:>12.5} м домен {:.2} × {:.2} м",
u.dx, spec.nx as R * u.dx, spec.ny as R * u.dx);
println!(" тело {:>12} {:.3} м ({:.0} ячеек), угол атаки {:.1}°",
- format!("{:?}", spec.body.kind).to_lowercase(), u.d_phys(), u.d_lat, cli.body_angle);
+ if spec.scene.bodies.len() > 1 {
+ format!("{} тел", spec.scene.bodies.len())
+ } else {
+ format!("{:?}", spec.scene.bodies[0].kind).to_lowercase()
+ },
+ u.d_phys(), u.d_lat, cli.body_angle);
println!(" Рейнольдс {:>12.1} вязкость {:.4e} м²/с", u.re, u.nu_phys());
println!(" Мах (реш.) {:>12.4} u_lat = {:.4}, τ = {:.4}", u.mach(), u.u_lat, tau);
println!(" блокировка канала {:>12.3} (габарит тела / ширина канала)", blockage);
@@ -739,7 +951,7 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) {
// Действительно ли субсеточная граница работает: если большинство линков сваливается в
// простой отскок, стенка де-факто ступенчатая, как её ни называй.
{
- let g = cpu::Geom::build(spec.nx, spec.ny, &spec.body);
+ let g = cpu::Geom::build(spec.nx, spec.ny, &spec.scene);
let w = g.wall_stats();
let pc = |k: usize| 100.0 * k as R / w.links.max(1) as R;
println!("
@@ -787,8 +999,11 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) {
print!(" — короче периода моды, будет звон!");
}
}
- if spec.sponge_len > 0 {
- println!(" + губка {} столбцов ×{:.0}", spec.sponge_len, spec.sponge_mult);
+ if spec.sponge_len > 0 || spec.sponge_in > 0 {
+ println!(
+ " + губка: вход {} / выход {} столбцов ×{:.0}",
+ spec.sponge_in, spec.sponge_len, spec.sponge_mult
+ );
} else {
println!("; губка ВЫКЛЮЧЕНА");
}
@@ -899,7 +1114,7 @@ fn final_report(
}
// ── установившийся режим ──
- let beta = spec.body.frontal_extent() / spec.ny as R;
+ let beta = spec.scene.frontal_extent() / spec.ny as R;
println!("\n── Установившийся режим (вторая половина ряда, шаги {}–{}) ──────────────────", h, n - 1);
println!(" {:<22}{:>10}{:>10}{:>10}{:>10}", "величина", "St", "⟨Cd⟩", "rms Cl", "⟨Cm⟩");
println!(" {:<22}{:>10.4}{:>10.4}{:>10.4}{:>10.5}", "как посчитано", st, cd_m, cl_rms, cm_m);
@@ -912,13 +1127,55 @@ fn final_report(
cl_rms * (1.0 - beta).powi(2),
"—"
);
- if spec.body.kind == ShapeKind::Cylinder && u.re > 100.0 && u.re < 200.0 {
+ if spec.scene.bodies.len() == 1
+ && spec.scene.bodies[0].kind == ShapeKind::Cylinder
+ && u.re > 100.0
+ && u.re < 200.0
+ {
println!(" {:<22}{:>10.4}{:>10.4}{:>10}{:>10.1}", "литература Re≈150", 0.183, 1.330, "~0.30", 0.0);
}
println!(" rms u_y в зонде = {uy_rms:.5} (узел {:?}); блокировка β = {beta:.3}", spec.probe);
println!(" поправка: St×(1−β), Cd и rms Cl ×(1−β)². Литература — безграничный цилиндр;");
println!(" ⟨Cm⟩≈0 для симметричного тела под нулевым углом — это контроль симметрии считывания.");
+ // ── силы по телам ──
+ // Для тандема и многоэлементных конфигураций литература даёт коэффициенты КАЖДОГО тела
+ // по отдельности, поэтому общей суммы недостаточно.
+ let nb = spec.scene.bodies.len().min(MAX_BODY_BUCKETS);
+ if nb > 1 {
+ println!("\n── Силы по телам (вторая половина ряда) ────────────────────────────────────");
+ println!(" {:>5}{:>12}{:>10}{:>10}{:>11}", "тело", "положение", "⟨Cd⟩", "rms Cl", "⟨Cm⟩");
+ for b in 0..nb {
+ let cdb: Vec =
+ recs[h..].iter().map(|r| 2.0 * r.body[b][0] / (u.u_lat * u.u_lat * d_ref)).collect();
+ let clb: Vec =
+ recs[h..].iter().map(|r| 2.0 * r.body[b][1] / (u.u_lat * u.u_lat * d_ref)).collect();
+ let cmb: Vec = recs[h..]
+ .iter()
+ .map(|r| 2.0 * r.body[b][2] / (u.u_lat * u.u_lat * d_ref * d_ref))
+ .collect();
+ let (m, _) = math::mean_std(&cdb);
+ let (_, sl) = math::mean_std(&clb);
+ let (mm, _) = math::mean_std(&cmb);
+ let pos = spec
+ .scene
+ .bodies
+ .get(b)
+ .map(|bd| format!("{:.0},{:.0}", bd.cx, bd.cy))
+ .unwrap_or_else(|| "—".into());
+ println!(" {:>5}{:>12}{:>10.4}{:>10.4}{:>11.5}", b + 1, pos, m, sl, mm);
+ }
+ if spec.scene.bodies.len() > MAX_BODY_BUCKETS {
+ println!(
+ " (тел {}, вёдер {} — тела с {}-го и дальше сведены в последнюю строку;",
+ spec.scene.bodies.len(),
+ MAX_BODY_BUCKETS,
+ MAX_BODY_BUCKETS
+ );
+ println!(" сумма по строкам при этом точна и равна общему ⟨Cd⟩ выше)");
+ }
+ }
+
// ── сходимость ──
let nw = cli.windows.max(2).min(n / 2);
let w = n / nw;
@@ -970,7 +1227,9 @@ fn final_report(
── Продольная пульсация (последняя десятая часть прогона) ───────────────────");
println!(" размах ⟨ρ⟩ = {:.5} вокруг {mean:.5} ⇒ амплитуда скорости u' = {:.2}% от U",
hi - lo, ratio * 100.0);
- if ratio > 0.05 {
+ if !ratio.is_finite() {
+ println!(" ⚠ ПУЛЬСАЦИЯ НЕ ОПРЕДЕЛЕНА: в рядах NaN/inf, счёт развалился.");
+ } else if ratio.abs() > 0.05 {
println!(" ⚠ ПОТОК ЗАМЕТНО ПУЛЬСИРУЕТ: это продольная мода канала, а не физика следа.");
println!(" Лечится стартом из однородного потока (--init uniform) и губкой");
println!(" (--sponge-len {}). Числа Cd/St на таком прогоне недостоверны.",
diff --git a/docs/theory/2d_solver/src/math.rs b/docs/theory/2d_solver/src/math.rs
index e1be810..6b9b187 100644
--- a/docs/theory/2d_solver/src/math.rs
+++ b/docs/theory/2d_solver/src/math.rs
@@ -238,6 +238,58 @@ pub fn collide_node_bgk(f: &mut [R; Q], beta: R) -> Kbc {
Kbc { gamma: 2.0, degenerate: false }
}
+/// Множитель скорости у стенки на СТАРТЕ, по знаковому расстоянию до тела.
+///
+/// Зачем. Однородный старт заливает набегающим потоком весь домен, включая то место, где
+/// стоит тело; на первом же шаге стенка резко останавливает эту жидкость, и рождается
+/// импульс сжатия. Замерено x–t диагностикой: на нулевом шаге |ρ−1| ≈ 1.3·10⁻² прямо у
+/// задней кромки, тогда как у входа и выхода 10⁻⁸. Дальше импульс уходит полосой через весь
+/// домен, а на высоких Re успевает раскачать неустойчивость.
+///
+/// Лечение — свести скорость к нулю на подходе к телу за `width` клеток. Разрыв исчезает,
+/// а вдали от тела поле остаётся ровно однородным.
+#[inline]
+pub fn wall_taper(phi: R, width: R) -> R {
+ if width <= 0.0 {
+ return if phi > 0.0 { 1.0 } else { 0.0 };
+ }
+ smoothstep(phi / width)
+}
+
+/// Профиль β по столбцам с поглощающими губками у входа и у выхода.
+///
+/// Губка — плавный (smoothstep) подъём вязкости к границе, в `mult` раз. Зачем она нужна:
+/// пара «вход по скорости / выход по давлению» акустически есть четвертьволновая труба, вход
+/// отражает продольные волны как жёсткий поршень, а собственное затухание моды идёт как
+/// ν(π/2Nx)², то есть на длинном домене его практически нет. Губка — единственное, что эту
+/// моду реально ест.
+///
+/// Там, где губки перекрываются, берётся более сильная — они не складываются.
+pub fn beta_profile(nx: usize, beta0: R, sponge_in: usize, sponge_out: usize, mult: R) -> Vec {
+ let nu = nu_of_beta(beta0);
+ (0..nx)
+ .map(|x| {
+ let a = if sponge_out > 0 {
+ let start = nx - 1 - sponge_out;
+ smoothstep((x as R - start as R) / sponge_out as R)
+ } else {
+ 0.0
+ };
+ let b = if sponge_in > 0 {
+ smoothstep((sponge_in as R - x as R) / sponge_in as R)
+ } else {
+ 0.0
+ };
+ let s = a.max(b);
+ if s <= 0.0 {
+ beta0
+ } else {
+ beta_of_nu(nu * (1.0 + (mult - 1.0) * s))
+ }
+ })
+ .collect()
+}
+
/// Кинематическая вязкость по β, формула (5): ν = c_s²(1/(2β) − 1/2).
#[inline]
pub fn nu_of_beta(beta: R) -> R {
@@ -357,6 +409,9 @@ pub enum ShapeKind {
Triangle,
/// Тонкая пластина: длина D, толщина ratio·D.
Plate,
+ /// Произвольный многоугольник, вершины задаются снаружи (в клетках, от центра тела).
+ /// Точный SDF уже есть — открывает клинья, зазубренные кромки, любые обводы.
+ Polygon,
}
impl ShapeKind {
@@ -369,11 +424,12 @@ impl ShapeKind {
"naca" | "airfoil" | "профиль" => ShapeKind::Naca,
"triangle" | "треугольник" => ShapeKind::Triangle,
"plate" | "пластина" => ShapeKind::Plate,
+ "polygon" | "poly" | "многоугольник" => ShapeKind::Polygon,
_ => return None,
})
}
- pub const ALL: [&'static str; 7] =
- ["cylinder", "square", "diamond", "ellipse", "naca", "triangle", "plate"];
+ pub const ALL: [&'static str; 8] =
+ ["cylinder", "square", "diamond", "ellipse", "naca", "triangle", "plate", "polygon"];
}
/// Тело: форма + положение + масштаб + поворот. Для многоугольных форм контур считается один раз.
@@ -389,26 +445,67 @@ pub struct Body {
pub angle: R,
/// Относительная толщина для Ellipse/Naca/Plate.
pub ratio: R,
+ /// Относительная кривизна средней линии профиля (первая цифра NACA/100).
+ pub camber: R,
+ /// Положение максимума кривизны в долях хорды (вторая цифра NACA/10).
+ pub camber_pos: R,
/// Контур в системе тела (без поворота), для многоугольных форм.
poly: Vec<[R; 2]>,
}
impl Body {
pub fn new(kind: ShapeKind, cx: R, cy: R, d: R, angle_deg: R, ratio: R) -> Self {
+ Self::full(kind, cx, cy, d, angle_deg, ratio, 0.0, 0.0, &[])
+ }
+
+ /// Полный конструктор: с кривизной средней линии (для NACA) и явными вершинами
+ /// (для Polygon; координаты в клетках, отсчитываются от центра тела).
+ #[allow(clippy::too_many_arguments)]
+ pub fn full(
+ kind: ShapeKind,
+ cx: R,
+ cy: R,
+ d: R,
+ angle_deg: R,
+ ratio: R,
+ camber: R,
+ camber_pos: R,
+ verts: &[[R; 2]],
+ ) -> Self {
let angle = angle_deg * PI / 180.0;
- let poly = build_polygon(kind, d, ratio);
- Body { kind, cx, cy, d, angle, ratio, poly }
+ let poly = if kind == ShapeKind::Polygon {
+ verts.to_vec()
+ } else {
+ build_polygon(kind, d, ratio, camber, camber_pos)
+ };
+ Body { kind, cx, cy, d, angle, ratio, camber, camber_pos, poly }
+ }
+
+ /// Разбор четырёхзначного обозначения NACA (например 4412) в кривизну и толщину.
+ pub fn naca_code(code: &str) -> Option<(R, R, R)> {
+ let b = code.trim().as_bytes();
+ if b.len() != 4 || !b.iter().all(|c| c.is_ascii_digit()) {
+ return None;
+ }
+ let m = (b[0] - b'0') as R / 100.0;
+ let pos = (b[1] - b'0') as R / 10.0;
+ let t = ((b[2] - b'0') * 10 + (b[3] - b'0')) as R / 100.0;
+ Some((m, pos, t))
}
/// Тот же контур, пересчитанный на уровень с измельчением `r` (координаты и размер ×r).
pub fn refined(&self, r: R, ox: R, oy: R) -> Body {
- Body::new(
+ let verts: Vec<[R; 2]> = self.poly.iter().map(|v| [v[0] * r, v[1] * r]).collect();
+ Body::full(
self.kind,
(self.cx - ox) * r,
(self.cy - oy) * r,
self.d * r,
self.angle * 180.0 / PI,
self.ratio,
+ self.camber,
+ self.camber_pos,
+ &verts,
)
}
@@ -453,11 +550,78 @@ impl Body {
}
}
+/// Сколько тел разносится по отдельным вёдрам силы. Больше — только суммарно.
+pub const MAX_BODY_BUCKETS: usize = 4;
+
+/// Несколько тел в домене.
+///
+/// SDF сцены — минимум по телам, поэтому маски, SDF-поле и Bouzidi-линки строятся ровно теми
+/// же процедурами, что и для одного тела: `min` из знаковых расстояний сам по себе есть
+/// знаковое расстояние до объединения (вне тел — точно, внутри — с занижением у стыков,
+/// что на границе, где живёт Bouzidi, не проявляется).
+#[derive(Clone, Debug)]
+pub struct Scene {
+ pub bodies: Vec,
+}
+
+impl Scene {
+ pub fn single(b: Body) -> Scene {
+ Scene { bodies: vec![b] }
+ }
+
+ #[inline]
+ pub fn sdf(&self, x: R, y: R) -> R {
+ let mut m = R::INFINITY;
+ for b in &self.bodies {
+ let v = b.sdf(x, y);
+ if v < m {
+ m = v;
+ }
+ }
+ m
+ }
+
+ /// Индекс ближайшего тела — им помечается линк, чтобы сила считалась по телам.
+ #[inline]
+ pub fn nearest(&self, x: R, y: R) -> usize {
+ let mut best = 0usize;
+ let mut m = R::INFINITY;
+ for (i, b) in self.bodies.iter().enumerate() {
+ let v = b.sdf(x, y);
+ if v < m {
+ m = v;
+ best = i;
+ }
+ }
+ best
+ }
+
+ /// Характерный размер для нормировки коэффициентов — размер ПЕРВОГО тела.
+ pub fn ref_size(&self) -> R {
+ self.bodies.first().map(|b| b.d).unwrap_or(1.0)
+ }
+
+ /// Центр первого тела: вокруг него считается момент.
+ pub fn center(&self) -> (R, R) {
+ self.bodies.first().map(|b| (b.cx, b.cy)).unwrap_or((0.0, 0.0))
+ }
+
+ /// Суммарный габарит поперёк потока — для оценки блокировки канала.
+ pub fn frontal_extent(&self) -> R {
+ self.bodies.iter().map(|b| b.frontal_extent()).fold(0.0, R::max)
+ }
+
+ pub fn refined(&self, r: R, ox: R, oy: R) -> Scene {
+ Scene { bodies: self.bodies.iter().map(|b| b.refined(r, ox, oy)).collect() }
+ }
+}
+
/// Контур тела в его собственной системе координат (центр в нуле, без поворота).
-fn build_polygon(kind: ShapeKind, d: R, ratio: R) -> Vec<[R; 2]> {
+fn build_polygon(kind: ShapeKind, d: R, ratio: R, camber: R, camber_pos: R) -> Vec<[R; 2]> {
let h = 0.5 * d;
match kind {
ShapeKind::Cylinder | ShapeKind::Ellipse => Vec::new(), // аналитические
+ ShapeKind::Polygon => Vec::new(), // вершины задаются снаружи
ShapeKind::Square => vec![[-h, -h], [h, -h], [h, h], [-h, h]],
ShapeKind::Diamond => vec![[h, 0.0], [0.0, h], [-h, 0.0], [0.0, -h]],
ShapeKind::Plate => {
@@ -475,36 +639,75 @@ fn build_polygon(kind: ShapeKind, d: R, ratio: R) -> Vec<[R; 2]> {
})
.collect()
}
- ShapeKind::Naca => naca_symmetric(d, ratio, 64),
+ ShapeKind::Naca => naca_profile(d, ratio, camber, camber_pos, 80),
}
}
-/// Симметричный профиль NACA00xx: y_t = 5t·c·(0.2969√ξ − 0.1260ξ − 0.3516ξ² + 0.2843ξ³ − 0.1036ξ⁴),
-/// ξ = x/c. Хорда направлена по +x тела, начало отсчёта смещено так, чтобы центр вращения
-/// (точка приложения угла атаки) был в c/4 — стандартная аэродинамическая четверть хорды.
-fn naca_symmetric(chord: R, t: R, n: usize) -> Vec<[R; 2]> {
+/// Профиль NACA четырёхзначной серии.
+///
+/// Толщина: y_t = 5t·c·(0.2969√ξ − 0.1260ξ − 0.3516ξ² + 0.2843ξ³ − 0.1036ξ⁴), ξ = x/c.
+/// Средняя линия (m — максимальная кривизна в долях хорды, p — её положение):
+/// y_c = (m/p²)(2pξ − ξ²) при ξ ≤ p
+/// y_c = (m/(1−p)²)((1−2p) + 2pξ − ξ²) при ξ > p
+/// Поверхности откладываются по НОРМАЛИ к средней линии, а не по вертикали — иначе у
+/// заметно изогнутого профиля толщина у носка получается завышенной.
+///
+/// При m = 0 получается симметричный 00xx: средняя линия вырождается в прямую, и деление
+/// на p² не выполняется вовсе.
+///
+/// Начало системы тела — в четверти хорды: это стандартная точка приложения угла атаки,
+/// вокруг неё и считается момент.
+fn naca_profile(chord: R, t: R, m: R, p: R, n: usize) -> Vec<[R; 2]> {
let yt = |xi: R| {
5.0 * t
* chord
* (0.2969 * xi.sqrt() - 0.1260 * xi - 0.3516 * xi * xi + 0.2843 * xi.powi(3)
- 0.1036 * xi.powi(4))
};
+ let cambered = m > 0.0 && p > 0.0 && p < 1.0;
+ let yc = |xi: R| -> R {
+ if !cambered {
+ return 0.0;
+ }
+ chord
+ * if xi <= p {
+ (m / (p * p)) * (2.0 * p * xi - xi * xi)
+ } else {
+ (m / ((1.0 - p) * (1.0 - p))) * ((1.0 - 2.0 * p) + 2.0 * p * xi - xi * xi)
+ }
+ };
+ let dyc = |xi: R| -> R {
+ if !cambered {
+ return 0.0;
+ }
+ if xi <= p {
+ (2.0 * m / (p * p)) * (p - xi)
+ } else {
+ (2.0 * m / ((1.0 - p) * (1.0 - p))) * (p - xi)
+ }
+ };
// косинусное сгущение к носку и хвосту
- let xs: Vec = (0..=n)
- .map(|k| 0.5 * (1.0 - (PI * k as R / n as R).cos()))
- .collect();
- let x0 = 0.25 * chord; // четверть хорды в нуле системы тела
+ let xs: Vec = (0..=n).map(|k| 0.5 * (1.0 - (PI * k as R / n as R).cos())).collect();
+ let x0 = 0.25 * chord;
let mut poly = Vec::with_capacity(2 * n);
+ let point = |xi: R, upper: bool| -> [R; 2] {
+ let th = dyc(xi).atan();
+ let (st, ct) = th.sin_cos();
+ let (tt, cc) = (yt(xi), yc(xi));
+ let (x, y) = if upper {
+ (xi * chord - tt * st, cc + tt * ct)
+ } else {
+ (xi * chord + tt * st, cc - tt * ct)
+ };
+ // носком против потока: профиль строится носком в −x, поэтому обе координаты
+ // отражаются по x относительно четверти хорды
+ [-(x - x0), y]
+ };
for &xi in xs.iter() {
- // нижняя поверхность, от носка к хвосту (обход по часовой → замкнётся против)
- poly.push([xi * chord - x0, -yt(xi)]);
+ poly.push(point(xi, false));
}
for &xi in xs.iter().rev().skip(1) {
- poly.push([xi * chord - x0, yt(xi)]);
- }
- // профиль строился носком в −x; развернём, чтобы носок смотрел против потока (+x → навстречу)
- for p in poly.iter_mut() {
- p[0] = -p[0];
+ poly.push(point(xi, true));
}
poly
}
@@ -617,8 +820,10 @@ pub struct Link {
pub kind: LinkKind,
/// Доля пересечения q = |x_f→стенка| / |x_f→x_solid| ∈ (0,1), из SDF, клип [0.02, 0.98].
pub q: R,
- /// Линк принадлежит обтекаемому телу (а не иной твёрдой поверхности) — по нему считается сила.
- pub body: bool,
+ /// Индекс тела, которого касается линк, зажатый в [0, MAX_BODY_BUCKETS). Сила считается
+ /// по каждому такому ведру отдельно; тела за пределом ведёрного набора сваливаются в
+ /// последнее, так что СУММА по вёдрам остаётся точной, а разбивка — только для первых.
+ pub body: u8,
}
/// Доля пересечения по SDF: φ_f / (φ_f − φ_solid). Клип отсекает вырожденные линки,