# Решётка 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) проектор на сдвиг