diff --git a/docs/theory/2d_solver/src/gpu.rs b/docs/theory/2d_solver/src/gpu.rs index fa54b11..12c0682 100644 --- a/docs/theory/2d_solver/src/gpu.rs +++ b/docs/theory/2d_solver/src/gpu.rs @@ -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, rho_out: f32, uy_out: f32) -> array { // ───── ядра ───── +// У 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, nwg: vec3) -> u32 { return gid.y * nwg.x * 64u + gid.x; } +fn wlin(wid: vec3, nwg: vec3) -> u32 { return wid.y * nwg.x + wid.x; } + @compute @workgroup_size(64) -fn k_collide(@builtin(global_invocation_id) gid: vec3) { - let nd = gid.x; +fn k_collide(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + let nd = lin(gid, nwg); if (nd >= P.n) { return; } var fv: array; load9(&fv, nd, P.n); @@ -327,8 +337,9 @@ fn k_collide(@builtin(global_invocation_id) gid: vec3) { } @compute @workgroup_size(64) -fn k_stream(@builtin(global_invocation_id) gid: vec3) { - let nd = gid.x; +fn k_stream(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + let nd = lin(gid, nwg); if (nd >= P.n) { return; } var cx = array(0, 1, 0, -1, 0, 1, -1, -1, 1); var cy = array(0, 0, 1, 0, -1, 1, 1, -1, -1); @@ -344,12 +355,18 @@ fn k_stream(@builtin(global_invocation_id) gid: vec3) { } @compute @workgroup_size(64) -fn k_bouzidi(@builtin(global_invocation_id) gid: vec3) { - let k = gid.x; +fn k_bouzidi(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + 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) { } @compute @workgroup_size(64) -fn k_walls(@builtin(global_invocation_id) gid: vec3) { - let x = gid.x; +fn k_walls(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + 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) { } @compute @workgroup_size(64) -fn k_channel(@builtin(global_invocation_id) gid: vec3) { - let y = gid.x; +fn k_channel(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + let y = lin(gid, nwg); if (y >= P.ny) { return; } // вход: скоростной Zou-He по ВСЕМУ столбцу, включая угловые узлы let a = y*P.nx; @@ -460,8 +479,9 @@ var wp: array; @compute @workgroup_size(64) fn k_stats1(@builtin(global_invocation_id) gid: vec3, @builtin(local_invocation_id) lid: vec3, - @builtin(workgroup_id) wid: vec3) { - let nd = gid.x; + @builtin(workgroup_id) wid: vec3, + @builtin(num_workgroups) nwg: vec3) { + 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, 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) { - let w = gid.x; +fn k_moment_wall(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + 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) { 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) { // ───── AMR ───── @compute @workgroup_size(64) -fn k_ghost(@builtin(global_invocation_id) gid: vec3) { - let k = gid.x; +fn k_ghost(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + 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) { } @compute @workgroup_size(64) -fn k_fill(@builtin(global_invocation_id) gid: vec3) { - let k = gid.x; +fn k_fill(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + 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) { } @compute @workgroup_size(64) -fn k_restrict(@builtin(global_invocation_id) gid: vec3) { - let k = gid.x; +fn k_restrict(@builtin(global_invocation_id) gid: vec3, + @builtin(num_workgroups) nwg: vec3) { + let k = lin(gid, nwg); if (k >= A.nrestrict) { return; } let pr = rest[k]; var ff: array; @@ -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 } diff --git a/docs/theory/2d_solver/src/main.rs b/docs/theory/2d_solver/src/main.rs index 62e1c1c..82d6e62 100644 --- a/docs/theory/2d_solver/src/main.rs +++ b/docs/theory/2d_solver/src/main.rs @@ -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 = "Тело")]