//! БЭКЕНД ПОД ПРОЦЕССОР. Раскладка памяти, обход сетки, AMR-связка уровней и параллелизм
//! через rayon. Вся физика берётся из `math` — здесь её нет, есть только организация счёта.
//!
//! Раскладка: AoS, `Vec<[R; Q]>`, узел = y·nx + x. Столкновение — самая дорогая часть шага
//! (sqrt на узел) и в AoS оно тривиально распараллеливается и идеально ложится в кэш;
//! перенос при этом собирает 9 значений с 9 разных узлов, что дешевле, чем кажется, т.к.
//! соседи по x лежат рядом. GPU-бэкенд использует SoA — там важнее коалесцированный доступ.

use rayon::prelude::*;

use crate::math::{self, Kbc, Link, LinkKind, Scene, WallModel, MAX_BODY_BUCKETS, R, CX, CY, OPP, Q};
use crate::{Case, Collision, FieldKind, Spec, StepRec};

// ─────────────────────────────────────────────────────────────────────────────
//  Стартовые поля
// ─────────────────────────────────────────────────────────────────────────────

/// Что нужно, чтобы разложить стартовое поле уровня.
pub struct Init<'a> {
    pub case: Case,
    pub nx: usize,
    pub ny: usize,
    /// Скорость набегающего потока (канал) либо амплитуда эталонного течения.
    pub u0: (R, R),
    pub beta: R,
    pub scene: Option<&'a Scene>,
    /// Ширина сведения скорости к нулю у тела; 0 — не сводить.
    pub taper: R,
}

/// Стартовое поле в виде популяций. Строится на хосте для ОБОИХ бэкендов: GPU просто
/// перекладывает результат в свою раскладку. Так эталонные течения из статей задаются
/// ровно одним куском кода.
pub fn initial_field(p: &Init) -> Vec<[R; Q]> {
    let (nx, ny) = (p.nx, p.ny);
    let n = nx * ny;
    let tau = 1.0 / (2.0 * p.beta);
    let nu = math::nu_of_beta(p.beta);
    let two_pi = 2.0 * std::f64::consts::PI;
    let mut f = vec![[0.0; Q]; n];

    match p.case {
        Case::Channel => {
            let plain = math::feq(1.0, p.u0.0, p.u0.1);
            for (k, cell) in f.iter_mut().enumerate() {
                *cell = match (p.scene, p.taper > 0.0) {
                    (Some(sc), true) => {
                        let (x, y) = ((k % nx) as R, (k / nx) as R);
                        let c = math::wall_taper(sc.sdf(x, y), p.taper);
                        if c < 1.0 {
                            math::feq(1.0, p.u0.0 * c, p.u0.1 * c)
                        } else {
                            plain
                        }
                    }
                    _ => plain,
                };
            }
        }

        // Разд. VI 2D-статьи: u = ∇×[(u₀/k₂)cos(k₁x)cos(k₂y)·exp(−ν(k₁²+k₂²)t)], k₁=1, k₂=4.
        // Стартуем приближением Града и с СОГЛАСОВАННЫМ давлением: течение несёт собственное
        // поле p ~ ρu₀², и старт с ρ ≡ 1 сбрасывает эту разницу в акустику, которая в
        // периодическом ящике почти не затухает.
        Case::TaylorGreen => {
            let (k1, k2) = (1.0, 4.0);
            let kk1 = two_pi * k1 / nx as R;
            let kk2 = two_pi * k2 / ny as R;
            let a = p.u0.0;
            for y in 0..ny {
                for x in 0..nx {
                    let (px, py) = (kk1 * x as R, kk2 * y as R);
                    let ux = -a * px.cos() * py.sin();
                    let uy = (k1 / k2) * a * px.sin() * py.cos();
                    let pr = -(a * a / 4.0)
                        * ((2.0 * px).cos() + (k1 * k1) / (k2 * k2) * (2.0 * py).cos());
                    let rho = 1.0 + pr / math::CS2;
                    let dxux = a * kk1 * px.sin() * py.sin();
                    let dyux = -a * kk2 * px.cos() * py.cos();
                    let dxuy = (k1 / k2) * a * kk1 * px.cos() * py.cos();
                    let dyuy = -(k1 / k2) * a * kk2 * px.sin() * py.sin();
                    let pref = -tau * rho * math::CS2;
                    f[y * nx + x] = math::grad_init(
                        rho,
                        ux,
                        uy,
                        rho * math::CS2 + rho * ux * ux + pref * 2.0 * dxux,
                        rho * ux * uy + pref * (dyux + dxuy),
                        rho * math::CS2 + rho * uy * uy + pref * 2.0 * dyuy,
                    );
                }
            }
        }

        // Разд. VII: дважды периодический сдвиговый слой, κ = 80, δ = 0.05.
        Case::ShearLayer => {
            let (kappa, delta) = (80.0, 0.05);
            let a = p.u0.0;
            for y in 0..ny {
                let yy = y as R / ny as R;
                let ux = if y as R <= ny as R / 2.0 {
                    a * (kappa * (yy - 0.25)).tanh()
                } else {
                    a * (kappa * (0.75 - yy)).tanh()
                };
                for x in 0..nx {
                    let uy = delta * a * (two_pi * (x as R / nx as R + 0.25)).sin();
                    f[y * nx + x] = math::feq(1.0, ux, uy);
                }
            }
        }

        // Затухающая турбулентность: функция тока собирается из мод в узкой полосе вокруг
        // k₀, фазы — детерминированный ГПСЧ (прогон обязан быть воспроизводимым).
        Case::DecayingTurbulence => {
            let k0 = 8.0;
            let a = p.u0.0;
            let mut seed = 0x2545_F491_4F6C_DD1Du64;
            let mut rnd = || {
                seed ^= seed << 13;
                seed ^= seed >> 7;
                seed ^= seed << 17;
                (seed >> 11) as R / (1u64 << 53) as R
            };
            let mut modes = Vec::new();
            for ky in -12i32..=12 {
                for kx in -12i32..=12 {
                    if kx == 0 && ky == 0 {
                        continue;
                    }
                    let km = ((kx * kx + ky * ky) as R).sqrt();
                    if !(k0 - 4.0..=k0 + 4.0).contains(&km) || kx < 0 {
                        continue;
                    }
                    let amp = (-(km - k0) * (km - k0) / 8.0).exp() / km;
                    modes.push((kx as R, ky as R, amp, two_pi * rnd()));
                }
            }
            let norm: R = modes.iter().map(|m| m.2 * m.2).sum::<R>().sqrt().max(1e-30);
            for y in 0..ny {
                for x in 0..nx {
                    let (fx, fy) = (x as R / nx as R, y as R / ny as R);
                    let (mut ux, mut uy) = (0.0, 0.0);
                    for &(kx, ky, amp, ph) in &modes {
                        let arg = two_pi * (kx * fx + ky * fy) + ph;
                        let c = arg.cos() * amp / norm;
                        // u = ∇×(ψ ẑ) для ψ = Σ amp·sin(2π(k·r)+φ)
                        ux += c * two_pi * ky;
                        uy -= c * two_pi * kx;
                    }
                    f[y * nx + x] = math::feq(1.0, a * ux, a * uy);
                }
            }
            // приводим к заданной среднеквадратичной скорости
            let rms = (f
                .iter()
                .map(|c| {
                    let (_, ux, uy) = math::macros(c);
                    ux * ux + uy * uy
                })
                .sum::<R>()
                / n as R)
                .sqrt()
                .max(1e-30);
            let scale = a / rms;
            for cell in f.iter_mut() {
                let (_, ux, uy) = math::macros(cell);
                *cell = math::feq(1.0, ux * scale, uy * scale);
            }
        }
    }
    let _ = nu;
    f
}

// ─────────────────────────────────────────────────────────────────────────────
//  Геометрия уровня
// ─────────────────────────────────────────────────────────────────────────────

/// Маски, SDF и предсобранные Bouzidi-линки одного уровня сетки.
pub struct Geom {
    pub nx: usize,
    pub ny: usize,
    pub solid: Vec<bool>,
    /// Линки «жидкий → твёрдый» с готовой геометрией пересечения.
    pub links: Vec<Link>,
    /// Индекс по граничным узлам: условие Града ставится сразу на весь узел, а не полинково.
    pub wall_nodes: Vec<WallNode>,
    /// Плечи для момента: координаты узлов относительно центра тела.
    pub body_cx: R,
    pub body_cy: R,
}

