# Popular explanation and simple implementation of PCISPH

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

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()

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()
dV = velocity[i] - velocity[j]
dA = acceleration[i] - acceleration[j] #a_nonp
sum1 += kernelFunc(r)

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)

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()

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()
dA = accP[i] - accP[j]

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()
dA = accP[i] - accP[j]

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:

Posted on Tue, 23 Nov 2021 10:27:19 -0500 by shinichi_nguyen