Батчинг GPU-чтения, компенсированные суммы и HRR-стенка на обоих бэкендах

БАТЧИНГ. Раньше после каждого шага делался map_async + poll(Wait) ради 48 байт итогов: на
сетке 240×120 счёт упирался в 863 шаг/с при том, что сам счёт занимал 0.27 мс из 1.16 — три
четверти времени машина стояла. Теперь k_stats2 пишет итоги в слот истории step % 128, а хост
читает пачкой. Контракт бэкенда изменён с пошагового step() на advance(&mut out) + flush():
записи дописываются в вектор, синхронизация происходит только там, где дальше нужны
актуальные данные (кадр гифки, живая строка отчёта).

Замер: 240×120 — 6715 шаг/с против 863 (7.8×); пропускная способность 195 MLUPS против прежних
105. Паритет CPU/GPU сохранён: Cd 1.798 против 1.799, ⟨ρ⟩ 1.00090 против 1.00091.

КОМПЕНСИРОВАННОЕ СУММИРОВАНИЕ (Кэхена–Ноймайера) в редукциях статистики и в сумме сил. Именно
там теряется основная точность f32: наивная сумма по 10^5…10^7 узлам съедает ~log2(N) бит.
Аппаратного f64 на целевом железе нет (в WGSL типа f64 не существует вовсе, а локальная Iris Xe
сообщает shaderFloat64 = false), поэтому компенсация — единственный доступный способ.

РАЗВАЛ СЧЁТА теперь ловится по уже посчитанным max|u| и ⟨ρ⟩, а не полным проходом по полю:
на сетке 4096×2048 такой проход тянет с устройства сотни мегабайт, и делать его регулярно
нельзя. Полная проверка осталась одна, в конце.

HRR-СТЕНКА (--wall hrr, теперь умолчание). Условие Града — это ряд Эрмита, оборванный на 2-м
порядке (ρ, u, Π). HRR продолжает его на третий, вычисляя коэффициенты не из популяций (их на
стенке как раз и не хватает), а рекурсивно из уже известных:
  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 и Π — соответствующие моменты весов обнуляются по симметрии, — что закреплено тестом.

Моментная стенка перенесена на GPU (раньше её там не было вовсе, бэкенд отказывался
запускаться): буфер индекса граничных узлов плюс ядро k_moment_wall. Скорость на момент t
берётся из пост-столкновительного поля — столкновение сохраняет ρ и ρu, поэтому отдельное
хранилище прошлого шага не нужно ни на одном бэкенде. Добавлена проверка лимита
storage-биндингов адаптера: моментной стенке нужно 9 против 8 гарантируемых.

Паритет на HRR: Cd 1.847 (CPU) против 1.848 (GPU). Три модели стенки на одной постановке дают
1.847 (hrr) / 1.845 (grad) / 1.856 (bouzidi).

33 теста, обе сборки чисты.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

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(),
         }
     }

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 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 };
+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<storage, read>       rest   : array<vec2<u32>>;
 @group(2) @binding(6) var<uniform>             A      : AmrParams;
 
+// (узел, смещение первого линка, число линков, —). Линки одного узла лежат подряд.
+@group(3) @binding(0) var<storage, read>       wnodes : array<vec4<u32>>;
+
 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<f32> {
+    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<f32>(t, cc);
+}
+
 // ───── физика: построчный перенос math.rs ─────
 
 // Энтропийное равновесие в product-form (точный максимизатор энтропии при заданных rho, rho*u)
