Semi Lagrangian Advection
“全拉格朗日”是：就去算一个个粒子，网格啥的一概没有。
“半拉格朗日”就“半”在，它每次用粒子方式Advect之后，就又回到了欧拉网格上。
具体地，第n+1步的A(i,j)等于第n步的A(x,y)，其中(x,y)是(i,j)这个点在第n步时应该在的位置。
MacCormack Scheme
考虑这个方程
[latex]\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0[/latex]
MacCormack方法分成两步：Predictor Step和Corrector Step.
Predictor Step：
把方程离散化，就变成了
[latex]\frac{u_i^{\overline{n+1}}-u_i^n}{\Delta t}+a\frac{u_{i+1}^n-u_i^n}{\Delta x}=0[/latex]
这个[latex]u_i^{\overline{n+1}}[/latex]表示的是“预估的第n+1步的值”的意思。把它写成explicit形式：
[latex]u_i^{\overline{n+1}}=u_i^n-a\frac{\Delta t}{\Delta x}\left(u_{i+1}^n-u_i^n\right)[/latex]
Corrector Step：
假设我们知道时间n+0.5的信息，那么就应该有
[latex]\frac{u_i^{n+1}-u_i^{n+\frac{1}{2}}}{1/2\Delta t}+a\frac{u_i^{\overline{n+1}}-u_{i-1}^{\overline{n+1}}}{\Delta x}=0[/latex]
注意，Predictor Step里面，关于x偏导那一项是i+1减i（前向差分），但Corrector Step里面就是后向差分了。这两个也可以交换，先后向差分再前向差分。但是必须不一样。这是达到二阶精度所必须做的。
Corrector Step关于x偏导是用“预估值”[latex]u^{\overline{n+1}}[/latex]做的，需要注意。
我们不知道[latex]u_i^{n+\frac{1}{2}}[/latex]是什么，就用“预估值”和[latex]u^n[/latex]进行插值：
[latex]u_i^{n+\frac{1}{2}}=\frac{1}{2}\left( u_i^n+u_i^{\overline{n+1}} \right)[/latex]
代入进去就变成了
[latex]u_i^{n+1}=\frac{u_i^n+u_i^{\overline{n+1}}}{2}-a\frac{\Delta t}{2\Delta x}\left( u_i^{\overline{n+1}}-u_{i-1}^{\overline{n+1}}\right)[/latex]
BFECC Scheme
BFECC是缩写，全拼是Back and Forth Error Compensation and Correction.
假设有场[latex]\phi[/latex].
先用“正向”Advect的算子[latex]A[/latex]算出一个预估值[latex]\hat{\phi}^{n+1}=A(\phi^n)[/latex]
然后再用“反向”的算子[latex]A^R[/latex]算出第n步的“估计”值[latex]\hat{\phi}^n=A^R(\hat{\phi}^{n+1})[/latex]
然后我们估计得到一个误差[latex]e=(\hat{\phi}^n-\phi^n)/2.[/latex]
然后用这个误差修正初始值：[latex]\bar{\phi}^n=\phi^n-e[/latex]
并重新Advect：[latex]\phi^{n+1}=A(\bar{\phi}^n)[/latex]
Projection
Advect并不保不可压性。假设Advect得到流场[latex]\textbf{v}[/latex]。
根据亥姆霍兹分解，我们可以设想，求一个“压力”场[latex]q[/latex]，使得[latex]\nabla^2q=\nabla\cdot\textbf{v}[/latex]
然后更新[latex]\textbf{v}_1=\textbf{v}-\nabla q[/latex]
我们设[latex]\textbf{v}=(u,v)[/latex].右边对[latex]\textbf{v}[/latex]求散度，由于拉普拉斯算子对称，所以这里也可以用中心差分：
[latex](\nabla\cdot\textbf{v})_{i,j}=(u_{i+1,j}-u_{i,j})/2\Delta x+(v_{i,j+1}-v_{i,j-1})/2\Delta y[/latex]
前面的拉普拉斯算子写出来就是
[latex](\nabla^2q)_{i,j}=(q_{i+1,j}-2q_{i,j}+q_{i-1,j})/\Delta x^2+(q_{i,j+1}-2q_{i,j}+q_{i,j-1})/\Delta y^2[/latex]
我们要解的是[latex]q[/latex].如果把它写成一个线性方程组的形式，那么对角线就是[latex]-2(1/\Delta x^2+1/\Delta y^2)[/latex]，(i,j)那一行里面，(i+1,j)和(i-1,j)对应项的元素是[latex]1/\Delta x^2[/latex]，(i,j-1)和(i,j+1)对应项都是[latex]1/\Delta y^2[/latex].方程右边的常数向量则是展开的[latex](\nabla\cdot\textbf{v})[/latex].
这里还要注意一个Boundary Condition的情况。假设我们的网格范围是(1,1)~(n-1,m-1)，外面一圈是“假”的边界值。如果采取“复制型”的边界条件，就是例如q(0,j)=q(1,j)这样。这时候(1,j)这一行就少了一项。