impl Geom {
    /// Собрать геометрию уровня по сцене. `solid` определяется знаком SDF сцены (минимум по
    /// телам): φ ≤ 0 — тело. Каждому линку проставляется индекс тела, которого он касается,
    /// чтобы сила считалась по телам отдельно.
    pub fn build(nx: usize, ny: usize, scene: &Scene) -> Self {
        let n = nx * ny;
        let mut phi = vec![0.0; n];
        let mut solid = vec![false; n];
        let mut owner = vec![0u8; n];
        for y in 0..ny {
            for x in 0..nx {
                let k = y * nx + x;
                let p = scene.sdf(x as R, y as R);
                phi[k] = p;
                solid[k] = p <= 0.0;
                if solid[k] {
                    owner[k] =
                        scene.nearest(x as R, y as R).min(MAX_BODY_BUCKETS - 1) as u8;
                }
            }
        }
        let links = build_links(nx, ny, &solid, &phi, &owner);
        let wall_nodes = build_wall_nodes(&links);
        let (bcx, bcy) = scene.center();
        Geom { nx, ny, solid, links, wall_nodes, body_cx: bcx, body_cy: bcy }
    }
}

/// Индекс соседа по направлению i с периодическим заворотом (перенос сделан через roll,
/// поэтому решётка периодична по обеим осям; ГУ снимают заворот там, где он не физичен).
#[inline]
fn nb(nx: usize, ny: usize, x: usize, y: usize, i: usize, sign: i32) -> usize {
    let xs = (x as i32 + sign * CX[i]).rem_euclid(nx as i32) as usize;
    let ys = (y as i32 + sign * CY[i]).rem_euclid(ny as i32) as usize;
    ys * nx + xs
}

/// Для каждого жидкого узла и каждого направления, упирающегося в тело, решаем, какая из трёх
/// формул Bouzidi применима, и считаем долю пересечения q из SDF. Делается один раз.
fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R], owner: &[u8]) -> Vec<Link> {
    let mut links = Vec::new();
    for y in 0..ny {
        for x in 0..nx {
            let node = y * nx + x;
            if solid[node] {
                continue;
            }
            for i in 1..Q {
                let s = nb(nx, ny, x, y, i, 1); // сосед по +c_i
                if !solid[s] {
                    continue;
                }
                let far = nb(nx, ny, x, y, i, -1); // «дальний» сосед x_f − c_i
                let q = math::bouzidi_q(phi[node], phi[s]);
                let kind = if q >= 0.5 {
                    LinkKind::Far
                } else if !solid[far] {
                    LinkKind::Near
                } else {
                    LinkKind::Simple
                };
                links.push(Link {
                    node: node as u32,
                    far: far as u32,
                    i: i as u8,
                    ib: OPP[i] as u8,
                    kind,
                    q,
                    body: owner[s],
                });
            }
        }
    }
    links
}

/// Сводка по границе тела: действительно ли работает субсеточная интерполяция.
#[derive(Clone, Copy, Debug, Default)]
pub struct WallStats {
    pub nodes: usize,
    pub links: usize,
    /// q < ½ и есть дальний жидкий сосед — полная интерполяционная формула.
    pub near: usize,
    /// q ≥ ½ — интерполяция по двум популяциям того же узла.
    pub far: usize,
    /// Дальнего жидкого соседа нет: откат на ПРОСТОЙ отскок, то есть ступенчатую стенку.
    pub simple: usize,
    pub q_min: R,
    pub q_max: R,
    pub q_mean: R,
    /// Линки, где q упёрлась в клип [0.02, 0.98] — узел практически лежит на поверхности.
    pub q_clipped: usize,
}

impl Geom {
    pub fn wall_stats(&self) -> WallStats {
        let mut w = WallStats { q_min: R::INFINITY, q_max: R::NEG_INFINITY, ..Default::default() };
        w.nodes = self.wall_nodes.len();
        w.links = self.links.len();
        for l in &self.links {
            match l.kind {
                LinkKind::Near => w.near += 1,
                LinkKind::Far => w.far += 1,
                LinkKind::Simple => w.simple += 1,
            }
            w.q_min = w.q_min.min(l.q);
            w.q_max = w.q_max.max(l.q);
            w.q_mean += l.q;
            if l.q <= 0.0201 || l.q >= 0.9799 {
                w.q_clipped += 1;
            }
        }
        if w.links > 0 {
            w.q_mean /= w.links as R;
        }
        w
    }
}

/// Жидкий узел, в который после переноса пришла бы хоть одна популяция из тела.
/// Линки одного узла в `Geom::links` лежат подряд (сборка идёт по y, потом x, потом
/// направлениям), поэтому достаточно смещения и количества.
#[derive(Clone, Copy)]
pub struct WallNode {
    pub node: u32,
    pub first: u32,
    pub count: u8,
}

fn build_wall_nodes(links: &[Link]) -> Vec<WallNode> {
    let mut out: Vec<WallNode> = Vec::new();
    for (k, l) in links.iter().enumerate() {
        match out.last_mut() {
            Some(w) if w.node == l.node => w.count += 1,
            _ => out.push(WallNode { node: l.node, first: k as u32, count: 1 }),
        }
    }
    out
}

// ─────────────────────────────────────────────────────────────────────────────
//  Уровень сетки
// ─────────────────────────────────────────────────────────────────────────────

/// Один уровень: поля популяций до/после столкновения, β и геометрия.
pub struct Level {
    pub nx: usize,
    pub ny: usize,
    /// Текущее состояние (после переноса и ГУ).
    pub f: Vec<[R; Q]>,
    /// Состояние после столкновения — источник переноса и вход Bouzidi/GMEM.
    pub post: Vec<[R; Q]>,
    /// Энтропийный стабилизатор поузлово (для диагностики и картинки).
    pub gamma: Vec<R>,
    /// β по столбцам x (губка делает его полем); длина nx.
    pub beta: Vec<R>,
    pub geom: Geom,
}

impl Level {
    /// `u0` — скорость, которой заполняется поле на старте. Заполнять сразу набегающим
    /// потоком принципиально: старт из покоя разгоняет весь столб жидкости и закачивает в
    /// канал продольную акустическую моду, которую вязкость потом почти не гасит.
    fn new(nx: usize, ny: usize, beta: Vec<R>, geom: Geom, f0: Vec<[R; Q]>) -> Self {
        let n = nx * ny;
        debug_assert_eq!(f0.len(), n);
        Level {
            nx,
            ny,
            post: f0.clone(),
            f: f0,
            gamma: vec![2.0; n],
            beta,
            geom,
        }
    }

    /// Столкновение по всем жидким узлам. Внутри тела не считаем: эти популяции фиктивны
    /// (ГУ Bouzidi перезаписывает всё, что могло бы прийти из тела в жидкость), а счёт там
    /// только жжёт такты и способен родить NaN при экстремальных режимах.
    fn collide(&mut self, op: Collision, model: math::KbcModel) -> Stats {
        let nx = self.nx;
        let beta = &self.beta;
        let solid = &self.geom.solid;
        let post = &mut self.post;
        let gamma = &mut self.gamma;
        let f = &self.f;

        post.par_iter_mut()
            .zip(gamma.par_iter_mut())
            .enumerate()
            .map(|(node, (p, g))| {
                *p = f[node];
                if solid[node] {
                    *g = 2.0;
                    return Stats::EMPTY;
                }
                let b = beta[node % nx];
                let k: Kbc = match op {
                    Collision::Kbc => math::collide_node(p, b, model),
                    Collision::Bgk => math::collide_node_bgk(p, b),
                };
                *g = k.gamma;
                Stats::of(k, b, model)
            })
            .reduce(|| Stats::EMPTY, Stats::merge)
    }

    /// Перенос f_i(x, t+1) = f_i(x − c_i, t), периодический по обеим осям.
    fn stream(&mut self) {
        let (nx, ny) = (self.nx, self.ny);
        let post = &self.post;
        self.f
            .par_chunks_mut(nx)
            .enumerate()
            .for_each(|(y, row)| {
                for (x, cell) in row.iter_mut().enumerate() {
                    for i in 0..Q {
                        let xs = (x as i32 - CX[i]).rem_euclid(nx as i32) as usize;
                        let ys = (y as i32 - CY[i]).rem_euclid(ny as i32) as usize;
                        cell[i] = post[ys * nx + xs][i];
                    }
                }
            });
        let _ = ny;
    }

    /// Простой отскок по всем линкам: q игнорируется. База сравнения для субсеточных моделей.
    fn apply_staircase(&mut self) {
        for l in &self.geom.links {
            let node = l.node as usize;
            self.f[node][l.ib as usize] = self.post[node][l.i as usize];
        }
    }

