Multithreading bug fix

Project.toml

@@ -1,7 +1,7 @@
 name = "Peridynamics"
 uuid = "4dc47793-80f3-4232-b30e-ca78ca9d621b"
 authors = ["Kai Friebertshäuser"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d"

docs/src/logo.md

@@ -75,7 +75,7 @@ matₛ = BondBasedMaterial(;
 
 All material points of the spheres have a initial velocity of $-20\, \mathrm{m} \, \mathrm{s}^{-1}$ in $z$-direction.
 ```julia
-icₛ = [VelocityIC(-20.0, 1:(pcₛ₁.n_points), 3)]
+icₛ = [VelocityIC(-20.0, 1:pcₛ₁.n_points, 3)]
 ```
 For the contact analysis, every body needs to be specified with a [`BodySetup`](@ref).
 Time can be saved by using only one sphere for the calculation of the stable time step.

examples/Logo.jl

@@ -48,7 +48,7 @@ matₛ = BondBasedMaterial(;
 )
 
 ##--   4. DEFINE INITIAL CONDITIONS   ------------------------------------------------------
-icₛ = [VelocityIC(-20.0, 1:(pcₛ₁.n_points), 3)]
+icₛ = [VelocityIC(-20.0, 1:pcₛ₁.n_points, 3)]
 
 ##--   5. DEFINE THE BODY-SETUP   ----------------------------------------------------------
 plate = BodySetup(pcₚ, matₚ)

src/BondBased.jl

@@ -99,7 +99,7 @@ function BondBasedBody(mat::BondBasedMaterial, pc::PointCloud)
     owned_bonds = defaultdist(n_bonds, n_threads)
     bond_failure = ones(Int, n_bonds)
     n_active_family_members = zeros(Int, (n_points, n_threads))
-    @threads :static for tid in 1:n_threads
+    @threads for tid in 1:n_threads
         for i in owned_points[tid]
             n_active_family_members[i, tid] = n_family_members[i]
         end
@@ -138,8 +138,7 @@ create_simmodel(mat::BondBasedMaterial, pc::PointCloud) = BondBasedBody(mat, pc)
 function compute_forcedensity!(body::BondBasedBody, mat::BondBasedMaterial)
     body.b_int .= 0.0
     body.n_active_family_members .= 0
-    @threads :static for _ in 1:(body.n_threads)
-        tid = threadid()
+    @threads for tid in 1:body.n_threads
         for current_bond in body.owned_bonds[tid]
             (i, j, ξ, failureflag) = body.bond_data[current_bond]
             ϑx = body.position[1, j] - body.position[1, i]

src/ContactBase.jl

@@ -97,15 +97,15 @@ function submit(sim::PDContactAnalysis)
         log_describe_sim(sim)
         bodies = [create_simmodel(ps.mat, ps.pc) for ps in sim.body_setup]
         log_describe_interactions(bodies, sim.es)
-        for i in 1:(sim.n_bodies)
+        for i in 1:sim.n_bodies
             for precrack in sim.body_setup[i].precracks
                 define_precrack!(bodies[i], precrack)
             end
             calc_damage!(bodies[i])
         end
         if sim.td.Δt < 0.0
             _Δt = Vector{Float64}(undef, 0)
-            for i in 1:(sim.n_bodies)
+            for i in 1:sim.n_bodies
                 if sim.body_setup[i].calc_timestep
                     push!(
                         _Δt,
@@ -121,11 +121,11 @@ function submit(sim::PDContactAnalysis)
             sim.td.Δt = minimum(_Δt)
         end
         log_simsetup(sim)
-        for i in 1:(sim.n_bodies)
+        for i in 1:sim.n_bodies
             apply_ics!(bodies[i], sim.body_setup[i].ics)
         end
         if sim.es.exportflag
-            for i in 1:(sim.n_bodies)
+            for i in 1:sim.n_bodies
                 export_results(bodies[i], string(sim.es.resultfile_prefix, "_b", i), 0, 0.0)
             end
         end
@@ -137,7 +137,7 @@ end
 
 function log_describe_sim(sim::PDContactAnalysis)
     msg = "Peridynamic multibody contact analysis: " * sim.name * "\nMaterial parameters:\n"
-    for i in 1:(sim.n_bodies)
+    for i in 1:sim.n_bodies
         msg *= @sprintf("Body %d:\n", i)
         msg *= @sprintf(
             "  - Number of material points [-]:      %30g\n",
@@ -209,21 +209,21 @@ function velocity_verlet!(body::Vector{<:AbstractPDBody}, sim::PDContactAnalysis
         enabled=!is_logging(stderr),
     )
     Δt½ = 0.5 * sim.td.Δt
-    for t in 1:(sim.td.n_timesteps)
+    for t in 1:sim.td.n_timesteps
         time = t * sim.td.Δt
-        for i in 1:(sim.n_bodies)
+        for i in 1:sim.n_bodies
             update_velhalf!(body[i], Δt½)
             apply_bcs!(body[i], sim.body_setup[i].bcs, time)
             update_disp_and_position!(body[i], sim.td.Δt)
             compute_forcedensity!(body[i], sim.body_setup[i].mat)
             calc_damage!(body[i])
         end
         compute_contactforcedensity!(body, sim)
-        for i in 1:(sim.n_bodies)
+        for i in 1:sim.n_bodies
             compute_equation_of_motion!(body[i], Δt½, sim.body_setup[i].mat.rho)
         end
         if mod(t, sim.es.exportfreq) == 0
-            for i in 1:(sim.n_bodies)
+            for i in 1:sim.n_bodies
                 export_results(body[i], string(sim.es.resultfile_prefix, "_b", i), t, time)
             end
         end
@@ -242,10 +242,9 @@ function compute_contactforcedensity!(
         body_b = bodies[jj]
         r = contact.search_radius
         C = 9 * contact.spring_constant / (π * sim.body_setup[ii].mat.δ^4)
-        @threads :static for _ in 1:nthreads()
-            tid = threadid()
+        @threads for tid in 1:nthreads()
             for i in body_a.owned_points[tid]
-                for j in 1:(body_b.n_points)
+                for j in 1:body_b.n_points
                     ϑx = body_b.position[1, j] - body_a.position[1, i]
                     ϑy = body_b.position[2, j] - body_a.position[2, i]
                     ϑz = body_b.position[3, j] - body_a.position[3, i]

src/PeridynamicsBase.jl

@@ -431,16 +431,15 @@ function find_bonds(pc::PointCloud, δ::Float64, owned_points::Vector{UnitRange{
         color=:normal,
         enabled=!is_logging(stderr),
     )
-    @threads :static for _ in 1:n_threads
-        tid = threadid()
+    @threads for tid in 1:n_threads
         local_bond_data = Vector{Tuple{Int,Int,Float64,Bool}}(undef, 0)
         sizehint!(local_bond_data, pc.n_points * 500)
         idst = 0.0
         num = 0
         fail = false
         for a in owned_points[tid]
             num = 0
-            for i in 1:(pc.n_points)
+            for i in 1:pc.n_points
                 if a !== i
                     idst = sqrt(
                         (pc.position[1, i] - pc.position[1, a])^2 +
@@ -475,16 +474,15 @@ function find_unique_bonds(pc::PointCloud, δ::Float64, owned_points::Vector{Uni
         color=:normal,
         enabled=!is_logging(stderr),
     )
-    @threads :static for _ in 1:n_threads
-        tid = threadid()
+    @threads for tid in 1:n_threads
         local_bonds_data = Vector{Tuple{Int,Int,Float64,Bool}}(undef, 0)
         sizehint!(local_bonds_data, pc.n_points * 500)
         idst = 0.0
         num = 0
         fail = false
         for a in owned_points[tid]
             num = 0
-            for i in 1:(pc.n_points)
+            for i in 1:pc.n_points
                 if a !== i
                     idst = sqrt(
                         (pc.position[1, i] - pc.position[1, a])^2 +
@@ -531,7 +529,7 @@ function defaultdist(sz::Int, nc::Int)
         return grid
     else
         sidx = [1:(sz + 1);]
-        grid = fill(0:0, nc)
+        grid = fill(0:-1, nc)
         for i in 1:sz
             grid[i] = sidx[i]:(sidx[i + 1] - 1)
         end
@@ -737,8 +735,7 @@ function log_closesim(part::AbstractPDBody, timingsummary::NamedTuple, es::Expor
 end
 
 function define_precrack!(body::AbstractPDBody, precrack::PreCrack)
-    @threads :static for _ in 1:(body.n_threads)
-        tid = threadid()
+    @threads for tid in 1:body.n_threads
         (@view body.n_active_family_members[:, tid]) .= 0
         for current_one_ni in body.owned_bonds[tid]
             i, j, _, _ = body.bond_data[current_one_ni]
@@ -760,8 +757,7 @@ end
 
 function calc_stable_timestep(body::AbstractPDBody, rho::Float64, K::Float64, δ::Float64)
     _Δt = zeros(Float64, body.n_threads)
-    @threads :static for _ in 1:(body.n_threads)
-        tid = threadid()
+    @threads for tid in 1:body.n_threads
         timesteps = zeros(Float64, body.n_points)
         dtsum = zeros(Float64, (body.n_points, body.n_threads))
         for current_one_ni in body.owned_bonds[tid]
@@ -796,8 +792,7 @@ function calc_stable_user_timestep(pc::PointCloud, mat::AbstractPDMaterial, Sf::
     bond_data, _ = find_bonds(pc, mat.δ, owned_points)
     owned_bonds = defaultdist(length(bond_data), nthreads())
     _Δt = zeros(Float64, nthreads())
-    @threads :static for _ in 1:nthreads()
-        tid = threadid()
+    @threads for tid in 1:nthreads()
         timesteps = zeros(Float64, pc.n_points)
         dtsum = zeros(Float64, (pc.n_points, nthreads()))
         for current_one_ni in owned_bonds[tid]
@@ -815,8 +810,8 @@ function calc_stable_user_timestep(pc::PointCloud, mat::AbstractPDMaterial, Sf::
 end
 
 function calc_damage!(body::AbstractPDBody)
-    @threads :static for _ in 1:(body.n_threads)
-        for a in body.owned_points[threadid()]
+    @threads for tid in 1:body.n_threads
+        for a in body.owned_points[tid]
             body.n_active_family_members[a, 1] = sum(
                 @view body.n_active_family_members[a, :]
             )
@@ -897,7 +892,7 @@ end
 
 function export_vtk(body::AbstractPDBody, expfile::String, timestep::Int, time::Float64)
     filename = expfile * "_t" * string(timestep)
-    cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:(body.n_points)]
+    cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (j,)) for j in 1:body.n_points]
     vtkfile = vtk_grid(filename, body.position, cells)
     vtkfile["Damage", VTKPointData()] = body.damage
     vtkfile["Force density", VTKPointData()] = @views body.b_int[:, :, 1]
@@ -923,7 +918,7 @@ function velocity_verlet!(body::AbstractPDBody, sim::PDSingleBodyAnalysis)
         enabled=!is_logging(stderr),
     )
     Δt½ = 0.5 * sim.td.Δt
-    for t in 1:(sim.td.n_timesteps)
+    for t in 1:sim.td.n_timesteps
         time = t * sim.td.Δt
         update_velhalf!(body, Δt½)
         apply_bcs!(body, sim.bcs, time)
@@ -956,7 +951,7 @@ function dynamic_relaxation_finite_difference!(
         color=:normal,
         enabled=!is_logging(stderr),
     )
-    for t in 1:(sim.td.n_timesteps)
+    for t in 1:sim.td.n_timesteps
         time = t * sim.td.Δt
         apply_bcs!(body, sim.bcs, time)
         compute_forcedensity!(body, sim.mat)
@@ -989,7 +984,7 @@ function calc_damping(
 )
     cn1 = 0.0
     cn2 = 0.0
-    for i in 1:(body.n_points)
+    for i in 1:body.n_points
         for d in 1:3
             body.b_int[d, i, 1] = sum(@view body.b_int[d, i, :])
             if velocity_half_old[d, i] !== 0.0
@@ -1022,7 +1017,7 @@ function finite_difference_first_step!(
     b_int_old::Matrix{Float64},
     Δt::Float64,
 )
-    @threads for i in 1:(body.n_points)
+    @threads for i in 1:body.n_points
         for d in 1:3
             body.velocity_half[d, i] =
                 0.5 * Δt / damping_matrix[d, i] * (body.b_int[d, i, 1] + body.b_ext[d, i])
@@ -1044,7 +1039,7 @@ function finite_difference!(
     Δt::Float64,
     cn::Float64,
 )
-    @threads for i in 1:(body.n_points)
+    @threads for i in 1:body.n_points
         for d in 1:3
             body.velocity_half[d, i] =
                 (
@@ -1063,7 +1058,7 @@ function finite_difference!(
 end
 
 function update_velhalf!(body::AbstractPDBody, Δt½::Float64)
-    @threads for i in 1:(body.n_points)
+    @threads for i in 1:body.n_points
         body.velocity_half[1, i] = body.velocity[1, i] + body.acceleration[1, i] * Δt½
         body.velocity_half[2, i] = body.velocity[2, i] + body.acceleration[2, i] * Δt½
         body.velocity_half[3, i] = body.velocity[3, i] + body.acceleration[3, i] * Δt½
@@ -1072,7 +1067,7 @@ function update_velhalf!(body::AbstractPDBody, Δt½::Float64)
 end
 
 function update_disp_and_position!(body::AbstractPDBody, Δt::Float64)
-    @threads for i in 1:(body.n_points)
+    @threads for i in 1:body.n_points
         body.displacement[1, i] += body.velocity_half[1, i] * Δt
         body.displacement[2, i] += body.velocity_half[2, i] * Δt
         body.displacement[3, i] += body.velocity_half[3, i] * Δt
@@ -1084,7 +1079,7 @@ function update_disp_and_position!(body::AbstractPDBody, Δt::Float64)
 end
 
 function compute_equation_of_motion!(body::AbstractPDBody, Δt½::Float64, rho::Float64)
-    @threads for i in 1:(body.n_points)
+    @threads for i in 1:body.n_points
         body.b_int[1, i, 1] = sum(@view body.b_int[1, i, :])
         body.b_int[2, i, 1] = sum(@view body.b_int[2, i, :])
         body.b_int[3, i, 1] = sum(@view body.b_int[3, i, :])

test/bondbased.jl

@@ -49,24 +49,32 @@ if Threads.nthreads() <= 2
     @test body.n_threads == Threads.nthreads()
     @test body.unique_bonds == true
     @test body.b_int == zeros(3, n_points, Threads.nthreads())
-    @test body.n_active_family_members == [1; 1;;]
+    if Threads.nthreads() == 1
+        @test body.n_active_family_members == [1; 1;;]
+    elseif Threads.nthreads() == 2
+        @test body.n_active_family_members == [1 0; 0 1]
+    end
     @test body.bond_failure == [1]
 
     # Boundary Condition:
     # Point 2 with ẋ_z = 1 m/s with Δt = 0.0015 s
     body.position[1, 2] = 1.0015
     Peridynamics.compute_forcedensity!(body, mat)
 
-    @test body.n_active_family_members == [0; 0;;]
+    if Threads.nthreads() == 1
+        @test body.n_active_family_members == [0; 0;;]
+    elseif Threads.nthreads() == 2
+        @test body.n_active_family_members == [0 0; 0 0]
+    end
     @test body.bond_failure == [0]
 
     b¹² = [
         18 * E / (3 * (1 - 2 * 0.25)) / (π * δ^4) * 1.0015 * 0.0015/1.0015 * 1.0
         0.0
         0.0
     ]
-    Threads.@threads for _ in 1:Threads.nthreads()
-        for i in body.owned_points[Threads.threadid()]
+    Threads.@threads for tid in 1:Threads.nthreads()
+        for i in body.owned_points[tid]
             body.b_int[1,i,1] = sum(@view body.b_int[1,i,:])
             body.b_int[2,i,1] = sum(@view body.b_int[2,i,:])
             body.b_int[3,i,1] = sum(@view body.b_int[3,i,:])

test/bonds.jl

@@ -37,14 +37,22 @@ if Threads.nthreads() <= 2
     body = Peridynamics.create_simmodel(mat, pc)
     precrack = PreCrack([1], [2])
 
-    @test body.n_active_family_members == [1; 1;;]
+    if Threads.nthreads() == 1
+        @test body.n_active_family_members == [1; 1;;]
+    elseif Threads.nthreads() == 2
+        @test body.n_active_family_members == [1 0; 0 1]
+    end
     @test body.bond_failure == [1]
     @test body.damage == [0, 0]
 
     Peridynamics.define_precrack!(body, precrack)
     Peridynamics.calc_damage!(body)
 
-    @test body.n_active_family_members == [0; 0;;]
+    if Threads.nthreads() == 1
+        @test body.n_active_family_members == [0; 0;;]
+    elseif Threads.nthreads() == 2
+        @test body.n_active_family_members == [0 0; 0 0]
+    end
     @test body.bond_failure == [0]
     @test body.damage == [1, 1]
 else

test/timeintegration.jl

@@ -91,7 +91,7 @@ if Threads.nthreads() <= 2
     @test body.velocity_half == [-0.09375 0.09375; 0.0 0.0; 0.0 0.0]
     @test body.displacement == [-0.09375 0.09375; 0.0 0.0; 0.0 0.0]
     @test body.position == [-0.09375 1.09375; 0.0 0.0; 0.0 0.0]
-    @test body.b_int == zeros(3, 2, 1)
+    @test body.b_int == zeros(3, 2, Threads.nthreads())
     @test body.acceleration == zeros(3, 2)
     @test body.velocity == [-0.046875 0.046875; 0.0 0.0; 0.0 0.0]
 else

test/utilities.jl

@@ -6,4 +6,4 @@ using Test
 
 @test Peridynamics.defaultdist(4, 2) == [1:2, 3:4]
 @test Peridynamics.defaultdist(3, 2) == [1:2, 3:3]
-@test Peridynamics.defaultdist(2, 3) == [1:1, 2:2, 0:0]
+@test Peridynamics.defaultdist(2, 3) == [1:1, 2:2, 0:-1]