Transfer between grid and particles

Transfer macros

@P2G grid=>i particles=>p mpvalues=>ip [space] begin
    equations...
end

@P2G grid=>i particles=>p mpvalues=>ip begin

    # Particle-to-grid transfer
    m[i]  = @∑ w[ip] * m[p]
    mv[i] = @∑ w[ip] * m[p] * v[p]
    f[i]  = @∑ -V[p] * σ[p] ⋅ ∇w[ip]

    # Calculation on grid
    vⁿ[i] = mv[i] / m[i]
    v[i]  = vⁿ[i] + (f[i] / m[i]) * Δt

end

# Reset grid properties
@. grid.m  = zero(grid.m)
@. grid.mv = zero(grid.mv)
@. grid.f  = zero(grid.f)

# Particle-to-grid transfer
for p in eachindex(particles)
    mp = mpvalues[p]
    nodeindices = neighboringnodes(mp)
    for ip in eachindex(nodeindices)
        i = nodeindices[ip]
        grid.m [i] += mp.w[ip] * particles.m[p]
        grid.mv[i] += mp.w[ip] * particles.m[p] * particles.v[p]
        grid.mv[i] += -particles.V[p] * particles.σ[p] ⋅ mp.∇w[ip]
    end
end

# Calculation on grid
@. grid.vⁿ = grid.mv / grid.m
@. grid.v  = grid.vⁿ + (grid.f / grid.m) * Δt

@G2P grid=>i particles=>p mpvalues=>ip begin
    equations...
end

@G2P grid=>i particles=>p mpvalues=>ip begin

    # Grid-to-particle transfer
    v[p] += @∑ w[ip] * (vⁿ[i] - v[i])
    ∇v[p] = @∑ v[i] ⊗ ∇w[ip]
    x[p] += @∑ w[ip] * v[i] * Δt

    # Calculation on particle
    Δϵₚ = symmetric(∇v[p]) * Δt
    F[p]  = (I + ∇v[p]*Δt) ⋅ F[p]
    V[p]  = V⁰[p] * det(F[p])
    σ[p] += λ*tr(Δϵₚ)*I + 2μ*Δϵₚ # Linear elastic material

end

# Grid-to-particle transfer
for p in eachindex(particles)
    mp = mpvalues[p]
    nodeindices = neighboringnodes(mp)
    Δvₚ = zero(eltype(particles.v))
    ∇vₚ = zero(eltype(particles.∇v))
    Δxₚ = zero(eltype(particles.x))
    for ip in eachindex(nodeindices)
        i = nodeindices[ip]
        Δvₚ += mp.w[ip] * (grid.vⁿ[i] - grid.v[i])
        ∇vₚ += grid.v[i] ⊗ mp.∇w[ip]
        Δxₚ += mp.w[ip] * grid.v[i] * Δt
    end
    particles.v[p] += Δvₚ
    particles.∇v[p] = ∇vₚ
    particles.x[p] += Δxₚ
end

# Calculation on particle
for p in eachindex(particles)
    Δϵₚ = symmetric(particles.∇v[p]) * Δt
    particles.F[p]  = (I + particles.∇v[p]*Δt) ⋅ particles.F[p]
    particles.V[p]  = particles.V⁰[p] * det(particles.F[p])
    particles.σ[p] += λ*tr(Δϵₚ)*I + 2μ*Δϵₚ # Linear elastic material
end

Code snippets

In this section, the following transfer schemes are presented:

• PIC–FLIP mixed transfer^[1]
• Affine PIC (APIC)^[2]
• Taylor PIC (TPIC)^[3]
• eXtended PIC (XPIC)^[4]

We assume that the grid and particle properties have following variables:

GridProp = @NamedTuple begin
    x   :: Vec{2, Float64} # Position
    m   :: Float64         # Mass
    mv  :: Vec{2, Float64} # Momentum
    v   :: Vec{2, Float64} # Velocity
    vⁿ  :: Vec{2, Float64} # Velocity at t = tⁿ
    # XPIC
    vᵣ★ :: Vec{2, Float64}
    v★  :: Vec{2, Float64}
end
ParticleProp = @NamedTuple begin
    x  :: Vec{2, Float64}                           # Position
    m  :: Float64                                   # Mass
    V⁰ :: Float64                                   # Initial volume
    V  :: Float64                                   # Volume
    v  :: Vec{2, Float64}                           # Velocity
    ∇v :: SecondOrderTensor{2, Float64, 4}          # Velocity gradient
    σ  :: SymmetricSecondOrderTensor{2, Float64, 3} # Cauchy stress
    F  :: SecondOrderTensor{2, Float64, 4}          # Deformation gradient
    # APIC
    B  :: SecondOrderTensor{2, Float64, 4}
    # XPIC
    vᵣ★ :: Vec{2, Float64}
    a★  :: Vec{2, Float64}
end

PIC–FLIP mixed transfer

\[\begin{aligned} m^n\bm{v}_i^n &= \sum_p w_{ip}^n m_p \bm{v}_p^n \\ \bm{v}_p^{n+1} &= \sum_i w_{ip}^n \left( (1-\alpha)\bm{v}_i^{n+1} + \alpha (\bm{v}_p^n + (\bm{v}_i^{n+1} - \bm{v}_i^n) \right)) \\ \end{aligned}\]

@P2G grid=>i particles=>p mpvalues=>ip begin
    mv[i] = @∑ w[ip] * m[p] * v[p]
end
@G2P grid=>i particles=>p mpvalues=>ip begin
    v[p] = @∑ w[ip] * ((1-α)*v[i] + α*(v[p] + (v[i]-vⁿ[i])))
end

Affine PIC (APIC)

\[\begin{aligned} m^n\bm{v}_i^n &= \sum_p w_{ip}^n m_p \left(\bm{v}_p^n + \bm{B}_p^n (\bm{D}_p^n)^{-1} (\bm{x}_i^n - \bm{x}_p^n) \right) \\ \bm{v}_p^{n+1} &= \sum_i w_{ip}^n \bm{v}_i^{n+1} \\ \bm{B}_p^{n+1} &= \sum_i w_{ip}^n \bm{v}_i^{n+1} \otimes (\bm{x}_i^n - \bm{x}_p^n) \end{aligned}\]

@P2G grid=>i particles=>p mpvalues=>ip begin
    mv[i] = @∑ w[ip] * m[p] * (v[p] + B[p] * Dₚ⁻¹ * (x[i] - x[p]))
end
@G2P grid=>i particles=>p mpvalues=>ip begin
    v[p] = @∑ w[ip] * v[i]
    B[p] = @∑ w[ip] * v[i] ⊗ (x[i] - x[p])
end

where Dₚ should be defined as Dₚ⁻¹ = inv(1/4 * h^2 * I) for BSpline(Quadratic()) (see the tutorial Transfer schemes).

Taylor PIC (TPIC)

\[\begin{aligned} m^n\bm{v}_i^n &= \sum_p w_{ip}^n m_p \left(\bm{v}_p^n + \nabla\bm{v}_p^n (\bm{x}_i^n - \bm{x}_p^n) \right) \\ \bm{v}_p^{n+1} &= \sum_i w_{ip}^n \bm{v}_i^{n+1} \\ \nabla\bm{v}_p^{n+1} &= \sum_i \bm{v}_i^{n+1} \otimes \nabla w_{ip}^n \\ \end{aligned}\]

@P2G grid=>i particles=>p mpvalues=>ip begin
    mv[i] = @∑ w[ip] * m[p] * (v[p] + ∇v[p] * (x[i] - x[p]))
end
@G2P grid=>i particles=>p mpvalues=>ip begin
    v[p]  = @∑ w[ip] * v[i]
    ∇v[p] = @∑ w[ip] * v[i] ⊗ (x[i] - x[p])
end

eXtended PIC (XPIC)

Overview of XPIC

