NotBigGhostandClaude Opus 5
11ff7b79b4
Состояние на момент заведения репозитория. C++ приложение (src/, shaders/, tests/) — минимальный редактор 3D-моделей на Vulkan 1.3: орбитальная камера, три опорные сетки через начало координат, загрузка .obj с режимами отображения. Весь Vulkan изолирован в src/vk/. Исследование (docs/) — оригинальные статьи по KBC (docs/origins) и Python-решатель D2Q9 KBC-N1 с AMR 2x и SDF+Bouzidi (docs/theory). В решателе перед коммитом исправлены дефекты, найденные сверкой с первоисточниками: относительный порог знаменателя энтропийного стабилизатора (абсолютный вырождал KBC в LBGK на 77-99% узлов), заворот вход/выход в углах домена, диагностика средней плотности по фиктивным узлам тела, зашитый refine=2. Подробности — docs/theory/solver_2x_sdf/README.md. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
53 lines
3.8 KiB
Python
53 lines
3.8 KiB
Python
# СИЛА и МОМЕНТ на теле — галилей-инвариантный обмен импульсом (GMEM, Wen et al. 2014).
|
||
#
|
||
# Для каждого линка жидкость→тело по направлению i (маска cl):
|
||
# ΔF = (c_i − u_w)·f_i^{после столкновения}(x_f) − (c_ī − u_w)·f_ī^{после стриминга}(x_f)
|
||
# где u_w — скорость стенки на линке (для неподвижного цилиндра = 0). Т.к. c_ī = −c_i, по компонентам:
|
||
# F_x += (Cx_i − u_wx)·f_i + (Cx_i + u_wx)·f_ī
|
||
# F_y += (Cy_i − u_wy)·f_i + (Cy_i + u_wy)·f_ī
|
||
#
|
||
# МОМЕНТ (торк) вокруг центра тела: T_z = Σ_links (r_w × ΔF)_z, где точка приложения —
|
||
# ПЕРЕСЕЧЕНИЕ линка со стенкой r_w = r_f + q·c_i (q — Bouzidi-доля, уже лежит в bc).
|
||
# Узел жидкости дал бы ошибку плеча до 1 ячейки; q бесплатен. Для кругового цилиндра
|
||
# ⟨Cm⟩ ≈ 0 — это контроль симметрии считывания; для подвижных/вращающихся тел (цель SimV4) —
|
||
# рабочая величина.
|
||
#
|
||
# Свойства:
|
||
# * f_i — популяция, уходящая в стенку (после столкновения); f_ī — вернувшаяся (после стриминга/Bouzidi).
|
||
# f_ī берётся из поля ПОСЛЕ стриминга — иначе ведущий симметричный член ~2ρw_i сокращается и Cd врёт.
|
||
# * При u_w=0 сводится к классическому Σ c_i (f_i + f_ī); члены с u_w делают оценку галилей-инвариантной
|
||
# и ГОТОВОЙ к подвижным/вращающимся телам (цель проекта SimV4) — туда передаётся u_w на линке.
|
||
#
|
||
# Независимый кросс-чек этой оценки — интегрирование тензора напряжений по контуру вокруг тела:
|
||
# см. forces_stress.py (σ = −p′·I + σ^v из неравновесных моментов; вызывается редко, вне CUDA-графа).
|
||
import backend as B
|
||
import lattice as L
|
||
|
||
xp = B.xp; ZERO = B.ZERO
|
||
Cx = L.Cx; Cy = L.Cy
|
||
|
||
def force_gmem(fpost, fnew, bc, rxy=None, u_wall=(0.0, 0.0)):
|
||
"""GMEM по линкам тела (маска cl). Возвращает (Fx, Fy, Tz) — device-скаляры.
|
||
rxy=(rx,ry) — координаты узлов решётки относительно центра тела (для торка); None → Tz=0.
|
||
u_wall — скорость стенки. Graph-safe: только where/арифметика с литералами."""
|
||
uwx, uwy = u_wall
|
||
Fx = ZERO; Fy = ZERO; Tz = ZERO
|
||
if rxy is not None:
|
||
rx, ry = rxy
|
||
for (i, ib, mask, q, xff, cl) in bc:
|
||
fi = fpost[i]; fib = fnew[ib]
|
||
dFx = xp.where(cl, (Cx[i] - uwx)*fi + (Cx[i] + uwx)*fib, ZERO)
|
||
dFy = xp.where(cl, (Cy[i] - uwy)*fi + (Cy[i] + uwy)*fib, ZERO)
|
||
Fx = Fx + dFx.sum(); Fy = Fy + dFy.sum()
|
||
if rxy is not None: # плечо до точки пересечения линка со стенкой
|
||
Tz = Tz + ((rx + q*Cx[i])*dFy - (ry + q*Cy[i])*dFx).sum()
|
||
return Fx, Fy, Tz
|
||
|
||
def coefficients(Fx, Fy, U, DL, Tz=None):
|
||
"""Безразмерные коэффициенты: Cd = 2Fx/(U²·DL), Cl = 2Fy/(U²·DL), Cm = 2Tz/(U²·DL²).
|
||
DL — диаметр у тела (на L1); сила/плечо тоже в единицах L1, так что нормировка согласована."""
|
||
Cd = 2.0*Fx/(U**2*DL); Cl = 2.0*Fy/(U**2*DL)
|
||
if Tz is None:
|
||
return Cd, Cl
|
||
return Cd, Cl, 2.0*Tz/(U**2*DL*DL)
|