学习教程来自:GAMES201：高级物理引擎实战指南2020 以下大部分图片来自教程PPT，仅作为笔记用于学习和分享，侵删
笔记内容大多为课程内容的翻译和转述，外加一些自己的理解，若有不正确的地方恳请大家交流和指正
笔记
1. 欧拉视角的计算方法概述
1.1 材料导数 Material Derivatives
物理量的变化 = 时间上的变化（欧拉视角下？） + 由于移动的变化
1.2 N-S方程 Navier–Stokes equations
速度的变化 = 压强 + 粘度（一般遗弃） + 重力
由于不可压缩性 速度的散度为0
1.3 求解
使用Operator splitting方法拆分
- 不考虑压强的重力的作用，先计算速度场的变化（欧拉空间下，每个节点固定不动，故每个time step后要重新计算当前节点的速度）
- 施加重力加速度（或其他外力）的作用，加速度为g
- 施加压强产生的加速度，并使得累加后，速度散度为0
2. Grid 网格方法计算
2.1 均匀网格下的存储
- cell-centered grids：都存在网格中心
- staggered grids：在不同的位置存x方向速度、y方向速度、压强，有利于计算有限元差分
2.2 网格中的插值
双线性插值 Bilinear interpolation：每个顶点的值x如图对应颜色的面积 累加
3. Advection 平流/对流阶段
3.1 Advection schemes
不同的方法在数值粘性viscosity、稳定性、性能、复杂程度做了取舍
Semi-Lagrangian advection：稳定，数值粘性高
MacCormack/BFECC：上边的升级版
BiMocq2
Particle advection (PIC/FLIP/APIC/PolyPIC)：在欧拉网格的速度场中，把物理量记录在例子上
3.2 Semi-Lagrangian advection
假设速度场不变（之前提到advection阶段的先不考虑压强的重力的作用，所以速度场没有变化，这也是平移的含义），则当前位置x的速度 = 沿速度场到上一时间步的位置（backtrace），速度场中插值求得速度
沿直线回溯Forward Euler (“RK1”)：
p -= dt * velocity(p) # 直接回溯但假设会导致误差，因为速度场实际是会变化的，所以沿直线回溯不一定能回到之前的位置
显式中点法Explicit Midpoint (“RK2”)：
p_mid = p - 0.5 * dt * velocity(p) # 先回溯一半，得到中点
p -= dt * velocity(p_mid) # 在中点的位置执行回溯RK3：比RK2更加精确，但是差别不大
v1 = velocity(p) # 取得一系列速度
p1 = p - 0.5 * dt * v1
v2 = velocity(p1)
p2 = p - 0.75 * dt * v2
v3 = velocity(p2)
p -= dt * (2 / 9 * v1 + 1 / 3 * v2 + 4 / 9 * v3) # 加权平均以上方法会在一段时间的运算后慢慢模糊，作为速度场可视作能量的衰减，流体变粘
下面这个方法会改善这种情况
3.3 MacCormack/BFECC
BFECC: Back and Forth Error Compensation and Correction
分别使用3.2的方法先后向前或则x*和x**，x error取1/2(x** - x)，x final = x* + x error
其次，需要视情况裁剪掉一些错误值避免overshooting
source_pos = backtrace(I, dt)
min_val = sample_min(x, source_pos)
max_val = sample_max(x, source_pos)
if new_x[I] < min_val or new_x[I] > max_val:
new_x[I] = sample_bilinear(x, source_pos) # 遗弃之前得到的x final，重新赋值4. Projection
求出一个标量场p，使速度场散度为0
4.1 求解过程（2维）
拉普拉斯算子（Laplace operator）
中心差分法近似拉普拉斯（这里使用了2维的five point stencil？）先求梯度再求散度
近似u的散度（流入和流出）
整个线性系统（线性方程组），其中A是稀疏矩阵
4.2 Solving large-scale linear systems 解大规模线性方程组
解Ax=b 将速度场迭代到无散度
4.2.1 一些解法
Direct solvers（直接求解，如PARDISO）：适用于scaler不是特别大的
Iterative solvers（迭代求解）：Gauss-Seidel、(Damped) Jacobi、(Preconditioned) Krylov-subspace solvers (e.g., conjugate gradients)
以上solver可以进行组合使用
4.2.2 A（如何存储）
稀疏的（sparse）、对称正定（SPD symmetric & positive-definite）的
存储方法：
- dense matrix：直接存，规模不大时
- sparse matrix：用稀疏矩阵的方式存，CSR、COO等
- Matrix-free：不存，究极解法
4.2.3 Krylov-subspace solvers
其中一个变形为共轭梯度法 conjugate gradients (CG)
一些其他在图形学中不常用的方法：CR、GMRES、BiCGStab等
一本关于共轭梯度教程的书：An Introduction to the Conjugate Gradient Method Without the Agonizing Pain5 by Jonathan Richard Shewchuk.
conjugate gradients算法过程：迭代更新x直到r足够小
4.2.4 Preconditioning 使解更快
condition number κ：一个评价收敛速度的值，越小收敛越快，等于SPD的最大特征值（max eigenvalue）/最小特征值
如何使得condition number变小：找到一个近似的矩阵M与A相近，且容易求逆（对角阵）。左右同乘M逆，此时condition number可能会变小
4.2.5 Multigrid methods
Multigrid preconditioned conjugate gradients (MGPCG)求解Poisson’s equation
5. 一些改进的Paper
- Restoring the missing vorticity inadvection-projection fluid solvers：修正速度场在经过advection和projection之后能量的损失，降低了流体的数值粘性
- An advection-reflection solver for detail-preserving fluid simulation：另一种降低能量损失的方法
5.1 一些扩展方向
- 从2D到3D模拟
- 精确的边界以及流体固体的耦合
- Two phase fluid simulation
- 处理自由表面，level sets方法
- 处理涡量守恒
