diff --git a/docs/theory/2d_solver/src/cpu.rs b/docs/theory/2d_solver/src/cpu.rs index c10b45b..da3b7f0 100644 --- a/docs/theory/2d_solver/src/cpu.rs +++ b/docs/theory/2d_solver/src/cpu.rs @@ -275,7 +275,9 @@ impl Level { } } - /// ГРАНИЧНОЕ УСЛОВИЕ ГРАДА на теле (Dorschner и др., JFM 801 (2016), прил. B). + /// ГРАНИЧНОЕ УСЛОВИЕ НА МОМЕНТАХ (Dorschner и др., JFM 801 (2016), прил. B). + /// + /// `third_order` — продолжать ли ряд Эрмита рекурсивно до третьего порядка (HRR). /// /// В отличие от Bouzidi, работающего полинково, здесь условие ставится сразу на весь /// граничный узел и не на популяции, а на моменты: @@ -286,7 +288,7 @@ impl Level { /// /// после чего недостающие популяции собираются приближением Града (2.13). Скорости /// соседей и градиенты берутся с прошлого шага — см. `u_prev`. - fn apply_grad_wall(&mut self) { + fn apply_moment_wall(&mut self, third_order: bool) { let nx = self.nx; // стенка неподвижна; для подвижного тела сюда пойдёт её скорость на линке, // и добавится динамическая часть плотности (B 4) @@ -325,7 +327,17 @@ impl Level { } let (dudx, dudy, dvdx, dvdy) = self.grad_u_at_t(node); - let g = math::grad_wall(rho, ux, uy, dudx, dudy, dvdx, dvdy, self.beta[node % nx]); + let g = math::moment_wall( + rho, + ux, + uy, + dudx, + dudy, + dvdx, + dvdy, + self.beta[node % nx], + third_order, + ); for i in 0..Q { if missing[i] { self.f[node][i] = g[i]; @@ -373,7 +385,8 @@ impl Level { fn apply_wall(&mut self, model: WallModel) { match model { WallModel::Bouzidi => self.apply_bouzidi(), - WallModel::Grad => self.apply_grad_wall(), + WallModel::Grad => self.apply_moment_wall(false), + WallModel::Hrr => self.apply_moment_wall(true), WallModel::Staircase => self.apply_staircase(), } } diff --git a/docs/theory/2d_solver/src/gpu.rs b/docs/theory/2d_solver/src/gpu.rs index 7063dc1..42824d9 100644 --- a/docs/theory/2d_solver/src/gpu.rs +++ b/docs/theory/2d_solver/src/gpu.rs @@ -25,6 +25,14 @@ use crate::{Collision, FieldKind, Spec, StepRec}; const WG: u32 = 64; +/// Сколько шагов копится в буфере итогов до одной синхронизации с устройством. +/// +/// Раньше после каждого шага делалось `map_async` + `poll(Wait)` ради 48 байт: на сетке +/// 240×120 счёт упирался в 863 шаг/с при том, что сам счёт занимал 0.27 мс из 1.16 — три +/// четверти времени машина стояла. Теперь `k_stats2` пишет итоги в слот `step % HIST`, а +/// хост читает всю пачку разом. Обязано совпадать с `HIST` в шейдере. +const HIST: usize = 128; + // ───────────────────────────────────────────────────────────────────────────── // Структуры, разделяемые с шейдером // ───────────────────────────────────────────────────────────────────────────── @@ -41,6 +49,9 @@ struct LevelParams { probe_node: u32, /// бит 0 — писать частичные суммы статистики; бит 1 — на этом уровне лежит зонд flags: u32, + /// число граничных узлов (для моментных моделей стенки) + nwall: u32, + _pad: [u32; 3], } #[repr(C)] @@ -55,6 +66,11 @@ struct Dyn { nparts: u32, refine: u32, collision: u32, + /// Слот истории, в который `k_stats2` кладёт итоги этого шага. + slot: u32, + /// 1 — продолжать ряд Эрмита до 3-го порядка (HRR), 0 — обрыв на Π (Град). + wall_hrr: u32, + _pad: [u32; 2], } #[repr(C)] @@ -126,8 +142,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 }; -struct Dyn { ux_in:f32, uy_in:f32, rho_out:f32, kbc_model:u32, outlet_extrap:u32, nparts:u32, refine:u32, collision:u32 }; +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 }; struct Substep { idx:u32, p0:u32, p1:u32, p2:u32 }; struct AmrParams { nghost:u32, nrestrict:u32, ccount:u32, fcount:u32, r01:f32, rfc:f32, w:f32, pad:f32 }; struct GLink { node:u32, far:u32, i:u32, ib:u32, kind:u32, q:f32, p0:u32, p1:u32 }; @@ -158,9 +174,28 @@ struct Partial { rho:f32, maxu:f32, gsum:f32, gmin:f32, gmax:f32, cnt:f32, degen @group(2) @binding(5) var rest : array>; @group(2) @binding(6) var A : AmrParams; +// (узел, смещение первого линка, число линков, —). Линки одного узла лежат подряд. +@group(3) @binding(0) var wnodes : array>; + const GREL: f32 = 1e-8; const CS2 : f32 = 0.3333333333; +// Раскладка буфера итогов. results[slot*12 .. +12) — история одного шага; за историей идут +// рабочие слоты силы по подшагам. История нужна, чтобы не синхронизироваться с устройством +// на каждом шаге: результаты копятся и читаются пачкой (см. HIST в gpu.rs). +const HIST: u32 = 128u; +const FORCE_BASE: u32 = 1536u; // = HIST*12u + +// Компенсированное сложение Кэхена–Ноймайера: возвращает (сумма, накопленная поправка). +// Наивная сумма по 10^5…10^7 значений в f32 съедает ~log2(N) бит; здесь потеря не копится, +// а итог берётся как s + c. Стоит несколько операций на элемент. +fn kadd(s: f32, c: f32, x: f32) -> vec2 { + let t = s + x; + var cc: f32; + if (abs(s) >= abs(x)) { cc = c + ((s - t) + x); } else { cc = c + ((x - t) + s); } + return vec2(t, cc); +} + // ───── физика: построчный перенос math.rs ───── // Энтропийное равновесие в product-form (точный максимизатор энтропии при заданных rho, rho*u) @@ -369,7 +404,9 @@ fn k_force(@builtin(local_invocation_id) lid: vec3) { let t = lid.x; var cx = array(0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0); var cy = array(0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0); - var sx = 0.0; var sy = 0.0; var sz = 0.0; + var sx = vec2(0.0, 0.0); + var sy = vec2(0.0, 0.0); + var sz = vec2(0.0, 0.0); var k = t; loop { if (k >= P.nlinks) { break; } @@ -378,15 +415,15 @@ fn k_force(@builtin(local_invocation_id) lid: vec3) { let fb = f[L.ib*P.n + L.node]; let dfx = cx[L.i]*fp + cx[L.i]*fb; let dfy = cy[L.i]*fp + cy[L.i]*fb; - sx = sx + dfx; - sy = sy + dfy; + sx = kadd(sx.x, sx.y, dfx); + sy = kadd(sy.x, sy.y, dfy); // плечо до точки пересечения линка со стенкой, а не до узла let rx = f32(L.node % P.nx) - P.bcx + L.q*cx[L.i]; let ry = f32(L.node / P.nx) - P.bcy + L.q*cy[L.i]; - sz = sz + rx*dfy - ry*dfx; + sz = kadd(sz.x, sz.y, rx*dfy - ry*dfx); k = k + 256u; } - wfx[t] = sx; wfy[t] = sy; wtz[t] = sz; + wfx[t] = sx.x + sx.y; wfy[t] = sy.x + sy.y; wtz[t] = sz.x + sz.y; workgroupBarrier(); var s = 128u; loop { @@ -400,9 +437,9 @@ fn k_force(@builtin(local_invocation_id) lid: vec3) { s = s >> 1u; } if (t == 0u) { - results[16u + S.idx*4u + 0u] = wfx[0]; - results[16u + S.idx*4u + 1u] = wfy[0]; - results[16u + S.idx*4u + 2u] = wtz[0]; + 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]; } } @@ -440,7 +477,7 @@ fn k_stats1(@builtin(global_invocation_id) gid: vec3, var fv: array; load9(&fv, nd, P.n); let m = macros9(fv); - results[3] = m.z; + results[D.slot*12u + 3u] = m.z; } wp[t] = p; workgroupBarrier(); @@ -469,60 +506,198 @@ fn k_stats1(@builtin(global_invocation_id) gid: vec3, @compute @workgroup_size(64) fn k_stats2(@builtin(local_invocation_id) lid: vec3) { 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; - p.cnt = 0.0; p.degen = 0.0; p.xineg = 0.0; + // слагаемых до 10^5 на поток — здесь и живёт основная потеря f32, поэтому все + // аддитивные накопители компенсированные + var rho = vec2(0.0, 0.0); + var gsum = vec2(0.0, 0.0); + var cnt = vec2(0.0, 0.0); + var degen = vec2(0.0, 0.0); + var xineg = vec2(0.0, 0.0); + var maxu = 0.0; var gmin = 1e30; var gmax = -1e30; var k = t; loop { if (k >= D.nparts) { break; } let o = parts[k]; - p.rho = p.rho + o.rho; - p.maxu = max(p.maxu, o.maxu); - p.gsum = p.gsum + o.gsum; - p.gmin = min(p.gmin, o.gmin); - p.gmax = max(p.gmax, o.gmax); - p.cnt = p.cnt + o.cnt; - p.degen = p.degen + o.degen; - p.xineg = p.xineg + o.xineg; + rho = kadd(rho.x, rho.y, o.rho); + gsum = kadd(gsum.x, gsum.y, o.gsum); + cnt = kadd(cnt.x, cnt.y, o.cnt); + degen = kadd(degen.x, degen.y, o.degen); + xineg = kadd(xineg.x, xineg.y, o.xineg); + maxu = max(maxu, o.maxu); + gmin = min(gmin, o.gmin); + gmax = max(gmax, o.gmax); k = k + 64u; } + var p: Partial; + p.rho = rho.x + rho.y; p.gsum = gsum.x + gsum.y; p.cnt = cnt.x + cnt.y; + p.degen = degen.x + degen.y; p.xineg = xineg.x + xineg.y; + p.maxu = maxu; p.gmin = gmin; p.gmax = gmax; wp[t] = p; workgroupBarrier(); // сведение 256 потоков через 64 ячейки делаем последовательно нулевым потоком: // объём данных крошечный, а корректность важнее пары микросекунд if (t == 0u) { - var q: Partial; - q.rho = 0.0; q.maxu = 0.0; q.gsum = 0.0; q.gmin = 1e30; q.gmax = -1e30; - q.cnt = 0.0; q.degen = 0.0; q.xineg = 0.0; + var qrho = vec2(0.0, 0.0); + var qgsum = vec2(0.0, 0.0); + var qcnt = vec2(0.0, 0.0); + var qdeg = vec2(0.0, 0.0); + var qxin = vec2(0.0, 0.0); + var qmaxu = 0.0; var qgmin = 1e30; var qgmax = -1e30; for (var j = 0u; j < 64u; j = j + 1u) { let o = wp[j]; - q.rho = q.rho + o.rho; - q.maxu = max(q.maxu, o.maxu); - q.gsum = q.gsum + o.gsum; - q.gmin = min(q.gmin, o.gmin); - q.gmax = max(q.gmax, o.gmax); - q.cnt = q.cnt + o.cnt; - q.degen = q.degen + o.degen; - q.xineg = q.xineg + o.xineg; + qrho = kadd(qrho.x, qrho.y, o.rho); + qgsum = kadd(qgsum.x, qgsum.y, o.gsum); + qcnt = kadd(qcnt.x, qcnt.y, o.cnt); + qdeg = kadd(qdeg.x, qdeg.y, o.degen); + qxin = kadd(qxin.x, qxin.y, o.xineg); + qmaxu = max(qmaxu, o.maxu); + qgmin = min(qgmin, o.gmin); + qgmax = max(qgmax, o.gmax); } var fx = 0.0; var fy = 0.0; var tz = 0.0; for (var s = 0u; s < D.refine; s = s + 1u) { - fx = fx + results[16u + s*4u + 0u]; - fy = fy + results[16u + s*4u + 1u]; - tz = tz + results[16u + s*4u + 2u]; + 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); - results[0] = fx*inv; - results[1] = fy*inv; - results[2] = tz*inv; - results[4] = q.rho; - results[5] = q.maxu; - results[6] = q.gsum; - results[7] = q.gmin; - results[8] = q.gmax; - results[9] = q.cnt; - results[10] = q.degen; - results[11] = q.xineg; + let o = D.slot*12u; + results[o + 0u] = fx*inv; + results[o + 1u] = fy*inv; + results[o + 2u] = tz*inv; + results[o + 4u] = qrho.x + qrho.y; + results[o + 5u] = qmaxu; + results[o + 6u] = qgsum.x + qgsum.y; + results[o + 7u] = qgmin; + results[o + 8u] = qgmax; + results[o + 9u] = qcnt.x + qcnt.y; + results[o + 10u] = qdeg.x + qdeg.y; + results[o + 11u] = qxin.x + qxin.y; + } +} + +// ───── граничное условие на моментах (Град / HRR) ───── + +// Скорость узла на момент t: берётся из ПОСТ-СТОЛКНОВИТЕЛЬНОГО поля. Столкновение сохраняет +// rho и rho*u точно, поэтому это ровно u(x, t) — то, что требует прил. B, без отдельного +// хранилища прошлого шага. +fn u_at_t(nd: u32) -> vec2 { + var ff: array; + for (var i = 0u; i < 9u; i = i + 1u) { ff[i] = post[i*P.n + nd]; } + let m = macros9(ff); + return vec2(m.y, m.z); +} + +// (du/dx, du/dy, dv/dx, dv/dy): центральная разность там, где оба соседа жидкие, +// односторонняя — где один твёрдый, ноль — если твёрдые оба. +fn grad_u_at_t(nd: u32) -> vec4 { + let x = nd % P.nx; + let y = nd / P.nx; + let xm = y*P.nx + (x + P.nx - 1u) % P.nx; + let xp = y*P.nx + (x + 1u) % P.nx; + let ym = ((y + P.ny - 1u) % P.ny)*P.nx + x; + let yp = ((y + 1u) % P.ny)*P.nx + x; + let uc = u_at_t(nd); + var dx = vec2(0.0, 0.0); + let am = solid[xm] == 0u; let ap = solid[xp] == 0u; + if (am && ap) { dx = 0.5*(u_at_t(xp) - u_at_t(xm)); } + else if (ap) { dx = u_at_t(xp) - uc; } + else if (am) { dx = uc - u_at_t(xm); } + var dy = vec2(0.0, 0.0); + let bm = solid[ym] == 0u; let bp = solid[yp] == 0u; + if (bm && bp) { dy = 0.5*(u_at_t(yp) - u_at_t(ym)); } + else if (bp) { dy = u_at_t(yp) - uc; } + else if (bm) { dy = uc - u_at_t(ym); } + return vec4(dx.x, dy.x, dx.y, dy.y); +} + +// Приближение Града, ур. (2.13): популяции по rho, u и полному тензору давлений. +fn grad_init9(rho: f32, ux: f32, uy: f32, pxx: f32, pxy: f32, pyy: f32) -> array { + let axx = pxx - rho*CS2; + let ayy = pyy - rho*CS2; + var w = array(4.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, + 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0); + var cx = array(0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0); + var cy = array(0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0); + var out: array; + for (var i = 0u; i < 9u; i = i + 1u) { + let quad = axx*(cx[i]*cx[i] - CS2) + 2.0*pxy*cx[i]*cy[i] + ayy*(cy[i]*cy[i] - CS2); + out[i] = w[i]*(rho + rho*(ux*cx[i] + uy*cy[i])/CS2 + quad/(2.0*CS2*CS2)); + } + return out; +} + +// third=1 — HRR: ряд Эрмита продолжается на 3-й порядок рекурсивно из неравновесных +// коэффициентов 2-го. Построчный перенос math::moment_wall. +fn moment_wall9(rho: f32, ux: f32, uy: f32, dudx: f32, dudy: f32, dvdx: f32, dvdy: f32, + b: f32, third: u32) -> array { + let pref = rho*CS2/(2.0*b); + let nxx = -pref*2.0*dudx; + let nyy = -pref*2.0*dvdy; + let nxy = -pref*(dudy + dvdx); + var out = grad_init9(rho, ux, uy, + rho*CS2 + rho*ux*ux + nxx, + rho*ux*uy + nxy, + rho*CS2 + rho*uy*uy + nyy); + if (third == 1u) { + let a3xxy = 2.0*ux*nxy + uy*nxx; + let a3xyy = 2.0*uy*nxy + ux*nyy; + let k = 1.0/(2.0*CS2*CS2*CS2); + var w = array(4.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, + 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0); + var cx = array(0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0); + var cy = array(0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0); + for (var i = 0u; i < 9u; i = i + 1u) { + out[i] = out[i] + w[i]*k*((cx[i]*cx[i] - CS2)*cy[i]*a3xxy + + cx[i]*(cy[i]*cy[i] - CS2)*a3xyy); + } + } + return out; +} + +@compute @workgroup_size(64) +fn k_moment_wall(@builtin(global_invocation_id) gid: vec3) { + let w = gid.x; + if (w >= P.nwall) { return; } + let wn = wnodes[w]; + let node = wn.x; + let first = wn.y; + let cnt = wn.z; + + // целевая скорость (B 1): интерполяция между стенкой и дальним жидким соседом + var ux = 0.0; + var uy = 0.0; + for (var k = 0u; k < cnt; k = k + 1u) { + let L = links[first + k]; + var fx = 0.0; + var fy = 0.0; // стенка неподвижна; дальнего соседа может не быть + if (solid[L.far] == 0u) { + let uf = u_at_t(L.far); + fx = uf.x; + fy = uf.y; + } + ux = ux + (L.q*fx) / (1.0 + L.q); + uy = uy + (L.q*fy) / (1.0 + L.q); + } + let inv = 1.0 / f32(cnt); + ux = ux * inv; + uy = uy * inv; + + // целевая плотность (B 3): известные популяции плюс отскок недостающих + var missing = array(0u,0u,0u,0u,0u,0u,0u,0u,0u); + for (var k = 0u; k < cnt; k = k + 1u) { missing[links[first + k].ib] = 1u; } + var opp = array(0u, 3u, 4u, 1u, 2u, 7u, 8u, 5u, 6u); + var rho = 0.0; + for (var i = 0u; i < 9u; i = i + 1u) { + if (missing[i] == 1u) { rho = rho + post[opp[i]*P.n + node]; } + else { rho = rho + f[i*P.n + node]; } + } + + 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); + for (var i = 0u; i < 9u; i = i + 1u) { + if (missing[i] == 1u) { f[i*P.n + node] = g[i]; } } } @@ -595,6 +770,9 @@ struct GpuLevel { /// Стабилизатор γ поузлово — нужен для картинки; забирается с устройства как есть. gam: wgpu::Buffer, bind: wgpu::BindGroup, + /// Группа 3: индекс граничных узлов для моментной стенки. + wall_bind: wgpu::BindGroup, + nwall: u32, /// число рабочих групп редукции статистики (оно же длина буфера частичных сумм) parts_count: u32, nlinks: u32, @@ -628,6 +806,7 @@ fn make_level( probe_node: u32, flags: u32, u0: (R, R), + wall_layout: &wgpu::BindGroupLayout, ) -> GpuLevel { let n = nx * ny; let init = soa_equilibrium(n, u0); @@ -668,6 +847,20 @@ fn make_level( let links_buf = storage_init(device, "links", bytemuck::cast_slice(&links)); let beta_buf = storage_init(device, "beta", bytemuck::cast_slice(beta)); + // индекс граничных узлов: (узел, смещение первого линка, число линков, —) + let wnodes: Vec<[u32; 4]> = geom + .wall_nodes + .iter() + .map(|w| [w.node, w.first, w.count as u32, 0]) + .collect(); + let nwall = wnodes.len() as u32; + let wall_buf = storage_init(device, "wall nodes", bytemuck::cast_slice(&wnodes)); + let wall_bind = device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("wall"), + layout: wall_layout, + entries: &[bind(0, &wall_buf)], + }); + let parts_count = ceil_div(n as u32, WG); let parts = storage(device, "parts", (parts_count as u64) * 32); let params = device.create_buffer_init(&wgpu::util::BufferInitDescriptor { @@ -681,6 +874,8 @@ fn make_level( bcy: geom.body_cy as f32, probe_node, flags, + nwall, + _pad: [0; 3], }), usage: wgpu::BufferUsages::UNIFORM, }); @@ -699,7 +894,7 @@ fn make_level( bind(7, ¶ms), ], }); - GpuLevel { nx, ny, n, f, gam, bind, parts_count, nlinks: links.len() as u32 } + GpuLevel { nx, ny, n, f, gam, bind, wall_bind, nwall, parts_count, nlinks: links.len() as u32 } } // ───────────────────────────────────────────────────────────────────────────── @@ -730,6 +925,9 @@ pub struct Sim { fluid_count: R, step_index: u64, d_ref: R, + /// Номера шагов, посчитанных на устройстве, но ещё не прочитанных на хост. + /// Длина = занятые слоты истории. + pending: Vec, } struct AmrRes { @@ -745,6 +943,7 @@ struct Pipes { collide: wgpu::ComputePipeline, stream: wgpu::ComputePipeline, bouzidi: wgpu::ComputePipeline, + moment_wall: wgpu::ComputePipeline, walls: wgpu::ComputePipeline, channel: wgpu::ComputePipeline, force: wgpu::ComputePipeline, @@ -784,11 +983,17 @@ impl Sim { )) .map_err(|e| format!("не удалось получить устройство: {e}"))?; - // Условие Града на GPU пока не перенесено. Молча считать по Bouzidi нельзя — это - // была бы другая физика под тем же ключом, поэтому отказываемся явно. - if spec.wall == math::WallModel::Grad { - return Err("--wall grad на GPU не реализован (перенесён только Bouzidi). Запусти с --backend cpu либо оставь --wall bouzidi" - .into()); + // Моментная стенка добавляет к связке уровня ещё один storage-биндинг (индекс + // граничных узлов), итого 9 против 8, гарантируемых лимитами по умолчанию. На + // десктопных адаптерах их обычно >= 16, но проверить надо явно. + if spec.wall.is_moment_based() { + let lim = adapter.limits().max_storage_buffers_per_shader_stage; + if lim < 9 { + return Err(format!( + "адаптер даёт только {lim} storage-биндингов на стадию, моментной стенке \ + нужно 9. Запусти с --wall bouzidi либо --backend cpu" + )); + } } // ── топология берётся из процессорного бэкенда, а не строится заново ── @@ -815,6 +1020,10 @@ impl Sim { let bgl_level = level_layout(&device); let bgl_dyn = dyn_layout(&device); let bgl_amr = amr_layout(&device); + let bgl_wall = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("wall"), + entries: &[ro(0)], + }); let bgl_empty = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { label: Some("empty"), entries: &[], @@ -825,6 +1034,11 @@ impl Sim { bind_group_layouts: &[&bgl_level, &bgl_dyn], push_constant_ranges: &[], }); + let pl_wall = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { + label: Some("moment wall"), + bind_group_layouts: &[&bgl_level, &bgl_dyn, &bgl_empty, &bgl_wall], + push_constant_ranges: &[], + }); let pl_amr = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { label: Some("amr"), bind_group_layouts: &[&bgl_empty, &bgl_empty, &bgl_amr], @@ -849,16 +1063,17 @@ impl Sim { usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST, mapped_at_creation: false, }); - // 16 слотов итогов + 8 подшагов по 4 числа на силу + // история на HIST шагов + рабочие слоты силы на 8 подшагов + let hist_bytes = (HIST * std::mem::size_of::()) as u64; let results_buf = device.create_buffer(&wgpu::BufferDescriptor { label: Some("results"), - size: (16 + 8 * 4) * 4, + size: hist_bytes + 8 * 4 * 4, usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, mapped_at_creation: false, }); let staging = device.create_buffer(&wgpu::BufferDescriptor { label: Some("staging"), - size: std::mem::size_of::() as u64, + size: hist_bytes, usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST, mapped_at_creation: false, }); @@ -906,6 +1121,7 @@ impl Sim { probe0 as u32, flags0, u0, + &bgl_wall, ); let pre = device.create_buffer(&wgpu::BufferDescriptor { @@ -934,7 +1150,8 @@ impl Sim { let probe1 = if probe_on_fine { ((py - ay) * r) * nfx + (px - ax) * r } else { 0 }; let flags1 = if probe_on_fine { 2 } else { 0 }; let lvl1 = - make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1, u0); + make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1, + u0, &bgl_wall); let ghosts: Vec = patch .ghosts() @@ -1021,6 +1238,7 @@ impl Sim { collide: mk(&pl_level, "k_collide"), stream: mk(&pl_level, "k_stream"), bouzidi: mk(&pl_level, "k_bouzidi"), + moment_wall: mk(&pl_wall, "k_moment_wall"), walls: mk(&pl_level, "k_walls"), channel: mk(&pl_level, "k_channel"), force: mk(&pl_level, "k_force"), @@ -1065,6 +1283,7 @@ impl Sim { readback, fluid_count, step_index: 0, + pending: Vec::with_capacity(HIST), d_ref, }) } @@ -1095,10 +1314,14 @@ impl Sim { (u * c, uy) } - pub fn step(&mut self) -> StepRec { + /// Посчитать один шаг. Результаты не читаются сразу: они копятся в истории на устройстве + /// и попадают в `out` пачкой — либо когда история заполнится, либо по явному `flush`. + pub fn advance(&mut self, out: &mut Vec) { let t = self.step_index; + let slot = self.pending.len() as u32; let (ux_in, uy_in) = self.inlet(t); let refine = self.spec.refine.max(1); + let moment_wall = self.spec.wall.is_moment_based(); self.queue.write_buffer( &self.dyn_buf, 0, @@ -1111,6 +1334,9 @@ impl Sim { nparts: self.l0.parts_count, refine: refine as u32, collision: (self.spec.collision == Collision::Bgk) as u32, + slot, + wall_hrr: (self.spec.wall == math::WallModel::Hrr) as u32, + _pad: [0; 2], }), ); @@ -1135,8 +1361,15 @@ impl Sim { p.dispatch_workgroups(ncell, 1, 1); p.set_pipeline(&self.pipes.stream); p.dispatch_workgroups(ncell, 1, 1); - p.set_pipeline(&self.pipes.bouzidi); - p.dispatch_workgroups(self.l0.link_groups(), 1, 1); + 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); + } else { + p.set_pipeline(&self.pipes.bouzidi); + p.dispatch_workgroups(self.l0.link_groups(), 1, 1); + } // порядок обязателен: стенки снимают заворот по y, затем Zou-He — по x, на весь столбец p.set_pipeline(&self.pipes.walls); p.dispatch_workgroups(ceil_div(self.l0.nx as u32, WG), 1, 1); @@ -1172,8 +1405,15 @@ impl Sim { p.dispatch_workgroups(ncell1, 1, 1); p.set_pipeline(&self.pipes.stream); p.dispatch_workgroups(ncell1, 1, 1); - p.set_pipeline(&self.pipes.bouzidi); - p.dispatch_workgroups(nlink1, 1, 1); + 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); + } else { + p.set_pipeline(&self.pipes.bouzidi); + p.dispatch_workgroups(nlink1, 1, 1); + } // силу снимаем на каждом подшаге, усредняется она в k_stats2 p.set_pipeline(&self.pipes.force); p.dispatch_workgroups(1, 1, 1); @@ -1231,18 +1471,50 @@ impl Sim { p.dispatch_workgroups(1, 1, 1); } - enc.copy_buffer_to_buffer( - &self.results_buf, - 0, - &self.staging, - 0, - std::mem::size_of::() as u64, - ); + self.queue.submit(Some(enc.finish())); + self.pending.push(t); + self.step_index += 1; + if self.pending.len() == HIST { + self.flush(out); + } + } + + /// Дочитать всё, что уже посчитано устройством. Обязательно вызывать перед тем, как + /// смотреть на поля (кадр гифки, проверка на NaN) — иначе записи отстанут от состояния. + pub fn flush(&mut self, out: &mut Vec) { + let n = self.pending.len(); + if n == 0 { + return; + } + let bytes = (n * std::mem::size_of::()) as u64; + let mut enc = self + .device + .create_command_encoder(&wgpu::CommandEncoderDescriptor { label: Some("flush") }); + enc.copy_buffer_to_buffer(&self.results_buf, 0, &self.staging, 0, bytes); self.queue.submit(Some(enc.finish())); - let res: Results = self.read_staging(); - self.step_index += 1; + let slice = self.staging.slice(..bytes); + let (tx, rx) = std::sync::mpsc::channel(); + slice.map_async(wgpu::MapMode::Read, move |r| { + let _ = tx.send(r); + }); + self.device.poll(wgpu::Maintain::Wait); + let got: Vec = match rx.recv() { + Ok(Ok(())) => { + let data = slice.get_mapped_range(); + bytemuck::cast_slice::(&data[..bytes as usize]).to_vec() + } + _ => vec![Results::default(); n], + }; + self.staging.unmap(); + for (i, &t) in self.pending.iter().enumerate() { + out.push(self.make_rec(t, &got[i])); + } + self.pending.clear(); + } + + fn make_rec(&self, t: u64, res: &Results) -> StepRec { let cnt = if res.cnt > 0.0 { res.cnt as R } else { self.fluid_count }; StepRec { step: t, @@ -1260,24 +1532,6 @@ impl Sim { } } - fn read_staging(&self) -> Results { - let slice = self.staging.slice(..); - let (tx, rx) = std::sync::mpsc::channel(); - slice.map_async(wgpu::MapMode::Read, move |r| { - let _ = tx.send(r); - }); - self.device.poll(wgpu::Maintain::Wait); - let out = match rx.recv() { - Ok(Ok(())) => { - let data = slice.get_mapped_range(); - *bytemuck::from_bytes::(&data[..std::mem::size_of::()]) - } - _ => Results::default(), - }; - self.staging.unmap(); - out - } - /// Скачать популяции L0 на хост (нужно для кадров и проверки на NaN). fn download_l0(&self) -> Vec { self.download(&self.l0.f, (9 * self.l0.n * 4) as u64) diff --git a/docs/theory/2d_solver/src/main.rs b/docs/theory/2d_solver/src/main.rs index d4072f3..d083d49 100644 --- a/docs/theory/2d_solver/src/main.rs +++ b/docs/theory/2d_solver/src/main.rs @@ -197,11 +197,13 @@ struct Cli { /// Оператор столкновения #[arg(long, default_value = "kbc", value_parser = ["kbc", "bgk"], help_heading = "Схема")] collision: String, - /// Модель стенки на теле: bouzidi — интерполированный отскок по линкам; grad — условие - /// Града (Dorschner и др., JFM 801 (2016)), где задаются не популяции, а целевые моменты - /// ρ, u и тензор давлений. Авторы метода предпочитают grad: интерполяционные схемы, по их - /// словам, «ограничены низкими Re, поскольку на границе возникают паразитные скачки». - #[arg(long, default_value = "bouzidi", value_parser = math::WallModel::ALL, + /// Модель стенки на теле. hrr (умолчание) — восстановление недостающих популяций по + /// целевым моментам с рекурсивной регуляризацией до 3-го порядка Эрмита; grad — то же, + /// но с обрывом на тензоре давлений; bouzidi — интерполированный отскок по каждому линку + /// со своей долей пересечения q; staircase — простой отскок, q игнорируется (не для + /// счёта, а как база сравнения). Замерено: по разрешению геометрии bouzidi точнее + /// моментных схем, зато у тех естественная форма для подвижных стенок — см. README. + #[arg(long, default_value = "hrr", value_parser = math::WallModel::ALL, help_heading = "Схема")] wall: String, /// Состав сдвиговой части KBC (табл. I 2D-статьи): n1 — только девиатор {N, Π_xy} @@ -451,11 +453,23 @@ enum Backend { } impl Backend { - fn step(&mut self) -> StepRec { + /// Посчитать шаг. Готовые записи ДОПИСЫВАЮТСЯ в `out` — их может быть ноль (GPU копит + /// итоги в истории и читает пачкой) или сразу много (когда история заполнилась). + fn advance(&mut self, out: &mut Vec) { match self { - Backend::Cpu(s) => s.step(), + Backend::Cpu(s) => out.push(s.step()), #[cfg(feature = "gpu")] - Backend::Gpu(s) => s.step(), + Backend::Gpu(s) => s.advance(out), + } + } + /// Дочитать всё посчитанное. Обязательно перед чтением полей: иначе записи отстанут + /// от состояния, которое покажет кадр. + fn flush(&mut self, out: &mut Vec) { + let _ = &out; // процессорный бэкенд отдаёт записи сразу, копить нечего + match self { + Backend::Cpu(_) => {} + #[cfg(feature = "gpu")] + Backend::Gpu(s) => s.flush(out), } } fn sample_field(&mut self, k: FieldKind) -> (Vec, Vec) { @@ -598,11 +612,18 @@ fn run(cli: Cli) -> Result<(), String> { .unwrap_or(0)) as f64; for t in 0..=spec.steps { - let rec = back.step(); - recs.push(rec); + back.advance(&mut recs); - if let Some(w) = writer.as_mut() { - if 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; + // и кадр, и живая строка требуют актуальных данных — синхронизируемся только здесь + if need_frame || need_report { + back.flush(&mut recs); + } + let rec = *recs.last().unwrap_or(&StepRec::default()); + + if need_frame { + if let Some(w) = writer.as_mut() { let (field, solid) = back.sample_field(field_kind); let (cd, _, _) = math::coefficients(rec.fx, rec.fy, rec.tz, spec.units.u_lat, d_ref); @@ -610,18 +631,25 @@ fn run(cli: Cli) -> Result<(), String> { .map_err(|e| format!("запись кадра: {e}"))?; } } - - if cli.report_every > 0 && !quiet && t % cli.report_every == 0 { + if need_report { live_line(&spec, &rec, d_ref, t_start.elapsed().as_secs_f64(), nodes_per_step, full); } - // развал счёта: проверяем редко, проверка стоит полного прохода по полю - if t % 500 == 0 && !back.is_finite() { + // Развал счёта ловится по УЖЕ ПОСЧИТАННЫМ величинам, а не полным проходом по полю: + // на сетке 4096×2048 такой проход тянет с устройства сотни мегабайт, и делать его + // регулярно нельзя. Полная проверка — один раз, для подтверждения. + if !rec.max_u.is_finite() || !rec.rho_mean.is_finite() { + back.flush(&mut recs); eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ на шаге {t}: в поле появились NaN/inf."); blew_up = true; break; } } + back.flush(&mut recs); + if !blew_up && !back.is_finite() { + eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ: итоговое поле содержит NaN/inf."); + blew_up = true; + } let wall = t_start.elapsed().as_secs_f64(); if let Some(w) = writer { @@ -729,6 +757,7 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) { match spec.wall { math::WallModel::Bouzidi => "Bouzidi", math::WallModel::Grad => "Град (моменты)", + math::WallModel::Hrr => "HRR (моменты, 3-й пор.)", math::WallModel::Staircase => "простой отскок", }); println!(" выход по u_y {:>12} губка {} столбцов ×{:.0}", diff --git a/docs/theory/2d_solver/src/math.rs b/docs/theory/2d_solver/src/math.rs index 35718bf..e1be810 100644 --- a/docs/theory/2d_solver/src/math.rs +++ b/docs/theory/2d_solver/src/math.rs @@ -672,8 +672,8 @@ pub enum WallModel { /// после чего недостающие популяции собираются приближением Града (ур. 2.13). Авторы /// метода предпочитают его интерполяционным схемам: те «ограничены низкими числами /// Рейнольдса, поскольку на границе возникают паразитные скачки» (разд. 2.1). - /// Заодно это естественная форма для подвижных стенок: скорость стенки входит в целевые - /// значения, а не в отдельную поправку. + /// Обрыв ряда Эрмита на 2-м порядке. Заодно естественная форма для подвижных стенок: + /// скорость стенки входит в целевые значения, а не в отдельную поправку. /// /// ВАЖНО про субсеточность: положение стенки входит сюда ТОЛЬКО через целевую скорость /// (B 1) — одну усреднённую по узлу величину. Целевая плотность (B 3) и тензор давлений @@ -683,6 +683,17 @@ pub enum WallModel { /// видно и по чувствительности к положению тела внутри клетки, и по точности на грубых /// сетках (см. README). Grad, + /// HRR — рекурсивная регуляризация высокого порядка (Malaspinas 2015; Coreixas и др., + /// PRE 96, 033306). То же восстановление по целевым моментам, что у Града, но ряд Эрмита + /// продолжен на третий порядок, причём коэффициенты берутся не из популяций (их на стенке + /// как раз и не хватает), а РЕКУРСИВНО из уже известных второго порядка: + /// + /// a₃_xxy = 2·u_x·a₂_xy + u_y·a₂_xx, a₃_xyy = 2·u_y·a₂_xy + u_x·a₂_yy + /// + /// В D2Q9 компоненты a₃_xxx и a₃_yyy решёткой не поддерживаются и отбрасываются. + /// Третий порядок не трогает ρ, ρu и Π — соответствующие моменты весов обнуляются по + /// симметрии, — поэтому целевые значения выполняются ровно так же, как у Града. + Hrr, /// ПРОСТОЙ ОТСКОК: q игнорируется, стенка по построению лежит ровно посередине между /// узлами — ступенчатая поверхность. Первый порядок. Держится не для счёта, а как база /// сравнения: показывает, сколько именно даёт субсеточность. @@ -694,11 +705,17 @@ impl WallModel { Some(match s.to_ascii_lowercase().as_str() { "bouzidi" | "bb" => WallModel::Bouzidi, "grad" => WallModel::Grad, + "hrr" => WallModel::Hrr, "staircase" | "step" => WallModel::Staircase, _ => return None, }) } - pub const ALL: [&'static str; 3] = ["bouzidi", "grad", "staircase"]; + pub const ALL: [&'static str; 4] = ["hrr", "bouzidi", "grad", "staircase"]; + + /// Ставится ли условие на моменты (в отличие от полинковых схем). + pub fn is_moment_based(&self) -> bool { + matches!(self, WallModel::Grad | WallModel::Hrr) + } } /// Недостающие популяции по целевым моментам — ур. (2.13) вместе с (2.14)–(2.16): @@ -710,7 +727,7 @@ impl WallModel { /// граничном узле ещё не определена — недостающие популяции как раз и вычисляются. #[allow(clippy::too_many_arguments)] #[inline] -pub fn grad_wall( +pub fn moment_wall( rho: R, ux: R, uy: R, @@ -719,12 +736,37 @@ pub fn grad_wall( dvdx: R, dvdy: R, beta: R, + third_order: bool, ) -> [R; Q] { + // неравновесная часть тензора давлений — ур. (2.16) let pref = rho * CS2 / (2.0 * beta); - let pxx = rho * CS2 + rho * ux * ux - pref * 2.0 * dudx; - let pyy = rho * CS2 + rho * uy * uy - pref * 2.0 * dvdy; - let pxy = rho * ux * uy - pref * (dudy + dvdx); - grad_init(rho, ux, uy, pxx, pxy, pyy) + let nxx = -pref * 2.0 * dudx; + let nyy = -pref * 2.0 * dvdy; + let nxy = -pref * (dudy + dvdx); + + let mut out = grad_init( + rho, + ux, + uy, + rho * CS2 + rho * ux * ux + nxx, + rho * ux * uy + nxy, + rho * CS2 + rho * uy * uy + nyy, + ); + if !third_order { + return out; + } + // HRR: коэффициенты 3-го порядка строятся рекурсивно из НЕРАВНОВЕСНЫХ коэффициентов 2-го + let a3xxy = 2.0 * ux * nxy + uy * nxx; + let a3xyy = 2.0 * uy * nxy + ux * nyy; + // H³_xxy = (c_x² − c_s²)c_y, H³_xyy = c_x(c_y² − c_s²); множитель 1/(2c_s⁶) вместо + // 1/(6c_s⁶) — это три перестановки каждого индекса + let k = 1.0 / (2.0 * CS2 * CS2 * CS2); + for i in 0..Q { + let cx = CX[i] as R; + let cy = CY[i] as R; + out[i] += W[i] * k * ((cx * cx - CS2) * cy * a3xxy + cx * (cy * cy - CS2) * a3xyy); + } + out } /// Целевая скорость граничного узла — ур. (B 1): линейная интерполяция между стенкой и @@ -1121,13 +1163,53 @@ mod tests { assert!((mxy - pxy).abs() < 1e-13, "Π_xy: {mxy} vs {pxy}"); } + /// Третий порядок HRR не имеет права трогать ρ, ρu и Π — иначе целевые значения, + /// ради которых условие и ставится на моменты, перестанут выполняться. + #[test] + fn hrr_third_order_preserves_target_moments() { + let (rho, ux, uy) = (1.03, 0.05, -0.02); + let (dudx, dudy, dvdx, dvdy) = (0.003, -0.004, 0.002, -0.003); + let g = moment_wall(rho, ux, uy, dudx, dudy, dvdx, dvdy, 0.8, false); + let h = moment_wall(rho, ux, uy, dudx, dudy, dvdx, dvdy, 0.8, true); + // третий порядок реально что-то добавил + let diff: R = (0..Q).map(|i| (h[i] - g[i]).abs()).sum(); + assert!(diff > 1e-9, "HRR не отличается от Града: {diff:e}"); + // но моменты те же + for (a, b) in [(macros(&g), macros(&h))] { + assert!((a.0 - b.0).abs() < 1e-13, "ρ разошлась"); + assert!((a.1 - b.1).abs() < 1e-13, "u_x разошлась"); + assert!((a.2 - b.2).abs() < 1e-13, "u_y разошлась"); + } + for (cxp, cyp) in [(2u32, 0u32), (1, 1), (0, 2)] { + let m = |f: &[R; Q]| -> R { + (0..Q) + .map(|i| { + (CX[i] as R).powi(cxp as i32) * (CY[i] as R).powi(cyp as i32) * f[i] + }) + .sum() + }; + assert!((m(&g) - m(&h)).abs() < 1e-13, "Π_{cxp}{cyp} разошёлся"); + } + } + + /// При нулевых градиентах неравновесной части нет, значит и рекурсивные коэффициенты + /// третьего порядка нулевые — HRR обязан совпасть с Градом побитово. + #[test] + fn hrr_equals_grad_without_gradients() { + let g = moment_wall(1.0, 0.04, 0.02, 0.0, 0.0, 0.0, 0.0, 0.9, false); + let h = moment_wall(1.0, 0.04, 0.02, 0.0, 0.0, 0.0, 0.0, 0.9, true); + for i in 0..Q { + assert_eq!(g[i], h[i], "i={i}"); + } + } + /// Граничная сборка (2.13)+(2.14)–(2.16): при нулевых градиентах неравновесной части нет, /// и результат обязан совпасть с равновесием по тем же ρ, u с точностью до порядка Ma³ /// (product-form равновесие не полиномиально, поэтому не побитово). #[test] fn grad_wall_without_gradients_is_near_equilibrium() { let (rho, ux, uy) = (1.0, 0.02, 0.01); - let g = grad_wall(rho, ux, uy, 0.0, 0.0, 0.0, 0.0, 0.9); + let g = moment_wall(rho, ux, uy, 0.0, 0.0, 0.0, 0.0, 0.9, false); let fe = feq(rho, ux, uy); for i in 0..Q { assert!((g[i] - fe[i]).abs() < 1e-4, "i={i}: {} vs {}", g[i], fe[i]);