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

Снять предел GPU на размер сетки и молчаливые подмены модели стенки

Browse Source 
Диспетчеризация была одномерной, а у GPU жёсткий предел 65535 рабочих групп
на измерение: при 64 узлах на группу это упирало сетку в 4.2 млн узлов, то
есть примерно 2048×2048. Всё, что крупнее, не считало медленно, а падало
ошибкой валидации — то есть вся группа сверхмелких сеток кампании (4096×2048
и 4096×4096) не запустилась бы вовсе. Теперь сетка рабочих групп двумерная, а
линейный индекс собирается в шейдере через lin()/wlin(); отображение
«группа → узлы» при этом остаётся ровно линейным, поэтому редукции и буфер
частичных сумм не потребовали изменений, кроме отсечения хвостовых групп.
Проверено до 4096×4096; паритет с CPU не сдвинулся ни в одной цифре
(Cd 2.4650 против 2.4649, как и до правки).

--wall staircase на GPU молча считался по Bouzidi: ветка выбиралась по
is_moment_based(), и ступенчатая модель попадала в ту же ветку, что
интерполированный отскок. Обнаружилось по тому, что две модели дали
побитово одинаковый Cd там, где обязаны были разойтись. Вместо булева
«третий порядок» в Dyn теперь режим стенки числом, и все четыре модели
живут на обоих бэкендах.

