From 010fecb0070f98d1a30a2da0da02f0c8389f3322 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 17 Aug 2026 23:16:36 +0200 Subject: [PATCH 1/5] Add parallel intitialization for TLSPH arrays to ensure NUMA awareness --- .../structure/total_lagrangian_sph/system.jl | 22 +++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index a1cd48fb54..dc710682bf 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -145,10 +145,24 @@ function TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing poisson_ratio_sorted = poisson_ratio end - initial_coordinates = copy(initial_condition_sorted.coordinates) - current_coordinates = copy(initial_condition_sorted.coordinates) - mass = copy(initial_condition_sorted.mass) - material_density = copy(initial_condition_sorted.density) + initial_coordinates = similar(initial_condition_sorted.coordinates) + current_coordinates = similar(initial_condition_sorted.coordinates) + mass = similar(initial_condition_sorted.mass) + material_density = similar(initial_condition_sorted.density) + + # Initialize the runtime arrays in parallel so that first-touch allocation places + # their memory close to the threads that will access it during the simulation. + # The semidiscretization backend is not available in the constructor yet, and using + # Polyester here is harmless for serial and GPU simulations. + parallelization_backend = PolyesterBackend() + copyto!(ThreadedBroadcastArray(initial_coordinates; parallelization_backend), + initial_condition_sorted.coordinates) + copyto!(ThreadedBroadcastArray(current_coordinates; parallelization_backend), + initial_condition_sorted.coordinates) + copyto!(ThreadedBroadcastArray(mass; parallelization_backend), + initial_condition_sorted.mass) + copyto!(ThreadedBroadcastArray(material_density; parallelization_backend), + initial_condition_sorted.density) correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) pk1_rho2 = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) deformation_grad = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) From fa834b374197b27c5ef0693af8668f9f04db43b9 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 17 Aug 2026 23:28:36 +0200 Subject: [PATCH 2/5] Fix --- .../structure/total_lagrangian_sph/system.jl | 17 ++++++++-------- src/util.jl | 20 +++++++++++++++++++ 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index dc710682bf..d64db95efc 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -154,15 +154,16 @@ function TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing # their memory close to the threads that will access it during the simulation. # The semidiscretization backend is not available in the constructor yet, and using # Polyester here is harmless for serial and GPU simulations. + # This makes the RHS significantly faster on large data center CPUs + # (see http://github.com/trixi-framework/TrixiParticles.jl/pull/1294). parallelization_backend = PolyesterBackend() - copyto!(ThreadedBroadcastArray(initial_coordinates; parallelization_backend), - initial_condition_sorted.coordinates) - copyto!(ThreadedBroadcastArray(current_coordinates; parallelization_backend), - initial_condition_sorted.coordinates) - copyto!(ThreadedBroadcastArray(mass; parallelization_backend), - initial_condition_sorted.mass) - copyto!(ThreadedBroadcastArray(material_density; parallelization_backend), - initial_condition_sorted.density) + copyto_threaded!(initial_coordinates, initial_condition_sorted.coordinates, + parallelization_backend) + copyto_threaded!(current_coordinates, initial_condition_sorted.coordinates, + parallelization_backend) + copyto_threaded!(mass, initial_condition_sorted.mass, parallelization_backend) + copyto_threaded!(material_density, initial_condition_sorted.density, + parallelization_backend) correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) pk1_rho2 = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) deformation_grad = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) diff --git a/src/util.jl b/src/util.jl index 650a2d0e88..03c93029b6 100644 --- a/src/util.jl +++ b/src/util.jl @@ -301,3 +301,23 @@ function Base.similar(::Broadcast.Broadcasted{ThreadedBroadcastStyle{P}}, # TODO we only have the type `P` here and just assume that we can do `P()` return ThreadedBroadcastArray(similar(Array{T}, dims), parallelization_backend=P()) end + +# Like `copyto!`, but parallelize over particles to use the same memory access pattern +# as the particle loops in the right-hand side. +@inline function copyto_threaded!(dest::Matrix, src, parallelization_backend) + @threaded parallelization_backend for particle in axes(dest, 2) + for dimension in axes(dest, 1) + @inbounds dest[dimension, particle] = src[dimension, particle] + end + end + + return dest +end + +@inline function copyto_threaded!(dest::Vector, src, parallelization_backend) + @threaded parallelization_backend for particle in eachindex(dest) + @inbounds dest[particle] = src[particle] + end + + return dest +end From 2a94489207939bb22b789b8d6d81e4eeb3e2c5f5 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 17 Aug 2026 23:29:35 +0200 Subject: [PATCH 3/5] Update comment --- src/schemes/structure/total_lagrangian_sph/system.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index d64db95efc..a0530de508 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -154,8 +154,8 @@ function TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing # their memory close to the threads that will access it during the simulation. # The semidiscretization backend is not available in the constructor yet, and using # Polyester here is harmless for serial and GPU simulations. - # This makes the RHS significantly faster on large data center CPUs - # (see http://github.com/trixi-framework/TrixiParticles.jl/pull/1294). + # This makes the RHS significantly faster on large data center CPUs with multiple + # NUMA domains (see http://github.com/trixi-framework/TrixiParticles.jl/pull/1294). parallelization_backend = PolyesterBackend() copyto_threaded!(initial_coordinates, initial_condition_sorted.coordinates, parallelization_backend) From f9ca7957cc6e8f43283f338486374e0e2e597d67 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 17 Aug 2026 23:34:21 +0200 Subject: [PATCH 4/5] Add NEWS.md entry --- NEWS.md | 5 +++++ src/schemes/structure/total_lagrangian_sph/system.jl | 5 +++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index ad3e44c715..81fc27af62 100644 --- a/NEWS.md +++ b/NEWS.md @@ -18,6 +18,11 @@ used in the Julia ecosystem. Notable changes will be documented in this file for - Added the number of split integration time steps to the `InfoCallback` output when a `SplitIntegrationCallback` is used (#1194). +### Performance + +- Improved multithreaded TLSPH performance on NUMA systems by initializing runtime + arrays in parallel (#1294). + ### Important Bugfixes - Fixed signed-distance constraints in `ParticlePackingSystem` when using a separate diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index a0530de508..339306a1e5 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -150,8 +150,9 @@ function TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing mass = similar(initial_condition_sorted.mass) material_density = similar(initial_condition_sorted.density) - # Initialize the runtime arrays in parallel so that first-touch allocation places - # their memory close to the threads that will access it during the simulation. + # Initialize the runtime arrays in parallel so that, on NUMA systems, the first-touch + # allocation policy places their memory close to the threads that will access it + # during the simulation. # The semidiscretization backend is not available in the constructor yet, and using # Polyester here is harmless for serial and GPU simulations. # This makes the RHS significantly faster on large data center CPUs with multiple From 104fe9a537cd7e7349fddf46ff919dcfd4989c19 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 17 Aug 2026 23:48:57 +0200 Subject: [PATCH 5/5] Fix bounds checks --- src/util.jl | 6 +++++- test/general/util.jl | 16 ++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/util.jl b/src/util.jl index 03c93029b6..6ca7c9ec6a 100644 --- a/src/util.jl +++ b/src/util.jl @@ -305,6 +305,10 @@ end # Like `copyto!`, but parallelize over particles to use the same memory access pattern # as the particle loops in the right-hand side. @inline function copyto_threaded!(dest::Matrix, src, parallelization_backend) + if axes(dest) != axes(src) + throw(DimensionMismatch("source and destination must have the same axes")) + end + @threaded parallelization_backend for particle in axes(dest, 2) for dimension in axes(dest, 1) @inbounds dest[dimension, particle] = src[dimension, particle] @@ -315,7 +319,7 @@ end end @inline function copyto_threaded!(dest::Vector, src, parallelization_backend) - @threaded parallelization_backend for particle in eachindex(dest) + @threaded parallelization_backend for particle in eachindex(dest, src) @inbounds dest[particle] = src[particle] end diff --git a/test/general/util.jl b/test/general/util.jl index 6f96607012..701b1314ea 100644 --- a/test/general/util.jl +++ b/test/general/util.jl @@ -14,6 +14,22 @@ @test visited == [(1, 3), (2, 4)] end +@testset verbose=true "copyto_threaded!" begin + backend = SerialBackend() + + matrix = zeros(2, 3) + TrixiParticles.copyto_threaded!(matrix, reshape(1:6, 2, 3), backend) + @test matrix == reshape(1:6, 2, 3) + @test_throws DimensionMismatch TrixiParticles.copyto_threaded!(matrix, zeros(2, 2), + backend) + + vector = zeros(3) + TrixiParticles.copyto_threaded!(vector, 1:3, backend) + @test vector == 1:3 + @test_throws DimensionMismatch TrixiParticles.copyto_threaded!(vector, zeros(2), + backend) +end + @testset verbose=true "ThreadedBroadcastArray" begin A = TrixiParticles.ThreadedBroadcastArray(ones(3, 3)) B = ones(3, 3)