PCISPH的通俗解释与简单实现

PCISPH是SPH的一个变种，针对压力求解部分进行了改进。

本文参照

1. Koschier, D., Bender, J., Solenthaler, B., & Teschner, M.
   (2020). Smoothed particle hydrodynamics techniques for the
   physics based simulation of fluids and solids. arXiv preprint
   arXiv:2009.06944.（主要参考）主要是P14-15
   https://interactivecomputergraphics.github.io/SPH-Tutorial/

伪代码

取自参考[1]

算法流程通俗解释

我们只关注PCISPH相较于WCSPH的改进部分，就是计算压力和压力梯度力的部分。其余部分完全一致。

WCSPH中直接通过状态方程解算出压力，而这里通过循环迭代求解。该循环就是上述伪代码中的算法。

1. 我们将非压力梯度力的部分，即粘性力和重力先求出来，加入到加速度中，并通过时间积分预测出一个位置速度（此时，未施加压力梯度力，相当于迭代初始时压力梯度力为0）。

2. 根据预测出的位置，使用核函数公式得到预测密度density_star。再减去静止密度得到密度差。

3. 通过密度差乘以一个系数k_PCI得到预测压力。

4. 根据预测压力，由核函数离散公式计算出预测压力梯度力

5. 根据预测的压力梯度力，加入到加速度中，进行时间积分后得到新的预测位置和速度，进而得到新的预测密度（此时称为校正密度）。这一步虽然是多步，但是可以写为一个公式。

6. 同样是密度差乘以k_PCI，得到校正压力

7. 假如密度差小于规定值，就跳出循环，否则，重复第4、5、6步，不断迭代压力梯度力以最小化密度差。

这里我们看出来，循环的自变量是压力梯度力，循环的跳出条件则是密度误差。通过不断更正压力梯度力，从而最小化密度误差。

公式和公式翻译成代码

主要参考[1]的公式（53）-（60）

所有的代码都是针对粒子i的，最外层还需要嵌套对粒子i的循环。

1. 计算k_PCI
   $$k_{PCI}=-\frac{0.5{\rho }_0^2}{\Delta t^2{m}^2}\cdot \frac{1}{{\Sigma }_j\nabla W_{ij}\cdot {\Sigma }_j\nabla W_{ij}+{\Sigma }_j(\nabla W_{ij}\cdot \nabla W_{ij})}$$

def formula_kPCI(i):
    sum1=ti.Vector([0.0, 0.0])
    sum2=ti.Vector([0.0, 0.0])
    for k in range(numNeighbor[i]):
        j = neighbor[i, k]      
        dir = (position[i]-position[j]).normalized()
        r = (position[i]-position[j]).norm()
        gradW= firstDW(r) * dir
        sum1 += gradW
        sum2 += gradW.dot(gradW)

    kPCI= - 0.5 * density0 **2 / timeStepSize **2 / mass**2 \
            * 1.0/ ( sum1.dot( sum1) + sum2 )

2. 计算density_star
   

def formula_densityStar(i):
    sum1 = ti.Vector([0.0, 0.0])
    sum2 = ti.Vector([0.0, 0.0])
    sum3 = ti.Vector([0.0, 0.0])
    for k in range(numNeighbor[i]):
        j = neighbor[i, k] 
        dir = (position[i]-position[j]).normalized()
        r = (position[i]-position[j]).norm()
        gradW = firstDW(r) * dir
        dV = velocity[i] - velocity[j]
        dA = acceleration[i] - acceleration[j] #a_nonp
        sum1 += kernelFunc(r)
        sum2 += dV.dot(gradW)
        sum3 += dA.dot(gradW)

    #densityStar
    densityStar=  sum1 * mass \
                + sum2 * mass * timeStepSize  \
                + sum3 * mass * timeStepSize * timeStepSize

3. 计算预测压力
   
   我们这里称之为pressureStar

def formula_pressureStar():
    pressureStar = kPCI* (densityStar-density0)

4. 计算压力梯度加速度
   

def formula_accP(i):
    sum1 = ti.Vector([0.0, 0.0])

    for k in range(numNeighbor[i]):
        j = neighbor[i, k] 
        dir = (position[i]-position[j]).normalized()
        r = (position[i]-position[j]).norm()
        gradW = firstDW(r) * dir
        sum1 += gradW

    accP= -1.0 * mass * 2 * pressureStar / (density0 *density0) \
                    * sum1

5. 计算校正密度

与densityStar的区别仅仅在于把accNonP换为了accP(accNonP用acceleration表示了）

def formula_densityP(i):
    sum1 = ti.Vector([0.0, 0.0])
    for k in range(numNeighbor[i]):
        j = neighbor[i, k] 
        dir = (position[i]-position[j]).normalized()
        r = (position[i]-position[j]).norm()
        gradW = firstDW(r) * dir
        dA = accP[i] - accP[j] 
        sum1 += dA.dot(gradW)
        
    densityP = sum1 * mass * timeStepSize * timeStepSize

6. 计算矫正压力
   
   由于要迭代，所以
   我们把该压力（等式左边的）称之为pressureNew
   右边的称之为pressureOld
   最开始，pressureOld=pressureStar
   后来不断迭代pressureOld与pressureNew

def formula_pressureNew(i):
    sum1 = ti.Vector([0.0, 0.0])
    for k in range(numNeighbor[i]):
        j = neighbor[i, k] 
        dir = (position[i]-position[j]).normalized()
        r = (position[i]-position[j]).norm()
        gradW = firstDW(r) * dir
        dA = accP[i] - accP[j] 
        sum1 += dA.dot(gradW)
        
    pressureNew =   pressureOld \
                    + kPCI * \
                    (
                        density0 - densityStar 
                    -   timeStepSize * timeStepSize * mass * sum1
                    )
    pressureOld = pressureNew

7 计算收敛条件

def formula_isConverge(i):
    judge = ( density0 - densityStar + densityP ) / density0
    if judge <1e-3:
        isConverge = True

代码的重整与合并

我们发现许多代码部分是重合的，所以可以进行合并，合并后的代码为：