    /// Интерполированный отскок Bouzidi по всем линкам тела.
    /// `f` — поле ПОСЛЕ переноса (его правим), `post` — ПОСЛЕ столкновения (до переноса).
    fn apply_bouzidi(&mut self) {
        for l in &self.geom.links {
            let node = l.node as usize;
            let i = l.i as usize;
            let ib = l.ib as usize;
            let fi = self.post[node][i];
            let v = match l.kind {
                LinkKind::Near => 2.0 * l.q * fi + (1.0 - 2.0 * l.q) * self.post[l.far as usize][i],
                LinkKind::Far => {
                    let h = 1.0 / (2.0 * l.q);
                    h * fi + (1.0 - h) * self.post[node][ib]
                }
                LinkKind::Simple => fi,
            };
            self.f[node][ib] = v;
        }
    }

    /// ГРАНИЧНОЕ УСЛОВИЕ НА МОМЕНТАХ (Dorschner и др., JFM 801 (2016), прил. B).
    ///
    /// `third_order` — продолжать ли ряд Эрмита рекурсивно до третьего порядка (HRR).
    ///
    /// В отличие от Bouzidi, работающего полинково, здесь условие ставится сразу на весь
    /// граничный узел и не на популяции, а на моменты:
    ///
    ///   целевая скорость  u_tgt = (1/n) Σ (q_i·u_f,i + u_w,i)/(1 + q_i)              (B 1)
    ///   целевая плотность ρ_tgt = Σ_известные f_i + Σ_недостающие f_i^отскок          (B 3)
    ///   тензор давлений   Π = ρc_s²I + ρu⊗u − (ρc_s²/2β)(∇u + ∇uᵀ)              (2.14)–(2.16)
    ///
    /// после чего недостающие популяции собираются приближением Града (2.13). Скорости
    /// соседей и градиенты берутся с прошлого шага — см. `u_prev`.
    fn apply_moment_wall(&mut self, third_order: bool) {
        let nx = self.nx;
        // стенка неподвижна; для подвижного тела сюда пойдёт её скорость на линке,
        // и добавится динамическая часть плотности (B 4)
        let (uwx, uwy) = (0.0, 0.0);
        for wi in 0..self.geom.wall_nodes.len() {
            let w = self.geom.wall_nodes[wi];
            let node = w.node as usize;
            let lo = w.first as usize;
            let hi = lo + w.count as usize;

            // ── целевая скорость (B 1) ──
            let (mut ux, mut uy) = (0.0, 0.0);
            for l in &self.geom.links[lo..hi] {
                let far = l.far as usize;
                let (fx, fy) = if self.geom.solid[far] {
                    // дальнего жидкого соседа нет (тело тоньше двух клеток) — остаётся стенка
                    (uwx, uwy)
                } else {
                    self.u_at_t(far)
                };
                ux += math::grad_target_velocity_term(l.q, fx, uwx);
                uy += math::grad_target_velocity_term(l.q, fy, uwy);
            }
            let inv = 1.0 / w.count as R;
            ux *= inv;
            uy *= inv;

            // ── целевая плотность (B 3): известные популяции плюс отскок недостающих ──
            let mut missing = [false; Q];
            for l in &self.geom.links[lo..hi] {
                missing[l.ib as usize] = true;
            }
            let mut rho = 0.0;
            for i in 0..Q {
                rho += if missing[i] { self.post[node][OPP[i]] } else { self.f[node][i] };
            }

            let (dudx, dudy, dvdx, dvdy) = self.grad_u_at_t(node);
            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];
                }
            }
        }
    }

    /// Скорость узла на момент t — берётся из ПОСТ-СТОЛКНОВИТЕЛЬНОГО поля.
    ///
    /// Столкновение сохраняет ρ и ρu точно, поэтому macros(post) даёт ровно ту же скорость,
    /// что было в f до переноса, то есть u(x, t) — именно то, что требует прил. B. Отдельное
    /// хранилище «поля предыдущего шага» при этом не нужно.
    #[inline]
    fn u_at_t(&self, node: usize) -> (R, R) {
        let (_, a, b) = math::macros(&self.post[node]);
        (a, b)
    }

    /// Градиент скорости на момент t: центральная разность там, где оба соседа жидкие,
    /// односторонняя — где один твёрдый, ноль — если твёрдые оба.
    /// Возвращает (∂u_x/∂x, ∂u_x/∂y, ∂u_y/∂x, ∂u_y/∂y).
    fn grad_u_at_t(&self, node: usize) -> (R, R, R, R) {
        let (nx, ny) = (self.nx, self.ny);
        let (x, y) = (node % nx, node / nx);
        let xm = y * nx + (x + nx - 1) % nx;
        let xp = y * nx + (x + 1) % nx;
        let ym = ((y + ny - 1) % ny) * nx + x;
        let yp = ((y + 1) % ny) * nx + x;
        let solid = &self.geom.solid;
        let uc = self.u_at_t(node);
        let comp = |t: (R, R), c: usize| if c == 0 { t.0 } else { t.1 };
        let d = |a: usize, b: usize, c: usize| -> R {
            match (!solid[a], !solid[b]) {
                (true, true) => 0.5 * (comp(self.u_at_t(b), c) - comp(self.u_at_t(a), c)),
                (false, true) => comp(self.u_at_t(b), c) - comp(uc, c),
                (true, false) => comp(uc, c) - comp(self.u_at_t(a), c),
                (false, false) => 0.0,
            }
        };
        (d(xm, xp, 0), d(ym, yp, 0), d(xm, xp, 1), d(ym, yp, 1))
    }

    /// Замкнуть недостающие популяции выбранной моделью стенки.
    fn apply_wall(&mut self, model: WallModel) {
        match model {
            WallModel::Bouzidi => self.apply_bouzidi(),
            WallModel::Grad => self.apply_moment_wall(false),
            WallModel::Hrr => self.apply_moment_wall(true),
            WallModel::Staircase => self.apply_staircase(),
        }
    }

    /// FREE-SLIP (зеркальные) стенки канала сверху и снизу: касательный импульс сохраняется,
    /// нормальный заворачивается, масса сохраняется. Применять ПОСЛЕ переноса и ПЕРЕД Zou–He.
    fn free_slip_walls(&mut self) {
        let (nx, ny) = (self.nx, self.ny);
        for x in 0..nx {
            let b = &mut self.f[x];
            b[2] = b[4];
            b[5] = b[8];
            b[6] = b[7];
        }
        let top = (ny - 1) * nx;
        for x in 0..nx {
            let t = &mut self.f[top + x];
            t[4] = t[2];
            t[7] = t[6];
            t[8] = t[5];
        }
    }

    /// Zou–He вход/выход по всему столбцу, включая угловые узлы.
    ///
    /// ПОРЯДОК КРИТИЧЕН: только ПОСЛЕ free_slip_walls и обязательно на ВЕСЬ столбец. Перенос
    /// сделан через roll и потому периодичен по обеим осям: заворот по y снимают стенки,
    /// заворот по x — этот ГУ. Если оставить ряды 0 и ny−1 без Zou–He, в них популяции с
    /// c_x = +1 на входе приходят прямо из столбца выхода, и вход с выходом оказываются
    /// физически связаны в четырёх узлах.
    fn channel_bc(&mut self, ux_in: R, uy_in: R, rho_out: R, outlet_extrapolate: bool) {
        let (nx, ny) = (self.nx, self.ny);
        for y in 0..ny {
            math::zou_he_inlet(&mut self.f[y * nx], ux_in, uy_in);
        }
        for y in 0..ny {
            let uy_out = if outlet_extrapolate {
                // нуль-градиент поперечной скорости: берём u_y с предвыходного столбца
                let c = &self.f[y * nx + nx - 2];
                let s: R = c.iter().sum();
                ((c[2] + c[5] + c[6]) - (c[4] + c[7] + c[8])) / s
            } else {
                0.0
            };
            math::zou_he_outlet(&mut self.f[y * nx + nx - 1], rho_out, uy_out);
        }
    }

    /// Сила и момент по GMEM, разложенные по телам: `out[b] = (F_x, F_y, T_z)` для тела b.
    /// Момент считается вокруг центра ПЕРВОГО тела — общей точки отсчёта для всей сцены.
    fn force(&self) -> [[R; 3]; MAX_BODY_BUCKETS] {
        let mut out = [[0.0; 3]; MAX_BODY_BUCKETS];
        let nx = self.nx;
        for l in &self.geom.links {
            let node = l.node as usize;
            let i = l.i as usize;
            let (dfx, dfy) =
                math::gmem_link(i, self.post[node][i], self.f[node][l.ib as usize], 0.0, 0.0);
            // плечо — до ТОЧКИ ПЕРЕСЕЧЕНИЯ линка со стенкой r_w = r_f + q·c_i, а не до узла:
            // узел дал бы ошибку плеча до целой ячейки, а q уже посчитан
            let rx = (node % nx) as R - self.geom.body_cx + l.q * CX[i] as R;
            let ry = (node / nx) as R - self.geom.body_cy + l.q * CY[i] as R;
            let b = (l.body as usize).min(MAX_BODY_BUCKETS - 1);
            out[b][0] += dfx;
            out[b][1] += dfy;
            out[b][2] += rx * dfy - ry * dfx;
        }
        out
    }

    #[inline]
    pub fn macros_at(&self, node: usize) -> (R, R, R) {
        math::macros(&self.f[node])
    }
}

