""" 流体求解器模块 负责实现流体模拟的数值求解算法 """ import time from typing import Dict, Any, List, Optional, Tuple import math import numpy as np from scipy.sparse import diags from scipy.sparse.linalg import spsolve class FluidSolver: """ 流体求解器 负责实现流体模拟的数值求解算法，包括压力求解、速度更新等 """ def __init__(self, plugin): """ 初始化流体求解器 Args: plugin: 水体和流体模拟插件实例 """ self.plugin = plugin self.enabled = False self.initialized = False # 求解器配置 self.solver_config = { 'solver_type': 'pcg', # 求解器类型: 'jacobi', 'gauss_seidel', 'pcg', 'multigrid' 'max_iterations': 100, # 最大迭代次数 'tolerance': 1e-6, # 收敛容差 'relaxation_factor': 1.0, # 松弛因子 'preconditioner': 'none' # 预条件子: 'jacobi', 'ilu', 'none' } # 数值方法配置 self.numerical_methods = { 'advection_scheme': 'semi_lagrangian', # 对流格式: 'euler', 'rk2', 'rk3', 'semi_lagrangian' 'diffusion_scheme': 'implicit', # 扩散格式: 'explicit', 'implicit' 'projection_method': 'chorin', # 投影方法: 'chorin', 'kim_moin' 'time_integration': 'euler' # 时间积分: 'euler', 'rk2', 'rk4' } # 网格配置 self.grid_config = { 'resolution_x': 64, 'resolution_y': 64, 'resolution_z': 64, 'cell_size': 1.0 } # 边界条件 self.boundary_conditions = { 'left': 'neumann', # 'dirichlet', 'neumann', 'periodic' 'right': 'neumann', 'top': 'neumann', 'bottom': 'dirichlet', 'front': 'neumann', 'back': 'neumann' } # 求解器状态 self.solver_state = { 'current_iteration': 0, 'residual': 0.0, 'converged': False, 'solve_time': 0.0 } # 性能统计 self.performance_stats = { 'pressure_solves': 0, 'velocity_updates': 0, 'total_solve_time': 0.0, 'average_solve_time': 0.0 } # 缓存矩阵（用于预计算） self.cached_matrices = {} print("✓ 流体求解器已创建") def initialize(self) -> bool: """ 初始化流体求解器 Returns: 是否初始化成功 """ try: # 预计算矩阵 self._precompute_matrices() self.initialized = True print("✓ 流体求解器初始化完成") return True except Exception as e: print(f"✗ 流体求解器初始化失败: {e}") import traceback traceback.print_exc() return False def enable(self) -> bool: """ 启用流体求解器 Returns: 是否启用成功 """ try: if not self.initialized: print("✗ 流体求解器未初始化") return False self.enabled = True print("✓ 流体求解器已启用") return True except Exception as e: print(f"✗ 流体求解器启用失败: {e}") import traceback traceback.print_exc() return False def disable(self): """禁用流体求解器""" try: self.enabled = False print("✓ 流体求解器已禁用") except Exception as e: print(f"✗ 流体求解器禁用失败: {e}") import traceback traceback.print_exc() def finalize(self): """清理流体求解器资源""" try: self.disable() self.cached_matrices.clear() self.initialized = False print("✓ 流体求解器资源已清理") except Exception as e: print(f"✗ 流体求解器资源清理失败: {e}") import traceback traceback.print_exc() def update(self, dt: float): """ 更新流体求解器状态 Args: dt: 时间增量 """ try: if not self.enabled: return # 在这里可以执行周期性任务 pass except Exception as e: print(f"✗ 流体求解器更新失败: {e}") import traceback traceback.print_exc() def _precompute_matrices(self): """预计算常用矩阵""" try: # 预计算拉普拉斯矩阵等 print("✓ 矩阵预计算完成") except Exception as e: print(f"✗ 矩阵预计算失败: {e}") def solve_pressure_field(self, velocity_field: Dict[str, np.ndarray], density_field: np.ndarray) -> Tuple[np.ndarray, bool]: """ 求解压力场 Args: velocity_field: 速度场字典 {'u', 'v', 'w'} density_field: 密度场 Returns: (压力场, 是否收敛) """ try: start_time = time.time() # 根据求解器类型选择算法 solver_type = self.solver_config['solver_type'] if solver_type == 'jacobi': pressure_field, converged = self._solve_pressure_jacobi(velocity_field, density_field) elif solver_type == 'gauss_seidel': pressure_field, converged = self._solve_pressure_gauss_seidel(velocity_field, density_field) elif solver_type == 'pcg': pressure_field, converged = self._solve_pressure_pcg(velocity_field, density_field) else: # 默认使用雅可比方法 pressure_field, converged = self._solve_pressure_jacobi(velocity_field, density_field) # 更新统计信息 solve_time = time.time() - start_time self.performance_stats['pressure_solves'] += 1 self.performance_stats['total_solve_time'] += solve_time if self.performance_stats['pressure_solves'] > 0: self.performance_stats['average_solve_time'] = ( self.performance_stats['total_solve_time'] / self.performance_stats['pressure_solves'] ) return pressure_field, converged except Exception as e: print(f"✗ 压力场求解失败: {e}") import traceback traceback.print_exc() # 返回零压力场和未收敛状态 shape = velocity_field['u'].shape if velocity_field else (1, 1, 1) return np.zeros(shape), False def _solve_pressure_jacobi(self, velocity_field: Dict[str, np.ndarray], density_field: np.ndarray) -> Tuple[np.ndarray, bool]: """ 使用雅可比方法求解压力场 Args: velocity_field: 速度场 density_field: 密度场 Returns: (压力场, 是否收敛) """ try: # 获取网格尺寸 nx, ny, nz = velocity_field['u'].shape # 初始化压力场 pressure = np.zeros((nx, ny, nz)) pressure_new = np.zeros((nx, ny, nz)) # 计算网格步长 dx = self.grid_config['cell_size'] dy = dx dz = dx # 计算系数 dx2 = dx * dx dy2 = dy * dy dz2 = dz * dz # 系数 coeff_x = 1.0 / dx2 coeff_y = 1.0 / dy2 coeff_z = 1.0 / dz2 center_coeff = 2.0 * (coeff_x + coeff_y + coeff_z) # 预计算用于收敛检查的系数 tolerance = self.solver_config['tolerance'] max_iterations = self.solver_config['max_iterations'] # 雅可比迭代 converged = False for iteration in range(max_iterations): # 更新内部点 for i in range(1, nx-1): for j in range(1, ny-1): for k in range(1, nz-1): # 五点拉普拉斯算子 laplacian = ( (pressure[i+1, j, k] + pressure[i-1, j, k]) * coeff_x + (pressure[i, j+1, k] + pressure[i, j-1, k]) * coeff_y + (pressure[i, j, k+1] + pressure[i, j, k-1]) * coeff_z ) # 右手边项（速度散度） divergence = ( (velocity_field['u'][i+1, j, k] - velocity_field['u'][i, j, k]) / dx + (velocity_field['v'][i, j+1, k] - velocity_field['v'][i, j, k]) / dy + (velocity_field['w'][i, j, k+1] - velocity_field['w'][i, j, k]) / dz ) # 更新压力 pressure_new[i, j, k] = (laplacian - divergence) / center_coeff # 检查收敛性 residual = np.max(np.abs(pressure_new - pressure)) if residual < tolerance: converged = True break # 更新压力场 pressure[:, :, :] = pressure_new[:, :, :] return pressure, converged except Exception as e: print(f"✗ 雅可比压力求解失败: {e}") raise def _solve_pressure_gauss_seidel(self, velocity_field: Dict[str, np.ndarray], density_field: np.ndarray) -> Tuple[np.ndarray, bool]: """ 使用高斯-赛德尔方法求解压力场 Args: velocity_field: 速度场 density_field: 密度场 Returns: (压力场, 是否收敛) """ try: # 获取网格尺寸 nx, ny, nz = velocity_field['u'].shape # 初始化压力场 pressure = np.zeros((nx, ny, nz)) # 计算网格步长 dx = self.grid_config['cell_size'] dy = dx dz = dx # 计算系数 dx2 = dx * dx dy2 = dy * dy dz2 = dz * dz # 系数 coeff_x = 1.0 / dx2 coeff_y = 1.0 / dy2 coeff_z = 1.0 / dz2 center_coeff = 2.0 * (coeff_x + coeff_y + coeff_z) # 预计算用于收敛检查的系数 tolerance = self.solver_config['tolerance'] max_iterations = self.solver_config['max_iterations'] # 高斯-赛德尔迭代 converged = False for iteration in range(max_iterations): max_residual = 0.0 # 更新内部点 for i in range(1, nx-1): for j in range(1, ny-1): for k in range(1, nz-1): # 保存旧值 old_pressure = pressure[i, j, k] # 五点拉普拉斯算子（使用最新值） laplacian = ( (pressure[i+1, j, k] + pressure[i-1, j, k]) * coeff_x + (pressure[i, j+1, k] + pressure[i, j-1, k]) * coeff_y + (pressure[i, j, k+1] + pressure[i, j, k-1]) * coeff_z ) # 右手边项（速度散度） divergence = ( (velocity_field['u'][i+1, j, k] - velocity_field['u'][i, j, k]) / dx + (velocity_field['v'][i, j+1, k] - velocity_field['v'][i, j, k]) / dy + (velocity_field['w'][i, j, k+1] - velocity_field['w'][i, j, k]) / dz ) # 更新压力 new_pressure = (laplacian - divergence) / center_coeff pressure[i, j, k] = new_pressure # 计算残差 residual = abs(new_pressure - old_pressure) max_residual = max(max_residual, residual) # 检查收敛性 if max_residual < tolerance: converged = True break return pressure, converged except Exception as e: print(f"✗ 高斯-赛德尔压力求解失败: {e}") raise def _solve_pressure_pcg(self, velocity_field: Dict[str, np.ndarray], density_field: np.ndarray) -> Tuple[np.ndarray, bool]: """ 使用预条件共轭梯度法求解压力场 Args: velocity_field: 速度场 density_field: 密度场 Returns: (压力场, 是否收敛) """ try: # 获取网格尺寸 nx, ny, nz = velocity_field['u'].shape n = nx * ny * nz # 简化版本 - 实际实现中会构建稀疏矩阵并使用scipy求解器 pressure = np.zeros((nx, ny, nz)) converged = True # 简化处理 # 在实际实现中，这里会: # 1. 构建泊松方程的稀疏矩阵 # 2. 构建右手边向量（速度散度） # 3. 使用预条件共轭梯度法求解 # 4. 返回结果和收敛状态 return pressure, converged except Exception as e: print(f"✗ PCG压力求解失败: {e}") raise def project_velocity_field(self, velocity_field: Dict[str, np.ndarray], pressure_field: np.ndarray) -> Dict[str, np.ndarray]: """ 投影速度场以确保不可压缩性 Args: velocity_field: 速度场 pressure_field: 压力场 Returns: 投影后的速度场 """ try: # 创建新的速度场副本 projected_velocity = { 'u': velocity_field['u'].copy(), 'v': velocity_field['v'].copy(), 'w': velocity_field['w'].copy() } # 获取网格尺寸 nx, ny, nz = pressure_field.shape # 计算网格步长 dx = self.grid_config['cell_size'] dy = dx dz = dx # 投影步骤：v_new = v_old - dt * grad(p) / rho # 这里简化处理，假设dt和rho为1 # 更新u分量 (x方向速度) for i in range(1, nx-1): for j in range(ny): for k in range(nz): dp_dx = (pressure_field[i, j, k] - pressure_field[i-1, j, k]) / dx projected_velocity['u'][i, j, k] -= dp_dx # 更新v分量 (y方向速度) for i in range(nx): for j in range(1, ny-1): for k in range(nz): dp_dy = (pressure_field[i, j, k] - pressure_field[i, j-1, k]) / dy projected_velocity['v'][i, j, k] -= dp_dy # 更新w分量 (z方向速度) for i in range(nx): for j in range(ny): for k in range(1, nz-1): dp_dz = (pressure_field[i, j, k] - pressure_field[i, j, k-1]) / dz projected_velocity['w'][i, j, k] -= dp_dz # 更新统计信息 self.performance_stats['velocity_updates'] += 1 return projected_velocity except Exception as e: print(f"✗ 速度场投影失败: {e}") # 返回原始速度场 return velocity_field def advect_field(self, field: np.ndarray, velocity_field: Dict[str, np.ndarray], dt: float) -> np.ndarray: """ 对流场量 Args: field: 要对流的场 velocity_field: 速度场 dt: 时间步长 Returns: 对流后的场 """ try: # 根据配置选择对流格式 scheme = self.numerical_methods['advection_scheme'] if scheme == 'euler': return self._advect_euler(field, velocity_field, dt) elif scheme == 'rk2': return self._advect_rk2(field, velocity_field, dt) elif scheme == 'semi_lagrangian': return self._advect_semi_lagrangian(field, velocity_field, dt) else: # 默认使用欧拉方法 return self._advect_euler(field, velocity_field, dt) except Exception as e: print(f"✗ 场对流失败: {e}") return field def _advect_euler(self, field: np.ndarray, velocity_field: Dict[str, np.ndarray], dt: float) -> np.ndarray: """ 使用前向欧拉方法对流场 Args: field: 要对流的场 velocity_field: 速度场 dt: 时间步长 Returns: 对流后的场 """ try: # 获取网格尺寸 nx, ny, nz = field.shape # 计算网格步长 dx = self.grid_config['cell_size'] dy = dx dz = dx # 创建新场 new_field = field.copy() # 前向欧拉对流 for i in range(1, nx-1): for j in range(1, ny-1): for k in range(1, nz-1): # 获取速度 u = velocity_field['u'][i, j, k] v = velocity_field['v'][i, j, k] w = velocity_field['w'][i, j, k] # 计算追踪点 x_trace = i - u * dt / dx y_trace = j - v * dt / dy z_trace = k - w * dt / dz # 线性插值获取追踪点的场值 new_field[i, j, k] = self._trilinear_interpolate(field, x_trace, y_trace, z_trace) return new_field except Exception as e: print(f"✗ 欧拉对流失败: {e}") return field def _advect_rk2(self, field: np.ndarray, velocity_field: Dict[str, np.ndarray], dt: float) -> np.ndarray: """ 使用二阶Runge-Kutta方法对流场 Args: field: 要对流的场 velocity_field: 速度场 dt: 时间步长 Returns: 对流后的场 """ try: # RK2实现 # 中点方法: k1 = f(x_n, t_n), k2 = f(x_n + dt/2 * k1, t_n + dt/2) # x_{n+1} = x_n + dt * k2 return field # 简化处理 except Exception as e: print(f"✗ RK2对流失败: {e}") return field def _advect_semi_lagrangian(self, field: np.ndarray, velocity_field: Dict[str, np.ndarray], dt: float) -> np.ndarray: """ 使用半拉格朗日方法对流场 Args: field: 要对流的场 velocity_field: 速度场 dt: 时间步长 Returns: 对流后的场 """ try: # 半拉格朗日方法实现 # 从终点追踪到起点，然后插值 return field # 简化处理 except Exception as e: print(f"✗ 半拉格朗日对流失败: {e}") return field def _trilinear_interpolate(self, field: np.ndarray, x: float, y: float, z: float) -> float: """ 三线性插值 Args: field: 场数据 x, y, z: 插值点坐标 Returns: 插值结果 """ try: # 获取整数部分和小数部分 i = int(x) j = int(y) k = int(z) fx = x - i fy = y - j fz = z - k # 边界检查 nx, ny, nz = field.shape if i < 0 or i >= nx-1 or j < 0 or j >= ny-1 or k < 0 or k >= nz-1: return 0.0 # 三线性插值 c00 = field[i, j, k] * (1 - fx) + field[i+1, j, k] * fx c01 = field[i, j, k+1] * (1 - fx) + field[i+1, j, k+1] * fx c10 = field[i, j+1, k] * (1 - fx) + field[i+1, j+1, k] * fx c11 = field[i, j+1, k+1] * (1 - fx) + field[i+1, j+1, k+1] * fx c0 = c00 * (1 - fy) + c10 * fy c1 = c01 * (1 - fy) + c11 * fy return c0 * (1 - fz) + c1 * fz except Exception as e: return 0.0 def diffuse_field(self, field: np.ndarray, diffusion_coeff: float, dt: float) -> np.ndarray: """ 扩散场量 Args: field: 要扩散的场 diffusion_coeff: 扩散系数 dt: 时间步长 Returns: 扩散后的场 """ try: # 根据配置选择扩散格式 scheme = self.numerical_methods['diffusion_scheme'] if scheme == 'explicit': return self._diffuse_explicit(field, diffusion_coeff, dt) elif scheme == 'implicit': return self._diffuse_implicit(field, diffusion_coeff, dt) else: # 默认使用显式方法 return self._diffuse_explicit(field, diffusion_coeff, dt) except Exception as e: print(f"✗ 场扩散失败: {e}") return field def _diffuse_explicit(self, field: np.ndarray, diffusion_coeff: float, dt: float) -> np.ndarray: """ 显式扩散 Args: field: 要扩散的场 diffusion_coeff: 扩散系数 dt: 时间步长 Returns: 扩散后的场 """ try: # 拉普拉斯算子的显式更新 nx, ny, nz = field.shape dx = self.grid_config['cell_size'] # 扩散系数 alpha = diffusion_coeff * dt / (dx * dx) # 确保稳定性 (CFL条件) if alpha > 0.166667: # 1/6 for 3D alpha = 0.166667 new_field = field.copy() for i in range(1, nx-1): for j in range(1, ny-1): for k in range(1, nz-1): laplacian = ( field[i+1, j, k] + field[i-1, j, k] + field[i, j+1, k] + field[i, j-1, k] + field[i, j, k+1] + field[i, j, k-1] - 6 * field[i, j, k] ) new_field[i, j, k] = field[i, j, k] + alpha * laplacian return new_field except Exception as e: print(f"✗ 显式扩散失败: {e}") return field def _diffuse_implicit(self, field: np.ndarray, diffusion_coeff: float, dt: float) -> np.ndarray: """ 隐式扩散 Args: field: 要扩散的场 diffusion_coeff: 扩散系数 dt: 时间步长 Returns: 扩散后的场 """ try: # 隐式扩散需要求解线性系统 # (I - alpha * L) * field_new = field_old return field # 简化处理 except Exception as e: print(f"✗ 隐式扩散失败: {e}") return field def set_solver_config(self, config: Dict[str, Any]): """ 设置求解器配置 Args: config: 求解器配置字典 """ try: self.solver_config.update(config) print(f"✓ 求解器配置已更新: {self.solver_config}") except Exception as e: print(f"✗ 求解器配置设置失败: {e}") def get_solver_config(self) -> Dict[str, Any]: """ 获取求解器配置 Returns: 求解器配置字典 """ return self.solver_config.copy() def set_numerical_methods(self, methods: Dict[str, str]): """ 设置数值方法 Args: methods: 数值方法字典 """ try: self.numerical_methods.update(methods) print(f"✓ 数值方法已更新: {self.numerical_methods}") except Exception as e: print(f"✗ 数值方法设置失败: {e}") def get_numerical_methods(self) -> Dict[str, str]: """ 获取数值方法 Returns: 数值方法字典 """ return self.numerical_methods.copy() def get_performance_stats(self) -> Dict[str, Any]: """ 获取性能统计信息 Returns: 性能统计字典 """ return self.performance_stats.copy() def reset_performance_stats(self): """重置性能统计信息""" try: self.performance_stats = { 'pressure_solves': 0, 'velocity_updates': 0, 'total_solve_time': 0.0, 'average_solve_time': 0.0 } print("✓ 性能统计信息已重置") except Exception as e: print(f"✗ 性能统计信息重置失败: {e}") def get_solver_state(self) -> Dict[str, Any]: """ 获取求解器状态 Returns: 求解器状态字典 """ return self.solver_state.copy()