thermo-mechanical coupling

[vavrines]

Hi @kaipartmann

We're working on the thermo-mechanical coupling model. A script from scratch is here (https://github.com/vavrines/KitPD.jl/blob/main/couple/script.jl).
We're just about to clean and refactor the codes. The idea is to reuse the structs and methods here as much as possible.
I'd like to know if you're interested in implementing this functionality. If things go well, we should be able to contribute codes, but your assistance will likely be needed.

[kaipartmann]

Hello @vavrines,
Thank you for reaching out!

I am very interested in implementing this into Peridynamics.jl.
Your code is very well structured and comfortable to read, so implementing it using your example should be feasible.
Do you have literature describing your thermo-mechanical coupling model?

Would you also be interested in us collaborating on a publication about this topic? We have good HPC access and can run large-scale simulations once the model works with the package.

[vavrines]

@kaipartmann
Thanks for the feedback!
Basically we follow Chapter 13 in the monograph Peridynamic Theory and Its Applications.
I'd love to work together with you. We'll first try to implement the thermal model. I'll let you know our process and there may be some beginner's questions in between.
Our long-term plan is to do fluid-solid interaction by coupling CFD and PD codes.

[kaipartmann]

@vavrines, this sounds like a good plan!
I recently talked with the developers of TrixiParticles.jl.
Maybe their repo is also worth looking for you.

Using your example, I hacked something together. I am unsure if it works correctly - the code is probably full of bugs/errors. The structure of Peridynamics.jl is currently not suited for coupled problems, so I had to override many functions, which resulted in a lot of complicated internal code. It also only works with the current state of the main branch and is likely to break.

I am still in the process of reworking the internal structure of the package. I will consider how a coupled model could be easily implemented without changing many internal functions for the next version.

However, you can probably try this out for your model development.
If you have questions, please do not hesitate to contact me.

Example simulation: heating of a bar

The best thing is that you start from an empty new project and add the current main branch of the package:

shell> mkdir pd_thermal_coupling

shell> cd pd_thermal_coupling

shell> julia

julia> ]

(@v1.10) pkg> activate .

(pd_thermal_coupling) pkg> add Peridynamics#main

Then you create the following file:

bb_thermal_coupling_model.jl

using Peridynamics
using Base.Threads

struct BBTMaterial <: Peridynamics.AbstractMaterial end

struct BBTPointParameters <: Peridynamics.AbstractPointParameters
    δ::Float64
    rho::Float64
    E::Float64
    nu::Float64
    G::Float64
    K::Float64
    λ::Float64
    μ::Float64
    Gc::Float64
    εc::Float64
    bc::Float64
    kc::Float64 # thermal conductivity
    kp::Float64 # microconductivity
    aph::Float64 # thermal expansion
    cv::Float64 # specific heat capacity
end

@inline Peridynamics.point_param_type(::BBTMaterial) = BBTPointParameters

@inline function Peridynamics.allowed_material_kwargs(::BBTMaterial)
    return (Peridynamics.DEFAULT_POINT_KWARGS..., :kc, :aph, :cv)
end

function Peridynamics.get_point_params(::BBTMaterial, p::Dict{Symbol,Any})
    δ = Peridynamics.get_horizon(p)
    rho = Peridynamics.get_density(p)
    if haskey(p, :nu)
        msg = "keyword `nu` is not allowed for BBTMaterial!\n"
        msg *= "Bond-based peridynamics has a limitation on the possion ratio.\n"
        msg *= "Therefore, when using BBTMaterial, `nu` is hardcoded as 1/4.\n"
        throw(ArgumentError(msg))
    else
        p[:nu] = 0.25
    end
    E, nu, G, K, λ, μ = Peridynamics.get_elastic_params(p)
    Gc, εc = Peridynamics.get_frac_params(p, δ, K)
    bc = 18 * K / (π * δ^4) # bond constant
    kc, kp, aph, cv = get_thermal_params(p, δ)
    return BBTPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc, kc, kp, aph, cv)
end

function get_thermal_params(p::Dict{Symbol,Any}, δ)
    haskey(p, :kc) || throw(UndefKeywordError(:kc))
    haskey(p, :aph) || throw(UndefKeywordError(:aph))
    haskey(p, :cv) || throw(UndefKeywordError(:cv))
    kc::Float64 = float(p[:kc])
    kc ≤ 0 && throw(ArgumentError("`kc` should be larger than zero!\n"))
    aph::Float64 = float(p[:aph])
    aph ≤ 0 && throw(ArgumentError("`aph` should be larger than zero!\n"))
    cv::Float64 = float(p[:cv])
    cv ≤ 0 && throw(ArgumentError("`cv` should be larger than zero!\n"))
    kp = 6 * kc / (π * δ^4)
    return kc, kp, aph, cv
end

@inline Peridynamics.discretization_type(::BBTMaterial) = Peridynamics.BondDiscretization

@inline function Peridynamics.init_discretization(body::Body{BBTMaterial}, args...)
    return Peridynamics.init_bond_discretization(body, args...)
end

struct BBTVerletStorage <: Peridynamics.AbstractStorage
    position::Matrix{Float64}
    displacement::Matrix{Float64}
    velocity::Matrix{Float64}
    velocity_half::Matrix{Float64}
    acceleration::Matrix{Float64}
    b_int::Matrix{Float64}
    b_ext::Matrix{Float64}
    damage::Vector{Float64}
    bond_active::Vector{Bool}
    n_active_bonds::Vector{Int}
    temperature::Vector{Float64}
    pflux::Vector{Float64}
end

const BBTStorage = Union{BBTVerletStorage}

@inline Peridynamics.storage_type(::BBTMaterial, ::VelocityVerlet) = BBTVerletStorage

function Peridynamics.init_storage(::Body{BBTMaterial}, ::VelocityVerlet,
                                   bd::Peridynamics.BondDiscretization,
                                   ch::Peridynamics.ChunkHandler)
    n_loc_points = length(ch.loc_points)
    n_points_with_halo = length(ch.point_ids)
    position = copy(bd.position)
    displacement = zeros(3, n_loc_points)
    velocity = zeros(3, n_points_with_halo)
    velocity_half = zeros(3, n_loc_points)
    acceleration = zeros(3, n_loc_points)
    b_int = zeros(3, n_loc_points)
    b_ext = zeros(3, n_loc_points)
    damage = zeros(n_loc_points)
    bond_active = ones(Bool, length(bd.bonds))
    n_active_bonds = copy(bd.n_neighbors)
    temperature = zeros(n_points_with_halo)
    pflux = zeros(n_loc_points)
    return BBTVerletStorage(position, displacement, velocity, velocity_half, acceleration,
                            b_int, b_ext, damage, bond_active, n_active_bonds, temperature,
                            pflux)