// ─────────────────────────────────────────────────────────────────────────────
//  Статистика KBC за шаг
// ─────────────────────────────────────────────────────────────────────────────

/// Сводка по стабилизатору γ за один проход столкновения. Собирается редукцией rayon.
#[derive(Clone, Copy, Debug)]
pub struct Stats {
    pub n: u64,
    pub gsum: R,
    pub gmin: R,
    pub gmax: R,
    /// Узлы, где ⟨Δh|Δh⟩ выродилось и γ подменён на 2 (локально — чистый LBGK).
    pub degenerate: u64,
    /// Узлы с отрицательной объёмной вязкостью ξ = c_s²(1/(γβ) − ½) — локальное антизатухание.
    pub xi_negative: u64,
}

impl Stats {
    pub const EMPTY: Stats = Stats {
        n: 0,
        gsum: 0.0,
        gmin: R::INFINITY,
        gmax: R::NEG_INFINITY,
        degenerate: 0,
        xi_negative: 0,
    };
    #[inline]
    fn of(k: Kbc, beta: R, model: math::KbcModel) -> Stats {
        Stats {
            n: 1,
            gsum: k.gamma,
            gmin: k.gamma,
            gmax: k.gamma,
            degenerate: k.degenerate as u64,
            xi_negative: (math::xi_of_gamma(k.gamma, beta, model) < 0.0) as u64,
        }
    }
    fn merge(a: Stats, b: Stats) -> Stats {
        Stats {
            n: a.n + b.n,
            gsum: a.gsum + b.gsum,
            gmin: a.gmin.min(b.gmin),
            gmax: a.gmax.max(b.gmax),
            degenerate: a.degenerate + b.degenerate,
            xi_negative: a.xi_negative + b.xi_negative,
        }
    }
    pub fn gamma_mean(&self) -> R {
        if self.n == 0 {
            0.0
        } else {
            self.gsum / self.n as R
        }
    }
    pub fn degenerate_frac(&self) -> R {
        if self.n == 0 {
            0.0
        } else {
            self.degenerate as R / self.n as R
        }
    }
    pub fn xi_negative_frac(&self) -> R {
        if self.n == 0 {
            0.0
        } else {
            self.xi_negative as R / self.n as R
        }
    }
}

// ─────────────────────────────────────────────────────────────────────────────
//  AMR-патч
// ─────────────────────────────────────────────────────────────────────────────

/// Узел рамки тонкого уровня с готовым билинейным стенсилем с грубого уровня.
/// Публичен: GPU-бэкенд переиспользует ту же топологию патча, а не строит её заново.
#[derive(Clone, Copy)]
pub struct Ghost {
    pub fine: u32,
    pub c00: u32,
    pub c10: u32,
    pub c01: u32,
    pub c11: u32,
    pub tx: R,
    pub ty: R,
}

/// Связка «грубый L0 ↔ тонкий L1»: на каждый шаг L0 тонкий уровень делает r подшагов.
/// Рамка патча заполняется интерполяцией с L0 (билинейно по пространству и линейно по
/// времени между состояниями «до» и «после» шага L0), внутренность после подшагов
/// проецируется обратно на L0.
///
/// Неравновесная часть при смене уровня масштабируется: f^neq ∝ τ·δt, поэтому
/// коэффициент грубый→тонкий равен R01 = τ_f/(r·τ_c), обратно — 1/R01.
pub struct Patch {
    pub ax: usize,
    pub bx: usize,
    pub ay: usize,
    pub by: usize,
    pub r: usize,
    pub nfx: usize,
    pub r01: R,
    ghosts: Vec<Ghost>,
    /// Пары (узел L0, узел L1) для рестрикции — только внутренние жидкие узлы перекрытия.
    restrict: Vec<(u32, u32)>,
}

impl Patch {
    /// Топология рамки и рестрикции. Публична: GPU-бэкенд загружает эти же списки в буферы,
    /// вместо того чтобы строить их заново. В сборке без GPU читателей у них нет.
    #[cfg_attr(not(feature = "gpu"), allow(dead_code))]
    pub fn ghosts(&self) -> &[Ghost] {
        &self.ghosts
    }
    #[cfg_attr(not(feature = "gpu"), allow(dead_code))]
    pub fn restrict_pairs(&self) -> &[(u32, u32)] {
        &self.restrict
    }

    pub fn new(spec: &Spec, coarse: &Geom, fine_solid: &[bool], r01: R) -> Patch {
        let (ax, bx, ay, by) = spec.patch.expect("патч запрошен без границ");
        let r = spec.refine;
        let nfx = r * (bx - ax) + 1;
        let nfy = r * (by - ay) + 1;
        let cnx = coarse.nx;

        // рамка: один ряд по периметру тонкого поля
        let mut ghosts = Vec::with_capacity(2 * (nfx + nfy));
        let push = |gx: usize, gy: usize, out: &mut Vec<Ghost>| {
            // координата тонкого узла в системе L0
            let fx = ax as R + gx as R / r as R;
            let fy = ay as R + gy as R / r as R;
            let x0 = fx.floor() as usize;
            let y0 = fy.floor() as usize;
            let x1 = (x0 + 1).min(cnx - 1);
            let y1 = (y0 + 1).min(coarse.ny - 1);
            out.push(Ghost {
                fine: (gy * nfx + gx) as u32,
                c00: (y0 * cnx + x0) as u32,
                c10: (y0 * cnx + x1) as u32,
                c01: (y1 * cnx + x0) as u32,
                c11: (y1 * cnx + x1) as u32,
                tx: fx - x0 as R,
                ty: fy - y0 as R,
            });
        };
        for gx in 0..nfx {
            push(gx, 0, &mut ghosts);
            push(gx, nfy - 1, &mut ghosts);
        }
        for gy in 1..nfy - 1 {
            push(0, gy, &mut ghosts);
            push(nfx - 1, gy, &mut ghosts);
        }

        // рестрикция: каждый r-й тонкий узел внутренности, только там, где на L0 жидкость
        let mut restrict = Vec::new();
        for cy in ay + 1..by {
            for cx in ax + 1..bx {
                let cnode = cy * cnx + cx;
                if coarse.solid[cnode] {
                    continue;
                }
                let fnode = ((cy - ay) * r) * nfx + (cx - ax) * r;
                if fine_solid[fnode] {
                    continue;
                }
                restrict.push((cnode as u32, fnode as u32));
            }
        }
        // рамка — ровно периметр тонкого поля: два ряда по nfx плюс два столбца без углов
        assert_eq!(ghosts.len(), 2 * nfx + 2 * (nfy - 2));
        Patch { ax, bx, ay, by, r, nfx, r01, ghosts, restrict }
    }

    /// Ghost-значения из грубого поля: равновесие по интерполированным ρ, u плюс
    /// масштабированная неравновесная часть.
    fn ghost_values(&self, coarse: &[[R; Q]], out: &mut Vec<[R; Q]>) {
        out.clear();
        out.reserve(self.ghosts.len());
        for g in &self.ghosts {
            let w00 = (1.0 - g.tx) * (1.0 - g.ty);
            let w10 = g.tx * (1.0 - g.ty);
            let w01 = (1.0 - g.tx) * g.ty;
            let w11 = g.tx * g.ty;
            // Интерполируются ОТДЕЛЬНО макропеременные и отдельно неравновесная часть:
            // рамка = feq(интерполированные ρ, u) + R01·интерполированная neq. Так сделано в
            // эталонном python-решателе, и это не то же самое, что интерполировать сами
            // популяции: u = interp(ρu)/interp(ρ) отличается от interp(u) во втором порядке,
            // и расщепление на eq/neq тогда тоже смещается. Разница мала поузлово, но рамка
            // задаёт весь обмен между уровнями каждый шаг, так что копится.
            let mut rho = 0.0;
            let mut ux = 0.0;
            let mut uy = 0.0;
            let mut neq = [0.0; Q];
            for (w, c) in [
                (w00, g.c00 as usize),
                (w10, g.c10 as usize),
                (w01, g.c01 as usize),
                (w11, g.c11 as usize),
            ] {
                let cf = &coarse[c];
                let (r, x, y) = math::macros(cf);
                rho += w * r;
                ux += w * x;
                uy += w * y;
                let fe = math::feq(r, x, y);
                for i in 0..Q {
                    neq[i] += w * (cf[i] - fe[i]);
                }
            }
            let fe = math::feq(rho, ux, uy);
            let mut v = [0.0; Q];
            for i in 0..Q {
                v[i] = fe[i] + self.r01 * neq[i];
            }
            out.push(v);
        }
    }

