高级物理引擎实战教程

课程简介

本教程旨在帮助读者深入理解现代物理引擎的核心算法与实现技术。内容涵盖从基础的弹簧质点系统到前沿的物质点法(MPM)，从理论推导到性能优化，为读者提供完整的物理仿真知识体系。

预备知识

• 高等数学（微积分、线性代数）
• 基础物理学（牛顿力学、连续介质力学）
• 编程基础（Python或任一编程语言）

学习目标

• 掌握拉格朗日与欧拉视角的仿真方法
• 理解有限元、有限差分等数值方法
• 熟练运用Taichi编程语言进行物理仿真
• 具备设计和优化物理引擎的能力

第一章：导论

1.1 基于物理的动画简介

• 1.1.1 物理引擎的定义与应用领域：计算机软件提供物理系统近似仿真，包括刚体动力学、软体动力学和流体动力学
• 1.1.2 计算机图形学中的物理仿真：从影视特效(如迪士尼动画)到游戏引擎(Angry Birds, Besieged)的演进历程
• 1.1.3 工程CAE与游戏物理的差异：精度要求、实时性约束、稳定性与艺术效果的平衡
• 1.1.4 实时性与精确性的权衡：时间步长选择、简化模型、LOD技术

1.2 Taichi编程语言简介

• 1.2.1 安装与环境配置：pip install taichi，支持Windows/Linux/OSX，CPU/GPU后端选择
• 1.2.2 基本数据类型与张量操作：ti.i32/f32/f64，ti.var/Vector/Matrix，shape参数定义维度
• 1.2.3 核函数(kernel)与函数(func)：@ti.kernel装饰器，类型提示，强制内联，编译时优化
• 1.2.4 并行for循环与原子操作：自动并行化，ti.atomic_add，避免竞态条件

1.3 Taichi的自动并行化

• 1.3.1 并行计算基础概念：数据并行vs任务并行，SIMD/SIMT执行模型
• 1.3.2 Range-for与Struct-for循环：ti.ndrange多维迭代，稀疏张量遍历，最外层作用域自动并行
• 1.3.3 线程安全与竞态条件：原子操作保证，增强赋值自动原子化(x[i] += 1)
• 1.3.4 性能分析与调优基础：ti.profiler，kernel时间测量，内存带宽分析

1.4 Taichi程序的调试技巧

• 1.4.1 调试模式与边界检查：ti.init(debug=True, arch=ti.cpu)启用边界检查器
• 1.4.2 打印调试与断言：kernel内print支持，assert语句，运行时错误捕获
• 1.4.3 常见错误与解决方案：越界访问、类型不匹配、未初始化变量
• 1.4.4 性能瓶颈识别：内存访问模式分析、核函数粒度调整

1.5 面向数据的编程(DOP)

• 1.5.1 AoS vs SoA布局：结构数组vs数组结构，缓存局部性影响
• 1.5.2 缓存友好的数据组织：连续内存访问，避免跨步访问
• 1.5.3 内存访问模式优化：合并访问、预取策略
• 1.5.4 向量化编程技巧：SIMD指令利用，数据对齐

1.6 面向对象的编程(OOP)

• 1.6.1 Taichi中的类与装饰器：@ti.data_oriented类装饰器，成员函数kernel化
• 1.6.2 ODOP设计模式：Objective Data-Oriented Programming，结合OOP与DOP优势
• 1.6.3 太阳系仿真案例：行星类设计，引力计算，轨道积分
• 1.6.4 代码复用与模块化：类继承、组合模式、接口设计

1.7 元编程(MP)

• 1.7.1 模板核函数：@ti.template()参数化kernel，编译时实例化
• 1.7.2 维度无关编程：ti.static(range(dim))，2D/3D代码统一
• 1.7.3 编译时分支与循环展开：ti.static(if/for)，零开销抽象
• 1.7.4 静态类型推导：编译时类型检查，自动类型转换规则

第二章：拉格朗日视角（1）

2.1 弹簧质点系统

• 2.1.1 胡克定律与牛顿第二定律：f_ij = -k(||x_i - x_j|| - l_ij)(x_i - x_j)/||x_i - x_j||，dv_i/dt = f_i/m_i
• 2.1.2 弹簧刚度与阻尼系数：刚度k决定系统响应速度，阻尼系数控制能量耗散，临界阻尼ζ=1
• 2.1.3 系统的能量守恒：动能T=½mv²，势能U=½kx²，数值积分的能量漂移问题
• 2.1.4 多质点系统的矩阵表示：质量矩阵M，刚度矩阵K，阻尼矩阵C，状态空间形式

2.2 布料模拟

• 2.2.1 结构弹簧、剪切弹簧与弯曲弹簧：结构弹簧维持基本形状，剪切弹簧抵抗剪切变形，弯曲弹簧模拟弯曲刚度
• 2.2.2 外力：重力mg，风力基于相对速度和法向量，碰撞力基于穿透深度
• 2.2.3 约束处理：固定点设置速度为0，滑动约束投影速度到切向量
• 2.2.4 自碰撞检测与响应：空间哈希加速，连续碰撞检测(CCD)，冲量响应

2.3 显式与隐式时间积分器

