NotBigGhost/CFDManager
1
0
Fork 0
Code Issues 6 Pull Requests 1 Actions Packages Projects Releases Wiki Activity

Несколько тел в домене, произвольная геометрия и разбор стартовой акустики

Browse Source 
НЕСКОЛЬКО ТЕЛ (--bodies). SDF сцены — минимум по телам, поэтому маски, поле SDF и
Bouzidi-линки строятся ровно теми же процедурами, что и для одного тела. Каждому линку
проставляется индекс тела, и сила считается по телам отдельно: для тандема и многоэлементных
конфигураций литература даёт коэффициенты каждого тела, общей суммы недостаточно. Разнос идёт
по четырём вёдрам; тела за пределом набора сваливаются в последнее, так что сумма по вёдрам
всегда точна. Синтаксис: "cylinder:d=24,x=120,y=120; cylinder:d=24,x=192,y=120".

Проверка на тандеме L/D=3 при Re=150 даёт физически правильную картину: передний цилиндр
Cd=1.22 при rms Cl=0.02, задний Cd=-0.10 при rms Cl=0.19. Отрицательное сопротивление заднего
— классический режим экранирования, когда сдвиговые слои переднего замыкаются на задний, а
восьмикратный рост rms Cl отвечает тому, что задний треплет след переднего.

ПРОИЗВОЛЬНЫЙ МНОГОУГОЛЬНИК (--shape polygon --poly) — точный SDF уже был, добавлен разбор.
Открывает клинья, зазубренные кромки, любые обводы. Проверено на клине: 100% линков идут по
интерполяционной формуле.

ПРОФИЛИ С ИЗГИБОМ (--naca 4412). Средняя линия по четырёхзначной серии, поверхности
откладываются по НОРМАЛИ к ней, а не по вертикали — иначе у заметно изогнутого профиля толщина
у носка завышается. При нулевой кривизне вырождается в прежний симметричный 00xx.

ПРОРЕЖИВАНИЕ РЯДОВ (--series-every): на 5·10^6 шагов полный ряд занимал бы сотни мегабайт.

РАЗБОР СТАРТОВОЙ АКУСТИКИ. Добавлена x-t диагностика (--xt): срез ⟨ρ⟩ и u_x вдоль осевой
линии, по строке на срез. Что она показала:

* возмущение рождается НА ТЕЛЕ, а не на входе. На нулевом шаге |ρ−1| ≈ 1.3·10^-2 у задней
  кромки при 10^-8 у входа и выхода; причина — поле стартует однородным потоком сквозь то
  место, где стоит тело;
* дальше импульс уходит полосой через весь домен — ровно то, что наблюдалось на низких Re;
* на Re ≥ 1000 при выключенной губке он раскачивает неустойчивость у ВЫХОДА, и счёт гибнет
  около шага 17000. С губкой (умолчание) Re=2000 доживает.

Попытка лечения сглаживанием стартовой скорости у тела (--init-taper) ЗАМЕРЕНА И ОТВЕРГНУТА:
ширина 24 клетки давит возмущение нулевого шага восьмикратно (2.54·10^-2 → 3.13·10^-3), но пик
за прогон при этом даже растёт (2.54·10^-2 → 3.22·10^-2). Возмущение просто переносится во
времени. Ключ оставлен выключенным, замер записан в его описании.