@@ -369,7 +404,9 @@ fn k_force(@builtin(local_invocation_id) lid: vec3<u32>) {
     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 cy = array<f32,9>(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<f32>(0.0, 0.0);
+    var sy = vec2<f32>(0.0, 0.0);
+    var sz = vec2<f32>(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<u32>) {
         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;
+        sx = kadd(sx.x, sx.y, dfx);
-        sy = sy + dfy;
+        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<u32>) {
         s = s >> 1u;
     }
     if (t == 0u) {
-        results[16u + S.idx*4u + 0u] = wfx[0];
+        results[FORCE_BASE + S.idx*4u + 0u] = wfx[0];
-        results[16u + S.idx*4u + 1u] = wfy[0];
+        results[FORCE_BASE + S.idx*4u + 1u] = wfy[0];
-        results[16u + S.idx*4u + 2u] = wtz[0];
+        results[FORCE_BASE + S.idx*4u + 2u] = wtz[0];
     }
 }
 
@@ -440,7 +477,7 @@ fn k_stats1(@builtin(global_invocation_id) gid: vec3<u32>,
         var fv: array<f32,9>;
         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<u32>,
 @compute @workgroup_size(64)
 fn k_stats2(@builtin(local_invocation_id) lid: vec3<u32>) {
     let t = lid.x;
-    var p: Partial;
+    // слагаемых до 10^5 на поток — здесь и живёт основная потеря f32, поэтому все
-    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;
+    var rho = vec2<f32>(0.0, 0.0);
+    var gsum = vec2<f32>(0.0, 0.0);
+    var cnt = vec2<f32>(0.0, 0.0);
+    var degen = vec2<f32>(0.0, 0.0);
+    var xineg = vec2<f32>(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;
+        rho   = kadd(rho.x, rho.y, o.rho);
-        p.maxu = max(p.maxu, o.maxu);
+        gsum  = kadd(gsum.x, gsum.y, o.gsum);
-        p.gsum = p.gsum + o.gsum;
+        cnt   = kadd(cnt.x, cnt.y, o.cnt);
-        p.gmin = min(p.gmin, o.gmin);
+        degen = kadd(degen.x, degen.y, o.degen);
-        p.gmax = max(p.gmax, o.gmax);
+        xineg = kadd(xineg.x, xineg.y, o.xineg);
-        p.cnt = p.cnt + o.cnt;
+        maxu = max(maxu, o.maxu);
-        p.degen = p.degen + o.degen;
+        gmin = min(gmin, o.gmin);
-        p.xineg = p.xineg + o.xineg;
+        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;
+        var qrho = vec2<f32>(0.0, 0.0);
-        q.rho = 0.0; q.maxu = 0.0; q.gsum = 0.0; q.gmin = 1e30; q.gmax = -1e30;
+        var qgsum = vec2<f32>(0.0, 0.0);
-        q.cnt = 0.0; q.degen = 0.0; q.xineg = 0.0;
+        var qcnt = vec2<f32>(0.0, 0.0);
+        var qdeg = vec2<f32>(0.0, 0.0);
+        var qxin = vec2<f32>(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;
+            qrho  = kadd(qrho.x, qrho.y, o.rho);
-            q.maxu = max(q.maxu, o.maxu);
+            qgsum = kadd(qgsum.x, qgsum.y, o.gsum);
-            q.gsum = q.gsum + o.gsum;
+            qcnt  = kadd(qcnt.x, qcnt.y, o.cnt);
-            q.gmin = min(q.gmin, o.gmin);
+            qdeg  = kadd(qdeg.x, qdeg.y, o.degen);
-            q.gmax = max(q.gmax, o.gmax);
+            qxin  = kadd(qxin.x, qxin.y, o.xineg);
-            q.cnt = q.cnt + o.cnt;
+            qmaxu = max(qmaxu, o.maxu);
-            q.degen = q.degen + o.degen;
+            qgmin = min(qgmin, o.gmin);
-            q.xineg = q.xineg + o.xineg;
+            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];
+            fx = fx + results[FORCE_BASE + s*4u + 0u];
-            fy = fy + results[16u + s*4u + 1u];
+            fy = fy + results[FORCE_BASE + s*4u + 1u];
-            tz = tz + results[16u + s*4u + 2u];
+            tz = tz + results[FORCE_BASE + s*4u + 2u];
         }
         let inv = 1.0 / f32(D.refine);
-        results[0] = fx*inv;
+        let o = D.slot*12u;
-        results[1] = fy*inv;
+        results[o + 0u] = fx*inv;
-        results[2] = tz*inv;
+        results[o + 1u] = fy*inv;
-        results[4] = q.rho;
+        results[o + 2u] = tz*inv;
-        results[5] = q.maxu;
+        results[o + 4u] = qrho.x + qrho.y;
-        results[6] = q.gsum;
+        results[o + 5u] = qmaxu;
-        results[7] = q.gmin;
+        results[o + 6u] = qgsum.x + qgsum.y;
-        results[8] = q.gmax;
+        results[o + 7u] = qgmin;
-        results[9] = q.cnt;
+        results[o + 8u] = qgmax;
-        results[10] = q.degen;
+        results[o + 9u] = qcnt.x + qcnt.y;
-        results[11] = q.xineg;
+        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<f32> {
+    var ff: array<f32,9>;
+    for (var i = 0u; i < 9u; i = i + 1u) { ff[i] = post[i*P.n + nd]; }
+    let m = macros9(ff);
+    return vec2<f32>(m.y, m.z);
+}
+
+// (du/dx, du/dy, dv/dx, dv/dy): центральная разность там, где оба соседа жидкие,
+// односторонняя — где один твёрдый, ноль — если твёрдые оба.
+fn grad_u_at_t(nd: u32) -> vec4<f32> {
+    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<f32>(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<f32>(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<f32>(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<f32,9> {
+    let axx = pxx - rho*CS2;
+    let ayy = pyy - rho*CS2;
+    var w = array<f32,9>(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<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 out: array<f32,9>;
+    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<f32,9> {
+    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<f32,9>(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<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);
+        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<u32>) {
+    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<u32,9>(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<u32,9>(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, &params),
         ],
     });
-    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<u64>,
 }
 
 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 нельзя — это
+        // Моментная стенка добавляет к связке уровня ещё один storage-биндинг (индекс
-        // была бы другая физика под тем же ключом, поэтому отказываемся явно.
+        // граничных узлов), итого 9 против 8, гарантируемых лимитами по умолчанию. На
-        if spec.wall == math::WallModel::Grad {
+        // десктопных адаптерах их обычно >= 16, но проверить надо явно.
-            return Err("--wall grad на GPU не реализован (перенесён только Bouzidi). Запусти с --backend cpu либо оставь --wall bouzidi"
+        if spec.wall.is_moment_based() {
-                .into());
+            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::<Results>()) 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::<Results>() 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<GGhost> = 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<StepRec>) {
         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);
+            if moment_wall {
-            p.dispatch_workgroups(self.l0.link_groups(), 1, 1);
+                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);
+                    if moment_wall {
-                    p.dispatch_workgroups(nlink1, 1, 1);
+                        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.queue.submit(Some(enc.finish()));
-            &self.results_buf,
+        self.pending.push(t);
-            0,
+        self.step_index += 1;
-            &self.staging,
+        if self.pending.len() == HIST {
-            0,
+            self.flush(out);
-            std::mem::size_of::<Results>() as u64,
+        }
-        );
+    }
+
+    /// Дочитать всё, что уже посчитано устройством. Обязательно вызывать перед тем, как
+    /// смотреть на поля (кадр гифки, проверка на NaN) — иначе записи отстанут от состояния.
+    pub fn flush(&mut self, out: &mut Vec<StepRec>) {
+        let n = self.pending.len();
+        if n == 0 {
+            return;
+        }
+        let bytes = (n * std::mem::size_of::<Results>()) 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();
+        let slice = self.staging.slice(..bytes);
-        self.step_index += 1;
+        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<Results> = match rx.recv() {
+            Ok(Ok(())) => {
+                let data = slice.get_mapped_range();
+                bytemuck::cast_slice::<u8, Results>(&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::<Results>(&data[..std::mem::size_of::<Results>()])
-            }
-            _ => Results::default(),
-        };
-        self.staging.unmap();
-        out
-    }
-
     /// Скачать популяции L0 на хост (нужно для кадров и проверки на NaN).
     fn download_l0(&self) -> Vec<f32> {
         self.download(&self.l0.f, (9 * self.l0.n * 4) as u64)

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 — условие
+    /// Модель стенки на теле. hrr (умолчание) — восстановление недостающих популяций по
-    /// Града (Dorschner и др., JFM 801 (2016)), где задаются не популяции, а целевые моменты
+    /// целевым моментам с рекурсивной регуляризацией до 3-го порядка Эрмита; grad — то же,
-    /// ρ, u и тензор давлений. Авторы метода предпочитают grad: интерполяционные схемы, по их
+    /// но с обрывом на тензоре давлений; bouzidi — интерполированный отскок по каждому линку
-    /// словам, «ограничены низкими Re, поскольку на границе возникают паразитные скачки».
+    /// со своей долей пересечения q; staircase — простой отскок, q игнорируется (не для
-    #[arg(long, default_value = "bouzidi", value_parser = math::WallModel::ALL,
+    /// счёта, а как база сравнения). Замерено: по разрешению геометрии 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<StepRec>) {
         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<StepRec>) {
+        let _ = &out; // процессорный бэкенд отдаёт записи сразу, копить нечего
+        match self {
+            Backend::Cpu(_) => {}
+            #[cfg(feature = "gpu")]
+            Backend::Gpu(s) => s.flush(out),
         }
     }
     fn sample_field(&mut self, k: FieldKind) -> (Vec<R>, Vec<bool>) {
@@ -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();
+        back.advance(&mut recs);
-        recs.push(rec);
 
-        if let Some(w) = writer.as_mut() {
+        let need_frame = writer.is_some() && t % plan.stride == 0;
-            if 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 need_report {
-        if cli.report_every > 0 && !quiet && t % cli.report_every == 0 {
             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}",

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 nxx = -pref * 2.0 * dudx;
-    let pyy = rho * CS2 + rho * uy * uy - pref * 2.0 * dvdy;
+    let nyy = -pref * 2.0 * dvdy;
-    let pxy = rho * ux * uy - pref * (dudy + dvdx);
+    let nxy = -pref * (dudy + dvdx);
-    grad_init(rho, ux, uy, pxx, pxy, pyy)
+
+    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]);