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>
35 lines
1.6 KiB
Python
35 lines
1.6 KiB
Python
# Решётка D2Q9 и KBC-проектор. Чистый CuPy; всё считается один раз при импорте.
|
|
#
|
|
# Скорости e_i и веса w_i (стандартная нумерация D2Q9):
|
|
# 0:( 0, 0) 1:(+1,0) 2:(0,+1) 3:(-1,0) 4:(0,-1)
|
|
# 5:(+1,+1) 6:(-1,+1) 7:(-1,-1) 8:(+1,-1)
|
|
# OPP[i] — индекс противоположного направления.
|
|
import backend as B
|
|
|
|
cp = B.cp; xp = B.xp; DTYPE = B.DTYPE
|
|
|
|
Q = 9
|
|
Cx = [0, 1, 0, -1, 0, 1, -1, -1, 1] # python-инты (для roll/индексации)
|
|
Cy = [0, 0, 1, 0, -1, 1, 1, -1, -1]
|
|
OPP = [0, 3, 4, 1, 2, 7, 8, 5, 6]
|
|
_Wlist = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
|
|
CS2 = 1.0 / 3.0 # квадрат скорости звука решётки
|
|
|
|
# device-векторы
|
|
Wa = cp.asarray(_Wlist, DTYPE)
|
|
CXa = cp.asarray(Cx, DTYPE)
|
|
CYa = cp.asarray(Cy, DTYPE)
|
|
|
|
def _build_projector():
|
|
"""Проектор Ps на сдвиг-моменты {cx²−cy², cx·cy} (ядро KBC-N1).
|
|
Ps = M⁻¹ · diag(оставить только эти 2 момента) · M, где M — моментный базис."""
|
|
cxv = cp.asarray(Cx, cp.float64); cyv = cp.asarray(Cy, cp.float64)
|
|
M = cp.zeros((Q, Q), cp.float64)
|
|
M[0] = 1.0; M[1] = cxv; M[2] = cyv; M[3] = 3.0*(cxv**2 + cyv**2) - 2.0
|
|
M[4] = cxv**2 - cyv**2; M[5] = cxv*cyv
|
|
M[6] = cxv**2*cyv; M[7] = cxv*cyv**2; M[8] = cxv**2*cyv**2
|
|
Dm = cp.zeros((Q, Q), cp.float64); Dm[4, 4] = Dm[5, 5] = 1.0
|
|
return (cp.linalg.inv(M) @ Dm @ M)
|
|
|
|
Ps = _build_projector().astype(DTYPE) # (Q,Q) проектор на сдвиг
|