    /// Записать рамку тонкого поля линейной комбинацией ghost-значений «до» и «после».
    fn fill(&self, fine: &mut [[R; Q]], old: &[[R; Q]], new: &[[R; Q]], w: R) {
        for (k, g) in self.ghosts.iter().enumerate() {
            let dst = &mut fine[g.fine as usize];
            for i in 0..Q {
                dst[i] = (1.0 - w) * old[k][i] + w * new[k][i];
            }
        }
    }

    /// Спроецировать внутренность тонкого уровня обратно на грубый, масштабируя neq на 1/R01.
    fn restrict_to(&self, fine: &[[R; Q]], coarse: &mut [[R; Q]]) {
        let rfc = 1.0 / self.r01;
        for &(cnode, fnode) in &self.restrict {
            let ff = &fine[fnode as usize];
            let (rho, ux, uy) = math::macros(ff);
            let fe = math::feq(rho, ux, uy);
            let dst = &mut coarse[cnode as usize];
            for i in 0..Q {
                dst[i] = fe[i] + rfc * (ff[i] - fe[i]);
            }
        }
    }
}

// ─────────────────────────────────────────────────────────────────────────────
//  Симуляция
// ─────────────────────────────────────────────────────────────────────────────

pub struct Sim {
    pub spec: Spec,
    pub l0: Level,
    pub l1: Option<Level>,
    pub patch: Option<Patch>,
    /// Состояние L0 до столкновения — «старый» край для временной интерполяции рамки.
    pre: Vec<[R; Q]>,
    gh_old: Vec<[R; Q]>,
    gh_new: Vec<[R; Q]>,
    fluid_count: R,
    probe_node: usize,
    probe_on_fine: bool,
    pub step_index: u64,
}

impl Sim {
    pub fn new(spec: Spec) -> Sim {
        let (nx, ny) = (spec.nx, spec.ny);
        let geom0 = Geom::build(nx, ny, &spec.scene);

        let beta0 = math::beta_profile(
            nx,
            spec.beta0,
            spec.sponge_in,
            spec.sponge_len,
            spec.sponge_mult,
        );

        // стартовое поле: либо сразу набегающий поток, либо покой
        let u0 = if spec.init_uniform {
            let (sn, cs) = spec.flow_angle.sin_cos();
            (spec.units.u_lat * cs, spec.units.u_lat * sn)
        } else {
            (0.0, 0.0)
        };
        let init0 = initial_field(&Init {
            case: spec.case,
            nx,
            ny,
            u0,
            beta: spec.beta0,
            scene: Some(&spec.scene),
            taper: if spec.init_uniform { spec.init_taper } else { 0.0 },
        });
        let l0 = Level::new(nx, ny, beta0, geom0, init0);
        let fluid_count = l0.geom.solid.iter().filter(|s| !**s).count() as R;

        let (l1, patch) = if spec.refine > 1 {
            let (ax, bx, ay, by) = spec.patch.expect("refine > 1 требует патч");
            let r = spec.refine;
            let scene1 = spec.scene.refined(r as R, ax as R, ay as R);
            let nfx = r * (bx - ax) + 1;
            let nfy = r * (by - ay) + 1;
            let geom1 = Geom::build(nfx, nfy, &scene1);
            // τ_f = r(τ_c − ½) + ½ ⇒ одинаковая ν на обоих уровнях
            let tau0 = 1.0 / (2.0 * spec.beta0);
            let tau1 = r as R * (tau0 - 0.5) + 0.5;
            let beta1 = 1.0 / (2.0 * tau1);
            let r01 = tau1 / (r as R * tau0);
            let patch = Patch::new(&spec, &l0.geom, &geom1.solid, r01);
            let init1 = initial_field(&Init {
                case: spec.case,
                nx: nfx,
                ny: nfy,
                u0,
                beta: 1.0 / (2.0 * tau1),
                scene: Some(&scene1),
                taper: if spec.init_uniform { spec.init_taper * r as R } else { 0.0 },
            });
            (Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1, init1)), Some(patch))
        } else {
            (None, None)
        };

        // зонд следа: берём с тонкой сетки, если точка внутри патча (меньше численного
        // размытия вихрей → чище спектр и St), иначе с грубой
        let (px, py) = spec.probe;
        let (probe_node, probe_on_fine) = match &patch {
            Some(p) if px >= p.ax && px <= p.bx && py >= p.ay && py <= p.by => {
                (((py - p.ay) * p.r) * p.nfx + (px - p.ax) * p.r, true)
            }
            _ => (py * nx + px, false),
        };

        let n = nx * ny;
        Sim {
            spec,
            l0,
            l1,
            patch,
            pre: vec![[0.0; Q]; n],
            gh_old: Vec::new(),
            gh_new: Vec::new(),
            fluid_count,
            probe_node,
            probe_on_fine,
            step_index: 0,
        }
    }

    /// Скорость на входе в момент t: разгон smoothstep плюс окно поперечного возмущения.
    ///
    /// При старте из однородного потока разгон не нужен и пропускается: поле и вход и так
    /// согласованы, а плавный разгон поверх согласованного поля сам стал бы рассогласованием.
    /// При старте из покоя разгон обязателен — мгновенное включение входа шлёт по домену
    /// ударную волну. Возмущение — короткий поперечный импульс, сбивающий симметрию: без него
    /// дорожка Кармана заводится только на численном шуме и стартует на порядок позже.
    fn inlet(&self, t: u64) -> (R, R) {
        let sp = &self.spec;
        let ramp = if sp.init_uniform {
            1.0
        } else {
            math::smoothstep(t as R / sp.ramp.max(1) as R)
        };
        let u = sp.units.u_lat * ramp;
        let (s, c) = sp.flow_angle.sin_cos();
        let mut uy = u * s;
        if sp.pert_dur > 0 && t >= sp.ramp && t < sp.ramp + sp.pert_dur {
            let ph = (t - sp.ramp) as R / sp.pert_dur as R;
            uy += sp.pert_amp * sp.units.u_lat * (std::f64::consts::PI * ph).sin();
        }
        (u * c, uy)
    }

    pub fn step(&mut self) -> StepRec {
        let t = self.step_index;
        let (ux_in, uy_in) = self.inlet(t);
        let sp_collision = self.spec.collision;
        let outlet_extrap = self.spec.outlet_extrapolate;

        // ── уровень 0 ──
        self.pre.copy_from_slice(&self.l0.f);
        let model = self.spec.kbc_model;
        let wall = self.spec.wall;
        let stats0 = self.l0.collide(sp_collision, model);
        self.l0.stream();
        // Эталонные течения статей периодичны по обеим осям и ГУ не имеют вовсе: перенос
        // уже периодичен, поэтому достаточно ничего не накладывать.
        if self.spec.case == Case::Channel {
            self.l0.apply_wall(wall);
            self.l0.free_slip_walls();
            self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap);
        }

        // ── уровень 1: r подшагов с временной интерполяцией рамки ──
        let mut fb = [[0.0 as R; 3]; MAX_BODY_BUCKETS];
        if let (Some(l1), Some(p)) = (self.l1.as_mut(), self.patch.as_ref()) {
            p.ghost_values(&self.pre, &mut self.gh_old);
            p.ghost_values(&self.l0.f, &mut self.gh_new);
            for s in 0..p.r {
                l1.collide(sp_collision, model);
                l1.stream();
                if self.spec.case == Case::Channel {
                    l1.apply_wall(wall);
                }
                // силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на
                // последнем подшаге даёт лишний шум в рядах при том же среднем
                let g = l1.force();
                for k in 0..MAX_BODY_BUCKETS {
                    for c in 0..3 {
                        fb[k][c] += g[k][c];
                    }
                }
                let w = (s + 1) as R / p.r as R;
                p.fill(&mut l1.f, &self.gh_old, &self.gh_new, w);
            }
            let inv = 1.0 / p.r as R;
            for k in 0..MAX_BODY_BUCKETS {
                for c in 0..3 {
                    fb[k][c] *= inv;
                }
            }
            p.restrict_to(&l1.f, &mut self.l0.f);
        } else {
            fb = self.l0.force();
        }

        // ── диагностика ──
        let probe = if self.probe_on_fine {
            self.l1.as_ref().unwrap().macros_at(self.probe_node)
        } else {
            self.l0.macros_at(self.probe_node)
        };

        let solid = &self.l0.geom.solid;
        let (rho_sum, max_u) = self
            .l0
            .f
            .par_iter()
            .enumerate()
            .filter(|(n, _)| !solid[*n])
            .map(|(_, c)| {
                let (r, ux, uy) = math::macros(c);
                (r, (ux * ux + uy * uy).sqrt())
            })
            .reduce(|| (0.0, 0.0), |a, b| (a.0 + b.0, a.1.max(b.1)));

        self.step_index += 1;
        let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0);
        for k in 0..MAX_BODY_BUCKETS {
            fx += fb[k][0];
            fy += fb[k][1];
            tz += fb[k][2];
        }
        StepRec {
            step: t,
            fx,
            fy,
            tz,
            body: fb,
            uy_probe: probe.2,
            rho_mean: rho_sum / self.fluid_count,
            max_u,
            gamma_mean: stats0.gamma_mean(),
            gamma_min: stats0.gmin,
            gamma_max: stats0.gmax,
            degenerate_frac: stats0.degenerate_frac(),
            xi_negative_frac: stats0.xi_negative_frac(),
        }
    }

    /// Характерный размер тела в единицах того уровня, где снимается сила.
    pub fn force_ref_size(&self) -> R {
        match &self.patch {
            Some(p) => self.spec.scene.ref_size() * p.r as R,
            None => self.spec.scene.ref_size(),
        }
    }

    /// Поле для картинки: (значение, маска тела) на сетке L0.
    pub fn sample_field(&self, kind: FieldKind) -> (Vec<R>, &[bool]) {
        let (nx, ny) = (self.l0.nx, self.l0.ny);
        let mut out = vec![0.0; nx * ny];
        match kind {
            FieldKind::Speed => {
                for n in 0..nx * ny {
                    let (_, ux, uy) = self.l0.macros_at(n);
                    out[n] = (ux * ux + uy * uy).sqrt();
                }
            }
            FieldKind::Vorticity => {
                // ω = ∂u_y/∂x − ∂u_x/∂y, центральные разности с заворотом по краям
                let mut ux = vec![0.0; nx * ny];
                let mut uy = vec![0.0; nx * ny];
                for n in 0..nx * ny {
                    let (_, a, b) = self.l0.macros_at(n);
                    ux[n] = a;
                    uy[n] = b;
                }
                for y in 0..ny {
                    for x in 0..nx {
                        let xp = (x + 1).min(nx - 1);
                        let xm = x.saturating_sub(1);
                        let yp = (y + 1).min(ny - 1);
                        let ym = y.saturating_sub(1);
                        let dvdx = (uy[y * nx + xp] - uy[y * nx + xm]) / (xp - xm).max(1) as R;
                        let dudy = (ux[yp * nx + x] - ux[ym * nx + x]) / (yp - ym).max(1) as R;
                        out[y * nx + x] = dvdx - dudy;
                    }
                }
            }
            FieldKind::Density => {
                for n in 0..nx * ny {
                    out[n] = self.l0.macros_at(n).0;
                }
            }
            FieldKind::Gamma => out.copy_from_slice(&self.l0.gamma),
        }
        (out, &self.l0.geom.solid)
    }

    /// Полное поле скорости уровня L0 — для метрик эталонных течений и радиуса влияния.
    pub fn sample_velocity(&self) -> (Vec<R>, Vec<R>) {
        let n = self.l0.nx * self.l0.ny;
        let mut ux = Vec::with_capacity(n);
        let mut uy = Vec::with_capacity(n);
        for k in 0..n {
            let (_, a, b) = self.l0.macros_at(k);
            ux.push(a);
            uy.push(b);
        }
        (ux, uy)
    }

    /// Срез вдоль осевой линии: (ρ, u_x) по каждому столбцу. Из последовательности таких
    /// срезов складывается x–t диаграмма, по которой видно, бежит возмущение со скоростью
    /// звука или конвекции и есть ли стоячие узлы.
    pub fn sample_centerline(&self) -> (Vec<R>, Vec<R>) {
        let y = self.l0.ny / 2;
        let nx = self.l0.nx;
        let mut rho = Vec::with_capacity(nx);
        let mut ux = Vec::with_capacity(nx);
        for x in 0..nx {
            let (r, a, _) = self.l0.macros_at(y * nx + x);
            rho.push(r);
            ux.push(a);
        }
        (rho, ux)
    }

    /// Есть ли в поле NaN/inf — признак развала счёта.
    pub fn is_finite(&self) -> bool {
        self.l0.f.par_iter().all(|c| c.iter().all(|v| v.is_finite()))
    }
}

