diff --git a/docs/origins/Grad's_aproximation.pdf b/docs/origins/Grad's_aproximation.pdf new file mode 100644 index 0000000..844b7d3 Binary files /dev/null and b/docs/origins/Grad's_aproximation.pdf differ diff --git a/docs/theory/2d_solver/README.md b/docs/theory/2d_solver/README.md index 948530e..210a19b 100644 --- a/docs/theory/2d_solver/README.md +++ b/docs/theory/2d_solver/README.md @@ -83,6 +83,19 @@ cargo test --release -- --include-ignored # плюс эталоны стат `--gif-speed 0.1` даёт замедление в 10 раз (тоже точно, через тот же накопитель). +## Шаг — не фиксированная порция времени + +Физическая длительность шага привязана к размеру клетки: δt = u_lat·δx/u_phys. Уменьшили клетку +вдвое — вдвое уменьшился и δt, и то же число шагов покроет вдвое меньше физического времени. +Отсюда `--time`: задаёте длительность в секундах, число шагов считается само. + +Измельчение стоит дважды: клеток становится (1/δx)², а шагов на ту же секунду — 1/δx, итого +работы ~(1/δx)³ в двумерии. И помните, что `--size` задаётся В КЛЕТКАХ: уменьшив клетку и не +тронув `--size`, вы уменьшите тело физически. + +Судить о длительности удобнее всего по **конвективным временам D/U** — их печатает шапка. Это +единственная мера, не зависящая ни от сетки, ни от выбора u_lat. + ## Что параметризовано - **поток**: скорость (м/с), направление, число Рейнольдса, решёточная скорость (число Маха); @@ -90,10 +103,12 @@ cargo test --release -- --include-ignored # плюс эталоны стат вложенного патча и его границы; - **тело**: семь форм на выбор — `cylinder`, `square`, `diamond`, `ellipse`, `naca`, `triangle`, `plate` — плюс характерный размер, относительная толщина, угол атаки и положение; -- **время**: число шагов, начальное поле (однородный поток либо покой с разгоном), длина - разгона, амплитуда и длительность стартового возмущения; -- **схема**: оператор столкновения (`kbc`/`bgk`), состав сдвиговой части (`n1`/`n2`), режим - выхода, поглощающая губка, бэкенд, число потоков; +- **время**: длительность прогона — либо числом шагов (`--steps`), либо прямо в СЕКУНДАХ + физического времени (`--time`, число шагов считается как time/δt); начальное поле (однородный + поток либо покой с разгоном), длина разгона, амплитуда и длительность стартового возмущения; +- **схема**: оператор столкновения (`kbc`/`bgk`), состав сдвиговой части (`n1`/`n2`), модель + стенки на теле (`bouzidi`/`grad`/`staircase`), режим выхода, поглощающая губка, бэкенд, + число потоков; - **анимация**: файл, поле (`speed`/`vorticity`/`density`/`gamma`), палитра, масштаб, шаг кадра, частота, скорость воспроизведения, диапазон нормировки; - **вывод**: период живых строк, три уровня подробности, число окон в отчёте о сходимости, CSV. @@ -210,6 +225,62 @@ c_s²(1/(γβ) − ½) и, поскольку измеренная ⟨γ⟩ ≈ постановке (старт из покоя, губка выключена) `n1` доживает до конца с пульсацией 75% от U, а `n2` разваливается. Поэтому умолчание — `n1`. +### Модель стенки на теле: насколько она субсеточная + +Ключ `--wall`: + +- `bouzidi` (умолчание) — интерполированный отскок: доля пересечения q входит в КАЖДУЮ + восстанавливаемую популяцию, полинково; +- `grad` — условие Града (Dorschner, Bösch, Chikatamarla, Boulouchos, Karlin, JFM 801 (2016), + разд. 2.1 и прил. B): задаются не популяции, а целевые моменты — ρ, u и тензор давлений, — + после чего недостающие популяции собираются приближением Града (2.13); +- `staircase` — простой отскок, q игнорируется. Не для счёта: это база сравнения, показывающая, + сколько именно даёт субсеточность. + +**Субсеточность — измеренная.** Прямой тест: сдвигаем тело внутри клетки и смотрим, насколько +поедет Cd. У по-настоящему субсеточной границы ответ не должен зависеть от того, где тело стоит +относительно узлов (Re = 20, D = 16, стационар, сдвиги 0…½ клетки): + +| модель | разброс Cd | Cd | +|---|---|---| +| `staircase` | 1.11% | 2.68–2.71 | +| `grad` | 0.64% | 2.64–2.66 | +| `bouzidi` | **0.19%** | 2.653–2.658 | + +Град оказывается ровно между ступенькой и Bouzidi, и это следует из его устройства: положение +стенки входит туда ТОЛЬКО через целевую скорость (B 1) — одну усреднённую по узлу величину. +Целевая плотность (B 3) — обычная сумма отскочивших и известных популяций, без q вовсе; тензор +давлений — конечные разности по решётке, тоже без q. Плюс все недостающие популяции узла +собираются из ОДНОГО набора моментов, так что полинковая направленность теряется. Bouzidi же +подставляет свою q в каждую популяцию отдельно. Ступенчатой поверхность у Града не становится, +но геометрия у него разрешена заметно грубее. + +**Сходимость по разрешению тела.** Физическая постановка фиксирована (домен 15D × 10D, +блокировка 0.1, Re = 20), меняется только число клеток на диаметр: + +| D | `bouzidi` | `grad` | +|---|---|---| +| 8 | 2.581 | 2.618 | +| 16 | 2.529 | 2.534 | +| 32 | **2.521** | **2.521** | + +Обе модели состоятельны и сходятся к одному пределу с наблюдаемым порядком ≈2.7; к D = 32 они +неразличимы. Но на грубой сетке Град заметно хуже: ошибка при D = 8 равна 0.097 против 0.060. +Для сравнения, `staircase` при D = 16 даёт 2.69 — то есть +6.7% к пределу, тогда как обе +субсеточные модели держатся в пределах +0.4%. + +**Зачем тогда Град.** Его преимущество в статье — не геометрическая точность, а устойчивость +на турбулентных режимах (авторы пишут, что интерполяционные схемы «ограничены низкими числами +Рейнольдса, поскольку на границе возникают паразитные скачки») и естественная форма для +подвижных стенок: скорость стенки входит в целевые значения, а не отдельной поправкой. +В здешней канальной постановке преимущества по устойчивости воспроизвести не удалось: при росте +Re обе модели теряют счёт на одном и том же значении (Re ≈ 5·10⁴ при теле в 16 клеток), то есть +ограничивает не стенка, а что-то другое — вероятнее всего Zou–He при τ → ½. Поэтому умолчание — +`bouzidi`, а `grad` стоит держать в виду для будущих подвижных тел. + +Ограничение реализации: `--wall grad` пока только на процессорном бэкенде; GPU при таком выборе +отказывается запускаться явно, а не считает молча по Bouzidi. + ### Паритет бэкендов и согласованность уровней CPU (f64) и GPU (f32) на одной постановке совпадают до 4–5 значащих цифр шаг в шаг: ⟨ρ⟩ diff --git a/docs/theory/2d_solver/src/cpu.rs b/docs/theory/2d_solver/src/cpu.rs index fc399e5..c10b45b 100644 --- a/docs/theory/2d_solver/src/cpu.rs +++ b/docs/theory/2d_solver/src/cpu.rs @@ -8,7 +8,7 @@ use rayon::prelude::*; -use crate::math::{self, Body, Kbc, Link, LinkKind, R, CX, CY, OPP, Q}; +use crate::math::{self, Body, Kbc, Link, LinkKind, WallModel, R, CX, CY, OPP, Q}; use crate::{Collision, FieldKind, Spec, StepRec}; // ───────────────────────────────────────────────────────────────────────────── @@ -22,6 +22,8 @@ pub struct Geom { pub solid: Vec, /// Линки «жидкий → твёрдый» с готовой геометрией пересечения. pub links: Vec, + /// Индекс по граничным узлам: условие Града ставится сразу на весь узел, а не полинково. + pub wall_nodes: Vec, /// Плечи для момента: координаты узлов относительно центра тела. pub body_cx: R, pub body_cy: R, @@ -41,7 +43,8 @@ impl Geom { } } let links = build_links(nx, ny, &solid, &phi); - Geom { nx, ny, solid, links, body_cx: body.cx, body_cy: body.cy } + let wall_nodes = build_wall_nodes(&links); + Geom { nx, ny, solid, links, wall_nodes, body_cx: body.cx, body_cy: body.cy } } } @@ -93,6 +96,70 @@ fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R]) -> Vec { links } +/// Сводка по границе тела: действительно ли работает субсеточная интерполяция. +#[derive(Clone, Copy, Debug, Default)] +pub struct WallStats { + pub nodes: usize, + pub links: usize, + /// q < ½ и есть дальний жидкий сосед — полная интерполяционная формула. + pub near: usize, + /// q ≥ ½ — интерполяция по двум популяциям того же узла. + pub far: usize, + /// Дальнего жидкого соседа нет: откат на ПРОСТОЙ отскок, то есть ступенчатую стенку. + pub simple: usize, + pub q_min: R, + pub q_max: R, + pub q_mean: R, + /// Линки, где q упёрлась в клип [0.02, 0.98] — узел практически лежит на поверхности. + pub q_clipped: usize, +} + +impl Geom { + pub fn wall_stats(&self) -> WallStats { + let mut w = WallStats { q_min: R::INFINITY, q_max: R::NEG_INFINITY, ..Default::default() }; + w.nodes = self.wall_nodes.len(); + w.links = self.links.len(); + for l in &self.links { + match l.kind { + LinkKind::Near => w.near += 1, + LinkKind::Far => w.far += 1, + LinkKind::Simple => w.simple += 1, + } + w.q_min = w.q_min.min(l.q); + w.q_max = w.q_max.max(l.q); + w.q_mean += l.q; + if l.q <= 0.0201 || l.q >= 0.9799 { + w.q_clipped += 1; + } + } + if w.links > 0 { + w.q_mean /= w.links as R; + } + w + } +} + +/// Жидкий узел, в который после переноса пришла бы хоть одна популяция из тела. +/// Линки одного узла в `Geom::links` лежат подряд (сборка идёт по y, потом x, потом +/// направлениям), поэтому достаточно смещения и количества. +#[derive(Clone, Copy)] +pub struct WallNode { + pub node: u32, + pub first: u32, + pub count: u8, +} + +fn build_wall_nodes(links: &[Link]) -> Vec { + let mut out: Vec = Vec::new(); + for (k, l) in links.iter().enumerate() { + match out.last_mut() { + Some(w) if w.node == l.node => w.count += 1, + _ => out.push(WallNode { node: l.node, first: k as u32, count: 1 }), + } + } + out +} + // ───────────────────────────────────────────────────────────────────────────── // Уровень сетки // ───────────────────────────────────────────────────────────────────────────── @@ -119,7 +186,15 @@ impl Level { fn new(nx: usize, ny: usize, beta: Vec, geom: Geom, u0: (R, R)) -> Self { let f0 = math::feq(1.0, u0.0, u0.1); let n = nx * ny; - Level { nx, ny, f: vec![f0; n], post: vec![f0; n], gamma: vec![2.0; n], beta, geom } + Level { + nx, + ny, + f: vec![f0; n], + post: vec![f0; n], + gamma: vec![2.0; n], + beta, + geom, + } } /// Столкновение по всем жидким узлам. Внутри тела не считаем: эти популяции фиктивны @@ -172,6 +247,14 @@ impl Level { let _ = ny; } + /// Простой отскок по всем линкам: q игнорируется. База сравнения для субсеточных моделей. + fn apply_staircase(&mut self) { + for l in &self.geom.links { + let node = l.node as usize; + self.f[node][l.ib as usize] = self.post[node][l.i as usize]; + } + } + /// Интерполированный отскок Bouzidi по всем линкам тела. /// `f` — поле ПОСЛЕ переноса (его правим), `post` — ПОСЛЕ столкновения (до переноса). fn apply_bouzidi(&mut self) { @@ -192,6 +275,109 @@ impl Level { } } + /// ГРАНИЧНОЕ УСЛОВИЕ ГРАДА на теле (Dorschner и др., JFM 801 (2016), прил. B). + /// + /// В отличие от Bouzidi, работающего полинково, здесь условие ставится сразу на весь + /// граничный узел и не на популяции, а на моменты: + /// + /// целевая скорость u_tgt = (1/n) Σ (q_i·u_f,i + u_w,i)/(1 + q_i) (B 1) + /// целевая плотность ρ_tgt = Σ_известные f_i + Σ_недостающие f_i^отскок (B 3) + /// тензор давлений Π = ρc_s²I + ρu⊗u − (ρc_s²/2β)(∇u + ∇uᵀ) (2.14)–(2.16) + /// + /// после чего недостающие популяции собираются приближением Града (2.13). Скорости + /// соседей и градиенты берутся с прошлого шага — см. `u_prev`. + fn apply_grad_wall(&mut self) { + let nx = self.nx; + // стенка неподвижна; для подвижного тела сюда пойдёт её скорость на линке, + // и добавится динамическая часть плотности (B 4) + let (uwx, uwy) = (0.0, 0.0); + for wi in 0..self.geom.wall_nodes.len() { + let w = self.geom.wall_nodes[wi]; + let node = w.node as usize; + let lo = w.first as usize; + let hi = lo + w.count as usize; + + // ── целевая скорость (B 1) ── + let (mut ux, mut uy) = (0.0, 0.0); + for l in &self.geom.links[lo..hi] { + let far = l.far as usize; + let (fx, fy) = if self.geom.solid[far] { + // дальнего жидкого соседа нет (тело тоньше двух клеток) — остаётся стенка + (uwx, uwy) + } else { + self.u_at_t(far) + }; + ux += math::grad_target_velocity_term(l.q, fx, uwx); + uy += math::grad_target_velocity_term(l.q, fy, uwy); + } + let inv = 1.0 / w.count as R; + ux *= inv; + uy *= inv; + + // ── целевая плотность (B 3): известные популяции плюс отскок недостающих ── + let mut missing = [false; Q]; + for l in &self.geom.links[lo..hi] { + missing[l.ib as usize] = true; + } + let mut rho = 0.0; + for i in 0..Q { + rho += if missing[i] { self.post[node][OPP[i]] } else { self.f[node][i] }; + } + + let (dudx, dudy, dvdx, dvdy) = self.grad_u_at_t(node); + let g = math::grad_wall(rho, ux, uy, dudx, dudy, dvdx, dvdy, self.beta[node % nx]); + for i in 0..Q { + if missing[i] { + self.f[node][i] = g[i]; + } + } + } + } + + /// Скорость узла на момент t — берётся из ПОСТ-СТОЛКНОВИТЕЛЬНОГО поля. + /// + /// Столкновение сохраняет ρ и ρu точно, поэтому macros(post) даёт ровно ту же скорость, + /// что было в f до переноса, то есть u(x, t) — именно то, что требует прил. B. Отдельное + /// хранилище «поля предыдущего шага» при этом не нужно. + #[inline] + fn u_at_t(&self, node: usize) -> (R, R) { + let (_, a, b) = math::macros(&self.post[node]); + (a, b) + } + + /// Градиент скорости на момент t: центральная разность там, где оба соседа жидкие, + /// односторонняя — где один твёрдый, ноль — если твёрдые оба. + /// Возвращает (∂u_x/∂x, ∂u_x/∂y, ∂u_y/∂x, ∂u_y/∂y). + fn grad_u_at_t(&self, node: usize) -> (R, R, R, R) { + let (nx, ny) = (self.nx, self.ny); + let (x, y) = (node % nx, node / nx); + let xm = y * nx + (x + nx - 1) % nx; + let xp = y * nx + (x + 1) % nx; + let ym = ((y + ny - 1) % ny) * nx + x; + let yp = ((y + 1) % ny) * nx + x; + let solid = &self.geom.solid; + let uc = self.u_at_t(node); + let comp = |t: (R, R), c: usize| if c == 0 { t.0 } else { t.1 }; + let d = |a: usize, b: usize, c: usize| -> R { + match (!solid[a], !solid[b]) { + (true, true) => 0.5 * (comp(self.u_at_t(b), c) - comp(self.u_at_t(a), c)), + (false, true) => comp(self.u_at_t(b), c) - comp(uc, c), + (true, false) => comp(uc, c) - comp(self.u_at_t(a), c), + (false, false) => 0.0, + } + }; + (d(xm, xp, 0), d(ym, yp, 0), d(xm, xp, 1), d(ym, yp, 1)) + } + + /// Замкнуть недостающие популяции выбранной моделью стенки. + fn apply_wall(&mut self, model: WallModel) { + match model { + WallModel::Bouzidi => self.apply_bouzidi(), + WallModel::Grad => self.apply_grad_wall(), + WallModel::Staircase => self.apply_staircase(), + } + } + /// FREE-SLIP (зеркальные) стенки канала сверху и снизу: касательный импульс сохраняется, /// нормальный заворачивается, масса сохраняется. Применять ПОСЛЕ переноса и ПЕРЕД Zou–He. fn free_slip_walls(&mut self) { @@ -635,9 +821,10 @@ impl Sim { // ── уровень 0 ── self.pre.copy_from_slice(&self.l0.f); let model = self.spec.kbc_model; + let wall = self.spec.wall; let stats0 = self.l0.collide(sp_collision, model); self.l0.stream(); - self.l0.apply_bouzidi(); + self.l0.apply_wall(wall); self.l0.free_slip_walls(); self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap); @@ -649,7 +836,7 @@ impl Sim { for s in 0..p.r { l1.collide(sp_collision, model); l1.stream(); - l1.apply_bouzidi(); + l1.apply_wall(wall); // силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на // последнем подшаге даёт лишний шум в рядах при том же среднем let (a, b, c) = l1.force(); @@ -782,6 +969,7 @@ mod tests { ny, solid: vec![false; nx * ny], links: Vec::new(), + wall_nodes: Vec::new(), body_cx: 0.0, body_cy: 0.0, }, @@ -843,6 +1031,44 @@ mod tests { } } + /// СУБСЕТОЧНОСТЬ ГРАНИЦЫ. Стенка обязана стоять не на узлах решётки: при сдвиге тела на + /// долю клетки доли пересечения q обязаны поехать, а у разрешённого тела не должно быть + /// ни одного отката на простой (ступенчатый) отскок. Если бы граница была ступенчатой, + /// средняя q не реагировала бы на смещение вовсе. + #[test] + fn wall_is_subgrid_not_staircase() { + let mut means = Vec::new(); + for off in [0.0, 0.25, 0.5] { + let body = Body::new(math::ShapeKind::Cylinder, 60.0, 40.0 + off, 24.0, 0.0, 0.3); + let g = Geom::build(120, 80, &body); + let w = g.wall_stats(); + assert!(w.links > 100, "смещение {off}: линков всего {}", w.links); + assert_eq!(w.simple, 0, "смещение {off}: {} ступенчатых линков", w.simple); + assert!(w.q_min > 0.0 && w.q_max < 1.0, "смещение {off}: q вышла за (0,1)"); + means.push(w.q_mean); + } + let hi = means.iter().cloned().fold(R::NEG_INFINITY, R::max); + let lo = means.iter().cloned().fold(R::INFINITY, R::min); + assert!(hi - lo > 0.02, "средняя q не отреагировала на сдвиг тела: {means:?}"); + } + + /// Индекс граничных узлов обязан в точности разбивать список линков: без пропусков, + /// без пересечений, в том же порядке. + #[test] + fn wall_node_index_partitions_links() { + let body = Body::new(math::ShapeKind::Naca, 60.0, 40.0, 30.0, 12.0, 0.25); + let g = Geom::build(120, 80, &body); + let mut covered = 0usize; + for w in &g.wall_nodes { + assert_eq!(w.first as usize, covered, "разрыв в индексе граничных узлов"); + for l in &g.links[covered..covered + w.count as usize] { + assert_eq!(l.node, w.node, "линк чужого узла в диапазоне"); + } + covered += w.count as usize; + } + assert_eq!(covered, g.links.len(), "индекс покрыл не все линки"); + } + /// Однородное равновесие — неподвижная точка схемы: не должно никуда уехать. #[test] fn uniform_equilibrium_is_fixed_point() { diff --git a/docs/theory/2d_solver/src/gpu.rs b/docs/theory/2d_solver/src/gpu.rs index 2d6fc6b..7063dc1 100644 --- a/docs/theory/2d_solver/src/gpu.rs +++ b/docs/theory/2d_solver/src/gpu.rs @@ -784,6 +784,13 @@ impl Sim { )) .map_err(|e| format!("не удалось получить устройство: {e}"))?; + // Условие Града на GPU пока не перенесено. Молча считать по Bouzidi нельзя — это + // была бы другая физика под тем же ключом, поэтому отказываемся явно. + if spec.wall == math::WallModel::Grad { + return Err("--wall grad на GPU не реализован (перенесён только Bouzidi). Запусти с --backend cpu либо оставь --wall bouzidi" + .into()); + } + // ── топология берётся из процессорного бэкенда, а не строится заново ── let geom0 = cpu::Geom::build(spec.nx, spec.ny, &spec.body); let fluid_count = geom0.solid.iter().filter(|s| !**s).count() as R; diff --git a/docs/theory/2d_solver/src/main.rs b/docs/theory/2d_solver/src/main.rs index 0dccff5..d4072f3 100644 --- a/docs/theory/2d_solver/src/main.rs +++ b/docs/theory/2d_solver/src/main.rs @@ -95,6 +95,8 @@ pub struct Spec { pub collision: Collision, /// Что входит в сдвиговую часть s (см. `math::KbcModel`). pub kbc_model: math::KbcModel, + /// Чем замыкаются недостающие популяции на теле (см. `math::WallModel`). + pub wall: math::WallModel, /// Узел зонда следа в координатах L0. pub probe: (usize, usize), } @@ -167,9 +169,14 @@ struct Cli { body_y: Option, // ── время ── - /// Число шагов симуляции - #[arg(long, default_value_t = 20000, help_heading = "Время")] - steps: u64, + /// Число шагов симуляции. Взаимоисключающе с --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, @@ -190,6 +197,13 @@ struct Cli { /// Оператор столкновения #[arg(long, default_value = "kbc", value_parser = ["kbc", "bgk"], help_heading = "Схема")] collision: String, + /// Модель стенки на теле: bouzidi — интерполированный отскок по линкам; grad — условие + /// Града (Dorschner и др., JFM 801 (2016)), где задаются не популяции, а целевые моменты + /// ρ, u и тензор давлений. Авторы метода предпочитают grad: интерполяционные схемы, по их + /// словам, «ограничены низкими Re, поскольку на границе возникают паразитные скачки». + #[arg(long, default_value = "bouzidi", value_parser = math::WallModel::ALL, + help_heading = "Схема")] + wall: String, /// Состав сдвиговой части KBC (табл. I 2D-статьи): n1 — только девиатор {N, Π_xy} /// (KBC D), объёмная вязкость гуляет вместе с γ и может стать отрицательной; n2 — /// девиатор со следом {N, Π_xy, T} (KBC C), объёмная вязкость фиксирована ξ = ν. @@ -375,6 +389,33 @@ fn build_spec(cli: &Cli) -> Result { (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, + }; + Ok(Spec { nx: cli.nx, ny: cli.ny, @@ -383,7 +424,7 @@ fn build_spec(cli: &Cli) -> Result { patch, units, beta0, - steps: cli.steps, + steps, ramp: cli.ramp.max(1), flow_angle: cli.flow_angle * std::f64::consts::PI / 180.0, pert_amp: cli.pert_amp, @@ -394,6 +435,7 @@ fn build_spec(cli: &Cli) -> Result { sponge_mult: cli.sponge_mult, 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, }) } @@ -666,6 +708,29 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) { } else { "LBGK" }); + // Действительно ли субсеточная граница работает: если большинство линков сваливается в + // простой отскок, стенка де-факто ступенчатая, как её ни называй. + { + let g = cpu::Geom::build(spec.nx, spec.ny, &spec.body); + 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::Staircase => "простой отскок", + }); println!(" выход по u_y {:>12} губка {} столбцов ×{:.0}", if spec.outlet_extrapolate { "extrapolate" } else { "zero" }, spec.sponge_len, spec.sponge_mult); diff --git a/docs/theory/2d_solver/src/math.rs b/docs/theory/2d_solver/src/math.rs index d95364f..35718bf 100644 --- a/docs/theory/2d_solver/src/math.rs +++ b/docs/theory/2d_solver/src/math.rs @@ -97,10 +97,6 @@ pub fn feq(rho: R, ux: R, uy: R) -> [R; Q] { /// /// По построению сохраняет ρ и ρu точно (третий момент весов обнуляется по симметрии). /// -/// В самом решателе пока не используется: канальная постановка стартует с однородного потока, -/// где Π⁽¹⁾ = 0 и приближение Града вырождается в равновесие. Нужна для эталонных течений, -/// у которых стартовое поле имеет ненулевые градиенты. -#[cfg_attr(not(test), allow(dead_code))] #[inline] pub fn grad_init(rho: R, ux: R, uy: R, pxx: R, pxy: R, pyy: R) -> [R; Q] { let axx = pxx - rho * CS2; @@ -662,6 +658,86 @@ pub fn zou_he_outlet(f: &mut [R; Q], rho_out: R, uy_out: R) { f[7] = f[5] + d - (1.0 / 6.0) * rho_out * ux - 0.5 * rho_out * uy_out; } +/// Чем замыкаются недостающие популяции на теле — те, что после переноса должны были прийти +/// из твёрдого узла. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum WallModel { + /// Интерполированный отскок Bouzidi: формула применяется к каждому линку по отдельности, + /// стенка стоит на доле q вдоль линка. Второй порядок по положению стенки, дёшево. + Bouzidi, + /// Граничное условие Града (Dorschner, Bösch, Chikatamarla, Boulouchos, Karlin, + /// J. Fluid Mech. 801 (2016) 623–651, разд. 2.1 и прил. B). + /// + /// Условие ставится не на сами популяции, а на МОМЕНТЫ — целевые ρ, u и тензор давлений, — + /// после чего недостающие популяции собираются приближением Града (ур. 2.13). Авторы + /// метода предпочитают его интерполяционным схемам: те «ограничены низкими числами + /// Рейнольдса, поскольку на границе возникают паразитные скачки» (разд. 2.1). + /// Заодно это естественная форма для подвижных стенок: скорость стенки входит в целевые + /// значения, а не в отдельную поправку. + /// + /// ВАЖНО про субсеточность: положение стенки входит сюда ТОЛЬКО через целевую скорость + /// (B 1) — одну усреднённую по узлу величину. Целевая плотность (B 3) и тензор давлений + /// строятся вообще без q, а все недостающие популяции узла собираются из одного набора + /// моментов, то есть полинковая направленность теряется. Bouzidi, наоборот, применяет + /// свою q к каждой популяции отдельно. Поэтому геометрия у Града разрешена грубее — это + /// видно и по чувствительности к положению тела внутри клетки, и по точности на грубых + /// сетках (см. README). + Grad, + /// ПРОСТОЙ ОТСКОК: q игнорируется, стенка по построению лежит ровно посередине между + /// узлами — ступенчатая поверхность. Первый порядок. Держится не для счёта, а как база + /// сравнения: показывает, сколько именно даёт субсеточность. + Staircase, +} + +impl WallModel { + pub fn from_str(s: &str) -> Option { + Some(match s.to_ascii_lowercase().as_str() { + "bouzidi" | "bb" => WallModel::Bouzidi, + "grad" => WallModel::Grad, + "staircase" | "step" => WallModel::Staircase, + _ => return None, + }) + } + pub const ALL: [&'static str; 3] = ["bouzidi", "grad", "staircase"]; +} + +/// Недостающие популяции по целевым моментам — ур. (2.13) вместе с (2.14)–(2.16): +/// +/// Π = Π^eq + Π^neq, Π^eq = ρc_s²I + ρu⊗u, Π^neq = −(ρc_s²/2β)(∇u + ∇uᵀ) +/// +/// Возвращает ПОЛНЫЙ набор из девяти популяций; вызывающий берёт из него только недостающие. +/// Градиенты берутся с поля предыдущего шага: в момент применения условия скорость в самом +/// граничном узле ещё не определена — недостающие популяции как раз и вычисляются. +#[allow(clippy::too_many_arguments)] +#[inline] +pub fn grad_wall( + rho: R, + ux: R, + uy: R, + dudx: R, + dudy: R, + dvdx: R, + dvdy: R, + beta: R, +) -> [R; Q] { + let pref = rho * CS2 / (2.0 * beta); + let pxx = rho * CS2 + rho * ux * ux - pref * 2.0 * dudx; + let pyy = rho * CS2 + rho * uy * uy - pref * 2.0 * dvdy; + let pxy = rho * ux * uy - pref * (dudy + dvdx); + grad_init(rho, ux, uy, pxx, pxy, pyy) +} + +/// Целевая скорость граничного узла — ур. (B 1): линейная интерполяция между стенкой и +/// «дальним» жидким соседом по каждому недостающему направлению, усреднённая по ним. +/// +/// u_tgt = (1/n) Σ (q_i·u_f,i + u_w,i)/(1 + q_i) +/// +/// Три точки лежат на одной прямой: стенка на −q от узла, сам узел в нуле, дальний сосед на +1. +#[inline] +pub fn grad_target_velocity_term(q: R, u_far: R, u_wall: R) -> R { + (q * u_far + u_wall) / (1.0 + q) +} + // ───────────────────────────────────────────────────────────────────────────── // Сила на теле // ───────────────────────────────────────────────────────────────────────────── @@ -1025,6 +1101,52 @@ mod tests { assert!((rho - 1.0).abs() < 1e-12, "ρ={rho}"); } + /// Приближение Града обязано ТОЧНО воспроизводить моменты, из которых собрано: + /// плотность, импульс и полный тензор давлений. На этом и держится идея граничного + /// условия — задаём моменты, а не популяции. + #[test] + fn grad_reproduces_target_moments() { + let (rho, ux, uy) = (1.02, 0.03, -0.01); + let (pxx, pxy, pyy) = (0.35, 0.004, 0.33); + let g = grad_init(rho, ux, uy, pxx, pxy, pyy); + let (r, x, y) = macros(&g); + assert!((r - rho).abs() < 1e-13, "ρ: {r}"); + assert!((x - ux).abs() < 1e-13, "ux: {x}"); + assert!((y - uy).abs() < 1e-13, "uy: {y}"); + let mxx: R = (0..Q).map(|i| (CX[i] * CX[i]) as R * g[i]).sum(); + let myy: R = (0..Q).map(|i| (CY[i] * CY[i]) as R * g[i]).sum(); + let mxy: R = (0..Q).map(|i| (CX[i] * CY[i]) as R * g[i]).sum(); + assert!((mxx - pxx).abs() < 1e-13, "Π_xx: {mxx} vs {pxx}"); + assert!((myy - pyy).abs() < 1e-13, "Π_yy: {myy} vs {pyy}"); + assert!((mxy - pxy).abs() < 1e-13, "Π_xy: {mxy} vs {pxy}"); + } + + /// Граничная сборка (2.13)+(2.14)–(2.16): при нулевых градиентах неравновесной части нет, + /// и результат обязан совпасть с равновесием по тем же ρ, u с точностью до порядка Ma³ + /// (product-form равновесие не полиномиально, поэтому не побитово). + #[test] + fn grad_wall_without_gradients_is_near_equilibrium() { + let (rho, ux, uy) = (1.0, 0.02, 0.01); + let g = grad_wall(rho, ux, uy, 0.0, 0.0, 0.0, 0.0, 0.9); + let fe = feq(rho, ux, uy); + for i in 0..Q { + assert!((g[i] - fe[i]).abs() < 1e-4, "i={i}: {} vs {}", g[i], fe[i]); + } + } + + /// Целевая скорость (B 1) — линейная интерполяция между стенкой и дальним соседом. + #[test] + fn grad_target_velocity_interpolates() { + // стенка вплотную к узлу (q→0) ⇒ берём скорость стенки + assert!((grad_target_velocity_term(0.0, 1.0, 0.0) - 0.0).abs() < 1e-15); + // стенка на целую ячейку (q=1) ⇒ ровно середина между стенкой и соседом + assert!((grad_target_velocity_term(1.0, 1.0, 0.0) - 0.5).abs() < 1e-15); + // неподвижная стенка, покоящийся сосед ⇒ ноль при любом q + for q in [0.1, 0.5, 0.9] { + assert!(grad_target_velocity_term(q, 0.0, 0.0).abs() < 1e-15); + } + } + /// SDF: знак внутри/снаружи и |∇φ| ≈ 1 у поверхности для каждой формы. #[test] fn sdf_sign_and_gradient() {