• 2.3.1 前向欧拉法：v_{t+1} = v_t + Δt·f_t/m，x_{t+1} = x_t + Δt·v_t，条件稳定Δt ≤ c√(m/k)
• 2.3.2 辛欧拉法：v_{t+1} = v_t + Δt·f_t/m，x_{t+1} = x_t + Δt·v_{t+1}，保持相空间体积
• 2.3.3 后向欧拉法：v_{t+1} = v_t + Δt·f(x_{t+1})/m，需要求解非线性系统，无条件稳定
• 2.3.4 中点法与RK4：β=0显式，β=0.5中点法，β=1后向欧拉，RK4四阶精度

2.4 隐式积分的线性化与求解

• 2.4.1 牛顿法线性化：f(x_{t+1}) ≈ f(x_t) + ∂f/∂x·Δx，一步Newton迭代
• 2.4.2 雅可比矩阵的构建：[I - Δt²M⁻¹∂f/∂x]v_{t+1} = v_t + ΔtM⁻¹f(x_t)
• 2.4.3 Jacobi与Gauss-Seidel迭代：Jacobi并行友好，G-S收敛更快，松弛因子ω选择
• 2.4.4 共轭梯度法简介：Krylov子空间方法，预条件加速收敛

2.5 光滑粒子流体动力学(SPH)

• 2.5.1 核函数与粒子近似：A(x) = Σ_j A_j·m_j/ρ_j·W(||x-x_j||, h)，Cubic Spline核
• 2.5.2 密度计算与状态方程：ρ_i = Σ_j m_j·W(||x_i-x_j||, h)，p = B((ρ/ρ_0)^γ - 1)
• 2.5.3 压力梯度与粘性力：∇p使用对称形式保证动量守恒，人工粘性稳定数值解
• 2.5.4 表面张力模型：基于颜色场梯度，CSF方法，表面法向量计算

2.6 基于位置的流体(PBF)

• 2.6.1 位置约束与拉格朗日乘子：C(x) = ρ/ρ_0 - 1 = 0，λ通过Newton迭代求解
• 2.6.2 不可压缩性约束：密度约束投影，位置修正Δx = -λ∇C
• 2.6.3 人工压力与涡量补偿：s_corr项防止粒子聚集，涡量约束恢复旋转运动
• 2.6.4 PCISPH与DFSPH方法：预测-校正策略，散度自由条件，压力求解优化

2.7 体素化：从三角网格生成粒子

• 2.7.1 点在多边形内测试：射线法，绕数(winding number)，角度和方法
• 2.7.2 射线-三角形相交：Möller-Trumbore算法，重心坐标判断
• 2.7.3 有符号距离场(SDF)：内部负值外部正值，快速行进法(FMM)构建
• 2.7.4 自适应采样策略：基于曲率的采样密度，蓝噪声分布，Poisson disk采样

2.8 快速邻居搜索

• 2.8.1 均匀网格法：O(1)查询复杂度，网格大小=支持半径，多重网格处理不同尺度
• 2.8.2 紧凑哈希(Compact Hashing)：Z-index空间哈希，处理稀疏分布，内存效率
• 2.8.3 Z-order与空间填充曲线：Morton编码，保持空间局部性，缓存友好
• 2.8.4 并行邻居搜索算法：原子操作构建，避免竞态，GPU实现优化

2.9 刚体模拟简介

• 2.9.1 刚体运动学：位置x，四元数姿态q，角速度ω，速度扭曲(twist)表示
• 2.9.2 刚体动力学：力F=ma，力矩τ=Iα，惯性张量I，Euler方程
• 2.9.3 碰撞检测：GJK距离查询，SAT分离轴定理，宽相位与窄相位
• 2.9.4 碰撞响应：冲量J计算，恢复系数e，摩擦锥约束，LCP求解

第三章：拉格朗日视角（2）：有限元仿真

3.1 弱形式与拉格朗日有限元入门

• 3.1.1 强形式vs弱形式：∇·∇u = 0 ⟺ ∀w, ∫∫ w(∇·∇u)dA = 0，将PDE转换为积分形式
• 3.1.2 分部积分与散度定理：∇w·∇u + w∇·∇u = ∇·(w∇u)，利用散度定理∫∫ ∇·(w∇u)dA = ∮ w∇u·dn
• 3.1.3 Galerkin方法：FEM属于Galerkin方法族，使用测试函数将连续PDE离散为线性系统
• 3.1.4 试函数与测试函数：使用分片线性基函数φᵢ(x)作为测试函数，u(x) = Σ uⱼφⱼ(x)

3.2 变形与弹性基础

• 3.2.1 变形梯度张量F：F := ∂x_deformed/∂x_rest，描述局部变形，具有平移不变性
• 3.2.2 格林应变与柯西应变：体积比J = det(F)，右柯西-格林张量C = F^T F
• 3.2.3 应力张量：第一Piola-Kirchhoff应力P(F) = ∂ψ/∂F，柯西应力σ，关系式τ = PF^T，σ = J₀F^(-T)P
• 3.2.4 本构关系与材料参数：杨氏模量E，泊松比ν∈[0, 0.5)，Lamé参数λ和μ，转换公式K = E/[3(1-2ν)]

3.3 超弹性材料模型

