diff --git a/docs/theory/2d_solver/src/cpu.rs b/docs/theory/2d_solver/src/cpu.rs index 9017e3a..b2e2fa0 100644 --- a/docs/theory/2d_solver/src/cpu.rs +++ b/docs/theory/2d_solver/src/cpu.rs @@ -9,7 +9,169 @@ use rayon::prelude::*; use crate::math::{self, Kbc, Link, LinkKind, Scene, WallModel, MAX_BODY_BUCKETS, R, CX, CY, OPP, Q}; -use crate::{Collision, FieldKind, Spec, StepRec}; +use crate::{Case, Collision, FieldKind, Spec, StepRec}; + +// ───────────────────────────────────────────────────────────────────────────── +// Стартовые поля +// ───────────────────────────────────────────────────────────────────────────── + +/// Что нужно, чтобы разложить стартовое поле уровня. +pub struct Init<'a> { + pub case: Case, + pub nx: usize, + pub ny: usize, + /// Скорость набегающего потока (канал) либо амплитуда эталонного течения. + pub u0: (R, R), + pub beta: R, + pub scene: Option<&'a Scene>, + /// Ширина сведения скорости к нулю у тела; 0 — не сводить. + pub taper: R, +} + +/// Стартовое поле в виде популяций. Строится на хосте для ОБОИХ бэкендов: GPU просто +/// перекладывает результат в свою раскладку. Так эталонные течения из статей задаются +/// ровно одним куском кода. +pub fn initial_field(p: &Init) -> Vec<[R; Q]> { + let (nx, ny) = (p.nx, p.ny); + let n = nx * ny; + let tau = 1.0 / (2.0 * p.beta); + let nu = math::nu_of_beta(p.beta); + let two_pi = 2.0 * std::f64::consts::PI; + let mut f = vec![[0.0; Q]; n]; + + match p.case { + Case::Channel => { + let plain = math::feq(1.0, p.u0.0, p.u0.1); + for (k, cell) in f.iter_mut().enumerate() { + *cell = match (p.scene, p.taper > 0.0) { + (Some(sc), true) => { + let (x, y) = ((k % nx) as R, (k / nx) as R); + let c = math::wall_taper(sc.sdf(x, y), p.taper); + if c < 1.0 { + math::feq(1.0, p.u0.0 * c, p.u0.1 * c) + } else { + plain + } + } + _ => plain, + }; + } + } + + // Разд. VI 2D-статьи: u = ∇×[(u₀/k₂)cos(k₁x)cos(k₂y)·exp(−ν(k₁²+k₂²)t)], k₁=1, k₂=4. + // Стартуем приближением Града и с СОГЛАСОВАННЫМ давлением: течение несёт собственное + // поле p ~ ρu₀², и старт с ρ ≡ 1 сбрасывает эту разницу в акустику, которая в + // периодическом ящике почти не затухает. + Case::TaylorGreen => { + let (k1, k2) = (1.0, 4.0); + let kk1 = two_pi * k1 / nx as R; + let kk2 = two_pi * k2 / ny as R; + let a = p.u0.0; + for y in 0..ny { + for x in 0..nx { + let (px, py) = (kk1 * x as R, kk2 * y as R); + let ux = -a * px.cos() * py.sin(); + let uy = (k1 / k2) * a * px.sin() * py.cos(); + let pr = -(a * a / 4.0) + * ((2.0 * px).cos() + (k1 * k1) / (k2 * k2) * (2.0 * py).cos()); + let rho = 1.0 + pr / math::CS2; + let dxux = a * kk1 * px.sin() * py.sin(); + let dyux = -a * kk2 * px.cos() * py.cos(); + let dxuy = (k1 / k2) * a * kk1 * px.cos() * py.cos(); + let dyuy = -(k1 / k2) * a * kk2 * px.sin() * py.sin(); + let pref = -tau * rho * math::CS2; + f[y * nx + x] = math::grad_init( + rho, + ux, + uy, + rho * math::CS2 + rho * ux * ux + pref * 2.0 * dxux, + rho * ux * uy + pref * (dyux + dxuy), + rho * math::CS2 + rho * uy * uy + pref * 2.0 * dyuy, + ); + } + } + } + + // Разд. VII: дважды периодический сдвиговый слой, κ = 80, δ = 0.05. + Case::ShearLayer => { + let (kappa, delta) = (80.0, 0.05); + let a = p.u0.0; + for y in 0..ny { + let yy = y as R / ny as R; + let ux = if y as R <= ny as R / 2.0 { + a * (kappa * (yy - 0.25)).tanh() + } else { + a * (kappa * (0.75 - yy)).tanh() + }; + for x in 0..nx { + let uy = delta * a * (two_pi * (x as R / nx as R + 0.25)).sin(); + f[y * nx + x] = math::feq(1.0, ux, uy); + } + } + } + + // Затухающая турбулентность: функция тока собирается из мод в узкой полосе вокруг + // k₀, фазы — детерминированный ГПСЧ (прогон обязан быть воспроизводимым). + Case::DecayingTurbulence => { + let k0 = 8.0; + let a = p.u0.0; + let mut seed = 0x2545_F491_4F6C_DD1Du64; + let mut rnd = || { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + (seed >> 11) as R / (1u64 << 53) as R + }; + let mut modes = Vec::new(); + for ky in -12i32..=12 { + for kx in -12i32..=12 { + if kx == 0 && ky == 0 { + continue; + } + let km = ((kx * kx + ky * ky) as R).sqrt(); + if !(k0 - 4.0..=k0 + 4.0).contains(&km) || kx < 0 { + continue; + } + let amp = (-(km - k0) * (km - k0) / 8.0).exp() / km; + modes.push((kx as R, ky as R, amp, two_pi * rnd())); + } + } + let norm: R = modes.iter().map(|m| m.2 * m.2).sum::().sqrt().max(1e-30); + for y in 0..ny { + for x in 0..nx { + let (fx, fy) = (x as R / nx as R, y as R / ny as R); + let (mut ux, mut uy) = (0.0, 0.0); + for &(kx, ky, amp, ph) in &modes { + let arg = two_pi * (kx * fx + ky * fy) + ph; + let c = arg.cos() * amp / norm; + // u = ∇×(ψ ẑ) для ψ = Σ amp·sin(2π(k·r)+φ) + ux += c * two_pi * ky; + uy -= c * two_pi * kx; + } + f[y * nx + x] = math::feq(1.0, a * ux, a * uy); + } + } + // приводим к заданной среднеквадратичной скорости + let rms = (f + .iter() + .map(|c| { + let (_, ux, uy) = math::macros(c); + ux * ux + uy * uy + }) + .sum::() + / n as R) + .sqrt() + .max(1e-30); + let scale = a / rms; + for cell in f.iter_mut() { + let (_, ux, uy) = math::macros(cell); + *cell = math::feq(1.0, ux * scale, uy * scale); + } + } + } + let _ = nu; + f +} // ───────────────────────────────────────────────────────────────────────────── // Геометрия уровня @@ -192,14 +354,14 @@ impl Level { /// `u0` — скорость, которой заполняется поле на старте. Заполнять сразу набегающим /// потоком принципиально: старт из покоя разгоняет весь столб жидкости и закачивает в /// канал продольную акустическую моду, которую вязкость потом почти не гасит. - fn new(nx: usize, ny: usize, beta: Vec, geom: Geom, u0: (R, R)) -> Self { - let f0 = math::feq(1.0, u0.0, u0.1); + fn new(nx: usize, ny: usize, beta: Vec, geom: Geom, f0: Vec<[R; Q]>) -> Self { let n = nx * ny; + debug_assert_eq!(f0.len(), n); Level { nx, ny, - f: vec![f0; n], - post: vec![f0; n], + post: f0.clone(), + f: f0, gamma: vec![2.0; n], beta, geom, @@ -472,20 +634,6 @@ impl Level { } } -/// Свести стартовую скорость к нулю на подходе к телу, чтобы на первом шаге не возникло -/// разрыва. Вдали от тела поле не трогается вовсе. -fn taper_start(lvl: &mut Level, scene: &Scene, u0: (R, R), width: R) { - let nx = lvl.nx; - for (n, cell) in lvl.f.iter_mut().enumerate() { - let (x, y) = ((n % nx) as R, (n / nx) as R); - let k = math::wall_taper(scene.sdf(x, y), width); - if k < 1.0 { - *cell = math::feq(1.0, u0.0 * k, u0.1 * k); - } - } - lvl.post.copy_from_slice(&lvl.f); -} - // ───────────────────────────────────────────────────────────────────────────── // Статистика KBC за шаг // ───────────────────────────────────────────────────────────────────────────── @@ -770,10 +918,16 @@ impl Sim { } else { (0.0, 0.0) }; - let mut l0 = Level::new(nx, ny, beta0, geom0, u0); - if spec.init_uniform { - taper_start(&mut l0, &spec.scene, u0, spec.init_taper); - } + let init0 = initial_field(&Init { + case: spec.case, + nx, + ny, + u0, + beta: spec.beta0, + scene: Some(&spec.scene), + taper: if spec.init_uniform { spec.init_taper } else { 0.0 }, + }); + let l0 = Level::new(nx, ny, beta0, geom0, init0); let fluid_count = l0.geom.solid.iter().filter(|s| !**s).count() as R; let (l1, patch) = if spec.refine > 1 { @@ -789,11 +943,16 @@ impl Sim { let beta1 = 1.0 / (2.0 * tau1); let r01 = tau1 / (r as R * tau0); let patch = Patch::new(&spec, &l0.geom, &geom1.solid, r01); - let mut lvl1 = Level::new(nfx, nfy, vec![beta1; nfx], geom1, u0); - if spec.init_uniform { - taper_start(&mut lvl1, &scene1, u0, spec.init_taper * r as R); - } - (Some(lvl1), Some(patch)) + let init1 = initial_field(&Init { + case: spec.case, + nx: nfx, + ny: nfy, + u0, + beta: 1.0 / (2.0 * tau1), + scene: Some(&scene1), + taper: if spec.init_uniform { spec.init_taper * r as R } else { 0.0 }, + }); + (Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1, init1)), Some(patch)) } else { (None, None) }; @@ -860,9 +1019,13 @@ impl Sim { let wall = self.spec.wall; let stats0 = self.l0.collide(sp_collision, model); self.l0.stream(); - self.l0.apply_wall(wall); - self.l0.free_slip_walls(); - self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap); + // Эталонные течения статей периодичны по обеим осям и ГУ не имеют вовсе: перенос + // уже периодичен, поэтому достаточно ничего не накладывать. + if self.spec.case == Case::Channel { + self.l0.apply_wall(wall); + self.l0.free_slip_walls(); + self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap); + } // ── уровень 1: r подшагов с временной интерполяцией рамки ── let mut fb = [[0.0 as R; 3]; MAX_BODY_BUCKETS]; @@ -872,7 +1035,9 @@ impl Sim { for s in 0..p.r { l1.collide(sp_collision, model); l1.stream(); - l1.apply_wall(wall); + if self.spec.case == Case::Channel { + l1.apply_wall(wall); + } // силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на // последнем подшаге даёт лишний шум в рядах при том же среднем let g = l1.force(); @@ -989,6 +1154,19 @@ impl Sim { (out, &self.l0.geom.solid) } + /// Полное поле скорости уровня L0 — для метрик эталонных течений и радиуса влияния. + pub fn sample_velocity(&self) -> (Vec, Vec) { + let n = self.l0.nx * self.l0.ny; + let mut ux = Vec::with_capacity(n); + let mut uy = Vec::with_capacity(n); + for k in 0..n { + let (_, a, b) = self.l0.macros_at(k); + ux.push(a); + uy.push(b); + } + (ux, uy) + } + /// Срез вдоль осевой линии: (ρ, u_x) по каждому столбцу. Из последовательности таких /// срезов складывается x–t диаграмма, по которой видно, бежит возмущение со скоростью /// звука или конвекции и есть ли стоячие узлы. @@ -1033,7 +1211,7 @@ mod tests { body_cx: 0.0, body_cy: 0.0, }, - (0.0, 0.0), + vec![math::feq(1.0, 0.0, 0.0); nx * ny], ) } diff --git a/docs/theory/2d_solver/src/gif.rs b/docs/theory/2d_solver/src/gif.rs index 72a373b..b880e29 100644 --- a/docs/theory/2d_solver/src/gif.rs +++ b/docs/theory/2d_solver/src/gif.rs @@ -369,6 +369,12 @@ pub struct Hud { pub struct GifWriter { enc: gif::Encoder>, pub plan: GifPlan, + /// Размеры ИСХОДНОЙ сетки. + src_nx: usize, + src_ny: usize, + /// Во сколько раз клетки усредняются в пиксель перед отрисовкой. + down: usize, + /// Размеры картинки после прореживания. nx: usize, ny: usize, scale: usize, @@ -388,6 +394,7 @@ impl GifWriter { nx: usize, ny: usize, scale: usize, + down: usize, cmap: ColorMap, range: Range, plan: GifPlan, @@ -395,6 +402,11 @@ impl GifWriter { patch: Option<(usize, usize, usize, usize)>, ) -> std::io::Result { let scale = scale.max(1); + let down = down.max(1); + let (src_nx, src_ny) = (nx, ny); + // прореживание с округлением вверх: последний блок может быть неполным + let nx = nx.div_ceil(down); + let ny = ny.div_ceil(down); let w = nx * scale; let h = ny * scale; let file = BufWriter::new(File::create(path)?); @@ -404,6 +416,9 @@ impl GifWriter { Ok(GifWriter { enc, plan, + src_nx, + src_ny, + down, nx, ny, scale, @@ -426,6 +441,34 @@ impl GifWriter { t_phys: R, cd: Option, ) -> std::io::Result<()> { + // Прореживание: блок down×down усредняется в один пиксель, а телом пиксель считается, + // если тело занимает хотя бы половину блока. Без этого кадр с сетки 4096×2048 весит + // столько, что гифка становится непригодной. + let (field, solid) = if self.down == 1 { + (field.to_vec(), solid.to_vec()) + } else { + let d = self.down; + let mut fv = vec![0.0; self.nx * self.ny]; + let mut sv = vec![false; self.nx * self.ny]; + for gy in 0..self.ny { + for gx in 0..self.nx { + let (mut acc, mut cnt, mut sol) = (0.0, 0usize, 0usize); + for yy in gy * d..((gy + 1) * d).min(self.src_ny) { + for xx in gx * d..((gx + 1) * d).min(self.src_nx) { + let k = yy * self.src_nx + xx; + acc += field[k]; + sol += solid[k] as usize; + cnt += 1; + } + } + let g = gy * self.nx + gx; + fv[g] = if cnt > 0 { acc / cnt as R } else { 0.0 }; + sv[g] = cnt > 0 && 2 * sol >= cnt; + } + } + (fv, sv) + }; + let (field, solid) = (&field[..], &solid[..]); let mut buf = vec![0u8; self.w * self.h]; // строка 0 изображения — это ВЕРХ, а y = 0 решётки — низ канала: переворачиваем for y in 0..self.ny { @@ -441,7 +484,8 @@ impl GifWriter { } } if let Some((ax, bx, ay, by)) = self.patch { - self.draw_patch_outline(&mut buf, ax, bx, ay, by); + let d = self.down; + self.draw_patch_outline(&mut buf, ax / d, bx / d, ay / d, by / d); } if self.hud.show { let line = match cd { diff --git a/docs/theory/2d_solver/src/gpu.rs b/docs/theory/2d_solver/src/gpu.rs index 50b96ea..fa54b11 100644 --- a/docs/theory/2d_solver/src/gpu.rs +++ b/docs/theory/2d_solver/src/gpu.rs @@ -21,7 +21,7 @@ use wgpu::util::DeviceExt; use crate::cpu; use crate::math::{self, R}; -use crate::{Collision, FieldKind, Spec, StepRec}; +use crate::{Case, Collision, FieldKind, Spec, StepRec}; const WG: u32 = 64; @@ -804,31 +804,14 @@ impl GpuLevel { } } -/// Стартовое поле в раскладке SoA. Если задана сцена, скорость сводится к нулю на подходе -/// к телу (см. `math::wall_taper`) — иначе на первом шаге возникает разрыв и импульс сжатия. -fn soa_equilibrium( - n: usize, - nx: usize, - u0: (R, R), - taper: Option<(&math::Scene, R)>, -) -> Vec { +/// Переложить готовое стартовое поле (построенное общим кодом в `cpu::initial_field`) +/// из AoS в раскладку SoA, которой пользуется GPU. +fn to_soa(f: &[[R; math::Q]]) -> Vec { + let n = f.len(); let mut v = vec![0.0f32; 9 * n]; - let plain = math::feq(1.0, u0.0, u0.1); - for k in 0..n { - let fe = match taper { - Some((scene, w)) => { - let (x, y) = ((k % nx) as R, (k / nx) as R); - let c = math::wall_taper(scene.sdf(x, y), w); - if c < 1.0 { - math::feq(1.0, u0.0 * c, u0.1 * c) - } else { - plain - } - } - None => plain, - }; + for (k, cell) in f.iter().enumerate() { for i in 0..9 { - v[i * n + k] = fe[i] as f32; + v[i * n + k] = cell[i] as f32; } } v @@ -844,12 +827,11 @@ fn make_level( beta: &[f32], probe_node: u32, flags: u32, - u0: (R, R), + init_field: &[[R; math::Q]], wall_layout: &wgpu::BindGroupLayout, - taper: Option<(&math::Scene, R)>, ) -> GpuLevel { let n = nx * ny; - let init = soa_equilibrium(n, nx, u0, taper); + let init = to_soa(init_field); let mkf = |label: &str| { device.create_buffer_init(&wgpu::util::BufferInitDescriptor { label: Some(label), @@ -885,7 +867,11 @@ fn make_level( _p: 0, }) .collect(); - let links_buf = storage_init(device, "links", bytemuck::cast_slice(&links)); + // Пустая сцена (эталонные течения) даёт нулевой список линков, а шейдер всё равно + // объявляет массив структур: буфер обязан вмещать хотя бы один элемент, иначе валидация + // ругается на несоответствие размера. Читать его при этом некому — nlinks = 0. + let links_pad = if links.is_empty() { vec![GLink::zeroed()] } else { links.clone() }; + let links_buf = storage_init(device, "links", bytemuck::cast_slice(&links_pad)); let beta_buf = storage_init(device, "beta", bytemuck::cast_slice(beta)); // индекс граничных узлов: (узел, смещение первого линка, число линков, —) @@ -895,7 +881,8 @@ fn make_level( .map(|w| [w.node, w.first, w.count as u32, 0]) .collect(); let nwall = wnodes.len() as u32; - let wall_buf = storage_init(device, "wall nodes", bytemuck::cast_slice(&wnodes)); + let wnodes_pad = if wnodes.is_empty() { vec![[0u32; 4]] } else { wnodes.clone() }; + let wall_buf = storage_init(device, "wall nodes", bytemuck::cast_slice(&wnodes_pad)); let wall_bind = device.create_bind_group(&wgpu::BindGroupDescriptor { label: Some("wall"), layout: wall_layout, @@ -1155,9 +1142,16 @@ impl Sim { &beta0, probe0 as u32, flags0, - u0, + &cpu::initial_field(&cpu::Init { + case: spec.case, + nx: spec.nx, + ny: spec.ny, + u0, + beta: spec.beta0, + scene: Some(&spec.scene), + taper: if spec.init_uniform { spec.init_taper } else { 0.0 }, + }), &bgl_wall, - if spec.init_uniform { Some((&spec.scene, spec.init_taper)) } else { None }, ); let pre = device.create_buffer(&wgpu::BufferDescriptor { @@ -1186,13 +1180,19 @@ impl Sim { let probe1 = if probe_on_fine { ((py - ay) * r) * nfx + (px - ax) * r } else { 0 }; let flags1 = if probe_on_fine { 2 } else { 0 }; let lvl1 = - make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1, - u0, &bgl_wall, - if spec.init_uniform { - Some((&scene1, spec.init_taper * r as R)) - } else { - None - }); + make_level( + &device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1, + &cpu::initial_field(&cpu::Init { + case: spec.case, + nx: nfx, + ny: nfy, + u0, + beta: 1.0 / (2.0 * tau1), + scene: Some(&scene1), + taper: if spec.init_uniform { spec.init_taper * r as R } else { 0.0 }, + }), + &bgl_wall, + ); let ghosts: Vec = patch .ghosts() @@ -1363,6 +1363,7 @@ impl Sim { let (ux_in, uy_in) = self.inlet(t); let refine = self.spec.refine.max(1); let moment_wall = self.spec.wall.is_moment_based(); + let channel = self.spec.case == Case::Channel; self.queue.write_buffer( &self.dyn_buf, 0, @@ -1402,20 +1403,23 @@ impl Sim { p.dispatch_workgroups(ncell, 1, 1); p.set_pipeline(&self.pipes.stream); p.dispatch_workgroups(ncell, 1, 1); - if moment_wall { - p.set_bind_group(2, &self.pipes.empty_bg, &[]); - p.set_bind_group(3, &self.l0.wall_bind, &[]); - p.set_pipeline(&self.pipes.moment_wall); - p.dispatch_workgroups(ceil_div(self.l0.nwall.max(1), WG), 1, 1); - } else { - p.set_pipeline(&self.pipes.bouzidi); - p.dispatch_workgroups(self.l0.link_groups(), 1, 1); + // Эталонные течения статей периодичны по обеим осям: ГУ не накладываются вовсе. + if channel { + if moment_wall { + p.set_bind_group(2, &self.pipes.empty_bg, &[]); + p.set_bind_group(3, &self.l0.wall_bind, &[]); + p.set_pipeline(&self.pipes.moment_wall); + p.dispatch_workgroups(ceil_div(self.l0.nwall.max(1), WG), 1, 1); + } else { + p.set_pipeline(&self.pipes.bouzidi); + p.dispatch_workgroups(self.l0.link_groups(), 1, 1); + } + // порядок обязателен: стенки снимают заворот по y, затем Zou-He — по x + p.set_pipeline(&self.pipes.walls); + p.dispatch_workgroups(ceil_div(self.l0.nx as u32, WG), 1, 1); + p.set_pipeline(&self.pipes.channel); + p.dispatch_workgroups(ceil_div(self.l0.ny as u32, WG), 1, 1); } - // порядок обязателен: стенки снимают заворот по y, затем Zou-He — по x, на весь столбец - p.set_pipeline(&self.pipes.walls); - p.dispatch_workgroups(ceil_div(self.l0.nx as u32, WG), 1, 1); - p.set_pipeline(&self.pipes.channel); - p.dispatch_workgroups(ceil_div(self.l0.ny as u32, WG), 1, 1); } if let (Some(l1), Some(a)) = (self.l1.as_ref(), self.amr.as_ref()) { @@ -1446,14 +1450,16 @@ impl Sim { p.dispatch_workgroups(ncell1, 1, 1); p.set_pipeline(&self.pipes.stream); p.dispatch_workgroups(ncell1, 1, 1); - if moment_wall { - p.set_bind_group(2, &self.pipes.empty_bg, &[]); - p.set_bind_group(3, &l1.wall_bind, &[]); - p.set_pipeline(&self.pipes.moment_wall); - p.dispatch_workgroups(ceil_div(l1.nwall.max(1), WG), 1, 1); - } else { - p.set_pipeline(&self.pipes.bouzidi); - p.dispatch_workgroups(nlink1, 1, 1); + if channel { + if moment_wall { + p.set_bind_group(2, &self.pipes.empty_bg, &[]); + p.set_bind_group(3, &l1.wall_bind, &[]); + p.set_pipeline(&self.pipes.moment_wall); + p.dispatch_workgroups(ceil_div(l1.nwall.max(1), WG), 1, 1); + } else { + p.set_pipeline(&self.pipes.bouzidi); + p.dispatch_workgroups(nlink1, 1, 1); + } } // силу снимаем на каждом подшаге, усредняется она в k_stats2 p.set_pipeline(&self.pipes.force); @@ -1612,6 +1618,21 @@ impl Sim { out } + /// Полное поле скорости уровня L0 — для метрик эталонных течений и радиуса влияния. + pub fn sample_velocity(&self) -> (Vec, Vec) { + let n = self.l0.n; + let raw = self.download_l0(); + let mut ux = Vec::with_capacity(n); + let mut uy = Vec::with_capacity(n); + for k in 0..n { + let c: [R; 9] = std::array::from_fn(|i| raw[i * n + k] as R); + let (_, a, b) = math::macros(&c); + ux.push(a); + uy.push(b); + } + (ux, uy) + } + /// Срез вдоль осевой линии: (ρ, u_x) по столбцам. На GPU это полное скачивание поля, /// поэтому x–t диагностика включается редким шагом и только там, где нужна. pub fn sample_centerline(&self) -> (Vec, Vec) { diff --git a/docs/theory/2d_solver/src/main.rs b/docs/theory/2d_solver/src/main.rs index 770ef23..62e1c1c 100644 --- a/docs/theory/2d_solver/src/main.rs +++ b/docs/theory/2d_solver/src/main.rs @@ -33,6 +33,38 @@ pub enum Collision { 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 { @@ -105,6 +137,8 @@ pub struct Spec { pub kbc_model: math::KbcModel, /// Чем замыкаются недостающие популяции на теле (см. `math::WallModel`). pub wall: math::WallModel, + /// Постановка задачи. + pub case: Case, /// Узел зонда следа в координатах L0. pub probe: (usize, usize), } @@ -140,6 +174,15 @@ struct Cli { /// Направление потока, градусы (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, ячеек @@ -292,6 +335,10 @@ struct Cli { /// Целочисленное увеличение картинки #[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, @@ -326,6 +373,13 @@ struct Cli { /// Период записи срезов 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, } // ───────────────────────────────────────────────────────────────────────────── @@ -402,6 +456,7 @@ fn parse_poly(v: &str) -> Result, String> { } 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()); } @@ -417,7 +472,12 @@ fn build_spec(cli: &Cli) -> Result { 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 = match &cli.bodies { + // В эталонных течениях тела нет вовсе: сцена пуста, маска твёрдого пуста, ГУ не + // накладываются. Иначе цилиндр по умолчанию стоял бы прямо посреди эталона. + 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()) { @@ -457,9 +517,13 @@ fn build_spec(cli: &Cli) -> Result { &verts, )) } + } }; - let units = Units::new(cli.dx, cli.u_phys, cli.u_lat, cli.re, scene.ref_size()); + // Характерный размер: калибр тела для канала и РАЗМЕР ДОМЕНА для эталонных течений — + // именно так статьи определяют 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 { @@ -468,6 +532,9 @@ fn build_spec(cli: &Cli) -> Result { )); } + if case.is_periodic() && cli.refine > 1 { + return Err("эталонные течения периодичны и патча измельчения не имеют: --refine 1".into()); + } // патч измельчения: по умолчанию охватывает тело и ближний след let patch = if cli.refine > 1 { let d = cli.size; @@ -577,6 +644,36 @@ fn build_spec(cli: &Cli) -> Result { (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, @@ -596,6 +693,7 @@ fn build_spec(cli: &Cli) -> Result { 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("неизвестная модель стенки")?, @@ -643,6 +741,13 @@ impl Backend { 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(), @@ -754,6 +859,7 @@ fn run(cli: Cli) -> Result<(), String> { spec.nx, spec.ny, cli.gif_scale, + cli.gif_downsample, cmap, range, plan, @@ -792,6 +898,7 @@ fn run(cli: Cli) -> Result<(), String> { }) .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(); @@ -823,6 +930,12 @@ fn run(cli: Cli) -> Result<(), String> { 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(); @@ -859,6 +972,54 @@ fn run(cli: Cli) -> Result<(), String> { } 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}"))?; @@ -874,10 +1035,103 @@ fn run(cli: Cli) -> Result<(), String> { } } - final_report(&spec, &cli, &recs, d_ref, wall, nodes_per_step, blew_up, &plan, full); + 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 }, @@ -908,10 +1162,10 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) { println!(" размер ячейки {:>12.5} м домен {:.2} × {:.2} м", u.dx, spec.nx as R * u.dx, spec.ny as R * u.dx); println!(" тело {:>12} {:.3} м ({:.0} ячеек), угол атаки {:.1}°", - if spec.scene.bodies.len() > 1 { - format!("{} тел", spec.scene.bodies.len()) - } else { - format!("{:?}", spec.scene.bodies[0].kind).to_lowercase() + 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()); @@ -1084,6 +1338,7 @@ fn final_report( blew_up: bool, plan: &gif::GifPlan, full: bool, + extras: &Extras, ) { let u = &spec.units; let n = recs.len(); @@ -1113,6 +1368,33 @@ fn final_report( 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); @@ -1127,8 +1409,8 @@ fn final_report( cl_rms * (1.0 - beta).powi(2), "—" ); - if spec.scene.bodies.len() == 1 - && spec.scene.bodies[0].kind == ShapeKind::Cylinder + 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 { @@ -1212,6 +1494,23 @@ fn final_report( 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. @@ -1298,3 +1597,104 @@ fn write_csv(path: &str, recs: &[StepRec], spec: &Spec, d_ref: R) -> std::io::Re } 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(()) +}