We assume that $\bm{\mathsf{S}}^+$ matrix maps particle velocities to the grid and $\bm{\mathsf{S}}$ matrix maps them back. In XPIC, a new effective acceleration for particles, $\mathbb{A}$, is employed to update the particle velocity $\bm{V}$ and position $\bm{X}$, respectively, as

\[\begin{aligned} \bm{V}^{(k+1)} &= \bm{V}^{(k)} + \mathbb{A}^{(k)} \Delta{t} \\ \bm{X}^{(k+1)} &= \bm{X}^{(k)} + \bm{\mathsf{S}} \bm{v}^{(k+)} \Delta{t} + \left( \frac{1}{2} \mathbb{A}^{(k)} - \bm{\mathsf{S}}\bm{a}^{(k)} \right) (\Delta{t})^2 \end{aligned}\]

where $\bm{v}^{(k+)}$ is the updated grid velocity:

\[\bm{v}^{(k+)} = \bm{v}^{(k)} + \bm{a}^{(k)} \Delta{t}\]

The effective acceleration $\mathbb{A}$ is represented in XPIC as follows:

\[\mathbb{A}^{(k)} \Delta{t} = (1-m) \bm{\mathsf{S}}\bm{a}^{(k)}\Delta{t} - \bm{V}^{(k)} + m\bm{\mathsf{S}}\bm{v}^{(k+)} - m\bm{\mathsf{S}}\bm{v}^{*}\]

where

\[\bm{v}^* = \sum_r^m (-1)^r \bm{v}_r^*\]

This new $\bm{v}_r^*$ term, which is unique to XPIC($m$) with $m>1$, can be evaluated by recursion:

\[\bm{v}_r^* = \frac{m-r+1}{r} \bm{\mathsf{S}}^+ \bm{\mathsf{S}} \bm{v}_{r-1}^*\]

starting with $\bm{v}_1^*=\bm{v}^{(k)}$.

Implementation using Tesserae

Based on the above equations, we can derive the following:

\[\begin{aligned} \bm{V}^{(k+1)} &= \bm{V}^{(k)} + \bm{\mathsf{S}}\left(\bm{v}^{(k+)}-\bm{v}^{(k)}\right) - \bm{A}^* \Delta{t} \\ \bm{X}^{(k+1)} &= \bm{X}^{(k)} + \frac{1}{2} \bm{\mathsf{S}} \left(\bm{v}^{(k+)} + \bm{v}^{(k)}\right) \Delta{t} - \frac{1}{2} \bm{A}^* (\Delta{t})^2 \end{aligned}\]

where

\[\bm{A}^* \Delta{t} = \bm{V}^{(k)} + m \bm{\mathsf{S}} \left( \bm{v}^* - \bm{v}^{(k)} \right)\]

# Set the initial values for the recursion:
@. grid.vᵣ★ = grid.vⁿ
@. grid.v★ = zero(grid.v★)

# The recursion process to calculate `v★`
for r in 2:m
    @G2P grid=>i particles=>p mpvalues=>ip begin
        vᵣ★[p] = @∑ w[ip] * vᵣ★[i]
    end
    @P2G grid=>i particles=>p mpvalues=>ip begin
        vᵣ★[i] = @∑ (m-r+1)/r * w[ip] * m[p] * vᵣ★[p] / m[i]
        v★[i] += (-1)^r * vᵣ★[i]
    end
end

# Grid-to-particle transfer in XPIC
@G2P grid=>i particles=>p mpvalues=>ip begin
    v[p] += @∑ w[ip] * (v[i] - vⁿ[i]) # same as FLIP
    x[p] += @∑ w[ip] * (v[i] + vⁿ[i]) * Δt / 2
    a★[p] = @∑ w[ip] * (v[p] + m*(v★[i] - vⁿ[i])) / Δt
    v[p] -= a★[p] * Δt
    x[p] -= a★[p] * Δt^2 / 2
end

where a★ represents $\bm{A}^*$.