• 3.3.1 应变能密度函数：ψ = ψ(F)，超弹性材料的应力-应变关系由应变能密度函数定义
• 3.3.2 Neo-Hookean模型：ψ(F) = μ/2 Σᵢ[(F^T F)ᵢᵢ - 1] - μlog(J) + λ/2 log²(J)，P(F) = μ(F - F^(-T)) + λlog(J)F^(-T)
• 3.3.3 Corotated模型：ψ(F) = μΣᵢ(σᵢ - 1)² + λ/2(J - 1)²，其中σᵢ为F的奇异值，P(F) = 2μ(F - R) + λ(J - 1)JF^(-T)
• 3.3.4 Mooney-Rivlin模型：用于橡胶类材料，包含不变量I₁和I₂的高阶项

3.4 基于六面体网格的拉格朗日有限元

• 3.4.1 等参单元与形函数：x_deformed = Fx_rest + p，线性单元内F为常数，使用三线性形函数
• 3.4.2 数值积分：高斯求积计算Kᵢⱼ = ∫∫ ∇φᵢ·∇φⱼ dA，通常使用2×2×2高斯积分点
• 3.4.3 单元刚度矩阵组装：8×8单元矩阵，每节点2个自由度(2D)或3个自由度(3D)
• 3.4.4 边界条件处理：Dirichlet边界直接设置uᵢ = f(xᵢ)，Neumann边界修改载荷向量f

3.5 基于四面体网格的拉格朗日有限元

• 3.5.1 线性四面体单元：F = DB，其中B = [a_rest-c_rest, b_rest-c_rest]^(-1)，D = [a_deformed-c_deformed, b_deformed-c_deformed]
• 3.5.2 常应变假设：四面体内应变场均匀，适合大变形但可能过于刚硬
• 3.5.3 体积锁定与缓解策略：不可压缩材料(ν→0.5)时出现，使用混合有限元或降阶积分缓解
• 3.5.4 网格生成与质量控制：Delaunay三角剖分，控制最小角度和纵横比

3.6 可逆(Invertible)有限元

• 3.6.1 单元翻转问题：当det(F) ≤ 0时单元发生翻转，导致数值不稳定和非物理结果
• 3.6.2 SVD分解与修复：F = UΣV^T，检测并修正负奇异值，保证det(U) = det(V) = 1
• 3.6.3 能量势垒方法：在应变能中添加势垒项防止单元翻转，确保J > 0
• 3.6.4 稳定性保证：使用投影方法将变形梯度约束到可行域内

3.7 拓扑优化

• 3.7.1 最小柔度问题：min L(ρ) = u^T K(ρ)u，s.t. K(ρ)u = f，Σρₑ ≤ cV，ρₑ ∈ [ρₘᵢₙ, 1]
• 3.7.2 SIMP方法：固体各向同性材料惩罚(Solid Isotropic Material with Penalization)，E(ρ) = ρ^p E₀
• 3.7.3 敏感度分析：计算目标函数对设计变量的梯度，∂L/∂ρₑ = -u^T ∂K/∂ρₑ u
• 3.7.4 优化准则法(OC)：基于拉格朗日乘子的迭代更新，ρₑ^(new) = max(ρₘᵢₙ, min(1, ρₑB_e^η))

3.8 高级Taichi特性（1）

• 3.8.1 可微编程基础：@ti.data_oriented类装饰器，使用ODOP模式结合OOP和DOP
• 3.8.2 自动微分：前向与反向模式，ti.template()参数化kernel，支持编译时实例化
• 3.8.3 梯度计算与优化：needs_grad=True声明梯度张量，使用ti.Tape(loss)计算梯度
• 3.8.4 物理过程的端到端优化：forces Fᵢ = -∂U(x)/∂xᵢ，利用AutoDiff避免手动求导

第四章：欧拉视角（1）

4.1 欧拉vs拉格朗日描述

• 4.1.1 物质导数与对流项：D/Dt := ∂/∂t + u·∇，欧拉视角的物质导数包含对流项u·∇
• 4.1.2 不可压Navier-Stokes方程：ρDu/Dt = -∇p + μ∇²u + ρg，∇·u = 0，v = μ/ρ为运动粘度
• 4.1.3 算子分裂方法：将N-S方程分解为对流、外力、投影三个步骤分别求解
• 4.1.4 稳定性与CFL条件：显式时间积分需满足CFL条件Δt ≤ CΔx/|u|_max，确保数值稳定

4.2 网格类型与离散化

• 4.2.1 同位网格vs交错网格(Staggered Grid)：同位网格所有量存储在网格中心，MAC网格速度分量存储在面中心
• 4.2.2 MAC网格的优势：自然满足散度自由条件，边界条件处理更简单，避免数值振荡
• 4.2.3 双线性插值：四个角点的加权平均，f(x,y) = f₁₁(1-α)(1-β) + f₂₁α(1-β) + f₁₂(1-α)β + f₂₂αβ
• 4.2.4 散度与梯度算子离散：中心差分(∇·u)ᵢ,ⱼ = (uᵢ₊₁,ⱼ - uᵢ,ⱼ)/Δx + (vᵢ,ⱼ₊₁ - vᵢ,ⱼ)/Δy

4.3 半拉格朗日输送(Semi-Lagrangian Advection)

• 4.3.1 反向粒子追踪：从网格点向后追踪找到上一时刻的位置，x_prev = x - Δt·u(x)
• 4.3.2 速度插值策略：使用RK2/RK3提高精度，p_mid = p - 0.5Δt·velocity(p)，p_prev = p - Δt·velocity(p_mid)
• 4.3.3 数值耗散问题：基本半拉格朗日方法具有较大的数值粘性，需要高阶方法改善
• 4.3.4 单调性保持：使用限制器防止非物理的过冲和欠冲，保持数值稳定

4.4 高阶输送格式

• 4.4.1 MacCormack方法：x = SL(x, Δt)，x = SL(x, -Δt)，x_final = x + 0.5(x - x*)，二阶精度
• 4.4.2 BFECC方法：前后误差补偿校正(Back and Forth Error Compensation and Correction)，估计误差并补偿
• 4.4.3 Runge-Kutta积分：RK3: v₁=velocity(p)，v₂=velocity(p-0.5Δt·v₁)，v₃=velocity(p-0.75Δt·v₂)
• 4.4.4 限制器与夹持(Clamping)：防止非物理的过冲值，对BFECC结果进行夹持到原始值范围内

4.5 Chorin式压力投影

• 4.5.1 Helmholtz-Hodge分解：任意向量场可分解为无散场和无旋场之和，u = u_div-free + ∇φ
• 4.5.2 压力泊松方程推导：∇·u^(n+1) = 0 ⇒ ∇·(u - Δt/ρ ∇p) = 0 ⇒ ∇²p = ρ/Δt ∇·u
• 4.5.3 离散化与边界条件：5点Laplace模板[-4 1; 1 -4 1; 1 -4]，Dirichlet和Neumann边界条件
• 4.5.4 速度修正步：u^(n+1) = u* - Δt/ρ ∇p，确保速度场满足不可压缩条件

4.6 固体边界条件

• 4.6.1 无滑移条件：v·n = 0，v·t = 0，流体在壁面速度为零(黏性)
• 4.6.2 自由滑移条件：v·n = 0，仅约束法向速度，切向可以自由滑动
• 4.6.3 浸入边界法：在流体中嵌入固体边界，使用分布函数表示边界力
• 4.6.4 切割单元法：将相交边界的网格单元切割，精确处理复杂几何边界

4.7 自由表面边界条件

• 4.7.1 运动学条件：表面上的粒子始终在表面上，Dφ/Dt = 0，即∂φ/∂t + u·∇φ = 0
• 4.7.2 动力学条件：表面上压力连续，考虑表面张力的影响σ·n = γκ n
• 4.7.3 表面张力处理：使用CSF模型(Continuum Surface Force)将表面张力转换为体积力
• 4.7.4 Ghost Fluid方法：在不同相中使用虚拟值来处理边界条件，保持数值精度

第五章：欧拉视角（2）：线性系统求解器

5.1 稀疏矩阵与零空间(Nullspaces)

• 5.1.1 稀疏矩阵存储格式：CSR(Compressed Sparse Row)、COO(Coordinate)格式，仅存储非零元素
• 5.1.2 兼容性条件：泊松方程有解的必要条件，右端项必须垂直于零空间
• 5.1.3 零空间投影：对于奇异系统，需要将解投影到零空间的正交补空间
• 5.1.4 奇异系统的处理：添加约束条件或使用广义逆矩阵(Moore-Penrose pseudoinverse)

5.2 Krylov子空间求解器

• 5.2.1 Krylov子空间理论：K_m(A,r₀) = span{r₀, Ar₀, A²r₀, ..., A^(m-1)r₀}，在子空间中寻找最优解
• 5.2.2 共轭梯度法(CG)：适用于对称正定矩阵(SPD)，α_k = r_k^T r_k / (p_k^T A p_k)，β_k = r_{k+1}^T r_{k+1} / (r_k^T r_k)
• 5.2.3 BiCGSTAB方法：适用于不对称矩阵，双共轭梯度稳定化方法，收敛性更稳定
• 5.2.4 收敛性分析：收敛速度与条件数κ(A) = λ_max/λ_min相关，条件数越小收敛越快

5.3 预条件(Preconditioning)

• 5.3.1 预条件的作用：找到接近A但更易求逆的矩阵M，使M^(-1)A有更小的条件数
• 5.3.2 Jacobi预条件：M = diag(A)，最简单的对角预条件，并行友好但效果有限
• 5.3.3 不完全Cholesky分解：A ≈ LL^T，但只保持A的稀疏模式，减少存储和计算量
• 5.3.4 修正不完全Cholesky(MIC)：保持行和不变的MIC(0)，改善数值稳定性

5.4 多重网格方法

• 5.4.1 误差的频率分析：高频误差由Jacobi/Gauss-Seidel快速消除，低频误差需要粗网格处理
• 5.4.2 限制与延拓算子：限制(Restriction)I_{2h}^h将细网格映射到粗网格，延拓(Prolongation)I_h^{2h}将粗网格映射到细网格
• 5.4.3 V-cycle与W-cycle：V-cycle单次递归，W-cycle双次递归，F-cycle可变递归，每种cycle有不同的效率和成本
• 5.4.4 粗网格修正：在粗网格上求解残差方程，然后延拓到细网格作为校正

5.5 几何多重网格

• 5.5.1 网格层次构建：通过粗化几何网格生成层次结构，通常每个维度粗化1/2
• 5.5.2 Galerkin粗化：A_{2h} = I_{2h}^h A_h I_h^{2h}，保证粗网格算子的数学性质
• 5.5.3 光滑器选择：Jacobi、Gauss-Seidel、Red-Black Gauss-Seidel，选择不同光滑器影响收敛性
• 5.5.4 并行多重网格：在多核/GPU上并行化实现，需要处理负载均衡和通信优化

