PCISPH2

一句话，PCISPH的核心思想是通过迭代压力梯度力最小化密度误差。

1. 核心思想的推导

这里有两个核心问题
第一个问题是：为什么要最小化密度误差？
第二个问题是：怎么把压力梯度力与密度联系起来？

1.1 为什么要最小化密度误差

第一个问题是：为什么要最小化密度误差？

我们知道，不可压缩性的核心在于密度不能改变，或者说体积不能改变。

那么就有控制密度误差和控制体积误差两种思路。

在连续性方程中，分别对应着左边的密度项与右边的速度散度项（速度的散度就是体积变化量，这点可以看本章附录2）。所以两种思路在实质上是等价的。其实质都是希望连续性方程为0，即

PCISPH选择的思路是最小化密度误差，即左边项。

1.2 压力与密度建立联系

第二个问题是：怎么把压力梯度力与密度联系起来？

压力要与密度建立联系，其核心就是使用压力泊松方程。关于压力泊松方程的推导，见本章附录。如不了解，请务必先看附录。

对压力泊松方程进行一系列的改造，我们可以建立，密度变化量与压力变化量的线性关系（不是那么严格的）。虽然这个关系式并不严格，但是却只需要一个系数即可：

左边的表示压力的变化量
右边的表示密度的变化量
delta就是那个系数

市面上绝大部分PCISPH的推导都围绕着如何推导出delta的公式，我们这里就不再过多赘述。

注意，~的意思是，这个量并不是那么严格的，而是加了某种简化的假设的。
这个假设就是假设迈出后半步（校正步）后，粒子所有邻域的压力和密度都是相等的，并且密度都等于静止密度。这显然是个理想状态，不可能发生。但是却能得到极大简化，让delta在当前时间步为一定值。之所以能进行如此假设，是因为我们并不清楚迈出后半步后的邻域是怎么样的。如果要严格来预测，必须对全局所有粒子同时进行，这就必须同时解所有粒子的校正步。而这种方法需要线性求解器，IISPH就是这种方法的其中一种。所以PCISPH的好处就在于，不需要线性求解器。

这个*的意思是预测量，也就是前半步所得到的值。
这个err的意思是代表与静止密度之间的误差。

附录1： 通过连续性方程推导压力泊松方程

已知连续性方程为：

这个公式的核心作用在于：把速度和密度相互转换

预测：只考虑非压力梯度力

先移动半步，即只考虑非压力梯度力得到预测速度位置与预测密度。

用星号表示预测量。

预测速度

把预测速度带入到连续性方程得到预测密度

校正：只考虑压力梯度力

再移动半步，即只考虑压力梯度力。
由压力梯度力除以质量得到压力加速度，即

再对压力加速度时间积分得到压力对应的速度变化量，即

这个速度变化量只考虑了压力的作用。

再把这个速度变化量带入到连续性方程，得到对应的密度变化量，
这个密度变化量只考虑了压力的作用。

抵消：校正量恰好抵消预测量

静止密度是真实密度。我们希望移动一步之后，密度恰好等于静止密度。

于是就让校正的密度变化量与预测的密度变化量恰好加起来为0。

此时第一项的rho0就是假设了后半步恰好走到静止密度。

这个式子也可以看成就是连续性方程，第一项就是迈步的目标，右边就是实际迈的步长。

核心思想在于，原本连续性方程是密度与速度相互转换，现在改造为了**密度与压力的相互转换**。

改造的核心技巧就在于对压力梯度力求时间积分，把加速度变为了速度。

整理

把第二项分子分母的rho(t)消掉，两项移动到等号两边，就得到了。

这是一个只关于压力的方程，因为右边的rho*是已知量，在前半步，也就是只考虑非压力梯度力的时候，就已经得到了。

这类左边是拉普拉斯项，右边是定值的二阶微分方程，在数学上称之为泊松方程。又由于是关于压力的，所以称之为压力泊松方程。

至此，我们推导出了压力泊松方程。

总结一下：压力泊松方程是只关于后半步的，也就是校正步的。不断迭代校正步，找到正确的步长，就是压力泊松方程的目的。

附录2：为什么速度的散度就是体积变化量

我们看到，速度可以分解为Vx Vy Vz三个方向，对于微元体来说，就是如上图。

而散度为：
$$\nabla \cdot \vec{V}=\frac{\partial V_x}{\partial x}+\frac{\partial V_y}{\partial y}+\frac{\partial V_z}{\partial z}$$

借助高斯公式，对体积的积分可以变为对边界面的积分

$\Omega$就是这个微元体的体积，
$\partial \Omega$就是其表面

可见，
dydz就是前后两个面的面积
dxdz就是左右两个面的面积
dydx就是上下两个面的面积

那么，右边的公式其实就是相应的速度乘以速度正对着的面积
恰好就是整个微元体的体积增量。

所以，速度的散度（对体积积分后）就是体积的变化量。

2. 算法流程

这里只大致阐述：

前置步骤

1. 1-3的邻域搜索不必阐述，是SPH常规流程。
2. 4-5在计算非压力梯度力，这也是SPH常规流程。
3. 6-7 始化压力和压力梯度力为0。

核心循环

4. 8-17就是核心的循环，跳出条件为压力误差小于给定值。下面详述：

其中9-11在利用非压力梯度力计算校正速度，这个v* 和x*和我们前面在
“附录1： 通过连续性方程推导压力泊松方程”中有所不同，是校正后的速度位置，即考虑了压力梯度力之后的。

其中12-15是利用v* 和x* 计算rho*，其公式为

这是连续性方程改造出来的公式。

然后用rho* 计算压力增量p~，其公式为

p~ 累加到原p上，就是校正后的p

16-17 就是利用校正后的p计算压力梯度力。这点与常规SPH无异。

后续步骤

4. 18-20 是计算新的速度和位置，这就是常规SPH中的时间积分。

结束

所以，PCISPH与常规SPH的核心区别就在于多了中间那个循环。

当然也多了循环的初始化。总结起来就是6-18步。

参考资料

B. Solenthaler and R. Pajarola. 2009. Predictive-corrective incompressible SPH. ACM Trans. Graph. 28, 3, Article 40 (August 2009), 6 pages. DOI:https://doi.org/10.1145/1531326.1531346