解决NumPy广播错误：在LBM CFD求解器中处理三维数组形状不匹配问题

def eq(geq,rho,ux,uy):
    # Calcul de la fonction d'équilibre
    geq[:, :, 0] = w[0] * rho * (1 - 0.5 * (c0**(-2)) * (ux**2 + uy**2))
    # 错误发生在此行
    geq[:, :, 1:9] = w[1:] * rho * (1 + (c0**(-2)) * (ca[1:9, 0]*ux + ca[1:9, 1]*uy) + 0.5* (c0**-4) * (ca[1:9, 0]*ux + ca[1:9, 1]*uy)**2 - 0.5 * (c0**(-2)) * (ux**2 + uy**2))

import numpy as np
import matplotlib.pyplot as plt
import time

# Parametres du modele D2Q9 (部分参数定义，用于代码完整性)
Re= 20
ca= np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1], [1, 1], [-1, 1], [-1, -1], [1, -1]])
w= np.array([4/9, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36])
c0= 1/np.sqrt(3)

D= 4
nx,ny= 20*D, 10*D
M0= 0.3

taug= (M0*D)/(c0*Re) + 1
nt= int((150*D)/(M0*c0))

mask=np.zeros((nx,ny));
cx1,cx2= int(8*D - 0.5*D), int(8*D + 0.5*D)
cy1,cy2= int(5*D - 0.5*D), int(5*D + 0.5*D)
mask[cx1:cx2,cy1:cy2]=1


def eq(geq, rho, ux, uy):
    """
    计算平衡态分布函数。
    通过扩展数组维度以实现NumPy广播兼容性。
    """
    # 扩展宏观变量的维度，使其形状从 (nx, ny) 变为 (nx, ny, 1)
    uxb = ux[:, :, None]
    uyb = uy[:, :, None]
    rhob = rho[:, :, None]

    # 扩展权重系数w的维度，使其形状从 (9,) 变为 (1, 1, 9)
    # 这样 w[..., 1:] 的形状变为 (1, 1, 8)
    wb = w[None, None, :]

    # 扩展离散速度ca的维度，使其形状从 (9, 2) 变为 (1, 1, 9, 2)
    # 这样 ca[..., 1:9, 0] 和 ca[..., 1:9, 1] 的形状都变为 (1, 1, 8)
    cab = ca[None, None, :, :] # 扩展所有维度，然后切片

    # 计算 geq[:, :, 0]
    geq[:, :, 0] = w[0] * rho * (1 - 0.5 * (c0**(-2)) * (ux**2 + uy**2))

    # 计算 geq[:, :, 1:9]
    # 使用扩展后的数组进行广播计算
    # 最终结果的形状将是 (nx, ny, 8)
    geq[:, :, 1:9] = wb[..., 1:] * (rhob * (1 + (c0**(-2)) * (cab[..., 1:9, 0]*uxb + cab[..., 1:9, 1]*uyb) + \
                                            0.5 * (c0**-4) * (cab[..., 1:9, 0]*uxb + cab[..., 1:9, 1]*uyb)**2 - \
                                            0.5 * (c0**(-2)) * (uxb**2 + uyb**2)))

# 其他函数（init, collide, propagate, macro, boundary, wall）保持不变
def init(M0):
    rho= np.ones((nx, ny)) * c0**2
    ux= np.full((nx, ny), M0 * c0)
    uy= np.zeros((nx, ny))
    geq=np.zeros((nx,ny,9))
    eq(geq,rho,ux,uy) # 确保在init中调用eq
    return geq,rho,ux,uy

def collide(gcoll,g,geq,taug):
    gcoll[:,:,:]= g[:,:,:] - (1/taug)*(g[:,:,:] - geq[:,:,:])

def propagate(g,gcoll):
    g[:, :, 0] = gcoll[:, :, 0]
    g[:, :, 1] = gcoll[:, :, 1]
    g[:, :, 2] = gcoll[:, :, 2]
    g[:, :, 3] = gcoll[:, :, 3]
    g[:, :, 4] = gcoll[:, :, 4]
    g[:, :, 5] = gcoll[:, :, 5]
    g[:, :, 6] = gcoll[:, :, 6]
    g[:, :, 7] = gcoll[:, :, 7]
    g[:, :, 8] = gcoll[:, :, 8]

def macro(g,rho,ux,uy):
    rho[:, :] = np