5.6 代数多重网格(AMG)

• 5.6.1 强连接与粗网格选择：基于矩阵系数大小判断强连接|a_{ij}| ≥ θ max_k |a_{ik}|，选择粗网格点
• 5.6.2 插值算子构造：构造延拓算子I_h^{2h}，使网格间传递信息最少失真
• 5.6.3 经典AMG与平滑聚合：经典AMG使用C/F分裂，平滑聚合(Smoothed Aggregation)使用聚合方法
• 5.6.4 自适应AMG：根据矩阵性质自动选择粗化策略和插值算子

5.7 无矩阵(Matrix-free)方法

• 5.7.1 矩阵-向量乘积的隐式计算：不存储矩阵A，直接计算Av，利用空间局部性清除缓存
• 5.7.2 内存带宽优化：现代架构中FLOPs免费但内存带宽昂贵，重新计算比读取更快
• 5.7.3 算子融合技术：将多个算子操作融合在一个kernel中，减少内存访问次数
• 5.7.4 GPU实现策略：充分利用GPU的并行计算能力和内存层次，优化线程占用率

5.8 泊松方程的快速解法

• 5.8.1 格林函数与基本解：φ(x) = -1/(4π||x-y||) ∫ρ(y)dy，泊松方程的基本解
• 5.8.2 快速多极子方法(FMM)：O(N)algirithm，使用多极子展开和局部展开(M2M, M2L, L2L)
• 5.8.3 PPPM方法：粒子-粒子粒子-网格(Particle-Particle Particle-Mesh)，近远场分解
• 5.8.4 FFT方法(周期边界)：利用快速傅里叶变换解周期边界条件下的泊松方程，O(N log N)

第六章：高级输送格式与等势面方法

6.1 高阶输送格式深入

• 6.1.1 WENO格式：5阶WENO使用3个子模板，权重基于光滑度指标β_k，避免跨越不连续处
• 6.1.2 TVD限制器：Total Variation Diminishing条件TV(u^{n+1}) ≤ TV(u^n)，保证无新的极值产生
• 6.1.3 通量限制器：minmod、superbee、van Leer限制器，φ(r)满足TVD条件0≤φ≤2
• 6.1.4 特征线方法：沿特征线dx/dt=u求解，拉格朗日观点在欧拉网格上的应用

6.2 有符号距离场(SDF)

• 6.2.1 距离场的定义与性质：φ(x)=±min||x-x_surface||，内部负外部正，|∇φ|=1
• 6.2.2 快速行进法(FMM)：解Eikonal方程|∇T|=1/F，使用堆结构O(N log N)复杂度
• 6.2.3 快速扫描法(FSM)：Gauss-Seidel迭代8个扫描方向，O(N)复杂度但常数大
• 6.2.4 重新初始化：∂φ/∂τ + sign(φ₀)(|∇φ|-1)=0，保持零等势面位置不变

6.3 等势面(Level Set)方法

• 6.3.1 界面演化方程：∂φ/∂t + u·∇φ = 0，Hamilton-Jacobi型方程，使用HJ-WENO格式
• 6.3.2 速度延拓：∂u_ext/∂τ + sign(φ)n·∇u_ext = 0，将界面速度延拓到整个域
• 6.3.3 质量守恒问题：Level Set方法不守恒，需要粒子校正或VOF耦合
• 6.3.4 粒子等势面(PLS)：在界面附近撒粒子，根据粒子位置修正level set

6.4 自由表面追踪

• 6.4.1 VOF方法对比：Volume of Fluid保证质量守恒，但界面重构复杂，难以计算曲率
• 6.4.2 界面重构技术：SLIC(Simple Line Interface Calculation)，PLIC(Piecewise Linear)
• 6.4.3 表面张力计算：κ=∇·n=∇·(∇φ/|∇φ|)，CSF模型F_st=σκδ(φ)∇φ
• 6.4.4 拓扑变化处理：Level Set自然处理合并分离，VOF需要特殊处理

6.5 渲染技术

• 6.5.1 路径追踪(Path Tracing)：蒙特卡洛光线追踪，渲染方程L_o=L_e+∫f_r·L_i·cosθdω
• 6.5.2 球面追踪(Sphere Tracing)：使用SDF加速光线行进，步长=到最近表面距离
• 6.5.3 行军立方体(Marching Cubes)：256种立方体配置，通过查表生成三角网格
• 6.5.4 运动模糊与景深：时间积分运动轨迹，光圈模拟实现景深效果

6.6 体素渲染

• 6.6.1 数字微分分析器(DDA)：3D Bresenham算法，高效遍历光线经过的体素
• 6.6.2 光线行进算法：Ray marching with adaptive步长，early termination优化
• 6.6.3 体积渲染方程：L(s)=∫τ(s,t)·σ_a(t)·L_e(t)dt，τ=exp(-∫σ_t dt)透射率
• 6.6.4 实时体渲染优化：空间跳跃(empty space skipping)，预积分传输函数，GPU加速

第七章：混合欧拉-拉格朗日视角（1）

7.1 混合方法概述

• 7.1.1 拉格朗日与欧拉的优缺点：拉格朗日善于对流但投影困难，欧拉善于投影但对流有耗散
• 7.1.2 物理量守恒问题：动量、角动量、体积(不可压缩性)、能量(低耗散)的守恒
• 7.1.3 数值耗散与噪声：欧拉方法的耗散vs拉格朗日方法的噪声，需要平衡权衡
• 7.1.4 混合策略设计：使用粒子追踪材料信息，使用网格计算力场(投影)

7.2 粒子元胞法(PIC)

• 7.2.1 P2G与G2P传输：P2G: grid_v[i] += weight * v[p]，G2P: new_v += weight * grid_v[i]，使用B样条核函数
• 7.2.2 核函数选择：线性、二次、三次B样条，高阶函数更平滑但成本更高
• 7.2.3 信息损失分析：18个网格自由度vs 2个粒子自由度，信息不匹配导致耗散
• 7.2.4 动量守恒性：使用加权速度传输v_i = Σ_p m_p v_p W_{ip} / m_i保证动量守恒

7.3 流体隐粒子(FLIP)

• 7.3.1 增量式速度更新：v_p^{n+1} = v_p^n + gather(v^{n+1} - v^n)，传输的是速度变化量而非绝对值
• 7.3.2 噪声vs耗散权衡：PIC过度耗散但稳定，FLIP保持能量但噪声大，需要权衡
• 7.3.3 PIC/FLIP混合：v_p^{n+1} = (1-α) * PIC + α * FLIP，通常α=0.95-0.99
• 7.3.4 自适应混合比例：根据局部流动特征动态调整混合系数α

7.4 仿射粒子元胞法(APIC)

• 7.4.1 仿射速度场：v(x) = v_p + C_p(x - x_p)，粒子周围保持线性速度分布
• 7.4.2 角动量守恒：APIC在P2G和G2P过程中保持角动量守恒，已有数学证明
• 7.4.3 C矩阵的物理意义：C_p表示粒子周围速度梯度，类似于变形梯度但用于速度场
• 7.4.4 APIC的稳定性：相比PIC/FLIP更稳定，不需要额外的速度存储，是现代MPM的基础

7.5 PolyPIC方法

• 7.5.1 高阶多项式表示：18模式=9网格点×2自由度，实现无损传输，包含二次项
• 7.5.2 Taylor展开传输：使用高阶Taylor展开表示粒子周围的速度场，最大化信息保留
• 7.5.3 精度与效率分析：高精度但计算成本高，内存开销大，适合对精度要求极高的应用
• 7.5.4 实现细节：需要处理大量的多项式系数和复杂的传输计算，实现复杂度较高

7.6 粒子-网格传输细节

• 7.6.1 B样条核函数：线性N(x)=1-|x|，二次N(x)=3/4-(x-1)²，三次N(x)=(2-|x|)³/6
• 7.6.2 二次vs三次插值：二次更平滑但计算量大，三次精度更高但支持域更大
• 7.6.3 核函数的紧支性：插值核在有限支持域内非零，影响计算效率和数值精度
• 7.6.4 并行传输算法：使用原子操作避免竞态条件，GPU实现需要优化内存访问模式

第八章：混合欧拉-拉格朗日视角（2）：物质点法

8.1 物质点法(MPM)基础

• 8.1.1 MPM的历史与发展：Sulsky & Schreyer 1996发明，2013年Stomakhin等引入图形学
• 8.1.2 与FEM的关系：MPM属于无单元Galerkin(EFG)，粒子对应FEM积分点
• 8.1.3 弱形式推导：使用弱形式推导方程，通常使用单点积分规则
• 8.1.4 空间离散化：粒子携带材料信息，网格用于碰撞处理和动量更新

8.2 经典MPM算法

• 8.2.1 应力更新(USL vs USF)：Update Stress Last vs Update Stress First，影响数值精度和稳定性
• 8.2.2 形函数选择：使用B样条基函数，通常三次或二次，影响光滑性和计算效率
• 8.2.3 积分点与背景网格：粒子作为积分点，背景网格为计算网格，可以是SPGrid或OpenVDB
• 8.2.4 边界条件处理：粘性边界(sticky)、滑动边界(slip)、分离边界(separate)

8.3 移动最小二乘MPM(MLS-MPM)

• 8.3.1 MLS插值理论：使用MLS形函数替代B样条，直接复用APIC传输方案
• 8.3.2 88行实现解析：F^{n+1} = (I+ΔtC^n)F^n，变形更新在P2G阶段，J[p] *= 1 + dt * new_C.trace()
• 8.3.3 性能优势分析：FLOPs减半因为直接使用C_p作为∇v的近似，不需计算∇W_{ip}
• 8.3.4 与APIC的关系：MLS-MPM基于APIC，在弹性力计算中融合了APIC动量传输

8.4 MPM中的本构模型

• 8.4.1 弹性固体(Neo-Hookean, Corotated)：使用超弹性材料模型，P(F) = 2μ(F-R) + λ(J-1)JF^{-T}
• 8.4.2 弱可压缩流体：状态方程 p = K(1-J)，直接维护J而非完整的F矩阵避免数值不稳定
• 8.4.3 弹塑性材料：F_p = F_{p,elastic} × F_{p,plastic}，弹性部分贡献势能
• 8.4.4 屈服准则与塑性流动：Box屈服准则、Drucker-Prager准则、Cam-clay模型，用于SVD夹持

8.5 MPM中的拉格朗日力

