学习教程来自:GAMES201：高级物理引擎实战指南2020
以下大部分图片来自教程PPT，仅作为笔记用于学习和分享，侵删


笔记内容大多为课程内容的翻译和转述，外加一些自己的理解，若有不正确的地方恳请大家交流和指正

笔记

1. 欧拉视角的计算方法概述

1.1 材料导数 Material Derivatives

物理量的变化 = 时间上的变化（欧拉视角下？） + 由于移动的变化


</br>

1.2 N-S方程 Navier–Stokes equations


</br>

速度的变化 = 压强 + 粘度（一般遗弃） + 重力

由于不可压缩性 速度的散度为0


</br>

1.3 求解

使用Operator splitting方法拆分

  1. 不考虑压强的重力的作用，先计算速度场的变化（欧拉空间下，每个节点固定不动，故每个time step后要重新计算当前节点的速度）


</br>

  1. 施加重力加速度（或其他外力）的作用，加速度为g


</br>

  1. 施加压强产生的加速度，并使得累加后，速度散度为0


</br>

2. Grid 网格方法计算

2.1 均匀网格下的存储

  1. cell-centered grids：都存在网格中心
  2. staggered grids：在不同的位置存x方向速度、y方向速度、压强，有利于计算有限元差分


</br>

2.2 网格中的插值

双线性插值 Bilinear interpolation：每个顶点的值x如图对应颜色的面积 累加


</br>

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），速度场中插值求得速度


</br>

沿直线回溯Forward Euler (“RK1”)：

p -= dt * velocity(p) # 直接回溯

但假设会导致误差，因为速度场实际是会变化的，所以沿直线回溯不一定能回到之前的位置


</br>

显式中点法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) # 加权平均

以上方法会在一段时间的运算后慢慢模糊，作为速度场可视作能量的衰减，流体变粘


</br>

下面这个方法会改善这种情况

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


</br> 

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


</br>

4.1 求解过程（2维）


</br>

拉普拉斯算子（Laplace operator）


</br>

中心差分法近似拉普拉斯（这里使用了2维的five point stencil？）先求梯度再求散度


</br>

近似u的散度（流入和流出）


</br>


</br>

整个线性系统（线性方程组），其中A是稀疏矩阵


</br>

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）的
存储方法：

  1. dense matrix：直接存，规模不大时
  2. sparse matrix：用稀疏矩阵的方式存，CSR、COO等
  3. 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足够小


</br>

4.2.4 Preconditioning 使解更快

condition number κ：一个评价收敛速度的值，越小收敛越快，等于SPD的最大特征值（max eigenvalue）/最小特征值


</br>

如何使得condition number变小：找到一个近似的矩阵M与A相近，且容易求逆（对角阵）。左右同乘M逆，此时condition number可能会变小


</br>

4.2.5 Multigrid methods

Multigrid preconditioned conjugate gradients (MGPCG)求解Poisson’s equation

5. 一些改进的Paper

  1. Restoring the missing vorticity inadvection-projection fluid solvers：修正速度场在经过advection和projection之后能量的损失，降低了流体的数值粘性
  2. An advection-reflection solver for detail-preserving fluid simulation：另一种降低能量损失的方法

5.1 一些扩展方向

  1. 从2D到3D模拟
  2. 精确的边界以及流体固体的耦合
  3. Two phase fluid simulation
  4. 处理自由表面，level sets方法
  5. 处理涡量守恒
赏