diff --git a/src/Peridynamics.jl b/src/Peridynamics.jl index a1258344..b60baeec 100644 --- a/src/Peridynamics.jl +++ b/src/Peridynamics.jl @@ -9,7 +9,7 @@ end export BBMaterial, CKIMaterial, NOSBMaterial, OSBMaterial, Body, point_set!, failure_permit!, material!, velocity_bc!, velocity_ic!, forcedensity_bc!, precrack!, VelocityVerlet, MultibodySetup, contact!, Job, read_vtk, uniform_box, submit, - process_each_export + process_each_export, mpi_isroot const MPI_INITIALIZED = Ref(false) const MPI_RANK = Ref(-1) @@ -25,6 +25,8 @@ const QUIET = Ref(false) QUIET[] = b return nothing end +@inline mpi_chunk_id() = mpi_rank() + 1 +@inline mpi_isroot() = mpi_rank() == 0 const TO = TimerOutput() @@ -113,6 +115,7 @@ include("core/job.jl") include("core/submit.jl") include("core/halo_exchange.jl") include("core/threads_data_handler.jl") +include("core/mpi_data_handler.jl") include("time_solvers/solve_velocity_verlet.jl") diff --git a/src/auxiliary/io.jl b/src/auxiliary/io.jl index a65b6292..4db4736b 100644 --- a/src/auxiliary/io.jl +++ b/src/auxiliary/io.jl @@ -76,7 +76,6 @@ function check_storage_fields(::Type{S}, fields::Vector{Symbol}) where {S} return nothing end - function export_results(dh::AbstractDataHandler, options::ExportOptions, chunk_id::Int, timestep::Int, time::Float64) options.exportflag || return nothing @@ -96,7 +95,8 @@ end function _export_results(b::AbstractBodyChunk, chunk_id::Int, n_chunks::Int, options::ExportOptions, n::Int, t::Float64) - filename = @sprintf("timestep_%05d", n) + # filename = @sprintf("timestep_%05d", n) + filename = joinpath(options.vtk, @sprintf("timestep_%05d", n)) position = get_loc_position(b) pvtk_grid(filename, position, b.cells; part=chunk_id, nparts=n_chunks) do vtk for field in options.fields diff --git a/src/core/halo_exchange.jl b/src/core/halo_exchange.jl index 93472f86..a9070ba6 100644 --- a/src/core/halo_exchange.jl +++ b/src/core/halo_exchange.jl @@ -1,64 +1,83 @@ +function exchange!(dest::Matrix{T}, src::Matrix{T}, dest_idxs::Vector{Int}, + src_idxs::Vector{Int}) where {T} + for i in eachindex(dest_idxs) + for d in axes(dest, 1) + @inbounds dest[d, dest_idxs[i]] = src[d, src_idxs[i]] + end + end + return nothing +end -struct HaloExchange - src_chunk_id::Int - dest_chunk_id::Int - src_idxs::Vector{Int} - dest_idxs::Vector{Int} +function exchange_from_buf!(dest::T, buf::T, dest_idxs::Vector{Int}) where {T<:Matrix} + for i in eachindex(dest_idxs) + for d in axes(dest, 1) + dest[d, dest_idxs[i]] = buf[d, i] + end + end + return nothing end -function find_halo_exchanges(chunks::Vector{B}) where {B<:AbstractBodyChunk} - has_read = has_read_halos(chunks) - has_write = has_write_halos(chunks) - read_halo_exs = Vector{Vector{HaloExchange}}(undef, length(chunks)) - write_halo_exs = Vector{Vector{HaloExchange}}(undef, length(chunks)) - @threads :static for chunk_id in eachindex(chunks) - read_exs, write_exs = find_exs(chunks, chunk_id, has_read, has_write) - read_halo_exs[chunk_id] = read_exs - write_halo_exs[chunk_id] = write_exs +function exchange_to_buf!(buf::T, src::T, src_idxs::Vector{Int}) where {T<:Matrix} + for i in eachindex(src_idxs) + for d in axes(buf, 1) + buf[d, i] = src[d, src_idxs[i]] + end end - reorder_write_exs!(write_halo_exs) - return read_halo_exs, write_halo_exs + return nothing end -function has_read_halos(chunk::AbstractBodyChunk) - hrfs = get_halo_read_fields(chunk.store) - isempty(hrfs) && return false - return true +function exchange!(dest::Vector{T}, src::Vector{T}, dest_idxs::Vector{Int}, + src_idxs::Vector{Int}) where {T} + for i in eachindex(dest_idxs) + @inbounds dest[dest_idxs[i]] = src[src_idxs[i]] + end + return nothing end -function has_read_halos(chunks::Vector{B}) where {B<:AbstractBodyChunk} - return has_read_halos(first(chunks)) +function exchange_from_buf!(dest::T, buf::T, dest_idxs::Vector{Int}) where {T<:Vector} + for i in eachindex(dest_idxs) + dest[dest_idxs[i]] = buf[i] + end + return nothing end -function has_write_halos(chunk::AbstractBodyChunk) - hwfs = get_halo_write_fields(chunk.store) - isempty(hwfs) && return false - return true +function exchange_to_buf!(buf::T, src::T, src_idxs::Vector{Int}) where {T<:Vector} + for i in eachindex(src_idxs) + buf[i] = src[src_idxs[i]] + end + return nothing end -function has_write_halos(chunks::Vector{B}) where {B<:AbstractBodyChunk} - return has_write_halos(first(chunks)) +function exchange_add!(dest::Matrix{T}, src::Matrix{T}, dest_idxs::Vector{Int}, + src_idxs::Vector{Int}) where {T<:Number} + for i in eachindex(dest_idxs) + for d in axes(dest, 1) + @inbounds dest[d, dest_idxs[i]] += src[d, src_idxs[i]] + end + end + return nothing end -function find_exs(chunks::Vector{B}, chunk_id::Int, hasread::Bool, - haswrite::Bool) where {B<:AbstractBodyChunk} - readexs = Vector{HaloExchange}() - writeexs = Vector{HaloExchange}() - chunk = chunks[chunk_id] - for (halo_chunk_id, idxs) in chunk.ch.halo_by_src - halo_chunk = chunks[halo_chunk_id] - halo_idxs = chunk.ch.point_ids[idxs] - localize!(halo_idxs, halo_chunk.ch.localizer) - hasread && push!(readexs, HaloExchange(halo_chunk_id, chunk_id, halo_idxs, idxs)) - haswrite && push!(writeexs, HaloExchange(chunk_id, halo_chunk_id, idxs, halo_idxs)) +function exchange_from_buf_add!(dest::T, buf::T, dest_idxs::Vector{Int}) where {T<:Matrix} + for i in eachindex(dest_idxs) + for d in axes(dest, 1) + dest[d, dest_idxs[i]] += buf[d, i] + end end - return readexs, writeexs + return nothing +end + +function exchange_add!(dest::Vector{T}, src::Vector{T}, dest_idxs::Vector{Int}, + src_idxs::Vector{Int}) where {T<:Number} + for i in eachindex(dest_idxs) + @inbounds dest[dest_idxs[i]] += src[src_idxs[i]] + end + return nothing end -function reorder_write_exs!(write_halo_exs::Vector{Vector{HaloExchange}}) - all_write_exs = reduce(vcat, write_halo_exs) - @threads :static for chunk_id in eachindex(write_halo_exs) - write_halo_exs[chunk_id] = filter(x -> x.dest_chunk_id == chunk_id, all_write_exs) +function exchange_from_buf_add!(dest::T, buf::T, dest_idxs::Vector{Int}) where {T<:Vector} + for i in eachindex(dest_idxs) + dest[dest_idxs[i]] += buf[i] end return nothing end diff --git a/src/core/job.jl b/src/core/job.jl index 85ed1e62..38ffffd7 100644 --- a/src/core/job.jl +++ b/src/core/job.jl @@ -15,6 +15,6 @@ end function Job(spatial_setup::S, time_solver::T; kwargs...) where {S,T} o = Dict{Symbol,Any}(kwargs) check_kwargs(o, JOB_KWARGS) - options = get_export_options(storage_type(spatial_setup.mat, time_solver), o) + options = get_export_options(storage_type(spatial_setup, time_solver), o) return Job(spatial_setup, time_solver, options) end diff --git a/src/core/mpi_data_handler.jl b/src/core/mpi_data_handler.jl new file mode 100644 index 00000000..231b5cda --- /dev/null +++ b/src/core/mpi_data_handler.jl @@ -0,0 +1,250 @@ +struct HaloExchangeBuf{T,K} + src_chunk_id::Int + dest_chunk_id::Int + src_idxs::Vector{Int} + dest_idxs::Vector{Int} + bufs::T + tags::K + function HaloExchangeBuf(src_chunk_id::Int, dest_chunk_id::Int, + src_idxs::AbstractVector{Int}, dest_idxs::AbstractVector{Int}, + bufs::T, tags::K) where {T,K} + return new{T,K}(src_chunk_id, dest_chunk_id, src_idxs, dest_idxs, bufs, tags) + end +end + +struct MPIHaloInfo + point_ids::Dict{Int,Vector{Int}} + halos_points::Dict{Int,Vector{Int}} + localizers::Dict{Int,Dict{Int,Int}} + hidxs_by_src::Dict{Int,Dict{Int,Vector{Int}}} +end + +struct MPIDataHandler{C<:AbstractBodyChunk,RB,WB} <: AbstractDataHandler + chunk::C + read_halo_exs::Vector{HaloExchangeBuf{RB}} + write_halo_exs::Vector{HaloExchangeBuf{WB}} +end + +function MPIDataHandler(body::Body, time_solver::AbstractTimeSolver, + point_decomp::PointDecomposition) + n_param = length(body.point_params) + v = n_param == 1 ? Val{1}() : Val{2}() + return MPIDataHandler(body, time_solver, point_decomp, v) +end + +function MPIDataHandler(body::Body, time_solver::AbstractTimeSolver, + point_decomp::PointDecomposition, v::Val{N}) where {N} + chunk = chop_body_mpi(body, time_solver, point_decomp, v) + halo_infos = get_halo_infos(chunk, point_decomp, body.n_points) + read_halo_exs, write_halo_exs = find_halo_exchange_bufs(chunk, halo_infos) + return MPIDataHandler(chunk, read_halo_exs, write_halo_exs) +end + +function MPIDataHandler(multibody::MultibodySetup, time_solver::AbstractTimeSolver, + point_decomp::PointDecomposition) + error("MultibodySetup not yet implemented!\n") +end + +function chop_body_mpi(body::Body{M,P}, ts::T, pd::PointDecomposition, + v::Val{N}) where {M,P,T,N} + body_chunks = _chop_body_mpi(body, ts, pd, v) + return body_chunks +end + +function _chop_body_mpi(body::Body{M,P}, ts::T, pd::PointDecomposition, + ::Val{1}) where {M,P,T} + chunk = BodyChunk(body, ts, pd, mpi_chunk_id()) + apply_precracks!(chunk, body) + apply_initial_conditions!(chunk, body) + return chunk +end + +function _chop_body_mpi(body::Body{M,P}, ts::T, pd::PointDecomposition, + ::Val{N}) where {M,P,T,N} + chunk = MultiParamBodyChunk(body, ts, pd, mpi_chunk_id()) + apply_precracks!(chunk, body) + apply_initial_conditions!(chunk, body) + return chunk +end + +function get_halo_infos(chunk::AbstractBodyChunk, pd::PointDecomposition, n_points::Int) + point_ids = get_all_point_ids(chunk, n_points) + halo_points = Dict{Int,Vector{Int}}() + localizers = Dict{Int,Dict{Int,Int}}() + hidxs_by_src = Dict{Int,Dict{Int,UnitRange{Int}}}() + for (chunk_id, loc_point_ids) in point_ids + n_loc_points = length(pd.decomp[chunk_id]) + halos = loc_point_ids[n_loc_points+1:end] + halo_sources = [pd.point_src[i] for i in halos] + halo_points[chunk_id] = halos + localizers[chunk_id] = find_localizer(loc_point_ids) + hidxs_by_src[chunk_id] = get_hidxs_by_src(halo_sources, n_loc_points) + end + return MPIHaloInfo(point_ids, halo_points, localizers, hidxs_by_src) +end + +function get_all_point_ids(chunk::AbstractBodyChunk, n_points::Int) + loc_point_ids = zeros(Int, n_points) + for i in eachindex(chunk.ch.point_ids) + loc_point_ids[i] = chunk.ch.point_ids[i] + end + _point_ids = MPI.Allgather(loc_point_ids, mpi_comm()) + point_ids = reshape(_point_ids, n_points, mpi_nranks()) + all_point_ids = Dict{Int,Vector{Int}}() + for chunk_id in axes(point_ids, 2) + all_point_ids[chunk_id] = filter(x -> x > 0, point_ids[:, chunk_id]) + end + return all_point_ids +end + +function get_halos_by_chunk_id(halo_chunk_ids::Vector{Int}) + halos_by_chunk_id = Dict{Int,Vector{Int}}() + unique_chunk_ids = unique(halo_chunk_ids) + for chunk_id in unique_chunk_ids + ids = findall(x -> x == chunk_id, halo_chunk_ids) + halos_by_chunk_id[chunk_id] = ids + end + return halos_by_chunk_id +end + +function find_halo_exchange_bufs(chunk::AbstractBodyChunk, halo_infos::MPIHaloInfo) + has_read = has_read_halos(chunk) + has_write = has_write_halos(chunk) + read_exs, write_exs = find_exs(chunk, halo_infos, has_read, has_write) + return read_exs, write_exs +end + +function has_read_halos(chunk::AbstractBodyChunk) + hrfs = get_halo_read_fields(chunk.store) + isempty(hrfs) && return false + return true +end + +function has_write_halos(chunk::AbstractBodyChunk) + hwfs = get_halo_write_fields(chunk.store) + isempty(hwfs) && return false + return true +end + +function find_exs(chunk::AbstractBodyChunk, halo_info::MPIHaloInfo, hasread::Bool, + haswrite::Bool) + RB = get_halo_read_buftype(chunk.store) + WB = get_halo_write_buftype(chunk.store) + readexs = Vector{HaloExchangeBuf{RB}}() + writeexs = Vector{HaloExchangeBuf{WB}}() + # current_chunk_id = mpi_chunk_id() + for (chunk_id, hidxs_by_src) in halo_info.hidxs_by_src + for (halo_chunk_id, idxs) in hidxs_by_src + halo_idxs = halo_info.point_ids[chunk_id][idxs] + localize!(halo_idxs, halo_info.localizers[halo_chunk_id]) + if hasread + src_id = halo_chunk_id + dest_id = chunk_id + src_idxs = halo_idxs + dest_idxs = idxs + bufs = get_halo_read_bufs(chunk.store, length(idxs)) + tags = get_tags(bufs, chunk_id * mpi_nranks()) + heb = HaloExchangeBuf(src_id, dest_id, src_idxs, dest_idxs, bufs, tags) + push!(readexs, heb) + end + if haswrite + src_id = chunk_id + dest_id = halo_chunk_id + src_idxs = idxs + dest_idxs = halo_idxs + bufs = get_halo_read_bufs(chunk.store, length(idxs)) + tags = get_tags(bufs, chunk_id * mpi_nranks()) + heb = HaloExchangeBuf(src_id, dest_id, src_idxs, dest_idxs, bufs, tags) + push!(writeexs, heb) + end + end + end + return readexs, writeexs +end + +@inline function get_halo_read_buftype(s::AbstractStorage) + return typeof(get_halo_read_fields(s)) +end + +@inline function get_halo_write_buftype(s::AbstractStorage) + return typeof(get_halo_write_fields(s)) +end + +@inline function get_halo_read_bufs(s::AbstractStorage, n_halo_points::Int) + return Tuple(get_buf(a, n_halo_points) for a in get_halo_read_fields(s)) +end + +@inline function get_halo_write_bufs(s::AbstractStorage, n_halo_points::Int) + return Tuple(get_buf(a, n_halo_points) for a in get_halo_write_fields(s)) +end + +@inline get_buf(a::Matrix{T}, n::Int) where {T} = zeros(T, size(a,1), n) +@inline get_buf(::Vector{T}, n::Int) where {T} = zeros(T, n) + +@inline function get_tags(bufs, i) + return Tuple(i + j for j in eachindex(bufs)) +end + +function calc_stable_timestep(dh::MPIDataHandler, safety_factor::Float64) + _Δt = calc_timestep(dh.chunk) + Δt = MPI.Allreduce(_Δt, MPI.MIN, mpi_comm()) + return Δt * safety_factor +end + +function exchange_read_fields!(dh::MPIDataHandler) + for he in dh.read_halo_exs + if mpi_chunk_id() == he.src_chunk_id + src_fields = get_halo_read_fields(dh.chunk.store) + for i in eachindex(src_fields, he.bufs) + exchange_to_buf!(he.bufs[i], src_fields[i], he.src_idxs) + MPI.Isend(he.bufs[i], mpi_comm(); dest=he.dest_chunk_id-1, tag=he.tags[i]) + end + end + end + for he in dh.read_halo_exs + if mpi_chunk_id() == he.dest_chunk_id + dest_fields = get_halo_read_fields(dh.chunk.store) + for i in eachindex(dest_fields, he.bufs) + MPI.Recv!(he.bufs[i], mpi_comm(); source=he.src_chunk_id-1, tag=he.tags[i]) + exchange_from_buf!(dest_fields[i], he.bufs[i], he.dest_idxs) + end + end + end + return nothing +end + +function exchange_write_fields!(dh::MPIDataHandler) + for he in dh.write_halo_exs + if mpi_chunk_id() == he.src_chunk_id + src_fields = get_halo_write_fields(dh.chunk.store) + for i in eachindex(src_fields, he.bufs) + exchange_to_buf!(he.bufs[i], src_fields[i], he.src_idxs) + MPI.Isend(he.bufs[i], mpi_comm(); dest=he.dest_chunk_id-1, tag=he.tags[i]) + end + end + end + for he in dh.write_halo_exs + if mpi_chunk_id() == he.dest_chunk_id + dest_fields = get_halo_write_fields(dh.chunk.store) + for i in eachindex(dest_fields, he.bufs) + MPI.Recv!(he.bufs[i], mpi_comm(); source=he.src_chunk_id-1, tag=he.tags[i]) + exchange_from_buf_add!(dest_fields[i], he.bufs[i], he.dest_idxs) + end + end + end + return nothing +end + +function export_results(dh::MPIDataHandler, options::ExportOptions, n::Int, t::Float64) + options.exportflag || return nothing + if mod(n, options.freq) == 0 + _export_results(dh.chunk, mpi_chunk_id(), mpi_nranks(), options, n, t) + end + return nothing +end + +function export_reference_results(dh::MPIDataHandler, options::ExportOptions) + options.exportflag || return nothing + _export_results(dh.chunk, mpi_chunk_id(), mpi_nranks(), options, 0, 0.0) + return nothing +end diff --git a/src/core/submit.jl b/src/core/submit.jl index 716c48ec..ed4d6430 100644 --- a/src/core/submit.jl +++ b/src/core/submit.jl @@ -26,9 +26,14 @@ function get_submit_options(o::Dict{Symbol,Any}) return quiet end -function submit_mpi(::Job) - error("functionality for MPI simulations not yet implemented!\n") - return nothing +function submit_mpi(job::Job) + point_decomp = PointDecomposition(job.spatial_setup, mpi_nranks()) + mdh = MPIDataHandler(job.spatial_setup, job.time_solver, point_decomp) + init_time_solver!(job.time_solver, mdh) + _pwd = init_export(job.options) + solve!(mdh, job.time_solver, job.options) + finish_export(job.options, _pwd) + return mdh end function submit_threads(job::Job, n_chunks::Int) @@ -45,11 +50,11 @@ function init_export(options) options.exportflag || return "" _pwd = pwd() mkpath(options.vtk) - cd(options.vtk) + # cd(options.vtk) return _pwd end function finish_export(options, _pwd) - options.exportflag && cd(_pwd) + # options.exportflag && cd(_pwd) return nothing end diff --git a/src/core/threads_data_handler.jl b/src/core/threads_data_handler.jl index 3c96e7e3..73a7b65d 100644 --- a/src/core/threads_data_handler.jl +++ b/src/core/threads_data_handler.jl @@ -1,3 +1,10 @@ +struct HaloExchange + src_chunk_id::Int + dest_chunk_id::Int + src_idxs::Vector{Int} + dest_idxs::Vector{Int} +end + struct ThreadsDataHandler{C<:AbstractBodyChunk} <: AbstractDataHandler n_chunks::Int chunks::Vector{C} @@ -8,7 +15,7 @@ end function ThreadsDataHandler(body::Body, time_solver::AbstractTimeSolver, point_decomp::PointDecomposition) n_param = length(body.point_params) - v = length(n_param) == 1 ? Val{1}() : Val{2}() + v = n_param == 1 ? Val{1}() : Val{2}() return ThreadsDataHandler(body, time_solver, point_decomp, v) end @@ -25,6 +32,87 @@ function ThreadsDataHandler(multibody::MultibodySetup, time_solver::AbstractTime error("MultibodySetup not yet implemented!\n") end +function chop_body_threads(body::Body{M,P}, ts::T, pd::PointDecomposition, + v::Val{N}) where {M,P,T,N} + D = discretization_type(body.mat) + S = storage_type(body.mat, ts) + body_chunks = _chop_body_threads(body, ts, pd, D, S, v) + return body_chunks +end + +function _chop_body_threads(body::Body{M,P}, ts::T, pd::PointDecomposition, ::Type{D}, + ::Type{S}, ::Val{1}) where {M,P,D,S,T} + body_chunks = Vector{BodyChunk{M,P,D,S}}(undef, pd.n_chunks) + + @threads :static for chunk_id in eachindex(pd.decomp) + body_chunk = BodyChunk(body, ts, pd, chunk_id) + apply_precracks!(body_chunk, body) + apply_initial_conditions!(body_chunk, body) + body_chunks[chunk_id] = body_chunk + end + + return body_chunks +end + +function _chop_body_threads(body::Body{M,P}, ts::T, pd::PointDecomposition, ::Type{D}, + ::Type{S}, ::Val{N}) where {M,P,D,S,T,N} + body_chunks = Vector{MultiParamBodyChunk{M,P,D,S}}(undef, pd.n_chunks) + + @threads :static for chunk_id in eachindex(pd.decomp) + body_chunk = MultiParamBodyChunk(body, ts, pd, chunk_id) + apply_precracks!(body_chunk, body) + apply_initial_conditions!(body_chunk, body) + body_chunks[chunk_id] = body_chunk + end + + return body_chunks +end + +function find_halo_exchanges(chunks::Vector{B}) where {B<:AbstractBodyChunk} + has_read = has_read_halos(chunks) + has_write = has_write_halos(chunks) + read_halo_exs = Vector{Vector{HaloExchange}}(undef, length(chunks)) + write_halo_exs = Vector{Vector{HaloExchange}}(undef, length(chunks)) + @threads :static for chunk_id in eachindex(chunks) + read_exs, write_exs = find_exs(chunks, chunk_id, has_read, has_write) + read_halo_exs[chunk_id] = read_exs + write_halo_exs[chunk_id] = write_exs + end + reorder_write_exs!(write_halo_exs) + return read_halo_exs, write_halo_exs +end + +function has_read_halos(chunks::Vector{B}) where {B<:AbstractBodyChunk} + return has_read_halos(first(chunks)) +end + +function has_write_halos(chunks::Vector{B}) where {B<:AbstractBodyChunk} + return has_write_halos(first(chunks)) +end + +function find_exs(chunks::Vector{B}, chunk_id::Int, hasread::Bool, + haswrite::Bool) where {B<:AbstractBodyChunk} + readexs = Vector{HaloExchange}() + writeexs = Vector{HaloExchange}() + chunk = chunks[chunk_id] + for (halo_chunk_id, idxs) in chunk.ch.hidxs_by_src + halo_chunk = chunks[halo_chunk_id] + halo_idxs = chunk.ch.point_ids[idxs] + localize!(halo_idxs, halo_chunk.ch.localizer) + hasread && push!(readexs, HaloExchange(halo_chunk_id, chunk_id, halo_idxs, idxs)) + haswrite && push!(writeexs, HaloExchange(chunk_id, halo_chunk_id, idxs, halo_idxs)) + end + return readexs, writeexs +end + +function reorder_write_exs!(write_halo_exs::Vector{Vector{HaloExchange}}) + all_write_exs = reduce(vcat, write_halo_exs) + @threads :static for chunk_id in eachindex(write_halo_exs) + write_halo_exs[chunk_id] = filter(x -> x.dest_chunk_id == chunk_id, all_write_exs) + end + return nothing +end + function calc_stable_timestep(dh::ThreadsDataHandler, safety_factor::Float64) Δt = zeros(length(dh.chunks)) @threads :static for chunk_id in eachindex(dh.chunks) @@ -57,38 +145,19 @@ function exchange_write_fields!(dh::ThreadsDataHandler, chunk_id::Int) return nothing end -function exchange!(dest::Matrix{T}, src::Matrix{T}, dest_idxs::Vector{Int}, - src_idxs::Vector{Int}) where {T} - for i in eachindex(dest_idxs) - for d in axes(dest, 1) - @inbounds dest[d, dest_idxs[i]] = src[d, src_idxs[i]] - end +function export_results(dh::ThreadsDataHandler, options::ExportOptions, chunk_id::Int, + timestep::Int, time::Float64) + options.exportflag || return nothing + if mod(timestep, options.freq) == 0 + _export_results(dh.chunks[chunk_id], chunk_id, dh.n_chunks, options, timestep, time) end return nothing end -function exchange!(dest::Vector{T}, src::Vector{T}, dest_idxs::Vector{Int}, - src_idxs::Vector{Int}) where {T} - for i in eachindex(dest_idxs) - @inbounds dest[dest_idxs[i]] = src[src_idxs[i]] - end - return nothing -end - -function exchange_add!(dest::Matrix{T}, src::Matrix{T}, dest_idxs::Vector{Int}, - src_idxs::Vector{Int}) where {T<:Number} - for i in eachindex(dest_idxs) - for d in axes(dest, 1) - @inbounds dest[d, dest_idxs[i]] += src[d, src_idxs[i]] - end - end - return nothing -end - -function exchange_add!(dest::Vector{T}, src::Vector{T}, dest_idxs::Vector{Int}, - src_idxs::Vector{Int}) where {T<:Number} - for i in eachindex(dest_idxs) - @inbounds dest[dest_idxs[i]] += src[src_idxs[i]] +function export_reference_results(dh::ThreadsDataHandler, options::ExportOptions) + options.exportflag || return nothing + @threads :static for chunk_id in eachindex(dh.chunks) + _export_results(dh.chunks[chunk_id], chunk_id, dh.n_chunks, options, 0, 0.0) end return nothing end diff --git a/src/discretizations/body.jl b/src/discretizations/body.jl index 6278842b..89312169 100644 --- a/src/discretizations/body.jl +++ b/src/discretizations/body.jl @@ -249,3 +249,5 @@ end @inline function get_point_param(b::Body, key::Symbol, i::Int) return getfield(b.point_params[b.params_map[i]], key) end + +@inline storage_type(b::Body, ts::AbstractTimeSolver) = storage_type(b.mat, ts) diff --git a/src/discretizations/body_chunk.jl b/src/discretizations/body_chunk.jl index 506d4a21..040be228 100644 --- a/src/discretizations/body_chunk.jl +++ b/src/discretizations/body_chunk.jl @@ -63,42 +63,6 @@ end return b.param end -function chop_body_threads(body::Body{M,P}, ts::T, pd::PointDecomposition, - v::Val{N}) where {M,P,T,N} - D = discretization_type(body.mat) - S = storage_type(body.mat, ts) - body_chunks = _chop_body_threads(body, ts, pd, D, S, v) - return body_chunks -end - -function _chop_body_threads(body::Body{M,P}, ts::T, pd::PointDecomposition, ::Type{D}, - ::Type{S}, ::Val{1}) where {M,P,D,S,T} - body_chunks = Vector{BodyChunk{M,P,D,S}}(undef, pd.n_chunks) - - @threads :static for chunk_id in eachindex(pd.decomp) - body_chunk = BodyChunk(body, ts, pd, chunk_id) - apply_precracks!(body_chunk, body) - apply_initial_conditions!(body_chunk, body) - body_chunks[chunk_id] = body_chunk - end - - return body_chunks -end - -function _chop_body_threads(body::Body{M,P}, ts::T, pd::PointDecomposition, ::Type{D}, - ::Type{S}, ::Val{N}) where {M,P,D,S,T,N} - body_chunks = Vector{MultiParamBodyChunk{M,P,D,S}}(undef, pd.n_chunks) - - @threads :static for chunk_id in eachindex(pd.decomp) - body_chunk = MultiParamBodyChunk(body, ts, pd, chunk_id) - apply_precracks!(body_chunk, body) - apply_initial_conditions!(body_chunk, body) - body_chunks[chunk_id] = body_chunk - end - - return body_chunks -end - @inline get_material_type(::MultiParamBodyChunk{M,P,D,S}) where {M,P,D,S} = M @inline get_material_type(::BodyChunk{M,P,D,S}) where {M,P,D,S} = M @inline get_material_type(::Vector{MultiParamBodyChunk{M,P,D,S}}) where {M,P,D,S} = M diff --git a/src/discretizations/chunk_handler.jl b/src/discretizations/chunk_handler.jl index a35a8568..a704549d 100644 --- a/src/discretizations/chunk_handler.jl +++ b/src/discretizations/chunk_handler.jl @@ -3,7 +3,7 @@ struct ChunkHandler point_ids::Vector{Int} loc_points::UnitRange{Int} halo_points::Vector{Int} - halo_by_src::Dict{Int,UnitRange{Int}} + hidxs_by_src::Dict{Int,UnitRange{Int}} localizer::Dict{Int,Int} end @@ -11,10 +11,10 @@ function ChunkHandler(bonds::Vector{Bond}, pd::PointDecomposition, chunk_id::Int loc_points = pd.decomp[chunk_id] n_loc_points = length(loc_points) halo_points = find_halo_points(bonds, loc_points) - halo_by_src = sort_halo_by_src!(halo_points, pd.point_src, length(loc_points)) + hidxs_by_src = sort_halo_by_src!(halo_points, pd.point_src, length(loc_points)) point_ids = vcat(loc_points, halo_points) localizer = find_localizer(point_ids) - return ChunkHandler(n_loc_points, point_ids, loc_points, halo_points, halo_by_src, + return ChunkHandler(n_loc_points, point_ids, loc_points, halo_points, hidxs_by_src, localizer) end @@ -35,18 +35,23 @@ function sort_halo_by_src!(halo_points::Vector{Int}, point_src::Dict{Int,Int}, idx_sorted = sortperm(halo_sources) halo_sources .= halo_sources[idx_sorted] halo_points .= halo_points[idx_sorted] + hidxs_by_src = get_hidxs_by_src(halo_sources, n_loc_points) + return hidxs_by_src +end - halo_by_src = Dict{Int,UnitRange{Int}}() +function get_hidxs_by_src(halo_sources::Vector{Int}, n_loc_points::Int) + @assert sort(halo_sources) == halo_sources + hidxs_by_src = Dict{Int,UnitRange{Int}}() unique_sources = unique(halo_sources) for source in unique_sources idxs = findall(x -> x == source, halo_sources) idx_begin = first(idxs) + n_loc_points idx_end = last(idxs) + n_loc_points - halo_by_src[source] = idx_begin:idx_end + hidxs_by_src[source] = idx_begin:idx_end end - return halo_by_src + return hidxs_by_src end function find_localizer(point_ids::Vector{Int}) diff --git a/src/discretizations/multibody_setup.jl b/src/discretizations/multibody_setup.jl index b1ae1888..fe30cfc8 100644 --- a/src/discretizations/multibody_setup.jl +++ b/src/discretizations/multibody_setup.jl @@ -55,3 +55,7 @@ function pre_submission_check(ms::MultibodySetup) #TODO: check if everything is defined for job submission! return nothing end + +@inline function storage_type(ms::MultibodySetup, ts::AbstractTimeSolver) + return storage_type(first(ms.bodies).mat, ts) +end diff --git a/src/time_solvers/solve_velocity_verlet.jl b/src/time_solvers/solve_velocity_verlet.jl index cd18b50d..4153753b 100644 --- a/src/time_solvers/solve_velocity_verlet.jl +++ b/src/time_solvers/solve_velocity_verlet.jl @@ -109,3 +109,34 @@ function solve_timestep!(dh::ThreadsDataHandler, options::ExportOptions, Δt::Fl end return nothing end + +function solve!(dh::MPIDataHandler, vv::VelocityVerlet, options::ExportOptions) + export_reference_results(dh, options) + Δt = vv.Δt + Δt½ = 0.5 * vv.Δt + if mpi_isroot() + p = Progress(vv.n_steps; dt=1, desc="solve...", color=:normal, barlen=20) + end + for n in 1:vv.n_steps + solve_timestep!(dh, options, Δt, Δt½, n) + mpi_isroot() && next!(p) + end + mpi_isroot() && finish!(p) + return dh +end + +function solve_timestep!(dh::MPIDataHandler, options::ExportOptions, Δt::Float64, + Δt½::Float64, n::Int) + t = n * Δt + chunk = dh.chunk + update_vel_half!(chunk, Δt½) + apply_bcs!(chunk, t) + update_disp_and_pos!(chunk, Δt) + exchange_read_fields!(dh) + calc_force_density!(chunk) + exchange_write_fields!(dh) + calc_damage!(chunk) + update_acc_and_vel!(chunk, Δt½) + export_results(dh, options, n, t) + return nothing +end diff --git a/test/discretizations/test_body_chunk.jl b/test/discretizations/test_body_chunk.jl index 4565822b..f6c18e35 100644 --- a/test/discretizations/test_body_chunk.jl +++ b/test/discretizations/test_body_chunk.jl @@ -34,7 +34,7 @@ @test bc.ch.point_ids == [1, 2, 3, 4] @test bc.ch.loc_points == pd.decomp[1] @test bc.ch.halo_points == [3, 4] - @test bc.ch.halo_by_src[2] == 3:4 + @test bc.ch.hidxs_by_src[2] == 3:4 for i in 1:4 @test bc.ch.localizer[i] == i end @@ -119,7 +119,7 @@ end @test b1.ch.point_ids == [1, 2, 3, 4] @test b1.ch.loc_points == 1:2 @test b1.ch.halo_points == [3, 4] - @test b1.ch.halo_by_src[2] == 3:4 + @test b1.ch.hidxs_by_src[2] == 3:4 for i in 1:4 @test b1.ch.localizer[i] == i end @@ -181,7 +181,7 @@ end @test b2.ch.point_ids == [3, 4, 1, 2] @test b2.ch.loc_points == [3, 4] @test b2.ch.halo_points == [1, 2] - @test b2.ch.halo_by_src[1] == 3:4 + @test b2.ch.hidxs_by_src[1] == 3:4 @test b2.ch.localizer[3] == 1 @test b2.ch.localizer[4] == 2 @test b2.ch.localizer[1] == 3 @@ -282,7 +282,7 @@ end @test b1.ch.point_ids == [1, 2, 3, 4] @test b1.ch.loc_points == 1:2 @test b1.ch.halo_points == [3, 4] - @test b1.ch.halo_by_src[2] == 3:4 + @test b1.ch.hidxs_by_src[2] == 3:4 for i in 1:4 @test b1.ch.localizer[i] == i end @@ -361,7 +361,7 @@ end @test b2.ch.point_ids == [3, 4, 1, 2] @test b2.ch.loc_points == [3, 4] @test b2.ch.halo_points == [1, 2] - @test b2.ch.halo_by_src[1] == 3:4 + @test b2.ch.hidxs_by_src[1] == 3:4 @test b2.ch.localizer[3] == 1 @test b2.ch.localizer[4] == 2 @test b2.ch.localizer[1] == 3 diff --git a/test/discretizations/test_bond_discretization.jl b/test/discretizations/test_bond_discretization.jl index c0e863ac..8f7f3f33 100644 --- a/test/discretizations/test_bond_discretization.jl +++ b/test/discretizations/test_bond_discretization.jl @@ -28,7 +28,7 @@ @test ch.point_ids == [1, 2, 3, 4] @test ch.loc_points == [1, 2] @test ch.halo_points == [3, 4] - @test ch.halo_by_src[2] == 3:4 + @test ch.hidxs_by_src[2] == 3:4 for i in 1:4 @test ch.localizer[i] == i @@ -53,7 +53,7 @@ @test ch.point_ids == [3, 4, 1, 2] @test ch.loc_points == [3, 4] @test ch.halo_points == [1, 2] - @test ch.halo_by_src[1] == 3:4 + @test ch.hidxs_by_src[1] == 3:4 @test ch.localizer[3] == 1 @test ch.localizer[4] == 2 @test ch.localizer[1] == 3 diff --git a/test/discretizations/test_chunk_handler.jl b/test/discretizations/test_chunk_handler.jl index f7059930..5096d251 100644 --- a/test/discretizations/test_chunk_handler.jl +++ b/test/discretizations/test_chunk_handler.jl @@ -36,10 +36,10 @@ end loc_points = 101:200 n_loc_points = length(loc_points) halo_points = Vector{Int}() - halo_by_src = Dict{Int,UnitRange{Int}}() + hidxs_by_src = Dict{Int,UnitRange{Int}}() localizer = Peridynamics.find_localizer(point_ids) ch = Peridynamics.ChunkHandler(n_loc_points, point_ids, loc_points, halo_points, - halo_by_src, localizer) + hidxs_by_src, localizer) point_set = [101, 110, 120, 210, 220] loc_point_set = Peridynamics.localize(point_set, ch) @@ -51,10 +51,10 @@ end loc_points = 101:200 n_loc_points = length(loc_points) halo_points = Vector{Int}() - halo_by_src = Dict{Int,UnitRange{Int}}() + hidxs_by_src = Dict{Int,UnitRange{Int}}() localizer = Peridynamics.find_localizer(point_ids) ch = Peridynamics.ChunkHandler(n_loc_points, point_ids, loc_points, halo_points, - halo_by_src, localizer) + hidxs_by_src, localizer) point_sets = Dict(:a => [101, 110, 120, 210, 220], :b => [1, 2, 3]) loc_point_sets = Peridynamics.localized_point_sets(point_sets, ch) @@ -94,7 +94,7 @@ end @test ch.point_ids == [1, 2, 3, 4] @test ch.loc_points == 1:2 @test ch.halo_points == [3, 4] - @test ch.halo_by_src[2] == 3:4 + @test ch.hidxs_by_src[2] == 3:4 @test ch.localizer[1] == 1 @test ch.localizer[2] == 2 @test ch.localizer[3] == 3 @@ -104,7 +104,7 @@ end @test ch.point_ids == [3, 4, 2, 1] @test ch.loc_points == 3:4 @test ch.halo_points == [2, 1] - @test ch.halo_by_src[1] == 3:4 + @test ch.hidxs_by_src[1] == 3:4 @test ch.localizer[1] == 4 @test ch.localizer[2] == 3 @test ch.localizer[3] == 1 @@ -115,9 +115,9 @@ end @test ch.point_ids == [1, 2, 3, 4] @test ch.loc_points == 1:1 @test ch.halo_points == [2, 3, 4] - @test ch.halo_by_src[2] == 2:2 - @test ch.halo_by_src[3] == 3:3 - @test ch.halo_by_src[4] == 4:4 + @test ch.hidxs_by_src[2] == 2:2 + @test ch.hidxs_by_src[3] == 3:3 + @test ch.hidxs_by_src[4] == 4:4 @test ch.localizer[1] == 1 @test ch.localizer[2] == 2 @test ch.localizer[3] == 3 @@ -127,9 +127,9 @@ end @test ch.point_ids == [2, 1, 3, 4] @test ch.loc_points == 2:2 @test ch.halo_points == [1, 3, 4] - @test ch.halo_by_src[1] == 2:2 - @test ch.halo_by_src[3] == 3:3 - @test ch.halo_by_src[4] == 4:4 + @test ch.hidxs_by_src[1] == 2:2 + @test ch.hidxs_by_src[3] == 3:3 + @test ch.hidxs_by_src[4] == 4:4 @test ch.localizer[1] == 2 @test ch.localizer[2] == 1 @test ch.localizer[3] == 3 diff --git a/test/integration/mpi_threads_comparison.jl b/test/integration/mpi_threads_comparison.jl new file mode 100644 index 00000000..1c66b9e4 --- /dev/null +++ b/test/integration/mpi_threads_comparison.jl @@ -0,0 +1,200 @@ +@testitem "MPI-Threads comparison BBMaterial" begin + root = joinpath(@__DIR__, "temp_mpi_threads_comparison") + path_threads = joinpath(root, "results_threads") + path_threads_vtk = joinpath(path_threads, "vtk") + path_mpi = joinpath(root, "results_mpi") + path_mpi_vtk = joinpath(path_mpi, "vtk") + + function sim_bb(N::Int, path::String) + l, Δx, δ, a = 1.0, 1/N, 3.015/N, 0.5 + pos, vol = uniform_box(l, l, 0.1l, Δx) + ids = sortperm(pos[2,:]) + b = Body(BBMaterial(), pos[:, ids], vol[ids]) + material!(b; horizon=3.015Δx, E=2.1e5, rho=8e-6, Gc=2.7) + point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, b, :set_a) + point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, b, :set_b) + precrack!(b, :set_a, :set_b) + point_set!(p -> p[2] > l/2-Δx, b, :set_top) + point_set!(p -> p[2] < -l/2+Δx, b, :set_bottom) + velocity_bc!(t -> -30, b, :set_bottom, :y) + velocity_bc!(t -> 30, b, :set_top, :y) + vv = VelocityVerlet(steps=100) + job = Job(b, vv; path=path, freq=50) + submit(job) + return nothing + end + sim_bb(30, path_threads) + + mpi_cmd = """ + using Peridynamics + function sim_bb(N::Int, path::String) + l, Δx, δ, a = 1.0, 1/N, 3.015/N, 0.5 + pos, vol = uniform_box(l, l, 0.1l, Δx) + ids = sortperm(pos[2,:]) + b = Body(BBMaterial(), pos[:, ids], vol[ids]) + material!(b; horizon=3.015Δx, E=2.1e5, rho=8e-6, Gc=2.7) + point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, b, :set_a) + point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, b, :set_b) + precrack!(b, :set_a, :set_b) + point_set!(p -> p[2] > l/2-Δx, b, :set_top) + point_set!(p -> p[2] < -l/2+Δx, b, :set_bottom) + velocity_bc!(t -> -30, b, :set_bottom, :y) + velocity_bc!(t -> 30, b, :set_top, :y) + vv = VelocityVerlet(steps=100) + job = Job(b, vv; path=path, freq=50) + submit(job) + return nothing + end + sim_bb(30, "$path_mpi") + """ + run(`$(Peridynamics.MPI.mpiexec()) -n 2 $(Base.julia_cmd()) --project -e $(mpi_cmd)`) + + @test isdir(path_threads_vtk) + @test isdir(path_mpi_vtk) + vtk_files_threads = Peridynamics.find_vtk_files(path_threads_vtk) + vtk_files_mpi = Peridynamics.find_vtk_files(path_mpi_vtk) + @test length(vtk_files_mpi) == length(vtk_files_threads) == 3 + for i in eachindex(vtk_files_threads, vtk_files_mpi) + res_threads = read_vtk(vtk_files_threads[i]) + res_mpi = read_vtk(vtk_files_mpi[i]) + for key in keys(res_threads) + @test res_threads[key] ≈ res_mpi[key] + end + end + + rm(root; recursive=true, force=true) +end + +@testitem "MPI-Threads comparison OSBMaterial" begin + root = joinpath(@__DIR__, "temp_mpi_threads_comparison") + path_threads = joinpath(root, "results_threads") + path_threads_vtk = joinpath(path_threads, "vtk") + path_mpi = joinpath(root, "results_mpi") + path_mpi_vtk = joinpath(path_mpi, "vtk") + + function sim_osb(N::Int, path::String) + l, Δx, δ, a = 1.0, 1/N, 3.015/N, 0.5 + pos, vol = uniform_box(l, l, 0.1l, Δx) + ids = sortperm(pos[2,:]) + b = Body(OSBMaterial(), pos[:, ids], vol[ids]) + material!(b; horizon=3.015Δx, E=2.1e5, nu=0.25, rho=8e-6, Gc=2.7) + point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, b, :set_a) + point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, b, :set_b) + precrack!(b, :set_a, :set_b) + point_set!(p -> p[2] > l/2-Δx, b, :set_top) + point_set!(p -> p[2] < -l/2+Δx, b, :set_bottom) + velocity_bc!(t -> -30, b, :set_bottom, :y) + velocity_bc!(t -> 30, b, :set_top, :y) + vv = VelocityVerlet(steps=100) + job = Job(b, vv; path=path, freq=50) + submit(job) + return nothing + end + sim_osb(30, path_threads) + + mpi_cmd = """ + using Peridynamics + function sim_osb(N::Int, path::String) + l, Δx, δ, a = 1.0, 1/N, 3.015/N, 0.5 + pos, vol = uniform_box(l, l, 0.1l, Δx) + ids = sortperm(pos[2,:]) + b = Body(OSBMaterial(), pos[:, ids], vol[ids]) + material!(b; horizon=3.015Δx, E=2.1e5, nu=0.25, rho=8e-6, Gc=2.7) + point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, b, :set_a) + point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, b, :set_b) + precrack!(b, :set_a, :set_b) + point_set!(p -> p[2] > l/2-Δx, b, :set_top) + point_set!(p -> p[2] < -l/2+Δx, b, :set_bottom) + velocity_bc!(t -> -30, b, :set_bottom, :y) + velocity_bc!(t -> 30, b, :set_top, :y) + vv = VelocityVerlet(steps=100) + job = Job(b, vv; path=path, freq=50) + submit(job) + return nothing + end + sim_osb(30, "$path_mpi") + """ + run(`$(Peridynamics.MPI.mpiexec()) -n 2 $(Base.julia_cmd()) --project -e $(mpi_cmd)`) + + @test isdir(path_threads_vtk) + @test isdir(path_mpi_vtk) + vtk_files_threads = Peridynamics.find_vtk_files(path_threads_vtk) + vtk_files_mpi = Peridynamics.find_vtk_files(path_mpi_vtk) + @test length(vtk_files_mpi) == length(vtk_files_threads) == 3 + for i in eachindex(vtk_files_threads, vtk_files_mpi) + res_threads = read_vtk(vtk_files_threads[i]) + res_mpi = read_vtk(vtk_files_mpi[i]) + for key in keys(res_threads) + @test res_threads[key] ≈ res_mpi[key] + end + end + + rm(root; recursive=true, force=true) +end + +@testitem "MPI-Threads comparison NOSBMaterial" begin + root = joinpath(@__DIR__, "temp_mpi_threads_comparison") + path_threads = joinpath(root, "results_threads") + path_threads_vtk = joinpath(path_threads, "vtk") + path_mpi = joinpath(root, "results_mpi") + path_mpi_vtk = joinpath(path_mpi, "vtk") + + function sim_nosb(N::Int, path::String) + l, Δx, δ, a = 1.0, 1/N, 3.015/N, 0.5 + pos, vol = uniform_box(l, l, 0.1l, Δx) + ids = sortperm(pos[2,:]) + b = Body(NOSBMaterial(), pos[:, ids], vol[ids]) + material!(b; horizon=3.015Δx, E=2.1e5, nu=0.25, rho=8e-6, Gc=2.7) + point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, b, :set_a) + point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, b, :set_b) + precrack!(b, :set_a, :set_b) + point_set!(p -> p[2] > l/2-Δx, b, :set_top) + point_set!(p -> p[2] < -l/2+Δx, b, :set_bottom) + velocity_bc!(t -> -30, b, :set_bottom, :y) + velocity_bc!(t -> 30, b, :set_top, :y) + vv = VelocityVerlet(steps=100) + job = Job(b, vv; path=path, freq=50) + submit(job) + return nothing + end + sim_nosb(30, path_threads) + + mpi_cmd = """ + using Peridynamics + function sim_nosb(N::Int, path::String) + l, Δx, δ, a = 1.0, 1/N, 3.015/N, 0.5 + pos, vol = uniform_box(l, l, 0.1l, Δx) + ids = sortperm(pos[2,:]) + b = Body(NOSBMaterial(), pos[:, ids], vol[ids]) + material!(b; horizon=3.015Δx, E=2.1e5, nu=0.25, rho=8e-6, Gc=2.7) + point_set!(p -> p[1] ≤ -l/2+a && 0 ≤ p[2] ≤ 2δ, b, :set_a) + point_set!(p -> p[1] ≤ -l/2+a && -2δ ≤ p[2] < 0, b, :set_b) + precrack!(b, :set_a, :set_b) + point_set!(p -> p[2] > l/2-Δx, b, :set_top) + point_set!(p -> p[2] < -l/2+Δx, b, :set_bottom) + velocity_bc!(t -> -30, b, :set_bottom, :y) + velocity_bc!(t -> 30, b, :set_top, :y) + vv = VelocityVerlet(steps=100) + job = Job(b, vv; path=path, freq=50) + submit(job) + return nothing + end + sim_nosb(30, "$path_mpi") + """ + run(`$(Peridynamics.MPI.mpiexec()) -n 2 $(Base.julia_cmd()) --project -e $(mpi_cmd)`) + + @test isdir(path_threads_vtk) + @test isdir(path_mpi_vtk) + vtk_files_threads = Peridynamics.find_vtk_files(path_threads_vtk) + vtk_files_mpi = Peridynamics.find_vtk_files(path_mpi_vtk) + @test length(vtk_files_mpi) == length(vtk_files_threads) == 3 + for i in eachindex(vtk_files_threads, vtk_files_mpi) + res_threads = read_vtk(vtk_files_threads[i]) + res_mpi = read_vtk(vtk_files_mpi[i]) + for key in keys(res_threads) + @test res_threads[key] ≈ res_mpi[key] + end + end + + rm(root; recursive=true, force=true) +end