PCISPH is a variant of SPH, which is improved for the pressure solution part.
This article refers to
 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 P1415
https://interactivecomputergraphics.github.io/SPHTutorial/
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.

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

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.

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

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

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 multistep, it can be written as a formula.

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

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.
 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 )
 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
 Calculate predicted pressure
We call it pressureStar here
def formula_pressureStar(): pressureStar = kPCI* (densityStardensity0)
 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
 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
 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 <1e3: isConverge = True
Code reorganization and merger
We found that many code parts coincide, so we can merge them. The merged code is: