Skip to content

Commit

Permalink
PC2 refactoring to new data structure
Browse files Browse the repository at this point in the history
  • Loading branch information
kaipartmann committed Sep 6, 2023
1 parent 87b03b0 commit bcf2ad7
Show file tree
Hide file tree
Showing 13 changed files with 434 additions and 379 deletions.
6 changes: 4 additions & 2 deletions src/Peridynamics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
58 changes: 36 additions & 22 deletions src/bond_based.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
155 changes: 15 additions & 140 deletions src/conditions.jl
Original file line number Diff line number Diff line change
@@ -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
2 changes: 1 addition & 1 deletion src/dynamic_relaxation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
48 changes: 5 additions & 43 deletions src/io.jl
Original file line number Diff line number Diff line change
@@ -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"
Expand Down Expand Up @@ -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
Expand Down
18 changes: 9 additions & 9 deletions src/jobs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading

0 comments on commit bcf2ad7

Please sign in to comment.