Transfer between grid and particles

Transfer macros describe how quantities move between particles and grid nodes. @P2G evaluates particle-to-grid operations, such as scattering mass, momentum, or internal force to the grid. @G2P evaluates grid-to-particle operations, such as updating particle velocity, position, or deformation state from grid values.

Inside a transfer macro, grid fields, particle fields, and basis weights are written with indexed notation. This notation keeps the code close to the transfer equations while letting Tesserae choose the execution backend, including sequential CPU execution, threaded CPU execution, and GPU kernels.

Transfer macros

@P2G grid=>i particles=>p weights=>ip [partition] begin
    equations...
end

@P2G grid=>i particles=>p weights=>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)
    bw = weights[p]
    nodeindices = supportnodes(bw)
    for ip in eachindex(nodeindices)
        i = nodeindices[ip]
        grid.m [i] += bw.w[ip] * particles.m[p]
        grid.mv[i] += bw.w[ip] * particles.m[p] * particles.v[p]
        grid.mv[i] += -particles.V[p] * particles.σ[p] * bw.∇w[ip]
    end
end

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

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

@G2P grid=>i particles=>p weights=>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)
    bw = weights[p]
    nodeindices = supportnodes(bw)
    Δvₚ = zero(eltype(particles.v))
    ∇vₚ = zero(eltype(particles.∇v))
    Δxₚ = zero(eltype(particles.x))
    for ip in eachindex(nodeindices)
        i = nodeindices[ip]
        Δvₚ += bw.w[ip] * (grid.vⁿ[i] - grid.v[i])
        ∇vₚ += grid.v[i] ⊗ bw.∇w[ip]
        Δxₚ += bw.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

@G2P2G grid=>i particles=>p weights=>ip [partition] begin
    equations...
end

@G2P2G grid=>i particles=>p weights=>ip begin
    # G2P
    ∇v[p] = @∑ v[i] ⊗ ∇w[ip]

    # Particle update
    F[p] = (I + ∇v[p]*Δt) * F[p]
    σ[p] = cauchy_stress(F[p])

    # P2G
    f[i] = @∑ -V[p] * σ[p] * ∇w[ip]
end

@foreach collection=>i begin
    statements...
end

@foreach collection[:,begin]=>i begin
    statements...
end

Backend-aware local loops

Use @foreach for local grid or particle updates that are not transfers. It follows the same collection=>i field-access convention as transfer macros:

@foreach grid=>i begin
    m⁻¹[i] = ifelse(iszero(m[i]), zero(m[i]), inv(m[i]))
    v[i] = mv[i] * m⁻¹[i]
end

@foreach particles=>p begin
    x[p] += v[p] * Δt
end

Index the collection argument to visit a slice. This is useful for boundary conditions and uses the same CPU, GPU, sparse-grid, and threading dispatch as the full loop:

@foreach grid[:,:,begin]=>i begin
    v[i] = v[i] .* (true, true, false)
end

On dense grids, @foreach visits all grid nodes. On SpGrid, it visits only active sparse nodes. When the collection is on GPU, the loop is dispatched as a GPU kernel. Prefix it with @threaded to parallelize CPU loops.

Inspecting transfer code

@explain @P2G ...
@explain @G2P ...
@explain @G2P2G ...
@explain @P2G_Matrix ...

ExplainedCode

Transfer macros hide the explicit loops over particles, support nodes, and grid nodes. To inspect that structure, prefix a transfer macro with @explain:

@explain @P2G grid=>i particles=>p weights=>ip begin
    m[i] = @∑ w[ip] * m[p]
    mv[i] = @∑ w[ip] * m[p] * v[p]
end

This prints readable CPU reference code for understanding and debugging. It is not the optimized lowering used by the transfer macro.

The same form works for @G2P, @G2P2G, and @P2G_Matrix. For threaded transfers, place @explain before @threaded:

@explain @threaded @P2G grid=>i particles=>p weights=>ip partition begin
    m[i] = @∑ w[ip] * m[p]
end

For scattering transfers such as @P2G, @G2P2G, and @P2G_Matrix, the threaded reference code follows the same ThreadPartition structure as the actual CPU transfer: it loops over threadsafe_groups(partition) and then over the particle indices assigned to each thread-safe region. Without a partition, the reference code shows the sequential fallback.

Code snippets

This section lists common velocity transfer schemes:

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

As a rough guide, PIC–FLIP is the minimal baseline, APIC improves angular momentum behavior by carrying an affine velocity field, TPIC carries a velocity gradient directly, and XPIC reduces transfer noise through a recursive correction.

The snippets below assume the following grid and particle fields:

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 weights=>ip begin
    mv[i] = @∑ w[ip] * m[p] * v[p]
end
@G2P grid=>i particles=>p weights=>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 weights=>ip begin
    mv[i] = @∑ w[ip] * m[p] * (v[p] + B[p] * Dₚ⁻¹ * (x[i] - x[p]))
end
@G2P grid=>i particles=>p weights=>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 weights=>ip begin
    mv[i] = @∑ w[ip] * m[p] * (v[p] + ∇v[p] * (x[i] - x[p]))
end
@G2P grid=>i particles=>p weights=>ip begin
    v[p]  = @∑ w[ip] * v[i]
    ∇v[p] = @∑ v[i] ⊗ ∇w[ip]
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 used to update the particle velocity $\bm{V}$ and position $\bm{X}$ 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

The equations above can be written as

\[\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 weights=>ip begin
        vᵣ★[p] = @∑ w[ip] * vᵣ★[i]
    end
    @P2G grid=>i particles=>p weights=>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 weights=>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}^*$.