• 8.5.1 弹簧与阻尼器：将MPM粒子视为FEM顶点，使用FEM势能模型，需要三角网格连接
• 8.5.2 薄壳与薄膜：处理低维结构在高维空间中的变形，需要特殊的投影和积分方法
• 8.5.3 刚体约束：通过约束残差实现刚体连接，需要求解约束优化问题
• 8.5.4 接触力模型：处理粒子间或粒子-边界间的接触，可能需要额外的碰撞检测

8.6 数值断裂

• 8.6.1 连续损伤力学(CDM)：通过损伤变量描述材料退化，在本构模型中引入损伤项
• 8.6.2 相场断裂模型：使用相场变量描述裂纹的演化，需要解额外的相场方程
• 8.6.3 CPIC方法：Compatible PIC，保持网格连接性同时允许断裂
• 8.6.4 断裂准则：基于应力、应变或能量的断裂判断条件，如最大主应力准则

8.7 隐式MPM

• 8.7.1 隐式时间积分：处理刚性材料时需要隐式方法保证稳定性，解非线性系统
• 8.7.2 切线刚度矩阵：计算Hessian矩阵∂²U/∂xᵢ∂xⱼ，需要考虑几何非线性
• 8.7.3 Newton-Raphson迭代：x^{k+1} = x^k - J^{-1}F(x^k)，需要线搜索保证收敛
• 8.7.4 大变形处理：处理几何非线性和材料非线性，可能需要更新参考配置

8.8 高级Taichi特性（2）

• 8.8.1 稀疏数据结构设计：SPGrid和OpenVDB类型的层次稀疏网格，优化存储和访问效率
• 8.8.2 动态内存分配：ti.field替代ti.var，支持动态shape调整和内存管理
• 8.8.3 层次化网格：多层级稀疏数据结构，支持自适应精化和粗化
• 8.8.4 自定义数据布局：AoS vs SoA布局选择，优化内存访问模式和缓存局部性

第九章：高性能计算

9.1 处理器微架构

• 9.1.1 CPU流水线与超标量：取指、译码、执行、回写四阶段，超标量允许多发射
• 9.1.2 SIMD指令集：SSE/AVX/AVX-512，单指令多数据，向量化加速数值计算
• 9.1.3 分支预测与投机执行：现代CPU预测分支走向减少流水线停滞
• 9.1.4 乱序执行：CPU可以打乱指令顺序执行以提高效率，但保证程序语义

9.2 内存层级

• 9.2.1 缓存结构与大小：L1(32KB)、L2(256KB)、L3(8MB)三级缓存，访问时间逐级增加
• 9.2.2 缓存行与伪共享：缓存行64字节，伪共享会导致性能下降，需要数据对齐
• 9.2.3 TLB与页表：转换查找缓冲器加速虚拟地址翻译，减少页表访问
• 9.2.4 NUMA架构：非统一内存访问，本地内存vs远程内存访问延迟不同

9.3 性能分析与优化

• 9.3.1 Roofline模型：性能上限 = min(峰值计算能力, 运算强度 × 带宽)，指导优化方向
• 9.3.2 带宽vs计算瓶颈：识别瓶颈是内存带宽限制还是计算能力限制
• 9.3.3 性能计数器：PMU(Performance Monitoring Unit)监控指令、缓存、内存事件
• 9.3.4 Profile工具使用：VTune、perf、Instruments等工具分析性能瓶颈

9.4 并行编程模型

• 9.4.1 共享内存并行(OpenMP)：#pragma omp parallel for简化多线程编程，支持循环并行化
• 9.4.2 分布式内存(MPI)：消息传递接口，MPI_Send/MPI_Recv实现进程间通信
• 9.4.3 任务并行vs数据并行：任务并行适用于异构任务，数据并行适用于相同操作
• 9.4.4 负载均衡策略：静态分配、动态分配、工作窃取(work stealing)等策略

9.5 GPU编程

• 9.5.1 GPU架构概述：与多核并行的众核并行，SIMT(Single Instruction Multiple Thread)架构
• 9.5.2 线程组织：Grid, Block, Warp：Grid包含Block数组，Block包含线程Warp(32线程)
• 9.5.3 内存合并访问：相邻线程访问相邻地址可合并为一次128字节事务
• 9.5.4 占用率优化：最大化活跃warp数量隐藏内存延迟，考虑寄存器和共享内存限制

9.6 稀疏数据结构

• 9.6.1 SPGrid与OpenVDB：SPGrid使用虚拟内存管理，OpenVDB使用B+树层次结构
• 9.6.2 分层稀疏网格：多层级结构平衡存储和访问效率，支持自适应精化
• 9.6.3 位掩码与指针结构：位掩码标记有效区域，指针结构实现快速访问
• 9.6.4 动态拓扑更新：追踪活跃区域变化，动态分配和释放内存块

9.7 算法优化技术

• 9.7.1 循环优化：循环展开减少分支开销(unroll factor 4-8)，循环阻塞提高缓存命中率(tile size=√cache)，循环融合减少内存访问
• 9.7.2 数据重排与预取：AoS转SoA提高SIMD效率，软件预取__builtin_prefetch隐藏内存延迟，数据打包减少传输
• 9.7.3 向量化技巧：使用intrinsics(_mm256_add_ps)，数据对齐(attribute((aligned(32))))，避免gather/scatter操作
• 9.7.4 通信隐藏：计算与通信重叠，使用非阻塞MPI(MPI_Isend/Irecv)，双缓冲技术pipeline并行

9.8 多GPU并行

• 9.8.1 域分解策略：1D/2D/3D分解选择，最小化表面积/体积比，考虑负载均衡和通信开销
• 9.8.2 Halo区域交换：Ghost cell同步，使用GPUDirect P2P传输，重叠计算和通信减少等待
• 9.8.3 异步传输与计算重叠：CUDA流(streams)管理，使用cudaMemcpyAsync，内核执行与数据传输并行
• 9.8.4 负载动态迁移：工作窃取算法，基于性能监控的动态负载均衡，考虑迁移开销vs收益

第十章：可微编程与机器学习

10.1 可微仿真基础

• 10.1.1 端到端优化：从输入参数到最终损失函数的完整求导链，支持复杂物理系统优化
• 10.1.2 参数敏感度分析：计算∂L/∂θ分析参数对结果的影响，指导参数调优
• 10.1.3 伴随方法：解伴随方程计算梯度，适用于大规模优化问题
• 10.1.4 Checkpointing技术：在前向过程中保存关键状态，减少反向传播的内存需求

10.2 自动微分实现

• 10.2.1 计算图构建：构建计算图记录前向传播的计算过程，支持复杂控制流
• 10.2.2 前向vs反向传播：前向模式空间复杂度O(n)，反向模式时间复杂度低
• 10.2.3 高阶导数：递归应用微分规则计算Hessian矩阵，支持任意阶导数
• 10.2.4 混合模式AD：结合前向和反向模式的优势，针对不同问题选择最优策略

10.3 优化算法

• 10.3.1 梯度下降变种：Adam、RMSprop、AdaGrad等自适应学习率算法
• 10.3.2 拟牛顿方法：BFGS、L-BFGS估计Hessian逆矩阵，二阶收敛性
• 10.3.3 约束优化：活跃集方法、内点法、罚函数方法处理约束
• 10.3.4 随机优化：随机梯度下降(SGD)、遗传算法、模拟退火等参数估计

10.4 逆问题求解

• 10.4.1 参数识别：从观测数据推断材料参数、边界条件等物理量
• 10.4.2 初始条件优化：优化初始速度、位置分布以匹配目标轨迹
• 10.4.3 控制问题：设计控制信号使系统达到期望状态，如机器人运动控制
• 10.4.4 正则化技术：L1/L2正则化防止过拟合，提高反演问题的稳定性

10.5 神经网络与物理仿真

• 10.5.1 Physics-Informed Neural Networks：损失函数L=L_data+λL_physics，L_physics=||∂u/∂t+N[u]||²强制满足PDE
• 10.5.2 图神经网络(GNN)：节点表示粒子，边表示相互作用，消息传递更新状态h_i^{t+1}=σ(W·aggregate(h_j^t))
• 10.5.3 神经常微分方程：dy/dt=f_θ(y,t)，使用伴随方法反向传播，连续时间动力学建模
• 10.5.4 混合方法：神经网络学习残差/修正项，物理求解器提供基础解，结合两者优势

10.6 强化学习应用

• 10.6.1 环境建模：使用Taichi构建可微物理环境，状态空间=粒子位置速度，动作空间=外力/控制参数
• 10.6.2 奖励设计：任务奖励+正则化项，避免reward hacking，curriculum learning逐步增加难度
• 10.6.3 策略梯度方法：PPO/TRPO稳定训练，利用物理梯度∂r/∂a加速收敛，actor-critic架构
• 10.6.4 仿真到现实迁移：域随机化(domain randomization)，系统辨识校准参数，添加噪声提高鲁棒性

10.7 案例研究

• 10.7.1 软体机器人控制：优化驱动器位置和控制序列，目标函数=到达目标位置，使用CMA-ES或梯度方法
• 10.7.2 流体形状优化：优化障碍物形状减小阻力，伴随方法计算形状导数，level set表示形状
• 10.7.3 材料设计：拓扑优化最大化刚度/重量比，SIMP方法ρ^p插值，敏感度滤波避免棋盘格
• 10.7.4 动画生成：优化初始速度场生成目标形状，时间反向积分，艺术风格控制

10.8 未来展望

• 10.8.1 可微渲染：结合物理仿真和渲染，端到端优化从参数到像素，逆向图形学应用
• 10.8.2 量子计算应用：量子算法加速线性系统求解(HHL算法)，量子退火优化问题
• 10.8.3 神经隐式表示：NeRF/DeepSDF表示几何和物理场，连续可微表示，内存效率高
• 10.8.4 大规模优化：分布式自动微分，混合精度训练，稀疏梯度处理，千万自由度优化

附录

A. 数学基础回顾

• 向量与矩阵运算
• 张量记号与爱因斯坦求和
• 微分几何基础
• 变分法与泛函分析

B. Taichi API参考

• 常用函数速查
• 数据结构模板
• 性能调优参数
• 调试工具集

C. 代码模板库

• 基础物理系统实现
• 求解器模板
• 可视化工具
• 测试框架

D. 扩展阅读

• 经典教材推荐
• 重要论文列表
• 开源项目资源
• 在线课程链接

术语表

[按照CLAUDE.md中的中英对照表整理]

参考文献

[根据各章内容整理相关文献]