// ─────────────────────────────────────────────────────────────────────────────
//  Тесты бэкенда
// ─────────────────────────────────────────────────────────────────────────────

#[cfg(test)]
mod tests {
    use super::*;

    fn bare_level(nx: usize, ny: usize) -> Level {
        Level::new(
            nx,
            ny,
            vec![0.5; nx],
            Geom {
                nx,
                ny,
                solid: vec![false; nx * ny],
                links: Vec::new(),
                wall_nodes: Vec::new(),
                body_cx: 0.0,
                body_cy: 0.0,
            },
            vec![math::feq(1.0, 0.0, 0.0); nx * ny],
        )
    }

    /// Голый периодический шаг: столкновение + перенос, без единого ГУ.
    /// Нужен, чтобы отделить ядро схемы от граничных условий.
    fn periodic_step(f: &mut [[R; Q]], tmp: &mut [[R; Q]], nx: usize, ny: usize, beta: R) {
        periodic_step_op(f, tmp, nx, ny, beta, Collision::Kbc)
    }

    fn periodic_step_op(
        f: &mut [[R; Q]],
        tmp: &mut [[R; Q]],
        nx: usize,
        ny: usize,
        beta: R,
        op: Collision,
    ) {
        for (n, c) in f.iter().enumerate() {
            tmp[n] = *c;
            match op {
                Collision::Kbc => math::collide_node(&mut tmp[n], beta, math::KbcModel::N1),
                Collision::Bgk => math::collide_node_bgk(&mut tmp[n], beta),
            };
        }
        for y in 0..ny {
            for x in 0..nx {
                for i in 0..Q {
                    let xs = (x as i32 - CX[i]).rem_euclid(nx as i32) as usize;
                    let ys = (y as i32 - CY[i]).rem_euclid(ny as i32) as usize;
                    f[y * nx + x][i] = tmp[ys * nx + xs][i];
                }
            }
        }
    }

    /// Перенос обязан сдвигать каждую популяцию ровно на её c_i, с заворотом.
    #[test]
    fn streaming_shifts_by_lattice_velocity() {
        let (nx, ny) = (7usize, 5usize);
        for i in 0..Q {
            let mut lvl = bare_level(nx, ny);
            for c in lvl.post.iter_mut() {
                *c = [0.0; Q];
            }
            lvl.post[2 * nx + 3][i] = 1.0; // дельта в (x=3, y=2)
            lvl.stream();
            let xd = (3 + CX[i]).rem_euclid(nx as i32) as usize;
            let yd = (2 + CY[i]).rem_euclid(ny as i32) as usize;
            assert!(
                (lvl.f[yd * nx + xd][i] - 1.0).abs() < 1e-15,
                "направление {i}: масса не пришла в ({xd},{yd})"
            );
            let total: R = lvl.f.iter().map(|c| c[i]).sum();
            assert!((total - 1.0).abs() < 1e-15, "направление {i}: масса не сохранилась");
        }
    }

    /// СУБСЕТОЧНОСТЬ ГРАНИЦЫ. Стенка обязана стоять не на узлах решётки: при сдвиге тела на
    /// долю клетки доли пересечения q обязаны поехать, а у разрешённого тела не должно быть
    /// ни одного отката на простой (ступенчатый) отскок. Если бы граница была ступенчатой,
    /// средняя q не реагировала бы на смещение вовсе.
    #[test]
    fn wall_is_subgrid_not_staircase() {
        let mut means = Vec::new();
        for off in [0.0, 0.25, 0.5] {
            let body = math::Body::new(math::ShapeKind::Cylinder, 60.0, 40.0 + off, 24.0, 0.0, 0.3);
            let g = Geom::build(120, 80, &Scene::single(body));
            let w = g.wall_stats();
            assert!(w.links > 100, "смещение {off}: линков всего {}", w.links);
            assert_eq!(w.simple, 0, "смещение {off}: {} ступенчатых линков", w.simple);
            assert!(w.q_min > 0.0 && w.q_max < 1.0, "смещение {off}: q вышла за (0,1)");
            means.push(w.q_mean);
        }
        let hi = means.iter().cloned().fold(R::NEG_INFINITY, R::max);
        let lo = means.iter().cloned().fold(R::INFINITY, R::min);
        assert!(hi - lo > 0.02, "средняя q не отреагировала на сдвиг тела: {means:?}");
    }