end

Peridynamics.storage_field(s::BBTStorage, ::Val{:temperature}) = getfield(s, :temperature)
Peridynamics.storage_field(s::BBTStorage, ::Val{:pflux}) = getfield(s, :pflux)

@inline function Peridynamics.get_halo_read_fields(s::BBTStorage)
    return (s.position, s.velocity, s.temperature)
end
@inline Peridynamics.get_halo_write_fields(::BBTStorage) = ()

const BBTBodyChunk = Union{Peridynamics.BodyChunk{BBTMaterial,P,D,S},
                           Peridynamics.MultiParamBodyChunk{BBTMaterial,P,D,S}} where {P,D,S}

function calc_pflux!(b::BBTBodyChunk)
    b.store.pflux .= 0.0
    for point_id in Peridynamics.each_point_idx(b)
        pflux_point!(b.store, b.dscr, b.mat, get_param(b, point_id), point_id)
    end
    return nothing
end

function pflux_point!(s::BBTStorage, bd::Peridynamics.BondDiscretization, ::BBTMaterial,
                      param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length
        Δxijx = s.position[1, j] - s.position[1, i]
        Δxijy = s.position[2, j] - s.position[2, i]
        Δxijz = s.position[3, j] - s.position[3, i]
        Δvijx = s.velocity[1, j] - s.velocity[1, i]
        Δvijy = s.velocity[2, j] - s.velocity[2, i]
        Δvijz = s.velocity[3, j] - s.velocity[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)
        ev = (Δxijx * Δvijx + Δxijy * Δvijy + Δxijz * Δvijz) / l
        Δtem = s.temperature[j] - s.temperature[i]
        s.pflux[i] += s.bond_active[bond_id] * (param.kp * Δtem / L -
                      ev * 0.5 * param.bc * param.aph) * bd.volume[j]
    end
    return nothing
end

@inline function update_temperature!(b::Peridynamics.AbstractBodyChunk, Δt::Float64)
    for point_id in Peridynamics.each_point_idx(b)
        param = get_param(b, point_id)
        k = Δt / (param.rho * param.cv)
        _update_temperature!(b.store.temperature, b.store.pflux, k, point_id)
    end
    return nothing
end

@inline function _update_temperature!(temperature, pflux, k, i)
    temperature[i] += pflux[i] * k
    return nothing
end

function Peridynamics.force_density_point!(s::BBTStorage,
                                           bd::Peridynamics.BondDiscretization,
                                           ::BBTMaterial, param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length

        # current bond length
        Δxijx = s.position[1, j] - s.position[1, i]
        Δxijy = s.position[2, j] - s.position[2, i]
        Δxijz = s.position[3, j] - s.position[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)

        # bond strain
        εm = (l - L) / L
        εt = 0.5 * (s.temperature[j] - s.temperature[i]) * param.aph
        ε = εm + εt

        # failure mechanism
        if ε > param.εc && bond.fail_permit
            s.bond_active[bond_id] = false
        end
        s.n_active_bonds[i] += s.bond_active[bond_id]

        # update of force density
        temp = s.bond_active[bond_id] * param.bc * ε / l * bd.volume[j]
        s.b_int[1, i] += temp * Δxijx
        s.b_int[2, i] += temp * Δxijy
        s.b_int[3, i] += temp * Δxijz
    end

    return nothing
end

# You have to override the time stepping schemes
function Peridynamics.solve_timestep!(dh::Peridynamics.ThreadsDataHandler{BBTBodyChunk},
                                      options::Peridynamics.ExportOptions, Δt::Float64,
                                      Δt½::Float64, n::Int)
    t = n * Δt
    @threads :static for chunk_id in eachindex(dh.chunks)
        chunk = dh.chunks[chunk_id]
        Peridynamics.update_vel_half!(chunk, Δt½)
        Peridynamics.apply_bcs!(chunk, t)
        Peridynamics.update_disp_and_pos!(chunk, Δt)
    end
    @threads :static for chunk_id in eachindex(dh.chunks)
        Peridynamics.exchange_read_fields!(dh, chunk_id)
        chunk = dh.chunks[chunk_id]
        calc_pflux!(chunk)
        update_temperature!(chunk, Δt)
    end
    @threads :static for chunk_id in eachindex(dh.chunks)
        Peridynamics.exchange_read_fields!(dh, chunk_id)
        Peridynamics.calc_force_density!(dh.chunks[chunk_id])
    end
    @threads :static for chunk_id in eachindex(dh.chunks)
        exchange_write_fields!(dh, chunk_id)
        chunk = dh.chunks[chunk_id]
        calc_damage!(chunk)
        update_acc_and_vel!(chunk, Δt½)
        export_results(dh, options, chunk_id, n, t)
    end
    return nothing
end

function Peridynamics.solve_timestep!(dh::Peridynamics.MPIDataHandler{BBTBodyChunk},
                                      options::Peridynamics.ExportOptions, Δt::Float64,
                                      Δt½::Float64, n::Int)
    t = n * Δt
    chunk = dh.chunk
    Peridynamics.@timeit_debug Peridynamics.TO "update_vel_half!" Peridynamics.update_vel_half!(chunk, Δt½)
    Peridynamics.@timeit_debug Peridynamics.TO "apply_bcs!" Peridynamics.apply_bcs!(chunk, t)
    Peridynamics.@timeit_debug Peridynamics.TO "update_disp_and_pos!" Peridynamics.update_disp_and_pos!(chunk, Δt)
    Peridynamics.@timeit_debug Peridynamics.TO "exchange_read_fields!" Peridynamics.exchange_read_fields!(dh)
    Peridynamics.@timeit_debug Peridynamics.TO "calc_pflux!" calc_pflux!(chunk)
    Peridynamics.@timeit_debug Peridynamics.TO "update_temperature!" update_temperature!(chunk, Δt)
    Peridynamics.@timeit_debug Peridynamics.TO "exchange_read_fields!" Peridynamics.exchange_read_fields!(dh)
    Peridynamics.@timeit_debug Peridynamics.TO "calc_force_density!" Peridynamics.calc_force_density!(chunk)
    Peridynamics.@timeit_debug Peridynamics.TO "exchange_write_fields!" Peridynamics.exchange_write_fields!(dh)
    Peridynamics.@timeit_debug Peridynamics.TO "calc_damage!" Peridynamics.calc_damage!(chunk)
    Peridynamics.@timeit_debug Peridynamics.TO "update_acc_and_vel!" Peridynamics.update_acc_and_vel!(chunk, Δt½)
    Peridynamics.@timeit_debug Peridynamics.TO "export_results" Peridynamics.export_results(dh, options, n, t)
    return nothing
end

function temperature_bc!(f::F, b::Body, name::Symbol) where {F<:Function}
    Peridynamics.check_if_set_is_defined(b.point_sets, name)
    type = Peridynamics.check_condition_function(f)
    if type === :sdbc
        sdbc = Peridynamics.SingleDimBC(f, :temperature, name, 0x1)
        Peridynamics._condition!(b.single_dim_bcs, sdbc)
    elseif type === :pdsdbc
        pdsdbc = Peridynamics.PosDepSingleDimBC(f, :temperature, name, 0x1)
        Peridynamics._condition!(b.posdep_single_dim_bcs, pdsdbc)
    end
    return nothing
end

function temperature_ic!(b::Body, name::Symbol, value::Real)
    Peridynamics.check_if_set_is_defined(b.point_sets, name)
    sdic = Peridynamics.SingleDimIC(convert(Float64, value), :temperature, name, 0x1)
    Peridynamics._condition!(b.single_dim_ics, sdic)
    return nothing
end

@inline function Peridynamics.apply_sdbc!(field::Vector{Float64}, value::Float64, ::UInt8,
                                          point_ids::Vector{Int})
    @simd for i in point_ids
        @inbounds field[i] = value
    end
    return nothing
end

@inline function Peridynamics.apply_pdsdbc!(field::Vector{Float64}, pos::Matrix{Float64},
                                            bc::Peridynamics.PosDepSingleDimBC{F},
                                            point_ids::Vector{Int}, t::Float64) where {F}
    @simd for i in point_ids
        value = bc(SVector{3}(pos[1,i], pos[2,i], pos[3,i]), t)
        if !isnan(value)
            field[i] = value
        end
    end
    return nothing
end

function Peridynamics.apply_ic!(b::BBTBodyChunk, ic::Peridynamics.SingleDimIC)
    for point_id in b.psets[ic.point_set]
        _apply_ic!(Peridynamics.get_storage_field(b.store, ic.field), ic.value, ic.dim,
                   point_id)
    end
    return nothing
end

function _apply_ic!(field::Matrix{Float64}, value::Float64, dim::UInt8, point_id::Int)
    setindex!(field, value, dim, point_id)
    return nothing
end

function _apply_ic!(field::Vector{Float64}, value::Float64, ::UInt8, point_id::Int)
    setindex!(field, value, point_id)
    return nothing
end

Then you can run a simulation with the new API of the package:

include("bb_thermal_coupling_model.jl")

l, Δx, δ = 1.0, 1/100, 3.015/100
pos, vol = uniform_box(l, 0.1l, 0.1l, Δx)
body = Body(BBTMaterial(), pos, vol)
# add material parameters as in your code
material!(body; horizon=3.015Δx, E=2.1e5, rho=8e-6, Gc=2.7, kc=51.9, aph=11.5e-6, cv=472)
point_set!(x -> x > 0.0, body, :heating_area)
point_set!(x -> x < -l/2+3Δx, body, :fixed_area)
# temperature boundary condition: increase temperature by 1e5 K/s
temperature_bc!(t -> 1e5 * t, body, :heating_area)
velocity_bc!(t -> 0.0, body, :fixed_area, :x)
vv = VelocityVerlet(steps=2000)
job = Job(body, vv; path=joinpath(@__DIR__, "thermal_coupling"))
@time submit(job);

How it works

I mainly defined a bond-based thermal material and added new internal methods for this material type.

struct BBTMaterial <: Peridynamics.AbstractMaterial end

The pflux (thermal loop) is calculated in the pflux_point! function:

function pflux_point!(s::BBTStorage, bd::Peridynamics.BondDiscretization, ::BBTMaterial,
                      param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length
        Δxijx = s.position[1, j] - s.position[1, i]
        Δxijy = s.position[2, j] - s.position[2, i]
        Δxijz = s.position[3, j] - s.position[3, i]
        Δvijx = s.velocity[1, j] - s.velocity[1, i]
        Δvijy = s.velocity[2, j] - s.velocity[2, i]
        Δvijz = s.velocity[3, j] - s.velocity[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)
        ev = (Δxijx * Δvijx + Δxijy * Δvijy + Δxijz * Δvijz) / l
        Δtem = s.temperature[j] - s.temperature[i]
        s.pflux[i] += s.bond_active[bond_id] * (param.kp * Δtem / L -
                      ev * 0.5 * param.bc * param.aph) * bd.volume[j]
    end
    return nothing
end

Then, the force density calculation of the plain BBMaterial is extended to account for temperature changes:

function Peridynamics.force_density_point!(s::BBTStorage,
                                           bd::Peridynamics.BondDiscretization,
                                           ::BBTMaterial, param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length

        # current bond length
        Δxijx = s.position[1, j] - s.position[1, i]
        Δxijy = s.position[2, j] - s.position[2, i]
        Δxijz = s.position[3, j] - s.position[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)

        # bond strain
        εm = (l - L) / L
        εt = 0.5 * (s.temperature[j] - s.temperature[i]) * param.aph
        ε = εm + εt

        # failure mechanism
        if ε > param.εc && bond.fail_permit
            s.bond_active[bond_id] = false
        end
        s.n_active_bonds[i] += s.bond_active[bond_id]

        # update of force density
        temp = s.bond_active[bond_id] * param.bc * ε / l * bd.volume[j]
        s.b_int[1, i] += temp * Δxijx
        s.b_int[2, i] += temp * Δxijy
        s.b_int[3, i] += temp * Δxijz
    end

    return nothing
end

The other methods are mostly internal functions.

[kaipartmann]

@vavrines I will introduce changes to the main branch in a few minutes that will break the provided code.
Therefore, I created a branch thermal_coupling that works as described with the example.

So in the local environment just do

(pd_thermal_coupling) pkg> rm Peridynamics

(pd_thermal_coupling) pkg> add Peridynamics#thermal_coupling

[ShiWeiHuH]

Thank you for your assistance. I am a member of vavrines's team and am currently working on this project. I am benefiting greatly from reviewing your code.

[ShiWeiHuH]

@vavrines I will introduce changes to the main branch in a few minutes that will break the provided code. Therefore, I created a branch thermal_coupling that works as described with the example.

So in the local environment just do

(pd_thermal_coupling) pkg> rm Peridynamics

(pd_thermal_coupling) pkg> add Peridynamics#thermal_coupling

My apologies for the interruption, may I ask you a few questions?

1. I encountered a bug when trying out the “bb_thermal_coupling_model.jl” script. It prompted an UndefVarError: discretization_type not defined at line 63 of bb_thermal_coupling_model.jl.
2. I would like to know how the surface and volume corrections are implemented in the code, specifically in which function.
3. During the thermal-structural coupling calculation, due to the presence of virtual layers, the domains for force and heat calculations should be different. Previously, we found that if this is not taken into account, there are significant calculation deviations, especially near the boundaries. In other words, during heat calculation, the displacement virtual layer needs to be excluded, and during force calculation, the temperature virtual layer needs to be excluded. Which specific function(s) need(s) to be modified for this setting?
   Thank you again for your assistance.
   

[kaipartmann]

Hello @ShiWeiHuH! I am happy that the example benefits you! :)

1. I encountered a bug when trying out the “bb_thermal_coupling_model.jl” script. It prompted an UndefVarError: discretization_type not defined at line 63 of bb_thermal_coupling_model.jl.

It seems that you're encountering an issue with the branch on GitHub. Please try in the package mode

pkg> status
[4dc47793] Peridynamics v0.2.0 `https://github.com/kaipartmann/Peridynamics.jl.git#thermal_coupling`

and check if it reads thermal_coupling at the end of the Peridynamics.jl entry.

2. I would like to know how the surface and volume corrections are implemented in the code, specifically in which function.

We currently do not have volume or surface corrections implemented. If you need that, I can add it. I tried to implement it before but did not finish to debug the code.

3. During the thermal-structural coupling calculation, due to the presence of virtual layers, the domains for force and heat calculations should be different. Previously, we found that if this is not taken into account, there are significant calculation deviations, especially near the boundaries. In other words, during heat calculation, the displacement virtual layer needs to be excluded, and during force calculation, the temperature virtual layer needs to be excluded. Which specific function(s) need(s) to be modified for this setting?

This is very interesting, and I did not know that. I am very new to thermal models with peridynamics, so I am unsure if I correctly understand the actual consequences of this regarding the code. So, is it a problem if a point has some mechanical and thermal boundary conditions simultaneously? Can you highlight your implementation of these effects in the script you provided for me?

[ShiWeiHuH]

Hi, @kaipartmann.Thank you very much for your response.
As for question 1: I'm very sorry, the bug was caused by my own operational mistake. It has been resolved now. Thank you for your guidance.
As for question 2: With a sufficient number of PD points, surface corrections and volume corrections can also yield good calculation results. I'm also trying to make some modifications to enhance my understanding of your code. Your code is very inspiring, and I'm learning from it.
3. Regarding the issue of thermal-mechanical coupling, I primarily addressed it in the init.jl file, where I also implemented surface corrections and volume corrections. I'll mark the relevant parts for you.
Once again, thank you for your assistance.

[ShiWeiHuH]

Hi, @kaipartmann. We have updated the init.jl file.
In the new init.jl, compared to the original version, the main content remains unchanged, but detailed comments have been added. We have annotated the codes for handling the thermal domain and the mechanical domain separately, while also marking the codes for surface correction and volume correction.
During thermal coupling calculations, it's not always the case that the mechanical domain and the thermal domain are different. If there are displacement or velocity boundary conditions, as well as temperature boundary conditions, it is necessary to add a fictitious layer. Only when the fictitious layers for thermal and deformation solutions are located at different positions, should the scenario of different domains for motion and heat transfer calculations be considered.
There are also cases where both temperature and displacement or velocity boundary conditions are set in the same fictitious layer, such as in thermal pulling. On the side receiving the pulling force, it's possible to simultaneously set temperature and displacement boundary conditions. However, for the fixed end on the other side, there might not be temperature boundary conditions, leading to a scenario where the thermal and mechanical domains are different.

[kaipartmann]

Hi @ShiWeiHuH,
Thank you for your updates in the init.jl file. I think I can now better understand what the fictitious layers do. I hope to get some time in the next few days to implement these into an updated version of the bb_thermal_coupling_model.jl script. I am unsure how long it will take, but I will ping you if I encounter problems or have any questions.

[ShiWeiHuH]

Hi @ShiWeiHuH, Thank you for your updates in the init.jl file. I think I can now better understand what the fictitious layers do. I hope to get some time in the next few days to implement these into an updated version of the bb_thermal_coupling_model.jl script. I am unsure how long it will take, but I will ping you if I encounter problems or have any questions.
Hello, @kaipartmann, I want to express my gratitude once again for your assistance. Over this period, I've been working on implementing surface and volume corrections based on your provided code, particularly focusing on the 'bond_discretization.jl' section. As a result, I've developed a 'modify.jl' file dedicated to this task. The basic functionality has been implemented. By executing the following three lines of code, we can derive the necessary correction coefficients.
Here is the modified code :

include("bb_thermal_coupling_model.jl")
include("modify.jl")
Peridynamics.get_git_info() = "*"
l, Δx, δ = 1.0, 1/100, 3.015/100
pos, vol = uniform_box(l+6*Δx, 0.1l, 0.1l, Δx)
body = Body(BBTMaterial(), pos, vol)
# add material parameters as in your code
material!(body; horizon=3.015*Δx, E=2.1e5, rho=8e-6, Gc=2.7, kc=51.9, aph=11.5e-6, cv=472)
point_set!(x -> x < -l/2, body, :tem_bc)
point_set!(x -> x >  l/2, body, :fix_bc)
# temperature boundary condition: increase temperature by 1e5 K/s
temperature_bc!(t -> 1e5 * t, body, :tem_bc)
velocity_bc!(t -> 0.0, body, :fix_bc, :x)
mp = mpm_init(Dict(:dx => Δx, :horizon => 3.015*Δx, :emod => 2.1e5, :th_con => 51.9)) ### creat parameters for modify
sbd = mbd_init(body, mp) ### creat the index and volume modify coefficient
mm = undate_mbd(sbd, body, mp) ### creat surface coefficient
vv = VelocityVerlet(steps=2000)
job = Job(body, vv; path=joinpath(@__DIR__, "thermal_coupling"), fields=(:displacement, :damage))
@time submit(job)

This part is newly added：

mp = mpm_init(Dict(:dx => Δx, :horizon => 3.015*Δx, :emod => 2.1e5, :th_con => 51.9)) ### creat parameters for modify
sbd = mbd_init(body, mp) ### creat the index and volume modify coefficient
mm = undate_mbd(sbd, body, mp) ### creat surface coefficient

The modify.jl

using Peridynamics
using NearestNeighbors

mutable struct Mbond
    neighbor::Int
    length::Float64
    vol_m::Float64    ### voleme modify coefficient 
    surf_th::Float64  ### surface modify coefficient of thermal zone
    surf_mh::Float64  ### surface modify coefficient of deformation zone 
end

struct MbondDiscretization <: Peridynamics.AbstractDiscretization
    bonds::Vector{Mbond}
    n_neighbors::Vector{Int}
    bond_ids::Vector{UnitRange{Int}}
    volume::Vector{Float64}
end

struct Mparameters 
dx::Float64
δ::Float64
kc::Float64 # thermal conductivity
kp::Float64 # microconductivity
E::Float64
bc::Float64
end

function mpm_init(p::Dict{Symbol,Float64})
    dx = p[:dx]
    δ = p[:horizon]
    kc = p[:th_con]
    kp = 6 * kc / (π * δ^4)
    E = p[:emod]
    bc = 12 * E / (π * δ^4)
    mpm = Mparameters(dx, δ, kc, kp, E, bc)
    return mpm
end

function mbd_init(body::Body, mpm::Mparameters)
    bonds, n_neighbors = mfind_bonds(body, mpm)
    bond_ids = mfind_bond_ids(n_neighbors)
    volume = body.volume
    mbd = MbondDiscretization(bonds, n_neighbors, bond_ids,volume)
    return mbd
end

function mfind_bonds(body::Body, mpm::Mparameters)
    balltree = BallTree(body.position)
    bonds = Vector{Mbond}()
    sizehint!(bonds, body.n_points * 100)
    n_neighbors = zeros(Int, body.n_points)
    for i in 1:body.n_points
        n_neighbors[i] = mfind_bonds!(bonds, balltree, body.position, mpm, i)
    end
    return bonds, n_neighbors
end

function mfind_bonds!(bonds::Vector{Mbond}, balltree::BallTree, position::Matrix{Float64}, mpm::Mparameters, i::Int)

    neigh_idxs_with_i = inrange(balltree, view(position, :, i), mpm.δ)
    neigh_idxs = filter(x -> x != i, neigh_idxs_with_i)
    for j in neigh_idxs
        L = sqrt((position[1, j] - position[1, i])^2 +
                 (position[2, j] - position[2, i])^2 +
                 (position[3, j] - position[3, i])^2)

                 if L < (mpm.δ - 0.5 * mpm.dx) 
                    vol_m = 1.0
                elseif L < (mpm.δ + 0.5 * mpm.dx)  
                    vol_m  = (mpm.δ + 0.5 * mpm.dx - L) / mpm.dx
                else
                    vol_m = 0.0
                end
                 
                surf_th = 1.0
                surf_mh = 1.0
        push!(bonds, Mbond(j, L, vol_m, surf_th, surf_mh))
    end
    return length(neigh_idxs)
end

function mfind_bond_ids(n_neighbors::Vector{Int})
    bond_ids = fill(0:0, length(n_neighbors))
    bonds_start, bonds_end = 1, 0
    for i in eachindex(n_neighbors)
        bonds_end = bonds_start + n_neighbors[i] - 1
        bond_ids[i] = bonds_start:bonds_end
        bonds_start += n_neighbors[i]
    end
    return bond_ids
end

function undate_mbd(mbd::MbondDiscretization, bo::Body, mpm::Mparameters)

    tem = 0.001.*[sum(bo.position[:, i]) for i in 1:bo.n_points]
    disx = copy(bo.position)
    disx[1, :] .= 1.001 .* disx[1, :]
    disy = copy(bo.position)
    disy[2, :] .= 1.001 .* disy[2, :]
    disz = copy(bo.position)
    disz[3, :] .= 1.001 .* disz[3, :]

    point_th = zeros(Float64, bo.n_points)
    point_mx = zeros(Float64, bo.n_points)
    point_my = zeros(Float64, bo.n_points)
    point_mz = zeros(Float64, bo.n_points)

    for i in 1:bo.n_points
        pd_th = 0
        pd_mechx = 0
        pd_mechy = 0
        pd_mechz = 0

        for bond_id in meach_bond_idx(mbd, i)
            bond = mbd.bonds[bond_id]
            j, L = bond.neighbor, bond.length

            if i in bo.point_sets[:fix_bc] || j in bo.point_sets[:fix_bc]
                mfth = 0
            else
                mfth = 1
            end

            if i in bo.point_sets[:tem_bc] || j in bo.point_sets[:tem_bc]
                mfmech = 0
            else
                mfmech = 1
            end

            nlengthx = sqrt((disx[1,j] - disx[1,i])^2 + (disx[2,j] - disx[2,i])^2 + (disx[3,j] - disx[3,i])^2)
            nlengthy = sqrt((disy[1,j] - disy[1,i])^2 + (disy[2,j] - disy[2,i])^2 + (disy[3,j] - disy[3,i])^2)
            nlengthz = sqrt((disz[1,j] - disz[1,i])^2 + (disz[2,j] - disz[2,i])^2 + (disz[3,j] - disz[3,i])^2)

            pd_th += 0.25 * mpm.kp * (tem[j] - tem[i])^2 / L * mbd.volume[j] * bond.vol_m * mfth
            pd_mechx += 0.25 * mpm.bc * ((nlengthx-L) / L)^2 * mbd.volume[j] * bond.vol_m * mfmech
            pd_mechy += 0.25 * mpm.bc * ((nlengthy-L) / L)^2 * mbd.volume[j] * bond.vol_m * mfmech
            pd_mechz += 0.25 * mpm.bc * ((nlengthz-L) / L)^2 * mbd.volume[j] * bond.vol_m * mfmech
        end 

        if pd_th != 0
            point_th[i,1] = 1e-6 *  mpm.kc * 1.5 / pd_th
        else
            point_th[i,1] = 0
        end

        if pd_mechx != 0
            point_mx[i,1] = 0.6 * 1e-6 * mpm.E / pd_mechx
        else
            point_mx[i,1] = 0
        end
        
        if pd_mechy != 0
            point_my[i,1] = 0.6 * 1e-6 * mpm.E / pd_mechy
        else
            point_my[i,1] = 0
        end
        
        if pd_mechz != 0
            point_mz[i,1] = 0.6 * 1e-6 * mpm.E / pd_mechz
        else
            point_mz[i,1] = 0
        end
    end

    for i in 1:bo.n_points
        for bond_id in meach_bond_idx(mbd, i)
            bond = mbd.bonds[bond_id]
            j, L = bond.neighbor, bond.length

            bond.surf_th = 0.5 * (point_th[i]+point_th[j])

            if abs(bo.position[3,j] - bo.position[3,i]) < 1e-10
                if abs(bo.position[2,j] - bo.position[2,i]) < 1e-10
                    θ = 0
                elseif abs(bo.position[1,j] - bo.position[1,i]) < 1e-10
                    θ = π/2
                else
                    θ = atan(abs((bo.position[2,j] - bo.position[2,i])/(bo.position[1,j] - bo.position[1,i])))
                end
                ϕ = π/2
                sx = (point_mx[i] + point_mx[j]) / 2
                sy = (point_my[i] + point_my[j]) / 2
                sz = (point_mz[i] + point_mz[j]) / 2
                bond.surf_mh = sqrt(1/((cos(θ)*sin(ϕ)/sx)^2 + (sin(θ)*sin(ϕ)/sy)^2 + (cos(ϕ)/sz)^2))
            elseif (abs(bo.position[2,j] - bo.position[2,i]) < 1e-10) && (abs(bo.position[1,j] - bo.position[1,i]) < 1e-10)
                sz = (point_mz[i] + point_mz[j]) / 2
                bond.surf_mh = sz
            else
                θ = atan(abs((bo.position[2,j] - bo.position[2,i])/(bo.position[1,j] - bo.position[1,i])))
                ϕ = acos((bo.position[3,j] - bo.position[3,i])/L)
                sx = (point_mx[i] + point_mx[j]) / 2
                sy = (point_my[i] + point_my[j]) / 2
                sz = (point_mz[i] + point_mz[j]) / 2
                bond.surf_mh = sqrt(1/((cos(θ)*sin(ϕ)/sx)^2 + (sin(θ)*sin(ϕ)/sy)^2 + (cos(ϕ)/sz)^2))
            end
        end  
    end  
    return mbd        
end
    
@inline meach_bond_idx(mbd::MbondDiscretization, point_id::Int) = mbd.bond_ids[point_id]

My next step is to pass the correction coefficients into 'bb_thermal_coupling_model.jl'.
Additionally, you mentioned that the time stepping schemes need to be rewritten. What are the main points to pay attention to in this part? I will try to proceed with it next.

[kaipartmann]

Hi @ShiWeiHuH,

Thank you very much for providing your code!

First, I'm sorry about the bug with the get_git_info() function. It has already been solved in the main branch.

I ran your code, and it should work fine. You are correct in bridging the parallel computation features with your newly defined MBondDiscretization, because they cause problems when calculating the surface correction factors.

I encountered the same issue and am working on it (#73). I'm now prioritizing the implementation of the surface correction (see #74), and your assistance is greatly appreciated.

I already tried to do something similar to what you did, but my code still has some bugs. You can take a look at the surface_correction branch in the file:
https://github.com/kaipartmann/Peridynamics.jl/blob/surface_correction/src/physics/bond_based_corrected.jl
These lines look very similar to your code:

Peridynamics.jl/src/physics/bond_based_corrected.jl

Lines 205 to 243 in 1cb4317

	function calc_scr(Δxijx, Δxijy, Δxijz, L, h, i, j)
	if abs(Δxijz) <= 1e-10
	if abs(Δxijy) <= 1e-10
	θ = 0.0
	elseif abs(Δxijx) <= 1e-10
	θ = 90 * π / 180
	else
	θ = atan(abs(Δxijy)/abs(Δxijx))
	end
	ϕ = 90 * π / 180
	scx = (h[1,i] + h[1,j]) / 2
	scy = (h[2,i] + h[2,j]) / 2
	scz = (h[3,i] + h[3,j]) / 2
	scr = sqrt(
	1/(
	(cos(θ) * sin(ϕ))^2 / scx^2 +
	(sin(θ) * sin(ϕ))^2 / scy^2 +
	cos(ϕ)^2 / scz^2
	)
	)
	elseif abs(Δxijx) <= 1e-10 && abs(Δxijy) <= 1e-10
	scz = (h[3,i] + h[3,j]) / 2
	scr = scz
	else
	θ = atan(abs(Δxijy)/abs(Δxijx))
	ϕ = acos(abs(Δxijz)/L)
	scx = (h[1,i] + h[1,j]) / 2
	scy = (h[2,i] + h[2,j]) / 2
	scz = (h[3,i] + h[3,j]) / 2
	scr = sqrt(
	1/(
	(cos(θ) * sin(ϕ))^2 / scx^2 +
	(sin(θ) * sin(ϕ))^2 / scy^2 +
	cos(ϕ)^2 / scz^2
	)
	)
	end
	return scr
	end

You can proceed by adding your MBondDiscretization as a discretization method with

@inline function Peridynamics.init_discretization(body::Body{BBTMaterial}, args...)
    return mbd_init(body, ...)
end

But your mbd_init function needs to return a ChunkHandler and has to be similar to:

Peridynamics.jl/src/discretizations/bond_discretization.jl

Lines 15 to 28 in a71541f

	function init_bond_discretization(body::Body, pd::PointDecomposition, chunk_id::Int)
	check_bond_discretization_compat(body.mat)
	bonds, n_neighbors = find_bonds(body, pd.decomp[chunk_id])
	ch = ChunkHandler(bonds, pd, chunk_id)
	localize!(bonds, ch.localizer)
	bond_ids = find_bond_ids(n_neighbors)
	position, volume = get_pos_and_vol_chunk(body, ch.point_ids)
	bd = BondDiscretization(position, volume, bonds, n_neighbors, bond_ids)
	return bd, ch
	end

However, please note that in the main branch, this process is already more straightforward when using the @storage macro. It might be easier to proceed by copying some of the code from there.

If I am ready with #73, I can provide a new script of the current state that includes your modify.jl with the thermal computation.

[vavrines]

Hi @kaipartmann ,

Thanks for the tips!
Besides of #73 and #74, please let us know when the main branch reaches a relatively stable time.
It might not be a bad idea to start pulling things before we go too deep in the obsolete codes.

[ShiWeiHuH]

@kaipartmann
Thank you for your reply. After reviewing this portion of surface_correction, It effectively modifies surface effects. 💯 🥇 👍
This approach reminds me of the method used in the 《Peridynamic Theory and Its Applications》. Volume and surface corrections are computed as coefficients for bond force calculation in each iteration. Even though the code in this book calculates correction coefficients with each iteration when computing forces, the calculations still rely on the coordinates of the initial PD points. Correction coefficients are a physical property that reflects the initial geometric characteristics, and once the geometric structure is generated, these coefficients remain fixed and do not change. Some subsequent articles have suggested that parameters such as correction coefficients should be precomputed, and then directly called in subsequent iterations, which can significantly reduce computational overhead. Inspired by the part of your code related to bond indexing and length handling, I attempted to precompute correction coefficients and then call them in each iteration.
Regarding the issue of different thermal dynamics in different regions, it can be addressed during surface correction. For example, in thermal dynamics computations, if any PD point of a bond lies within the virtual layer of displacement constraints, the correction coefficient is set to 0, thus eliminating the influence of the virtual layer. I will proceed with experimentation as you suggested.
Thanks again for your assistance. 👍

[kaipartmann]

@ShiWeiHuH thank you for taking the time to review the code. Your assistance is greatly appreciated!

I wanted to do what you suggested but needed #73 solved for that.
A halo exchange in the initialization phase is needed to calculate the bond correction factors, and I previously designed the exchange functions so that they would only work during the time integration. Therefore, I took this quick and dirty approach to calculate the factors during each time step because then it would be possible to exchange the values to the halo points.

I will try to add the surface correction as part of the BondSystem. If I am ready, you can probably copy the code in bond_system.jl and create a ThermalBondSystem that handles thermal correction factors.

Unfortunately, the surface correction still has some bugs; I will explain this closer in #74 and ping you there.

[ShiWeiHuH]

@kaipartmann
Thank you. I understand what you mean now. My current approach is somewhat makeshift. Considering that the calculation of correction factors only occurs once at the beginning, I bypassed the issue of parallelization and simply calculated the correction factors separately and left them as is. It would be great if we could handle them uniformly within your parallel framework.
Today, I managed to integrate the computed correction factors into the coupling program via function calls. Similarly, it's not the most elegant solution, haha. Below is the code I modified. Take a look when you have time. Lastly, looking forward to your new code. Many thanks, as always.
I made just minor adjustments to bb_thermal_coupling_model.jl.
I just added a function call in the heat flux and force calculations to incorporate the correction factors calculated in th.jl and modify.jl into the computations. (th.jl and modify.jl are both the same as yesterday, no changes.）
thermal flux

function pflux_point!(s::BBTStorage, bd::Peridynamics.BondDiscretization, ::BBTMaterial,
                      param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length
     
############## introduced datas into the calculation #####newly added
        mbond = get_init_mbd(mm, bond_id) 
        volm = mbond.vol_m
        mof_th = mbond.surf_th
##############

        Δxijx = s.position[1, j] - s.position[1, i]
        Δxijy = s.position[2, j] - s.position[2, i]
        Δxijz = s.position[3, j] - s.position[3, i]
        Δvijx = s.velocity[1, j] - s.velocity[1, i]
        Δvijy = s.velocity[2, j] - s.velocity[2, i]
        Δvijz = s.velocity[3, j] - s.velocity[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)
        ev = (Δxijx * Δvijx + Δxijy * Δvijy + Δxijz * Δvijz) / l
        Δtem = s.temperature[j] - s.temperature[i]
        s.pflux[i] += s.bond_active[bond_id] * (param.kp * Δtem / L -
                      ev * 0.5 * param.bc * param.aph) * bd.volume[j] * volm * mof_th   ### modifys coefficient
    end
    return nothing
end

force density

function Peridynamics.force_density_point!(s::BBTStorage,
                                           bd::Peridynamics.BondDiscretization,
                                           ::BBTMaterial, param::BBTPointParameters, i::Int)
    for bond_id in Peridynamics.each_bond_idx(bd, i)
        bond = bd.bonds[bond_id]
        j, L = bond.neighbor, bond.length

######### introduced datas into the calculation#####newly added
        mbond = get_init_mbd(mm, bond_id)
        volm = mbond.vol_m
        mof_mh = mbond.surf_mh
#########

        # current bond length
        Δxijx = s.position[1, j] - s.position[1, i]
        Δxijy = s.position[2, j] - s.position[2, i]
        Δxijz = s.position[3, j] - s.position[3, i]
        l = sqrt(Δxijx * Δxijx + Δxijy * Δxijy + Δxijz * Δxijz)

        # bond strain
        εm = (l - L) / L
        εt = 0.5 * (s.temperature[j] - s.temperature[i]) * param.aph
        ε = εm + εt

        # failure mechanism
        if ε > param.εc && bond.fail_permit
            s.bond_active[bond_id] = false
        end
        s.n_active_bonds[i] += s.bond_active[bond_id]

        # update of force density
        temp = s.bond_active[bond_id] * param.bc * ε / l * bd.volume[j] * volm * mof_mh  ### modifys coefficient
        s.b_int[1, i] += temp * Δxijx
        s.b_int[2, i] += temp * Δxijy
        s.b_int[3, i] += temp * Δxijz
    end

    return nothing
end

Function

#####newly added
@inline get_init_mbd(modifys_coefficients::MbondDiscretization, point_id::Int64) = modifys_coefficients.bonds[point_id]

[kaipartmann]

@ShiWeiHuH @vavrines
I have a question regarding your surface correction factors.
In your code, you used:

0.6 * 1e-6 * mpm.E

How did you derive these analytical strain energy density?
The surface correction works correctly with this value, but I would like to better understand the assumptions behind it.

[ShiWeiHuH]

Hi, @kaipartmann
You can take a look at this simple derivation.
In addition, I would like to ask you what is the command to execute multi-threading or multi-node parallelism.
Also, I see that you have corrected the surface effect, and it seems that the volume has not been corrected. After all, the volume at the boundary of the near-field domain is not necessarily complete.

[kaipartmann]

Hi @ShiWeiHuH, Thank you very much for your detailed description! I added an issue (#94) because the OSBMaterial allows for the use of arbitrary Poisson's ratios.

On the parallelism:

The parallelism functionality is implemented at a very high level. With the PointDecomposition, the Body is chopped into multiple BodyChunks, where each chunk contains all relevant information for the simulation. Depending on the type of parallelism, either a ThreadsDataHandler or a MPIDataHandler is created that contain all the body chunks and all information for the simulation.

Then, everything is done on the body chunk level. For multithreading, each thread works on its chunk:

Peridynamics.jl/src/time_solvers/velocity_verlet.jl

Lines 212 to 233 in 9bb65ae

	function solve_timestep!(dh::AbstractThreadsDataHandler, options::AbstractOptions,
	Δt::Float64, Δt½::Float64, n::Int)
	t = n * Δt
	@threads :static for chunk_id in eachindex(dh.chunks)
	chunk = dh.chunks[chunk_id]
	update_vel_half!(chunk, Δt½)
	apply_bcs!(chunk, t)
	update_disp_and_pos!(chunk, Δt)
	end
	@threads :static for chunk_id in eachindex(dh.chunks)
	exchange_loc_to_halo!(dh, chunk_id)
	calc_force_density!(dh.chunks[chunk_id])
	end
	@threads :static for chunk_id in eachindex(dh.chunks)
	exchange_halo_to_loc!(dh, chunk_id)
	chunk = dh.chunks[chunk_id]
	calc_damage!(chunk)
	update_acc_and_vel!(chunk, Δt½)
	export_results(dh, options, chunk_id, n, t)
	end
	return nothing
	end

For MPI, each process uses its chunk, and then the processes need to communicate to exchange the information of points at the chunk border (halo points):

Peridynamics.jl/src/time_solvers/velocity_verlet.jl

Lines 250 to 264 in 9bb65ae

	function solve_timestep!(dh::AbstractMPIDataHandler, options::AbstractOptions, Δt::Float64,
	Δt½::Float64, n::Int)
	t = n * Δt
	chunk = dh.chunk
	@timeit_debug TO "update_vel_half!" update_vel_half!(chunk, Δt½)
	@timeit_debug TO "apply_bcs!" apply_bcs!(chunk, t)
	@timeit_debug TO "update_disp_and_pos!" update_disp_and_pos!(chunk, Δt)
	@timeit_debug TO "exchange_loc_to_halo!" exchange_loc_to_halo!(dh)
	@timeit_debug TO "calc_force_density!" calc_force_density!(chunk)
	@timeit_debug TO "exchange_halo_to_loc!" exchange_halo_to_loc!(dh)
	@timeit_debug TO "calc_damage!" calc_damage!(chunk)
	@timeit_debug TO "update_acc_and_vel!" update_acc_and_vel!(chunk, Δt½)
	@timeit_debug TO "export_results" export_results(dh, options, n, t)
	return nothing
	end

If you have further questions regarding the parallelism, please do not hesitate to contact me.

On volume correction:

I did not include the volume correction because it needs a uniformly distributed point cloud to be correct, which is not guaranteed with the current approach of our package. If someone uses an arbitrary mesh with Abaqus and converts that into a point cloud, the volume correction will introduce errors, and the point spacing is unknown.

However, from our observations, the correspondence formulation is also very precise with arbitrary meshes, and the volume error had only a very negligible effect on the simulation results, so we did not prioritize implementing it.

Do you have observations on how much the volume error influences the simulation results in your thermo-mechanical coupled simulations? If you use the correspondence formulation instead of the bond-based formulation, all the corrections could theoretically be dismissed.

However, I have some ideas about how to guarantee that a point cloud is uniformly distributed with a given point spacing. This would again break the current API but make it possible to use volume correction if a point cloud is guaranteed to be uniform and to ignore it otherwise. I will open an issue and ping you if my idea is clearer.

[ShiWeiHuH]

Hi, @kaipartmann It's been a while!

I noticed that PD has a new version with many great features. That's fantastic!

I have a couple of questions:

Currently, I'm using single-thread computation with one region. If I want to switch to using 64 threads, should I set the environment variable NUM=64 first? Will this automatically divide the region into 64 blocks?

Can computation be resumed if it shuts down unexpectedly?

[kaipartmann]

Hi @ShiWeiHuH, I am very happy you like the new features and changes we made in the last weeks!

Threads

To use multiple threads you can use a few different approaches, as also mentioned [here].(https://docs.julialang.org/en/v1/manual/multi-threading/#Starting-Julia-with-multiple-threads)

Use Julia command line options

julia -t 64
julia --threads 64

Specify the number of threads in the VSCode settings
If you use VSCode, then you can specify the number of threads for a Julia REPL in the settings.

Setting an environment variable

JULIA_NUM_THREADS=64

You can check how many threads you use with

Threads.nthreads()

The body is then divided in different blocks with distribute_equally

Peridynamics.jl/src/discretization/decomposition.jl

Lines 18 to 44 in bc81309

	function distribute_equally(n_elems::Int, n_chunks::Int)
	if n_elems >= n_chunks
	chunk_size = div(n_elems, n_chunks)
	remainder = rem(n_elems, n_chunks)
	sidx = zeros(Int64, n_chunks + 1)
	for i in 1:(n_chunks + 1)
	sidx[i] += (i - 1) * chunk_size + 1
	if i <= remainder
	sidx[i] += i - 1
	else
	sidx[i] += remainder
	end
	end
	decomp = fill(0:0, n_chunks)
	for i in 1:n_chunks
	decomp[i] = sidx[i]:(sidx[i + 1] - 1)
	end
	return decomp
	else
	sidx = [1:(n_elems + 1);]
	decomp = fill(0:-1, n_chunks)
	for i in 1:n_elems
	decomp[i] = sidx[i]:(sidx[i + 1] - 1)
	end
	return decomp
	end
	end

You can check the different regions with

body = ...
n_chunks = 64
point_decomp = Peridynamics.PointDecomposition(body, n_chunks)
point_decomp.decomp # Vector with n_chunks elements and a range of the point ids for each region

Resume simulations

For now, it is not possible to resume simulations. The current approach is to capture as many bugs as possible before submitting the simulation job. If a simulation stops unexpectedly, all exported data up to this point is preserved.

However, the main reason I have not considered it for now is performance. It would be necessary to export the current state of the data handler regularly, which would then take a significant fraction of the simulation time (I estimate up to 50% for some setups).

I made an issue #144, so maybe there is a way to do it efficiently, and it will be included in a future version!