clap отвергал --body-angle -8 как неизвестный ключ, из-за чего не
запускался прогон на зеркальную симметрию профиля.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
NotBigGhost
2026-08-15 17:36:18 +03:00
co-authored by Claude Opus 5
parent 9120affab5
commit 26a8475396
2 changed files with 90 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/theory/2d_solver/src/gpu.rs
+88 -42
View File
@@ -68,8 +68,8 @@ struct Dyn {
collision: u32,
/// Слот истории, в который `k_stats2` кладёт итоги этого шага.
slot: u32,
/// 1 — продолжать ряд Эрмита до 3-го порядка (HRR), 0 — обрыв на Π (Град).
wall_hrr: u32,
/// Режим стенки: 0 — Bouzidi, 1 — Град, 2 — HRR, 3 — простой отскок.
wall_mode: u32,
_pad: [u32; 2],
}
@@ -146,7 +146,8 @@ struct Results {
const SHADER: &str = r#"
// ───── структуры (обязаны совпадать с gpu.rs) ─────
struct LevelParams { n:u32, nx:u32, ny:u32, nlinks:u32, bcx:f32, bcy:f32, probe_node:u32, flags:u32, nwall:u32, lp1:u32, lp2:u32, lp3: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 };
// wall_mode: 0 — Bouzidi, 1 — Град, 2 — HRR, 3 — простой отскок
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_mode: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, body:u32, p1:u32 };
@@ -309,9 +310,18 @@ fn zou_he_outlet(fin: array<f32,9>, rho_out: f32, uy_out: f32) -> array<f32,9> {
// ───── ядра ─────
// У GPU жёсткий предел на число рабочих групп в ОДНОМ измерении (65535). Сетка 4096×2048 —
// это 131072 группы по 64 узла, то есть вдвое больше предела, и запуск падает с ошибкой
// валидации, а не считает медленно. Поэтому диспетчеризация двумерная, а линейный индекс
// собирается обратно вручную. Сборка точная: gid.x = wid.x*64 + lid.x, поэтому
// lin() = (wid.y*nwg.x + wid.x)*64 + lid.x — то есть номер группы тоже остаётся линейным.
fn lin(gid: vec3<u32>, nwg: vec3<u32>) -> u32 { return gid.y * nwg.x * 64u + gid.x; }
fn wlin(wid: vec3<u32>, nwg: vec3<u32>) -> u32 { return wid.y * nwg.x + wid.x; }
@compute @workgroup_size(64)
fn k_collide(@builtin(global_invocation_id) gid: vec3<u32>) {
let nd = gid.x;
fn k_collide(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let nd = lin(gid, nwg);
if (nd >= P.n) { return; }
var fv: array<f32,9>;
load9(&fv, nd, P.n);
@@ -327,8 +337,9 @@ fn k_collide(@builtin(global_invocation_id) gid: vec3<u32>) {
}
@compute @workgroup_size(64)
fn k_stream(@builtin(global_invocation_id) gid: vec3<u32>) {
let nd = gid.x;
fn k_stream(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let nd = lin(gid, nwg);
if (nd >= P.n) { return; }
var cx = array<i32,9>(0, 1, 0, -1, 0, 1, -1, -1, 1);
var cy = array<i32,9>(0, 0, 1, 0, -1, 1, 1, -1, -1);
@@ -344,12 +355,18 @@ fn k_stream(@builtin(global_invocation_id) gid: vec3<u32>) {
}
@compute @workgroup_size(64)
fn k_bouzidi(@builtin(global_invocation_id) gid: vec3<u32>) {
let k = gid.x;
fn k_bouzidi(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let k = lin(gid, nwg);
if (k >= P.nlinks) { return; }
let L = links[k];
let fi = post[L.i*P.n + L.node];
var v = fi;
// ступенчатая модель игнорирует долю пересечения: стенка ровно посередине между узлами
if (D.wall_mode == 3u) {
f[L.ib*P.n + L.node] = fi;
return;
}
if (L.kind == 0u) { // q < 1/2, есть дальний жидкий сосед
v = 2.0*L.q*fi + (1.0 - 2.0*L.q)*post[L.i*P.n + L.far];
} else if (L.kind == 1u) { // q >= 1/2
@@ -360,8 +377,9 @@ fn k_bouzidi(@builtin(global_invocation_id) gid: vec3<u32>) {
}
@compute @workgroup_size(64)
fn k_walls(@builtin(global_invocation_id) gid: vec3<u32>) {
let x = gid.x;
fn k_walls(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let x = lin(gid, nwg);
if (x >= P.nx) { return; }
// зеркальное отражение: касательный импульс сохраняется, нормальный заворачивается
f[2u*P.n + x] = f[4u*P.n + x];
@@ -374,8 +392,9 @@ fn k_walls(@builtin(global_invocation_id) gid: vec3<u32>) {
}
@compute @workgroup_size(64)
fn k_channel(@builtin(global_invocation_id) gid: vec3<u32>) {
let y = gid.x;
fn k_channel(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let y = lin(gid, nwg);
if (y >= P.ny) { return; }
// вход: скоростной Zou-He по ВСЕМУ столбцу, включая угловые узлы
let a = y*P.nx;
@@ -460,8 +479,9 @@ var<workgroup> wp: array<Partial, 64>;
@compute @workgroup_size(64)
fn k_stats1(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(local_invocation_id) lid: vec3<u32>,
@builtin(workgroup_id) wid: vec3<u32>) {
let nd = gid.x;
@builtin(workgroup_id) wid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let nd = lin(gid, nwg);
let t = lid.x;
var p: Partial;
p.rho = 0.0; p.maxu = 0.0; p.gsum = 0.0; p.gmin = 1e30; p.gmax = -1e30;
@@ -509,7 +529,10 @@ fn k_stats1(@builtin(global_invocation_id) gid: vec3<u32>,
workgroupBarrier();
s = s >> 1u;
}
if (t == 0u && (P.flags & 1u) != 0u) { parts[wid.x] = wp[0]; }
// хвостовые группы двумерной сетки лежат за пределами поля и в parts не пишут:
// отображение группа→узлы осталось линейным, поэтому wl >= nparts ⇒ узлы >= n
let wl = wlin(wid, nwg);
if (t == 0u && wl < D.nparts && (P.flags & 1u) != 0u) { parts[wl] = wp[0]; }
}
// Размер группы обязан совпадать с длиной wp: при 256 потоках на 64 ячейки четыре потока
@@ -676,8 +699,9 @@ fn moment_wall9(rho: f32, ux: f32, uy: f32, dudx: f32, dudy: f32, dvdx: f32, dvd
}
@compute @workgroup_size(64)
fn k_moment_wall(@builtin(global_invocation_id) gid: vec3<u32>) {
let w = gid.x;
fn k_moment_wall(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let w = lin(gid, nwg);
if (w >= P.nwall) { return; }
let wn = wnodes[w];
let node = wn.x;
@@ -715,7 +739,7 @@ fn k_moment_wall(@builtin(global_invocation_id) gid: vec3<u32>) {
let g4 = grad_u_at_t(node);
var g = moment_wall9(rho, ux, uy, g4.x, g4.y, g4.z, g4.w,
beta[node % P.nx], D.wall_hrr);
beta[node % P.nx], select(0u, 1u, D.wall_mode == 2u));
for (var i = 0u; i < 9u; i = i + 1u) {
if (missing[i] == 1u) { f[i*P.n + node] = g[i]; }
}
@@ -724,8 +748,9 @@ fn k_moment_wall(@builtin(global_invocation_id) gid: vec3<u32>) {
// ───── AMR ─────
@compute @workgroup_size(64)
fn k_ghost(@builtin(global_invocation_id) gid: vec3<u32>) {
let k = gid.x;
fn k_ghost(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let k = lin(gid, nwg);
if (k >= A.nghost) { return; }
let g = ghosts[k];
let w00 = (1.0 - g.tx)*(1.0 - g.ty);
@@ -753,8 +778,9 @@ fn k_ghost(@builtin(global_invocation_id) gid: vec3<u32>) {
}
@compute @workgroup_size(64)
fn k_fill(@builtin(global_invocation_id) gid: vec3<u32>) {
let k = gid.x;
fn k_fill(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let k = lin(gid, nwg);
if (k >= A.nghost) { return; }
let g = ghosts[k];
// временная интерполяция рамки между состояниями L0 «до» и «после» шага
@@ -764,8 +790,9 @@ fn k_fill(@builtin(global_invocation_id) gid: vec3<u32>) {
}
@compute @workgroup_size(64)
fn k_restrict(@builtin(global_invocation_id) gid: vec3<u32>) {
let k = gid.x;
fn k_restrict(@builtin(global_invocation_id) gid: vec3<u32>,
@builtin(num_workgroups) nwg: vec3<u32>) {
let k = lin(gid, nwg);
if (k >= A.nrestrict) { return; }
let pr = rest[k];
var ff: array<f32,9>;
@@ -1377,7 +1404,12 @@ impl Sim {
refine: refine as u32,
collision: (self.spec.collision == Collision::Bgk) as u32,
slot,
wall_hrr: (self.spec.wall == math::WallModel::Hrr) as u32,
wall_mode: match self.spec.wall {
math::WallModel::Bouzidi => 0,
math::WallModel::Grad => 1,
math::WallModel::Hrr => 2,
math::WallModel::Staircase => 3,
},
_pad: [0; 2],
}),
);
@@ -1400,25 +1432,25 @@ impl Sim {
p.set_bind_group(1, &self.dyn_binds[0], &[]);
let ncell = ceil_div(self.l0.n as u32, WG);
p.set_pipeline(&self.pipes.collide);
p.dispatch_workgroups(ncell, 1, 1);
dispatch(&mut p, ncell);
p.set_pipeline(&self.pipes.stream);
p.dispatch_workgroups(ncell, 1, 1);
dispatch(&mut p, ncell);
// Эталонные течения статей периодичны по обеим осям: ГУ не накладываются вовсе.
if channel {
if moment_wall {
p.set_bind_group(2, &self.pipes.empty_bg, &[]);
p.set_bind_group(3, &self.l0.wall_bind, &[]);
p.set_pipeline(&self.pipes.moment_wall);
p.dispatch_workgroups(ceil_div(self.l0.nwall.max(1), WG), 1, 1);
dispatch(&mut p, ceil_div(self.l0.nwall.max(1), WG));
} else {
p.set_pipeline(&self.pipes.bouzidi);
p.dispatch_workgroups(self.l0.link_groups(), 1, 1);
dispatch(&mut p, self.l0.link_groups());
}
// порядок обязателен: стенки снимают заворот по y, затем Zou-He — по x
p.set_pipeline(&self.pipes.walls);
p.dispatch_workgroups(ceil_div(self.l0.nx as u32, WG), 1, 1);
dispatch(&mut p, ceil_div(self.l0.nx as u32, WG));
p.set_pipeline(&self.pipes.channel);
p.dispatch_workgroups(ceil_div(self.l0.ny as u32, WG), 1, 1);
dispatch(&mut p, ceil_div(self.l0.ny as u32, WG));
}
}
@@ -1432,9 +1464,9 @@ impl Sim {
p.set_bind_group(1, &self.pipes.empty_bg, &[]);
p.set_pipeline(&self.pipes.ghost);
p.set_bind_group(2, &a.bg_ghost_old, &[]);
p.dispatch_workgroups(ceil_div(a.nghost, WG), 1, 1);
dispatch(&mut p, ceil_div(a.nghost, WG));
p.set_bind_group(2, &a.bg_ghost_new, &[]);
p.dispatch_workgroups(ceil_div(a.nghost, WG), 1, 1);
dispatch(&mut p, ceil_div(a.nghost, WG));
}
let ncell1 = ceil_div(l1.n as u32, WG);
let nlink1 = l1.link_groups();
@@ -1447,18 +1479,18 @@ impl Sim {
p.set_bind_group(0, &l1.bind, &[]);
p.set_bind_group(1, &self.dyn_binds[s], &[]);
p.set_pipeline(&self.pipes.collide);
p.dispatch_workgroups(ncell1, 1, 1);
dispatch(&mut p, ncell1);
p.set_pipeline(&self.pipes.stream);
p.dispatch_workgroups(ncell1, 1, 1);
dispatch(&mut p, ncell1);
if channel {
if moment_wall {
p.set_bind_group(2, &self.pipes.empty_bg, &[]);
p.set_bind_group(3, &l1.wall_bind, &[]);
p.set_pipeline(&self.pipes.moment_wall);
p.dispatch_workgroups(ceil_div(l1.nwall.max(1), WG), 1, 1);
dispatch(&mut p, ceil_div(l1.nwall.max(1), WG));
} else {
p.set_pipeline(&self.pipes.bouzidi);
p.dispatch_workgroups(nlink1, 1, 1);
dispatch(&mut p, nlink1);
}
}
// силу снимаем на каждом подшаге, усредняется она в k_stats2
@@ -1474,7 +1506,7 @@ impl Sim {
p.set_bind_group(1, &self.pipes.empty_bg, &[]);
p.set_bind_group(2, &a.bg_fill[s], &[]);
p.set_pipeline(&self.pipes.fill);
p.dispatch_workgroups(ceil_div(a.nghost, WG), 1, 1);
dispatch(&mut p, ceil_div(a.nghost, WG));
}
}
{
@@ -1486,7 +1518,7 @@ impl Sim {
p.set_bind_group(1, &self.pipes.empty_bg, &[]);
p.set_bind_group(2, &a.bg_restrict, &[]);
p.set_pipeline(&self.pipes.restrict);
p.dispatch_workgroups(ceil_div(a.nrestrict, WG), 1, 1);
dispatch(&mut p, ceil_div(a.nrestrict, WG));
}
} else {
let mut p = enc.begin_compute_pass(&wgpu::ComputePassDescriptor {
@@ -1507,11 +1539,11 @@ impl Sim {
p.set_bind_group(1, &self.dyn_binds[0], &[]);
p.set_bind_group(0, &self.l0.bind, &[]);
p.set_pipeline(&self.pipes.stats1);
p.dispatch_workgroups(self.l0.parts_count, 1, 1);
dispatch(&mut p, self.l0.parts_count);
if let Some(l1) = self.l1.as_ref() {
// на тонком уровне stats1 нужен только чтобы снять зонд (флаг записи сумм снят)
p.set_bind_group(0, &l1.bind, &[]);
p.dispatch_workgroups(ceil_div(l1.n as u32, WG), 1, 1);
dispatch(&mut p, ceil_div(l1.n as u32, WG));
p.set_bind_group(0, &self.l0.bind, &[]);
}
p.set_pipeline(&self.pipes.stats2);
@@ -1712,6 +1744,20 @@ impl Sim {
// Мелкие помощники
// ─────────────────────────────────────────────────────────────────────────────
/// Предел числа рабочих групп на измерение — 65535 (`maxComputeWorkgroupsPerDimension`).
/// Всё, что не влезло, раскладывается по второму измерению; шейдеры собирают линейный
/// индекс через lin()/wlin(). Хвост за пределами поля отсекается обычной проверкой границ.
fn wg_grid(groups: u32) -> (u32, u32) {
const MAX: u32 = 65535;
if groups <= MAX { (groups.max(1), 1) } else { (MAX, ceil_div(groups, MAX)) }
}
/// Запустить ядро на `groups` рабочих групп, разложив их по двум измерениям.
fn dispatch(p: &mut wgpu::ComputePass<'_>, groups: u32) {
let (gx, gy) = wg_grid(groups);
p.dispatch_workgroups(gx, gy, 1);
}
fn ceil_div(a: u32, b: u32) -> u32 {
(a + b - 1) / b
}
docs/theory/2d_solver/src/main.rs
+2 -2
View File
@@ -172,7 +172,7 @@ struct Cli {
#[arg(long, default_value_t = 150.0, help_heading = "Физика")]
re: R,
/// Направление потока, градусы (0 — вдоль канала)
#[arg(long, default_value_t = 0.0, help_heading = "Физика")]
#[arg(long, default_value_t = 0.0, help_heading = "Физика", allow_hyphen_values = true)]
flow_angle: R,
/// Постановка. channel — обтекание тела; остальные три периодичны по обеим осям, тела и
/// граничных условий не имеют и служат эталонами из статей: taylor-green (разд. VI, есть
@@ -210,7 +210,7 @@ struct Cli {
#[arg(long, default_value_t = 0.3, help_heading = "Тело")]
thickness: R,
/// Угол атаки тела, градусы
#[arg(long, default_value_t = 0.0, help_heading = "Тело")]
#[arg(long, default_value_t = 0.0, help_heading = "Тело", allow_hyphen_values = true)]
body_angle: R,
/// Положение центра тела по x, ячеек (по умолчанию nx/4)
#[arg(long, help_heading = "Тело")]