Multi-threading

Multi-threading in Tesserae parallelizes CPU work over particles, grid nodes, and particle-grid transfers. The @threaded macro adds this parallelism while keeping transfer expressions close to their sequential form. It behaves similarly to Julia's built-in Threads.@threads, but is designed to work with particle-grid transfer macros such as @G2P, @P2G, @G2P2G, and @P2G_Matrix.

Usage guidelines

Particle-grid transfers have two directions: gathering and scattering. In a gathering transfer, each particle reads values from nearby grid nodes, so the operation can be threaded directly. In a scattering transfer, particles write contributions to grid nodes. If multiple threads update the same grid node at the same time, this is a data race; see Julia's discussion of data races between threads. Threaded scattering therefore uses a ThreadPartition.

Gathering (@G2P)

To parallelize @G2P, simply prefix it with @threaded.

@threaded @G2P grid=>i particles=>p weights=>ip begin
    # your code here
end

Scattering (@P2G, @G2P2G and @P2G_Matrix)

For scattering operations, prefix @P2G with @threaded and use ThreadPartition to avoid data races on the grid.

partition = ThreadPartition(mesh)
update!(partition, particles.x) # CartesianMesh only
@threaded @P2G grid=>i particles=>p weights=>ip partition begin
    # your code here
end

For FEMesh and IGAMesh, the partition is built from the fixed cell connectivity, so it does not need an update! call. The same partitioning applies to @G2P2G and @P2G_Matrix.

Updating basis weights

To update basis weights, either use the update! function, or simply:

@threaded for p in eachindex(particles)
    update!(weights[p], particles.x[p], mesh)
end

Reordering particles

For @P2G and related scattering operations, using reorder_particles! together with ThreadPartition can significantly improve cache efficiency and thread scaling:

partition = ThreadPartition(mesh)
update!(partition, particles.x)
reorder_particles!(particles, partition)

Reordering ensures that particles within the same grid block are stored contiguously in memory, reducing random memory access during parallel execution. When reordering is checked every step, use an adaptive threshold such as threshold=0.85:

update!(partition, particles.x)
reorder_particles!(particles, partition; threshold=0.85)

For 0 ≤ threshold ≤ 1, larger values reorder more often. Particles are reordered when Tesserae.block_ordered_particle_contiguity is below threshold.

At the endpoints, threshold=0 never reorders and threshold=1 always reorders.

Multi-threading API

@threaded [scheduler] for ...
@threaded [scheduler] @P2G ...
@threaded @P2G ...

# Parallel loop
@threaded for i in 1:100
    println(i)
end

# Grid-to-particle transfer
@threaded @G2P grid=>i particles=>p weights=>ip begin
    v[p] = @∑ w[ip] * v[i]
end

ThreadPartition(::CartesianMesh)
ThreadPartition(::FEMesh)
ThreadPartition(::IGAMesh)

# Construct ThreadPartition
partition = ThreadPartition(mesh)

# Update partition using current particle positions
update!(partition, particles.x) # Required only for `CartesianMesh`.

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

reorder_particles!(particles, partition; threshold=1)

Tesserae.block_ordered_particle_contiguity(partition)