Добавлена губка у входа (--sponge-in); профиль β вынесен в общую math::beta_profile, раньше
он дублировался в обоих бэкендах. Исправлен вердикт о пульсации: при развале счёта отношение
уходило в минус и проверка «> 0.05» молча не срабатывала.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
NotBigGhost
2026-08-15 12:20:07 +03:00
co-authored by Claude Opus 5
parent 039e5f4806
commit b97686c8db
4 changed files with 737 additions and 160 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/theory/2d_solver/src/cpu.rs
+105 -58
View File
@@ -8,7 +8,7 @@
use rayon::prelude::*; 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}; use crate::{Collision, FieldKind, Spec, StepRec};
// ───────────────────────────────────────────────────────────────────────────── // ─────────────────────────────────────────────────────────────────────────────
@@ -30,21 +30,30 @@ pub struct Geom {
} }
impl Geom { impl Geom {
/// Собрать геометрию уровня по телу. `solid` определяется знаком SDF: φ ≤ 0 — тело. /// Собрать геометрию уровня по сцене. `solid` определяется знаком SDF сцены (минимум по
pub fn build(nx: usize, ny: usize, body: &Body) -> Self { /// телам): φ ≤ 0 — тело. Каждому линку проставляется индекс тела, которого он касается,
/// чтобы сила считалась по телам отдельно.
pub fn build(nx: usize, ny: usize, scene: &Scene) -> Self {
let n = nx * ny; let n = nx * ny;
let mut phi = vec![0.0; n]; let mut phi = vec![0.0; n];
let mut solid = vec![false; n]; let mut solid = vec![false; n];
let mut owner = vec![0u8; n];
for y in 0..ny { for y in 0..ny {
for x in 0..nx { for x in 0..nx {
let p = body.sdf(x as R, y as R); let k = y * nx + x;
phi[y * nx + x] = p; let p = scene.sdf(x as R, y as R);
solid[y * nx + x] = p <= 0.0; 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); 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. Делается один раз. /// формул Bouzidi применима, и считаем долю пересечения q из SDF. Делается один раз.
fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R]) -> Vec<Link> { fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R], owner: &[u8]) -> Vec<Link> {
let mut links = Vec::new(); let mut links = Vec::new();
for y in 0..ny { for y in 0..ny {
for x in 0..nx { for x in 0..nx {
@@ -88,7 +97,7 @@ fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R]) -> Vec<Link> {
ib: OPP[i] as u8, ib: OPP[i] as u8,
kind, kind,
q, q,
body: true, body: owner[s],
}); });
} }
} }
@@ -435,26 +444,26 @@ impl Level {
} }
} }
/// Сила и момент на теле по GMEM, суммой по линкам тела. /// Сила и момент по GMEM, разложенные по телам: `out[b] = (F_x, F_y, T_z)` для тела b.
fn force(&self) -> (R, R, R) { /// Момент считается вокруг центра ПЕРВОГО тела — общей точки отсчёта для всей сцены.
let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0); fn force(&self) -> [[R; 3]; MAX_BODY_BUCKETS] {
let mut out = [[0.0; 3]; MAX_BODY_BUCKETS];
let nx = self.nx; let nx = self.nx;
for l in &self.geom.links { for l in &self.geom.links {
if !l.body {
continue;
}
let node = l.node as usize; let node = l.node as usize;
let i = l.i 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); let (dfx, dfy) =
fx += dfx; math::gmem_link(i, self.post[node][i], self.f[node][l.ib as usize], 0.0, 0.0);
fy += dfy;
// плечо — до ТОЧКИ ПЕРЕСЕЧЕНИЯ линка со стенкой r_w = r_f + q·c_i, а не до узла: // плечо — до ТОЧКИ ПЕРЕСЕЧЕНИЯ линка со стенкой r_w = r_f + q·c_i, а не до узла:
// узел дал бы ошибку плеча до целой ячейки, а q уже посчитан // узел дал бы ошибку плеча до целой ячейки, а q уже посчитан
let rx = (node % nx) as R - self.geom.body_cx + l.q * CX[i] as R; 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; 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] #[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 за шаг // Статистика KBC за шаг
// ───────────────────────────────────────────────────────────────────────────── // ─────────────────────────────────────────────────────────────────────────────
@@ -730,22 +753,15 @@ pub struct Sim {
impl Sim { impl Sim {
pub fn new(spec: Spec) -> Sim { pub fn new(spec: Spec) -> Sim {
let (nx, ny) = (spec.nx, spec.ny); 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 = math::beta_profile(
// «скорость-вход + давление-выход» — недодемпфированный акустический резонатор, nx,
// губка гасит и вихри, и акустику до прихода на выход. spec.beta0,
let beta0: Vec<R> = (0..nx) spec.sponge_in,
.map(|x| { spec.sponge_len,
if spec.sponge_len == 0 { spec.sponge_mult,
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 u0 = if spec.init_uniform { let u0 = if spec.init_uniform {
@@ -754,23 +770,30 @@ impl Sim {
} else { } else {
(0.0, 0.0) (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 fluid_count = l0.geom.solid.iter().filter(|s| !**s).count() as R;
let (l1, patch) = if spec.refine > 1 { let (l1, patch) = if spec.refine > 1 {
let (ax, bx, ay, by) = spec.patch.expect("refine > 1 требует патч"); let (ax, bx, ay, by) = spec.patch.expect("refine > 1 требует патч");
let r = spec.refine; 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 nfx = r * (bx - ax) + 1;
let nfy = r * (by - ay) + 1; let nfy = r * (by - ay) + 1;
let geom1 = Geom::build(nfx, nfy, &body1); let geom1 = Geom::build(nfx, nfy, &scene1);
// τ_f = r(τ_c − ½) + ½ ⇒ одинаковая ν на обоих уровнях // τ_f = r(τ_c − ½) + ½ ⇒ одинаковая ν на обоих уровнях
let tau0 = 1.0 / (2.0 * spec.beta0); let tau0 = 1.0 / (2.0 * spec.beta0);
let tau1 = r as R * (tau0 - 0.5) + 0.5; let tau1 = r as R * (tau0 - 0.5) + 0.5;
let beta1 = 1.0 / (2.0 * tau1); let beta1 = 1.0 / (2.0 * tau1);
let r01 = tau1 / (r as R * tau0); let r01 = tau1 / (r as R * tau0);
let patch = Patch::new(&spec, &l0.geom, &geom1.solid, r01); 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 { } else {
(None, None) (None, None)
}; };
@@ -842,7 +865,7 @@ impl Sim {
self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap); self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap);
// ── уровень 1: r подшагов с временной интерполяцией рамки ── // ── уровень 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()) { 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.pre, &mut self.gh_old);
p.ghost_values(&self.l0.f, &mut self.gh_new); p.ghost_values(&self.l0.f, &mut self.gh_new);
@@ -852,23 +875,24 @@ impl Sim {
l1.apply_wall(wall); l1.apply_wall(wall);
// силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на // силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на
// последнем подшаге даёт лишний шум в рядах при том же среднем // последнем подшаге даёт лишний шум в рядах при том же среднем
let (a, b, c) = l1.force(); let g = l1.force();
fx += a; for k in 0..MAX_BODY_BUCKETS {
fy += b; for c in 0..3 {
tz += c; fb[k][c] += g[k][c];
}
}
let w = (s + 1) as R / p.r as R; let w = (s + 1) as R / p.r as R;
p.fill(&mut l1.f, &self.gh_old, &self.gh_new, w); p.fill(&mut l1.f, &self.gh_old, &self.gh_new, w);
} }
let inv = 1.0 / p.r as R; let inv = 1.0 / p.r as R;
fx *= inv; for k in 0..MAX_BODY_BUCKETS {
fy *= inv; for c in 0..3 {
tz *= inv; fb[k][c] *= inv;
}
}
p.restrict_to(&l1.f, &mut self.l0.f); p.restrict_to(&l1.f, &mut self.l0.f);
} else { } else {
let (a, b, c) = self.l0.force(); fb = self.l0.force();
fx = a;
fy = b;
tz = c;
} }
// ── диагностика ── // ── диагностика ──
@@ -892,11 +916,18 @@ impl Sim {
.reduce(|| (0.0, 0.0), |a, b| (a.0 + b.0, a.1.max(b.1))); .reduce(|| (0.0, 0.0), |a, b| (a.0 + b.0, a.1.max(b.1)));
self.step_index += 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 { StepRec {
step: t, step: t,
fx, fx,
fy, fy,
tz, tz,
body: fb,
uy_probe: probe.2, uy_probe: probe.2,
rho_mean: rho_sum / self.fluid_count, rho_mean: rho_sum / self.fluid_count,
max_u, max_u,
@@ -911,8 +942,8 @@ impl Sim {
/// Характерный размер тела в единицах того уровня, где снимается сила. /// Характерный размер тела в единицах того уровня, где снимается сила.
pub fn force_ref_size(&self) -> R { pub fn force_ref_size(&self) -> R {
match &self.patch { match &self.patch {
Some(p) => self.spec.body.d * p.r as R, Some(p) => self.spec.scene.ref_size() * p.r as R,
None => self.spec.body.d, None => self.spec.scene.ref_size(),
} }
} }
@@ -958,6 +989,22 @@ impl Sim {
(out, &self.l0.geom.solid) (out, &self.l0.geom.solid)
} }
/// Срез вдоль осевой линии: (ρ, u_x) по каждому столбцу. Из последовательности таких
/// срезов складывается x–t диаграмма, по которой видно, бежит возмущение со скоростью
/// звука или конвекции и есть ли стоячие узлы.
pub fn sample_centerline(&self) -> (Vec<R>, Vec<R>) {
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 — признак развала счёта. /// Есть ли в поле NaN/inf — признак развала счёта.
pub fn is_finite(&self) -> bool { pub fn is_finite(&self) -> bool {
self.l0.f.par_iter().all(|c| c.iter().all(|v| v.is_finite())) 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() { fn wall_is_subgrid_not_staircase() {
let mut means = Vec::new(); let mut means = Vec::new();
for off in [0.0, 0.25, 0.5] { 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 body = math::Body::new(math::ShapeKind::Cylinder, 60.0, 40.0 + off, 24.0, 0.0, 0.3);
let g = Geom::build(120, 80, &body); let g = Geom::build(120, 80, &Scene::single(body));
let w = g.wall_stats(); let w = g.wall_stats();
assert!(w.links > 100, "смещение {off}: линков всего {}", w.links); assert!(w.links > 100, "смещение {off}: линков всего {}", w.links);
assert_eq!(w.simple, 0, "смещение {off}: {} ступенчатых линков", w.simple); assert_eq!(w.simple, 0, "смещение {off}: {} ступенчатых линков", w.simple);
@@ -1069,8 +1116,8 @@ mod tests {
/// без пересечений, в том же порядке. /// без пересечений, в том же порядке.
#[test] #[test]
fn wall_node_index_partitions_links() { fn wall_node_index_partitions_links() {
let body = Body::new(math::ShapeKind::Naca, 60.0, 40.0, 30.0, 12.0, 0.25); let body = math::Body::new(math::ShapeKind::Naca, 60.0, 40.0, 30.0, 12.0, 0.25);
let g = Geom::build(120, 80, &body); let g = Geom::build(120, 80, &Scene::single(body));
let mut covered = 0usize; let mut covered = 0usize;
for w in &g.wall_nodes { for w in &g.wall_nodes {
assert_eq!(w.first as usize, covered, "разрыв в индексе граничных узлов"); assert_eq!(w.first as usize, covered, "разрыв в индексе граничных узлов");
docs/theory/2d_solver/src/gpu.rs
+125 -59
View File
@@ -102,7 +102,8 @@ struct GLink {
ib: u32, ib: u32,
kind: u32, kind: u32,
q: f32, q: f32,
_p: [u32; 2], body: u32,
_p: u32,
} }
#[repr(C)] #[repr(C)]
@@ -134,6 +135,8 @@ struct Results {
cnt: f32, cnt: f32,
degen: f32, degen: f32,
xineg: 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 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 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 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 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 }; 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). // на каждом шаге: результаты копятся и читаются пачкой (см. HIST в gpu.rs).
const HIST: u32 = 128u; 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) бит; здесь потеря не копится, // Наивная сумма по 10^5…10^7 значений в f32 съедает ~log2(N) бит; здесь потеря не копится,
@@ -404,9 +409,10 @@ fn k_force(@builtin(local_invocation_id) lid: vec3<u32>) {
let t = lid.x; let t = lid.x;
var cx = array<f32,9>(0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0); var cx = array<f32,9>(0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0);
var cy = array<f32,9>(0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0); var cy = array<f32,9>(0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0);
var sx = vec2<f32>(0.0, 0.0); // накопители на каждое тело держатся в регистрах, редукция потом идёт по телам подряд
var sy = vec2<f32>(0.0, 0.0); var pfx = array<f32,4>(0.0, 0.0, 0.0, 0.0);
var sz = vec2<f32>(0.0, 0.0); var pfy = array<f32,4>(0.0, 0.0, 0.0, 0.0);
var ptz = array<f32,4>(0.0, 0.0, 0.0, 0.0);
var k = t; var k = t;
loop { loop {
if (k >= P.nlinks) { break; } if (k >= P.nlinks) { break; }
@@ -415,31 +421,36 @@ fn k_force(@builtin(local_invocation_id) lid: vec3<u32>) {
let fb = f[L.ib*P.n + L.node]; let fb = f[L.ib*P.n + L.node];
let dfx = cx[L.i]*fp + cx[L.i]*fb; let dfx = cx[L.i]*fp + cx[L.i]*fb;
let dfy = cy[L.i]*fp + cy[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 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]; 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; k = k + 256u;
} }
wfx[t] = sx.x + sx.y; wfy[t] = sy.x + sy.y; wtz[t] = sz.x + sz.y; for (var b = 0u; b < MAXB; b = b + 1u) {
workgroupBarrier(); wfx[t] = pfx[b]; wfy[t] = pfy[b]; wtz[t] = ptz[b];
var s = 128u; workgroupBarrier();
loop { var s = 128u;
if (s == 0u) { break; } loop {
if (t < s) { if (s == 0u) { break; }
wfx[t] = wfx[t] + wfx[t + s]; if (t < s) {
wfy[t] = wfy[t] + wfy[t + s]; wfx[t] = wfx[t] + wfx[t + s];
wtz[t] = wtz[t] + wtz[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(); 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<u32>,
var fv: array<f32,9>; var fv: array<f32,9>;
load9(&fv, nd, P.n); load9(&fv, nd, P.n);
let m = macros9(fv); let m = macros9(fv);
results[D.slot*12u + 3u] = m.z; results[D.slot*SLOT + 3u] = m.z;
} }
wp[t] = p; wp[t] = p;
workgroupBarrier(); workgroupBarrier();
@@ -554,17 +565,26 @@ fn k_stats2(@builtin(local_invocation_id) lid: vec3<u32>) {
qgmin = min(qgmin, o.gmin); qgmin = min(qgmin, o.gmin);
qgmax = max(qgmax, o.gmax); qgmax = max(qgmax, o.gmax);
} }
var fx = 0.0; var fy = 0.0; var tz = 0.0; // сила усредняется по подшагам L1 и раскладывается по телам; общая — их сумма
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];
}
let inv = 1.0 / f32(D.refine); let inv = 1.0 / f32(D.refine);
let o = D.slot*12u; let o = D.slot*SLOT;
results[o + 0u] = fx*inv; var tfx = 0.0; var tfy = 0.0; var ttz = 0.0;
results[o + 1u] = fy*inv; for (var b = 0u; b < MAXB; b = b + 1u) {
results[o + 2u] = tz*inv; 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 + 4u] = qrho.x + qrho.y;
results[o + 5u] = qmaxu; results[o + 5u] = qmaxu;
results[o + 6u] = qgsum.x + qgsum.y; results[o + 6u] = qgsum.x + qgsum.y;
@@ -784,11 +804,30 @@ impl GpuLevel {
} }
} }
fn soa_equilibrium(n: usize, u0: (R, R)) -> Vec<f32> { /// Стартовое поле в раскладке SoA. Если задана сцена, скорость сводится к нулю на подходе
let fe = math::feq(1.0, u0.0, u0.1); /// к телу (см. `math::wall_taper`) — иначе на первом шаге возникает разрыв и импульс сжатия.
fn soa_equilibrium(
n: usize,
nx: usize,
u0: (R, R),
taper: Option<(&math::Scene, R)>,
) -> Vec<f32> {
let mut v = vec![0.0f32; 9 * n]; let mut v = vec![0.0f32; 9 * n];
for i in 0..9 { let plain = math::feq(1.0, u0.0, u0.1);
for k in 0..n { 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; v[i * n + k] = fe[i] as f32;
} }
} }
@@ -807,9 +846,10 @@ fn make_level(
flags: u32, flags: u32,
u0: (R, R), u0: (R, R),
wall_layout: &wgpu::BindGroupLayout, wall_layout: &wgpu::BindGroupLayout,
taper: Option<(&math::Scene, R)>,
) -> GpuLevel { ) -> GpuLevel {
let n = nx * ny; let n = nx * ny;
let init = soa_equilibrium(n, u0); let init = soa_equilibrium(n, nx, u0, taper);
let mkf = |label: &str| { let mkf = |label: &str| {
device.create_buffer_init(&wgpu::util::BufferInitDescriptor { device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
label: Some(label), label: Some(label),
@@ -841,7 +881,8 @@ fn make_level(
math::LinkKind::Simple => 2, math::LinkKind::Simple => 2,
}, },
q: l.q as f32, q: l.q as f32,
_p: [0; 2], body: l.body as u32,
_p: 0,
}) })
.collect(); .collect();
let links_buf = storage_init(device, "links", bytemuck::cast_slice(&links)); 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 fluid_count = geom0.solid.iter().filter(|s| !**s).count() as R;
let beta0: Vec<f32> = (0..spec.nx) let beta0: Vec<f32> =
.map(|x| { math::beta_profile(spec.nx, spec.beta0, spec.sponge_in, spec.sponge_len, spec.sponge_mult)
if spec.sponge_len == 0 { .into_iter()
return spec.beta0 as f32; .map(|v| v as f32)
} .collect();
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 module = device.create_shader_module(wgpu::ShaderModuleDescriptor { let module = device.create_shader_module(wgpu::ShaderModuleDescriptor {
label: Some("kbc2d.wgsl"), label: Some("kbc2d.wgsl"),
@@ -1122,6 +1157,7 @@ impl Sim {
flags0, flags0,
u0, u0,
&bgl_wall, &bgl_wall,
if spec.init_uniform { Some((&spec.scene, spec.init_taper)) } else { None },
); );
let pre = device.create_buffer(&wgpu::BufferDescriptor { let pre = device.create_buffer(&wgpu::BufferDescriptor {
@@ -1137,10 +1173,10 @@ impl Sim {
if spec.refine > 1 { if spec.refine > 1 {
let (ax, bx, ay, by) = spec.patch.ok_or("refine > 1 требует патч")?; let (ax, bx, ay, by) = spec.patch.ok_or("refine > 1 требует патч")?;
let r = spec.refine; 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 nfx = r * (bx - ax) + 1;
let nfy = r * (by - ay) + 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 tau0 = 1.0 / (2.0 * spec.beta0);
let tau1 = r as R * (tau0 - 0.5) + 0.5; let tau1 = r as R * (tau0 - 0.5) + 0.5;
let beta1 = vec![(1.0 / (2.0 * tau1)) as f32; nfx]; 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 flags1 = if probe_on_fine { 2 } else { 0 };
let lvl1 = let lvl1 =
make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1, 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<GGhost> = patch let ghosts: Vec<GGhost> = patch
.ghosts() .ghosts()
@@ -1262,8 +1303,8 @@ impl Sim {
}); });
let d_ref = match spec.patch { let d_ref = match spec.patch {
Some(_) if spec.refine > 1 => spec.body.d * spec.refine as R, Some(_) if spec.refine > 1 => spec.scene.ref_size() * spec.refine as R,
_ => spec.body.d, _ => spec.scene.ref_size(),
}; };
Ok(Sim { Ok(Sim {
@@ -1516,11 +1557,18 @@ impl Sim {
fn make_rec(&self, t: u64, res: &Results) -> StepRec { 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 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 { StepRec {
step: t, step: t,
fx: res.fx as R, fx: res.fx as R,
fy: res.fy as R, fy: res.fy as R,
tz: res.tz as R, tz: res.tz as R,
body,
uy_probe: res.uy_probe as R, uy_probe: res.uy_probe as R,
rho_mean: res.rho_sum as R / self.fluid_count, rho_mean: res.rho_sum as R / self.fluid_count,
max_u: res.max_u as R, max_u: res.max_u as R,
@@ -1564,6 +1612,24 @@ impl Sim {
out out
} }
/// Срез вдоль осевой линии: (ρ, u_x) по столбцам. На GPU это полное скачивание поля,
/// поэтому x–t диагностика включается редким шагом и только там, где нужна.
pub fn sample_centerline(&self) -> (Vec<R>, Vec<R>) {
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 { pub fn is_finite(&self) -> bool {
self.download_l0().iter().all(|v| v.is_finite()) 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) (out, solid)
} }
} }
docs/theory/2d_solver/src/main.rs
+278 -19
View File
@@ -19,7 +19,7 @@ use std::time::Instant;
use clap::Parser; 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)] #[derive(Clone, Copy, Debug, Default)]
pub struct StepRec { pub struct StepRec {
pub step: u64, pub step: u64,
/// Суммарная сила и момент по всей сцене.
pub fx: R, pub fx: R,
pub fy: R, pub fy: R,
pub tz: R, pub tz: R,
/// Она же по телам: `body[b] = (F_x, F_y, T_z)`. Тела за пределом набора вёдер
/// сваливаются в последнее, поэтому сумма по вёдрам всегда точна.
pub body: [[R; 3]; MAX_BODY_BUCKETS],
pub uy_probe: R, pub uy_probe: R,
pub rho_mean: R, pub rho_mean: R,
pub max_u: R, pub max_u: R,
@@ -75,7 +79,7 @@ pub struct StepRec {
pub struct Spec { pub struct Spec {
pub nx: usize, pub nx: usize,
pub ny: usize, pub ny: usize,
pub body: Body, pub scene: Scene,
pub refine: usize, pub refine: usize,
/// Границы патча измельчения в координатах L0: (ax, bx, ay, by). /// Границы патча измельчения в координатах L0: (ax, bx, ay, by).
pub patch: Option<(usize, usize, usize, usize)>, pub patch: Option<(usize, usize, usize, usize)>,
@@ -90,7 +94,11 @@ pub struct Spec {
pub outlet_extrapolate: bool, pub outlet_extrapolate: bool,
/// Начальное поле уже несёт набегающий поток (иначе — покой с разгоном входа). /// Начальное поле уже несёт набегающий поток (иначе — покой с разгоном входа).
pub init_uniform: bool, pub init_uniform: bool,
/// Ширина сведения стартовой скорости к нулю у тела, клеток.
pub init_taper: R,
pub sponge_len: usize, pub sponge_len: usize,
/// Длина губки после входа, столбцов.
pub sponge_in: usize,
pub sponge_mult: R, pub sponge_mult: R,
pub collision: Collision, pub collision: Collision,
/// Что входит в сдвиговую часть s (см. `math::KbcModel`). /// Что входит в сдвиговую часть s (см. `math::KbcModel`).
@@ -167,6 +175,20 @@ struct Cli {
/// Положение центра тела по y, ячеек (по умолчанию ny/2) /// Положение центра тела по y, ячеек (по умолчанию ny/2)
#[arg(long, help_heading = "Тело")] #[arg(long, help_heading = "Тело")]
body_y: Option<R>, body_y: Option<R>,
/// Четырёхзначное обозначение профиля NACA, например 4412 (кривизна 4%, её максимум на
/// 40% хорды, толщина 12%). Действует при --shape naca и перекрывает --thickness.
#[arg(long, help_heading = "Тело")]
naca: Option<String>,
/// Вершины произвольного многоугольника при --shape polygon: "x1,y1;x2,y2;..." в клетках,
/// отсчёт от центра тела. Обход любой, SDF точный.
#[arg(long, allow_hyphen_values = true, help_heading = "Тело")]
poly: Option<String>,
/// НЕСКОЛЬКО ТЕЛ вместо одного. Список через точку с запятой, каждое — форма и параметры:
/// "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<String>,
// ── время ── // ── время ──
/// Число шагов симуляции. Взаимоисключающе с --time. /// Число шагов симуляции. Взаимоисключающе с --time.
@@ -192,6 +214,17 @@ struct Cli {
#[arg(long, default_value = "uniform", value_parser = ["uniform", "rest"], #[arg(long, default_value = "uniform", value_parser = ["uniform", "rest"],
help_heading = "Время")] help_heading = "Время")]
init: String, 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 — выключить. Губка гасит продольную моду, которую вязкость сама не гасит. /// 0 — выключить. Губка гасит продольную моду, которую вязкость сама не гасит.
#[arg(long, help_heading = "Схема")] #[arg(long, help_heading = "Схема")]
sponge_len: Option<usize>, sponge_len: Option<usize>,
/// Длина губки ПОСЛЕ входа, столбцов (0 — выключена). Вход по скорости отражает продольные
/// волны как жёсткий поршень; губка у входа гасит их до отражения.
#[arg(long, default_value_t = 0, help_heading = "Схема")]
sponge_in: usize,
/// Во сколько раз губка поднимает вязкость /// Во сколько раз губка поднимает вязкость
#[arg(long, default_value_t = 30.0, help_heading = "Схема")] #[arg(long, default_value_t = 30.0, help_heading = "Схема")]
sponge_mult: R, sponge_mult: R,
@@ -273,9 +310,22 @@ struct Cli {
/// Число окон в отчёте о сходимости /// Число окон в отчёте о сходимости
#[arg(long, default_value_t = 10, help_heading = "Вывод")] #[arg(long, default_value_t = 10, help_heading = "Вывод")]
windows: usize, windows: usize,
/// Прореживание временных рядов: хранить каждую N-ю запись. Статистика и спектр от этого
/// не страдают, пока N много меньше периода схода вихрей, зато сверхдлинные прогоны
/// перестают съедать память (на 5·10⁶ шагов полный ряд — сотни мегабайт).
#[arg(long, default_value_t = 1, help_heading = "Вывод")]
series_every: u64,
/// Выгрузить временные ряды в CSV /// Выгрузить временные ряды в CSV
#[arg(long, help_heading = "Вывод")] #[arg(long, help_heading = "Вывод")]
csv: Option<String>, csv: Option<String>,
/// Файл x–t диаграммы: срез ⟨ρ⟩ и u_x вдоль осевой линии, по строке на срез. По наклону
/// полос видно, бежит возмущение со скоростью звука или конвекции, а стоячие узлы
/// проявляются сразу. Основной инструмент разбора стартовой акустики.
#[arg(long, help_heading = "Вывод")]
xt: Option<String>,
/// Период записи срезов 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)) Ok((a, b))
} }
/// Разбор одного тела из мини-синтаксиса `форма:ключ=значение,...`.
fn parse_body(spec: &str, def_x: R, def_y: R) -> Result<Body, String> {
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::<R>().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<Vec<[R; 2]>, 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<Spec, String> { fn build_spec(cli: &Cli) -> Result<Spec, String> {
if cli.nx < 16 || cli.ny < 16 { if cli.nx < 16 || cli.ny < 16 {
return Err("сетка меньше 16×16 не имеет смысла".into()); return Err("сетка меньше 16×16 не имеет смысла".into());
@@ -306,12 +415,51 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
return Err("--refine обязан быть ≥ 1".into()); 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 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 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 beta0 = math::beta_of_nu(units.nu_lat);
let tau0 = 1.0 / (2.0 * beta0); let tau0 = 1.0 / (2.0 * beta0);
if tau0 <= 0.5 { if tau0 <= 0.5 {
@@ -364,6 +512,17 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
// как открытый конец. Собственная мода затухает как ν(π/2Nx)², то есть практически не // как открытый конец. Собственная мода затухает как ν(π/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 { let sponge_len = match cli.sponge_len {
Some(v) => { Some(v) => {
if v > 0 { if v > 0 {
@@ -421,7 +580,7 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
Ok(Spec { Ok(Spec {
nx: cli.nx, nx: cli.nx,
ny: cli.ny, ny: cli.ny,
body, scene,
refine: cli.refine, refine: cli.refine,
patch, patch,
units, units,
@@ -433,7 +592,9 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
pert_dur: cli.pert_dur, pert_dur: cli.pert_dur,
outlet_extrapolate: cli.outlet == "extrapolate", outlet_extrapolate: cli.outlet == "extrapolate",
init_uniform: cli.init == "uniform", init_uniform: cli.init == "uniform",
init_taper: cli.init_taper.max(0.0),
sponge_len, sponge_len,
sponge_in: cli.sponge_in,
sponge_mult: cli.sponge_mult, sponge_mult: cli.sponge_mult,
collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc }, collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc },
kbc_model: math::KbcModel::from_str(&cli.kbc_model).ok_or("неизвестная модель 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), Backend::Gpu(s) => s.sample_field(k),
} }
} }
fn sample_centerline(&mut self) -> (Vec<R>, Vec<R>) {
match self {
Backend::Cpu(s) => s.sample_centerline(),
#[cfg(feature = "gpu")]
Backend::Gpu(s) => s.sample_centerline(),
}
}
fn is_finite(&mut self) -> bool { fn is_finite(&mut self) -> bool {
match self { match self {
Backend::Cpu(s) => s.is_finite(), Backend::Cpu(s) => s.is_finite(),
@@ -598,6 +766,19 @@ fn run(cli: Cli) -> Result<(), String> {
None => None, 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 d_ref = back.force_ref_size();
let mut recs: Vec<StepRec> = Vec::with_capacity(spec.steps as usize + 1); let mut recs: Vec<StepRec> = Vec::with_capacity(spec.steps as usize + 1);
@@ -611,16 +792,24 @@ fn run(cli: Cli) -> Result<(), String> {
}) })
.unwrap_or(0)) as f64; .unwrap_or(0)) as f64;
let series_every = cli.series_every.max(1);
let mut scratch: Vec<StepRec> = Vec::with_capacity(256);
let mut rec = StepRec::default();
for t in 0..=spec.steps { 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_frame = writer.is_some() && t % plan.stride == 0;
let need_report = cli.report_every > 0 && !quiet && t % cli.report_every == 0; let need_report = cli.report_every > 0 && !quiet && t % cli.report_every == 0;
// и кадр, и живая строка требуют актуальных данных — синхронизируемся только здесь // и кадр, и живая строка требуют актуальных данных — синхронизируемся только здесь
if need_frame || need_report { 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 need_frame {
if let Some(w) = writer.as_mut() { if let Some(w) = writer.as_mut() {
@@ -634,18 +823,36 @@ fn run(cli: Cli) -> Result<(), String> {
if need_report { if need_report {
live_line(&spec, &rec, d_ref, t_start.elapsed().as_secs_f64(), nodes_per_step, full); 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 такой проход тянет с устройства сотни мегабайт, и делать его // на сетке 4096×2048 такой проход тянет с устройства сотни мегабайт, и делать его
// регулярно нельзя. Полная проверка — один раз, для подтверждения. // регулярно нельзя. Полная проверка — один раз, для подтверждения.
if !rec.max_u.is_finite() || !rec.rho_mean.is_finite() { if !rec.max_u.is_finite() || !rec.rho_mean.is_finite() {
back.flush(&mut recs); back.flush(&mut scratch);
eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ на шаге {t}: в поле появились NaN/inf."); eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ на шаге {t}: в поле появились NaN/inf.");
blew_up = true; blew_up = true;
break; 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() { if !blew_up && !back.is_finite() {
eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ: итоговое поле содержит NaN/inf."); eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ: итоговое поле содержит NaN/inf.");
blew_up = true; 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) { fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) {
let u = &spec.units; let u = &spec.units;
let tau = 1.0 / (2.0 * spec.beta0); 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!("╔══════════════════════════════════════════════════════════════════════════╗");
println!("║ KBC-2D — обтекание тела в канале, D2Q9 + энтропийное столкновение KBC-D ║"); 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} м", println!(" размер ячейки {:>12.5} м домен {:.2} × {:.2} м",
u.dx, spec.nx as R * u.dx, spec.ny as R * u.dx); u.dx, spec.nx as R * u.dx, spec.ny as R * u.dx);
println!(" тело {:>12} {:.3} м ({:.0} ячеек), угол атаки {:.1}°", 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.1} вязкость {:.4e} м²/с", u.re, u.nu_phys());
println!(" Мах (реш.) {:>12.4} u_lat = {:.4}, τ = {:.4}", u.mach(), u.u_lat, tau); println!(" Мах (реш.) {:>12.4} u_lat = {:.4}, τ = {:.4}", u.mach(), u.u_lat, tau);
println!(" блокировка канала {:>12.3} (габарит тела / ширина канала)", blockage); 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 w = g.wall_stats();
let pc = |k: usize| 100.0 * k as R / w.links.max(1) as R; let pc = |k: usize| 100.0 * k as R / w.links.max(1) as R;
println!(" println!("
@@ -787,8 +999,11 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) {
print!(" — короче периода моды, будет звон!"); print!(" — короче периода моды, будет звон!");
} }
} }
if spec.sponge_len > 0 { if spec.sponge_len > 0 || spec.sponge_in > 0 {
println!(" + губка {} столбцов ×{:.0}", spec.sponge_len, spec.sponge_mult); println!(
" + губка: вход {} / выход {} столбцов ×{:.0}",
spec.sponge_in, spec.sponge_len, spec.sponge_mult
);
} else { } else {
println!("; губка ВЫКЛЮЧЕНА"); 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!("\n── Установившийся режим (вторая половина ряда, шаги {}–{}) ──────────────────", h, n - 1);
println!(" {:<22}{:>10}{:>10}{:>10}{:>10}", "величина", "St", "⟨Cd⟩", "rms Cl", "⟨Cm⟩"); 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); 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), 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!(" {:<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!(" rms u_y в зонде = {uy_rms:.5} (узел {:?}); блокировка β = {beta:.3}", spec.probe);
println!(" поправка: St×(1−β), Cd и rms Cl ×(1−β)². Литература — безграничный цилиндр;"); println!(" поправка: St×(1−β), Cd и rms Cl ×(1−β)². Литература — безграничный цилиндр;");
println!(" ⟨Cm⟩≈0 для симметричного тела под нулевым углом — это контроль симметрии считывания."); 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<R> =
recs[h..].iter().map(|r| 2.0 * r.body[b][0] / (u.u_lat * u.u_lat * d_ref)).collect();
let clb: Vec<R> =
recs[h..].iter().map(|r| 2.0 * r.body[b][1] / (u.u_lat * u.u_lat * d_ref)).collect();
let cmb: Vec<R> = 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 nw = cli.windows.max(2).min(n / 2);
let w = n / nw; let w = n / nw;
@@ -970,7 +1227,9 @@ fn final_report(
── Продольная пульсация (последняя десятая часть прогона) ───────────────────"); ── Продольная пульсация (последняя десятая часть прогона) ───────────────────");
println!(" размах ⟨ρ⟩ = {:.5} вокруг {mean:.5} ⇒ амплитуда скорости u' = {:.2}% от U", println!(" размах ⟨ρ⟩ = {:.5} вокруг {mean:.5} ⇒ амплитуда скорости u' = {:.2}% от U",
hi - lo, ratio * 100.0); hi - lo, ratio * 100.0);
if ratio > 0.05 { if !ratio.is_finite() {
println!(" ⚠ ПУЛЬСАЦИЯ НЕ ОПРЕДЕЛЕНА: в рядах NaN/inf, счёт развалился.");
} else if ratio.abs() > 0.05 {
println!(" ⚠ ПОТОК ЗАМЕТНО ПУЛЬСИРУЕТ: это продольная мода канала, а не физика следа."); println!(" ⚠ ПОТОК ЗАМЕТНО ПУЛЬСИРУЕТ: это продольная мода канала, а не физика следа.");
println!(" Лечится стартом из однородного потока (--init uniform) и губкой"); println!(" Лечится стартом из однородного потока (--init uniform) и губкой");
println!(" (--sponge-len {}). Числа Cd/St на таком прогоне недостоверны.", println!(" (--sponge-len {}). Числа Cd/St на таком прогоне недостоверны.",
docs/theory/2d_solver/src/math.rs
+229 -24
View File
@@ -238,6 +238,58 @@ pub fn collide_node_bgk(f: &mut [R; Q], beta: R) -> Kbc {
Kbc { gamma: 2.0, degenerate: false } 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<R> {
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). /// Кинематическая вязкость по β, формула (5): ν = c_s²(1/(2β) − 1/2).
#[inline] #[inline]
pub fn nu_of_beta(beta: R) -> R { pub fn nu_of_beta(beta: R) -> R {
@@ -357,6 +409,9 @@ pub enum ShapeKind {
Triangle, Triangle,
/// Тонкая пластина: длина D, толщина ratio·D. /// Тонкая пластина: длина D, толщина ratio·D.
Plate, Plate,
/// Произвольный многоугольник, вершины задаются снаружи (в клетках, от центра тела).
/// Точный SDF уже есть — открывает клинья, зазубренные кромки, любые обводы.
Polygon,
} }
impl ShapeKind { impl ShapeKind {
@@ -369,11 +424,12 @@ impl ShapeKind {
"naca" | "airfoil" | "профиль" => ShapeKind::Naca, "naca" | "airfoil" | "профиль" => ShapeKind::Naca,
"triangle" | "треугольник" => ShapeKind::Triangle, "triangle" | "треугольник" => ShapeKind::Triangle,
"plate" | "пластина" => ShapeKind::Plate, "plate" | "пластина" => ShapeKind::Plate,
"polygon" | "poly" | "многоугольник" => ShapeKind::Polygon,
_ => return None, _ => return None,
}) })
} }
pub const ALL: [&'static str; 7] = pub const ALL: [&'static str; 8] =
["cylinder", "square", "diamond", "ellipse", "naca", "triangle", "plate"]; ["cylinder", "square", "diamond", "ellipse", "naca", "triangle", "plate", "polygon"];
} }
/// Тело: форма + положение + масштаб + поворот. Для многоугольных форм контур считается один раз. /// Тело: форма + положение + масштаб + поворот. Для многоугольных форм контур считается один раз.
@@ -389,26 +445,67 @@ pub struct Body {
pub angle: R, pub angle: R,
/// Относительная толщина для Ellipse/Naca/Plate. /// Относительная толщина для Ellipse/Naca/Plate.
pub ratio: R, pub ratio: R,
/// Относительная кривизна средней линии профиля (первая цифра NACA/100).
pub camber: R,
/// Положение максимума кривизны в долях хорды (вторая цифра NACA/10).
pub camber_pos: R,
/// Контур в системе тела (без поворота), для многоугольных форм. /// Контур в системе тела (без поворота), для многоугольных форм.
poly: Vec<[R; 2]>, poly: Vec<[R; 2]>,
} }
impl Body { impl Body {
pub fn new(kind: ShapeKind, cx: R, cy: R, d: R, angle_deg: R, ratio: R) -> Self { 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 angle = angle_deg * PI / 180.0;
let poly = build_polygon(kind, d, ratio); let poly = if kind == ShapeKind::Polygon {
Body { kind, cx, cy, d, angle, ratio, poly } 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). /// Тот же контур, пересчитанный на уровень с измельчением `r` (координаты и размер ×r).
pub fn refined(&self, r: R, ox: R, oy: R) -> Body { 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.kind,
(self.cx - ox) * r, (self.cx - ox) * r,
(self.cy - oy) * r, (self.cy - oy) * r,
self.d * r, self.d * r,
self.angle * 180.0 / PI, self.angle * 180.0 / PI,
self.ratio, 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<Body>,
}
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; let h = 0.5 * d;
match kind { match kind {
ShapeKind::Cylinder | ShapeKind::Ellipse => Vec::new(), // аналитические ShapeKind::Cylinder | ShapeKind::Ellipse => Vec::new(), // аналитические
ShapeKind::Polygon => Vec::new(), // вершины задаются снаружи
ShapeKind::Square => vec![[-h, -h], [h, -h], [h, h], [-h, h]], 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::Diamond => vec![[h, 0.0], [0.0, h], [-h, 0.0], [0.0, -h]],
ShapeKind::Plate => { ShapeKind::Plate => {
@@ -475,36 +639,75 @@ fn build_polygon(kind: ShapeKind, d: R, ratio: R) -> Vec<[R; 2]> {
}) })
.collect() .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ξ⁴), /// Профиль NACA четырёхзначной серии.
/// ξ = x/c. Хорда направлена по +x тела, начало отсчёта смещено так, чтобы центр вращения ///
/// (точка приложения угла атаки) был в c/4 — стандартная аэродинамическая четверть хорды. /// Толщина: y_t = 5t·c·(0.2969√ξ − 0.1260ξ − 0.3516ξ² + 0.2843ξ³ − 0.1036ξ⁴), ξ = x/c.
fn naca_symmetric(chord: R, t: R, n: usize) -> Vec<[R; 2]> { /// Средняя линия (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| { let yt = |xi: R| {
5.0 * t 5.0 * t
* chord * chord
* (0.2969 * xi.sqrt() - 0.1260 * xi - 0.3516 * xi * xi + 0.2843 * xi.powi(3) * (0.2969 * xi.sqrt() - 0.1260 * xi - 0.3516 * xi * xi + 0.2843 * xi.powi(3)
- 0.1036 * xi.powi(4)) - 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<R> = (0..=n) let xs: Vec<R> = (0..=n).map(|k| 0.5 * (1.0 - (PI * k as R / n as R).cos())).collect();
.map(|k| 0.5 * (1.0 - (PI * k as R / n as R).cos())) let x0 = 0.25 * chord;
.collect();
let x0 = 0.25 * chord; // четверть хорды в нуле системы тела
let mut poly = Vec::with_capacity(2 * n); 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() { for &xi in xs.iter() {
// нижняя поверхность, от носка к хвосту (обход по часовой → замкнётся против) poly.push(point(xi, false));
poly.push([xi * chord - x0, -yt(xi)]);
} }
for &xi in xs.iter().rev().skip(1) { for &xi in xs.iter().rev().skip(1) {
poly.push([xi * chord - x0, yt(xi)]); poly.push(point(xi, true));
}
// профиль строился носком в −x; развернём, чтобы носок смотрел против потока (+x → навстречу)
for p in poly.iter_mut() {
p[0] = -p[0];
} }
poly poly
} }
@@ -617,8 +820,10 @@ pub struct Link {
pub kind: LinkKind, pub kind: LinkKind,
/// Доля пересечения q = |x_f→стенка| / |x_f→x_solid| ∈ (0,1), из SDF, клип [0.02, 0.98]. /// Доля пересечения q = |x_f→стенка| / |x_f→x_solid| ∈ (0,1), из SDF, клип [0.02, 0.98].
pub q: R, pub q: R,
/// Линк принадлежит обтекаемому телу (а не иной твёрдой поверхности) — по нему считается сила. /// Индекс тела, которого касается линк, зажатый в [0, MAX_BODY_BUCKETS). Сила считается
pub body: bool, /// по каждому такому ведру отдельно; тела за пределом ведёрного набора сваливаются в
/// последнее, так что СУММА по вёдрам остаётся точной, а разбивка — только для первых.
pub body: u8,
} }
/// Доля пересечения по SDF: φ_f / (φ_f − φ_solid). Клип отсекает вырожденные линки, /// Доля пересечения по SDF: φ_f / (φ_f − φ_solid). Клип отсекает вырожденные линки,