    /// Индекс граничных узлов обязан в точности разбивать список линков: без пропусков,
    /// без пересечений, в том же порядке.
    #[test]
    fn wall_node_index_partitions_links() {
        let body = math::Body::new(math::ShapeKind::Naca, 60.0, 40.0, 30.0, 12.0, 0.25);
        let g = Geom::build(120, 80, &Scene::single(body));
        let mut covered = 0usize;
        for w in &g.wall_nodes {
            assert_eq!(w.first as usize, covered, "разрыв в индексе граничных узлов");
            for l in &g.links[covered..covered + w.count as usize] {
                assert_eq!(l.node, w.node, "линк чужого узла в диапазоне");
            }
            covered += w.count as usize;
        }
        assert_eq!(covered, g.links.len(), "индекс покрыл не все линки");
    }

    /// Однородное равновесие — неподвижная точка схемы: не должно никуда уехать.
    #[test]
    fn uniform_equilibrium_is_fixed_point() {
        let (nx, ny) = (16usize, 16usize);
        let mut f = vec![math::feq(1.0, 0.03, -0.01); nx * ny];
        let mut tmp = f.clone();
        let f0 = f[0];
        for _ in 0..50 {
            periodic_step(&mut f, &mut tmp, nx, ny, 0.9);
        }
        for c in &f {
            for i in 0..Q {
                assert!((c[i] - f0[i]).abs() < 1e-12, "однородное равновесие поехало");
            }
        }
    }

    /// Затухание сдвиговой волны обязано идти с ν = c_s²(1/(2β) − ½), формула (5).
    /// Это прямая проверка того, что стабилизатор γ НЕ трогает вязкость.
    #[test]
    fn shear_wave_decays_at_prescribed_viscosity() {
        let n = 48usize;
        for &tau in &[0.6, 1.0] {
            let beta = 1.0 / (2.0 * tau);
            let nu = math::nu_of_beta(beta);
            let k = 2.0 * std::f64::consts::PI / n as R;
            let amp = 0.01;
            let mut f = vec![[0.0; Q]; n * n];
            for y in 0..n {
                for x in 0..n {
                    f[y * n + x] = math::feq(1.0, amp * (k * y as R).sin(), 0.0);
                }
            }
            let mut tmp = f.clone();
            let mode = |f: &[[R; Q]]| -> R {
                let mut s = 0.0;
                for y in 0..n {
                    for x in 0..n {
                        s += math::macros(&f[y * n + x]).1 * (k * y as R).sin();
                    }
                }
                (2.0 * s / (n * n) as R).abs()
            };
            let a0 = mode(&f);
            let steps = 1500;
            for _ in 0..steps {
                periodic_step(&mut f, &mut tmp, n, n, beta);
            }
            let a1 = mode(&f);
            let nu_measured = -(a1 / a0).ln() / (k * k * steps as R);
            let err = (nu_measured - nu).abs() / nu;
            assert!(err < 1e-2, "τ={tau}: ν измеренная {nu_measured:.6e} vs заданная {nu:.6e}");
        }
    }

    /// Диагностика: из чего складывается ошибка Тейлора–Грина. Сетка и вязкость ЗАФИКСИРОВАНЫ,
    /// меняется только амплитуда u₀. Дискретизационная ошибка от u₀ почти не зависит (задача
    /// в этом пределе линейна), а сжимаемостная идёт как Ma² ∝ u₀². Наклон и покажет, что
    /// доминирует.
    #[test]
    #[ignore = "диагностика; запуск: cargo test --release -- --ignored taylor_green_error_scaling --nocapture"]
    fn taylor_green_error_scaling() {
        let n = 128usize;
        let nu = 0.0384;
        let beta = math::beta_of_nu(nu);
        println!("  N={n}, ν={nu} (τ={:.4}) — меняем только u₀", 1.0 / (2.0 * beta));
        let mut prev: Option<(R, R)> = None;
        for &u0 in &[0.04, 0.02, 0.01, 0.005] {
            let e = taylor_green_error(n, u0, nu, beta);
            let slope = prev.map(|(p_u, p_e): (R, R)| (p_e / e.shape).ln() / (p_u / u0).ln());
            match slope {
                Some(s) => println!(
                    "  u₀={u0:<7} Ma={:.4}  форма {:.3e}  ампл {:.3e}   наклон формы по u₀: {s:.2}",
                    u0 / math::CS2.sqrt(), e.shape, e.amplitude),
                None => println!("  u₀={u0:<7} Ma={:.4}  форма {:.3e}  ампл {:.3e}",
                    u0 / math::CS2.sqrt(), e.shape, e.amplitude),
            }
            prev = Some((u0, e.shape));
        }
        println!("  наклон ≈2 ⇒ правит сжимаемость (Ma²); ≈0 ⇒ правит дискретизация");
    }

    /// Один прогон вихря Тейлора–Грина до полураспада; возвращает Σ|u_x−точн|/Σ|точн|.
    ///
    /// Начальное состояние ставится ПОЛНОСТЬЮ согласованным: скорость, давление и неравновесная
    /// часть. Давление здесь не константа — течение несёт собственное поле порядка ρu₀²,
    /// которое находится из ∇²p = 2ρ(ψ_xx·ψ_yy − ψ_xy²):
    ///
    ///   p = −(ρu₀²/4)[cos(2k₁x) + (k₁²/k₂²)·cos(2k₂y)],   δρ = p/c_s².
    ///
    /// Если стартовать с ρ ≡ 1, эта разница уходит в акустику, которая в периодическом ящике
    /// почти не затухает и садится полкой на ошибку скорости, ломая порядок сходимости.
    fn taylor_green_error(n: usize, u0: R, nu: R, beta: R) -> TgError {
        taylor_green_error_op(n, u0, nu, beta, Collision::Kbc)
    }

    fn taylor_green_error_op(n: usize, u0: R, nu: R, beta: R, op: Collision) -> TgError {
        let (k1, k2) = (1.0, 4.0);
        let tau = 1.0 / (2.0 * beta);
        let kk1 = 2.0 * std::f64::consts::PI * k1 / n as R;
        let kk2 = 2.0 * std::f64::consts::PI * k2 / n as R;
        let decay = nu * (kk1 * kk1 + kk2 * kk2);
        let tc = ((2.0_f64).ln() / decay).round() as usize;
        let exact = |x: usize, y: usize, e: R| -> (R, R) {
            let (a, b) = (kk1 * x as R, kk2 * y as R);
            (-u0 * a.cos() * b.sin() * e, (k1 / k2) * u0 * a.sin() * b.cos() * e)
        };
        let mut f = vec![[0.0; Q]; n * n];
        for y in 0..n {
            for x in 0..n {
                let (ux, uy) = exact(x, y, 1.0);
                let (a, b) = (kk1 * x as R, kk2 * y as R);
                let p = -(u0 * u0 / 4.0)
                    * ((2.0 * a).cos() + (k1 * k1) / (k2 * k2) * (2.0 * b).cos());
                let rho = 1.0 + p / math::CS2;
                let dxux = u0 * kk1 * a.sin() * b.sin();
                let dyux = -u0 * kk2 * a.cos() * b.cos();
                let dxuy = (k1 / k2) * u0 * kk1 * a.cos() * b.cos();
                let dyuy = -(k1 / k2) * u0 * kk2 * a.sin() * b.sin();
                // Π = ρc_s²δ + ρuu + Π⁽¹⁾, Π⁽¹⁾ = −τρc_s²(∂_αu_β + ∂_βu_α) — ур. (53)
                let pref = -tau * rho * math::CS2;
                let pxx = rho * math::CS2 + rho * ux * ux + pref * 2.0 * dxux;
                let pyy = rho * math::CS2 + rho * uy * uy + pref * 2.0 * dyuy;
                let pxy = rho * ux * uy + pref * (dyux + dxuy);
                f[y * n + x] = math::grad_init(rho, ux, uy, pxx, pxy, pyy);
            }
        }
        let mut tmp = f.clone();
        for _ in 0..tc {
            periodic_step_op(&mut f, &mut tmp, n, n, beta, op);
        }
        let e = (-decay * tc as R).exp();
        let (mut num, mut den) = (0.0, 0.0);
        // заодно раскладываем ошибку: наилучшая подгонка амплитуды к точной форме
        let (mut dot, mut nrm2) = (0.0, 0.0);
        for y in 0..n {
            for x in 0..n {
                let got = math::macros(&f[y * n + x]).1;
                let want = exact(x, y, e).0;
                num += (got - want).abs();
                den += want.abs();
                dot += got * want;
                nrm2 += want * want;
            }
        }
        let amp = dot / nrm2; // 1.0 = амплитуда совпала
        let (mut snum, mut sden) = (0.0, 0.0);
        for y in 0..n {
            for x in 0..n {
                let got = math::macros(&f[y * n + x]).1;
                let want = amp * exact(x, y, e).0;
                snum += (got - want).abs();
                sden += want.abs();
            }
        }
        TgError { total: num / den, amplitude: (amp - 1.0).abs(), shape: snum / sden }
    }

