Сверка со статьями: эталон Тейлора–Грина, выбор модели KBC, лечение продольной акустики
Проверка переориентирована на первоисточники из docs/origins вместо сравнения с
python-прототипом.
ЭТАЛОН ТЕЙЛОРА–ГРИНА (разд. VI 2D-статьи) — единственное течение статьи с точным
аналитическим решением. Второй порядок сходимости воспроизведён: 2.05 и 2.01 на
сетках 64/128/256, отдельно по амплитуде и по форме. Два места потребовали разбора:
* начальное давление. Течение несёт собственное поле p ~ ρu₀², находимое из
∇²p = 2ρ(ψ_xx·ψ_yy − ψ_xy²). Старт с ρ ≡ 1 сбрасывает разницу в акустику, которая
в периодическом ящике не затухает и садится полкой на ошибку. Работа 2024 года
делает то же самое явно, решая ∂ρ/∂t + ∇·(ρu₀) = D∇²ρ до стационара;
* способ измельчения. При фиксированном u₀ ошибка упирается в полку O(Ma²) (измерено:
относительная ошибка формы ∝ u₀^1.07). Порядок виден целиком только при диффузионном
измельчении, ν = const и u₀ ∝ 1/N. Это свойство слабо-сжимаемого метода, а не
реализации: LBGK на том же тесте даёт ту же полку, что согласуется с утверждением
статьи о практически одинаковом поведении всех моделей.
Добавлено приближение Града (ур. 58) для согласованного старта эталонов.
СВЕРКА ОПЕРАТОРА с пошаговым листингом работы 2024 года — совпадает дословно, включая
Δh = f − f^eq − Δs и f′ = f − β(2Δs + γΔh). Та же работа подтверждает относительный
порог вырожденности γ: стабилизатор «далеко не постоянен», а MRT с γ = const не
достигает той же устойчивости.
ВЫБОР СОСТАВА СДВИГОВОЙ ЧАСТИ (табл. I) — ключ --kbc-model: n1 = {N, Π_xy} (KBC D),
n2 = {N, Π_xy, T} (KBC C). По точности неразличимы, как и заявляет статья. Разница в
объёмной вязкости (ур. 57): у n2 она фиксирована ξ = ν, у n1 равна c_s²(1/(γβ) − ½) и
при измеренной ⟨γ⟩ ≈ 1.73 < 2 в среднем вчетверо больше ν. Поэтому вопреки ожиданию
именно n1 сильнее демпфирует продольную акустику и оставлен умолчанием.
ПРОДОЛЬНАЯ АКУСТИКА КАНАЛА — найдена причина «поршневого» поведения потока. Вход по
скорости акустически есть жёсткий поршень, выход по давлению — открытый конец; канал
работает четвертьволновым резонатором. Измерено: период пульсации 3332 шага против
расчётных 4·Nx/c_s = 3325, первый ноль автокорреляции ровно на Nx/c_s, затухание за
30000 шагов — 3.8%, то есть мода не гаснет. Возбуждал её сам старт из покоя.
Исправлено умолчаниями: --init uniform (домен сразу заполнен потоком) и автоподбор
губки по длине домена. Пульсация упала с 75% до 1.0% от скорости потока, а числа
выправились сами: St 0.183 против прежних 0.194 (литература 0.183), rms Cl 0.38
против 0.84, ⟨ρ⟩ = 1.0000. Отчёт теперь печатает период моды, время её вязкого
затухания и измеренную пульсацию с предупреждением при превышении 5%.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5
@@ -113,8 +113,11 @@ pub struct Level {
|
||||
}
|
||||
|
||||
impl Level {
|
||||
fn new(nx: usize, ny: usize, beta: Vec<R>, geom: Geom) -> Self {
|
||||
let f0 = math::feq(1.0, 0.0, 0.0);
|
||||
/// `u0` — скорость, которой заполняется поле на старте. Заполнять сразу набегающим
|
||||
/// потоком принципиально: старт из покоя разгоняет весь столб жидкости и закачивает в
|
||||
/// канал продольную акустическую моду, которую вязкость потом почти не гасит.
|
||||
fn new(nx: usize, ny: usize, beta: Vec<R>, 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 }
|
||||
}
|
||||
@@ -122,7 +125,7 @@ impl Level {
|
||||
/// Столкновение по всем жидким узлам. Внутри тела не считаем: эти популяции фиктивны
|
||||
/// (ГУ Bouzidi перезаписывает всё, что могло бы прийти из тела в жидкость), а счёт там
|
||||
/// только жжёт такты и способен родить NaN при экстремальных режимах.
|
||||
fn collide(&mut self, op: Collision) -> Stats {
|
||||
fn collide(&mut self, op: Collision, model: math::KbcModel) -> Stats {
|
||||
let nx = self.nx;
|
||||
let beta = &self.beta;
|
||||
let solid = &self.geom.solid;
|
||||
@@ -141,11 +144,11 @@ impl Level {
|
||||
}
|
||||
let b = beta[node % nx];
|
||||
let k: Kbc = match op {
|
||||
Collision::Kbc => math::collide_node(p, b),
|
||||
Collision::Kbc => math::collide_node(p, b, model),
|
||||
Collision::Bgk => math::collide_node_bgk(p, b),
|
||||
};
|
||||
*g = k.gamma;
|
||||
Stats::of(k, b)
|
||||
Stats::of(k, b, model)
|
||||
})
|
||||
.reduce(|| Stats::EMPTY, Stats::merge)
|
||||
}
|
||||
@@ -288,14 +291,14 @@ impl Stats {
|
||||
xi_negative: 0,
|
||||
};
|
||||
#[inline]
|
||||
fn of(k: Kbc, beta: R) -> Stats {
|
||||
fn of(k: Kbc, beta: R, model: math::KbcModel) -> Stats {
|
||||
Stats {
|
||||
n: 1,
|
||||
gsum: k.gamma,
|
||||
gmin: k.gamma,
|
||||
gmax: k.gamma,
|
||||
degenerate: k.degenerate as u64,
|
||||
xi_negative: (math::xi_of_gamma(k.gamma, beta) < 0.0) as u64,
|
||||
xi_negative: (math::xi_of_gamma(k.gamma, beta, model) < 0.0) as u64,
|
||||
}
|
||||
}
|
||||
fn merge(a: Stats, b: Stats) -> Stats {
|
||||
@@ -545,7 +548,14 @@ impl Sim {
|
||||
})
|
||||
.collect();
|
||||
|
||||
let l0 = Level::new(nx, ny, beta0, geom0);
|
||||
// стартовое поле: либо сразу набегающий поток, либо покой
|
||||
let u0 = if spec.init_uniform {
|
||||
let (sn, cs) = spec.flow_angle.sin_cos();
|
||||
(spec.units.u_lat * cs, spec.units.u_lat * sn)
|
||||
} else {
|
||||
(0.0, 0.0)
|
||||
};
|
||||
let l0 = Level::new(nx, ny, beta0, geom0, u0);
|
||||
let fluid_count = l0.geom.solid.iter().filter(|s| !**s).count() as R;
|
||||
|
||||
let (l1, patch) = if spec.refine > 1 {
|
||||
@@ -561,7 +571,7 @@ 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);
|
||||
(Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1)), Some(patch))
|
||||
(Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1, u0)), Some(patch))
|
||||
} else {
|
||||
(None, None)
|
||||
};
|
||||
@@ -594,12 +604,19 @@ impl Sim {
|
||||
|
||||
/// Скорость на входе в момент t: разгон smoothstep плюс окно поперечного возмущения.
|
||||
///
|
||||
/// Разгон гасит импульсный старт (мгновенное включение входа шлёт по домену ударную волну).
|
||||
/// Возмущение — короткий поперечный импульс, сбивающий симметрию: без него дорожка Кармана
|
||||
/// заводится только на численном шуме и стартует на порядок позже.
|
||||
/// При старте из однородного потока разгон не нужен и пропускается: поле и вход и так
|
||||
/// согласованы, а плавный разгон поверх согласованного поля сам стал бы рассогласованием.
|
||||
/// При старте из покоя разгон обязателен — мгновенное включение входа шлёт по домену
|
||||
/// ударную волну. Возмущение — короткий поперечный импульс, сбивающий симметрию: без него
|
||||
/// дорожка Кармана заводится только на численном шуме и стартует на порядок позже.
|
||||
fn inlet(&self, t: u64) -> (R, R) {
|
||||
let sp = &self.spec;
|
||||
let u = sp.units.u_lat * math::smoothstep(t as R / sp.ramp.max(1) as R);
|
||||
let ramp = if sp.init_uniform {
|
||||
1.0
|
||||
} else {
|
||||
math::smoothstep(t as R / sp.ramp.max(1) as R)
|
||||
};
|
||||
let u = sp.units.u_lat * ramp;
|
||||
let (s, c) = sp.flow_angle.sin_cos();
|
||||
let mut uy = u * s;
|
||||
if sp.pert_dur > 0 && t >= sp.ramp && t < sp.ramp + sp.pert_dur {
|
||||
@@ -617,7 +634,8 @@ impl Sim {
|
||||
|
||||
// ── уровень 0 ──
|
||||
self.pre.copy_from_slice(&self.l0.f);
|
||||
let stats0 = self.l0.collide(sp_collision);
|
||||
let model = self.spec.kbc_model;
|
||||
let stats0 = self.l0.collide(sp_collision, model);
|
||||
self.l0.stream();
|
||||
self.l0.apply_bouzidi();
|
||||
self.l0.free_slip_walls();
|
||||
@@ -629,7 +647,7 @@ impl Sim {
|
||||
p.ghost_values(&self.pre, &mut self.gh_old);
|
||||
p.ghost_values(&self.l0.f, &mut self.gh_new);
|
||||
for s in 0..p.r {
|
||||
l1.collide(sp_collision);
|
||||
l1.collide(sp_collision, model);
|
||||
l1.stream();
|
||||
l1.apply_bouzidi();
|
||||
// силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на
|
||||
@@ -767,15 +785,30 @@ mod tests {
|
||||
body_cx: 0.0,
|
||||
body_cy: 0.0,
|
||||
},
|
||||
(0.0, 0.0),
|
||||
)
|
||||
}
|
||||
|
||||
/// Голый периодический шаг: столкновение + перенос, без единого ГУ.
|
||||
/// Нужен, чтобы отделить ядро схемы от граничных условий.
|
||||
fn periodic_step(f: &mut [[R; Q]], tmp: &mut [[R; Q]], nx: usize, ny: usize, beta: R) {
|
||||
periodic_step_op(f, tmp, nx, ny, beta, Collision::Kbc)
|
||||
}
|
||||
|
||||
fn periodic_step_op(
|
||||
f: &mut [[R; Q]],
|
||||
tmp: &mut [[R; Q]],
|
||||
nx: usize,
|
||||
ny: usize,
|
||||
beta: R,
|
||||
op: Collision,
|
||||
) {
|
||||
for (n, c) in f.iter().enumerate() {
|
||||
tmp[n] = *c;
|
||||
math::collide_node(&mut tmp[n], beta);
|
||||
match op {
|
||||
Collision::Kbc => math::collide_node(&mut tmp[n], beta, math::KbcModel::N1),
|
||||
Collision::Bgk => math::collide_node_bgk(&mut tmp[n], beta),
|
||||
};
|
||||
}
|
||||
for y in 0..ny {
|
||||
for x in 0..nx {
|
||||
@@ -865,6 +898,183 @@ mod tests {
|
||||
}
|
||||
}
|
||||
|
||||
/// Диагностика: из чего складывается ошибка Тейлора–Грина. Сетка и вязкость ЗАФИКСИРОВАНЫ,
|
||||
/// меняется только амплитуда u₀. Дискретизационная ошибка от u₀ почти не зависит (задача
|
||||
/// в этом пределе линейна), а сжимаемостная идёт как Ma² ∝ u₀². Наклон и покажет, что
|
||||
/// доминирует.
|
||||
#[test]
|
||||
#[ignore = "диагностика; запуск: cargo test --release -- --ignored taylor_green_error_scaling --nocapture"]
|
||||
fn taylor_green_error_scaling() {
|
||||
let n = 128usize;
|
||||
let nu = 0.0384;
|
||||
let beta = math::beta_of_nu(nu);
|
||||
println!(" N={n}, ν={nu} (τ={:.4}) — меняем только u₀", 1.0 / (2.0 * beta));
|
||||
let mut prev: Option<(R, R)> = None;
|
||||
for &u0 in &[0.04, 0.02, 0.01, 0.005] {
|
||||
let e = taylor_green_error(n, u0, nu, beta);
|
||||
let slope = prev.map(|(p_u, p_e): (R, R)| (p_e / e.shape).ln() / (p_u / u0).ln());
|
||||
match slope {
|
||||
Some(s) => println!(
|
||||
" u₀={u0:<7} Ma={:.4} форма {:.3e} ампл {:.3e} наклон формы по u₀: {s:.2}",
|
||||
u0 / math::CS2.sqrt(), e.shape, e.amplitude),
|
||||
None => println!(" u₀={u0:<7} Ma={:.4} форма {:.3e} ампл {:.3e}",
|
||||
u0 / math::CS2.sqrt(), e.shape, e.amplitude),
|
||||
}
|
||||
prev = Some((u0, e.shape));
|
||||
}
|
||||
println!(" наклон ≈2 ⇒ правит сжимаемость (Ma²); ≈0 ⇒ правит дискретизация");
|
||||
}
|
||||
|
||||
/// Один прогон вихря Тейлора–Грина до полураспада; возвращает Σ|u_x−точн|/Σ|точн|.
|
||||
///
|
||||
/// Начальное состояние ставится ПОЛНОСТЬЮ согласованным: скорость, давление и неравновесная
|
||||
/// часть. Давление здесь не константа — течение несёт собственное поле порядка ρu₀²,
|
||||
/// которое находится из ∇²p = 2ρ(ψ_xx·ψ_yy − ψ_xy²):
|
||||
///
|
||||
/// p = −(ρu₀²/4)[cos(2k₁x) + (k₁²/k₂²)·cos(2k₂y)], δρ = p/c_s².
|
||||
///
|
||||
/// Если стартовать с ρ ≡ 1, эта разница уходит в акустику, которая в периодическом ящике
|
||||
/// почти не затухает и садится полкой на ошибку скорости, ломая порядок сходимости.
|
||||
fn taylor_green_error(n: usize, u0: R, nu: R, beta: R) -> TgError {
|
||||
taylor_green_error_op(n, u0, nu, beta, Collision::Kbc)
|
||||
}
|
||||
|
||||
fn taylor_green_error_op(n: usize, u0: R, nu: R, beta: R, op: Collision) -> TgError {
|
||||
let (k1, k2) = (1.0, 4.0);
|
||||
let tau = 1.0 / (2.0 * beta);
|
||||
let kk1 = 2.0 * std::f64::consts::PI * k1 / n as R;
|
||||
let kk2 = 2.0 * std::f64::consts::PI * k2 / n as R;
|
||||
let decay = nu * (kk1 * kk1 + kk2 * kk2);
|
||||
let tc = ((2.0_f64).ln() / decay).round() as usize;
|
||||
let exact = |x: usize, y: usize, e: R| -> (R, R) {
|
||||
let (a, b) = (kk1 * x as R, kk2 * y as R);
|
||||
(-u0 * a.cos() * b.sin() * e, (k1 / k2) * u0 * a.sin() * b.cos() * e)
|
||||
};
|
||||
let mut f = vec![[0.0; Q]; n * n];
|
||||
for y in 0..n {
|
||||
for x in 0..n {
|
||||
let (ux, uy) = exact(x, y, 1.0);
|
||||
let (a, b) = (kk1 * x as R, kk2 * y as R);
|
||||
let p = -(u0 * u0 / 4.0)
|
||||
* ((2.0 * a).cos() + (k1 * k1) / (k2 * k2) * (2.0 * b).cos());
|
||||
let rho = 1.0 + p / math::CS2;
|
||||
let dxux = u0 * kk1 * a.sin() * b.sin();
|
||||
let dyux = -u0 * kk2 * a.cos() * b.cos();
|
||||
let dxuy = (k1 / k2) * u0 * kk1 * a.cos() * b.cos();
|
||||
let dyuy = -(k1 / k2) * u0 * kk2 * a.sin() * b.sin();
|
||||
// Π = ρc_s²δ + ρuu + Π⁽¹⁾, Π⁽¹⁾ = −τρc_s²(∂_αu_β + ∂_βu_α) — ур. (53)
|
||||
let pref = -tau * rho * math::CS2;
|
||||
let pxx = rho * math::CS2 + rho * ux * ux + pref * 2.0 * dxux;
|
||||
let pyy = rho * math::CS2 + rho * uy * uy + pref * 2.0 * dyuy;
|
||||
let pxy = rho * ux * uy + pref * (dyux + dxuy);
|
||||
f[y * n + x] = math::grad_init(rho, ux, uy, pxx, pxy, pyy);
|
||||
}
|
||||
}
|
||||
let mut tmp = f.clone();
|
||||
for _ in 0..tc {
|
||||
periodic_step_op(&mut f, &mut tmp, n, n, beta, op);
|
||||
}
|
||||
let e = (-decay * tc as R).exp();
|
||||
let (mut num, mut den) = (0.0, 0.0);
|
||||
// заодно раскладываем ошибку: наилучшая подгонка амплитуды к точной форме
|
||||
let (mut dot, mut nrm2) = (0.0, 0.0);
|
||||
for y in 0..n {
|
||||
for x in 0..n {
|
||||
let got = math::macros(&f[y * n + x]).1;
|
||||
let want = exact(x, y, e).0;
|
||||
num += (got - want).abs();
|
||||
den += want.abs();
|
||||
dot += got * want;
|
||||
nrm2 += want * want;
|
||||
}
|
||||
}
|
||||
let amp = dot / nrm2; // 1.0 = амплитуда совпала
|
||||
let (mut snum, mut sden) = (0.0, 0.0);
|
||||
for y in 0..n {
|
||||
for x in 0..n {
|
||||
let got = math::macros(&f[y * n + x]).1;
|
||||
let want = amp * exact(x, y, e).0;
|
||||
snum += (got - want).abs();
|
||||
sden += want.abs();
|
||||
}
|
||||
}
|
||||
TgError { total: num / den, amplitude: (amp - 1.0).abs(), shape: snum / sden }
|
||||
}
|
||||
|
||||
/// Разложение ошибки эталона: полная, вклад амплитуды (скорость затухания) и вклад формы.
|
||||
struct TgError {
|
||||
total: R,
|
||||
amplitude: R,
|
||||
shape: R,
|
||||
}
|
||||
|
||||
/// ВИХРЬ ТЕЙЛОРА–ГРИНА — первый эталон 2D-статьи (разд. VI). Единственное из трёх течений
|
||||
/// статьи, у которого есть ТОЧНОЕ аналитическое решение, поэтому проверяется не «похоже на
|
||||
/// чужой прогон», а прямое совпадение с формулой и заявленный статьёй ВТОРОЙ ПОРЯДОК
|
||||
/// сходимости (рис. 1a: Re = 100, u₀ = 0.03, N ∈ {64, 128, 256}).
|
||||
///
|
||||
/// Постановка ровно по статье:
|
||||
/// u = ∇×[(u₀/k₂)cos(k₁x)cos(k₂y)·exp(−ν(k₁²+k₂²)t)], k₁ = 1, k₂ = 4,
|
||||
/// область 0 < x,y < 2π на сетке N×N, Re = u₀N/ν, полураспад t_c = ln2/[ν(k₁²+k₂²)].
|
||||
/// Волновые числа переводятся в решёточные: K = 2πk/N (узел — единица длины).
|
||||
/// Старт — приближением Града (ур. 58), как в статье.
|
||||
/// Метрика — та же, что на рис. 1: Σ|u_x − u_x^точн| / Σ|u_x^точн| в момент t_c.
|
||||
#[test]
|
||||
#[ignore = "долгий (до 256², t_c растёт как N²); запуск: cargo test --release -- --ignored"]
|
||||
fn taylor_green_converges_at_second_order() {
|
||||
// ДИФФУЗИОННОЕ ИЗМЕЛЬЧЕНИЕ: ν фиксирована, u₀ ∝ 1/N. Тогда Re = u₀N/ν сохраняется
|
||||
// (все три сетки считают ОДНО И ТО ЖЕ течение), а число Маха падает как 1/N — вместе
|
||||
// с ним падает и сжимаемостная ошибка метода. Только при таком измельчении второй
|
||||
// порядок виден целиком; при фиксированном u₀ ошибка упирается в полку O(Ma²), которая
|
||||
// от сетки не зависит вовсе (см. taylor_green_error_scaling).
|
||||
let nu = 0.0192;
|
||||
let beta = math::beta_of_nu(nu);
|
||||
let mut errs = Vec::new();
|
||||
for &n in &[64usize, 128, 256] {
|
||||
let u0 = 0.03 * 64.0 / n as R;
|
||||
let e = taylor_green_error(n, u0, nu, beta);
|
||||
println!(
|
||||
" N={n:>4} u₀={u0:.5} Re={:.0} полная {:.3e} амплитуда {:.3e} форма {:.3e}",
|
||||
u0 * n as R / nu, e.total, e.amplitude, e.shape
|
||||
);
|
||||
errs.push((n, e));
|
||||
}
|
||||
for w in errs.windows(2) {
|
||||
let (n0, e0) = (&w[0].0, &w[0].1);
|
||||
let (n1, e1) = (&w[1].0, &w[1].1);
|
||||
let ord = |a: R, b: R| (a / b).ln() / (*n1 as R / *n0 as R).ln();
|
||||
let (pt, pa, ps) = (
|
||||
ord(e0.total, e1.total),
|
||||
ord(e0.amplitude, e1.amplitude),
|
||||
ord(e0.shape, e1.shape),
|
||||
);
|
||||
println!(" порядок {n0}→{n1}: полная {pt:.2} амплитуда {pa:.2} форма {ps:.2}");
|
||||
assert!(
|
||||
(1.7..2.6).contains(&pt),
|
||||
"порядок {n0}→{n1} = {pt:.2}, статья (разд. VI) заявляет второй"
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// Диагностика: тот же вихрь Тейлора–Грина оператором LBGK. Статья (рис. 1) утверждает,
|
||||
/// что на этом течении LBGK и все варианты KBC идут практически одинаково; если полка по
|
||||
/// форме есть у обоих — она от схемы D2Q9, а не от энтропийного стабилизатора.
|
||||
#[test]
|
||||
#[ignore = "диагностика; запуск: cargo test --release -- --ignored taylor_green_kbc_vs_bgk --nocapture"]
|
||||
fn taylor_green_kbc_vs_bgk() {
|
||||
let (u0, re) = (0.03, 100.0);
|
||||
for &n in &[64usize, 128, 256] {
|
||||
let nu = u0 * n as R / re;
|
||||
let beta = math::beta_of_nu(nu);
|
||||
let k = taylor_green_error_op(n, u0, nu, beta, Collision::Kbc);
|
||||
let b = taylor_green_error_op(n, u0, nu, beta, Collision::Bgk);
|
||||
println!(
|
||||
" N={n:>4} KBC: форма {:.3e} ампл {:.3e} LBGK: форма {:.3e} ампл {:.3e}",
|
||||
k.shape, k.amplitude, b.shape, b.amplitude
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
/// СКВОЗНАЯ СВЕРКА С ЭТАЛОНОМ. Дважды периодический сдвиговый слой — один из трёх
|
||||
/// бенчмарков 2D-статьи (разд. VIII). Постановка: N=128, Re=30000, u0=0.04, κ=80, δ=0.05,
|
||||
/// одно конвективное время t_c = N/u0 шагов.
|
||||
|
||||
@@ -49,7 +49,8 @@ struct Dyn {
|
||||
ux_in: f32,
|
||||
uy_in: f32,
|
||||
rho_out: f32,
|
||||
_pad0: f32,
|
||||
/// 0 — сдвиговая часть только девиатор (N1), 1 — девиатор со следом (N2)
|
||||
kbc_model: u32,
|
||||
outlet_extrap: u32,
|
||||
nparts: u32,
|
||||
refine: u32,
|
||||
@@ -126,7 +127,7 @@ struct Results {
|
||||
const SHADER: &str = r#"
|
||||
// ───── структуры (обязаны совпадать с gpu.rs) ─────
|
||||
struct LevelParams { n:u32, nx:u32, ny:u32, nlinks:u32, bcx:f32, bcy:f32, probe_node:u32, flags:u32 };
|
||||
struct Dyn { ux_in:f32, uy_in:f32, rho_out:f32, p0:f32, outlet_extrap:u32, nparts:u32, refine:u32, collision:u32 };
|
||||
struct Dyn { ux_in:f32, uy_in:f32, rho_out:f32, kbc_model:u32, outlet_extrap:u32, nparts:u32, refine:u32, collision:u32 };
|
||||
struct Substep { idx:u32, p0:u32, p1:u32, p2:u32 };
|
||||
struct AmrParams { nghost:u32, nrestrict:u32, ccount:u32, fcount:u32, r01:f32, rfc:f32, w:f32, pad:f32 };
|
||||
struct GLink { node:u32, far:u32, i:u32, ib:u32, kind:u32, q:f32, p0:u32, p1:u32 };
|
||||
@@ -186,11 +187,18 @@ fn macros9(fv: array<f32,9>) -> vec3<f32> {
|
||||
return vec3<f32>(rho, mx/rho, my/rho);
|
||||
}
|
||||
|
||||
// Проекция неравновесия на сдвиговую часть модели KBC D: моменты N = M20-M02 и Pi_xy = M11
|
||||
fn shift9(d: array<f32,9>) -> array<f32,9> {
|
||||
// Проекция неравновесия на сдвиговую часть. model=0: {N, Pi_xy} (N1); model=1: плюс след T (N2)
|
||||
fn shift9(d: array<f32,9>, model: u32) -> array<f32,9> {
|
||||
let a = 0.25*(d[1] + d[3] - d[2] - d[4]);
|
||||
let b = 0.25*(d[5] + d[7] - d[6] - d[8]);
|
||||
return array<f32,9>(0.0, a, -a, a, -a, b, -b, b, -b);
|
||||
var s = array<f32,9>(0.0, a, -a, a, -a, b, -b, b, -b);
|
||||
if (model == 1u) {
|
||||
let dt = d[1] + d[2] + d[3] + d[4] + 2.0*(d[5] + d[6] + d[7] + d[8]);
|
||||
let q = 0.25*dt;
|
||||
s[0] = s[0] - dt;
|
||||
s[1] = s[1] + q; s[2] = s[2] + q; s[3] = s[3] + q; s[4] = s[4] + q;
|
||||
}
|
||||
return s;
|
||||
}
|
||||
|
||||
fn load9(base: ptr<function, array<f32,9>>, off: u32, n: u32) {
|
||||
@@ -200,7 +208,7 @@ fn load9(base: ptr<function, array<f32,9>>, off: u32, n: u32) {
|
||||
// Возвращает (пост-столкновительные популяции, gamma, признак вырождения)
|
||||
struct CollOut { fv: array<f32,9>, gamma: f32, degen: f32 };
|
||||
|
||||
fn collide9(fin: array<f32,9>, b: f32, op: u32) -> CollOut {
|
||||
fn collide9(fin: array<f32,9>, b: f32, op: u32, model: u32) -> CollOut {
|
||||
var out: CollOut;
|
||||
// WGSL разрешает переменный индекс только по памяти (var), а не по значению (let/параметр),
|
||||
// поэтому всё, что индексируется в цикле, кладётся в var
|
||||
@@ -215,7 +223,7 @@ fn collide9(fin: array<f32,9>, b: f32, op: u32) -> CollOut {
|
||||
}
|
||||
var d: array<f32,9>;
|
||||
for (var i = 0u; i < 9u; i = i + 1u) { d[i] = fv[i] - fe[i]; }
|
||||
var ds = shift9(d);
|
||||
var ds = shift9(d, model);
|
||||
var num = 0.0; var den = 0.0; var nrm = 0.0;
|
||||
for (var i = 0u; i < 9u; i = i + 1u) {
|
||||
let inv = 1.0 / fe[i];
|
||||
@@ -273,7 +281,7 @@ fn k_collide(@builtin(global_invocation_id) gid: vec3<u32>) {
|
||||
gam[nd] = 2.0;
|
||||
return;
|
||||
}
|
||||
var r = collide9(fv, beta[nd % P.nx], D.collision);
|
||||
var r = collide9(fv, beta[nd % P.nx], D.collision, D.kbc_model);
|
||||
for (var i = 0u; i < 9u; i = i + 1u) { post[i*P.n + nd] = r.fv[i]; }
|
||||
gam[nd] = r.gamma;
|
||||
}
|
||||
@@ -422,7 +430,9 @@ fn k_stats1(@builtin(global_invocation_id) gid: vec3<u32>,
|
||||
p.degen = select(0.0, 1.0, g == 2.0);
|
||||
// объёмная вязкость модели D: xi = cs^2 (1/(gamma*beta) - 1/2)
|
||||
let b = beta[nd % P.nx];
|
||||
let xi = CS2*(1.0/(g*b) - 0.5);
|
||||
// при следе в сдвиговой части (N2) объёмная вязкость равна сдвиговой и всегда > 0
|
||||
var xi = CS2*(1.0/(b + b) - 0.5);
|
||||
if (D.kbc_model == 0u) { xi = CS2*(1.0/(g*b) - 0.5); }
|
||||
p.xineg = select(0.0, 1.0, xi < 0.0);
|
||||
}
|
||||
// зонд следа снимается с того уровня, на котором он лежит
|
||||
@@ -596,8 +606,8 @@ impl GpuLevel {
|
||||
}
|
||||
}
|
||||
|
||||
fn soa_equilibrium(n: usize) -> Vec<f32> {
|
||||
let fe = math::feq(1.0, 0.0, 0.0);
|
||||
fn soa_equilibrium(n: usize, u0: (R, R)) -> Vec<f32> {
|
||||
let fe = math::feq(1.0, u0.0, u0.1);
|
||||
let mut v = vec![0.0f32; 9 * n];
|
||||
for i in 0..9 {
|
||||
for k in 0..n {
|
||||
@@ -617,9 +627,10 @@ fn make_level(
|
||||
beta: &[f32],
|
||||
probe_node: u32,
|
||||
flags: u32,
|
||||
u0: (R, R),
|
||||
) -> GpuLevel {
|
||||
let n = nx * ny;
|
||||
let init = soa_equilibrium(n);
|
||||
let init = soa_equilibrium(n, u0);
|
||||
let mkf = |label: &str| {
|
||||
device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
|
||||
label: Some(label),
|
||||
@@ -861,6 +872,14 @@ impl Sim {
|
||||
})
|
||||
.collect();
|
||||
|
||||
// стартовое поле: либо сразу набегающий поток, либо покой (см. cpu.rs)
|
||||
let u0 = if spec.init_uniform {
|
||||
let (sn, cs) = spec.flow_angle.sin_cos();
|
||||
(spec.units.u_lat * cs, spec.units.u_lat * sn)
|
||||
} else {
|
||||
(0.0, 0.0)
|
||||
};
|
||||
|
||||
// ── зонд ──
|
||||
let (px, py) = spec.probe;
|
||||
let probe_on_fine = matches!(spec.patch, Some((ax, bx, ay, by))
|
||||
@@ -879,6 +898,7 @@ impl Sim {
|
||||
&beta0,
|
||||
probe0 as u32,
|
||||
flags0,
|
||||
u0,
|
||||
);
|
||||
|
||||
let pre = device.create_buffer(&wgpu::BufferDescriptor {
|
||||
@@ -907,7 +927,7 @@ 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);
|
||||
make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1, u0);
|
||||
|
||||
let ghosts: Vec<GGhost> = patch
|
||||
.ghosts()
|
||||
@@ -1053,7 +1073,12 @@ impl Sim {
|
||||
|
||||
fn inlet(&self, t: u64) -> (R, R) {
|
||||
let sp = &self.spec;
|
||||
let u = sp.units.u_lat * math::smoothstep(t as R / sp.ramp.max(1) as R);
|
||||
let ramp = if sp.init_uniform {
|
||||
1.0
|
||||
} else {
|
||||
math::smoothstep(t as R / sp.ramp.max(1) as R)
|
||||
};
|
||||
let u = sp.units.u_lat * ramp;
|
||||
let (s, c) = sp.flow_angle.sin_cos();
|
||||
let mut uy = u * s;
|
||||
if sp.pert_dur > 0 && t >= sp.ramp && t < sp.ramp + sp.pert_dur {
|
||||
@@ -1074,7 +1099,7 @@ impl Sim {
|
||||
ux_in: ux_in as f32,
|
||||
uy_in: uy_in as f32,
|
||||
rho_out: 1.0,
|
||||
_pad0: 0.0,
|
||||
kbc_model: (self.spec.kbc_model == math::KbcModel::N2) as u32,
|
||||
outlet_extrap: self.spec.outlet_extrapolate as u32,
|
||||
nparts: self.l0.parts_count,
|
||||
refine: refine as u32,
|
||||
|
||||
@@ -88,9 +88,13 @@ pub struct Spec {
|
||||
pub pert_amp: R,
|
||||
pub pert_dur: u64,
|
||||
pub outlet_extrapolate: bool,
|
||||
/// Начальное поле уже несёт набегающий поток (иначе — покой с разгоном входа).
|
||||
pub init_uniform: bool,
|
||||
pub sponge_len: usize,
|
||||
pub sponge_mult: R,
|
||||
pub collision: Collision,
|
||||
/// Что входит в сдвиговую часть s (см. `math::KbcModel`).
|
||||
pub kbc_model: math::KbcModel,
|
||||
/// Узел зонда следа в координатах L0.
|
||||
pub probe: (usize, usize),
|
||||
}
|
||||
@@ -175,20 +179,32 @@ struct Cli {
|
||||
/// Длительность стартового импульса, шагов
|
||||
#[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,
|
||||
|
||||
// ── численная схема ──
|
||||
/// Оператор столкновения
|
||||
#[arg(long, default_value = "kbc", value_parser = ["kbc", "bgk"], help_heading = "Схема")]
|
||||
collision: 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, default_value_t = 0, help_heading = "Схема")]
|
||||
sponge_len: usize,
|
||||
/// Длина поглощающей губки перед выходом, столбцов. Не задана — подбирается по домену;
|
||||
/// 0 — выключить. Губка гасит продольную моду, которую вязкость сама не гасит.
|
||||
#[arg(long, help_heading = "Схема")]
|
||||
sponge_len: Option<usize>,
|
||||
/// Во сколько раз губка поднимает вязкость
|
||||
#[arg(long, default_value_t = 30.0, help_heading = "Схема")]
|
||||
sponge_mult: R,
|
||||
@@ -322,17 +338,38 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
if !(2 <= ax && ax < bx && bx <= cli.nx - 2) {
|
||||
return Err(format!("патч по x [{ax},{bx}] обязан лежать внутри (1,{})", cli.nx - 1));
|
||||
}
|
||||
if cli.sponge_len > 0 && cli.nx - 1 - cli.sponge_len <= bx {
|
||||
return Err(format!(
|
||||
"губка (последние {} столбцов) накрывает патч (bx={bx})",
|
||||
cli.sponge_len
|
||||
));
|
||||
}
|
||||
Some((ax, bx, ay, by))
|
||||
} else {
|
||||
None
|
||||
};
|
||||
|
||||
// Губка перед выходом. Задача канала «вход по скорости + выход по давлению» акустически
|
||||
// есть четвертьволновая труба: вход отражает продольные волны как жёсткий поршень, выход —
|
||||
// как открытый конец. Собственная мода затухает как ν(π/2Nx)², то есть практически не
|
||||
// затухает, и любой стартовый удар остаётся в домене навсегда. Губка — единственное, что
|
||||
// её реально ест, поэтому по умолчанию она включена и подобрана под длину домена.
|
||||
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),
|
||||
@@ -352,9 +389,11 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
pert_amp: cli.pert_amp,
|
||||
pert_dur: cli.pert_dur,
|
||||
outlet_extrapolate: cli.outlet == "extrapolate",
|
||||
sponge_len: cli.sponge_len,
|
||||
init_uniform: cli.init == "uniform",
|
||||
sponge_len,
|
||||
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")?,
|
||||
probe,
|
||||
})
|
||||
}
|
||||
@@ -618,21 +657,52 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) {
|
||||
}
|
||||
None => println!(" L1 {:>12} (измельчение выключено)", "нет"),
|
||||
}
|
||||
println!(" оператор {:>12}", if spec.collision == Collision::Kbc { "KBC-D" } else { "LBGK" });
|
||||
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"
|
||||
});
|
||||
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 {
|
||||
println!(" + губка {} столбцов ×{:.0}", spec.sponge_len, spec.sponge_mult);
|
||||
} else {
|
||||
println!("; губка ВЫКЛЮЧЕНА");
|
||||
}
|
||||
}
|
||||
|
||||
// Пара «вход по скорости / выход по давлению» — недодемпфированный акустический резонатор:
|
||||
// затухание продольной моды ~ ν(π/Nx)², то есть на длинном домене её почти ничто не гасит,
|
||||
// и стартовый транзиент способен раскачать её до смещённого режима или до развала счёта.
|
||||
if spec.sponge_len == 0 && spec.nx >= 250 {
|
||||
println!("
|
||||
⚠ домен длинный ({} столбцов), губка выключена. Продольная акустическая", spec.nx);
|
||||
println!(" мода затухает как ν(π/Nx)² и на такой длине почти не гасится: возможен уход");
|
||||
println!(" ⟨ρ⟩ от единицы (режим смещён, числа несопоставимы) вплоть до развала счёта.");
|
||||
println!(" Смотрите ⟨ρ⟩ в отчёте о сходимости; лечится ключом --sponge-len 32.");
|
||||
}
|
||||
|
||||
if cli.gif.is_some() {
|
||||
println!("\n── Анимация (синхронизация с физическим временем) ───────────────────────────");
|
||||
@@ -791,6 +861,31 @@ fn final_report(
|
||||
println!("rms Cl насыщен ({:+.4}); дрейф-устойчивый Cd = {cdn1:.3}", clr1 - clr0);
|
||||
}
|
||||
|
||||
// ── измеренная продольная пульсация ──
|
||||
// Колебание средней плотности — это и есть та самая «поршневая» мода. Переводим его в
|
||||
// амплитуду скорости: для бегущей акустической волны 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::<R>() / 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 > 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::<Vec<_>>());
|
||||
|
||||
@@ -86,6 +86,33 @@ pub fn feq(rho: R, ux: R, uy: R) -> [R; Q] {
|
||||
]
|
||||
}
|
||||
|
||||
/// Приближение Града — ур. (58) 2D-статьи. Восстанавливает популяции по ρ, u и ПОЛНОМУ
|
||||
/// тензору давлений Π (то есть вместе с неравновесной частью Π⁽¹⁾):
|
||||
///
|
||||
/// f*_i = W_i [ ρ + ρu_α v_iα/c_s² + (1/2c_s⁴)(Π_αβ − ρc_s²δ_αβ)(v_iα v_iβ − c_s²δ_αβ) ]
|
||||
///
|
||||
/// Нужно для корректного СТАРТА: инициализация одним равновесием обнуляет Π⁽¹⁾, схема тратит
|
||||
/// несколько шагов на её восстановление, и этот стартовый дефект портит наблюдаемый порядок
|
||||
/// сходимости. Статья инициализирует эталонные течения именно так.
|
||||
///
|
||||
/// По построению сохраняет ρ и ρ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;
|
||||
let ayy = pyy - rho * CS2;
|
||||
std::array::from_fn(|i| {
|
||||
let cx = CX[i] as R;
|
||||
let cy = CY[i] as R;
|
||||
let quad = axx * (cx * cx - CS2) + 2.0 * pxy * cx * cy + ayy * (cy * cy - CS2);
|
||||
W[i] * (rho + rho * (ux * cx + uy * cy) / CS2 + quad / (2.0 * CS2 * CS2))
|
||||
})
|
||||
}
|
||||
|
||||
/// Макропеременные: ρ = Σ f_i, u = (Σ c_i f_i)/ρ.
|
||||
#[inline]
|
||||
pub fn macros(f: &[R; Q]) -> (R, R, R) {
|
||||
@@ -99,22 +126,58 @@ pub fn macros(f: &[R; Q]) -> (R, R, R) {
|
||||
// KBC: проектор на сдвиг, энтропийный стабилизатор, столкновение
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
|
||||
/// Δs = Ps·Δ — проекция неравновесия на сдвиговую часть модели KBC D: натуральные моменты
|
||||
/// N = M20 − M02 и Π_xy = M11.
|
||||
/// Состав сдвиговой части s — то, что релаксирует жёстко с 2β. Всё остальное попадает в h и
|
||||
/// отдаётся энтропийному стабилизатору γ. Выбор берётся из табл. I 2D-статьи; в 3D-работах
|
||||
/// те же варианты называются N1 и N2.
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
|
||||
pub enum KbcModel {
|
||||
/// s = {N, Π_xy} — только девиатор напряжений. 2D-статья: **KBC D**; 3D: **KBC-N1**.
|
||||
/// След T уходит в h, поэтому объёмная вязкость ξ = c_s²(1/(γβ) − ½) гуляет вместе с γ
|
||||
/// и способна уходить в минус — локальное антизатухание звука (ур. 57).
|
||||
N1,
|
||||
/// s = {N, Π_xy, T} — девиатор плюс след. 2D-статья: **KBC C**; 3D: **KBC-N2**.
|
||||
/// Объёмная вязкость фиксирована и положительна: ξ = ν (ур. 57). Полезно там, где в
|
||||
/// задаче есть слабозатухающая акустика.
|
||||
N2,
|
||||
}
|
||||
|
||||
impl KbcModel {
|
||||
pub fn from_str(s: &str) -> Option<Self> {
|
||||
Some(match s.to_ascii_lowercase().as_str() {
|
||||
"n1" | "d" => KbcModel::N1,
|
||||
"n2" | "c" => KbcModel::N2,
|
||||
_ => return None,
|
||||
})
|
||||
}
|
||||
pub const ALL: [&'static str; 2] = ["n1", "n2"];
|
||||
}
|
||||
|
||||
/// Δs = Ps·Δ — проекция неравновесия на сдвиговую часть выбранной модели.
|
||||
///
|
||||
/// Аналитическая форма, прямо из представления популяций через натуральные моменты (10):
|
||||
/// вклад N сидит только в f(σ,0) (+ρN/4) и f(0,λ) (−ρN/4), вклад Π_xy — только в f(σ,λ)
|
||||
/// (+σλ·ρΠ_xy/4). Поскольку n_i и p_i линейны по своим моментам, а ρ у f и feq одинакова,
|
||||
/// разности берутся покомпонентно. Эквивалентно матричному Ps = M⁻¹·diag(0,0,0,0,1,1,0,0,0)·M
|
||||
/// (проверяется тестом `projector_matches_moment_matrix`), но без матрицы 9×9.
|
||||
/// Аналитическая форма, прямо из представления популяций через натуральные моменты (ур. 10).
|
||||
/// Вклад N = M20 − M02 сидит только в f(σ,0) (+ρN/4) и f(0,λ) (−ρN/4); вклад Π_xy = M11 —
|
||||
/// только в f(σ,λ) (+σλ·ρΠ_xy/4); вклад следа T = M20 + M02 — в f(0,0) (−ρT) и в f(σ,0),
|
||||
/// f(0,λ) (по +ρT/4). Все три вклада линейны по своим моментам, а ρ у f и feq одинакова,
|
||||
/// поэтому разности берутся покомпонентно. Эквивалентно матричному
|
||||
/// Ps = M⁻¹·diag(оставить нужные моменты)·M, но без матрицы 9×9.
|
||||
#[inline]
|
||||
pub fn project_shift(d: &[R; Q]) -> [R; Q] {
|
||||
pub fn project_shift(d: &[R; Q], model: KbcModel) -> [R; Q] {
|
||||
// ρΔN = Σ (c_x² − c_y²)Δ_i и ρΔΠ_xy = Σ c_x c_y Δ_i
|
||||
let dn = d[1] + d[3] - d[2] - d[4];
|
||||
let dp = d[5] + d[7] - d[6] - d[8];
|
||||
let a = 0.25 * dn;
|
||||
let b = 0.25 * dp;
|
||||
[0.0, a, -a, a, -a, b, -b, b, -b]
|
||||
let mut s = [0.0, a, -a, a, -a, b, -b, b, -b];
|
||||
if let KbcModel::N2 = model {
|
||||
// ρΔT = Σ (c_x² + c_y²)Δ_i
|
||||
let dt = d[1] + d[2] + d[3] + d[4] + 2.0 * (d[5] + d[6] + d[7] + d[8]);
|
||||
s[0] -= dt;
|
||||
let q = 0.25 * dt;
|
||||
for v in s.iter_mut().take(5).skip(1) {
|
||||
*v += q;
|
||||
}
|
||||
}
|
||||
s
|
||||
}
|
||||
|
||||
/// Что вернуло столкновение помимо новых популяций — нужно для диагностики KBC.
|
||||
@@ -137,7 +200,7 @@ pub struct Kbc {
|
||||
/// При γ = 2 совпадает с LBGK. Сдвиговые моменты релаксируют с точным 2β = 1/τ при любом γ,
|
||||
/// поэтому кинематическая вязкость (5) от стабилизатора не зависит.
|
||||
#[inline]
|
||||
pub fn collide_node(f: &mut [R; Q], beta: R) -> Kbc {
|
||||
pub fn collide_node(f: &mut [R; Q], beta: R, model: KbcModel) -> Kbc {
|
||||
let (rho, ux, uy) = macros(f);
|
||||
let fe = feq(rho, ux, uy);
|
||||
|
||||
@@ -145,7 +208,7 @@ pub fn collide_node(f: &mut [R; Q], beta: R) -> Kbc {
|
||||
for i in 0..Q {
|
||||
d[i] = f[i] - fe[i];
|
||||
}
|
||||
let ds = project_shift(&d);
|
||||
let ds = project_shift(&d, model);
|
||||
|
||||
let (mut num, mut den, mut nrm) = (0.0, 0.0, 0.0);
|
||||
for i in 0..Q {
|
||||
@@ -197,8 +260,13 @@ pub fn beta_of_nu(nu: R) -> R {
|
||||
/// Отрицательная ξ означает локальное антизатухание акустической моды — поэтому в отчёте
|
||||
/// печатается доля узлов с ξ < 0 (см. report в main.rs).
|
||||
#[inline]
|
||||
pub fn xi_of_gamma(gamma: R, beta: R) -> R {
|
||||
CS2 * (1.0 / (gamma * beta) - 0.5)
|
||||
pub fn xi_of_gamma(gamma: R, beta: R, model: KbcModel) -> R {
|
||||
match model {
|
||||
// след в h — ξ зависит от стабилизатора и может стать отрицательной
|
||||
KbcModel::N1 => CS2 * (1.0 / (gamma * beta) - 0.5),
|
||||
// след в s — ξ равна сдвиговой вязкости и всегда положительна
|
||||
KbcModel::N2 => nu_of_beta(beta),
|
||||
}
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
@@ -737,17 +805,18 @@ mod tests {
|
||||
use super::*;
|
||||
|
||||
/// Проектор Δs, построенный аналитически, обязан совпадать с матричным
|
||||
/// Ps = M⁻¹·diag(…,1,1,…)·M, где M — базис натуральных моментов из (6)–(7).
|
||||
/// Ps = M⁻¹·D·M, где M — базис натуральных моментов (6)–(7), а D оставляет ровно те
|
||||
/// моменты, которые входят в сдвиговую часть выбранной модели (табл. I).
|
||||
#[test]
|
||||
fn projector_matches_moment_matrix() {
|
||||
// строки M: 1, cx, cy, 3(cx²+cy²)−2, cx²−cy², cx·cy, cx²cy, cx cy², cx²cy²
|
||||
// строки M: 1, cx, cy, T=cx²+cy², N=cx²−cy², Π_xy=cx·cy, cx²cy, cx cy², cx²cy²
|
||||
let mut m = [[0.0f64; Q]; Q];
|
||||
for i in 0..Q {
|
||||
let (x, y) = (CX[i] as f64, CY[i] as f64);
|
||||
m[0][i] = 1.0;
|
||||
m[1][i] = x;
|
||||
m[2][i] = y;
|
||||
m[3][i] = 3.0 * (x * x + y * y) - 2.0;
|
||||
m[3][i] = x * x + y * y;
|
||||
m[4][i] = x * x - y * y;
|
||||
m[5][i] = x * y;
|
||||
m[6][i] = x * x * y;
|
||||
@@ -755,29 +824,62 @@ mod tests {
|
||||
m[8][i] = x * x * y * y;
|
||||
}
|
||||
let minv = invert9(&m);
|
||||
// Ps = M⁻¹ · D · M, D оставляет строки 4 (N) и 5 (Π_xy)
|
||||
let mut ps = [[0.0f64; Q]; Q];
|
||||
for a in 0..Q {
|
||||
for b in 0..Q {
|
||||
ps[a][b] = minv[a][4] * m[4][b] + minv[a][5] * m[5][b];
|
||||
}
|
||||
}
|
||||
// сверяем на случайных (детерминированных) векторах
|
||||
let mut seed = 12345u64;
|
||||
let mut rnd = || {
|
||||
seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
|
||||
((seed >> 33) as f64 / (1u64 << 31) as f64) - 1.0
|
||||
};
|
||||
for _ in 0..200 {
|
||||
let d: [f64; Q] = std::array::from_fn(|_| rnd());
|
||||
let got = project_shift(&d);
|
||||
for (model, keep) in
|
||||
[(KbcModel::N1, vec![4usize, 5]), (KbcModel::N2, vec![3usize, 4, 5])]
|
||||
{
|
||||
let mut ps = [[0.0f64; Q]; Q];
|
||||
for a in 0..Q {
|
||||
let want: f64 = (0..Q).map(|b| ps[a][b] * d[b]).sum();
|
||||
assert!((got[a] - want).abs() < 1e-13, "i={a}: {} vs {}", got[a], want);
|
||||
for b in 0..Q {
|
||||
ps[a][b] = keep.iter().map(|&k| minv[a][k] * m[k][b]).sum();
|
||||
}
|
||||
}
|
||||
let mut seed = 12345u64;
|
||||
let mut rnd = || {
|
||||
seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
|
||||
((seed >> 33) as f64 / (1u64 << 31) as f64) - 1.0
|
||||
};
|
||||
for _ in 0..200 {
|
||||
// проектор сравниваем на ФИЗИЧЕСКИ допустимых Δ: у них моменты массы и
|
||||
// импульса нулевые (и f, и feq несут одни и те же ρ, ρu)
|
||||
let mut d: [f64; Q] = std::array::from_fn(|_| rnd());
|
||||
let mass: f64 = d.iter().sum();
|
||||
let mx: f64 = (0..Q).map(|i| CX[i] as f64 * d[i]).sum();
|
||||
let my: f64 = (0..Q).map(|i| CY[i] as f64 * d[i]).sum();
|
||||
for i in 0..Q {
|
||||
d[i] -= W[i] * (mass + 3.0 * (CX[i] as f64 * mx + CY[i] as f64 * my));
|
||||
}
|
||||
let got = project_shift(&d, model);
|
||||
for a in 0..Q {
|
||||
let want: f64 = (0..Q).map(|b| ps[a][b] * d[b]).sum();
|
||||
assert!(
|
||||
(got[a] - want).abs() < 1e-13,
|
||||
"{model:?} i={a}: {} vs {want}",
|
||||
got[a]
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Объёмная вязкость по ур. (57): у модели со следом в сдвиговой части она равна
|
||||
/// сдвиговой и не зависит от стабилизатора, у модели без следа — зависит.
|
||||
#[test]
|
||||
fn bulk_viscosity_follows_model() {
|
||||
let beta = 1.0 / (2.0 * 0.7);
|
||||
let nu = nu_of_beta(beta);
|
||||
for &g in &[0.5, 1.0, 1.7, 3.0] {
|
||||
assert!(
|
||||
(xi_of_gamma(g, beta, KbcModel::N2) - nu).abs() < 1e-15,
|
||||
"N2: ξ обязана равняться ν при любой γ"
|
||||
);
|
||||
let want = CS2 * (1.0 / (g * beta) - 0.5);
|
||||
assert!((xi_of_gamma(g, beta, KbcModel::N1) - want).abs() < 1e-15);
|
||||
}
|
||||
// у N1 ξ уходит в минус при достаточно большой γ — то самое антизатухание звука
|
||||
let g_neg = 2.0 / (beta * CS2) * CS2;
|
||||
assert!(xi_of_gamma(g_neg * 2.0, beta, KbcModel::N1) < 0.0);
|
||||
}
|
||||
|
||||
/// Ps обязан быть идемпотентным (это проектор) и не трогать сохраняющиеся моменты.
|
||||
#[test]
|
||||
fn projector_is_idempotent_and_conservative() {
|
||||
@@ -788,8 +890,9 @@ mod tests {
|
||||
};
|
||||
for _ in 0..100 {
|
||||
let d: [f64; Q] = std::array::from_fn(|_| rnd());
|
||||
let s1 = project_shift(&d);
|
||||
let s2 = project_shift(&s1);
|
||||
for model in [KbcModel::N1, KbcModel::N2] {
|
||||
let s1 = project_shift(&d, model);
|
||||
let s2 = project_shift(&s1, model);
|
||||
for i in 0..Q {
|
||||
assert!((s1[i] - s2[i]).abs() < 1e-14);
|
||||
}
|
||||
@@ -798,6 +901,7 @@ mod tests {
|
||||
let mx: f64 = (0..Q).map(|i| CX[i] as f64 * s1[i]).sum();
|
||||
let my: f64 = (0..Q).map(|i| CY[i] as f64 * s1[i]).sum();
|
||||
assert!(mx.abs() < 1e-13 && my.abs() < 1e-13);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -823,7 +927,7 @@ mod tests {
|
||||
let (rho, ux, uy) = macros(&f);
|
||||
let fe = feq(rho, ux, uy);
|
||||
let d: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]);
|
||||
let ds = project_shift(&d);
|
||||
let ds = project_shift(&d, KbcModel::N1);
|
||||
let dh: [f64; Q] = std::array::from_fn(|i| d[i] - ds[i]);
|
||||
|
||||
// невязка условия (15)
|
||||
@@ -835,7 +939,7 @@ mod tests {
|
||||
.sum()
|
||||
};
|
||||
let mut fc = f;
|
||||
let info = collide_node(&mut fc, beta);
|
||||
let info = collide_node(&mut fc, beta, KbcModel::N1);
|
||||
let r0 = resid(info.gamma).abs();
|
||||
let r1 = resid(info.gamma + 1.0).abs().min(resid(info.gamma - 1.0).abs());
|
||||
assert!(r0 < 0.05 * r1, "невязка (15) при γ*: {r0:.3e}, на ±1 от него: {r1:.3e}");
|
||||
@@ -850,10 +954,10 @@ mod tests {
|
||||
let (rho, ux, uy) = macros(&f);
|
||||
let fe = feq(rho, ux, uy);
|
||||
let d: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]);
|
||||
let ds = project_shift(&d);
|
||||
let ds = project_shift(&d, KbcModel::N1);
|
||||
let mut pure_shear: [f64; Q] = std::array::from_fn(|i| fe[i] + ds[i]);
|
||||
let mut bgk = pure_shear;
|
||||
collide_node(&mut pure_shear, beta);
|
||||
collide_node(&mut pure_shear, beta, KbcModel::N1);
|
||||
collide_node_bgk(&mut bgk, beta);
|
||||
for i in 0..Q {
|
||||
assert!((pure_shear[i] - bgk[i]).abs() < 1e-12, "i={i}");
|
||||
@@ -864,20 +968,26 @@ mod tests {
|
||||
/// что вязкость (5) задаётся только β, а стабилизатор её не трогает.
|
||||
#[test]
|
||||
fn shear_relaxes_at_two_beta_regardless_of_gamma() {
|
||||
for &beta in &[0.3, 0.6, 0.95] {
|
||||
let mut f = perturbed_state(0.8);
|
||||
let (rho, ux, uy) = macros(&f);
|
||||
let fe = feq(rho, ux, uy);
|
||||
let d0: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]);
|
||||
let s0 = project_shift(&d0);
|
||||
collide_node(&mut f, beta);
|
||||
let (r1, x1, y1) = macros(&f);
|
||||
let fe1 = feq(r1, x1, y1);
|
||||
let d1: [f64; Q] = std::array::from_fn(|i| f[i] - fe1[i]);
|
||||
let s1 = project_shift(&d1);
|
||||
for i in 0..Q {
|
||||
let want = (1.0 - 2.0 * beta) * s0[i];
|
||||
assert!((s1[i] - want).abs() < 1e-10, "β={beta} i={i}: {} vs {want}", s1[i]);
|
||||
for model in [KbcModel::N1, KbcModel::N2] {
|
||||
for &beta in &[0.3, 0.6, 0.95] {
|
||||
let mut f = perturbed_state(0.8);
|
||||
let (rho, ux, uy) = macros(&f);
|
||||
let fe = feq(rho, ux, uy);
|
||||
let d0: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]);
|
||||
let s0 = project_shift(&d0, model);
|
||||
collide_node(&mut f, beta, model);
|
||||
let (r1, x1, y1) = macros(&f);
|
||||
let fe1 = feq(r1, x1, y1);
|
||||
let d1: [f64; Q] = std::array::from_fn(|i| f[i] - fe1[i]);
|
||||
let s1 = project_shift(&d1, model);
|
||||
for i in 0..Q {
|
||||
let want = (1.0 - 2.0 * beta) * s0[i];
|
||||
assert!(
|
||||
(s1[i] - want).abs() < 1e-10,
|
||||
"{model:?} β={beta} i={i}: {} vs {want}",
|
||||
s1[i]
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -888,7 +998,7 @@ mod tests {
|
||||
for &beta in &[0.25, 0.5, 0.99] {
|
||||
let mut f = perturbed_state(0.7);
|
||||
let before = macros(&f);
|
||||
collide_node(&mut f, beta);
|
||||
collide_node(&mut f, beta, KbcModel::N1);
|
||||
let after = macros(&f);
|
||||
assert!((before.0 - after.0).abs() < 1e-13);
|
||||
assert!((before.1 - after.1).abs() < 1e-12);
|
||||
|
||||
Reference in New Issue
Block a user