//! ОСНОВНОЙ ФАЙЛ: разбор параметров, сборка задачи, цикл по шагам, вывод. //! //! Здесь же живёт контракт между оркестратором и бэкендами: `Spec` (что считать), //! `StepRec` (что бэкенд отдаёт за шаг) и `FieldKind` (что рисовать). Бэкенды (cpu.rs, //! gpu.rs) реализуют один и тот же интерфейс и взаимозаменяемы флагом `--backend`. //! //! Раскладка задачи — канал с обтекаемым телом: вход слева (Zou–He по скорости), выход //! справа (Zou–He по давлению), сверху и снизу зеркальные (free-slip) стенки, тело — //! no-slip через SDF + интерполированный отскок Bouzidi. mod cpu; mod gif; mod math; #[cfg(feature = "gpu")] mod gpu; use std::io::Write as _; use std::time::Instant; use clap::Parser; use math::{Body, Scene, ShapeKind, Units, MAX_BODY_BUCKETS, R}; // ───────────────────────────────────────────────────────────────────────────── // Контракт оркестратор ↔ бэкенд // ───────────────────────────────────────────────────────────────────────────── #[derive(Clone, Copy, Debug, PartialEq, Eq)] pub enum Collision { /// Энтропийный KBC-D (он же N1) — основной оператор модели. Kbc, /// Классический LBGK — для сравнительных прогонов; это ровно KBC при γ ≡ 2. Bgk, } /// Постановка задачи. Канал — рабочая; остальные три периодичны по обеим осям и не имеют /// граничных условий вовсе: это эталонные течения из статей, на которых схема проверяется /// в чистом виде, без вклада стенок. #[derive(Clone, Copy, Debug, PartialEq, Eq)] pub enum Case { /// Обтекание тела в канале: вход Zou–He по скорости, выход по давлению, стенки зеркальные. Channel, /// Вихрь Тейлора–Грина (разд. VI 2D-статьи) — единственное течение с ТОЧНЫМ решением. TaylorGreen, /// Дважды периодический сдвиговый слой (разд. VII). ShearLayer, /// Затухающая двумерная турбулентность. DecayingTurbulence, } impl Case { pub fn from_str(v: &str) -> Option { Some(match v { "channel" => Case::Channel, "taylor-green" => Case::TaylorGreen, "shear-layer" => Case::ShearLayer, "decaying-turbulence" => Case::DecayingTurbulence, _ => return None, }) } pub const ALL: [&'static str; 4] = ["channel", "taylor-green", "shear-layer", "decaying-turbulence"]; pub fn is_periodic(&self) -> bool { !matches!(self, Case::Channel) } } /// Что рисовать в анимации. #[derive(Clone, Copy, Debug, PartialEq, Eq)] pub enum FieldKind { Speed, Vorticity, Density, Gamma, } impl FieldKind { pub fn label(&self) -> &'static str { match self { FieldKind::Speed => "SPEED", FieldKind::Vorticity => "VORTICITY", FieldKind::Density => "DENSITY", FieldKind::Gamma => "GAMMA", } } } /// Запись за один шаг — сырьё для всего отчёта. #[derive(Clone, Copy, Debug, Default)] pub struct StepRec { pub step: u64, /// Суммарная сила и момент по всей сцене. pub fx: R, pub fy: R, pub tz: R, /// Она же по телам: `body[b] = (F_x, F_y, T_z)`. Тела за пределом набора вёдер /// сваливаются в последнее, поэтому сумма по вёдрам всегда точна. pub body: [[R; 3]; MAX_BODY_BUCKETS], pub uy_probe: R, pub rho_mean: R, pub max_u: R, pub gamma_mean: R, pub gamma_min: R, pub gamma_max: R, pub degenerate_frac: R, pub xi_negative_frac: R, } /// Полностью разрешённая постановка задачи в решёточных единицах. #[derive(Clone)] pub struct Spec { pub nx: usize, pub ny: usize, pub scene: Scene, pub refine: usize, /// Границы патча измельчения в координатах L0: (ax, bx, ay, by). pub patch: Option<(usize, usize, usize, usize)>, pub units: Units, pub beta0: R, pub steps: u64, pub ramp: u64, /// Направление набегающего потока, рад (0 — вдоль +x). pub flow_angle: R, pub pert_amp: R, pub pert_dur: u64, pub outlet_extrapolate: bool, /// Начальное поле уже несёт набегающий поток (иначе — покой с разгоном входа). pub init_uniform: bool, /// Ширина сведения стартовой скорости к нулю у тела, клеток. pub init_taper: R, pub sponge_len: usize, /// Длина губки после входа, столбцов. pub sponge_in: usize, pub sponge_mult: R, pub collision: Collision, /// Что входит в сдвиговую часть s (см. `math::KbcModel`). pub kbc_model: math::KbcModel, /// Чем замыкаются недостающие популяции на теле (см. `math::WallModel`). pub wall: math::WallModel, /// Постановка задачи. pub case: Case, /// Узел зонда следа в координатах L0. pub probe: (usize, usize), } // ───────────────────────────────────────────────────────────────────────────── // Параметры запуска // ───────────────────────────────────────────────────────────────────────────── #[derive(Parser, Debug)] #[command( name = "kbc2d", about = "Двумерный решатель LBM D2Q9 с энтропийным столкновением KBC (модель D/N1)", long_about = "Обтекание тела в канале: KBC-D по Bösch/Chikatamarla/Karlin, SDF+Bouzidi на \ теле, Zou–He вход/выход, вложенный патч измельчения. Анимация синхронизируется с ФИЗИЧЕСКИМ \ временем потока, а не со скоростью счёта." )] struct Cli { // ── физика потока ── /// Скорость набегающего потока, м/с #[arg(long, default_value_t = 30.0, help_heading = "Физика")] u_phys: R, /// Размер ячейки (плотность сетки), м #[arg(long, default_value_t = 0.1, help_heading = "Физика")] dx: R, /// Решёточная скорость потока (ячеек за шаг). Задаёт δt = u_lat·dx/u_phys. /// 0.05 ⇒ Ma≈0.09. Значение 1.0 воспроизводит арифметику «30 м/с ÷ 0.1 м = 300 шагов/с», /// но это Ma≈1.73 — за пределами применимости LBM. #[arg(long, default_value_t = 0.05, help_heading = "Физика")] u_lat: R, /// Число Рейнольдса по характерному размеру тела #[arg(long, default_value_t = 150.0, help_heading = "Физика")] re: R, /// Направление потока, градусы (0 — вдоль канала) #[arg(long, default_value_t = 0.0, help_heading = "Физика")] flow_angle: R, /// Постановка. channel — обтекание тела; остальные три периодичны по обеим осям, тела и /// граничных условий не имеют и служат эталонами из статей: taylor-green (разд. VI, есть /// точное решение), shear-layer (разд. VII), decaying-turbulence. #[arg(long, default_value = "channel", value_parser = Case::ALL, help_heading = "Физика")] case: String, /// Период съёма метрик эталонных течений (энергия, энстрофия, палинстрофия), шагов. /// Каждый съём тянет поле с устройства, поэтому редкий. #[arg(long, default_value_t = 200, help_heading = "Физика")] case_every: u64, // ── сетка ── /// Размер домена по x, ячеек #[arg(long, default_value_t = 480, help_heading = "Сетка")] nx: usize, /// Размер домена по y, ячеек #[arg(long, default_value_t = 240, help_heading = "Сетка")] ny: usize, /// Коэффициент измельчения вложенного патча (1 — патча нет) #[arg(long, default_value_t = 2, help_heading = "Сетка")] refine: usize, /// Границы патча в координатах L0: ax,bx,ay,by (по умолчанию строятся вокруг тела) #[arg(long, help_heading = "Сетка")] patch: Option, // ── тело ── /// Форма обтекаемого тела #[arg(long, default_value = "cylinder", value_parser = ShapeKind::ALL, help_heading = "Тело")] shape: String, /// Характерный размер тела (диаметр/сторона/хорда), ячеек #[arg(long, default_value_t = 24.0, help_heading = "Тело")] size: R, /// Относительная толщина для ellipse / naca / plate #[arg(long, default_value_t = 0.3, help_heading = "Тело")] thickness: R, /// Угол атаки тела, градусы #[arg(long, default_value_t = 0.0, help_heading = "Тело")] body_angle: R, /// Положение центра тела по x, ячеек (по умолчанию nx/4) #[arg(long, help_heading = "Тело")] body_x: Option, /// Положение центра тела по y, ячеек (по умолчанию ny/2) #[arg(long, help_heading = "Тело")] body_y: Option, /// Четырёхзначное обозначение профиля NACA, например 4412 (кривизна 4%, её максимум на /// 40% хорды, толщина 12%). Действует при --shape naca и перекрывает --thickness. #[arg(long, help_heading = "Тело")] naca: Option, /// Вершины произвольного многоугольника при --shape polygon: "x1,y1;x2,y2;..." в клетках, /// отсчёт от центра тела. Обход любой, SDF точный. #[arg(long, allow_hyphen_values = true, help_heading = "Тело")] poly: Option, /// НЕСКОЛЬКО ТЕЛ вместо одного. Список через точку с запятой, каждое — форма и параметры: /// "cylinder:d=24,x=120,y=120; cylinder:d=24,x=200,y=120". Ключи: d — размер, x/y — центр, /// a — угол атаки (град), t — относительная толщина, naca — код профиля, /// p — вершины многоугольника через |. Перекрывает --shape и прочие одиночные ключи. #[arg(long, allow_hyphen_values = true, help_heading = "Тело")] bodies: Option, // ── время ── /// Число шагов симуляции. Взаимоисключающе с --time. #[arg(long, conflicts_with = "time", help_heading = "Время")] steps: Option, /// Длительность прогона в СЕКУНДАХ физического времени. Число шагов считается как /// time/δt, где δt = u_lat·dx/u_phys — то есть зависит и от размера клетки, и от /// решёточной скорости. Взаимоисключающе с --steps. #[arg(long, help_heading = "Время")] time: Option, /// Длина smoothstep-разгона входа, шагов #[arg(long, default_value_t = 1000, help_heading = "Время")] ramp: u64, /// Амплитуда стартового поперечного импульса (доля от U); 0 — выключить #[arg(long, default_value_t = 0.1, help_heading = "Время")] pert_amp: R, /// Длительность стартового импульса, шагов #[arg(long, default_value_t = 300, help_heading = "Время")] pert_dur: u64, /// Начальное поле: uniform — домен сразу заполнен набегающим потоком и вход включён на /// полную (умолчание); rest — покой, вход разгоняется за --ramp шагов. Старт из покоя /// закачивает в канал продольную акустическую моду, которая почти не затухает. #[arg(long, default_value = "uniform", value_parser = ["uniform", "rest"], help_heading = "Время")] init: String, /// Ширина полосы у тела, на которой стартовая скорость сводится к нулю, клеток. /// /// ЗАМЕРЕНО, что это НЕ ПОМОГАЕТ, и ключ оставлен только для повторной проверки. /// Однородный старт заливает потоком место, где стоит тело, и на первом шаге у стенки /// рождается возмущение плотности ~2.5·10⁻². Сглаживание монотонно давит именно его /// (при ширине 24 клетки — до 3.1·10⁻³, восьмикратно), но ПИК ЗА ПРОГОН при этом даже /// подрастает: 2.54·10⁻² → 3.22·10⁻². Возмущение просто переносится во времени — поток /// всё равно обязан разогнаться вокруг тела, и энергия этого переходного процесса задана /// физикой, а не гладкостью начального поля. Продольную моду ест губка, а не это. #[arg(long, default_value_t = 0.0, help_heading = "Время")] init_taper: R, // ── численная схема ── /// Оператор столкновения #[arg(long, default_value = "kbc", value_parser = ["kbc", "bgk"], help_heading = "Схема")] collision: String, /// Модель стенки на теле. hrr (умолчание) — восстановление недостающих популяций по /// целевым моментам с рекурсивной регуляризацией до 3-го порядка Эрмита; grad — то же, /// но с обрывом на тензоре давлений; bouzidi — интерполированный отскок по каждому линку /// со своей долей пересечения q; staircase — простой отскок, q игнорируется (не для /// счёта, а как база сравнения). Замерено: по разрешению геометрии bouzidi точнее /// моментных схем, зато у тех естественная форма для подвижных стенок — см. README. #[arg(long, default_value = "hrr", value_parser = math::WallModel::ALL, help_heading = "Схема")] wall: String, /// Состав сдвиговой части KBC (табл. I 2D-статьи): n1 — только девиатор {N, Π_xy} /// (KBC D), объёмная вязкость гуляет вместе с γ и может стать отрицательной; n2 — /// девиатор со следом {N, Π_xy, T} (KBC C), объёмная вязкость фиксирована ξ = ν. #[arg(long, default_value = "n1", value_parser = math::KbcModel::ALL, help_heading = "Схема")] kbc_model: String, /// Поперечная скорость на выходе: extrapolate — нуль-градиент (умолчание: в факторном /// исследовании эталонного решателя это лучший вариант по всем метрикам), zero — жёсткий /// ноль (классический Zou–He; отражает вихри дорожки назад к телу и завышает rms Cl) #[arg(long, default_value = "extrapolate", value_parser = ["zero", "extrapolate"], help_heading = "Схема")] outlet: String, /// Длина поглощающей губки перед выходом, столбцов. Не задана — подбирается по домену; /// 0 — выключить. Губка гасит продольную моду, которую вязкость сама не гасит. #[arg(long, help_heading = "Схема")] sponge_len: Option, /// Длина губки ПОСЛЕ входа, столбцов (0 — выключена). Вход по скорости отражает продольные /// волны как жёсткий поршень; губка у входа гасит их до отражения. #[arg(long, default_value_t = 0, help_heading = "Схема")] sponge_in: usize, /// Во сколько раз губка поднимает вязкость #[arg(long, default_value_t = 30.0, help_heading = "Схема")] sponge_mult: R, /// Бэкенд вычислений #[arg(long, default_value = "cpu", value_parser = ["cpu", "gpu"], help_heading = "Схема")] backend: String, /// Число потоков CPU (0 — по числу ядер) #[arg(long, default_value_t = 0, help_heading = "Схема")] threads: usize, // ── анимация ── /// Файл гифки (не указан — анимация не пишется) #[arg(long, help_heading = "Анимация")] gif: Option, /// Шагов на кадр: число либо auto (подобрать под --gif-fps в реальном времени) #[arg(long, default_value = "auto", help_heading = "Анимация")] gif_every: String, /// Целевая частота кадров для режима auto #[arg(long, default_value_t = 30.0, help_heading = "Анимация")] gif_fps: R, /// Скорость воспроизведения: 1 — реальное время, 0.1 — замедление в 10 раз #[arg(long, default_value_t = 1.0, help_heading = "Анимация")] gif_speed: R, /// Что рисовать #[arg(long, default_value = "vorticity", value_parser = ["speed", "vorticity", "density", "gamma"], help_heading = "Анимация")] gif_field: String, /// Палитра #[arg(long, default_value = "coolwarm", value_parser = gif::ColorMap::ALL, help_heading = "Анимация")] gif_cmap: String, /// Целочисленное увеличение картинки #[arg(long, default_value_t = 2, help_heading = "Анимация")] gif_scale: usize, /// Усреднять блок k×k клеток в один пиксель перед отрисовкой. Обязательно на крупных /// сетках: кадр с 4096×2048 иначе неподъёмен. #[arg(long, default_value_t = 1, help_heading = "Анимация")] gif_downsample: usize, /// Диапазон нормировки цвета: lo,hi (по умолчанию — по скорости потока) #[arg(long, help_heading = "Анимация")] gif_range: Option, /// Не рисовать служебную надпись на кадре #[arg(long, help_heading = "Анимация")] no_hud: bool, // ── вывод ── /// Период строк живого отчёта, шагов (0 — выключить) #[arg(long, default_value_t = 1000, help_heading = "Вывод")] report_every: u64, /// Подробность: quiet — только итог, normal — живой отчёт, full — плюс разбор KBC #[arg(long, default_value = "normal", value_parser = ["quiet", "normal", "full"], help_heading = "Вывод")] verbose: String, /// Число окон в отчёте о сходимости #[arg(long, default_value_t = 10, help_heading = "Вывод")] windows: usize, /// Прореживание временных рядов: хранить каждую N-ю запись. Статистика и спектр от этого /// не страдают, пока N много меньше периода схода вихрей, зато сверхдлинные прогоны /// перестают съедать память (на 5·10⁶ шагов полный ряд — сотни мегабайт). #[arg(long, default_value_t = 1, help_heading = "Вывод")] series_every: u64, /// Выгрузить временные ряды в CSV #[arg(long, help_heading = "Вывод")] csv: Option, /// Файл x–t диаграммы: срез ⟨ρ⟩ и u_x вдоль осевой линии, по строке на срез. По наклону /// полос видно, бежит возмущение со скоростью звука или конвекции, а стоячие узлы /// проявляются сразу. Основной инструмент разбора стартовой акустики. #[arg(long, help_heading = "Вывод")] xt: Option, /// Период записи срезов x–t, шагов #[arg(long, default_value_t = 200, help_heading = "Вывод")] xt_every: u64, /// Машиночитаемая сводка прогона в JSON: все ключевые метрики одним файлом. Без неё /// разбор кампании из десятков прогонов пришлось бы вести глазами. #[arg(long, help_heading = "Вывод")] summary: Option, /// Ряд метрик эталонного течения в CSV (шаг, время, энергия, энстрофия, палинстрофия). #[arg(long, help_heading = "Вывод")] case_csv: Option, } // ───────────────────────────────────────────────────────────────────────────── // Сборка постановки // ───────────────────────────────────────────────────────────────────────────── fn parse_pair(s: &str, what: &str) -> Result<(R, R), String> { let p: Vec<&str> = s.split(',').collect(); if p.len() != 2 { return Err(format!("{what}: ожидалось два числа через запятую, получено «{s}»")); } let a = p[0].trim().parse::().map_err(|e| format!("{what}: {e}"))?; let b = p[1].trim().parse::().map_err(|e| format!("{what}: {e}"))?; Ok((a, b)) } /// Разбор одного тела из мини-синтаксиса `форма:ключ=значение,...`. fn parse_body(spec: &str, def_x: R, def_y: R) -> Result { let (name, rest) = match spec.split_once(':') { Some((a, b)) => (a.trim(), b), None => (spec.trim(), ""), }; let kind = ShapeKind::from_str(name).ok_or(format!("неизвестная форма «{name}»"))?; let (mut d, mut x, mut y, mut a, mut t) = (24.0, def_x, def_y, 0.0, 0.3); let (mut camber, mut cpos) = (0.0, 0.0); let mut verts: Vec<[R; 2]> = Vec::new(); for kv in rest.split(',').filter(|z| !z.trim().is_empty()) { let (k, v) = kv .split_once('=') .ok_or(format!("в «{}» ожидалось ключ=значение", kv.trim()))?; let (k, v) = (k.trim(), v.trim()); let num = |v: &str| v.parse::().map_err(|e| format!("{k}: {e}")); match k { "d" | "size" | "c" => d = num(v)?, "x" => x = num(v)?, "y" => y = num(v)?, "a" | "angle" => a = num(v)?, "t" | "thickness" => t = num(v)?, "naca" => { let (m, pp, tt) = Body::naca_code(v).ok_or(format!("naca: ожидалось четыре цифры, «{v}»"))?; camber = m; cpos = pp; t = tt; } "p" | "poly" => { let n: Vec<&str> = v.split('|').collect(); if n.len() < 6 || n.len() % 2 != 0 { return Err("p: нужно чётное число координат, минимум три вершины".into()); } for pair in n.chunks(2) { verts.push([num(pair[0])?, num(pair[1])?]); } } _ => return Err(format!("неизвестный ключ «{k}»")), } } if kind == ShapeKind::Polygon && verts.len() < 3 { return Err("для polygon нужен ключ p с вершинами".into()); } Ok(Body::full(kind, x, y, d, a, t, camber, cpos, &verts)) } fn parse_poly(v: &str) -> Result, String> { let mut out = Vec::new(); for pt in v.split(';').filter(|z| !z.trim().is_empty()) { let (a, b) = parse_pair(pt, "--poly")?; out.push([a, b]); } if out.len() < 3 { return Err("--poly: нужно минимум три вершины".into()); } Ok(out) } fn build_spec(cli: &Cli) -> Result { let case = Case::from_str(&cli.case).ok_or("неизвестная постановка")?; if cli.nx < 16 || cli.ny < 16 { return Err("сетка меньше 16×16 не имеет смысла".into()); } if cli.u_lat <= 0.0 || cli.u_lat >= 1.0 { return Err(format!("--u-lat обязан лежать в (0,1), получено {}", cli.u_lat)); } if cli.u_phys <= 0.0 || cli.dx <= 0.0 || cli.re <= 0.0 { return Err("--u-phys, --dx и --re обязаны быть положительными".into()); } if cli.refine == 0 { return Err("--refine обязан быть ≥ 1".into()); } let cx = cli.body_x.unwrap_or(cli.nx as R / 4.0); let cy = cli.body_y.unwrap_or(cli.ny as R / 2.0); // В эталонных течениях тела нет вовсе: сцена пуста, маска твёрдого пуста, ГУ не // накладываются. Иначе цилиндр по умолчанию стоял бы прямо посреди эталона. let scene = if case.is_periodic() { Scene { bodies: Vec::new() } } else { match &cli.bodies { Some(list) => { let mut v = Vec::new(); for one in list.split(';').filter(|z| !z.trim().is_empty()) { v.push(parse_body(one, cx, cy)?); } if v.is_empty() { return Err("--bodies: список пуст".into()); } Scene { bodies: v } } None => { let kind = ShapeKind::from_str(&cli.shape).ok_or("неизвестная форма")?; let (mut camber, mut cpos, mut thick) = (0.0, 0.0, cli.thickness); if let Some(code) = &cli.naca { let (m, pp, tt) = Body::naca_code(code) .ok_or("--naca: ожидалось четыре цифры, например 4412")?; camber = m; cpos = pp; thick = tt; } let verts = match &cli.poly { Some(v) => parse_poly(v)?, None => Vec::new(), }; if kind == ShapeKind::Polygon && verts.is_empty() { return Err("--shape polygon требует --poly с вершинами".into()); } Scene::single(Body::full( kind, cx, cy, cli.size, cli.body_angle, thick, camber, cpos, &verts, )) } } }; // Характерный размер: калибр тела для канала и РАЗМЕР ДОМЕНА для эталонных течений — // именно так статьи определяют Re = u₀N/ν. let ref_len = if case.is_periodic() { cli.nx as R } else { scene.ref_size() }; let units = Units::new(cli.dx, cli.u_phys, cli.u_lat, cli.re, ref_len); let beta0 = math::beta_of_nu(units.nu_lat); let tau0 = 1.0 / (2.0 * beta0); if tau0 <= 0.5 { return Err(format!( "τ = {tau0:.4} ≤ ½: вязкость неположительна. Подними --u-lat или опусти --re" )); } if case.is_periodic() && cli.refine > 1 { return Err("эталонные течения периодичны и патча измельчения не имеют: --refine 1".into()); } // патч измельчения: по умолчанию охватывает тело и ближний след let patch = if cli.refine > 1 { let d = cli.size; let (ax, bx, ay, by) = match &cli.patch { Some(s) => { let p: Vec<&str> = s.split(',').collect(); if p.len() != 4 { return Err("--patch: ожидалось ax,bx,ay,by".into()); } let v: Result, _> = p.iter().map(|t| t.trim().parse::()).collect(); let v = v.map_err(|e| format!("--patch: {e}"))?; (v[0], v[1], v[2], v[3]) } None => ( (cx - 1.5 * d).round().max(2.0) as usize, (cx + 6.25 * d).round().min(cli.nx as R - 3.0) as usize, (cy - 2.0 * d).round().max(1.0) as usize, (cy + 2.0 * d).round().min(cli.ny as R - 2.0) as usize, ), }; // патч обязан лежать строго внутри домена и НЕ накрывать зеркальные ряды стенок: // на них ГУ работает по всему ряду, а рамка патча их бы перезаписала if !(1 <= ay && ay < by && by <= cli.ny - 2) { return Err(format!( "патч по y [{ay},{by}] обязан лежать внутри (0,{}) и не трогать ряды стенок; \ при теле {d} ячеек минимальное ny ≈ {}", cli.ny - 1, (4.0 * d) as usize + 6 )); } if !(2 <= ax && ax < bx && bx <= cli.nx - 2) { return Err(format!("патч по x [{ax},{bx}] обязан лежать внутри (1,{})", cli.nx - 1)); } Some((ax, bx, ay, by)) } else { None }; // Губка перед выходом. Задача канала «вход по скорости + выход по давлению» акустически // есть четвертьволновая труба: вход отражает продольные волны как жёсткий поршень, выход — // как открытый конец. Собственная мода затухает как ν(π/2Nx)², то есть практически не // затухает, и любой стартовый удар остаётся в домене навсегда. Губка — единственное, что // её реально ест, поэтому по умолчанию она включена и подобрана под длину домена. // губка входа не должна доставать до патча измельчения if cli.sponge_in > 0 { if let Some((ax, _, _, _)) = patch { if cli.sponge_in + 2 >= ax { return Err(format!( "губка входа ({} столбцов) достаёт до патча (ax={ax})", cli.sponge_in )); } } } let sponge_len = match cli.sponge_len { Some(v) => { if v > 0 { if let Some((_, bx, _, _)) = patch { if cli.nx - 1 - v <= bx { return Err(format!("губка (последние {v} столбцов) накрывает патч (bx={bx})")); } } } v } None => { let want = (cli.nx / 12).max(16); // губка не должна доставать до патча (и вообще до тела) let limit = match patch { Some((_, bx, _, _)) => cli.nx.saturating_sub(bx + 2), None => cli.nx / 4, }; want.min(limit) } }; let probe = ( ((cx + 3.0 * cli.size).round() as usize).min(cli.nx - 2), (cy.round() as usize).min(cli.ny - 2), ); // Длительность задаётся либо в шагах, либо в секундах. Шаг не есть фиксированная порция // времени: δt привязан к размеру клетки, поэтому одно и то же число шагов на разных // сетках покрывает разное физическое время. let steps = match (cli.steps, cli.time) { (Some(s), _) => s, (None, Some(t)) => { if !(t > 0.0) { return Err("--time обязан быть положительным".into()); } let n = (t / units.dt).round(); if n < 1.0 { return Err(format!( "--time {t} с при шаге δt = {:.4e} с даёт меньше одного шага", units.dt )); } if n > 2e9 { return Err(format!( "--time {t} с при шаге δt = {:.4e} с требует {n:.3e} шагов — это заведомо неподъёмно. Подними --u-lat или --dx, либо задай --steps явно", units.dt )); } n as u64 } (None, None) => 20_000, }; // В периодической постановке губок нет по определению: они поднимают вязкость у границ, // а границ здесь нет. Иначе эталон считался бы с завышенным затуханием. if case.is_periodic() { return Ok(Spec { nx: cli.nx, ny: cli.ny, scene, refine: 1, patch: None, units, beta0, steps, ramp: cli.ramp.max(1), flow_angle: cli.flow_angle * std::f64::consts::PI / 180.0, pert_amp: 0.0, pert_dur: 0, outlet_extrapolate: false, init_uniform: true, init_taper: 0.0, sponge_len: 0, sponge_in: 0, sponge_mult: 1.0, collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc }, kbc_model: math::KbcModel::from_str(&cli.kbc_model).ok_or("неизвестная модель KBC")?, wall: math::WallModel::from_str(&cli.wall).ok_or("неизвестная модель стенки")?, case, probe: (cli.nx / 2, cli.ny / 2), }); } Ok(Spec { nx: cli.nx, ny: cli.ny, scene, refine: cli.refine, patch, units, beta0, steps, ramp: cli.ramp.max(1), flow_angle: cli.flow_angle * std::f64::consts::PI / 180.0, pert_amp: cli.pert_amp, pert_dur: cli.pert_dur, outlet_extrapolate: cli.outlet == "extrapolate", init_uniform: cli.init == "uniform", init_taper: cli.init_taper.max(0.0), sponge_len, sponge_in: cli.sponge_in, sponge_mult: cli.sponge_mult, case, collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc }, kbc_model: math::KbcModel::from_str(&cli.kbc_model).ok_or("неизвестная модель KBC")?, wall: math::WallModel::from_str(&cli.wall).ok_or("неизвестная модель стенки")?, probe, }) } // ───────────────────────────────────────────────────────────────────────────── // Бэкенд // ───────────────────────────────────────────────────────────────────────────── enum Backend { Cpu(Box), #[cfg(feature = "gpu")] Gpu(Box), } impl Backend { /// Посчитать шаг. Готовые записи ДОПИСЫВАЮТСЯ в `out` — их может быть ноль (GPU копит /// итоги в истории и читает пачкой) или сразу много (когда история заполнилась). fn advance(&mut self, out: &mut Vec) { match self { Backend::Cpu(s) => out.push(s.step()), #[cfg(feature = "gpu")] Backend::Gpu(s) => s.advance(out), } } /// Дочитать всё посчитанное. Обязательно перед чтением полей: иначе записи отстанут /// от состояния, которое покажет кадр. fn flush(&mut self, out: &mut Vec) { let _ = &out; // процессорный бэкенд отдаёт записи сразу, копить нечего match self { Backend::Cpu(_) => {} #[cfg(feature = "gpu")] Backend::Gpu(s) => s.flush(out), } } fn sample_field(&mut self, k: FieldKind) -> (Vec, Vec) { match self { Backend::Cpu(s) => { let (v, m) = s.sample_field(k); (v, m.to_vec()) } #[cfg(feature = "gpu")] Backend::Gpu(s) => s.sample_field(k), } } fn sample_velocity(&mut self) -> (Vec, Vec) { match self { Backend::Cpu(s) => s.sample_velocity(), #[cfg(feature = "gpu")] Backend::Gpu(s) => s.sample_velocity(), } } fn sample_centerline(&mut self) -> (Vec, Vec) { match self { Backend::Cpu(s) => s.sample_centerline(), #[cfg(feature = "gpu")] Backend::Gpu(s) => s.sample_centerline(), } } fn is_finite(&mut self) -> bool { match self { Backend::Cpu(s) => s.is_finite(), #[cfg(feature = "gpu")] Backend::Gpu(s) => s.is_finite(), } } fn force_ref_size(&self) -> R { match self { Backend::Cpu(s) => s.force_ref_size(), #[cfg(feature = "gpu")] Backend::Gpu(s) => s.force_ref_size(), } } fn name(&self) -> &'static str { match self { Backend::Cpu(_) => "CPU (rayon, f64)", #[cfg(feature = "gpu")] Backend::Gpu(s) => s.name(), } } } // ───────────────────────────────────────────────────────────────────────────── // main // ───────────────────────────────────────────────────────────────────────────── fn main() { let cli = Cli::parse(); if let Err(e) = run(cli) { eprintln!("\nОШИБКА: {e}"); std::process::exit(1); } } fn run(cli: Cli) -> Result<(), String> { if cli.threads > 0 { rayon::ThreadPoolBuilder::new() .num_threads(cli.threads) .build_global() .map_err(|e| e.to_string())?; } let spec = build_spec(&cli)?; let quiet = cli.verbose == "quiet"; let full = cli.verbose == "full"; // ── анимация: план тайминга считается ДО прогона и печатается в шапке ── let stride = if cli.gif_every == "auto" { None } else { Some(cli.gif_every.parse::().map_err(|e| format!("--gif-every: {e}"))?) }; let plan = gif::GifPlan::new(&spec.units, stride, cli.gif_fps, cli.gif_speed); let field_kind = match cli.gif_field.as_str() { "speed" => FieldKind::Speed, "density" => FieldKind::Density, "gamma" => FieldKind::Gamma, _ => FieldKind::Vorticity, }; if !quiet { print_header(&spec, &cli, &plan, field_kind); } // ── бэкенд ── let mut back = match cli.backend.as_str() { "gpu" => { #[cfg(feature = "gpu")] { Backend::Gpu(Box::new(gpu::Sim::new(spec.clone())?)) } #[cfg(not(feature = "gpu"))] { return Err("сборка без поддержки GPU (--no-default-features)".into()); } } _ => Backend::Cpu(Box::new(cpu::Sim::new(spec.clone()))), }; if !quiet { println!("бэкенд: {}", back.name()); } // ── писатель гифки ── let mut writer = match &cli.gif { Some(path) => { let cmap = gif::ColorMap::from_str(&cli.gif_cmap).ok_or("неизвестная палитра")?; let range = match &cli.gif_range { Some(s) => { let (lo, hi) = parse_pair(s, "--gif-range")?; gif::Range { lo, hi } } None => default_range(field_kind, spec.units.u_lat, spec.units.d_lat), }; let hud = gif::Hud { u_phys: spec.units.u_phys, field: field_kind.label(), show: !cli.no_hud, }; Some( gif::GifWriter::create( path, spec.nx, spec.ny, cli.gif_scale, cli.gif_downsample, cmap, range, plan, hud, spec.patch, ) .map_err(|e| format!("не удалось создать {path}: {e}"))?, ) } None => None, }; // ── x–t диагностика ── let mut xt_file = match &cli.xt { Some(path) => { let mut f = std::io::BufWriter::new( std::fs::File::create(path).map_err(|e| format!("не создать {path}: {e}"))?, ); writeln!(f, "# срез вдоль осевой линии, столбцы: step,t_phys_s,field,v0..v{}", spec.nx - 1) .map_err(|e| e.to_string())?; Some(f) } None => None, }; // ── цикл ── let d_ref = back.force_ref_size(); let mut recs: Vec = Vec::with_capacity(spec.steps as usize + 1); let t_start = Instant::now(); let mut blew_up = false; let nodes_per_step = (spec.nx * spec.ny + spec .patch .map(|(ax, bx, ay, by)| { spec.refine * (spec.refine * (bx - ax) + 1) * (spec.refine * (by - ay) + 1) }) .unwrap_or(0)) as f64; let mut extras = Extras::default(); let series_every = cli.series_every.max(1); let mut scratch: Vec = Vec::with_capacity(256); let mut rec = StepRec::default(); for t in 0..=spec.steps { back.advance(&mut scratch); let need_frame = writer.is_some() && t % plan.stride == 0; let need_report = cli.report_every > 0 && !quiet && t % cli.report_every == 0; // и кадр, и живая строка требуют актуальных данных — синхронизируемся только здесь if need_frame || need_report { back.flush(&mut scratch); } for r in scratch.drain(..) { if r.step % series_every == 0 { recs.push(r); } rec = r; } if need_frame { if let Some(w) = writer.as_mut() { let (field, solid) = back.sample_field(field_kind); let (cd, _, _) = math::coefficients(rec.fx, rec.fy, rec.tz, spec.units.u_lat, d_ref); w.push(&field, &solid, spec.units.time_of_step(t), Some(cd)) .map_err(|e| format!("запись кадра: {e}"))?; } } if need_report { live_line(&spec, &rec, d_ref, t_start.elapsed().as_secs_f64(), nodes_per_step, full); } if spec.case.is_periodic() && t % cli.case_every.max(1) == 0 { back.flush(&mut scratch); let (vx, vy) = back.sample_velocity(); let (e, om, pa) = periodic_metrics(&vx, &vy, spec.nx, spec.ny); extras.case_series.push([t as R, spec.units.time_of_step(t), e, om, pa]); } if let Some(f) = xt_file.as_mut() { if t % cli.xt_every.max(1) == 0 { let (rho, ux) = back.sample_centerline(); let tp = spec.units.time_of_step(t); for (name, v) in [("rho", &rho), ("ux", &ux)] { write!(f, "{t},{tp:.9e},{name}").map_err(|e| e.to_string())?; for z in v { write!(f, ",{z:.7e}").map_err(|e| e.to_string())?; } writeln!(f).map_err(|e| e.to_string())?; } } } // Развал счёта ловится по УЖЕ ПОСЧИТАННЫМ величинам, а не полным проходом по полю: // на сетке 4096×2048 такой проход тянет с устройства сотни мегабайт, и делать его // регулярно нельзя. Полная проверка — один раз, для подтверждения. if !rec.max_u.is_finite() || !rec.rho_mean.is_finite() { back.flush(&mut scratch); eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ на шаге {t}: в поле появились NaN/inf."); blew_up = true; break; } } back.flush(&mut scratch); for r in scratch.drain(..) { if r.step % series_every == 0 { recs.push(r); } } if !blew_up && !back.is_finite() { eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ: итоговое поле содержит NaN/inf."); blew_up = true; } let wall = t_start.elapsed().as_secs_f64(); // Разовые итоги: радиус влияния (канал) и сверка с точным решением (Тейлор–Грин). if !blew_up { let (vx, vy) = back.sample_velocity(); if spec.case == Case::Channel { let (bcx, bcy) = spec.scene.center(); extras.influence = Some(influence_extent( &vx, &vy, spec.nx, spec.ny, spec.units.u_lat, bcx, bcy, spec.scene.ref_size(), 0.01, )); } if spec.case == Case::TaylorGreen { // точное решение: u_x = −u₀·cos(K₁x)·sin(K₂y)·exp(−ν(K₁²+K₂²)t) let (k1, k2) = (1.0, 4.0); let kk1 = 2.0 * std::f64::consts::PI * k1 / spec.nx as R; let kk2 = 2.0 * std::f64::consts::PI * k2 / spec.ny as R; let nu = math::nu_of_beta(spec.beta0); let e = (-nu * (kk1 * kk1 + kk2 * kk2) * spec.steps as R).exp(); let (mut num, mut den) = (0.0, 0.0); for y in 0..spec.ny { for x in 0..spec.nx { let want = -spec.units.u_lat * (kk1 * x as R).cos() * (kk2 * y as R).sin() * e; num += (vx[y * spec.nx + x] - want).abs(); den += want.abs(); } } extras.tg_error = Some(if den > 0.0 { num / den } else { R::NAN }); } } if let Some(path) = &cli.case_csv { let mut f = std::io::BufWriter::new( std::fs::File::create(path).map_err(|e| format!("не создать {path}: {e}"))?, ); writeln!(f, "step,t_phys_s,energy,enstrophy,palinstrophy").map_err(|e| e.to_string())?; for r in &extras.case_series { writeln!(f, "{:.0},{:.9e},{:.9e},{:.9e},{:.9e}", r[0], r[1], r[2], r[3], r[4]) .map_err(|e| e.to_string())?; } } if let Some(w) = writer { let frames = w.frames; w.finish().map_err(|e| format!("закрытие гифки: {e}"))?; if !quiet { println!("\nгифка: {} — кадров {frames}", cli.gif.as_ref().unwrap()); } } if let Some(path) = &cli.csv { write_csv(path, &recs, &spec, d_ref).map_err(|e| format!("CSV: {e}"))?; if !quiet { println!("ряды: {path}"); } } final_report( &spec, &cli, &recs, d_ref, wall, nodes_per_step, blew_up, &plan, full, &extras, ); if let Some(path) = &cli.summary { write_summary(path, &spec, &cli, &recs, d_ref, wall, nodes_per_step, blew_up, &extras) .map_err(|e| format!("сводка: {e}"))?; if !quiet { println!("сводка: {path}"); } } Ok(()) } /// Итоги, которые считаются не каждый шаг, а разово или редкой выборкой. #[derive(Default)] struct Extras { /// Ряд метрик эталонного течения: [шаг, t, энергия, энстрофия, палинстрофия]. case_series: Vec<[R; 5]>, /// Радиус влияния тела в калибрах: вверх по потоку, вбок, вниз по потоку. influence: Option<(R, R, R)>, /// Относительная ошибка поля u_x против точного решения Тейлора–Грина. tg_error: Option, } /// Интегральные метрики периодического течения: энергия, энстрофия, палинстрофия. /// /// E = ½⟨|u|²⟩, Ω = ½⟨ω²⟩, P = ½⟨|∇ω|²⟩; производные — центральные разности с заворотом, /// как и положено в периодическом ящике. По ним статьи и сравнивают модели: энстрофия ловит /// разрешение мелких вихрей, палинстрофия — их градиентов. fn periodic_metrics(ux: &[R], uy: &[R], nx: usize, ny: usize) -> (R, R, R) { let n = nx * ny; let idx = |x: usize, y: usize| y * nx + x; let mut w = vec![0.0; n]; let mut energy = 0.0; for y in 0..ny { for x in 0..nx { let (xp, xm) = ((x + 1) % nx, (x + nx - 1) % nx); let (yp, ym) = ((y + 1) % ny, (y + ny - 1) % ny); w[idx(x, y)] = 0.5 * (uy[idx(xp, y)] - uy[idx(xm, y)]) - 0.5 * (ux[idx(x, yp)] - ux[idx(x, ym)]); energy += ux[idx(x, y)] * ux[idx(x, y)] + uy[idx(x, y)] * uy[idx(x, y)]; } } let mut enst = 0.0; let mut pal = 0.0; for y in 0..ny { for x in 0..nx { let (xp, xm) = ((x + 1) % nx, (x + nx - 1) % nx); let (yp, ym) = ((y + 1) % ny, (y + ny - 1) % ny); enst += w[idx(x, y)] * w[idx(x, y)]; let gx = 0.5 * (w[idx(xp, y)] - w[idx(xm, y)]); let gy = 0.5 * (w[idx(x, yp)] - w[idx(x, ym)]); pal += gx * gx + gy * gy; } } let inv = 1.0 / n as R; (0.5 * energy * inv, 0.5 * enst * inv, 0.5 * pal * inv) } /// Насколько далеко тело возмущает поток: расстояния до самой дальней точки, где скорость /// отклоняется от набегающей больше чем на `thresh`, в калибрах тела. /// /// Отвечает на практический вопрос «какой домен достаточен»: если возмущение достаёт до /// границы, домен мал и результат зависит от его размера, а не от физики. fn influence_extent( ux: &[R], uy: &[R], nx: usize, ny: usize, u0: R, cx: R, cy: R, d: R, thresh: R, ) -> (R, R, R) { let lim = thresh * u0; let (mut up, mut lat, mut down): (R, R, R) = (0.0, 0.0, 0.0); for y in 0..ny { for x in 0..nx { let k = y * nx + x; let dev = ((ux[k] - u0).powi(2) + uy[k] * uy[k]).sqrt(); if dev <= lim { continue; } let (fx, fy) = (x as R, y as R); if fx < cx { up = up.max(cx - fx); } else { down = down.max(fx - cx); } lat = lat.max((fy - cy).abs()); } } let d = d.max(1e-9); (up / d, lat / d, down / d) } fn default_range(k: FieldKind, u_lat: R, d_lat: R) -> gif::Range { match k { FieldKind::Speed => gif::Range { lo: 0.0, hi: 1.7 * u_lat }, FieldKind::Vorticity => { let w = 5.0 * u_lat / d_lat; gif::Range { lo: -w, hi: w } } FieldKind::Density => gif::Range { lo: 0.985, hi: 1.015 }, FieldKind::Gamma => gif::Range { lo: 0.0, hi: 4.0 }, } } // ───────────────────────────────────────────────────────────────────────────── // Вывод // ───────────────────────────────────────────────────────────────────────────── fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) { let u = &spec.units; let tau = 1.0 / (2.0 * spec.beta0); let blockage = spec.scene.frontal_extent() / spec.ny as R; println!("╔══════════════════════════════════════════════════════════════════════════╗"); println!("║ KBC-2D — обтекание тела в канале, D2Q9 + энтропийное столкновение KBC-D ║"); println!("╚══════════════════════════════════════════════════════════════════════════╝"); println!("\n── Физика ──────────────────────────────────────────────────────────────────"); println!(" скорость потока {:>12.4} м/с направление {:>6.1}°", u.u_phys, cli.flow_angle); println!(" размер ячейки {:>12.5} м домен {:.2} × {:.2} м", u.dx, spec.nx as R * u.dx, spec.ny as R * u.dx); println!(" тело {:>12} {:.3} м ({:.0} ячеек), угол атаки {:.1}°", match spec.scene.bodies.len() { 0 => "нет (периодич.)".to_string(), 1 => format!("{:?}", spec.scene.bodies[0].kind).to_lowercase(), k => format!("{k} тел"), }, u.d_phys(), u.d_lat, cli.body_angle); println!(" Рейнольдс {:>12.1} вязкость {:.4e} м²/с", u.re, u.nu_phys()); println!(" Мах (реш.) {:>12.4} u_lat = {:.4}, τ = {:.4}", u.mach(), u.u_lat, tau); println!(" блокировка канала {:>12.3} (габарит тела / ширина канала)", blockage); println!("\n── Время ───────────────────────────────────────────────────────────────────"); println!(" шаг по времени {:>12.6e} с {:.1} шагов на секунду физ. времени", u.dt, u.steps_per_second()); println!(" всего шагов {:>12} = {:.4} с физ. времени = {:.1} конв. времён D/U", spec.steps, u.time_of_step(spec.steps), u.convective_times(spec.steps)); println!(" разгон / возмущение {:>12} {} шагов импульса амплитудой {:.2}·U", spec.ramp, spec.pert_dur, spec.pert_amp); println!("\n── Сетка ───────────────────────────────────────────────────────────────────"); println!(" L0 {:>12} {} узлов", format!("{}×{}", spec.nx, spec.ny), spec.nx * spec.ny); match spec.patch { Some((ax, bx, ay, by)) => { let (nfx, nfy) = (spec.refine * (bx - ax) + 1, spec.refine * (by - ay) + 1); let tau1 = spec.refine as R * (tau - 0.5) + 0.5; println!(" L1 (×{}) {:>12} область x∈[{ax},{bx}] y∈[{ay},{by}], τ₁={:.4}", spec.refine, format!("{nfx}×{nfy}"), tau1); println!(" {:>12} {} подшагов на шаг L0", "", spec.refine); } None => println!(" L1 {:>12} (измельчение выключено)", "нет"), } println!(" оператор {:>12}", if spec.collision == Collision::Kbc { match spec.kbc_model { math::KbcModel::N1 => "KBC-N1 (D)", math::KbcModel::N2 => "KBC-N2 (C)", } } else { "LBGK" }); // Действительно ли субсеточная граница работает: если большинство линков сваливается в // простой отскок, стенка де-факто ступенчатая, как её ни называй. { let g = cpu::Geom::build(spec.nx, spec.ny, &spec.scene); let w = g.wall_stats(); let pc = |k: usize| 100.0 * k as R / w.links.max(1) as R; println!(" ── Граница тела ────────────────────────────────────────────────────────────"); println!(" граничных узлов {:>12} линков {}", w.nodes, w.links); println!(" интерполяция {:>11.1}% (ближняя {:.1}% + дальняя {:.1}%)", pc(w.near) + pc(w.far), pc(w.near), pc(w.far)); println!(" простой отскок {:>11.1}% {}", pc(w.simple), if w.simple == 0 { "— ступенчатых линков нет" } else { "← ЭТИ ЛИНКИ СТУПЕНЧАТЫЕ" }); println!(" доля пересечения q {:>12} мин {:.3}, средн {:.3}, макс {:.3}; зажато {:.1}%", "", w.q_min, w.q_mean, w.q_max, pc(w.q_clipped)); } println!(" стенка тела {:>12}", match spec.wall { math::WallModel::Bouzidi => "Bouzidi", math::WallModel::Grad => "Град (моменты)", math::WallModel::Hrr => "HRR (моменты, 3-й пор.)", math::WallModel::Staircase => "простой отскок", }); println!(" выход по u_y {:>12} губка {} столбцов ×{:.0}", if spec.outlet_extrapolate { "extrapolate" } else { "zero" }, spec.sponge_len, spec.sponge_mult); // Продольная акустика — не мелочь оформления, а то, что видно глазом: при заметной // амплитуде поток «дышит» рывками. Печатаем период моды и то, чем она гасится. { let cs = math::CS2.sqrt(); let period = 4.0 * spec.nx as R / cs; // четвертьволновая труба: lambda = 4L let nu = math::nu_of_beta(spec.beta0); let k = std::f64::consts::PI / (2.0 * spec.nx as R); let efold = 1.0 / (nu * k * k); println!(" ── Продольная акустика канала ──────────────────────────────────────────────"); println!(" вход по скорости отражает звук как жёсткий поршень, выход по давлению — как"); println!(" открытый конец: канал работает четвертьволновым резонатором."); println!(" период основной моды {:>10.0} шагов = {:.4} с", period, period * spec.units.dt); println!(" затухание вязкостью {:>10.0} шагов на е-кратное — то есть сама не гаснет", efold); print!(" чем гасится: "); if spec.init_uniform { print!("старт из однородного потока (моду не возбуждаем)"); } else { print!("РАЗГОН ИЗ ПОКОЯ за {} шагов", spec.ramp); if (spec.ramp as R) < period { print!(" — короче периода моды, будет звон!"); } } if spec.sponge_len > 0 || spec.sponge_in > 0 { println!( " + губка: вход {} / выход {} столбцов ×{:.0}", spec.sponge_in, spec.sponge_len, spec.sponge_mult ); } else { println!("; губка ВЫКЛЮЧЕНА"); } } // Пара «вход по скорости / выход по давлению» — недодемпфированный акустический резонатор: // затухание продольной моды ~ ν(π/Nx)², то есть на длинном домене её почти ничто не гасит, // и стартовый транзиент способен раскачать её до смещённого режима или до развала счёта. if cli.gif.is_some() { println!("\n── Анимация (синхронизация с физическим временем) ───────────────────────────"); println!(" поле {:>12} палитра {}", field.label().to_lowercase(), cli.gif_cmap); println!(" шаг кадра {:>12} {}", plan.stride, if plan.auto { "подобран автоматически" } else { "задан явно" }); println!(" между кадрами {:>12.6e} с физ. времени", plan.frame_dt); println!(" задержка кадра {:>12.4} сотых {:.2} кадр/с{}", plan.exact_delay_cs, plan.fps, if plan.dithered { " (задержки чередуются — точное время без ухода)" } else { "" }); println!(" воспроизведение {:>12.4}× от реального времени (запрошено {:.4}×)", plan.actual_playback, plan.requested_playback); match plan.advice(&spec.units) { Some(a) => println!(" ⚠ {a}"), None if plan.is_synced() => { println!(" ✓ синхронно с физическим временем (расхождение {:.4}%)", (plan.sync_error() - 1.0).abs() * 100.0) } None => println!(" ⚠ расхождение с реальным временем {:.3}×", plan.sync_error()), } } println!(); } fn live_line(spec: &Spec, r: &StepRec, d_ref: R, wall: f64, nodes: f64, full: bool) { let (cd, cl, _) = math::coefficients(r.fx, r.fy, r.tz, spec.units.u_lat, d_ref); let done = (r.step + 1) as f64; let sps = done / wall.max(1e-9); let eta = (spec.steps + 1 - r.step.min(spec.steps)) as f64 / sps.max(1e-9); let mlups = sps * nodes / 1e6; // во сколько раз счёт быстрее реального времени let rt = sps / spec.units.steps_per_second(); let head = format!( "шаг {:>8} t={:>9.4}с ⟨ρ⟩={:.5} max|u|={:.4} Cd={:>7.3} Cl={:>7.3}", r.step, spec.units.time_of_step(r.step), r.rho_mean, r.max_u, cd, cl ); if full { println!( "{head}\n ⟨γ⟩={:.3} [{:.2}…{:.2}] вырожд={:.2}% ξ<0={:.2}% \ {:.0} шаг/с ({:.1} MLUPS, {:.2}× реального) ETA {:.0}с", r.gamma_mean, r.gamma_min, r.gamma_max, r.degenerate_frac * 100.0, r.xi_negative_frac * 100.0, sps, mlups, rt, eta ); } else { println!("{head} {:.0} шаг/с ETA {:.0}с", sps, eta); } std::io::stdout().flush().ok(); } #[allow(clippy::too_many_arguments)] fn final_report( spec: &Spec, cli: &Cli, recs: &[StepRec], d_ref: R, wall: f64, nodes: f64, blew_up: bool, plan: &gif::GifPlan, full: bool, extras: &Extras, ) { let u = &spec.units; let n = recs.len(); if n < 8 { println!("слишком мало шагов для отчёта"); return; } let h = n / 2; // установившийся режим — вторая половина ряда let cd: Vec = recs.iter().map(|r| 2.0 * r.fx / (u.u_lat * u.u_lat * d_ref)).collect(); let cl: Vec = recs.iter().map(|r| 2.0 * r.fy / (u.u_lat * u.u_lat * d_ref)).collect(); let cm: Vec = recs.iter().map(|r| 2.0 * r.tz / (u.u_lat * u.u_lat * d_ref * d_ref)).collect(); let uy: Vec = recs.iter().map(|r| r.uy_probe).collect(); let rho: Vec = recs.iter().map(|r| r.rho_mean).collect(); let (cd_m, _) = math::mean_std(&cd[h..]); let (_, cl_rms) = math::mean_std(&cl[h..]); let (cm_m, _) = math::mean_std(&cm[h..]); let (_, uy_rms) = math::mean_std(&uy[h..]); let (st, _, _) = math::strouhal(&uy, u.d_lat, u.u_lat, 8192); println!("\n╔══════════════════════════════════════════════════════════════════════════╗"); println!("║ ИТОГОВЫЙ ОТЧЁТ ║"); println!("╚══════════════════════════════════════════════════════════════════════════╝"); if blew_up { println!("⚠ ПРОГОН ОБОРВАН: счёт развалился. Числа ниже относятся к тому, что успело сойтись."); } // ── эталонные течения: свои метрики вместо Cd/St ── if spec.case.is_periodic() { println!("\n── Эталонное течение: {} ─────────────────────────────────────────", cli.case); if let (Some(a), Some(b)) = (extras.case_series.first(), extras.case_series.last()) { println!(" {:>10}{:>14}{:>14}{:>16}", "момент", "энергия", "энстрофия", "палинстрофия"); println!(" {:>10}{:>14.6e}{:>14.6e}{:>16.6e}", "старт", a[2], a[3], a[4]); println!(" {:>10}{:>14.6e}{:>14.6e}{:>16.6e}", "конец", b[2], b[3], b[4]); if a[2] > 0.0 && a[3] > 0.0 { println!( " отношение к начальному: энергия {:.4}, энстрофия {:.4}", b[2] / a[2], b[3] / a[3] ); } } if let Some(e) = extras.tg_error { println!(" относительная ошибка u_x против ТОЧНОГО решения: {e:.4e}"); println!(" (метрика рис. 1 статьи: Σ|u_x − u_x^точн| / Σ|u_x^точн|)"); } println!("\n── Производительность ──────────────────────────────────────────────────────"); let sps = n as f64 / wall.max(1e-9); println!(" время счёта {wall:.1} с {sps:.0} шаг/с {:.1} MLUPS", sps * nodes / 1e6); println!(); return; } // ── установившийся режим ── let beta = spec.scene.frontal_extent() / spec.ny as R; println!("\n── Установившийся режим (вторая половина ряда, шаги {}–{}) ──────────────────", h, n - 1); println!(" {:<22}{:>10}{:>10}{:>10}{:>10}", "величина", "St", "⟨Cd⟩", "rms Cl", "⟨Cm⟩"); println!(" {:<22}{:>10.4}{:>10.4}{:>10.4}{:>10.5}", "как посчитано", st, cd_m, cl_rms, cm_m); // поправка на стеснение канала: St кинематическое ⇒ ×(1−β); Cd и Cl ~ скорость² ⇒ ×(1−β)² println!( " {:<22}{:>10.4}{:>10.4}{:>10.4}{:>10}", "с поправ. на блокир.", st * (1.0 - beta), cd_m * (1.0 - beta).powi(2), cl_rms * (1.0 - beta).powi(2), "—" ); if spec.scene.bodies.first().map(|b| b.kind) == Some(ShapeKind::Cylinder) && spec.scene.bodies.len() == 1 && u.re > 100.0 && u.re < 200.0 { println!(" {:<22}{:>10.4}{:>10.4}{:>10}{:>10.1}", "литература Re≈150", 0.183, 1.330, "~0.30", 0.0); } println!(" rms u_y в зонде = {uy_rms:.5} (узел {:?}); блокировка β = {beta:.3}", spec.probe); println!(" поправка: St×(1−β), Cd и rms Cl ×(1−β)². Литература — безграничный цилиндр;"); println!(" ⟨Cm⟩≈0 для симметричного тела под нулевым углом — это контроль симметрии считывания."); // ── силы по телам ── // Для тандема и многоэлементных конфигураций литература даёт коэффициенты КАЖДОГО тела // по отдельности, поэтому общей суммы недостаточно. let nb = spec.scene.bodies.len().min(MAX_BODY_BUCKETS); if nb > 1 { println!("\n── Силы по телам (вторая половина ряда) ────────────────────────────────────"); println!(" {:>5}{:>12}{:>10}{:>10}{:>11}", "тело", "положение", "⟨Cd⟩", "rms Cl", "⟨Cm⟩"); for b in 0..nb { let cdb: Vec = recs[h..].iter().map(|r| 2.0 * r.body[b][0] / (u.u_lat * u.u_lat * d_ref)).collect(); let clb: Vec = recs[h..].iter().map(|r| 2.0 * r.body[b][1] / (u.u_lat * u.u_lat * d_ref)).collect(); let cmb: Vec = recs[h..] .iter() .map(|r| 2.0 * r.body[b][2] / (u.u_lat * u.u_lat * d_ref * d_ref)) .collect(); let (m, _) = math::mean_std(&cdb); let (_, sl) = math::mean_std(&clb); let (mm, _) = math::mean_std(&cmb); let pos = spec .scene .bodies .get(b) .map(|bd| format!("{:.0},{:.0}", bd.cx, bd.cy)) .unwrap_or_else(|| "—".into()); println!(" {:>5}{:>12}{:>10.4}{:>10.4}{:>11.5}", b + 1, pos, m, sl, mm); } if spec.scene.bodies.len() > MAX_BODY_BUCKETS { println!( " (тел {}, вёдер {} — тела с {}-го и дальше сведены в последнюю строку;", spec.scene.bodies.len(), MAX_BODY_BUCKETS, MAX_BODY_BUCKETS ); println!(" сумма по строкам при этом точна и равна общему ⟨Cd⟩ выше)"); } } // ── сходимость ── let nw = cli.windows.max(2).min(n / 2); let w = n / nw; println!("\n── Сходимость по времени (окна по {w} шагов) ────────────────────────────────"); println!(" {:>4}{:>16}{:>10}{:>10}{:>10}{:>10}{:>10}", "окно", "шаги", "⟨ρ⟩", "⟨Cd⟩", "⟨Cd⟩/ρ", "rms Cl", "rms u_y"); let mut rows = Vec::new(); for k in 0..nw { let s = k * w; let e = if k + 1 == nw { n } else { (k + 1) * w }; let (rm, _) = math::mean_std(&rho[s..e]); let (cdm, _) = math::mean_std(&cd[s..e]); let (_, clr) = math::mean_std(&cl[s..e]); let (_, uyr) = math::mean_std(&uy[s..e]); let cdn = if rm != 0.0 { cdm / rm } else { R::NAN }; rows.push((rm, cdm, cdn, clr, uyr)); println!(" {:>4}{:>16}{:>10.5}{:>10.4}{:>10.4}{:>10.4}{:>10.5}", k + 1, format!("{s}–{e}"), rm, cdm, cdn, clr, uyr); } let (rm1, _, cdn1, clr1, _) = rows[nw - 1]; let (rm0, _, _, clr0, _) = rows[nw - 2]; print!(" вердикт: "); if (rm1 - 1.0).abs() > 0.02 { // «стабильно, но не там»: система села на смещённую ветвь (акустическая накачка массы), // на ней Cd/Cl несопоставимы с литературой print!("⟨ρ⟩={rm1:.4} — РЕЖИМ СМЕЩЁН, прогон невалиден для сравнения; "); } else if (rm1 - rm0).abs() > 1e-3 { print!("⟨ρ⟩={rm1:.4}, Δ за окно {:+.5} — ДРЕЙФ МАССЫ (правь ГУ выхода); ", rm1 - rm0); } else { print!("масса стабильна (⟨ρ⟩={rm1:.4}, Δ {:+.5}); ", rm1 - rm0); } if clr1 - clr0 > 1e-3 { println!("rms Cl ещё растёт ({:+.4}) — не насыщено; дрейф-устойчивый Cd = {cdn1:.3}", clr1 - clr0); } else { println!("rms Cl насыщен ({:+.4}); дрейф-устойчивый Cd = {cdn1:.3}", clr1 - clr0); } if let Some((up, lat, down)) = extras.influence { println!("\n── Радиус влияния тела (отклонение скорости больше 1% от U) ────────────────"); println!(" вверх по потоку {up:.1} калибра, вбок {lat:.1}, вниз по потоку {down:.1}"); let (bcx, _) = spec.scene.center(); let d = spec.scene.ref_size().max(1e-9); let (room_up, room_lat, room_down) = ( bcx / d, (spec.ny as R / 2.0) / d, (spec.nx as R - bcx) / d, ); let tight = up > 0.9 * room_up || lat > 0.9 * room_lat || down > 0.9 * room_down; println!( " до границ домена: {room_up:.1} / {room_lat:.1} / {room_down:.1} калибра{}", if tight { " ⚠ ВОЗМУЩЕНИЕ ДОСТАЁТ ДО ГРАНИЦЫ — домен мал" } else { "" } ); } // ── измеренная продольная пульсация ── // Колебание средней плотности — это и есть та самая «поршневая» мода. Переводим его в // амплитуду скорости: для бегущей акустической волны u' = (δρ/ρ)·c_s. { let cs = math::CS2.sqrt(); let last = &rho[n - n / 10..]; let lo = last.iter().cloned().fold(R::INFINITY, R::min); let hi = last.iter().cloned().fold(R::NEG_INFINITY, R::max); let mean = last.iter().sum::() / last.len() as R; let u_ac = 0.5 * (hi - lo) / mean * cs; let ratio = u_ac / u.u_lat; println!(" ── Продольная пульсация (последняя десятая часть прогона) ───────────────────"); println!(" размах ⟨ρ⟩ = {:.5} вокруг {mean:.5} ⇒ амплитуда скорости u' = {:.2}% от U", hi - lo, ratio * 100.0); if !ratio.is_finite() { println!(" ⚠ ПУЛЬСАЦИЯ НЕ ОПРЕДЕЛЕНА: в рядах NaN/inf, счёт развалился."); } else if ratio.abs() > 0.05 { println!(" ⚠ ПОТОК ЗАМЕТНО ПУЛЬСИРУЕТ: это продольная мода канала, а не физика следа."); println!(" Лечится стартом из однородного потока (--init uniform) и губкой"); println!(" (--sponge-len {}). Числа Cd/St на таком прогоне недостоверны.", (spec.nx / 12).max(16)); } else { println!(" пульсация в пределах нормы — продольная мода не раскачана"); } } // ── диагностика KBC ── if spec.collision == Collision::Kbc { let (g_m, g_s) = math::mean_std(&recs[h..].iter().map(|r| r.gamma_mean).collect::>()); let gmin = recs[h..].iter().map(|r| r.gamma_min).fold(R::INFINITY, R::min); let gmax = recs[h..].iter().map(|r| r.gamma_max).fold(R::NEG_INFINITY, R::max); let (deg, _) = math::mean_std(&recs[h..].iter().map(|r| r.degenerate_frac).collect::>()); let (xin, _) = math::mean_std(&recs[h..].iter().map(|r| r.xi_negative_frac).collect::>()); println!("\n── Энтропийный стабилизатор γ (формула (17) статьи) ─────────────────────────"); println!(" ⟨γ⟩ = {g_m:.4} ± {g_s:.4} размах по узлам [{gmin:.3} … {gmax:.3}]"); println!(" вырожденных узлов (⟨Δh|Δh⟩→0, γ подменён на 2 = локальный LBGK): {:.3}%", deg * 100.0); println!(" узлов с отрицательной объёмной вязкостью ξ = c_s²(1/(γβ) − ½): {:.2}%", xin * 100.0); if full { println!(" Пояснение. γ ≡ 2 — это в точности LBGK, поэтому доля вырожденных узлов есть"); println!(" доля решётки, где KBC ничего не даёт; при относительном пороге она обязана быть"); println!(" ~0 (абсолютный порог загонял её в 77–99% и молча превращал схему в LBGK)."); println!(" ξ по формуле (57) у моделей B и D зависит от γ, то есть от точки и времени;"); println!(" ξ < 0 означает локальное антизатухание акустики. У моделей A и C ξ = ν всегда."); } } // ── производительность и синхронность ── let sps = n as f64 / wall.max(1e-9); println!("\n── Производительность ──────────────────────────────────────────────────────"); println!(" время счёта {wall:.1} с {sps:.0} шаг/с {:.1} MLUPS", sps * nodes / 1e6); println!(" физического времени просчитано {:.4} с ⇒ счёт идёт в {:.3}× от реального времени", u.time_of_step(n as u64 - 1), sps / u.steps_per_second()); if cli.gif.is_some() { println!(" гифка воспроизводится в {:.4}× от реального времени — это НЕ зависит от скорости", plan.actual_playback); println!(" счёта: задержка кадра берётся из δt, а не из wall-clock."); } println!(); } fn write_csv(path: &str, recs: &[StepRec], spec: &Spec, d_ref: R) -> std::io::Result<()> { let mut f = std::io::BufWriter::new(std::fs::File::create(path)?); writeln!(f, "step,t_phys_s,cd,cl,cm,uy_probe,rho_mean,max_u,gamma_mean,gamma_min,gamma_max,degenerate_frac,xi_negative_frac")?; let s = spec.units.u_lat * spec.units.u_lat * d_ref; for r in recs { writeln!( f, "{},{:.9e},{:.6e},{:.6e},{:.6e},{:.6e},{:.9},{:.6e},{:.6},{:.6},{:.6},{:.6},{:.6}", r.step, spec.units.time_of_step(r.step), 2.0 * r.fx / s, 2.0 * r.fy / s, 2.0 * r.tz / (s * d_ref), r.uy_probe, r.rho_mean, r.max_u, r.gamma_mean, r.gamma_min, r.gamma_max, r.degenerate_frac, r.xi_negative_frac )?; } Ok(()) } /// Машиночитаемая сводка прогона. Пишется вручную, без serde: полей немного, а лишняя /// зависимость в решателе не нужна. #[allow(clippy::too_many_arguments)] fn write_summary( path: &str, spec: &Spec, cli: &Cli, recs: &[StepRec], d_ref: R, wall: f64, nodes: f64, blew_up: bool, extras: &Extras, ) -> std::io::Result<()> { let u = &spec.units; let n = recs.len(); let h = n / 2; let q = |v: R| if v.is_finite() { format!("{v:.9e}") } else { "null".into() }; let mut f = std::io::BufWriter::new(std::fs::File::create(path)?); writeln!(f, "{{")?; writeln!(f, " \"case\": \"{}\",", cli.case)?; writeln!(f, " \"backend\": \"{}\",", cli.backend)?; writeln!(f, " \"wall\": \"{}\", \"kbc_model\": \"{}\", \"collision\": \"{}\",", cli.wall, cli.kbc_model, cli.collision)?; writeln!(f, " \"nx\": {}, \"ny\": {}, \"refine\": {},", spec.nx, spec.ny, spec.refine)?; writeln!(f, " \"bodies\": {}, \"body_size_cells\": {},", spec.scene.bodies.len(), q(spec.scene.ref_size()))?; writeln!(f, " \"re\": {}, \"u_lat\": {}, \"mach\": {}, \"tau\": {},", q(u.re), q(u.u_lat), q(u.mach()), q(1.0 / (2.0 * spec.beta0)))?; writeln!(f, " \"steps\": {}, \"time_phys_s\": {}, \"convective_times\": {},", spec.steps, q(u.time_of_step(spec.steps)), q(u.convective_times(spec.steps)))?; writeln!(f, " \"blew_up\": {},", blew_up)?; let sps = n as f64 / wall.max(1e-9); writeln!(f, " \"wall_time_s\": {:.3}, \"steps_per_s\": {:.1}, \"mlups\": {:.2},", wall, sps, sps * nodes / 1e6)?; if n >= 8 { if spec.case.is_periodic() { if let (Some(a), Some(b)) = (extras.case_series.first(), extras.case_series.last()) { writeln!(f, " \"energy_start\": {}, \"energy_end\": {},", q(a[2]), q(b[2]))?; writeln!(f, " \"enstrophy_start\": {}, \"enstrophy_end\": {},", q(a[3]), q(b[3]))?; writeln!(f, " \"palinstrophy_end\": {},", q(b[4]))?; writeln!(f, " \"enstrophy_ratio\": {},", q(if a[3] > 0.0 { b[3] / a[3] } else { R::NAN }))?; } writeln!(f, " \"taylor_green_error\": {},", extras.tg_error.map(q).unwrap_or_else(|| "null".into()))?; } else { let sc = u.u_lat * u.u_lat * d_ref; let cd: Vec = recs[h..].iter().map(|r| 2.0 * r.fx / sc).collect(); let cl: Vec = recs[h..].iter().map(|r| 2.0 * r.fy / sc).collect(); let cm: Vec = recs[h..].iter().map(|r| 2.0 * r.tz / (sc * d_ref)).collect(); let uy: Vec = recs.iter().map(|r| r.uy_probe).collect(); let rho: Vec = recs[h..].iter().map(|r| r.rho_mean).collect(); let (st, _, _) = math::strouhal(&uy, u.d_lat, u.u_lat, 8192); let (cdm, _) = math::mean_std(&cd); let (_, clr) = math::mean_std(&cl); let (cmm, _) = math::mean_std(&cm); let (rm, _) = math::mean_std(&rho); let bl = spec.scene.frontal_extent() / spec.ny as R; writeln!(f, " \"strouhal\": {}, \"cd\": {}, \"cl_rms\": {}, \"cm\": {},", q(st), q(cdm), q(clr), q(cmm))?; writeln!(f, " \"blockage\": {}, \"strouhal_corr\": {}, \"cd_corr\": {}, \"cl_rms_corr\": {},", q(bl), q(st * (1.0 - bl)), q(cdm * (1.0 - bl).powi(2)), q(clr * (1.0 - bl).powi(2)))?; writeln!(f, " \"rho_mean\": {},", q(rm))?; // по телам write!(f, " \"bodies_cd\": [")?; let nb = spec.scene.bodies.len().min(MAX_BODY_BUCKETS); for b in 0..nb { let v: Vec = recs[h..].iter().map(|r| 2.0 * r.body[b][0] / sc).collect(); let (m, _) = math::mean_std(&v); write!(f, "{}{}", if b > 0 { ", " } else { "" }, q(m))?; } writeln!(f, "],")?; let (up, lat, down) = extras.influence.unwrap_or((R::NAN, R::NAN, R::NAN)); writeln!(f, " \"influence_up\": {}, \"influence_lat\": {}, \"influence_down\": {},", q(up), q(lat), q(down))?; } let cs = math::CS2.sqrt(); let last = &recs[n - n / 10..]; let lo = last.iter().map(|r| r.rho_mean).fold(R::INFINITY, R::min); let hi = last.iter().map(|r| r.rho_mean).fold(R::NEG_INFINITY, R::max); let mean = last.iter().map(|r| r.rho_mean).sum::() / last.len() as R; writeln!(f, " \"pulsation_u_over_U\": {},", q(0.5 * (hi - lo) / mean * cs / u.u_lat))?; let (gm, _) = math::mean_std(&recs[h..].iter().map(|r| r.gamma_mean).collect::>()); let (dg, _) = math::mean_std(&recs[h..].iter().map(|r| r.degenerate_frac).collect::>()); let (xn, _) = math::mean_std(&recs[h..].iter().map(|r| r.xi_negative_frac).collect::>()); writeln!(f, " \"gamma_mean\": {}, \"degenerate_frac\": {}, \"xi_negative_frac\": {}", q(gm), q(dg), q(xn))?; } else { writeln!(f, " \"note\": \"слишком мало шагов для статистики\"")?; } writeln!(f, "}}")?; Ok(()) }