From bcf2ad770897a2e5bdd70a31e08f6f663692530b Mon Sep 17 00:00:00 2001 From: Kai Partmann Date: Wed, 6 Sep 2023 10:11:54 +0200 Subject: [PATCH] PC2 refactoring to new data structure --- src/Peridynamics.jl | 6 +- src/bond_based.jl | 58 ++++--- src/conditions.jl | 155 ++--------------- src/dynamic_relaxation.jl | 2 +- src/io.jl | 48 +----- src/jobs.jl | 18 +- src/pdproblem.jl | 122 +++++++++----- src/precrack.jl | 29 ++++ src/spatial_discretization.jl | 50 ------ src/structs.jl | 187 +++++++++++++++++++++ src/utilities.jl | 21 +-- src/velocity_verlet.jl | 115 ++++++------- test/integration_tests/b_int_bond_based.jl | 2 +- 13 files changed, 434 insertions(+), 379 deletions(-) create mode 100644 src/precrack.jl create mode 100644 src/structs.jl diff --git a/src/Peridynamics.jl b/src/Peridynamics.jl index 1e55aa87..a62012d7 100644 --- a/src/Peridynamics.jl +++ b/src/Peridynamics.jl @@ -25,11 +25,13 @@ pinthreads(:cores; force=false) # Basics for peridynamics and contact problems include("abstract_types.jl") +include("structs.jl") include("multi_material.jl") include("spatial_discretization.jl") -include("conditions.jl") -include("io.jl") include("pdproblem.jl") +include("io.jl") +include("conditions.jl") +include("precrack.jl") include("velocity_verlet.jl") # include("dynamic_relaxation.jl") include("jobs.jl") diff --git a/src/bond_based.jl b/src/bond_based.jl index db1d9945..697e16a2 100644 --- a/src/bond_based.jl +++ b/src/bond_based.jl @@ -185,34 +185,48 @@ end return BondBasedBody(mat, pc) end -@timeit TO function compute_forcedensity!(body::BondBasedBody, mat::PDMaterial{BondBasedMaterial}) - body.b_int .= 0.0 - body.n_active_family_members .= 0 - @inbounds @threads :static 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] - ϑy = body.position[2, j] - body.position[2, i] - ϑz = body.position[3, j] - body.position[3, i] +@timeit TO function compute_forcedensity!(pdp::PDProblem) + @threads :static for tid in 1:pdp.sp.n_threads + tls = pdp.tls[tid] + tls.b_int .= 0 + tls.n_active_family_members .= 0 + for current_bond in pdp.sp.owned_bonds[tid] + (i, j, ξ, failureflag) = pdp.sp.bond_data[current_bond] + li = tls.pointmap[i] + lj = tls.pointmap[j] + ϑx = pdp.gs.position[1, j] - pdp.gs.position[1, i] + ϑy = pdp.gs.position[2, j] - pdp.gs.position[2, i] + ϑz = pdp.gs.position[3, j] - pdp.gs.position[3, i] η = sqrt(ϑx * ϑx + ϑy * ϑy + ϑz * ϑz) ε = (η - ξ) / ξ - temp = body.bond_failure[current_bond] * ε / η - tempi = temp * mat[i].bc * body.volume[j] - tempj = temp * mat[j].bc * body.volume[i] - body.b_int[1, i, tid] += tempi * ϑx - body.b_int[2, i, tid] += tempi * ϑy - body.b_int[3, i, tid] += tempi * ϑz - body.b_int[1, j, tid] -= tempj * ϑx - body.b_int[2, j, tid] -= tempj * ϑy - body.b_int[3, j, tid] -= tempj * ϑz - if ε > mat[i].εc || ε > mat[j].εc + temp = pdp.gs.bond_failure[current_bond] * ε / η + tempi = temp * pdp.sp.mat.bc * pdp.sp.pc.volume[j] + tempj = temp * pdp.sp.mat.bc * pdp.sp.pc.volume[i] + # tls.b_int[1, li] += tempi * ϑx + # tls.b_int[2, li] += tempi * ϑy + # tls.b_int[3, li] += tempi * ϑz + # tls.b_int[1, lj] -= tempj * ϑx + # tls.b_int[2, lj] -= tempj * ϑy + # tls.b_int[3, lj] -= tempj * ϑz + update_tls_b_int!(tls.b_int, 1, li, tempi * ϑx) + update_tls_b_int!(tls.b_int, 2, li, tempi * ϑy) + update_tls_b_int!(tls.b_int, 3, li, tempi * ϑz) + update_tls_b_int!(tls.b_int, 1, lj, tempj * ϑx) + update_tls_b_int!(tls.b_int, 2, lj, tempj * ϑy) + update_tls_b_int!(tls.b_int, 3, lj, tempj * ϑz) + + if ε > pdp.sp.mat.εc || ε > pdp.sp.mat.εc if failureflag - body.bond_failure[current_bond] = 0 + pdp.gs.bond_failure[current_bond] = 0 end end - body.n_active_family_members[i, tid] += body.bond_failure[current_bond] - body.n_active_family_members[j, tid] += body.bond_failure[current_bond] + tls.n_active_family_members[li] += pdp.gs.bond_failure[current_bond] + tls.n_active_family_members[lj] += pdp.gs.bond_failure[current_bond] end end return nothing end + +function update_tls_b_int!(b_int, d, li, val) + b_int[d, li] += val +end diff --git a/src/conditions.jl b/src/conditions.jl index ad602025..769047e6 100644 --- a/src/conditions.jl +++ b/src/conditions.jl @@ -1,173 +1,48 @@ -function Base.show(io::IO, ::MIME"text/plain", bc::T) where {T <: AbstractBC} - msg = string(length(bc.point_id_set), "-points ", T) - if bc.dim == 1 - msg *= " in x-direction (dim=1)" - elseif bc.dim == 2 - msg *= " in y-direction (dim=2)" - elseif bc.dim == 3 - msg *= " in z-direction (dim=3)" +@timeit TO function apply_bcs!(pdp::PDProblem, t::Float64) + @sync for bc in pdp.sp.bcs + Threads.@spawn apply_boundarycondition!(pdp.gs, bc, t) end - print(io, msg) return nothing end -function Base.show(io::IO, ::MIME"text/plain", ic::T) where {T <: AbstractIC} - msg = string(length(ic.point_id_set), "-points ", T) - if ic.dim == 1 - msg *= " in x-direction (dim=1)" - elseif ic.dim == 2 - msg *= " in y-direction (dim=2)" - elseif ic.dim == 3 - msg *= " in z-direction (dim=3)" - end - print(io, msg) - return nothing -end - -@doc raw""" - VelocityBC <: AbstractBC - -Velocity boundary condition. The value of the velocity is calculated every time step with -the function `fun` and applied to the dimension `dim` of the points specified by -`point_id_set`. - -# Fields -- `fun::Function`: function `f(t)` with current time `t` as argument, that calculates the -value of the velocity for each timestep -- `point_id_set::Vector{Int}`: point-id set with all points for which the boundary condition -gets applied every timestep -- `dim::Int`: dimension on which the boundary condition gets applied to. Possible values: -- x-direction: `dim=1` -- y-direction: `dim=2` -- z-direction: `dim=3` - -# Examples - -The constant velocity $v = 0.1$ in $y$-direction gets applied to the first $10$ points: -```julia-repl -julia> VelocityBC(t -> 0.1, 1:10, 2) -VelocityBC(var"#7#8"(), [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 2) -``` -""" -struct VelocityBC <: AbstractBC - fun::Function - point_id_set::Vector{Int} - dim::Int -end - -@doc raw""" - ForceDensityBC <: AbstractBC - -Force density boundary condition. The value of the force density is calculated every time -step with the function `fun` and applied to the dimension `dim` of the points specified by -`point_id_set`. - -# Fields -- `fun::Function`: function `f(t)` with current time `t` as argument, that calculates the -value of the force density for each timestep -- `point_id_set::Vector{Int}`: point-id set with all points for which the boundary condition -gets applied every timestep -- `dim::Int`: dimension on which the boundary condition gets applied to. Possible values: -- x-direction: `dim=1` -- y-direction: `dim=2` -- z-direction: `dim=3` -""" -struct ForceDensityBC <: AbstractBC - fun::Function - point_id_set::Vector{Int} - dim::Int -end - -@doc raw""" - VelocityIC <: AbstractIC - -Velocity initial condition. The value `val` of the velocity is applied as initial condition -to the dimension `dim` of the points specified by `point_id_set`. - -# Fields -- `val::Float64`: value of the velocity -- `point_id_set::Vector{Int}`: point-id set with all points for which the initial condition -gets applied to -- `dim::Int`: dimension on which the initial condition gets applied to. Possible values: -- x-direction: `dim=1` -- y-direction: `dim=2` -- z-direction: `dim=3` -""" -struct VelocityIC <: AbstractIC - val::Float64 - point_id_set::Vector{Int} - dim::Int -end - -@doc raw""" - PosDepVelBC <: AbstractBC - -Position dependent velocity boundary condition. The value of the force density is calculated -every time step with the function `fun` and applied to the dimension `dim` of the points -specified by `point_id_set`. - -# Fields -- `fun::Function`: function `f(x, y, z, t)` with x-, y-, z-position of the point and current -time `t` as arguments, calculates the value of the force density for each timestep -dependent of the point position -- `point_id_set::Vector{Int}`: point-id set with all points for which the boundary condition -gets applied to -- `dim::Int`: dimension on which the initial condition gets applied to. Possible values: -- x-direction: `dim=1` -- y-direction: `dim=2` -- z-direction: `dim=3` -""" -struct PosDepVelBC <: AbstractBC - fun::Function - point_id_set::Vector{Int} - dim::Int -end - -@timeit TO function apply_bcs!(body::AbstractPDBody, bcs::Vector{<:AbstractBC}, t::Float64) - @sync for bc in bcs - Threads.@spawn apply_boundarycondition!(body, bc, t) - end - return nothing -end - -function apply_boundarycondition!(body::AbstractPDBody, bc::VelocityBC, t::Float64) +function apply_boundarycondition!(gs::GlobalStorage, bc::VelocityBC, t::Float64) value = bc.fun(t) dim = bc.dim @inbounds @simd for i in bc.point_id_set - body.velocity_half[dim, i] = value + gs.velocity_half[dim, i] = value end return nothing end -function apply_boundarycondition!(body::AbstractPDBody, bc::ForceDensityBC, t::Float64) +function apply_boundarycondition!(gs::GlobalStorage, bc::ForceDensityBC, t::Float64) value = bc.fun(t) dim = bc.dim @inbounds @simd for i in bc.point_id_set - body.b_ext[dim, i] = value + gs.b_ext[dim, i] = value end return nothing end -function apply_boundarycondition!(body::AbstractPDBody, bc::PosDepVelBC, t::Float64) +function apply_boundarycondition!(gs::GlobalStorage, bc::PosDepVelBC, t::Float64) dim = bc.dim @inbounds @simd for i in bc.point_id_set - body.velocity_half[dim, i] = bc.fun(body.position[1, i], body.position[2, i], - body.position[3, i], t) + gs.velocity_half[dim, i] = bc.fun(gs.position[1, i], gs.position[2, i], + gs.position[3, i], t) end return nothing end -@timeit TO function apply_ics!(body::AbstractPDBody, ics::Vector{<:AbstractIC}) - @sync for ic in ics - Threads.@spawn apply_initialcondition!(body, ic) +@timeit TO function apply_ics!(pdp::PDProblem) + @sync for ic in pdp.sp.ics + Threads.@spawn apply_initialcondition!(pdp.gs, ic) end return nothing end -function apply_initialcondition!(body::AbstractPDBody, ic::VelocityIC) +function apply_initialcondition!(gs::GlobalStorage, ic::VelocityIC) dim = ic.dim @inbounds @simd for i in ic.point_id_set - body.velocity[dim, i] = ic.val + gs.velocity[dim, i] = ic.val end return nothing end diff --git a/src/dynamic_relaxation.jl b/src/dynamic_relaxation.jl index ba7b6d1c..bc7ca195 100644 --- a/src/dynamic_relaxation.jl +++ b/src/dynamic_relaxation.jl @@ -48,7 +48,7 @@ function time_loop!(body::AbstractPDBody, dr::DynamicRelaxation, mat::PDMaterial apply_bcs!(body, bcs, time) update_disp_and_position!(body, dr.Δt) compute_forcedensity!(body, mat) - update_thread_cache!(body) + reduce_tls_to_gs!(body) calc_damage!(body) cn = calc_damping(body, damping_matrix, velocity_half_old, b_int_old, dr.Δt) if t == 1 diff --git a/src/io.jl b/src/io.jl index 90eb36c0..fe137ca7 100644 --- a/src/io.jl +++ b/src/io.jl @@ -1,42 +1,4 @@ -""" - ExportSettings -Export settings. - -# Fields -- `path::String`: path where results will be saved -- `exportfreq::Int`: export frequency, will export every `exportfreq`-th timestep -- `resultfile_prefix::String`: prefix of the result-filename -- `logfile::String`: name of logfile -- `exportflag::Bool`: disable export for a simulation where saved results are not needed - ---- -```julia -ExportSettings([path::String, freq::Int]) -``` - -Create `ExportSettings` only by `path` and `freq`. If no arguments are specified, the -`exportflag` will be set to `false` and export disabled. - -# Arguments -- `path::String`: path where results will be saved -- `freq::Int`: export frequency -""" -mutable struct ExportSettings - path::String - exportfreq::Int - resultfile_prefix::String - logfile::String - exportflag::Bool -end - -ExportSettings() = ExportSettings("", 0, "", "", false) -ExportSettings(path::String, freq::Int) = ExportSettings(path, freq, "", "", true) - -function Base.show(io::IO, ::MIME"text/plain", es::ExportSettings) - print(io, typeof(es)) - return nothing -end function log_header(es::ExportSettings) msg = "="^70 * "\n" @@ -172,12 +134,12 @@ function log_displacement_and_damage(body::AbstractPDBody) return msg end -@timeit TO function export_vtk(body::AbstractPDBody, expfile::String, timestep::Int, time::Float64) - filename = @sprintf("%s_t%04d", expfile, timestep) - vtk_grid(filename, body.position, body.cells) do vtk - vtk["Damage", VTKPointData()] = body.damage +@timeit TO function export_vtk(pdp::PDProblem, timestep::Int, time::Float64) + filename = @sprintf("%s_t%04d", pdp.sp.es.resultfile_prefix, timestep) + vtk_grid(filename, pdp.sp.pc.position, pdp.sp.cells) do vtk + vtk["Damage", VTKPointData()] = pdp.gs.damage # vtk["ForceDensity", VTKPointData()] = @views body.b_int[:, :, 1] - vtk["Displacement", VTKPointData()] = body.displacement + vtk["Displacement", VTKPointData()] = pdp.gs.displacement # vtk["Velocity", VTKPointData()] = body.velocity vtk["Time", VTKFieldData()] = time end diff --git a/src/jobs.jl b/src/jobs.jl index 043fd0b9..405185dd 100644 --- a/src/jobs.jl +++ b/src/jobs.jl @@ -61,23 +61,23 @@ function submit(sim::PDSingleBodyAnalysis) @timeit TO "init time loop" begin log_header(sim.es) log_describe_sim(sim.name, sim.pc, sim.mat, sim.es) - body = init_body(sim.mat, sim.pc) - log_describe_interactions(body, sim.es) + pdp = init_pdp(sim) + # log_describe_interactions(body, sim.es) @timeit TO "define cracks" begin for precrack in sim.precracks - define_precrack!(body, precrack) + define_precrack!(pdp, precrack) end - update_thread_cache!(body) - calc_damage!(body) + reduce_tls_to_gs!(pdp) + calc_damage!(pdp) end @timeit TO "init time discretization" begin - init_time_discretization!(sim.td, body, sim.mat) + init_time_discretization!(pdp) end log_simsetup(sim) end - @timeit TO "time loop" time_loop!(body, sim.td, sim.mat, sim.bcs, sim.ics, sim.es) + @timeit TO "time loop" time_loop!(pdp) end - log_closesim(body, timingsummary, sim.es) + # log_closesim(body, timingsummary, sim.es) log_timers(sim.es) - return body + return pdp end diff --git a/src/pdproblem.jl b/src/pdproblem.jl index ed2b1a2a..061df065 100644 --- a/src/pdproblem.jl +++ b/src/pdproblem.jl @@ -7,13 +7,14 @@ struct SimParameters{M <: AbstractPDMaterial, T <: AbstractTimeDiscretization} name::String pc::PointCloud unique_bonds::Bool + n_threads::Int n_bonds::Int bond_data::Vector{Tuple{Int, Int, Float64, Bool}} n_family_members::Vector{Int} owned_points::Vector{UnitRange{Int}} owned_bonds::Vector{UnitRange{Int}} - single_tids::Vector{Tuple{Int, Int}} - multi_tids::Vector{Tuple{Int, Vector{Int}}} + # single_tids::Vector{Tuple{Int, Int}} + # multi_tids::Vector{Tuple{Int, Vector{Int}}} mat::M cells::Vector{MeshCell{VTKCellType, Tuple{Int64}}} precracks::Vector{PreCrack} @@ -27,25 +28,38 @@ end All things that are calculated and change globally ??? """ struct GlobalStorage - owned_points::Vector{UnitRange{Int}} - owned_bonds::Vector{UnitRange{Int}} - single_tids::Vector{Tuple{Int, Int}} - multi_tids::Vector{Tuple{Int, Vector{Int}}} position::Matrix{Float64} displacement::Matrix{Float64} velocity::Matrix{Float64} velocity_half::Matrix{Float64} acceleration::Matrix{Float64} - damage::Matrix{Float64} + damage::Vector{Float64} bond_failure::Vector{Int} + b_int::Matrix{Float64} + b_ext::Matrix{Float64} + n_active_family_members::Vector{Int} +end + +function GlobalStorage(n_points::Int, n_bonds::Int) + position = zeros(3, n_points) + displacement = zeros(3, n_points) + velocity = zeros(3, n_points) + velocity_half = zeros(3, n_points) + acceleration = zeros(3, n_points) + damage = zeros(Int, n_points) + bond_failure = ones(Int, n_bonds) + b_int = zeros(3, n_points) + b_ext = zeros(3, n_points) + n_active_family_members = zeros(n_points) + return GlobalStorage(position, displacement, velocity, velocity_half, acceleration, + damage, bond_failure, b_int, b_ext, n_active_family_members) end """ All things that are thread locally """ struct ThreadLocalStorage - pointidmap::Dict{Int, Int} - bondidmap::Dict{Int, Int} + pointmap::Dict{Int, Int} b_int::Matrix{Float64} n_active_family_members::Vector{Int} end @@ -53,13 +67,13 @@ end """ """ -struct PDProblem - sp::SimParameters +struct PDProblem{M <: AbstractPDMaterial, T <: AbstractTimeDiscretization} + sp::SimParameters{M,T} gs::GlobalStorage tls::Vector{ThreadLocalStorage} end -function init_pdproblem(sim::AbstractPDAnalysis) +function init_pdp(sim::AbstractPDAnalysis) n_threads = nthreads() name = sim.name pc = sim.pc @@ -71,42 +85,72 @@ function init_pdproblem(sim::AbstractPDAnalysis) precracks = sim.precracks unique_bonds = true owned_points = defaultdist(pc.n_points, n_threads) - point_mappings = find_mappings(owned_points) + # point_mappings = find_mappings(owned_points) # find the bonds bond_data, n_family_members = find_unique_bonds(pc, mat, owned_points) n_bonds = length(bond_data) owned_bonds = defaultdist(n_bonds, n_threads) - bond_mappings = find_mappings(owned_bonds) - n_active_family_members = zeros(Int, (pc.n_points, n_threads)) - _sum_tids = zeros(Bool, (n_points, n_threads)) - _sum_tids .= false + # n_active_family_members = zeros(Int, (pc.n_points, n_threads)) + # _sum_tids = zeros(Bool, (pc.n_points, n_threads)) + # _sum_tids .= false + # @threads for tid in 1:n_threads + # for current_bond in owned_bonds[tid] + # (i, j, _, _) = bond_data[current_bond] + # n_active_family_members[i, tid] += 1 + # n_active_family_members[j, tid] += 1 + # _sum_tids[i, tid] = true + # _sum_tids[j, tid] = true + # end + # end + # sum_tids = [findall(row) for row in eachrow(_sum_tids)] + # single_tids, multi_tids = find_tids(sum_tids) + + cells = get_cells(pc.n_points) + # position = pc.position + # displacement = zeros(Float64, (3, pc.n_points)) + # velocity = zeros(Float64, (3, pc.n_points)) + # velocity_half = zeros(Float64, (3, pc.n_points)) + # acceleration = zeros(Float64, (3, pc.n_points)) + # b_ext = zeros(Float64, (3, pc.n_points)) + # b_int = zeros(Float64, (3, pc.n_points)) + # damage = zeros(Int, pc.n_points) + + sp = SimParameters(name, pc, unique_bonds, n_threads, n_bonds, bond_data, n_family_members, + owned_points, owned_bonds, mat, cells, precracks, bcs, ics, td, es) + gs = GlobalStorage(pc.n_points, n_bonds) + + # @assert n_threads == length(point_mappings) + tls = Vector{ThreadLocalStorage}(undef, n_threads) @threads for tid in 1:n_threads - for current_bond in owned_bonds[tid] - (i, j, _, _) = bond_data[current_bond] - n_active_family_members[i, tid] += 1 - n_active_family_members[j, tid] += 1 - _sum_tids[i, tid] = true - _sum_tids[j, tid] = true + + all_points_of_tid = Set{Int}() + for cb in owned_bonds[tid] + i, j, _, _ = bond_data[cb] + push!(all_points_of_tid, i) + push!(all_points_of_tid, j) + end + + pointmap = Dict{Int,Int}() + n_local_points = 0 + for i in all_points_of_tid + n_local_points += 1 + pointmap[i] = n_local_points + end + @assert length(pointmap) == n_local_points + + b_int_local = zeros(3, n_local_points) + n_active_family_members_local = zeros(n_local_points) + for val in values(pointmap) + if val > n_local_points + @info "li > n_local_points" tid val n_local_points + end end + tls[tid] = ThreadLocalStorage(pointmap, b_int_local, n_active_family_members_local) end - sum_tids = [findall(row) for row in eachrow(_sum_tids)] - single_tids, multi_tids = find_tids(sum_tids) - - cells = get_cells(n_points) - position = pc.position - displacement = zeros(Float64, (3, n_points)) - velocity = zeros(Float64, (3, n_points)) - velocity_half = zeros(Float64, (3, n_points)) - acceleration = zeros(Float64, (3, n_points)) - b_ext = zeros(Float64, (3, n_points)) - b_int = zeros(Float64, (3, n_points)) - damage = zeros(Int, n_points) - pdparams = SimParameters(name, pc, unique_bonds, n_bonds, bond_data, n_family_members, - owned_points, owned_bonds, single_tids, multi_tids, mat, cells, - precracks, bcs, ics, td, es) + pdp = PDProblem(sp, gs, tls) - return pdparams + return pdp end diff --git a/src/precrack.jl b/src/precrack.jl new file mode 100644 index 00000000..960c3182 --- /dev/null +++ b/src/precrack.jl @@ -0,0 +1,29 @@ +@timeit TO function define_precrack!(pdp::PDProblem, precrack::PreCrack) + @inbounds @threads :static for tid in 1:pdp.sp.n_threads + pdp.tls[tid].n_active_family_members .= 0 + for current_one_ni in pdp.sp.owned_bonds[tid] + i, j, _, _ = pdp.sp.bond_data[current_one_ni] + i_is_in_set_a = in(i, precrack.point_id_set_a) + i_is_in_set_b = in(i, precrack.point_id_set_b) + j_is_in_set_a = in(j, precrack.point_id_set_a) + j_is_in_set_b = in(j, precrack.point_id_set_b) + if i_is_in_set_a && j_is_in_set_b || i_is_in_set_b && j_is_in_set_a + pdp.gs.bond_failure[current_one_ni] = 0 + end + li = pdp.tls[tid].pointmap[i] + pdp.tls[tid].n_active_family_members[li] += pdp.gs.bond_failure[current_one_ni] + if pdp.sp.unique_bonds + lj = pdp.tls[tid].pointmap[j] + pdp.tls[tid].n_active_family_members[lj] += pdp.gs.bond_failure[current_one_ni] + end + end + end + return nothing +end + +@timeit TO function calc_damage!(pdp::PDProblem) + @inbounds @threads :static for i in 1:pdp.sp.pc.n_points + pdp.gs.damage[i] = 1 - pdp.gs.n_active_family_members[i] / pdp.sp.n_family_members[i] + end + return nothing +end diff --git a/src/spatial_discretization.jl b/src/spatial_discretization.jl index e121ccba..3478a039 100644 --- a/src/spatial_discretization.jl +++ b/src/spatial_discretization.jl @@ -273,56 +273,6 @@ function pcmerge(v::Vector{PointCloud}) return pc end -""" - PreCrack(point_id_set_a::Vector{Int}, point_id_set_b::Vector{Int}) - -Definition of an preexisting crack in the model. Points in `point_id_set_a` cannot have -interactions with points in `point_id_set_b`. - -# Fields -- `point_id_set_a::Vector{Int}`: first point-id set -- `point_id_set_b::Vector{Int}`: second point-id set -""" -struct PreCrack - point_id_set_a::Vector{Int} - point_id_set_b::Vector{Int} -end - -function Base.show(io::IO, ::MIME"text/plain", precrack::PreCrack) - println(io, typeof(precrack), ":") - println(io, " ", length(precrack.point_id_set_a), " points in set a") - print(io, " ", length(precrack.point_id_set_b), " points in set b") - return nothing -end - -@timeit TO function define_precrack!(body::AbstractPDBody, precrack::PreCrack) - @inbounds @threads :static for tid in 1:body.n_threads - body.n_active_family_members[:, tid] .= 0 - for current_one_ni in body.owned_bonds[tid] - i, j, _, _ = body.bond_data[current_one_ni] - i_is_in_set_a = in(i, precrack.point_id_set_a) - i_is_in_set_b = in(i, precrack.point_id_set_b) - j_is_in_set_a = in(j, precrack.point_id_set_a) - j_is_in_set_b = in(j, precrack.point_id_set_b) - if i_is_in_set_a && j_is_in_set_b || i_is_in_set_b && j_is_in_set_a - body.bond_failure[current_one_ni] = 0 - end - body.n_active_family_members[i, tid] += body.bond_failure[current_one_ni] - if body.unique_bonds - body.n_active_family_members[j, tid] += body.bond_failure[current_one_ni] - end - end - end - return nothing -end - -@timeit TO function calc_damage!(body::AbstractPDBody) - @inbounds @threads :static for i in 1:body.n_points - body.damage[i] = 1 - body.n_active_family_members[i, 1] / body.n_family_members[i] - end - return nothing -end - @timeit TO function find_bonds(pc::PointCloud, mat::PDMaterial, owned_points::Vector{UnitRange{Int}}) n_threads = nthreads() _bond_data = fill([(0, 0, 0.0, true)], n_threads) diff --git a/src/structs.jl b/src/structs.jl new file mode 100644 index 00000000..09d6572c --- /dev/null +++ b/src/structs.jl @@ -0,0 +1,187 @@ + +""" + PreCrack(point_id_set_a::Vector{Int}, point_id_set_b::Vector{Int}) + +Definition of an preexisting crack in the model. Points in `point_id_set_a` cannot have +interactions with points in `point_id_set_b`. + +# Fields +- `point_id_set_a::Vector{Int}`: first point-id set +- `point_id_set_b::Vector{Int}`: second point-id set +""" +struct PreCrack + point_id_set_a::Vector{Int} + point_id_set_b::Vector{Int} +end + +function Base.show(io::IO, ::MIME"text/plain", precrack::PreCrack) + println(io, typeof(precrack), ":") + println(io, " ", length(precrack.point_id_set_a), " points in set a") + print(io, " ", length(precrack.point_id_set_b), " points in set b") + return nothing +end + +function Base.show(io::IO, ::MIME"text/plain", bc::T) where {T <: AbstractBC} + msg = string(length(bc.point_id_set), "-points ", T) + if bc.dim == 1 + msg *= " in x-direction (dim=1)" + elseif bc.dim == 2 + msg *= " in y-direction (dim=2)" + elseif bc.dim == 3 + msg *= " in z-direction (dim=3)" + end + print(io, msg) + return nothing +end + +function Base.show(io::IO, ::MIME"text/plain", ic::T) where {T <: AbstractIC} + msg = string(length(ic.point_id_set), "-points ", T) + if ic.dim == 1 + msg *= " in x-direction (dim=1)" + elseif ic.dim == 2 + msg *= " in y-direction (dim=2)" + elseif ic.dim == 3 + msg *= " in z-direction (dim=3)" + end + print(io, msg) + return nothing +end + +@doc raw""" + VelocityBC <: AbstractBC + +Velocity boundary condition. The value of the velocity is calculated every time step with +the function `fun` and applied to the dimension `dim` of the points specified by +`point_id_set`. + +# Fields +- `fun::Function`: function `f(t)` with current time `t` as argument, that calculates the +value of the velocity for each timestep +- `point_id_set::Vector{Int}`: point-id set with all points for which the boundary condition +gets applied every timestep +- `dim::Int`: dimension on which the boundary condition gets applied to. Possible values: +- x-direction: `dim=1` +- y-direction: `dim=2` +- z-direction: `dim=3` + +# Examples + +The constant velocity $v = 0.1$ in $y$-direction gets applied to the first $10$ points: +```julia-repl +julia> VelocityBC(t -> 0.1, 1:10, 2) +VelocityBC(var"#7#8"(), [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 2) +``` +""" +struct VelocityBC <: AbstractBC + fun::Function + point_id_set::Vector{Int} + dim::Int +end + +@doc raw""" + ForceDensityBC <: AbstractBC + +Force density boundary condition. The value of the force density is calculated every time +step with the function `fun` and applied to the dimension `dim` of the points specified by +`point_id_set`. + +# Fields +- `fun::Function`: function `f(t)` with current time `t` as argument, that calculates the +value of the force density for each timestep +- `point_id_set::Vector{Int}`: point-id set with all points for which the boundary condition +gets applied every timestep +- `dim::Int`: dimension on which the boundary condition gets applied to. Possible values: +- x-direction: `dim=1` +- y-direction: `dim=2` +- z-direction: `dim=3` +""" +struct ForceDensityBC <: AbstractBC + fun::Function + point_id_set::Vector{Int} + dim::Int +end + +@doc raw""" + VelocityIC <: AbstractIC + +Velocity initial condition. The value `val` of the velocity is applied as initial condition +to the dimension `dim` of the points specified by `point_id_set`. + +# Fields +- `val::Float64`: value of the velocity +- `point_id_set::Vector{Int}`: point-id set with all points for which the initial condition +gets applied to +- `dim::Int`: dimension on which the initial condition gets applied to. Possible values: +- x-direction: `dim=1` +- y-direction: `dim=2` +- z-direction: `dim=3` +""" +struct VelocityIC <: AbstractIC + val::Float64 + point_id_set::Vector{Int} + dim::Int +end + +@doc raw""" + PosDepVelBC <: AbstractBC + +Position dependent velocity boundary condition. The value of the force density is calculated +every time step with the function `fun` and applied to the dimension `dim` of the points +specified by `point_id_set`. + +# Fields +- `fun::Function`: function `f(x, y, z, t)` with x-, y-, z-position of the point and current +time `t` as arguments, calculates the value of the force density for each timestep +dependent of the point position +- `point_id_set::Vector{Int}`: point-id set with all points for which the boundary condition +gets applied to +- `dim::Int`: dimension on which the initial condition gets applied to. Possible values: +- x-direction: `dim=1` +- y-direction: `dim=2` +- z-direction: `dim=3` +""" +struct PosDepVelBC <: AbstractBC + fun::Function + point_id_set::Vector{Int} + dim::Int +end + +""" + ExportSettings + +Export settings. + +# Fields +- `path::String`: path where results will be saved +- `exportfreq::Int`: export frequency, will export every `exportfreq`-th timestep +- `resultfile_prefix::String`: prefix of the result-filename +- `logfile::String`: name of logfile +- `exportflag::Bool`: disable export for a simulation where saved results are not needed + +--- +```julia +ExportSettings([path::String, freq::Int]) +``` + +Create `ExportSettings` only by `path` and `freq`. If no arguments are specified, the +`exportflag` will be set to `false` and export disabled. + +# Arguments +- `path::String`: path where results will be saved +- `freq::Int`: export frequency +""" +mutable struct ExportSettings + path::String + exportfreq::Int + resultfile_prefix::String + logfile::String + exportflag::Bool +end + +ExportSettings() = ExportSettings("", 0, "", "", false) +ExportSettings(path::String, freq::Int) = ExportSettings(path, freq, "", "", true) + +function Base.show(io::IO, ::MIME"text/plain", es::ExportSettings) + print(io, typeof(es)) + return nothing +end diff --git a/src/utilities.jl b/src/utilities.jl index 06ed08e6..eafdf9f4 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -42,21 +42,16 @@ end return single_tids, multi_tids end -@timeit TO function update_thread_cache!(body::AbstractPDBody) - @threads :static for (i, tid) in body.single_tids - body.b_int[1, i, 1] = body.b_int[1, i, tid] - body.b_int[2, i, 1] = body.b_int[2, i, tid] - body.b_int[3, i, 1] = body.b_int[3, i, tid] - body.n_active_family_members[i, 1] = body.n_active_family_members[i, tid] - end - @threads :static for (i, tids) in body.multi_tids - for tid in tids - body.b_int[1, i, 1] += body.b_int[1, i, tid] - body.b_int[2, i, 1] += body.b_int[2, i, tid] - body.b_int[3, i, 1] += body.b_int[3, i, tid] - body.n_active_family_members[i, 1] += body.n_active_family_members[i, tid] +@timeit TO function reduce_tls_to_gs!(pdp::PDProblem) + for tid in 1:pdp.sp.n_threads + for (gi, li) in pdp.tls[tid].pointmap + pdp.gs.b_int[1, gi] = pdp.tls[tid].b_int[1, li] + pdp.gs.b_int[2, gi] = pdp.tls[tid].b_int[2, li] + pdp.gs.b_int[3, gi] = pdp.tls[tid].b_int[3, li] + pdp.gs.n_active_family_members[gi] = pdp.tls[tid].n_active_family_members[li] end end + return nothing end is_logging(io) = isa(io, Base.TTY) == false || (get(ENV, "CI", nothing) == "true") diff --git a/src/velocity_verlet.jl b/src/velocity_verlet.jl index fdcc4a86..4f6274da 100644 --- a/src/velocity_verlet.jl +++ b/src/velocity_verlet.jl @@ -25,91 +25,88 @@ mutable struct VelocityVerlet <: AbstractTimeDiscretization end end -function init_time_discretization!(vv::VelocityVerlet, body::AbstractPDBody, - mat::PDMaterial) - if vv.Δt < 0 - Δt = calc_stable_timestep(body, mat, vv.safety_factor) - vv.Δt = Δt +function init_time_discretization!(pdp::PDProblem) + if pdp.sp.td.Δt < 0 + Δt = calc_stable_timestep(pdp) + pdp.sp.td.Δt = Δt end return nothing end -function calc_stable_timestep(body::AbstractPDBody, mat::PDMaterial, safety_factor::Float64) - _Δt = zeros(Float64, body.n_threads) - @inbounds @threads :static 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] - (a, i, L, _) = body.bond_data[current_one_ni] - dtsum[a, tid] += body.volume[i] * 1 / L * 18 * mat[a].K / (π * mat[a].δ^4) +function calc_stable_timestep(pdp::PDProblem) + _Δt = zeros(Float64, pdp.sp.n_threads) + @inbounds @threads :static for tid in 1:pdp.sp.n_threads + timesteps = zeros(Float64, pdp.sp.pc.n_points) + dtsum = zeros(Float64, (pdp.sp.pc.n_points, pdp.sp.n_threads)) + for current_one_ni in pdp.sp.owned_bonds[tid] + (a, i, L, _) = pdp.sp.bond_data[current_one_ni] + dtsum[a, tid] += pdp.sp.pc.volume[i] * 1 / L * 18 * pdp.sp.mat[a].K / (π * pdp.sp.mat[a].δ^4) end - for a in body.owned_points[tid] + for a in pdp.sp.owned_points[tid] dtsum[a, 1] = sum(@view dtsum[a, :]) - timesteps[a] = √(2 * mat[a].rho / dtsum[a, 1]) + timesteps[a] = √(2 * pdp.sp.mat[a].rho / dtsum[a, 1]) end - _Δt[tid] = safety_factor * minimum(timesteps[timesteps .> 0]) + _Δt[tid] = pdp.sp.td.safety_factor * minimum(timesteps[timesteps .> 0]) end Δt = minimum(_Δt) return Δt end -function time_loop!(body::AbstractPDBody, vv::VelocityVerlet, mat::PDMaterial, - bcs::Vector{<:AbstractBC}, ics::Vector{<:AbstractIC}, - es::ExportSettings) - apply_ics!(body, ics) - if es.exportflag - export_vtk(body, es.resultfile_prefix, 0, 0.0) +function time_loop!(pdp::PDProblem) + apply_ics!(pdp) + if pdp.sp.es.exportflag + export_vtk(pdp, 0, 0.0) end - # p = Progress(vv.n_steps; dt = 1, desc = "Time integration... ", barlen = 30, - # color = :normal, enabled = !is_logging(stderr)) - Δt½ = 0.5 * vv.Δt - for t in 1:(vv.n_steps) - time = t * vv.Δt - update_velhalf!(body, Δt½) - apply_bcs!(body, bcs, time) - update_disp_and_position!(body, vv.Δt) - compute_forcedensity!(body, mat) - update_thread_cache!(body) #TODO: BAD NAME! - calc_damage!(body) - compute_equation_of_motion!(body, Δt½, mat) - if mod(t, es.exportfreq) == 0 - export_vtk(body, es.resultfile_prefix, t, time) + p = Progress(pdp.sp.td.n_steps; dt = 1, desc = "Time integration... ", barlen = 30, + color = :normal, enabled = !is_logging(stderr)) + Δt½ = 0.5 * pdp.sp.td.Δt + for t in 1:pdp.sp.td.n_steps + time = t * pdp.sp.td.Δt + update_velhalf!(pdp.gs, Δt½) + apply_bcs!(pdp, time) + update_disp_and_position!(pdp.gs, pdp.sp.td.Δt) + compute_forcedensity!(pdp) + reduce_tls_to_gs!(pdp) + calc_damage!(pdp) + compute_equation_of_motion!(pdp, Δt½) + if mod(t, pdp.sp.es.exportfreq) == 0 + export_vtk(pdp, t, time) end - # next!(p) + next!(p) end - # finish!(p) + finish!(p) return nothing end -@timeit TO function update_velhalf!(body::AbstractPDBody, Δt½::Float64) - @inbounds @threads :static 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½ +@timeit TO function update_velhalf!(gs::GlobalStorage, Δt½::Float64) + @inbounds @threads :static for i in axes(gs.velocity_half, 2) + gs.velocity_half[1, i] = gs.velocity[1, i] + gs.acceleration[1, i] * Δt½ + gs.velocity_half[2, i] = gs.velocity[2, i] + gs.acceleration[2, i] * Δt½ + gs.velocity_half[3, i] = gs.velocity[3, i] + gs.acceleration[3, i] * Δt½ end return nothing end -@timeit TO function update_disp_and_position!(body::AbstractPDBody, Δt::Float64) - @inbounds @threads :static 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 - body.position[1, i] += body.velocity_half[1, i] * Δt - body.position[2, i] += body.velocity_half[2, i] * Δt - body.position[3, i] += body.velocity_half[3, i] * Δt +@timeit TO function update_disp_and_position!(gs::GlobalStorage, Δt::Float64) + @inbounds @threads :static for i in axes(gs.displacement, 2) + gs.displacement[1, i] += gs.velocity_half[1, i] * Δt + gs.displacement[2, i] += gs.velocity_half[2, i] * Δt + gs.displacement[3, i] += gs.velocity_half[3, i] * Δt + gs.position[1, i] += gs.velocity_half[1, i] * Δt + gs.position[2, i] += gs.velocity_half[2, i] * Δt + gs.position[3, i] += gs.velocity_half[3, i] * Δt end return nothing end -@timeit TO function compute_equation_of_motion!(body::AbstractPDBody, Δt½::Float64, mat::PDMaterial) - @inbounds @threads :static for i in 1:body.n_points - body.acceleration[1, i] = (body.b_int[1, i, 1] + body.b_ext[1, i]) / mat[i].rho - body.acceleration[2, i] = (body.b_int[2, i, 1] + body.b_ext[2, i]) / mat[i].rho - body.acceleration[3, i] = (body.b_int[3, i, 1] + body.b_ext[3, i]) / mat[i].rho - body.velocity[1, i] = body.velocity_half[1, i] + body.acceleration[1, i] * Δt½ - body.velocity[2, i] = body.velocity_half[2, i] + body.acceleration[2, i] * Δt½ - body.velocity[3, i] = body.velocity_half[3, i] + body.acceleration[3, i] * Δt½ +@timeit TO function compute_equation_of_motion!(pdp::PDProblem, Δt½::Float64) + @inbounds @threads :static for i in 1:pdp.sp.pc.n_points + pdp.gs.acceleration[1, i] = (pdp.gs.b_int[1, i, 1] + pdp.gs.b_ext[1, i]) / pdp.sp.mat[i].rho + pdp.gs.acceleration[2, i] = (pdp.gs.b_int[2, i, 1] + pdp.gs.b_ext[2, i]) / pdp.sp.mat[i].rho + pdp.gs.acceleration[3, i] = (pdp.gs.b_int[3, i, 1] + pdp.gs.b_ext[3, i]) / pdp.sp.mat[i].rho + pdp.gs.velocity[1, i] = pdp.gs.velocity_half[1, i] + pdp.gs.acceleration[1, i] * Δt½ + pdp.gs.velocity[2, i] = pdp.gs.velocity_half[2, i] + pdp.gs.acceleration[2, i] * Δt½ + pdp.gs.velocity[3, i] = pdp.gs.velocity_half[3, i] + pdp.gs.acceleration[3, i] * Δt½ end return nothing end diff --git a/test/integration_tests/b_int_bond_based.jl b/test/integration_tests/b_int_bond_based.jl index e482531b..999610d7 100644 --- a/test/integration_tests/b_int_bond_based.jl +++ b/test/integration_tests/b_int_bond_based.jl @@ -44,7 +44,7 @@ if Threads.nthreads() <= 2 0.0 0.0 ] - Peridynamics.update_thread_cache!(body) + Peridynamics.reduce_tls_to_gs!(body) @test body.b_int[:,:,1] ≈ [b¹² -b¹²] else @warn "Test omitted! Threads.nthreads() should be <= 2"