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

use rayon::prelude::*;

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

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

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

impl Geom {
    /// Собрать геометрию уровня по телу. `solid` определяется знаком SDF: φ ≤ 0 — тело.
    pub fn build(nx: usize, ny: usize, body: &Body) -> Self {
        let n = nx * ny;
        let mut phi = vec![0.0; n];
        let mut solid = vec![false; n];
        for y in 0..ny {
            for x in 0..nx {
                let p = body.sdf(x as R, y as R);
                phi[y * nx + x] = p;
                solid[y * nx + x] = p <= 0.0;
            }
        }
        let links = build_links(nx, ny, &solid, &phi);
        Geom { nx, ny, solid, links, body_cx: body.cx, body_cy: body.cy }
    }
}

/// Индекс соседа по направлению 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]) -> 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: true,
                });
            }
        }
    }
    links
}

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

/// Один уровень: поля популяций до/после столкновения, β и геометрия.
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 {
    fn new(nx: usize, ny: usize, beta: Vec<R>, geom: Geom) -> Self {
        let f0 = math::feq(1.0, 0.0, 0.0);
        let n = nx * ny;
        Level { nx, ny, f: vec![f0; n], post: vec![f0; n], gamma: vec![2.0; n], beta, geom }
    }

    /// Столкновение по всем жидким узлам. Внутри тела не считаем: эти популяции фиктивны
    /// (ГУ Bouzidi перезаписывает всё, что могло бы прийти из тела в жидкость), а счёт там
    /// только жжёт такты и способен родить NaN при экстремальных режимах.
    fn collide(&mut self, op: Collision) -> 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),
                    Collision::Bgk => math::collide_node_bgk(p, b),
                };
                *g = k.gamma;
                Stats::of(k, b)
            })
            .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;
    }

    /// Интерполированный отскок 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;
        }
    }

    /// 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, суммой по линкам тела.
    fn force(&self) -> (R, R, R) {
        let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0);
        let nx = self.nx;
        for l in &self.geom.links {
            if !l.body {
                continue;
            }
            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);
            fx += dfx;
            fy += dfy;
            // плечо — до ТОЧКИ ПЕРЕСЕЧЕНИЯ линка со стенкой 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;
            tz += rx * dfy - ry * dfx;
        }
        (fx, fy, tz)
    }

    #[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) -> 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) < 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.body);

        // губка: плавный рост вязкости в последних sponge_len столбцах. Канал
        // «скорость-вход + давление-выход» — недодемпфированный акустический резонатор,
        // губка гасит и вихри, и акустику до прихода на выход.
        let beta0: Vec<R> = (0..nx)
            .map(|x| {
                if spec.sponge_len == 0 {
                    return spec.beta0;
                }
                let start = nx - 1 - spec.sponge_len;
                let s = math::smoothstep((x as R - start as R) / spec.sponge_len as R);
                let nu = math::nu_of_beta(spec.beta0);
                math::beta_of_nu(nu * (1.0 + (spec.sponge_mult - 1.0) * s))
            })
            .collect();

        let l0 = Level::new(nx, ny, beta0, geom0);
        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 body1 = spec.body.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, &body1);
            // τ_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);
            (Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1)), 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 u = sp.units.u_lat * math::smoothstep(t as R / sp.ramp.max(1) as R);
        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 stats0 = self.l0.collide(sp_collision);
        self.l0.stream();
        self.l0.apply_bouzidi();
        self.l0.free_slip_walls();
        self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap);

        // ── уровень 1: r подшагов с временной интерполяцией рамки ──
        let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0);
        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);
                l1.stream();
                l1.apply_bouzidi();
                // силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на
                // последнем подшаге даёт лишний шум в рядах при том же среднем
                let (a, b, c) = l1.force();
                fx += a;
                fy += b;
                tz += 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;
            fx *= inv;
            fy *= inv;
            tz *= inv;
            p.restrict_to(&l1.f, &mut self.l0.f);
        } else {
            let (a, b, c) = self.l0.force();
            fx = a;
            fy = b;
            tz = c;
        }

        // ── диагностика ──
        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;
        StepRec {
            step: t,
            fx,
            fy,
            tz,
            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.body.d * p.r as R,
            None => self.spec.body.d,
        }
    }

    /// Поле для картинки: (значение, маска тела) на сетке 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)
    }

    /// Есть ли в поле 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(),
                body_cx: 0.0,
                body_cy: 0.0,
            },
        )
    }

    /// Голый периодический шаг: столкновение + перенос, без единого ГУ.
    /// Нужен, чтобы отделить ядро схемы от граничных условий.
    fn periodic_step(f: &mut [[R; Q]], tmp: &mut [[R; Q]], nx: usize, ny: usize, beta: R) {
        for (n, c) in f.iter().enumerate() {
            tmp[n] = *c;
            math::collide_node(&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}: масса не сохранилась");
        }
    }

    /// Однородное равновесие — неподвижная точка схемы: не должно никуда уехать.
    #[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}");
        }
    }

    /// СКВОЗНАЯ СВЕРКА С ЭТАЛОНОМ. Дважды периодический сдвиговый слой — один из трёх
    /// бенчмарков 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)"
        );
    }
}