diff --git a/docs/src/DomainBuffers/Setup.md b/docs/src/DomainBuffers/Setup.md index 3c63ed2b..fc1187d0 100644 --- a/docs/src/DomainBuffers/Setup.md +++ b/docs/src/DomainBuffers/Setup.md @@ -26,14 +26,26 @@ set_time_increment!(::FerriteAssembly.DomainBuffers, ::Any) ## Coupled simulations The `Simulation` type contains an abstract domain buffer, along with (optionally) the global degree of freedom values, which are used to get the local values for each item. -The main purpose is to conveniently collect these when passing into [`work!`](@ref), -especially in the case of `CoupledSimulations`. -The idea behind the coupled simulation setup is to give access to values from a different simulation -at the item level. For example, when solving two separate problems in parallel, and using staggered -iterations. See the [Phase-field fracture tutorial](@ref Phase-field-fracture) for an example. +A [`CoupledSimulations`](@ref) group is built, once, from a set of named `Simulation`s, and +gives access to values from other simulations at the item level (e.g. state variables and +local dof-values) via [`get_coupled_buffer`](@ref). For example, when solving two separate +problems in parallel using staggered iterations. See the +[Phase-field fracture tutorial](@ref Phase-field-fracture) for an example. Coupling is +resolved entirely at group-construction time; `work!`ing a group member (`work!(worker, +group.member_name)`) never re-discovers or rebuilds the coupling. + +!!! warning "Concurrency contract" + Coupled buffers reference the partner's *actual* mutable storage — no copies are made. + `work!` calls that share any of that storage must therefore not run concurrently with + each other. This includes: working two members of the *same* group at the same time; + working a member of a group at the same time as its own original (pre-group) source + `Simulation`; working members of *two different* groups that were built from the same + source `Simulation`(s) (e.g. a group and a later `replace_material`-built group that + still shares some members' storage by reference); and re-entrant `work!` calls that + would reuse the same scratch. Ordinary staggered iteration — working one member, then + another, in sequence — is safe; it is *simultaneous* access to shared scratch that is not. ```@docs Simulation -couple_buffers CoupledSimulations ``` diff --git a/docs/src/literate_tutorials/phasefield_fracture.jl b/docs/src/literate_tutorials/phasefield_fracture.jl index adfc33b1..780041ec 100644 --- a/docs/src/literate_tutorials/phasefield_fracture.jl +++ b/docs/src/literate_tutorials/phasefield_fracture.jl @@ -158,8 +158,11 @@ db_d_uc, Kd, rd, ndofs_d = setup(PhaseFieldFracture{:d}(mbase), grid, :d; ip_quad = Lagrange{RefQuadrilateral, 2}() ) -sim_u = Simulation(couple_buffers(db_u_uc; d = db_d_uc), zeros(ndofs_u), zeros(ndofs_u)) -sim_d = Simulation(couple_buffers(db_d_uc; u = db_u_uc), zeros(ndofs_d), zeros(ndofs_d)); +g = CoupledSimulations(( + u = Simulation(db_u_uc, zeros(ndofs_u), zeros(ndofs_u)), + d = Simulation(db_d_uc, zeros(ndofs_d), zeros(ndofs_d)), + )) +sim_u, sim_d = g.u, g.d; # Setup loading and boundary conditions load_function(t) = 1e-4 * t @@ -177,8 +180,8 @@ function get_reaction_dofs(dh) end; # ## Solving -# Write function to solve one simulation part, given the other as input. -function solve_single_part(sim, coupled, K, r, ch; firsttol = 1e-5, tol = 1e-6, maxiter = 100) +# Common function to solve each simulation +function solve_single_part(sim, K, r, ch; firsttol = 1e-5, tol = 1e-6, maxiter = 100) if ch !== nothing # Displacement part reaction_dofs = get_reaction_dofs(FerriteAssembly.get_dofhandler(sim)) else @@ -186,7 +189,7 @@ function solve_single_part(sim, coupled, K, r, ch; firsttol = 1e-5, tol = 1e-6, end for i in 1:maxiter assembler = start_assemble(K, r) - work!(assembler, sim, coupled) + work!(assembler, sim) rf = sum(i -> r[i], reaction_dofs; init = zero(eltype(r))) ch === nothing || apply_zero!(K, r, ch) res = norm(r) @@ -212,9 +215,9 @@ function solve(sim_u, sim_d, Ku, ru, Kd, rd, ch_u, grid) max_staggered = 2500 for iter in 1:max_staggered num = iter - u_converged, rf = solve_single_part(sim_u, CoupledSimulations(d = sim_d), Ku, ru, ch_u) + u_converged, rf = solve_single_part(sim_u, Ku, ru, ch_u) u_converged && break # Displacement was converged without updating - d_converged, _ = solve_single_part(sim_d, CoupledSimulations(u = sim_u), Kd, rd, nothing) + d_converged, _ = solve_single_part(sim_d, Kd, rd, nothing) d_converged && break # Damage was converged without updating iter ≥ max_staggered && error("Did not converge in staggered iterations") end diff --git a/src/Autodiff/autodiff.jl b/src/Autodiff/autodiff.jl index 94a827a8..8845c1ee 100644 --- a/src/Autodiff/autodiff.jl +++ b/src/Autodiff/autodiff.jl @@ -1,4 +1,6 @@ -mutable struct ElementResidual{S,M,CV,B<:CellBuffer} <: Function +const AnyCellBuffer = Union{CellBuffer,CoupledCellBuffer} + +mutable struct ElementResidual{S,M,CV,B<:AnyCellBuffer} <: Function state::S material::M cellvalues::CV @@ -23,19 +25,19 @@ function create_jacobian_config(er::ElementResidual) return ForwardDiff.JacobianConfig(er, re, ae, ForwardDiff.Chunk{length(ae)}()) end -struct AutoDiffCellBuffer{CB<:CellBuffer,ER<:ElementResidual,JC} <: AbstractCellBuffer +struct AutoDiffCellBuffer{CB<:AnyCellBuffer,ER<:ElementResidual,JC} <: AbstractCellBuffer cb::CB er::ER cfg::JC # JacobianConfig end -include("autodiff_unwrap.jl") # Experimental feature, include to remove large docstring from src here. +include("autodiff_unwrap.jl") # Experimental feature, include to remove large docstring from src here. """ - AutoDiffCellBuffer(cb::CellBuffer) + AutoDiffCellBuffer(cb::Union{CellBuffer,CoupledCellBuffer}) """ -function AutoDiffCellBuffer(cb::CellBuffer) +function AutoDiffCellBuffer(cb::AnyCellBuffer) cellstate = deepcopy(get_old_state(cb)) # to be safe, copy shouldn't be required. material = unwrap_material_for_ad(get_material(cb)) cellvalues = get_values(cb) @@ -55,7 +57,7 @@ reinit_buffer!(cb::AutoDiffCellBuffer, args...; kwargs...) = reinit_buffer!(cb.c set_time_increment!(c::AutoDiffCellBuffer, Δt) = set_time_increment!(c.cb, Δt) function _replace_material_with(ad_cb::AutoDiffCellBuffer{CB}, new_material) where CB - cb = setproperties(ad_cb.cb; material = new_material) + cb = _replace_material_with(ad_cb.cb, new_material) if isa(cb, CB) # If type didn't change, no need to recalculate autodiff buffers return setproperties(ad_cb; cb) else @@ -63,10 +65,6 @@ function _replace_material_with(ad_cb::AutoDiffCellBuffer{CB}, new_material) whe end end -function couple_buffers(cb::AutoDiffCellBuffer; kwargs...) - return AutoDiffCellBuffer(couple_buffers(cb.cb; kwargs...)) -end - function create_local(c::AutoDiffCellBuffer) cb = create_local(c.cb) AutoDiffCellBuffer(cb, deepcopy(c.er), deepcopy(c.cfg)) @@ -93,7 +91,7 @@ end # Standard method if no AutoDiffCellBuffer is defined. Should be no need to use, but good to keep for # benchmarks if desired. -function element_routine_ad!(Ke, re, state, ae, material, cellvalues, buffer::CellBuffer) +function element_routine_ad!(Ke, re, state, ae, material, cellvalues, buffer::AnyCellBuffer) rf!(re_, ae_) = element_residual!(re_, state, ae_, material, cellvalues, buffer) try # Setting Chunk explicitly to solve https://github.com/KnutAM/FerriteAssembly.jl/issues/9 diff --git a/src/Coupling.jl b/src/Coupling.jl new file mode 100644 index 00000000..29f6a96b --- /dev/null +++ b/src/Coupling.jl @@ -0,0 +1,272 @@ +# CoupledSimulations group construction and the cell-buffer-dependent parts of coupling. +# +# CoupledSimulation itself (struct, accessors, iteration, the scatter! overloads that push a +# threaded partner's current base state into its task-locals before work!) lives in +# Simulation.jl since it needs no cell-buffer types. Everything here references CellBuffer, +# CoupledCellBuffer, or AutoDiffCellBuffer, so this file must be included after those exist +# (see FerriteAssembly.jl's include order). + +""" + reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum::Int) + +Reinitialize the reader's own `cb.primary` against `sim.sim`, then reinitialize each partner +buffer in `cb.partner_buffers` against its own partner `Simulation` stored in `sim.partners`. +Partner reinitialization does not recurse: partner buffers are plain `CellBuffer`s, so no +further coupling initialization happens. +""" +function reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum::Int) + reinit_buffer!(cb.primary, sim.sim, cellnum) + map((b, s) -> (reinit_buffer!(b, s, cellnum); nothing), cb.partner_buffers, sim.partners) + return nothing +end + +struct CoupledSimulations{P<:NamedTuple, R<:NamedTuple, M<:NamedTuple} + primaries::P + refs::R + members::M +end + +""" + CoupledSimulations(primaries::NamedTuple; refs::NamedTuple = NamedTuple()) + +Build a group of mutually-wired simulations from `primaries` (members that read partners +and are worked via the group) and, optionally, `refs` (members with no outgoing +dependencies, still accessible/workable through the group but never rewired themselves). + +Each primary reads every other primary and every ref (excluded: itself). Names must be +unique across `primaries` and `refs`. Member access is direct/nonrecursive: `g.a`'s view of +`g.b` exposes `b`'s own local values, not `b`'s further coupling. + +```julia +g = CoupledSimulations((a = sima, b = simb, c = simc)) # mutual +g = CoupledSimulations((a = sima,); refs = (b = simb,)) # one-way: a reads b +g = CoupledSimulations((a = sima, b = simb); refs = (c = simc,)) # mixed + +work!(worker_a, g.a) +``` + +See the package documentation for the full setup-validation and replacement contract. +""" +function CoupledSimulations(primaries::NamedTuple; refs::NamedTuple = NamedTuple()) + isempty(primaries) && throw(ArgumentError("`primaries` must be a nonempty named tuple of `Simulation`s")) + all(v -> v isa Simulation, primaries) || throw(ArgumentError("`primaries` values must be `Simulation`s")) + all(v -> v isa Simulation, refs) || throw(ArgumentError("`refs` values must be `Simulation`s")) + overlap = intersect(keys(primaries), keys(refs)) + isempty(overlap) || throw(ArgumentError("primary and ref names must be unique, overlap: $overlap")) + reserved = intersect(union(keys(primaries), keys(refs)), (:primaries, :refs, :members)) + isempty(reserved) || throw(ArgumentError( + "member name(s) $reserved are reserved and would shadow `CoupledSimulations` internals")) + + all_members = merge(primaries, refs) + validate_storage_identity(all_members) + validate_task_counts_positive(all_members) + members = NamedTuple{keys(primaries)}( + Tuple(build_coupled_simulation(name, sim, all_members) for (name, sim) in pairs(primaries)) + ) + return CoupledSimulations(primaries, refs, members) +end + +function Base.getproperty(cs::CoupledSimulations, name::Symbol) + name in (:primaries, :refs, :members) && return getfield(cs, name) + members = getfield(cs, :members) + haskey(members, name) && return members[name] + refs = getfield(cs, :refs) + haskey(refs, name) && return refs[name] + throw(ArgumentError("CoupledSimulations has no member named `$name`")) +end + +Base.propertynames(cs::CoupledSimulations) = (:primaries, :refs, :members, keys(getfield(cs, :primaries))..., keys(getfield(cs, :refs))...) + +unwrap_cb(cb::CellBuffer) = cb +unwrap_cb(ad::AutoDiffCellBuffer) = ad.cb +unwrap_cb(ib) = throw(ArgumentError( + "coupling only supports `CellBuffer`/autodiff cell buffers, got $(typeof(ib))")) + +select_partner(p::TaskLocals, i::Int) = get_local(p, i) +select_partner(p, ::Int) = p + +_is_autodiff(ib::AutoDiffCellBuffer) = true +_is_autodiff(ib) = ib isa TaskLocals && get_base(ib) isa AutoDiffCellBuffer + +function build_coupled_itembuffer(reader_ibuf, partner_containers::NamedTuple) + autodiff = _is_autodiff(reader_ibuf) + wrap(primary_cb, partners_nt) = autodiff ? + AutoDiffCellBuffer(CoupledCellBuffer(primary_cb, partners_nt)) : + CoupledCellBuffer(primary_cb, partners_nt) + if reader_ibuf isa TaskLocals + n = length(get_locals(reader_ibuf)) + base = wrap(unwrap_cb(get_base(reader_ibuf)), map(unwrap_cb ∘ get_base, partner_containers)) + locals = [wrap(unwrap_cb(get_local(reader_ibuf, i)), + map(c -> unwrap_cb(select_partner(c, i)), partner_containers)) for i in 1:n] + return TaskLocals(base, locals) + else + return wrap(unwrap_cb(reader_ibuf), map(unwrap_cb ∘ get_base, partner_containers)) + end +end + +function validate_domain_pair(reader_db::AbstractDomainBuffer, partner_db::AbstractDomainBuffer, partner_name::Symbol) + get_grid(reader_db) === get_grid(partner_db) || throw(ArgumentError( + "coupling partner `$partner_name` uses a different grid than the reader")) + for (role, db) in ((:reader, reader_db), (Symbol(partner_name), partner_db)) + ib = get_base(get_itembuffer(db)) + (ib isa CellBuffer || ib isa AutoDiffCellBuffer) || throw(ArgumentError( + "coupling only supports `CellBuffer`/autodiff cell buffers, got $(typeof(ib)) for `$role`")) + end + issubset(getset(reader_db), getset(partner_db)) || throw(ArgumentError( + "coupling partner `$partner_name` does not cover all cells read by the reader")) + reader_threaded = reader_db isa ThreadedDomainBuffer + if reader_threaded + reader_tasks = get_num_tasks(reader_db) + reader_tasks > 0 || throw(ArgumentError("task count must be positive")) + partner_tasks = partner_db isa ThreadedDomainBuffer ? get_num_tasks(partner_db) : 1 + reader_tasks == partner_tasks || throw(ArgumentError( + "threaded reader with $reader_tasks tasks requires coupling partner `$partner_name` to provide " * + "$reader_tasks task-local buffers (a sequential partner counts as 1 slot); got $partner_tasks")) + end + return nothing +end + +# Returns (new_db, partner_sims::NamedTuple): the rebuilt domain buffer with coupled +# itembuffer(s), and the resolved per-domain partner `Simulation`s (for the caller to store +# on the owning `CoupledSimulation`, not on the itembuffer itself). +function build_coupled_domain(reader_db::AbstractDomainBuffer, partners::NamedTuple) + for (pname, p) in pairs(partners) + validate_domain_pair(reader_db, p.db, pname) + end + reader_ibuf = get_itembuffer(reader_db) + partner_containers = map(p -> get_itembuffer(p.db), partners) + coupled_ibuf = build_coupled_itembuffer(reader_ibuf, partner_containers) + new_db = setproperties(reader_db; itembuffer = coupled_ibuf) + partner_sims = map(p -> p.sim, partners) + return new_db, partner_sims +end + +# Resolve, for a single reader domain (named `dname` when the reader is a `Dict`, or +# `nothing` for a single-domain reader), the single-domain `Simulation` of a partner +# (sharing the partner's own `a`/`aold`). Errors if a partner does not provide a required +# domain, or if reader/partner shapes are mixed. +function partner_domain_sim(dname::Union{Nothing,String}, partner_name::Symbol, partner_sim::Simulation) + pdb = partner_sim.db + if dname === nothing + pdb isa DomainBuffers && throw(ArgumentError( + "mixed single-domain/dictionary coupling is not supported (reader is single-domain, " * + "partner `$partner_name` is a domain dictionary)")) + return partner_sim + else + pdb isa DomainBuffers || throw(ArgumentError( + "mixed single-domain/dictionary coupling is not supported (reader is a domain dictionary, " * + "partner `$partner_name` is single-domain)")) + haskey(pdb, dname) || throw(ArgumentError( + "coupling partner `$partner_name` does not supply required domain \"$dname\"")) + return Simulation(pdb[dname], partner_sim.a, partner_sim.aold) + end +end + +function build_coupled_simulation(name::Symbol, reader_sim::Simulation, all_members::NamedTuple) + partner_names = Tuple(k for k in keys(all_members) if k != name) + partner_sims = NamedTuple{partner_names}(Tuple(all_members[k] for k in partner_names)) + reader_db = reader_sim.db + if reader_db isa DomainBuffers + if isempty(reader_db) + new_db = reader_db # nothing to couple; preserves the original (correctly-typed) empty Dict + partners_by_domain = Dict{String, NamedTuple}() + else + built = Any[] + partners_by_domain = Dict{String, Any}() + for (dname, rdb) in reader_db + partners = NamedTuple{partner_names}(Tuple( + let psim_dom = partner_domain_sim(dname, pname, psim) + (sim = psim_dom, db = psim_dom.db) + end for (pname, psim) in pairs(partner_sims) + )) + ndb, dpartner_sims = build_coupled_domain(rdb, partners) + push!(built, dname => ndb) + partners_by_domain[dname] = dpartner_sims + end + new_db = Dict(built...) # infers the narrowest common concrete value type, matching MultiDomain(Threaded)Sim dispatch + end + new_sim = Simulation(new_db, reader_sim.a, reader_sim.aold) + return CoupledSimulation(new_sim, partners_by_domain) + else + partners = NamedTuple{partner_names}(Tuple( + let psim_dom = partner_domain_sim(nothing, pname, psim) + (sim = psim_dom, db = psim_dom.db) + end for (pname, psim) in pairs(partner_sims) + )) + new_db, dpartner_sims = build_coupled_domain(reader_db, partners) + new_sim = Simulation(new_db, reader_sim.a, reader_sim.aold) + return CoupledSimulation(new_sim, dpartner_sims) + end +end + +_scratch_identity(cb::CellBuffer) = cb.ae # survives replace_material's setproperties (fields copied by reference) + +function _domain_entries(name::Symbol, sim::Simulation) + db = sim.db + db isa DomainBuffers && return [(name, dname, _scratch_identity(unwrap_cb(get_base(get_itembuffer(d))))) for (dname, d) in db] + return [(name, "", _scratch_identity(unwrap_cb(get_base(get_itembuffer(db)))))] +end + +function validate_storage_identity(all_members::NamedTuple) + entries = reduce(vcat, (_domain_entries(name, sim) for (name, sim) in pairs(all_members))) + for i in eachindex(entries), j in (i+1):length(entries) + if entries[i][3] === entries[j][3] + throw(ArgumentError( + "members `$(entries[i][1])` (domain \"$(entries[i][2])\") and `$(entries[j][1])` " * + "(domain \"$(entries[j][2])\") alias the same underlying item-buffer storage; " * + "each member must own distinct mutable scratch")) + end + end + return nothing +end + +# Every member's own domain(s) must have a positive task count, independent of whether that +# member is ever paired as a "reader" against a partner (e.g. a ref, or a sole primary with +# no partners, would otherwise never be checked). +function _validate_own_task_count(name::Symbol, db::AbstractDomainBuffer) + db isa ThreadedDomainBuffer || return nothing + get_num_tasks(db) > 0 || throw(ArgumentError("member `$name` has a nonpositive task count")) + return nothing +end + +function validate_task_counts_positive(all_members::NamedTuple) + for (name, sim) in pairs(all_members) + db = sim.db + if db isa DomainBuffers + for (_, d) in db + _validate_own_task_count(name, d) + end + else + _validate_own_task_count(name, db) + end + end + return nothing +end + +""" + replace_material(g::CoupledSimulations, member::Symbol, f; domain = nothing) + +Return a new `CoupledSimulations` group in which `member`'s material has been replaced by +`f` (applied as `f(old_material)`), either for the whole member (`domain = nothing`) or only +for the named domain of a multi-domain member. Rebuilds the whole group (rerunning +constructor validation and autodiff configuration construction); other members are reused by +reference. Previously obtained handles (from the old group) keep their prior configuration. +""" +function replace_material(g::CoupledSimulations, member::Symbol, f; domain::Union{Nothing,String} = nothing) + primaries = getfield(g, :primaries) + refs = getfield(g, :refs) + is_primary = haskey(primaries, member) + is_primary || haskey(refs, member) || throw(ArgumentError("unknown member `$member`")) + old_sim = is_primary ? primaries[member] : refs[member] + if domain === nothing + new_db = replace_material(old_sim.db, f) + else + old_sim.db isa DomainBuffers || throw(ArgumentError( + "`domain` selector requires member `$member` to be a domain dictionary")) + new_db = replace_material(old_sim.db, domain, f) + end + new_sim = Simulation(new_db, old_sim.a, old_sim.aold) + new_primaries = is_primary ? merge(primaries, NamedTuple{(member,)}((new_sim,))) : primaries + new_refs = is_primary ? refs : merge(refs, NamedTuple{(member,)}((new_sim,))) + return CoupledSimulations(new_primaries; refs = new_refs) +end diff --git a/src/DomainBuffers.jl b/src/DomainBuffers.jl index 4d8e332f..d426608d 100644 --- a/src/DomainBuffers.jl +++ b/src/DomainBuffers.jl @@ -154,22 +154,15 @@ function replace_material(dbs::DomainBuffers, replacement_function) end """ - couple_buffers(dbs::Dict{String, <:AbstractDomainBuffer}; kwargs::Dict{String, <:AbstractDomainBuffer}...) - couple_buffers(db::AbstractDomainBuffer; kwargs::AbstractDomainBuffer...) + replace_material(dbs::Dict{String,AbstractDomainBuffer}, domain::String, replacement_function) -Return new buffer(s) that are coupled with the buffers provided as keyword arguments. The key is used in -[`get_coupled_buffer`](@ref) to get the coupled itembuffer, such that its values may be queried. - -!!! note - This functionality assumes that each setup has the same grid, and in case of multiple domains, these should also - match. +Return a new instance of `dbs` where as much as possible is copied by reference, and +where the material, `m`, of `dbs[domain]` is replaced by `replacement_function(m)`. +Other domains are copied by reference, unchanged. """ -function couple_buffers(dbs::DomainBuffers; kwargs...) - return Dict( - key => (all(haskey(v, key) for (_, v) in kwargs) ? - couple_buffers(db; (k => v[key] for (k, v) in kwargs)...) : - db) for (key, db) in dbs) - #return Dict(key => couple_buffers(db; (k => v[key] for (k, v) in kwargs)...) for (key, db) in dbs) +function replace_material(dbs::DomainBuffers, domain::String, replacement_function) + haskey(dbs, domain) || throw(ArgumentError("domain \"$domain\" not found in $(collect(keys(dbs)))")) + return Dict(key => (key == domain ? replace_material(db, replacement_function) : db) for (key, db) in dbs) end """ @@ -243,20 +236,7 @@ function replace_material(db::ThreadedDomainBuffer, replacement_function) return setproperties(db; itembuffer = TaskLocals(base_ibuf, task_ibuf)) end -function couple_buffers(db::DomainBuffer; kwargs...) - itembuffer = couple_buffers(db.itembuffer; (k => v.itembuffer for (k, v) in kwargs)...) - return setproperties(db; itembuffer) -end - -function couple_buffers(db::ThreadedDomainBuffer; kwargs...) - base_ibuf = couple_buffers(get_base(db.itembuffer); (k => get_base(v.itembuffer) for (k, v) in kwargs)...) - task_ibuf = map(enumerate(get_locals(db.itembuffer))) do (i, ibuf) - couple_buffers(ibuf; (k => get_local(v.itembuffer, i) for (k, v) in kwargs)...) - end - return setproperties(db; itembuffer = TaskLocals(base_ibuf, task_ibuf)) -end - -# Experimental: Insert new states, allows reusing the buffer for multiple simulations with same +# Experimental: Insert new states, allows reusing the buffer for multiple simulations with same # initial state (grid, dh, etc.), but which experience different loading. Typically for RVE simulations. function replace_states!(dbs::Dict{String, <:AbstractDomainBuffer}, states::Dict{String, <:StateVariables}) keys(dbs) == keys(states) || throw(ArgumentError("keys of dictionaries don't match")) diff --git a/src/FerriteAssembly.jl b/src/FerriteAssembly.jl index 0e706ea4..c4618898 100644 --- a/src/FerriteAssembly.jl +++ b/src/FerriteAssembly.jl @@ -17,9 +17,11 @@ include("Simulation.jl") include("setup.jl") include("ItemBuffers/CellBuffer.jl") +include("ItemBuffers/CoupledCellBuffer.jl") include("ItemBuffers/FacetBuffer.jl") include("Autodiff/autodiff.jl") +include("Coupling.jl") include("work.jl") include("Workers/Assemblers.jl") include("Workers/Integrators.jl") @@ -29,7 +31,7 @@ include("LoadHandler/LoadHandler.jl") # Setup export DomainSpec, setup_domainbuffer, setup_domainbuffers -export Simulation, CoupledSimulations, couple_buffers +export Simulation, CoupledSimulations # Main functions to use during simulations export work!, update_states!, revert_states!, set_time_increment! # Workers diff --git a/src/ItemBuffers/AbstractItemBuffer.jl b/src/ItemBuffers/AbstractItemBuffer.jl index ab76eaad..5cf3692d 100644 --- a/src/ItemBuffers/AbstractItemBuffer.jl +++ b/src/ItemBuffers/AbstractItemBuffer.jl @@ -58,10 +58,10 @@ function get_user_cache end """ get_coupled_buffer(b::AbstractItemBuffer, key::Symbol) -Get the coupled buffer `key` from `b`. To enable this, use [`couple_buffers`](@ref) on the -domain buffers. The coupled buffer can be queried just like a normal item buffer, -e.g. by calling `get_state(coupled_buffer)`. -""" +Get the coupled buffer `key` from `b`. To enable this, build a [`CoupledSimulations`](@ref) +group and `work!` its member handles. The coupled buffer can be queried just like a normal +item buffer, e.g. by calling `get_state(coupled_buffer)`. +""" @inline get_coupled_buffer(b::AbstractItemBuffer, key::Symbol) = getfield(get_coupled_buffers(b), key) """ diff --git a/src/ItemBuffers/CellBuffer.jl b/src/ItemBuffers/CellBuffer.jl index 211df395..af066946 100644 --- a/src/ItemBuffers/CellBuffer.jl +++ b/src/ItemBuffers/CellBuffer.jl @@ -9,10 +9,10 @@ Each worker that supports a cellbuffer should overload this function. """ function work_single_cell! end -mutable struct CellBuffer{T,CC,CV,DR,MT,ST,UD,UC,CB} <: AbstractCellBuffer +mutable struct CellBuffer{T,CC,CV,DR,MT,ST,UD,UC} <: AbstractCellBuffer const ae_old::Vector{T} # Old element dof values const ae::Vector{T} # Current element dof values - const re::Vector{T} # Residual/force vector + const re::Vector{T} # Residual/force vector const Ke::Matrix{T} # Element stiffness matrix const dofs::Vector{Int} # celldofs const coords::CC # cellcoords (or what is required to reinit cellvalues) @@ -26,7 +26,6 @@ mutable struct CellBuffer{T,CC,CV,DR,MT,ST,UD,UC,CB} <: AbstractCellBuffer old_state::ST # Old state variables for the cell (updated in reinit!) const user_data::UD # User data for the cell (used for additional information) const user_cache::UC # Cache for the cell (user type) (deepcopy for each thread) - const coupled_buffers::CB # nothing or NamedTuple with staggered coupled `CellBuffer`s. end """ @@ -49,9 +48,9 @@ function CellBuffer(numdofs::Int, coords, cellvalues, material, state, dofrange, cellid = -1 cache = allocate_cell_cache(material, cellvalues) return CellBuffer( - zeros(numdofs), zeros(numdofs), zeros(numdofs), zeros(numdofs,numdofs), - zeros(Int, numdofs), coords, - cellvalues, Δt, cellid, dofrange, material, state, state, user_data, cache, nothing) + zeros(numdofs), zeros(numdofs), zeros(numdofs), zeros(numdofs,numdofs), + zeros(Int, numdofs), coords, + cellvalues, Δt, cellid, dofrange, material, state, state, user_data, cache) end setup_cellbuffer(ad::Bool, args...; kwargs...) = setup_cellbuffer(Val(ad), args...; kwargs...) @@ -61,10 +60,6 @@ function setup_cellbuffer(::Val{false}, sdh, cv, material, cell_state, dofrange, return CellBuffer(numdofs, coords, cv, material, cell_state, dofrange, user_data) end -function couple_buffers(cb::CellBuffer; kwargs...) - return setproperties(cb; coupled_buffers = NamedTuple{keys(kwargs)}(values(kwargs))) -end - function setup_cellbuffer(::Val{true}, args...) return AutoDiffCellBuffer(setup_cellbuffer(Val(false), args...)) end @@ -73,7 +68,7 @@ end # TaskLocals interface (only `create_local` required for other `AbstractCellBuffer`s) (unless gather! is req.) function create_local(cb::CellBuffer) dcpy = map(deepcopy, (cb.ae_old, cb.ae, cb.re, cb.Ke, cb.dofs, cb.coords, cb.cellvalues, cb.Δt, cb.cellid, cb.dofrange, cb.material, cb.state, cb.old_state)) - return CellBuffer(dcpy..., cb.user_data, deepcopy(cb.user_cache), create_local(cb.coupled_buffers)) + return CellBuffer(dcpy..., cb.user_data, deepcopy(cb.user_cache)) end set_time_increment!(cb::CellBuffer, Δt) = (cb.Δt=Δt) @@ -106,8 +101,6 @@ Ferrite.getfieldnames(cb::CellBuffer) = keys(cb.dofrange) @inline get_user_cache(cb::CellBuffer) = cb.user_cache -@inline get_coupled_buffers(cb::CellBuffer) = cb.coupled_buffers - """ FerriteAssembly.allocate_cell_cache(material, cellvalues) @@ -118,15 +111,15 @@ used to reduce allocations. Returns `nothing` by default. allocate_cell_cache(::Any, ::Any) = nothing """ - reinit_buffer!(cb::CellBuffer, sim::Simulation, coupled, cellnum::Int) + reinit_buffer!(cb::CellBuffer, sim::Simulation, cellnum::Int) Reinitialize the `cb::CellBuffer` for cell number `cellnum`. The global degree of freedom vectors `a` (current) and `aold` are used to update the cell degree of freedom vectors in `c`. If the global vectors are not included in `sim`, the corresponding local vectors are set to `NaN` -The element stiffness, `cb.Ke`, and residual, `cb.re`, are also zeroed. +The element stiffness, `cb.Ke`, and residual, `cb.re`, are also zeroed. """ -function reinit_buffer!(cb::CellBuffer, sim::Simulation, coupled, cellnum::Int) +function reinit_buffer!(cb::CellBuffer, sim::Simulation, cellnum::Int) dh = get_dofhandler(sim) grid = dh.grid cb.cellid = cellnum @@ -139,19 +132,7 @@ function reinit_buffer!(cb::CellBuffer, sim::Simulation, coupled, cellnum::Int) _copydofs!(cb.ae_old, sim.aold, cb.dofs) # ae_old .= a_old[dofs] fill!(cb.Ke, 0) fill!(cb.re, 0) - reinit_coupled!(cb.coupled_buffers, coupled, cellnum) - return nothing # Ferrite's reinit! doesn't return -end - -# No coupled buffer, no coupled simulation -reinit_coupled!(::Nothing, coupled::CoupledSimulations{@NamedTuple{}}, cellnum::Int) = nothing - -function reinit_coupled!(coupled_buffers::NamedTuple, coupled::CoupledSimulations, cellnum::Int) - if length(coupled_buffers) != length(coupled.sims) - throw(ArgumentError("When using coupled simulations, the coupled buffers must match the coupled simulations")) - end - tuple((reinit_buffer!(cb, coupled.sims[k], CoupledSimulations(), cellnum) for (k, cb) in pairs(coupled_buffers))...) - return nothing + return nothing # Ferrite's reinit! doesn't return end function _replace_material_with(cb::CellBuffer, new_material) diff --git a/src/ItemBuffers/CoupledCellBuffer.jl b/src/ItemBuffers/CoupledCellBuffer.jl new file mode 100644 index 00000000..f0ee9fd2 --- /dev/null +++ b/src/ItemBuffers/CoupledCellBuffer.jl @@ -0,0 +1,42 @@ +""" + CoupledCellBuffer(primary::CellBuffer, partner_buffers::NamedTuple) + +Wraps a reader's own `primary::CellBuffer` together with references to its coupling +partners' plain `CellBuffer`s (`partner_buffers`, returned by [`get_coupled_buffer`](@ref)). + +`partner_buffers` are never themselves `CoupledCellBuffer`s or `AutoDiffCellBuffer`s: +coupling is direct and nonrecursive, so a partner's own coupling (if any) is not exposed. + +Constructed once per task at [`CoupledSimulations`](@ref) setup time; never rebuilt per cell +or per `work!` call. Does *not* hold the partner `Simulation`s needed to reinitialize +`partner_buffers`: those live once on the [`CoupledSimulation`](@ref) passed into +[`reinit_buffer!`](@ref) (defined in `Coupling.jl`, once that type exists), rather than being +duplicated into every task-local copy of this buffer. +""" +struct CoupledCellBuffer{CB<:CellBuffer, PB<:NamedTuple} <: AbstractCellBuffer + primary::CB + partner_buffers::PB +end + +for op = (:get_Ke, :get_re, :get_ae, :get_material, :get_values, :get_time_increment, + :get_aeold, :get_state, :get_old_state, :get_user_data, :get_user_cache) + eval(quote + @inline $op(cb::CoupledCellBuffer) = $op(cb.primary) + end) +end + +get_coupled_buffers(cb::CoupledCellBuffer) = cb.partner_buffers + +set_time_increment!(cb::CoupledCellBuffer, Δt) = set_time_increment!(cb.primary, Δt) + +for op = (:celldofs, :getcoordinates, :getfieldnames, :cellid) + eval(quote + Ferrite.$op(cb::CoupledCellBuffer, args...) = Ferrite.$op(cb.primary, args...) + end) +end +Ferrite.dof_range(cb::CoupledCellBuffer, name::Symbol) = Ferrite.dof_range(cb.primary, name) + +function _replace_material_with(cb::CoupledCellBuffer, new_material) + new_primary = _replace_material_with(cb.primary, new_material) + return CoupledCellBuffer(new_primary, cb.partner_buffers) +end diff --git a/src/ItemBuffers/FacetBuffer.jl b/src/ItemBuffers/FacetBuffer.jl index 115cfb15..17adebea 100644 --- a/src/ItemBuffers/FacetBuffer.jl +++ b/src/ItemBuffers/FacetBuffer.jl @@ -100,7 +100,7 @@ allocations. Returns `nothing` by default. """ allocate_facet_cache(::Any, ::Any) = nothing -function reinit_buffer!(fb::FacetBuffer, sim::Simulation, #=coupled=#_, fi::FacetIndex) +function reinit_buffer!(fb::FacetBuffer, sim::Simulation, fi::FacetIndex) cellnum, facetnr = fi dh = get_dofhandler(sim) fb.cellid = cellnum diff --git a/src/Simulation.jl b/src/Simulation.jl index 7d8f7350..13f42f44 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -1,76 +1,112 @@ +abstract type AbstractSimulation{DB} end + +const AbstractSingleDomainSim = AbstractSimulation{<:DomainBuffer} +const AbstractMultiDomainSim = AbstractSimulation{<:Dict{String, <:DomainBuffer}} +const AbstractSingleDomainThreadedSim = AbstractSimulation{<:ThreadedDomainBuffer} +const AbstractMultiDomainThreadedSim = AbstractSimulation{<:Dict{String, <:ThreadedDomainBuffer}} + +# Must be defined +""" + get_domainbuffer(sim::AbstractSimulation) + +Accessor for the automatic forwarding for domainbuffer methods to work +""" +function get_domainbuffer end + +# Forwarding for public API +get_material(sim::AbstractSimulation, args::Vararg{Any, N}) where N = get_material(get_domainbuffer(sim), args...) +get_dofhandler(sim::AbstractSimulation) = get_dofhandler(get_domainbuffer(sim)) +get_grid(sim::AbstractSimulation) = get_grid(get_domainbuffer(sim)) +get_state(sim::AbstractSimulation, args::Vararg{Any, N}) where N = get_state(get_domainbuffer(sim), args...) +get_old_state(sim::AbstractSimulation, args::Vararg{Any, N}) where N = get_old_state(get_domainbuffer(sim), args...) +getset(sim::AbstractSimulation, args::Vararg{Any, N}) where N = getset(get_domainbuffer(sim), args...) +update_states!(sim::AbstractSimulation; kwargs...) = update_states!(get_domainbuffer(sim); kwargs...) +set_time_increment!(sim::AbstractSimulation, Δt) = set_time_increment!(get_domainbuffer(sim), Δt) +revert_states!(sim::AbstractSimulation) = revert_states!(get_domainbuffer(sim)) + +# Forwarding for internal API +get_num_tasks(sim::AbstractSimulation) = get_num_tasks(get_domainbuffer(sim)) +get_chunks(sim::AbstractSimulation{<:AbstractDomainBuffer}) = get_chunks(get_domainbuffer(sim)) +get_itembuffer(sim::AbstractSimulation, args::Vararg{Any, N}) where {N} = get_itembuffer(get_domainbuffer(sim), args...) + + """ Simulation(db, a = nothing, aold = nothing) -A `Simulation` is a collection of the simulation domain(s) `db`, and the -global degree of freedom vectors, `a` and `aold`. +A `Simulation` is a collection of the simulation domain(s) `db`, and the +global degree of freedom vectors, `a` and `aold`. -**Note:** +**Note:** If `a` or `aold` are not provided, the local vectors will have `NaN` values. """ struct Simulation{ - DB <: Union{DomainBuffers, AbstractDomainBuffer}, - TA <: Union{Nothing, AbstractVector}, + DB <: Union{DomainBuffers, AbstractDomainBuffer}, + TA <: Union{Nothing, AbstractVector}, TAO <: Union{Nothing, AbstractVector} - } + } <: AbstractSimulation{DB} db::DB a::TA aold::TAO end Simulation(db::Union{DomainBuffers, AbstractDomainBuffer}, a = nothing, aold = nothing) = Simulation(db, a, aold) -const SingleDomainSim = Simulation{<:DomainBuffer} -const MultiDomainSim = Simulation{<:Dict{String, <:DomainBuffer}} -const SingleDomainThreadedSim = Simulation{<:ThreadedDomainBuffer} -const MultiDomainThreadedSim = Simulation{<:Dict{String, <:ThreadedDomainBuffer}} +get_domainbuffer(sim::Simulation) = sim.db -# Forwarding for public API -get_material(sim::Simulation, args::Vararg{Any, N}) where N = get_material(sim.db, args...) -get_dofhandler(sim::Simulation) = get_dofhandler(sim.db) -get_grid(sim::Simulation) = get_grid(sim.db) -get_state(sim::Simulation, args::Vararg{Any, N}) where N = get_state(sim.db, args...) -get_old_state(sim::Simulation, args::Vararg{Any, N}) where N = get_old_state(sim.db, args...) -getset(sim::Simulation, args::Vararg{Any, N}) where N = getset(sim.db, args...) -update_states!(sim::Simulation; kwargs...) = update_states!(sim.db; kwargs...) -set_time_increment!(sim::Simulation, Δt) = set_time_increment!(sim.db, Δt) -revert_states!(sim::Simulation) = revert_states!(sim.db) - -# Forwarding for internal API -get_num_tasks(sim::Simulation) = get_num_tasks(sim.db) -get_chunks(sim::Simulation{<:AbstractDomainBuffer}) = get_chunks(sim.db) -get_itembuffer(sim::Simulation, args::Vararg{Any, N}) where {N} = get_itembuffer(sim.db, args...) - -# Internal API -get_domain_simulation(sim::Simulation{<:DomainBuffers}, name::String) = Simulation(sim.db[name], sim.a, sim.aold) ## Iterator interface @inline function _iterate(sim::Simulation{<:DomainBuffers}, iter) iter === nothing && return nothing - ((name, db), state) = iter + ((name, db), state) = iter return ((name, Simulation(db, sim.a, sim.aold)), state) end Base.iterate(sim::Simulation{<:DomainBuffers}) = _iterate(sim, iterate(sim.db)) Base.iterate(sim::Simulation{<:DomainBuffers}, iter) = _iterate(sim, iterate(sim.db, iter)) +scatter!(::AbstractSimulation) = nothing # only a threaded sim has task-local buffers to scatter into +scatter!(sim::Simulation{<:ThreadedDomainBuffer}) = scatter!(get_itembuffer(sim)) + """ - CoupledSimulations(; key1 = sim1::Simulation, key2 = sim2::Simulation, ...) + CoupledSimulation(sim, partners) -Setup the collection of coupled simulations to allow values (such as state variables and -local dof-values from these simulations to be available when `work!`ing another simulation, -if the buffers have been coupled with [`couple_buffers`](@ref). -The coupled itembuffer on the local level is accessed with [`get_coupled_buffer`](@ref). +A handle to one primary member of a [`CoupledSimulations`](@ref) group (e.g. `group.a`). +`sim` is a [`Simulation`](@ref) whose domain buffer(s) have been rebuilt with coupled +itembuffers. `partners` holds the resolved partner `Simulation`s this member reads from: a +`NamedTuple{name}` of partner `Simulation`s for a single-domain member, or a +`Dict{String,<:NamedTuple}` (one `NamedTuple` of partner `Simulation`s per domain name) for a +multi-domain member. This is the single, canonical copy of that information — passed into +[`reinit_buffer!`](@ref) at call time rather than duplicated into every task-local buffer. + +Forwards the ordinary [`Simulation`](@ref) accessor API (`.a`, `.aold`, `.db`, +`get_dofhandler`, `get_state`, `set_time_increment!`, `update_states!`, etc.). """ -struct CoupledSimulations{NT <: NamedTuple{<:Any, <:NTuple{<:Any, Simulation}}} - sims::NT +struct CoupledSimulation{DB, S <: Simulation{DB}, P} <: AbstractSimulation{DB} + sim::S + partners::P +end + +function Base.getproperty(csim::CoupledSimulation, name::Symbol) + name === :sim && return getfield(csim, :sim) + name === :partners && return getfield(csim, :partners) + return getproperty(getfield(csim, :sim), name) +end + +# Include the forwarded Simulation properties (`.a`, `.aold`, `.db`) so they tab-complete. +Base.propertynames(csim::CoupledSimulation) = (:sim, :partners, propertynames(getfield(csim, :sim))...) + +get_domainbuffer(sim::CoupledSimulation) = get_domainbuffer(sim.sim) + +function scatter!(sim::CoupledSimulation{<:ThreadedDomainBuffer}) + scatter!(sim.sim) + map(scatter!, sim.partners) +end + +replace_material(::CoupledSimulation, args...; kwargs...) = throw(ArgumentError( + "replace_material on a CoupledSimulations member is not supported; use " * + "replace_material(group, member_name, f) to rebuild the whole group instead.")) + +@inline function _iterate(csim::CoupledSimulation{<:DomainBuffers}, iter) + iter === nothing && return nothing + ((name, sim), state) = iter + return ((name, CoupledSimulation(sim, csim.partners[name])), state) end -CoupledSimulations(; kwargs...) = CoupledSimulations(NamedTuple{keys(kwargs)}(values(kwargs))) - -function get_domain_simulation(cs::CoupledSimulations, name::String) - # Need to return a named tuple with only the simulations that have a domain called `name` - sims = Pair{Symbol, Simulation}[] - for (key, sim) in zip(keys(cs.sims), values(cs.sims)) - if haskey(sim.db, name) - push!(sims, key => get_domain_simulation(sim, name)) - end - end - return CoupledSimulations(NamedTuple(sims)) -# return CoupledSimulations(map(s -> get_domain_simulation(s, name), cs.sims)) -end \ No newline at end of file +Base.iterate(sim::CoupledSimulation{<:DomainBuffers}) = _iterate(sim, iterate(sim.sim)) +Base.iterate(sim::CoupledSimulation{<:DomainBuffers}, iter) = _iterate(sim, iterate(sim.sim, iter)) diff --git a/src/work.jl b/src/work.jl index cd8a6d0d..f58fbf37 100644 --- a/src/work.jl +++ b/src/work.jl @@ -1,107 +1,104 @@ -function work!(worker, buffer::Union{AbstractDomainBuffer, DomainBuffers}; a = nothing, aold = nothing) - return work!(worker, Simulation(buffer, a, aold)) -end - -""" - work!(worker, sim::Simulation, [coupled_simulations::CoupledSimulations]) - -Perform the work according to `worker` over the domain(s) in `sim`. - -**Advance usage:** By passing the optional `coupled_simulations`, values from those simulations -(e.g. state variables and local dof-values) become available on the local level via -[`get_coupled_buffer`](@ref). This requires that the domainbuffer(s) in `sim` has been coupled -using [`couple_buffers`](@ref). - - work!(worker, db::Union{AbstractDomainBuffer, Dict}; a = nothing, aold = nothing) - -Simplified interface that doesn't support coupled simulations, directly forwarded to -`work!(worker, Simulation(db, a, aold))`. The global degree of freedom vectors, `a` and `aold`, -make their corresponding local values available. If not passed, the local values are `NaN`s. -""" -function work!(worker, multisim::MultiDomainSim, coupled_simulations = CoupledSimulations()) - for (name, sim) in multisim - skip_this_domain(worker, name) && continue - coupled = get_domain_simulation(coupled_simulations, name) - work_domain_sequential!(worker, sim, coupled) - end -end -function work!(worker, sim::SingleDomainSim, coupled_simulations = CoupledSimulations()) - work_domain_sequential!(worker, sim, coupled_simulations) -end -function work!(worker, multisim::MultiDomainThreadedSim, coupled_simulations = CoupledSimulations()) - if can_thread(worker) - workers = TaskLocals(worker, num_tasks = get_num_tasks(multisim)) - for (name, sim) in multisim - skip_this_domain(worker, name) && continue - coupled = get_domain_simulation(coupled_simulations, name) - work_domain_threaded!(workers, sim, coupled) - end - else - for (name, sim) in multisim - skip_this_domain(worker, name) && continue - coupled = get_domain_simulation(coupled_simulations, name) - work_domain_sequential!(worker, sim, coupled) - end - end -end -function work!(worker, sim::SingleDomainThreadedSim, coupled_simulations = CoupledSimulations()) - if can_thread(worker) - workers = TaskLocals(worker; num_tasks = get_num_tasks(sim)) - work_domain_threaded!(workers, sim, coupled_simulations) - else - work_domain_sequential!(worker, sim, coupled_simulations) - end -end - -function work_domain_sequential!(worker, sim::Simulation{<:AbstractDomainBuffer}, coupled) - itembuffer = get_base(get_itembuffer(sim)) # get_base if threaded buffer - for itemnr in getset(sim) - reinit_buffer!(itembuffer, sim, coupled, itemnr) - work_single!(worker, itembuffer) - end -end - -function work_domain_threaded!(workers, sim::SingleDomainThreadedSim, coupled) - itembuffers = get_itembuffer(sim) #::TaskLocals - scatter!(itembuffers) - scatter!(workers) - num_tasks = get_num_tasks(sim) # Default to Threads.nthreads() - for chunk_vector in get_chunks(sim) - taskchunks = TaskChunks(chunk_vector) - Base.Experimental.@sync begin - for taskid in 1:num_tasks - itembuffer = get_local(itembuffers, taskid) - worker = get_local(workers, taskid) - Threads.@spawn begin - while true - taskchunk = get_chunk(taskchunks) # Union{Vector{Int}, Nothing} - taskchunk === nothing && break - for itemnr in taskchunk - reinit_buffer!(itembuffer, sim, coupled, itemnr) - work_single!(worker, itembuffer) - end # itemnr - end #chunk - end #spawn - end #taskid - end #sync - end #chunk_vectors - gather!(itembuffers) - gather!(workers) -end - -# Worker interface -""" - can_thread(worker)::Bool - -Does the worker support multithreaded work? Defaults to `false`. -If this returns `true`, the worker must support the `TaskLocals` interface. -""" -can_thread(::Any) = false - -""" - skip_this_domain(worker, name::String) - -Should the domain with key `name` be skipped during work? Defaults to `false`. -Can be used to e.g. only loop over parts of a domain. -""" -skip_this_domain(::Any, ::String) = false # opt-in to skip domains (used for integration) +function work!(worker, buffer::Union{AbstractDomainBuffer, DomainBuffers}; a = nothing, aold = nothing) + return work!(worker, Simulation(buffer, a, aold)) +end + +""" + work!(worker, sim::Simulation) + +Perform the work according to `worker` over the domain(s) in `sim`. + +**Coupled simulations:** To make values from other simulations (e.g. state variables and +local dof-values) available on the local level via [`get_coupled_buffer`](@ref), build a +[`CoupledSimulations`](@ref) group and call `work!(worker, group.member_name)` instead; +the member handle already carries its resolved coupling. + + work!(worker, db::Union{AbstractDomainBuffer, Dict}; a = nothing, aold = nothing) + +Simplified interface, directly forwarded to `work!(worker, Simulation(db, a, aold))`. +The global degree of freedom vectors, `a` and `aold`, make their corresponding local values +available. If not passed, the local values are `NaN`s. +""" +function work!(worker, multisim::AbstractMultiDomainSim) + for (name, sim) in multisim + skip_this_domain(worker, name) && continue + work_domain_sequential!(worker, sim) + end +end +function work!(worker, sim::AbstractSingleDomainSim) + work_domain_sequential!(worker, sim) +end +function work!(worker, multisim::AbstractMultiDomainThreadedSim) + if can_thread(worker) + workers = TaskLocals(worker, num_tasks = get_num_tasks(multisim)) + for (name, sim) in multisim + skip_this_domain(worker, name) && continue + work_domain_threaded!(workers, sim) + end + else + for (name, sim) in multisim + skip_this_domain(worker, name) && continue + work_domain_sequential!(worker, sim) + end + end +end +function work!(worker, sim::AbstractSingleDomainThreadedSim) + if can_thread(worker) + workers = TaskLocals(worker; num_tasks = get_num_tasks(sim)) + work_domain_threaded!(workers, sim) + else + work_domain_sequential!(worker, sim) + end +end + +function work_domain_sequential!(worker, sim::AbstractSimulation{<:AbstractDomainBuffer}) + itembuffer = get_base(get_itembuffer(sim)) # get_base if threaded buffer + for itemnr in getset(sim) + reinit_buffer!(itembuffer, sim, itemnr) + work_single!(worker, itembuffer) + end +end + +function work_domain_threaded!(workers, sim::AbstractSingleDomainThreadedSim) + scatter!(sim) # Includes scatter of the `itembuffers` + itembuffers = get_itembuffer(sim) #::TaskLocals + scatter!(workers) + num_tasks = get_num_tasks(sim) # Default to Threads.nthreads() + for chunk_vector in get_chunks(sim) + taskchunks = TaskChunks(chunk_vector) + Base.Experimental.@sync begin + for taskid in 1:num_tasks + itembuffer = get_local(itembuffers, taskid) + worker = get_local(workers, taskid) + Threads.@spawn begin + while true + taskchunk = get_chunk(taskchunks) # Union{Vector{Int}, Nothing} + taskchunk === nothing && break + for itemnr in taskchunk + reinit_buffer!(itembuffer, sim, itemnr) + work_single!(worker, itembuffer) + end # itemnr + end #chunk + end #spawn + end #taskid + end #sync + end #chunk_vectors + gather!(itembuffers) + gather!(workers) +end + +# Worker interface +""" + can_thread(worker)::Bool + +Does the worker support multithreaded work? Defaults to `false`. +If this returns `true`, the worker must support the `TaskLocals` interface. +""" +can_thread(::Any) = false + +""" + skip_this_domain(worker, name::String) + +Should the domain with key `name` be skipped during work? Defaults to `false`. +Can be used to e.g. only loop over parts of a domain. +""" +skip_this_domain(::Any, ::String) = false # opt-in to skip domains (used for integration) diff --git a/test/assemblers.jl b/test/assemblers.jl index 5117c390..d4b923ff 100644 --- a/test/assemblers.jl +++ b/test/assemblers.jl @@ -105,4 +105,13 @@ @test rb ≈ ra end end + + @testset "can_thread/skip_this_domain defaults" begin + # Every worker in the package overrides both traits, so a worker relying purely on + # the generic `::Any` fallback (e.g. a user-defined worker that never opts in to + # threading or domain skipping) must still get sensible defaults. + struct FA_DummyWorker end + @test !FA.can_thread(FA_DummyWorker()) + @test !FA.skip_this_domain(FA_DummyWorker(), "somedomain") + end end \ No newline at end of file diff --git a/test/coupled_simulations.jl b/test/coupled_simulations.jl new file mode 100644 index 00000000..b7848600 --- /dev/null +++ b/test/coupled_simulations.jl @@ -0,0 +1,498 @@ +@testset "CoupledSimulations" begin + # `@inferred` can't check constant propagation for property access (`.a` etc.): its macro + # only accepts call expressions, and it infers based on the *runtime type* of arguments + # passed to `getproperty`, not the literal property name baked into `x.a` syntax at the + # call site (the case that actually matters, since that's how these are used everywhere). + # `Base.return_types` on a closure containing the literal dot-access captures that. + is_concrete_inferred(f, argtypes...) = begin + rt = Base.return_types(f, argtypes) + length(rt) == 1 && isconcretetype(rt[1]) + end + + grid = generate_grid(Quadrilateral, (2,2)) + addcellset!(grid, "left", x -> x[1] < eps()) + addcellset!(grid, "right", setdiff(1:getncells(grid), getcellset(grid, "left"))) + ip = Lagrange{RefQuadrilateral,1}() + dh1 = close!(add!(DofHandler(grid), :u, ip)) + dh2 = close!(add!(DofHandler(grid), :v, ip^2)) + qr = QuadratureRule{RefQuadrilateral}(2) + cvu = CellValues(qr, ip, ip) + cvv = CellValues(qr, ip^2, ip) + + struct CS_MA end + struct CS_MB end + struct CS_MB2 end # distinct material type, used to exercise replace_material changing type + # aold will be same in both cases (for both components in the case of MB) + # a will be 3 times larger for first component in MB, and 5 times for second component + # State will be 6 times larger for MB, obtained by multiplying the function values by factor 2 + FerriteAssembly.create_cell_state(::CS_MA, cv, x, ae, args...) = [function_value(cv, i, ae) for i in 1:getnquadpoints(cv)] + FerriteAssembly.create_cell_state(::CS_MB, cv, x, ae, args...) = [2 * function_value(cv, i, ae)[1] for i in 1:getnquadpoints(cv)] + FerriteAssembly.create_cell_state(::CS_MB2, cv, x, ae, args...) = nothing + + # Set to a Float64 (not NaN) by the BUG-003 regression below to check that the partner's + # (`:b`'s) *own* current time increment is observed, independent of the reader's own Δt. + expected_b_dt = Ref(NaN) + expected_b_material = Ref{DataType}(CS_MB) + function FerriteAssembly.element_routine!(Ke, re, state, ae, m::CS_MA, cv, buffer) + cb_b = FerriteAssembly.get_coupled_buffer(buffer, :b) + @test FerriteAssembly.get_material(cb_b) isa expected_b_material[] + if expected_b_material[] === CS_MB + @test 3 * ae ≈ FerriteAssembly.get_ae(cb_b)[1:2:end] + @test 5 * ae ≈ FerriteAssembly.get_ae(cb_b)[2:2:end] + @test FerriteAssembly.get_aeold(buffer) ≈ FerriteAssembly.get_aeold(cb_b)[1:2:end] + @test FerriteAssembly.get_aeold(buffer) ≈ FerriteAssembly.get_aeold(cb_b)[2:2:end] + @test 6 * state ≈ FerriteAssembly.get_state(cb_b) + end + isnan(expected_b_dt[]) || @test FerriteAssembly.get_time_increment(cb_b) == expected_b_dt[] + # Present only in the 3-member mutual-coupling test below; checks that :b and :c are + # not positionally swapped when a reader has two distinct partners. + if haskey(FerriteAssembly.get_coupled_buffers(buffer), :c) + cb_c = FerriteAssembly.get_coupled_buffer(buffer, :c) + @test 7 * ae ≈ FerriteAssembly.get_ae(cb_c) + end + end + function FerriteAssembly.element_routine!(Ke, re, state, ae, m::CS_MB, cv, buffer) + nothing # Only assembled from `:a`'s perspective in these tests + end + function FerriteAssembly.element_routine!(Ke, re, state, ae, m::CS_MB2, cv, buffer) + nothing + end + + a1 = rand(ndofs(dh1)) + a2 = zeros(ndofs(dh2)) + @assert length(a1) * 2 == length(a2) + a2[1:2:end] = 3 * a1 + a2[2:2:end] = 5 * a1 + aold1 = rand(ndofs(dh1)) + aold2 = zeros(ndofs(dh2)) + aold2[1:2:end] = aold1 + aold2[2:2:end] = aold1 + + @testset "threading=$threading, autodiffbuffer=$autodiffbuffer, singledomain=$singledomain" for + threading in (false, true), autodiffbuffer in (false, true), singledomain in (true, false) + if singledomain + d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1, threading, autodiffbuffer) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2, threading, autodiffbuffer) + else + sets = Dict(k => getcellset(grid, k) for k in ("left", "right")) + d1 = setup_domainbuffers(Dict(k => DomainSpec(dh1, CS_MA(), cvu; set) for (k, set) in sets); a = a1, threading, autodiffbuffer) + d2 = setup_domainbuffers(Dict(k => DomainSpec(dh2, CS_MB(), cvv; set) for (k, set) in sets); a = a2, threading, autodiffbuffer) + end + sim1 = Simulation(d1, a1, aold1) + sim2 = Simulation(d2, a2, aold2) + g = CoupledSimulations((a = sim1,); refs = (b = sim2,)) + @test g.a isa FerriteAssembly.CoupledSimulation + @test g.b === sim2 # refs are the plain source Simulation + + # A CoupledSimulation member forwards ordinary Simulation property access: `.a`/ + # `.aold` are the same global vectors (shared by reference, not copied), `.db` is the + # *rebuilt* (coupled) domain buffer, not the original source `d1`. + @test g.a.a === a1 + @test g.a.aold === aold1 + @test g.a.db === g.a.sim.db + @test g.a.db !== d1 + + # Property access must constant-propagate to a single concrete type (not a Union + # across every forwarding branch) for both the member handle's own getproperty + # override and the group's. + @test is_concrete_inferred(csim -> csim.a, typeof(g.a)) + @test is_concrete_inferred(csim -> csim.aold, typeof(g.a)) + @test is_concrete_inferred(csim -> csim.db, typeof(g.a)) + @test is_concrete_inferred(csim -> csim.sim, typeof(g.a)) + @test is_concrete_inferred(csim -> csim.partners, typeof(g.a)) + @test is_concrete_inferred(grp -> grp.a, typeof(g)) + + # Forwarded properties must tab-complete, not just the two real struct fields. + propnames = propertynames(g.a) + @test :a in propnames + @test :aold in propnames + @test :db in propnames + @test :sim in propnames + @test :partners in propnames + + # The group itself (`CoupledSimulations`) also forwards property access to its + # members/refs, must tab-complete accordingly, and must reject unknown names. + group_propnames = propertynames(g) + @test :primaries in group_propnames + @test :refs in group_propnames + @test :members in group_propnames + @test :a in group_propnames + @test :b in group_propnames + @test_throws ArgumentError g.nonexistent_member + + K = allocate_matrix(dh1) + r = zeros(ndofs(dh1)) + assembler = start_assemble(K, r) + expected_b_dt[] = NaN + work!(assembler, g.a) # Runs the @test's inside element_routine! + + # BUG-003 regression: changing a ref's Δt between two work! calls must be observed + # by every reader task, not just task 1, without re-working the ref itself. + set_time_increment!(g.b, 1.23) + expected_b_dt[] = 1.23 + work!(assembler, g.a) + set_time_increment!(g.b, 4.56) + expected_b_dt[] = 4.56 + work!(assembler, g.a) # element_routine! asserts Δt equality itself + expected_b_dt[] = NaN + + # Stable buffer/config identity across repeated work! calls (no rebuild per call) + get_ib() = singledomain ? FerriteAssembly.get_itembuffer(g.a) : FerriteAssembly.get_itembuffer(g.a, "left") + ib1 = FerriteAssembly.get_base(get_ib()) + work!(assembler, g.a) + ib2 = FerriteAssembly.get_base(get_ib()) + @test ib1 === ib2 + + # Coupling does not scale allocations with cell count (allow generous fixed overhead + # for task-spawn/chunk machinery; this is a smoke check, not a scaling sweep). + work!(assembler, g.a) # warm up (compile) + work!(assembler, g.a) # warm up again to be safe against any first-use effects + nalloc = @allocated work!(assembler, g.a) + @test nalloc < 2_000_000 + end + + @testset "coupling allocations do not scale with cell count" begin + # The single-mesh check above only bounds allocations against a fixed ceiling on one + # 4-cell grid; it cannot detect a small per-cell allocation (e.g. a reintroduced + # per-cell wrapper/config construction) that would still be far below that ceiling. + # Compare a much larger mesh against a tiny one instead: coupling-specific overhead + # (wrapper/config construction, task-spawn/chunk machinery) is paid once per `work!` + # call, not per cell, so it must not grow materially with cell count. + struct CS_AllocA end + struct CS_AllocB end + struct CS_AllocA0 end # uncoupled baseline: same per-cell work, no partner access + FerriteAssembly.create_cell_state(::CS_AllocA, args...) = nothing + FerriteAssembly.create_cell_state(::CS_AllocB, args...) = nothing + FerriteAssembly.create_cell_state(::CS_AllocA0, args...) = nothing + function FerriteAssembly.element_routine!(Ke, re, state, ae, ::CS_AllocA, cv, buffer) + cb = FerriteAssembly.get_coupled_buffer(buffer, :b) + ae_p = FerriteAssembly.get_ae(cb) + @inbounds for i in eachindex(re) + re[i] += ae_p[i] + end + return nothing + end + FerriteAssembly.element_routine!(Ke, re, state, ae, ::CS_AllocB, cv, buffer) = nothing + function FerriteAssembly.element_routine!(Ke, re, state, ae, ::CS_AllocA0, cv, buffer) + @inbounds for i in eachindex(re) + re[i] += ae[i] + end + return nothing + end + + function build_alloc_group(n) + grid_ = generate_grid(Quadrilateral, (n, n)) + ip_ = Lagrange{RefQuadrilateral,1}() + dhA = close!(add!(DofHandler(grid_), :u, ip_)) + dhB = close!(add!(DofHandler(grid_), :v, ip_^2)) + cvA = CellValues(qr, ip_, ip_) + cvB = CellValues(qr, ip_^2, ip_) + aA = zeros(ndofs(dhA)) + aB = zeros(ndofs(dhB)) + dA = setup_domainbuffer(DomainSpec(dhA, CS_AllocA(), cvA); a = aA) + dB = setup_domainbuffer(DomainSpec(dhB, CS_AllocB(), cvB); a = aB) + simA = Simulation(dA, aA, zeros(ndofs(dhA))) + simB = Simulation(dB, aB, zeros(ndofs(dhB))) + return CoupledSimulations((a = simA,); refs = (b = simB,)), dhA + end + function build_baseline(n) + grid_ = generate_grid(Quadrilateral, (n, n)) + ip_ = Lagrange{RefQuadrilateral,1}() + dhA = close!(add!(DofHandler(grid_), :u, ip_)) + cvA = CellValues(qr, ip_, ip_) + aA = zeros(ndofs(dhA)) + dA = setup_domainbuffer(DomainSpec(dhA, CS_AllocA0(), cvA); a = aA) + return Simulation(dA, aA, zeros(ndofs(dhA))), dhA + end + function measure_alloc(sim_or_group, dh_) + K = allocate_matrix(dh_) + r = zeros(ndofs(dh_)) + asm = start_assemble(K, r) + work!(asm, sim_or_group) # warm up (compile) + work!(asm, sim_or_group) # warm up again + return @allocated work!(asm, sim_or_group) + end + g_small, dh_small = build_alloc_group(2) + g_large, dh_large = build_alloc_group(20) # 100x the cells of g_small + nalloc_small = measure_alloc(g_small.a, dh_small) + nalloc_large = measure_alloc(g_large.a, dh_large) + + # Coupling-specific overhead relative to an uncoupled baseline doing equivalent + # per-cell work: both must be exactly zero (ordinary sequential assembly is + # allocation-free), so even a small per-cell allocation reintroduced by coupling + # would be caught, not just growth that outpaces cell count. + base_small = measure_alloc(build_baseline(2)...) + base_large = measure_alloc(build_baseline(20)...) + @test nalloc_small - base_small == 0 + @test nalloc_large - base_large == 0 + end + + @testset "threaded reader (1 task) with sequential partner" begin + # validate_domain_pair explicitly allows this (a sequential partner counts as 1 slot, + # matching a threaded reader with exactly 1 task); work! must not throw when scattering + # partners before dispatch, even though the partner has no task-local buffers to + # scatter into. + expected_b_dt[] = NaN + expected_b_material[] = CS_MB + d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1, threading = true, num_tasks = 1) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2, threading = false) + sim1 = Simulation(d1, a1, aold1) + sim2 = Simulation(d2, a2, aold2) + g = CoupledSimulations((a = sim1,); refs = (b = sim2,)) + K = allocate_matrix(dh1) + r = zeros(ndofs(dh1)) + work!(start_assemble(K, r), g.a) + end + + @testset "mutual coupling (3 members)" begin + struct CS_MC end + FerriteAssembly.create_cell_state(::CS_MC, cv, x, ae, args...) = nothing + function FerriteAssembly.element_routine!(Ke, re, state, ae, ::CS_MC, cv, buffer) + nothing + end + ip3 = Lagrange{RefQuadrilateral,1}() + dh3 = close!(add!(DofHandler(grid), :w, ip3)) + cv3 = CellValues(qr, ip3, ip3) + # Same dof ordering as dh1 (same grid/interpolation/single scalar field), so a + # component-wise multiple of a1 gives an independently checkable per-cell value, + # matching the existing a2/a1 pattern used for the :b partner above. + a3 = 7 * a1 + aold3 = zeros(ndofs(dh3)) + expected_b_material[] = CS_MB + for threading in (false, true) + d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1, threading) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2, threading) + d3 = setup_domainbuffer(DomainSpec(dh3, CS_MC(), cv3); a = a3, threading) + sim1 = Simulation(d1, a1, aold1) + sim2 = Simulation(d2, a2, aold2) + sim3 = Simulation(d3, a3, aold3) + g = CoupledSimulations((a = sim1, b = sim2, c = sim3)) + @test g.a isa FerriteAssembly.CoupledSimulation + @test g.b isa FerriteAssembly.CoupledSimulation + @test g.c isa FerriteAssembly.CoupledSimulation + cb_b = FerriteAssembly.get_coupled_buffers(FerriteAssembly.get_base(FerriteAssembly.get_itembuffer(g.a))) + @test haskey(cb_b, :b) && haskey(cb_b, :c) + # b and c views from a are nonrecursive: plain CellBuffer, no further coupling + @test !hasmethod(FerriteAssembly.get_coupled_buffers, Tuple{typeof(cb_b.b)}) + + # Actually work! the reader with two distinct partners: CS_MA's element_routine! + # checks both :b (2x/2y-scaled dof/state values) and :c (7x-scaled dof values), + # so a positional mix-up between the two partner NamedTuples (buffers vs. + # simulations) would fail here even though it could pass a construction-only check. + expected_b_dt[] = NaN + K = allocate_matrix(dh1) + r = zeros(ndofs(dh1)) + work!(start_assemble(K, r), g.a) + end + end + + @testset "autodiff through coupling: numerical agreement" begin + # `CS_AD_Reader` defines only `element_residual!` (no `element_routine!`), so its Ke + # is genuinely computed via ForwardDiff through the coupled partner buffer, not a + # hand-written Ke. `CS_AD_Reader_manual` computes the analytically-known Ke (c*I, + # since the residual is elementwise linear) directly via `element_routine!`, coupled + # to the same partner. Assembling both through the same mesh/assembly machinery and + # comparing the results verifies the AD-through-coupling path numerically, for both + # sequential and threaded execution. + struct CS_AD_Reader + c::Float64 + end + struct CS_AD_Reader_manual + c::Float64 + end + struct CS_AD_Partner end + FerriteAssembly.create_cell_state(::CS_AD_Reader, args...) = nothing + FerriteAssembly.create_cell_state(::CS_AD_Reader_manual, args...) = nothing + FerriteAssembly.create_cell_state(::CS_AD_Partner, args...) = nothing + + function FerriteAssembly.element_residual!(re, state, ae, m::CS_AD_Reader, cv, buffer) + ae_p = FerriteAssembly.get_ae(FerriteAssembly.get_coupled_buffer(buffer, :p)) + re .= m.c .* ae .- ae_p + return nothing + end + function FerriteAssembly.element_routine!(Ke, re, state, ae, m::CS_AD_Reader_manual, cv, buffer) + ae_p = FerriteAssembly.get_ae(FerriteAssembly.get_coupled_buffer(buffer, :p)) + re .= m.c .* ae .- ae_p + fill!(Ke, 0) + for i in axes(Ke, 1) + Ke[i, i] = m.c + end + return nothing + end + function FerriteAssembly.element_routine!(Ke, re, state, ae, ::CS_AD_Partner, cv, buffer) + nothing + end + + ipr = Lagrange{RefQuadrilateral,1}() + dhr = close!(add!(DofHandler(grid), :r, ipr)) + cvr = CellValues(qr, ipr, ipr) + ar = rand(ndofs(dhr)) + ap = rand(ndofs(dhr)) + aold_dummy = zeros(ndofs(dhr)) + c = 2.5 + + for threading in (false, true) + dr_ad = setup_domainbuffer(DomainSpec(dhr, CS_AD_Reader(c), cvr); a = ar, threading, autodiffbuffer=true) + dr_man = setup_domainbuffer(DomainSpec(dhr, CS_AD_Reader_manual(c), cvr); a = ar, threading) + dp = setup_domainbuffer(DomainSpec(dhr, CS_AD_Partner(), cvr); a = ap, threading) + simr_ad = Simulation(dr_ad, ar, aold_dummy) + simr_man = Simulation(dr_man, ar, aold_dummy) + simp = Simulation(dp, ap, aold_dummy) + g_ad = CoupledSimulations((r = simr_ad,); refs = (p = simp,)) + g_man = CoupledSimulations((r = simr_man,); refs = (p = simp,)) + K_ad = allocate_matrix(dhr); r_ad = zeros(ndofs(dhr)) + K_man = allocate_matrix(dhr); r_man = zeros(ndofs(dhr)) + work!(start_assemble(K_ad, r_ad), g_ad.r) + work!(start_assemble(K_man, r_man), g_man.r) + @test Matrix(K_ad) ≈ Matrix(K_man) + @test r_ad ≈ r_man + end + end + + @testset "replace_material through group" begin + d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2) + sim1 = Simulation(d1, a1, aold1) + sim2 = Simulation(d2, a2, aold2) + g = CoupledSimulations((a = sim1,); refs = (b = sim2,)) + K = allocate_matrix(dh1); r = zeros(ndofs(dh1)) + expected_b_dt[] = NaN + expected_b_material[] = CS_MB + work!(start_assemble(K, r), g.a) # exercise before replacement, observes CS_MB + + g2 = FerriteAssembly.replace_material(g, :b, m -> CS_MB2()) # changes the material TYPE + @test g2.a isa FerriteAssembly.CoupledSimulation + @test FerriteAssembly.get_material(g2.b) isa CS_MB2 + @test FerriteAssembly.get_material(g.b) isa CS_MB # old group/handle untouched + @test g2.a !== g.a + + # New handle actually works and observes the new material through coupling + expected_b_material[] = CS_MB2 + work!(start_assemble(K, r), g2.a) + # Old handle, worked again, still observes the original material (not silently switched) + expected_b_material[] = CS_MB + work!(start_assemble(K, r), g.a) + expected_b_material[] = CS_MB + + @test_throws ArgumentError FerriteAssembly.replace_material(g, :nope, identity) + @test_throws ArgumentError FerriteAssembly.replace_material(g.a, identity) + end + + @testset "replace_material with domain selector" begin + # Both members are multi-domain (matching "left"/"right" keys), as required by + # coupling's mixed single-domain/dictionary validation. + sets = Dict(k => getcellset(grid, k) for k in ("left", "right")) + d1m = setup_domainbuffers(Dict(k => DomainSpec(dh1, CS_MA(), cvu; set) for (k, set) in sets); a = a1) + d2m = setup_domainbuffers(Dict(k => DomainSpec(dh2, CS_MB(), cvv; set) for (k, set) in sets); a = a2) + sim1m = Simulation(d1m, a1, aold1) + sim2m = Simulation(d2m, a2, aold2) + g = CoupledSimulations((a = sim1m,); refs = (b = sim2m,)) + + # Success path: only the named domain's material is swapped, the other domain and + # the old group/handle are untouched. + g2 = FerriteAssembly.replace_material(g, :a, m -> CS_MB2(); domain = "left") + @test FerriteAssembly.get_material(g2.a.db["left"]) isa CS_MB2 + @test FerriteAssembly.get_material(g2.a.db["right"]) isa CS_MA + @test FerriteAssembly.get_material(g.a.db["left"]) isa CS_MA + + # Error path: `domain` given for a member whose `db` is not a domain dictionary + # (uses the plain single-domain group from the "replace_material through group" case). + d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2) + gs = CoupledSimulations((a = Simulation(d1, a1, aold1),); refs = (b = Simulation(d2, a2, aold2),)) + @test_throws ArgumentError FerriteAssembly.replace_material(gs, :a, identity; domain = "left") + end + + @testset "CoupledCellBuffer forwarding (dof_range, direct replace_material)" begin + d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2) + sim1 = Simulation(d1, a1, aold1) + sim2 = Simulation(d2, a2, aold2) + g = CoupledSimulations((a = sim1,); refs = (b = sim2,)) + cb = FerriteAssembly.get_base(FerriteAssembly.get_itembuffer(g.a)) + @test cb isa FerriteAssembly.CoupledCellBuffer + @test Ferrite.dof_range(cb, :u) == Ferrite.dof_range(cb.primary, :u) + + # Internal-plumbing check, not a demonstration of a supported user workflow: calling + # `replace_material` directly on a `.db` that already wraps `CoupledCellBuffer`s (as + # opposed to the documented, group-level `replace_material(group, member, f)`) must + # still dispatch correctly and preserve the partner buffers by reference. The result + # is not itself re-workable through `work!`, since the partner *Simulations* needed + # for reinitialization live on the owning `CoupledSimulation`, not on `.db` alone. + new_db = FerriteAssembly.replace_material(g.a.db, m -> CS_MB2()) + @test FerriteAssembly.get_material(new_db) isa CS_MB2 + new_cb = FerriteAssembly.get_base(FerriteAssembly.get_itembuffer(new_db)) + @test new_cb isa FerriteAssembly.CoupledCellBuffer + @test FerriteAssembly.get_coupled_buffers(new_cb) === FerriteAssembly.get_coupled_buffers(cb) + end + + @testset "empty multi-domain reader" begin + # Construction-only edge case: a primary member whose own `db` is an empty domain + # dictionary must succeed (nothing to couple), rather than erroring while trying to + # iterate it. This does not imply such a member is workable via `work!` afterwards + # (its `Dict` value type is the abstract `AbstractDomainBuffer`, which does not match + # the concrete-eltype bound `work!`'s multi-domain dispatch requires) - only that + # `CoupledSimulations` construction itself tolerates it. + empty_db = Dict{String, FerriteAssembly.AbstractDomainBuffer}() + sim_empty = Simulation(empty_db, Float64[], Float64[]) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2) + sim2 = Simulation(d2, a2, aold2) + g = CoupledSimulations((a = sim_empty,); refs = (b = sim2,)) + @test g.a isa FerriteAssembly.CoupledSimulation + @test g.a.db isa Dict + @test isempty(g.a.db) + end + + @testset "validation errors" begin + d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2) + sim1 = Simulation(d1, a1, aold1) + sim2 = Simulation(d2, a2, aold2) + + @test_throws ArgumentError CoupledSimulations(NamedTuple()) # empty primaries + @test_throws ArgumentError CoupledSimulations((a = sim1,); refs = (a = sim2,)) # duplicate name + + # Different grid + grid2 = generate_grid(Quadrilateral, (2,2)) + ip2 = Lagrange{RefQuadrilateral,1}() + dh2b = close!(add!(DofHandler(grid2), :v, ip2^2)) + a2b = zeros(ndofs(dh2b)) + d2b = setup_domainbuffer(DomainSpec(dh2b, CS_MB(), CellValues(qr, ip2^2, ip2)); a = a2b) + sim2b = Simulation(d2b, a2b, zeros(ndofs(dh2b))) + @test_throws ArgumentError CoupledSimulations((a = sim1,); refs = (b = sim2b,)) + + # Missing domain coverage (multi-domain reader, partner missing a domain) + sets = Dict(k => getcellset(grid, k) for k in ("left", "right")) + d1m = setup_domainbuffers(Dict(k => DomainSpec(dh1, CS_MA(), cvu; set) for (k, set) in sets); a = a1) + d2m_partial = Dict("left" => setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv; set=sets["left"]); a = a2)) + sim1m = Simulation(d1m, a1, aold1) + sim2m_partial = Simulation(d2m_partial, a2, aold2) + @test_throws ArgumentError CoupledSimulations((a = sim1m,); refs = (b = sim2m_partial,)) + + # Mixed single/dictionary coupling + @test_throws ArgumentError CoupledSimulations((a = sim1m,); refs = (b = sim2,)) + @test_throws ArgumentError CoupledSimulations((a = sim1,); refs = (b = Simulation(d2m_partial, a2, aold2),)) + + # Incompatible task counts + d1t = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1, threading=true, num_tasks=2) + d2t = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2, threading=true, num_tasks=3) + sim1t = Simulation(d1t, a1, aold1) + sim2t = Simulation(d2t, a2, aold2) + @test_throws ArgumentError CoupledSimulations((a = sim1t,); refs = (b = sim2t,)) + + # Duplicate storage: same domain buffer object used under two member names + @test_throws ArgumentError CoupledSimulations((a = sim1,); refs = (b = sim2, c = sim2)) + + # Unsupported buffer kind (facet buffers): rejected with an actionable ArgumentError, + # not a MethodError, both as a partner and as a sole primary with no partners at all + # (storage-identity validation runs for every member, regardless of pairing). + struct CS_MFacet end + dh_f = close!(add!(DofHandler(grid), :f, ip)) + fv = FacetValues(FacetQuadratureRule{RefQuadrilateral}(2), ip) + d_facet = setup_domainbuffer(DomainSpec(dh_f, CS_MFacet(), fv; set=getfacetset(grid, "left"))) + a_f = zeros(ndofs(dh_f)) + sim_facet = Simulation(d_facet, a_f, zeros(ndofs(dh_f))) + @test_throws ArgumentError CoupledSimulations((a = sim_facet,)) + @test_throws ArgumentError CoupledSimulations((a = sim1,); refs = (b = sim_facet,)) + end +end diff --git a/test/replacements.jl b/test/replacements.jl index 65d7e061..74979f75 100644 --- a/test/replacements.jl +++ b/test/replacements.jl @@ -1,108 +1,47 @@ -@testset "replace_material" begin - m_el = EE.LinearElastic(;E=1.0, ν=0.4) - m_elx2 = EE.LinearElastic(;E=2.0, ν=0.4) - f_repl1(::EE.LinearElastic) = m_elx2 - m_pl = EE.J2Plasticity(;E=1.0, ν=0.4, σ0=0.2, H=1.0) - f_repl2(::EE.LinearElastic) = m_pl - f_repl3(::EE.J2Plasticity) = m_el - grid = generate_grid(Quadrilateral, (2,2)) - ip = Lagrange{RefQuadrilateral,1}()^2 - dh = DofHandler(grid); add!(dh, :u, ip); close!(dh) - qr = QuadratureRule{RefQuadrilateral}(2) - cv = CellValues(qr, ip, ip) - dspec = DomainSpec(dh, m_el, cv) - buffer = setup_domainbuffer(dspec) - ad_buffer = setup_domainbuffer(dspec; autodiffbuffer=true) - td_buffer = setup_domainbuffer(dspec; threading=true) - - for b0 in (buffer, ad_buffer, td_buffer) - @test FerriteAssembly.get_material(b0) === m_el - b1 = FerriteAssembly.replace_material(b0, f_repl1) - @test FerriteAssembly.get_material(b1) === m_elx2 - b2 = FerriteAssembly.replace_material(b1, f_repl2) - @test FerriteAssembly.get_material(b2) === m_pl - b3 = FerriteAssembly.replace_material(b2, f_repl3) - @test FerriteAssembly.get_material(b3) === m_el - end - - n_half = getncells(grid)÷2 - buffers = setup_domainbuffers(Dict( - "a" => DomainSpec(dh, m_el, cv; set=1:(n_half-1)), - "b" => DomainSpec(dh, m_pl, cv; set=n_half:getncells(grid)) - )) - @test FerriteAssembly.get_material(buffers, "a") === m_el - @test FerriteAssembly.get_material(buffers, "b") === m_pl - f_repl(::EE.LinearElastic) = m_elx2 - f_repl(m::EE.J2Plasticity) = m - bs2 = FerriteAssembly.replace_material(buffers, f_repl) - @test FerriteAssembly.get_material(bs2, "a") === m_elx2 - @test FerriteAssembly.get_material(bs2, "b") === m_pl -end - -@testset "couple_buffers" begin - grid = generate_grid(Quadrilateral, (2,2)) - addcellset!(grid, "left", x -> x[1] < eps()) - addcellset!(grid, "right", setdiff(1:getncells(grid), getcellset(grid, "left"))) - ip = Lagrange{RefQuadrilateral,1}() - dh1 = close!(add!(DofHandler(grid), :u, ip)) - dh2 = close!(add!(DofHandler(grid), :v, ip^2)) - qr = QuadratureRule{RefQuadrilateral}(2) - cvu = CellValues(qr, ip, ip) - cvv = CellValues(qr, ip^2, ip) - - struct MA end - struct MB end - # We will test with the following dof value differences - # aold will be same in both cases (for both components in the case of MB) - # a will be 3 times larger for first component in MB, and 5 times for second component - # State will be 6 times larger for MB, obtained by multiplying the function values by factor 2 - FerriteAssembly.create_cell_state(::MA, cv, x, ae, args...) = [function_value(cv, i, ae) for i in 1:getnquadpoints(cv)] - FerriteAssembly.create_cell_state(::MB, cv, x, ae, args...) = [2 * function_value(cv, i, ae)[1] for i in 1:getnquadpoints(cv)] - - # Test case to check that values have been updated correctly - function FerriteAssembly.element_routine!(Ke, re, state, ae, m::MA, cv, buffer) - cb_b = FerriteAssembly.get_coupled_buffer(buffer, :b) - # Check that correct material has been set - @test FerriteAssembly.get_material(cb_b) isa MB - # Check that dofs have been updated - @test 3 * ae ≈ FerriteAssembly.get_ae(cb_b)[1:2:end] # 1st component - @test 5 * ae ≈ FerriteAssembly.get_ae(cb_b)[2:2:end] # 2nd component - # Check that old dofs have been updated - @test FerriteAssembly.get_aeold(buffer) ≈ FerriteAssembly.get_aeold(cb_b)[1:2:end] - @test FerriteAssembly.get_aeold(buffer) ≈ FerriteAssembly.get_aeold(cb_b)[2:2:end] - # Check that state variables have been updated - @test 6 * state ≈ FerriteAssembly.get_state(cb_b) - end - - a1 = rand(ndofs(dh1)) - a2 = zeros(ndofs(dh2)) - @assert length(a1) * 2 == length(a2) - a2[1:2:end] = 3 * a1 - a2[2:2:end] = 5 * a1 - aold1 = rand(ndofs(dh1)) - aold2 = zeros(ndofs(dh2)) - aold2[1:2:end] = aold1; - aold2[2:2:end] = aold1; - - for threading in (false, true) - for autodiffbuffer in (false, true) - for singledomain in (true, false) - if singledomain - d1 = setup_domainbuffer(DomainSpec(dh1, MA(), cvu); a = a1, threading, autodiffbuffer) - d2 = setup_domainbuffer(DomainSpec(dh2, MB(), cvv); a = a2, threading, autodiffbuffer) - else - sets = Dict(k => getcellset(grid, k) for k in ("left", "right")) - d1 = setup_domainbuffers(Dict(k => DomainSpec(dh1, MA(), cvu; set) for (k, set) in sets); a = a1, threading, autodiffbuffer) - d2 = setup_domainbuffers(Dict(k => DomainSpec(dh2, MB(), cvv; set) for (k, set) in sets); a = a2, threading, autodiffbuffer) - end - d1 = couple_buffers(d1; b = d2) - sim1 = Simulation(d1, a1, aold1) - sim2 = Simulation(d2, a2, aold2) - K = allocate_matrix(dh1) - r = zeros(ndofs(dh1)) - assembler = start_assemble(K, r) - work!(assembler, sim1, CoupledSimulations(b = sim2)) # Test - end - end - end -end +@testset "replace_material" begin + m_el = EE.LinearElastic(;E=1.0, ν=0.4) + m_elx2 = EE.LinearElastic(;E=2.0, ν=0.4) + f_repl1(::EE.LinearElastic) = m_elx2 + m_pl = EE.J2Plasticity(;E=1.0, ν=0.4, σ0=0.2, H=1.0) + f_repl2(::EE.LinearElastic) = m_pl + f_repl3(::EE.J2Plasticity) = m_el + grid = generate_grid(Quadrilateral, (2,2)) + ip = Lagrange{RefQuadrilateral,1}()^2 + dh = DofHandler(grid); add!(dh, :u, ip); close!(dh) + qr = QuadratureRule{RefQuadrilateral}(2) + cv = CellValues(qr, ip, ip) + dspec = DomainSpec(dh, m_el, cv) + buffer = setup_domainbuffer(dspec) + ad_buffer = setup_domainbuffer(dspec; autodiffbuffer=true) + td_buffer = setup_domainbuffer(dspec; threading=true) + + for b0 in (buffer, ad_buffer, td_buffer) + @test FerriteAssembly.get_material(b0) === m_el + b1 = FerriteAssembly.replace_material(b0, f_repl1) + @test FerriteAssembly.get_material(b1) === m_elx2 + b2 = FerriteAssembly.replace_material(b1, f_repl2) + @test FerriteAssembly.get_material(b2) === m_pl + b3 = FerriteAssembly.replace_material(b2, f_repl3) + @test FerriteAssembly.get_material(b3) === m_el + end + + n_half = getncells(grid)÷2 + buffers = setup_domainbuffers(Dict( + "a" => DomainSpec(dh, m_el, cv; set=1:(n_half-1)), + "b" => DomainSpec(dh, m_pl, cv; set=n_half:getncells(grid)) + )) + @test FerriteAssembly.get_material(buffers, "a") === m_el + @test FerriteAssembly.get_material(buffers, "b") === m_pl + f_repl(::EE.LinearElastic) = m_elx2 + f_repl(m::EE.J2Plasticity) = m + bs2 = FerriteAssembly.replace_material(buffers, f_repl) + @test FerriteAssembly.get_material(bs2, "a") === m_elx2 + @test FerriteAssembly.get_material(bs2, "b") === m_pl + + # Domain-selective replacement + bs3 = FerriteAssembly.replace_material(buffers, "a", f_repl1) + @test FerriteAssembly.get_material(bs3, "a") === m_elx2 + @test FerriteAssembly.get_material(bs3, "b") === m_pl # unchanged, f_repl1 not applied here + @test bs3["b"] === buffers["b"] # copied by reference + @test_throws ArgumentError FerriteAssembly.replace_material(buffers, "c", f_repl1) +end diff --git a/test/runtests.jl b/test/runtests.jl index 32035017..97c48f3e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,7 +9,8 @@ import MechanicalMaterialModels as MMM using Logging include("replacements.jl") -include("states.jl") +include("coupled_simulations.jl") +include("states.jl") include("threading_utils.jl") include("heatequation.jl") include("example_elements.jl") diff --git a/test/setup.jl b/test/setup.jl index 40b757f5..762f422d 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -21,7 +21,7 @@ aold_value = rand() aold = ones(ndofs(dh))*aold_value _getdomain(dbs::Dict, key::String) = dbs[key] - _getdomain(sim::Simulation, key) = FerriteAssembly.get_domain_simulation(sim, key) + _getdomain(sim::Simulation, key) = Simulation(sim.db[key], sim.a, sim.aold) for container in (buffers, buffers_ad, Simulation(buffers, nothing, aold), Simulation(buffers_ad, nothing, aold)) # Basic access functions @test FerriteAssembly.get_dofhandler(container) === dh @@ -36,7 +36,7 @@ cell_id = first(cellset) cb1 = FerriteAssembly.get_itembuffer(cont1) sim = isa(container, Simulation) ? cont1 : Simulation(cont1, nothing, aold) - FerriteAssembly.reinit_buffer!(cb1, sim, CoupledSimulations(), cell_id) + FerriteAssembly.reinit_buffer!(cb1, sim, cell_id) @test FerriteAssembly.get_user_data(cb1) === userdata @test FerriteAssembly.get_user_cache(cb1) == [1.0] ae_old = FerriteAssembly.get_aeold(cb1) @@ -82,7 +82,7 @@ aold = ones(ndofs(dh))*aold_value facetbuffer = FerriteAssembly.get_itembuffer(buffer) facet_id = first(FerriteAssembly.getset(buffer)) - FerriteAssembly.reinit_buffer!(facetbuffer, Simulation(buffer, zeros(ndofs(dh)), aold), CoupledSimulations(), facet_id) + FerriteAssembly.reinit_buffer!(facetbuffer, Simulation(buffer, zeros(ndofs(dh)), aold), facet_id) @test FerriteAssembly.get_user_data(facetbuffer) === userdata @test FerriteAssembly.get_user_cache(facetbuffer) == [1.0] @test FerriteAssembly.get_user_cache(facetbuffer) !== FerriteAssembly.get_user_cache(FerriteAssembly.get_itembuffer(buffers["right"]))