    /// Разложение ошибки эталона: полная, вклад амплитуды (скорость затухания) и вклад формы.
    struct TgError {
        total: R,
        amplitude: R,
        shape: R,
    }

    /// ВИХРЬ ТЕЙЛОРА–ГРИНА — первый эталон 2D-статьи (разд. VI). Единственное из трёх течений
    /// статьи, у которого есть ТОЧНОЕ аналитическое решение, поэтому проверяется не «похоже на
    /// чужой прогон», а прямое совпадение с формулой и заявленный статьёй ВТОРОЙ ПОРЯДОК
    /// сходимости (рис. 1a: Re = 100, u₀ = 0.03, N ∈ {64, 128, 256}).
    ///
    /// Постановка ровно по статье:
    ///   u = ∇×[(u₀/k₂)cos(k₁x)cos(k₂y)·exp(−ν(k₁²+k₂²)t)],  k₁ = 1, k₂ = 4,
    ///   область 0 < x,y < 2π на сетке N×N, Re = u₀N/ν, полураспад t_c = ln2/[ν(k₁²+k₂²)].
    /// Волновые числа переводятся в решёточные: K = 2πk/N (узел — единица длины).
    /// Старт — приближением Града (ур. 58), как в статье.
    /// Метрика — та же, что на рис. 1: Σ|u_x − u_x^точн| / Σ|u_x^точн| в момент t_c.
    #[test]
    #[ignore = "долгий (до 256², t_c растёт как N²); запуск: cargo test --release -- --ignored"]
    fn taylor_green_converges_at_second_order() {
        // ДИФФУЗИОННОЕ ИЗМЕЛЬЧЕНИЕ: ν фиксирована, u₀ ∝ 1/N. Тогда Re = u₀N/ν сохраняется
        // (все три сетки считают ОДНО И ТО ЖЕ течение), а число Маха падает как 1/N — вместе
        // с ним падает и сжимаемостная ошибка метода. Только при таком измельчении второй
        // порядок виден целиком; при фиксированном u₀ ошибка упирается в полку O(Ma²), которая
        // от сетки не зависит вовсе (см. taylor_green_error_scaling).
        let nu = 0.0192;
        let beta = math::beta_of_nu(nu);
        let mut errs = Vec::new();
        for &n in &[64usize, 128, 256] {
            let u0 = 0.03 * 64.0 / n as R;
            let e = taylor_green_error(n, u0, nu, beta);
            println!(
                "  N={n:>4} u₀={u0:.5} Re={:.0}  полная {:.3e}  амплитуда {:.3e}  форма {:.3e}",
                u0 * n as R / nu, e.total, e.amplitude, e.shape
            );
            errs.push((n, e));
        }
        for w in errs.windows(2) {
            let (n0, e0) = (&w[0].0, &w[0].1);
            let (n1, e1) = (&w[1].0, &w[1].1);
            let ord = |a: R, b: R| (a / b).ln() / (*n1 as R / *n0 as R).ln();
            let (pt, pa, ps) = (
                ord(e0.total, e1.total),
                ord(e0.amplitude, e1.amplitude),
                ord(e0.shape, e1.shape),
            );
            println!("  порядок {n0}→{n1}: полная {pt:.2}  амплитуда {pa:.2}  форма {ps:.2}");
            assert!(
                (1.7..2.6).contains(&pt),
                "порядок {n0}→{n1} = {pt:.2}, статья (разд. VI) заявляет второй"
            );
        }
    }

    /// Диагностика: тот же вихрь Тейлора–Грина оператором LBGK. Статья (рис. 1) утверждает,
    /// что на этом течении LBGK и все варианты KBC идут практически одинаково; если полка по
    /// форме есть у обоих — она от схемы D2Q9, а не от энтропийного стабилизатора.
    #[test]
    #[ignore = "диагностика; запуск: cargo test --release -- --ignored taylor_green_kbc_vs_bgk --nocapture"]
    fn taylor_green_kbc_vs_bgk() {
        let (u0, re) = (0.03, 100.0);
        for &n in &[64usize, 128, 256] {
            let nu = u0 * n as R / re;
            let beta = math::beta_of_nu(nu);
            let k = taylor_green_error_op(n, u0, nu, beta, Collision::Kbc);
            let b = taylor_green_error_op(n, u0, nu, beta, Collision::Bgk);
            println!(
                "  N={n:>4}  KBC: форма {:.3e} ампл {:.3e}   LBGK: форма {:.3e} ампл {:.3e}",
                k.shape, k.amplitude, b.shape, b.amplitude
            );
        }
    }

    /// СКВОЗНАЯ СВЕРКА С ЭТАЛОНОМ. Дважды периодический сдвиговый слой — один из трёх
    /// бенчмарков 2D-статьи (разд. VIII). Постановка: N=128, Re=30000, u0=0.04, κ=80, δ=0.05,
    /// одно конвективное время t_c = N/u0 шагов.
    ///
    /// Эталон 0.6035 — отношение энстрофии к начальной на fp64-версии python-решателя
    /// (docs/theory/solver_2x_sdf) после исправления порога вырожденности γ. Тот же прогон на
    /// испорченном (абсолютном) пороге давал 0.6599, то есть +9.3%, а чистый LBGK при этих
    /// параметрах разваливается. Тест ловит и «схема молча выродилась в LBGK», и «схема
    /// считает не то».
    #[test]
    #[ignore = "долгий (3200 шагов на 128²); запуск: cargo test --release -- --ignored"]
    fn doubly_periodic_shear_layer_matches_reference() {
        let n = 128usize;
        let (u0, kappa, delta, re) = (0.04, 80.0, 0.05, 30000.0);
        let nu = u0 * n as R / re;
        let beta = math::beta_of_nu(nu);
        let tc = (n as R / u0) as usize;

        let mut f = vec![[0.0; Q]; n * n];
        for y in 0..n {
            let yy = y as R / n as R;
            let ux = if y as R <= n as R / 2.0 {
                u0 * (kappa * (yy - 0.25)).tanh()
            } else {
                u0 * (kappa * (0.75 - yy)).tanh()
            };
            for x in 0..n {
                let uy = delta * u0 * (2.0 * std::f64::consts::PI * (x as R / n as R + 0.25)).sin();
                f[y * n + x] = math::feq(1.0, ux, uy);
            }
        }
        let mut tmp = f.clone();

        let enstrophy = |f: &[[R; Q]]| -> R {
            let mut ux = vec![0.0; n * n];
            let mut uy = vec![0.0; n * n];
            for k in 0..n * n {
                let (_, a, b) = math::macros(&f[k]);
                ux[k] = a;
                uy[k] = b;
            }
            let mut s = 0.0;
            for y in 0..n {
                for x in 0..n {
                    let (xp, xm) = ((x + 1) % n, (x + n - 1) % n);
                    let (yp, ym) = ((y + 1) % n, (y + n - 1) % n);
                    let w = 0.5 * (uy[y * n + xp] - uy[y * n + xm])
                        - 0.5 * (ux[yp * n + x] - ux[ym * n + x]);
                    s += w * w;
                }
            }
            s / (n * n) as R
        };

        let e0 = enstrophy(&f);
        for _ in 0..=tc {
            periodic_step(&mut f, &mut tmp, n, n, beta);
        }
        assert!(f.iter().all(|c| c.iter().all(|v| v.is_finite())), "счёт развалился");
        let ratio = enstrophy(&f) / e0;
        assert!(
            (ratio - 0.6035).abs() < 0.012,
            "энстрофия/начальная = {ratio:.4}, эталон fp64 python-решателя 0.6035 \
             (испорченный порог γ давал 0.6599)"
        );
    }
}