PCISPH is a variant of SPH, which is improved for the pressure solution part.

This article refers to

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. (main reference) mainly P14-15
   https://interactivecomputergraphics.github.io/SPH-Tutorial/

Pseudo code

From reference [1]

Popular explanation of algorithm flow

We only focus on the improvement of PCISPH compared with WCSPH, that is, the calculation of pressure and pressure gradient force. The rest is exactly the same.

In WCSPH, the pressure is calculated directly by solving the equation of state, and here it is solved by cyclic iteration. The loop is the algorithm in the above pseudo code.

1. We first calculate the non pressure gradient force, i.e. viscous force and gravity, add it to the acceleration, and predict a position velocity through time integration (at this time, no pressure gradient force is applied, which is equivalent to that the pressure gradient force is 0 at the beginning of the iteration).

2. According to the predicted position, the predicted density is obtained by using the kernel function formula_ star. Then subtract the static density to get the density difference.

3. Multiply the density difference by a factor k_PCI gets the predicted pressure.

4. According to the predicted pressure, the predicted pressure gradient force is calculated by the discrete formula of kernel function

5. According to the predicted pressure gradient force, it is added to the acceleration. After time integration, a new predicted position and velocity are obtained, and then a new predicted density (called corrected density at this time). Although this step is multi-step, it can be written as a formula.

6. It's also the density difference times k_PCI, get the corrected pressure

7. If the density difference is less than the specified value, jump out of the cycle. Otherwise, repeat steps 4, 5 and 6 and continuously iterate the pressure gradient force to minimize the density difference.

Here we can see that the independent variable of the cycle is the pressure gradient force, and the jump out condition of the cycle is the density error. The density error is minimized by continuously correcting the pressure gradient force.

Formula and formula translation into code

Mainly refer to formula (53) - (60) of [1]

All codes are for particle i, and the outermost layer also needs to nest the loop to particle i.

1. Calculate k_PCI
   k P C I = − 0.5 ρ 0 2 Δ t 2 m 2 ⋅ 1 Σ j ∇ W i j ⋅ Σ j ∇ W i j + Σ j ( ∇ W i j ⋅ ∇ W i j ) 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})} kPCI​=−Δt2m20.5ρ02​​⋅Σj​∇Wij​⋅Σj​∇Wij​+Σj​(∇Wij​⋅∇Wij​)1​

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 * gradW

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

1. Calculate 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)

    sum1 *= mass
    sum2 *= mass * timeStepSize
    sum3 *= mass * timeStepSize * timeStepSize

    densityStar=   sum1 + sum2 + sum3

1. Calculate predicted pressure
   
   We call it pressureStar here

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

1. Calculate pressure gradient acceleration
   

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

1. Calculate correction density

The only difference from densityStar is that accNonP is replaced by accP(accNonP is represented by 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)
        
    sum1 *= mass * timeStepSize * timeStepSize

    densityStar=   sum1

1. Calculate correction pressure
   
   Because of iteration, so
   We call this pressure (on the left side of the equation) pressureNew
   The one on the right is called pressureOld
   Initially, pressureOld=pressureStar
   Later, we iterated over pressureOld and 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 calculation convergence conditions

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

Code reorganization and merger

We found that many code parts coincide, so we can merge them. The merged code is: