From d9f6ca3aa2160c91e4a12345eb9f2299c2db50be Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Mon, 14 Sep 2026 12:18:52 -0400 Subject: [PATCH 01/16] Replace implicit coupled-buffer API with explicit setup-time CoupledSimulations MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Remove couple_buffers, the keyword-only CoupledSimulations partner container, and the 3-arg coupled work! interface. Replace with a setup-time CoupledSimulations((primaries...); refs=(...)) group constructor that validates and wires coupling once, producing CoupledMember handles whose work! needs no extra arguments. Fixes BUG-003 (stale partner Δt) at the root: coupled cell buffers now hold direct, nonrecursive references into the partner's own mutable CellBuffer and Simulation (never deepcopies), and work!(worker, ::CoupledMember, ...) scatters each declared partner's task-local buffers from its base before dispatch, so a threaded reader always observes the partner's current state (time increment included) without re-working the partner. Key pieces: - New CoupledCellBuffer (src/ItemBuffers/CoupledCellBuffer.jl) wraps a reader's own CellBuffer plus per-partner (CellBuffer, Simulation) pairs; reinit_buffer! reinitializes the primary and every direct partner, with no recursion into partners' own coupling. - New src/Coupling.jl: CoupledSimulations/CoupledMember construction, setup-time validation (matching grid, single/dictionary domain shape, supported buffer kinds, cell coverage, per-domain task-count compatibility, distinct mutable scratch across all members including replace_material'd buffers), task-count-aware partner binding for threaded readers, and group-level replace_material (whole-member or single named domain). - CellBuffer/AutoDiffCellBuffer/FacetBuffer/work.jl: removed coupling fields, couple_buffers methods, and the coupled argument from reinit_buffer!/work!; AutoDiffCellBuffer and its residual/Jacobian construction now also accept CoupledCellBuffer, built fresh after coupling is wired. - Migrated the phase-field fracture tutorial and coupling tests to the new API; added test/coupled_simulations.jl covering mutual/one-way/mixed groups, the BUG-003 Δt regression, autodiff-through-coupling numerical agreement against a hand-differentiated reference, replace_material through a group (including type changes), and setup-validation failures. Validation: `Pkg.test()` passes in full (CoupledSimulations testset: 1515 assertions, 0 failures; whole suite green). All literate tutorials/howtos run clean, including phasefield_fracture.jl (the only coupled-buffer user); full `docs/make.jl` HTML build completes without error. Reviewed with Codex (dual-review skill) at the plan stage (8 findings, all incorporated before implementation: partner scatter-before-dispatch, strict threaded task-count matching, per-domain execution records, facet/autodiff signature updates, partner-buffer unwrapping, cross-member storage-identity checks, empty-domain handling) and against the final diff (6 findings: scratch-array-identity storage check surviving replace_material, empty multi-domain Dict typing, positive-task-count validation independent of pairing, reserved member names, documented concurrency contract, strengthened autodiff/replace_material tests — all fixed and re-verified). Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- docs/src/DomainBuffers/Setup.md | 24 +- .../literate_tutorials/phasefield_fracture.jl | 17 +- src/Autodiff/autodiff.jl | 20 +- src/Coupling.jl | 301 ++++++++++++++++++ src/DomainBuffers.jl | 36 +-- src/FerriteAssembly.jl | 4 +- src/ItemBuffers/AbstractItemBuffer.jl | 8 +- src/ItemBuffers/CellBuffer.jl | 39 +-- src/ItemBuffers/CoupledCellBuffer.jl | 60 ++++ src/ItemBuffers/FacetBuffer.jl | 2 +- src/Simulation.jl | 24 -- src/work.jl | 211 ++++++------ test/coupled_simulations.jl | 261 +++++++++++++++ test/replacements.jl | 155 +++------ test/runtests.jl | 3 +- test/setup.jl | 4 +- 16 files changed, 840 insertions(+), 329 deletions(-) create mode 100644 src/Coupling.jl create mode 100644 src/ItemBuffers/CoupledCellBuffer.jl create mode 100644 test/coupled_simulations.jl 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..2b0a146b 100644 --- a/docs/src/literate_tutorials/phasefield_fracture.jl +++ b/docs/src/literate_tutorials/phasefield_fracture.jl @@ -158,8 +158,10 @@ 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)); +sim_u0 = Simulation(db_u_uc, zeros(ndofs_u), zeros(ndofs_u)) +sim_d0 = Simulation(db_d_uc, zeros(ndofs_d), zeros(ndofs_d)) +g = CoupledSimulations((u = sim_u0, d = sim_d0)) # `:u` and `:d` mutually read each other +sim_u, sim_d = g.u, g.d; # Setup loading and boundary conditions load_function(t) = 1e-4 * t @@ -177,8 +179,9 @@ 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) +# Write function to solve one simulation part. The coupling to the other part is already +# wired into `sim` (a `CoupledSimulations` member handle), so `work!` needs no extra input. +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..11e35576 --- /dev/null +++ b/src/Coupling.jl @@ -0,0 +1,301 @@ +# Explicit, setup-time coupling between simulations. +# +# A `CoupledSimulations` group is built once from a set of named `Simulation`s. Each +# primary member's domain buffer(s) are rebuilt to hold a `CoupledCellBuffer` (or +# `AutoDiffCellBuffer{<:CoupledCellBuffer}`) itembuffer that references the *actual* +# mutable buffers/simulations of its declared partners (primaries and refs alike, except +# itself). No coupling is resolved or discovered during `work!`. + +""" + CoupledMember(sim, partner_containers) + +A handle to one primary member of a [`CoupledSimulations`](@ref) group. `sim` is a +[`Simulation`](@ref) whose domain buffer(s) have been rebuilt with coupled itembuffers. +`partner_containers` are the raw itembuffer containers (`TaskLocals` or plain) of every +partner this member reads from, scattered once at the start of each `work!` call so that +threaded partners' task-local buffers observe the partner's current base state (e.g. its +time increment) even if the partner itself has not been `work!`ed since it last changed. + +Forwards the ordinary [`Simulation`](@ref) accessor API (`.a`, `.aold`, `.db`, +`get_dofhandler`, `get_state`, `set_time_increment!`, `update_states!`, etc.). +""" +struct CoupledMember{S<:Simulation, PT<:Tuple} + sim::S + partner_containers::PT +end + +function Base.getproperty(m::CoupledMember, name::Symbol) + name === :sim && return getfield(m, :sim) + name === :partner_containers && return getfield(m, :partner_containers) + return getproperty(getfield(m, :sim), name) +end + +get_material(m::CoupledMember, args::Vararg{Any,N}) where N = get_material(getfield(m, :sim), args...) +get_dofhandler(m::CoupledMember) = get_dofhandler(getfield(m, :sim)) +get_grid(m::CoupledMember) = get_grid(getfield(m, :sim)) +get_state(m::CoupledMember, args::Vararg{Any,N}) where N = get_state(getfield(m, :sim), args...) +get_old_state(m::CoupledMember, args::Vararg{Any,N}) where N = get_old_state(getfield(m, :sim), args...) +getset(m::CoupledMember, args::Vararg{Any,N}) where N = getset(getfield(m, :sim), args...) +update_states!(m::CoupledMember; kwargs...) = update_states!(getfield(m, :sim); kwargs...) +set_time_increment!(m::CoupledMember, Δt) = set_time_increment!(getfield(m, :sim), Δt) +revert_states!(m::CoupledMember) = revert_states!(getfield(m, :sim)) +get_itembuffer(m::CoupledMember, args::Vararg{Any,N}) where N = get_itembuffer(getfield(m, :sim), args...) +get_num_tasks(m::CoupledMember) = get_num_tasks(getfield(m, :sim)) +get_chunks(m::CoupledMember) = get_chunks(getfield(m, :sim)) + +replace_material(::CoupledMember, 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.")) + +_scatter_partner!(c::TaskLocals) = scatter!(c) +_scatter_partner!(::Any) = nothing + +function work!(worker, m::CoupledMember, args...; kwargs...) + foreach(_scatter_partner!, getfield(m, :partner_containers)) + return work!(worker, getfield(m, :sim), args...; kwargs...) +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_member(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 + +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, partner_sims::NamedTuple) + autodiff = _is_autodiff(reader_ibuf) + wrap(primary_cb, partners_nt) = autodiff ? + AutoDiffCellBuffer(CoupledCellBuffer(primary_cb, partners_nt, partner_sims)) : + CoupledCellBuffer(primary_cb, partners_nt, partner_sims) + 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 + rn = get_num_tasks(reader_db) + rn > 0 || throw(ArgumentError("task count must be positive")) + pn = partner_db isa ThreadedDomainBuffer ? get_num_tasks(partner_db) : 1 + rn == pn || throw(ArgumentError( + "threaded reader with $rn tasks requires coupling partner `$partner_name` to provide $rn " * + "task-local buffers (a sequential partner counts as 1 slot); got $pn")) + end + return nothing +end + +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) + partner_sims = map(p -> p.sim, partners) + coupled_ibuf = build_coupled_itembuffer(reader_ibuf, partner_containers, partner_sims) + new_db = setproperties(reader_db; itembuffer = coupled_ibuf) + return new_db, values(partner_containers) +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_member(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 + containers = Any[] + else + containers = Any[] + built = 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, conts = build_coupled_domain(rdb, partners) + push!(built, dname => ndb) + append!(containers, conts) + end + new_db = Dict(built...) # infers the narrowest common concrete value type, matching MultiDomain(Threaded)Sim dispatch + end + 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, conts = build_coupled_domain(reader_db, partners) + containers = collect(conts) + end + new_sim = Simulation(new_db, reader_sim.a, reader_sim.aold) + return CoupledMember(new_sim, Tuple(unique(containers))) +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..cba8f6f4 100644 --- a/src/FerriteAssembly.jl +++ b/src/FerriteAssembly.jl @@ -17,10 +17,12 @@ include("Simulation.jl") include("setup.jl") include("ItemBuffers/CellBuffer.jl") +include("ItemBuffers/CoupledCellBuffer.jl") include("ItemBuffers/FacetBuffer.jl") include("Autodiff/autodiff.jl") include("work.jl") +include("Coupling.jl") include("Workers/Assemblers.jl") include("Workers/Integrators.jl") include("Workers/QuadPointEvaluator.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..68fc9b86 --- /dev/null +++ b/src/ItemBuffers/CoupledCellBuffer.jl @@ -0,0 +1,60 @@ +""" + CoupledCellBuffer(primary::CellBuffer, partner_buffers::NamedTuple, partner_sims::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)) +and the partner `Simulation`s required to reinitialize them (`partner_sims`, internal only). + +`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. +""" +struct CoupledCellBuffer{CB<:CellBuffer, PB<:NamedTuple, PS<:NamedTuple} <: AbstractCellBuffer + primary::CB + partner_buffers::PB + partner_sims::PS +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) + +""" + reinit_buffer!(cb::CoupledCellBuffer, sim::Simulation, cellnum::Int) + +Reinitialize the reader's own `cb.primary` against `sim`, then reinitialize each +partner buffer in `cb.partner_buffers` against its own partner `Simulation` stored in +`cb.partner_sims`. Partner reinitialization does not recurse: partner buffers are plain +`CellBuffer`s, so no further coupling initialization happens. +""" +function reinit_buffer!(cb::CoupledCellBuffer, sim::Simulation, cellnum::Int) + reinit_buffer!(cb.primary, sim, cellnum) + reinit_partners!(cb.partner_buffers, cb.partner_sims, cellnum) + return nothing +end + +function reinit_partners!(buffers::NamedTuple, sims::NamedTuple, cellnum::Int) + map((b, s) -> (reinit_buffer!(b, s, cellnum); nothing), buffers, sims) + return nothing +end + +function _replace_material_with(cb::CoupledCellBuffer, new_material) + new_primary = _replace_material_with(cb.primary, new_material) + return CoupledCellBuffer(new_primary, cb.partner_buffers, cb.partner_sims) +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..5acee157 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -50,27 +50,3 @@ end Base.iterate(sim::Simulation{<:DomainBuffers}) = _iterate(sim, iterate(sim.db)) Base.iterate(sim::Simulation{<:DomainBuffers}, iter) = _iterate(sim, iterate(sim.db, iter)) -""" - CoupledSimulations(; key1 = sim1::Simulation, key2 = sim2::Simulation, ...) - -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). -""" -struct CoupledSimulations{NT <: NamedTuple{<:Any, <:NTuple{<:Any, Simulation}}} - sims::NT -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 diff --git a/src/work.jl b/src/work.jl index cd8a6d0d..bc73aa41 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::MultiDomainSim) + for (name, sim) in multisim + skip_this_domain(worker, name) && continue + work_domain_sequential!(worker, sim) + end +end +function work!(worker, sim::SingleDomainSim) + work_domain_sequential!(worker, sim) +end +function work!(worker, multisim::MultiDomainThreadedSim) + 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::SingleDomainThreadedSim) + 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::Simulation{<: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::SingleDomainThreadedSim) + 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, 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/coupled_simulations.jl b/test/coupled_simulations.jl new file mode 100644 index 00000000..508cbcbb --- /dev/null +++ b/test/coupled_simulations.jl @@ -0,0 +1,261 @@ +@testset "CoupledSimulations" 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 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[] + 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.CoupledMember + @test g.b === sim2 # refs are the plain source Simulation + + 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 "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) + a3 = zeros(ndofs(dh3)) + aold3 = zeros(ndofs(dh3)) + d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1) + d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2) + d3 = setup_domainbuffer(DomainSpec(dh3, CS_MC(), cv3); a = a3) + 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.CoupledMember + @test g.b isa FerriteAssembly.CoupledMember + @test g.c isa FerriteAssembly.CoupledMember + 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)}) + 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.CoupledMember + @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 "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)) + 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..53a4e596 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -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"])) From 0ed0a06638637a211bc97cf4ec2c1249ec7b29ba Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Mon, 14 Sep 2026 13:02:10 -0400 Subject: [PATCH 02/16] Move partner Simulations off CoupledCellBuffer onto a new CoupledSimulation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CoupledCellBuffer previously carried a partner_sims::NamedTuple field so reinit_buffer! could reinitialize each partner's CellBuffer. Since CoupledCellBuffer is built once per task, this duplicated an identical copy of partner_sims into every task-local buffer of a threaded coupled domain. Rename CoupledMember to CoupledSimulation and give it a `partners` field (NamedTuple for a single-domain member, Dict{String,<:NamedTuple} per domain for a multi-domain member) holding the single, canonical copy of the resolved partner Simulations. CoupledCellBuffer now only holds partner_buffers (the CellBuffers themselves, still needed per task/per domain for the correct threaded binding); reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum) reads partner Simulations from `sim.partners` at call time instead. work! gained CoupledSimulation-specific work_domain_sequential!/ work_domain_threaded! methods (mirroring work.jl's plain-Simulation ones) so that reinit_buffer! receives the full CoupledSimulation, not just the inner plain Simulation; scatter-before-dispatch (the BUG-003 fix) is now driven directly from `partners` instead of a separately-tracked container tuple. Multi-domain members gained Base.iterate on CoupledSimulation, pairing each per-domain Simulation with its own partner slice. Pure internal-representation refactor: no public API or numerical-behavior change. Reviewed with Codex (dual-review skill) at the plan stage (one finding: the new reinit_buffer! method needed to live in Coupling.jl, after CoupledSimulation is defined, not in CoupledCellBuffer.jl which is included first — fixed) and against the final diff (no findings). Pkg.test() passes in full (CoupledSimulations: 1515/1515). phasefield_fracture.jl and the full docs/make.jl build both run clean. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- src/Coupling.jl | 225 +++++++++++++++++++++------ src/ItemBuffers/CoupledCellBuffer.jl | 30 ++-- test/coupled_simulations.jl | 10 +- 3 files changed, 190 insertions(+), 75 deletions(-) diff --git a/src/Coupling.jl b/src/Coupling.jl index 11e35576..0e060a8c 100644 --- a/src/Coupling.jl +++ b/src/Coupling.jl @@ -3,56 +3,177 @@ # A `CoupledSimulations` group is built once from a set of named `Simulation`s. Each # primary member's domain buffer(s) are rebuilt to hold a `CoupledCellBuffer` (or # `AutoDiffCellBuffer{<:CoupledCellBuffer}`) itembuffer that references the *actual* -# mutable buffers/simulations of its declared partners (primaries and refs alike, except -# itself). No coupling is resolved or discovered during `work!`. +# mutable buffers of its declared partners (primaries and refs alike, except itself). The +# partner *Simulation*s needed to reinitialize those buffers are not duplicated into every +# task-local `CoupledCellBuffer`; they live once on the `CoupledSimulation` handle and are +# passed down to `reinit_buffer!` at call time. No coupling is resolved or discovered during +# `work!`. """ - CoupledMember(sim, partner_containers) + CoupledSimulation(sim, partners) -A handle to one primary member of a [`CoupledSimulations`](@ref) group. `sim` is a -[`Simulation`](@ref) whose domain buffer(s) have been rebuilt with coupled itembuffers. -`partner_containers` are the raw itembuffer containers (`TaskLocals` or plain) of every -partner this member reads from, scattered once at the start of each `work!` call so that -threaded partners' task-local buffers observe the partner's current base state (e.g. its -time increment) even if the partner itself has not been `work!`ed since it last changed. +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 CoupledMember{S<:Simulation, PT<:Tuple} +struct CoupledSimulation{S<:Simulation, P} sim::S - partner_containers::PT + partners::P end -function Base.getproperty(m::CoupledMember, name::Symbol) - name === :sim && return getfield(m, :sim) - name === :partner_containers && return getfield(m, :partner_containers) - return getproperty(getfield(m, :sim), name) +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 -get_material(m::CoupledMember, args::Vararg{Any,N}) where N = get_material(getfield(m, :sim), args...) -get_dofhandler(m::CoupledMember) = get_dofhandler(getfield(m, :sim)) -get_grid(m::CoupledMember) = get_grid(getfield(m, :sim)) -get_state(m::CoupledMember, args::Vararg{Any,N}) where N = get_state(getfield(m, :sim), args...) -get_old_state(m::CoupledMember, args::Vararg{Any,N}) where N = get_old_state(getfield(m, :sim), args...) -getset(m::CoupledMember, args::Vararg{Any,N}) where N = getset(getfield(m, :sim), args...) -update_states!(m::CoupledMember; kwargs...) = update_states!(getfield(m, :sim); kwargs...) -set_time_increment!(m::CoupledMember, Δt) = set_time_increment!(getfield(m, :sim), Δt) -revert_states!(m::CoupledMember) = revert_states!(getfield(m, :sim)) -get_itembuffer(m::CoupledMember, args::Vararg{Any,N}) where N = get_itembuffer(getfield(m, :sim), args...) -get_num_tasks(m::CoupledMember) = get_num_tasks(getfield(m, :sim)) -get_chunks(m::CoupledMember) = get_chunks(getfield(m, :sim)) - -replace_material(::CoupledMember, args...; kwargs...) = throw(ArgumentError( +get_material(csim::CoupledSimulation, args::Vararg{Any,N}) where N = get_material(getfield(csim, :sim), args...) +get_dofhandler(csim::CoupledSimulation) = get_dofhandler(getfield(csim, :sim)) +get_grid(csim::CoupledSimulation) = get_grid(getfield(csim, :sim)) +get_state(csim::CoupledSimulation, args::Vararg{Any,N}) where N = get_state(getfield(csim, :sim), args...) +get_old_state(csim::CoupledSimulation, args::Vararg{Any,N}) where N = get_old_state(getfield(csim, :sim), args...) +getset(csim::CoupledSimulation, args::Vararg{Any,N}) where N = getset(getfield(csim, :sim), args...) +update_states!(csim::CoupledSimulation; kwargs...) = update_states!(getfield(csim, :sim); kwargs...) +set_time_increment!(csim::CoupledSimulation, Δt) = set_time_increment!(getfield(csim, :sim), Δt) +revert_states!(csim::CoupledSimulation) = revert_states!(getfield(csim, :sim)) +get_itembuffer(csim::CoupledSimulation, args::Vararg{Any,N}) where N = get_itembuffer(getfield(csim, :sim), args...) +get_num_tasks(csim::CoupledSimulation) = get_num_tasks(getfield(csim, :sim)) +get_chunks(csim::CoupledSimulation) = get_chunks(getfield(csim, :sim)) + +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.")) -_scatter_partner!(c::TaskLocals) = scatter!(c) -_scatter_partner!(::Any) = nothing +# Per-domain iteration for a multi-domain member, mirroring `Simulation{<:DomainBuffers}`'s +# own iteration but pairing each per-domain `Simulation` with its own slice of `partners`. +function Base.iterate(csim::CoupledSimulation{<:Simulation{<:DomainBuffers}}) + it = iterate(getfield(csim, :sim)) + it === nothing && return nothing + ((name, dsim), st) = it + return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st) +end +function Base.iterate(csim::CoupledSimulation{<:Simulation{<:DomainBuffers}}, st) + it = iterate(getfield(csim, :sim), st) + it === nothing && return nothing + ((name, dsim), st2) = it + return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st2) +end + +_scatter_partner_container!(c::TaskLocals) = scatter!(c) +_scatter_partner_container!(::Any) = nothing + +_flatten_partner_sims(partners::NamedTuple) = values(partners) +_flatten_partner_sims(partners_by_domain::Dict) = (psim for nt in values(partners_by_domain) for psim in values(nt)) + +# Scatter every reachable partner's task-local buffers from its base once, before any +# per-cell work, so a threaded reader always observes the partner's *current* state (e.g. its +# time increment) even if the partner itself has not been `work!`ed since it last changed. +function _scatter_all_partners!(csim::CoupledSimulation) + for psim in _flatten_partner_sims(getfield(csim, :partners)) + _scatter_partner_container!(get_itembuffer(psim.db)) + end + return nothing +end + +const CoupledSingleDomainSim = CoupledSimulation{<:SingleDomainSim} +const CoupledMultiDomainSim = CoupledSimulation{<:MultiDomainSim} +const CoupledSingleDomainThreadedSim = CoupledSimulation{<:SingleDomainThreadedSim} +const CoupledMultiDomainThreadedSim = CoupledSimulation{<:MultiDomainThreadedSim} + +function work!(worker, csim::CoupledMultiDomainSim) + _scatter_all_partners!(csim) + for (name, dcsim) in csim + skip_this_domain(worker, name) && continue + work_domain_sequential!(worker, dcsim) + end +end +function work!(worker, csim::CoupledSingleDomainSim) + _scatter_all_partners!(csim) + work_domain_sequential!(worker, csim) +end +function work!(worker, csim::CoupledMultiDomainThreadedSim) + _scatter_all_partners!(csim) + if can_thread(worker) + workers = TaskLocals(worker, num_tasks = get_num_tasks(csim)) + for (name, dcsim) in csim + skip_this_domain(worker, name) && continue + work_domain_threaded!(workers, dcsim) + end + else + for (name, dcsim) in csim + skip_this_domain(worker, name) && continue + work_domain_sequential!(worker, dcsim) + end + end +end +function work!(worker, csim::CoupledSingleDomainThreadedSim) + _scatter_all_partners!(csim) + if can_thread(worker) + workers = TaskLocals(worker; num_tasks = get_num_tasks(csim)) + work_domain_threaded!(workers, csim) + else + work_domain_sequential!(worker, csim) + end +end + +# Mirror `work.jl`'s plain-`Simulation` `work_domain_sequential!`/`work_domain_threaded!`, +# dispatching on `CoupledSimulation` instead so that `reinit_buffer!` receives the full +# `CoupledSimulation` (and thereby its `.partners`), not just the plain inner `Simulation`. +function work_domain_sequential!(worker, sim::CoupledSimulation) + itembuffer = get_base(get_itembuffer(sim)) + for itemnr in getset(sim) + reinit_buffer!(itembuffer, sim, itemnr) + work_single!(worker, itembuffer) + end +end -function work!(worker, m::CoupledMember, args...; kwargs...) - foreach(_scatter_partner!, getfield(m, :partner_containers)) - return work!(worker, getfield(m, :sim), args...; kwargs...) +function work_domain_threaded!(workers, sim::CoupledSimulation) + itembuffers = get_itembuffer(sim) #::TaskLocals + scatter!(itembuffers) + scatter!(workers) + num_tasks = get_num_tasks(sim) + 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 + +""" + 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, getfield(sim, :sim), cellnum) + reinit_partners!(cb.partner_buffers, getfield(sim, :partners), cellnum) + return nothing end struct CoupledSimulations{P<:NamedTuple, R<:NamedTuple, M<:NamedTuple} @@ -96,7 +217,7 @@ function CoupledSimulations(primaries::NamedTuple; refs::NamedTuple = NamedTuple validate_storage_identity(all_members) validate_task_counts_positive(all_members) members = NamedTuple{keys(primaries)}( - Tuple(build_coupled_member(name, sim, all_members) for (name, sim) in pairs(primaries)) + Tuple(build_coupled_simulation(name, sim, all_members) for (name, sim) in pairs(primaries)) ) return CoupledSimulations(primaries, refs, members) end @@ -121,11 +242,11 @@ 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, partner_sims::NamedTuple) +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, partner_sims)) : - CoupledCellBuffer(primary_cb, partners_nt, partner_sims) + 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)) @@ -159,16 +280,19 @@ function validate_domain_pair(reader_db::AbstractDomainBuffer, partner_db::Abstr 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) - partner_sims = map(p -> p.sim, partners) - coupled_ibuf = build_coupled_itembuffer(reader_ibuf, partner_containers, partner_sims) + coupled_ibuf = build_coupled_itembuffer(reader_ibuf, partner_containers) new_db = setproperties(reader_db; itembuffer = coupled_ibuf) - return new_db, values(partner_containers) + 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 @@ -192,40 +316,41 @@ function partner_domain_sim(dname::Union{Nothing,String}, partner_name::Symbol, end end -function build_coupled_member(name::Symbol, reader_sim::Simulation, all_members::NamedTuple) +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 - containers = Any[] + partners_by_domain = Dict{String, NamedTuple}() else - containers = Any[] 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, conts = build_coupled_domain(rdb, partners) + ndb, dpartner_sims = build_coupled_domain(rdb, partners) push!(built, dname => ndb) - append!(containers, conts) + 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, conts = build_coupled_domain(reader_db, partners) - containers = collect(conts) + 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 - new_sim = Simulation(new_db, reader_sim.a, reader_sim.aold) - return CoupledMember(new_sim, Tuple(unique(containers))) end _scratch_identity(cb::CellBuffer) = cb.ae # survives replace_material's setproperties (fields copied by reference) diff --git a/src/ItemBuffers/CoupledCellBuffer.jl b/src/ItemBuffers/CoupledCellBuffer.jl index 68fc9b86..306b10f3 100644 --- a/src/ItemBuffers/CoupledCellBuffer.jl +++ b/src/ItemBuffers/CoupledCellBuffer.jl @@ -1,20 +1,21 @@ """ - CoupledCellBuffer(primary::CellBuffer, partner_buffers::NamedTuple, partner_sims::NamedTuple) + 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)) -and the partner `Simulation`s required to reinitialize them (`partner_sims`, internal only). +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. +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, PS<:NamedTuple} <: AbstractCellBuffer +struct CoupledCellBuffer{CB<:CellBuffer, PB<:NamedTuple} <: AbstractCellBuffer primary::CB partner_buffers::PB - partner_sims::PS end for op = (:get_Ke, :get_re, :get_ae, :get_material, :get_values, :get_time_increment, @@ -35,19 +36,8 @@ for op = (:celldofs, :getcoordinates, :getfieldnames, :cellid) end Ferrite.dof_range(cb::CoupledCellBuffer, name::Symbol) = Ferrite.dof_range(cb.primary, name) -""" - reinit_buffer!(cb::CoupledCellBuffer, sim::Simulation, cellnum::Int) - -Reinitialize the reader's own `cb.primary` against `sim`, then reinitialize each -partner buffer in `cb.partner_buffers` against its own partner `Simulation` stored in -`cb.partner_sims`. Partner reinitialization does not recurse: partner buffers are plain -`CellBuffer`s, so no further coupling initialization happens. -""" -function reinit_buffer!(cb::CoupledCellBuffer, sim::Simulation, cellnum::Int) - reinit_buffer!(cb.primary, sim, cellnum) - reinit_partners!(cb.partner_buffers, cb.partner_sims, cellnum) - return nothing -end +# `reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum::Int)` is defined in +# Coupling.jl, after `CoupledSimulation` exists (this file is included before Coupling.jl). function reinit_partners!(buffers::NamedTuple, sims::NamedTuple, cellnum::Int) map((b, s) -> (reinit_buffer!(b, s, cellnum); nothing), buffers, sims) @@ -56,5 +46,5 @@ end function _replace_material_with(cb::CoupledCellBuffer, new_material) new_primary = _replace_material_with(cb.primary, new_material) - return CoupledCellBuffer(new_primary, cb.partner_buffers, cb.partner_sims) + return CoupledCellBuffer(new_primary, cb.partner_buffers) end diff --git a/test/coupled_simulations.jl b/test/coupled_simulations.jl index 508cbcbb..1b003afc 100644 --- a/test/coupled_simulations.jl +++ b/test/coupled_simulations.jl @@ -65,7 +65,7 @@ sim1 = Simulation(d1, a1, aold1) sim2 = Simulation(d2, a2, aold2) g = CoupledSimulations((a = sim1,); refs = (b = sim2,)) - @test g.a isa FerriteAssembly.CoupledMember + @test g.a isa FerriteAssembly.CoupledSimulation @test g.b === sim2 # refs are the plain source Simulation K = allocate_matrix(dh1) @@ -117,9 +117,9 @@ sim2 = Simulation(d2, a2, aold2) sim3 = Simulation(d3, a3, aold3) g = CoupledSimulations((a = sim1, b = sim2, c = sim3)) - @test g.a isa FerriteAssembly.CoupledMember - @test g.b isa FerriteAssembly.CoupledMember - @test g.c isa FerriteAssembly.CoupledMember + @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 @@ -201,7 +201,7 @@ 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.CoupledMember + @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 From e82c69e62f72d124e2db32da839fc989bfd5b584 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Mon, 14 Sep 2026 13:14:22 -0400 Subject: [PATCH 03/16] Rename rn/pn task-count locals to fix typos CI false positive crate-ci/typos flagged the short local `pn` as a misspelling of `on` in validate_domain_pair's task-count check. Renamed to reader_tasks/partner_tasks (also just more readable than rn/pn). No behavior change. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- src/Coupling.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Coupling.jl b/src/Coupling.jl index 0e060a8c..9f002b2b 100644 --- a/src/Coupling.jl +++ b/src/Coupling.jl @@ -270,12 +270,12 @@ function validate_domain_pair(reader_db::AbstractDomainBuffer, partner_db::Abstr "coupling partner `$partner_name` does not cover all cells read by the reader")) reader_threaded = reader_db isa ThreadedDomainBuffer if reader_threaded - rn = get_num_tasks(reader_db) - rn > 0 || throw(ArgumentError("task count must be positive")) - pn = partner_db isa ThreadedDomainBuffer ? get_num_tasks(partner_db) : 1 - rn == pn || throw(ArgumentError( - "threaded reader with $rn tasks requires coupling partner `$partner_name` to provide $rn " * - "task-local buffers (a sequential partner counts as 1 slot); got $pn")) + 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 From 446b8cf724db4964f4e4984c7be1e888add4ff97 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Mon, 14 Sep 2026 17:27:09 -0400 Subject: [PATCH 04/16] Unify work!/work_domain_*! dispatch for Simulation and CoupledSimulation Coupling.jl duplicated all six work!/work_domain_sequential!/ work_domain_threaded! methods from work.jl, each body identical to the plain-Simulation version except dispatching on CoupledSimulation instead. CoupledSimulation already forwards every accessor these functions call (get_itembuffer, getset, get_num_tasks, get_chunks, Base.iterate for multi-domain) to its wrapped Simulation, so the duplication added nothing. Replace the six duplicated methods with a single set in work.jl, dispatching on Union type aliases (AnySingleDomainSim, AnyMultiDomainSim, etc.) that cover both a plain Simulation and a CoupledSimulation wrapping the same domain-buffer shape. Each work! method gains one line, _prepare_work!(sim), a tiny extension point: a no-op for Any (defined in work.jl, so plain Simulation calls are unaffected) with a single override in Coupling.jl, _prepare_work!(csim::CoupledSimulation) = _scatter_all_partners!(csim), replacing the scatter call that used to be hand-inlined into every duplicated work! method. This requires Coupling.jl (which defines CoupledSimulation) to be included before work.jl (whose new Union aliases reference it); swapped their order in FerriteAssembly.jl. Coupling.jl calls nothing from work.jl, so this has no other effect. Pure dispatch-layer refactor: no behavior or public API change. Reviewed with Codex (dual-review skill) at the plan stage and against the final diff: no findings either time. Pkg.test() passes in full, unchanged pass counts (CoupledSimulations: 1515/1515). phasefield_fracture.jl and the full docs/make.jl build both run clean. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- src/Coupling.jl | 83 ++---------------------------------------- src/FerriteAssembly.jl | 2 +- src/work.jl | 30 ++++++++++++--- 3 files changed, 29 insertions(+), 86 deletions(-) diff --git a/src/Coupling.jl b/src/Coupling.jl index 9f002b2b..db8cad7f 100644 --- a/src/Coupling.jl +++ b/src/Coupling.jl @@ -82,85 +82,10 @@ function _scatter_all_partners!(csim::CoupledSimulation) return nothing end -const CoupledSingleDomainSim = CoupledSimulation{<:SingleDomainSim} -const CoupledMultiDomainSim = CoupledSimulation{<:MultiDomainSim} -const CoupledSingleDomainThreadedSim = CoupledSimulation{<:SingleDomainThreadedSim} -const CoupledMultiDomainThreadedSim = CoupledSimulation{<:MultiDomainThreadedSim} - -function work!(worker, csim::CoupledMultiDomainSim) - _scatter_all_partners!(csim) - for (name, dcsim) in csim - skip_this_domain(worker, name) && continue - work_domain_sequential!(worker, dcsim) - end -end -function work!(worker, csim::CoupledSingleDomainSim) - _scatter_all_partners!(csim) - work_domain_sequential!(worker, csim) -end -function work!(worker, csim::CoupledMultiDomainThreadedSim) - _scatter_all_partners!(csim) - if can_thread(worker) - workers = TaskLocals(worker, num_tasks = get_num_tasks(csim)) - for (name, dcsim) in csim - skip_this_domain(worker, name) && continue - work_domain_threaded!(workers, dcsim) - end - else - for (name, dcsim) in csim - skip_this_domain(worker, name) && continue - work_domain_sequential!(worker, dcsim) - end - end -end -function work!(worker, csim::CoupledSingleDomainThreadedSim) - _scatter_all_partners!(csim) - if can_thread(worker) - workers = TaskLocals(worker; num_tasks = get_num_tasks(csim)) - work_domain_threaded!(workers, csim) - else - work_domain_sequential!(worker, csim) - end -end - -# Mirror `work.jl`'s plain-`Simulation` `work_domain_sequential!`/`work_domain_threaded!`, -# dispatching on `CoupledSimulation` instead so that `reinit_buffer!` receives the full -# `CoupledSimulation` (and thereby its `.partners`), not just the plain inner `Simulation`. -function work_domain_sequential!(worker, sim::CoupledSimulation) - itembuffer = get_base(get_itembuffer(sim)) - for itemnr in getset(sim) - reinit_buffer!(itembuffer, sim, itemnr) - work_single!(worker, itembuffer) - end -end - -function work_domain_threaded!(workers, sim::CoupledSimulation) - itembuffers = get_itembuffer(sim) #::TaskLocals - scatter!(itembuffers) - scatter!(workers) - num_tasks = get_num_tasks(sim) - 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 +# Hook hit once at the start of every top-level `work!` call (see `work.jl`); a no-op for a +# plain `Simulation`, overridden here so a `CoupledSimulation`'s partners are scattered before +# any per-cell work, without `work.jl` needing to know coupling exists. +_prepare_work!(csim::CoupledSimulation) = _scatter_all_partners!(csim) """ reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum::Int) diff --git a/src/FerriteAssembly.jl b/src/FerriteAssembly.jl index cba8f6f4..c4618898 100644 --- a/src/FerriteAssembly.jl +++ b/src/FerriteAssembly.jl @@ -21,8 +21,8 @@ include("ItemBuffers/CoupledCellBuffer.jl") include("ItemBuffers/FacetBuffer.jl") include("Autodiff/autodiff.jl") -include("work.jl") include("Coupling.jl") +include("work.jl") include("Workers/Assemblers.jl") include("Workers/Integrators.jl") include("Workers/QuadPointEvaluator.jl") diff --git a/src/work.jl b/src/work.jl index bc73aa41..54bbce08 100644 --- a/src/work.jl +++ b/src/work.jl @@ -2,6 +2,20 @@ function work!(worker, buffer::Union{AbstractDomainBuffer, DomainBuffers}; a = n return work!(worker, Simulation(buffer, a, aold)) end +# Hit once at the start of every top-level `work!` call, before any per-cell work. A no-op for +# a plain `Simulation`; `Coupling.jl` overrides this for `CoupledSimulation` to scatter its +# partners' task-local buffers, without this file needing to know coupling exists. +_prepare_work!(::Any) = nothing + +const AnySingleDomainSim = Union{SingleDomainSim, CoupledSimulation{<:SingleDomainSim}} +const AnyMultiDomainSim = Union{MultiDomainSim, CoupledSimulation{<:MultiDomainSim}} +const AnySingleDomainThreadedSim = Union{SingleDomainThreadedSim, CoupledSimulation{<:SingleDomainThreadedSim}} +const AnyMultiDomainThreadedSim = Union{MultiDomainThreadedSim, CoupledSimulation{<:MultiDomainThreadedSim}} +# Deliberately broad, mirroring the two concrete cases `work_domain_sequential!` is called +# for: a genuinely sequential domain, and the sequential fallback for a threaded domain when +# the worker doesn't support threading. +const AnyDomainSim = Union{Simulation{<:AbstractDomainBuffer}, CoupledSimulation{<:Simulation{<:AbstractDomainBuffer}}} + """ work!(worker, sim::Simulation) @@ -18,16 +32,19 @@ Simplified interface, directly forwarded to `work!(worker, Simulation(db, a, aol 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) +function work!(worker, multisim::AnyMultiDomainSim) + _prepare_work!(multisim) for (name, sim) in multisim skip_this_domain(worker, name) && continue work_domain_sequential!(worker, sim) end end -function work!(worker, sim::SingleDomainSim) +function work!(worker, sim::AnySingleDomainSim) + _prepare_work!(sim) work_domain_sequential!(worker, sim) end -function work!(worker, multisim::MultiDomainThreadedSim) +function work!(worker, multisim::AnyMultiDomainThreadedSim) + _prepare_work!(multisim) if can_thread(worker) workers = TaskLocals(worker, num_tasks = get_num_tasks(multisim)) for (name, sim) in multisim @@ -41,7 +58,8 @@ function work!(worker, multisim::MultiDomainThreadedSim) end end end -function work!(worker, sim::SingleDomainThreadedSim) +function work!(worker, sim::AnySingleDomainThreadedSim) + _prepare_work!(sim) if can_thread(worker) workers = TaskLocals(worker; num_tasks = get_num_tasks(sim)) work_domain_threaded!(workers, sim) @@ -50,7 +68,7 @@ function work!(worker, sim::SingleDomainThreadedSim) end end -function work_domain_sequential!(worker, sim::Simulation{<:AbstractDomainBuffer}) +function work_domain_sequential!(worker, sim::AnyDomainSim) itembuffer = get_base(get_itembuffer(sim)) # get_base if threaded buffer for itemnr in getset(sim) reinit_buffer!(itembuffer, sim, itemnr) @@ -58,7 +76,7 @@ function work_domain_sequential!(worker, sim::Simulation{<:AbstractDomainBuffer} end end -function work_domain_threaded!(workers, sim::SingleDomainThreadedSim) +function work_domain_threaded!(workers, sim::AnySingleDomainThreadedSim) itembuffers = get_itembuffer(sim) #::TaskLocals scatter!(itembuffers) scatter!(workers) From 29d3e2ef4c8fe528d7fbf772897c31ed1d2479a1 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 15 Sep 2026 09:59:55 -0400 Subject: [PATCH 05/16] WIP introduce AbstractSimulation --- src/Coupling.jl | 351 ------------------------------------- src/FerriteAssembly.jl | 1 - src/Simulation.jl | 386 ++++++++++++++++++++++++++++++++++++++--- test/setup.jl | 2 +- 4 files changed, 364 insertions(+), 376 deletions(-) delete mode 100644 src/Coupling.jl diff --git a/src/Coupling.jl b/src/Coupling.jl deleted file mode 100644 index db8cad7f..00000000 --- a/src/Coupling.jl +++ /dev/null @@ -1,351 +0,0 @@ -# Explicit, setup-time coupling between simulations. -# -# A `CoupledSimulations` group is built once from a set of named `Simulation`s. Each -# primary member's domain buffer(s) are rebuilt to hold a `CoupledCellBuffer` (or -# `AutoDiffCellBuffer{<:CoupledCellBuffer}`) itembuffer that references the *actual* -# mutable buffers of its declared partners (primaries and refs alike, except itself). The -# partner *Simulation*s needed to reinitialize those buffers are not duplicated into every -# task-local `CoupledCellBuffer`; they live once on the `CoupledSimulation` handle and are -# passed down to `reinit_buffer!` at call time. No coupling is resolved or discovered during -# `work!`. - -""" - CoupledSimulation(sim, partners) - -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 CoupledSimulation{S<:Simulation, P} - 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 - -get_material(csim::CoupledSimulation, args::Vararg{Any,N}) where N = get_material(getfield(csim, :sim), args...) -get_dofhandler(csim::CoupledSimulation) = get_dofhandler(getfield(csim, :sim)) -get_grid(csim::CoupledSimulation) = get_grid(getfield(csim, :sim)) -get_state(csim::CoupledSimulation, args::Vararg{Any,N}) where N = get_state(getfield(csim, :sim), args...) -get_old_state(csim::CoupledSimulation, args::Vararg{Any,N}) where N = get_old_state(getfield(csim, :sim), args...) -getset(csim::CoupledSimulation, args::Vararg{Any,N}) where N = getset(getfield(csim, :sim), args...) -update_states!(csim::CoupledSimulation; kwargs...) = update_states!(getfield(csim, :sim); kwargs...) -set_time_increment!(csim::CoupledSimulation, Δt) = set_time_increment!(getfield(csim, :sim), Δt) -revert_states!(csim::CoupledSimulation) = revert_states!(getfield(csim, :sim)) -get_itembuffer(csim::CoupledSimulation, args::Vararg{Any,N}) where N = get_itembuffer(getfield(csim, :sim), args...) -get_num_tasks(csim::CoupledSimulation) = get_num_tasks(getfield(csim, :sim)) -get_chunks(csim::CoupledSimulation) = get_chunks(getfield(csim, :sim)) - -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.")) - -# Per-domain iteration for a multi-domain member, mirroring `Simulation{<:DomainBuffers}`'s -# own iteration but pairing each per-domain `Simulation` with its own slice of `partners`. -function Base.iterate(csim::CoupledSimulation{<:Simulation{<:DomainBuffers}}) - it = iterate(getfield(csim, :sim)) - it === nothing && return nothing - ((name, dsim), st) = it - return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st) -end -function Base.iterate(csim::CoupledSimulation{<:Simulation{<:DomainBuffers}}, st) - it = iterate(getfield(csim, :sim), st) - it === nothing && return nothing - ((name, dsim), st2) = it - return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st2) -end - -_scatter_partner_container!(c::TaskLocals) = scatter!(c) -_scatter_partner_container!(::Any) = nothing - -_flatten_partner_sims(partners::NamedTuple) = values(partners) -_flatten_partner_sims(partners_by_domain::Dict) = (psim for nt in values(partners_by_domain) for psim in values(nt)) - -# Scatter every reachable partner's task-local buffers from its base once, before any -# per-cell work, so a threaded reader always observes the partner's *current* state (e.g. its -# time increment) even if the partner itself has not been `work!`ed since it last changed. -function _scatter_all_partners!(csim::CoupledSimulation) - for psim in _flatten_partner_sims(getfield(csim, :partners)) - _scatter_partner_container!(get_itembuffer(psim.db)) - end - return nothing -end - -# Hook hit once at the start of every top-level `work!` call (see `work.jl`); a no-op for a -# plain `Simulation`, overridden here so a `CoupledSimulation`'s partners are scattered before -# any per-cell work, without `work.jl` needing to know coupling exists. -_prepare_work!(csim::CoupledSimulation) = _scatter_all_partners!(csim) - -""" - 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, getfield(sim, :sim), cellnum) - reinit_partners!(cb.partner_buffers, getfield(sim, :partners), cellnum) - 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 - -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/FerriteAssembly.jl b/src/FerriteAssembly.jl index c4618898..d69e57af 100644 --- a/src/FerriteAssembly.jl +++ b/src/FerriteAssembly.jl @@ -21,7 +21,6 @@ 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") diff --git a/src/Simulation.jl b/src/Simulation.jl index 5acee157..6c97d615 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -1,3 +1,35 @@ +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::Simulation) = get_num_tasks(get_domainbuffer(sim)) +get_chunks(sim::Simulation{<:AbstractDomainBuffer}) = get_chunks(get_domainbuffer(sim)) +get_itembuffer(sim::Simulation, args::Vararg{Any, N}) where {N} = get_itembuffer(get_domainbuffer(sim), args...) + + """ Simulation(db, a = nothing, aold = nothing) @@ -11,36 +43,15 @@ struct Simulation{ 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}} - -# 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...) +get_domainbuffer(sim::Simulation) = sim.db -# 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 @@ -50,3 +61,332 @@ end Base.iterate(sim::Simulation{<:DomainBuffers}) = _iterate(sim, iterate(sim.db)) Base.iterate(sim::Simulation{<:DomainBuffers}, iter) = _iterate(sim, iterate(sim.db, iter)) +""" + CoupledSimulation(sim, partners) + +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 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 + +get_domainbuffer(sim::CoupledSimulation) = get_domainbuffer(getfield(sim, :sim)) + +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.")) + +# Per-domain iteration for a multi-domain member, mirroring `Simulation{<:DomainBuffers}`'s +# own iteration but pairing each per-domain `Simulation` with its own slice of `partners`. +function Base.iterate(csim::CoupledSimulation{<:Simulation{<:DomainBuffers}}) + it = iterate(getfield(csim, :sim)) + it === nothing && return nothing + ((name, dsim), st) = it + return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st) +end +function Base.iterate(csim::CoupledSimulation{<:Simulation{<:DomainBuffers}}, st) + it = iterate(getfield(csim, :sim), st) + it === nothing && return nothing + ((name, dsim), st2) = it + return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st2) +end + +_scatter_partner_container!(c::TaskLocals) = scatter!(c) +_scatter_partner_container!(::Any) = nothing + +_flatten_partner_sims(partners::NamedTuple) = values(partners) +_flatten_partner_sims(partners_by_domain::Dict) = (psim for nt in values(partners_by_domain) for psim in values(nt)) + +# Scatter every reachable partner's task-local buffers from its base once, before any +# per-cell work, so a threaded reader always observes the partner's *current* state (e.g. its +# time increment) even if the partner itself has not been `work!`ed since it last changed. +function _scatter_all_partners!(csim::CoupledSimulation) + for psim in _flatten_partner_sims(getfield(csim, :partners)) + _scatter_partner_container!(get_itembuffer(psim.db)) + end + return nothing +end + +# Hook hit once at the start of every top-level `work!` call (see `work.jl`); a no-op for a +# plain `Simulation`, overridden here so a `CoupledSimulation`'s partners are scattered before +# any per-cell work, without `work.jl` needing to know coupling exists. +_prepare_work!(csim::CoupledSimulation) = _scatter_all_partners!(csim) + +""" + 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, getfield(sim, :sim), cellnum) + reinit_partners!(cb.partner_buffers, getfield(sim, :partners), cellnum) + 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 + +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/test/setup.jl b/test/setup.jl index 53a4e596..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 From d113939056b59981eb26379e8fd48bbd739fd0d1 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Tue, 15 Sep 2026 10:15:12 -0400 Subject: [PATCH 06/16] Finish AbstractSimulation: split CoupledSimulation from the CellBuffer-dependent coupling machinery The prior WIP commit (29d3e2ef) introduced AbstractSimulation{DB} and moved all of Coupling.jl's content into Simulation.jl, but didn't precompile: Simulation.jl sits early in the include order (before CellBuffer.jl, CoupledCellBuffer.jl, Autodiff/autodiff.jl), while the moved-in coupling machinery (reinit_buffer! for CoupledCellBuffer, build_coupled_itembuffer, validate_domain_pair, the CoupledSimulations group constructor, ...) references CellBuffer/CoupledCellBuffer/AutoDiffCellBuffer, none of which exist yet at that point. work.jl (untouched by the WIP) also broke, since it still needed the concrete SingleDomainSim/MultiDomainSim/... aliases the WIP had removed from Simulation.jl. Clarified with the user: keep AbstractSimulation + CoupledSimulation's core (struct, accessors, iteration, the partner-scatter _prepare_work! hook) in Simulation.jl, since none of that needs cell-buffer types. Move only the genuinely CellBuffer-dependent coupling machinery (reinit_buffer! for CoupledCellBuffer, the CoupledSimulations struct/constructor/validation, build_coupled_itembuffer/build_coupled_domain/build_coupled_simulation, partner_domain_sim, the group-level replace_material) into a restored Coupling.jl, included after Autodiff/autodiff.jl and before work.jl, same position it held before the WIP. This also lets work.jl's dispatch simplify beyond the previous round's Union-based approach: since Simulation and CoupledSimulation share AbstractSimulation{DB} with the same DB parameter, work! and work_domain_sequential!/work_domain_threaded! now dispatch directly on AbstractSingleDomainSim/AbstractMultiDomainSim/.../AbstractSimulation{<:AbstractDomainBuffer} bounds - no Union{Simulation{<:X}, CoupledSimulation{<:Simulation{<:X}}} plumbing needed at all. Also finished generalizing get_num_tasks/get_chunks/ get_itembuffer from ::Simulation-typed to ::AbstractSimulation-typed, which the WIP had done for the other nine forwarding methods but not these three. Fixed one bug caught by Codex's plan review: CoupledSimulation's per-domain Base.iterate methods still matched on the old 2-type-param layout (CoupledSimulation{<:Simulation{<:DomainBuffers}}); with the new DB-first layout the correct bound is CoupledSimulation{<:DomainBuffers}, or a multi-domain coupled primary's work! throws MethodError on iterate. Pkg.test() passes in full, unchanged pass counts (CoupledSimulations: 1515/1515). phasefield_fracture.jl and the full docs/make.jl build both run clean. Reviewed with Codex (dual-review skill) at the plan stage (1 finding, fixed as above) and against the final diff (no findings). Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- src/Coupling.jl | 269 ++++++++++++++++++++++++++++++++++++++ src/FerriteAssembly.jl | 1 + src/Simulation.jl | 285 ++--------------------------------------- src/work.jl | 21 +-- 4 files changed, 287 insertions(+), 289 deletions(-) create mode 100644 src/Coupling.jl diff --git a/src/Coupling.jl b/src/Coupling.jl new file mode 100644 index 00000000..fa3d0a5e --- /dev/null +++ b/src/Coupling.jl @@ -0,0 +1,269 @@ +# CoupledSimulations group construction and the cell-buffer-dependent parts of coupling. +# +# CoupledSimulation itself (struct, accessors, iteration, partner-scatter hook) 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, getfield(sim, :sim), cellnum) + reinit_partners!(cb.partner_buffers, getfield(sim, :partners), cellnum) + 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 + +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/FerriteAssembly.jl b/src/FerriteAssembly.jl index d69e57af..c4618898 100644 --- a/src/FerriteAssembly.jl +++ b/src/FerriteAssembly.jl @@ -21,6 +21,7 @@ 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") diff --git a/src/Simulation.jl b/src/Simulation.jl index 6c97d615..c2d7f8cb 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -25,23 +25,23 @@ set_time_increment!(sim::AbstractSimulation, Δt) = set_time_increment!(get_doma revert_states!(sim::AbstractSimulation) = revert_states!(get_domainbuffer(sim)) # Forwarding for internal API -get_num_tasks(sim::Simulation) = get_num_tasks(get_domainbuffer(sim)) -get_chunks(sim::Simulation{<:AbstractDomainBuffer}) = get_chunks(get_domainbuffer(sim)) -get_itembuffer(sim::Simulation, args::Vararg{Any, N}) where {N} = get_itembuffer(get_domainbuffer(sim), args...) +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 @@ -55,7 +55,7 @@ get_domainbuffer(sim::Simulation) = sim.db ## 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)) @@ -94,13 +94,13 @@ replace_material(::CoupledSimulation, args...; kwargs...) = throw(ArgumentError( # Per-domain iteration for a multi-domain member, mirroring `Simulation{<:DomainBuffers}`'s # own iteration but pairing each per-domain `Simulation` with its own slice of `partners`. -function Base.iterate(csim::CoupledSimulation{<:Simulation{<:DomainBuffers}}) +function Base.iterate(csim::CoupledSimulation{<:DomainBuffers}) it = iterate(getfield(csim, :sim)) it === nothing && return nothing ((name, dsim), st) = it return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st) end -function Base.iterate(csim::CoupledSimulation{<:Simulation{<:DomainBuffers}}, st) +function Base.iterate(csim::CoupledSimulation{<:DomainBuffers}, st) it = iterate(getfield(csim, :sim), st) it === nothing && return nothing ((name, dsim), st2) = it @@ -127,266 +127,3 @@ end # plain `Simulation`, overridden here so a `CoupledSimulation`'s partners are scattered before # any per-cell work, without `work.jl` needing to know coupling exists. _prepare_work!(csim::CoupledSimulation) = _scatter_all_partners!(csim) - -""" - 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, getfield(sim, :sim), cellnum) - reinit_partners!(cb.partner_buffers, getfield(sim, :partners), cellnum) - 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 - -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/work.jl b/src/work.jl index 54bbce08..f6babd0d 100644 --- a/src/work.jl +++ b/src/work.jl @@ -7,15 +7,6 @@ end # partners' task-local buffers, without this file needing to know coupling exists. _prepare_work!(::Any) = nothing -const AnySingleDomainSim = Union{SingleDomainSim, CoupledSimulation{<:SingleDomainSim}} -const AnyMultiDomainSim = Union{MultiDomainSim, CoupledSimulation{<:MultiDomainSim}} -const AnySingleDomainThreadedSim = Union{SingleDomainThreadedSim, CoupledSimulation{<:SingleDomainThreadedSim}} -const AnyMultiDomainThreadedSim = Union{MultiDomainThreadedSim, CoupledSimulation{<:MultiDomainThreadedSim}} -# Deliberately broad, mirroring the two concrete cases `work_domain_sequential!` is called -# for: a genuinely sequential domain, and the sequential fallback for a threaded domain when -# the worker doesn't support threading. -const AnyDomainSim = Union{Simulation{<:AbstractDomainBuffer}, CoupledSimulation{<:Simulation{<:AbstractDomainBuffer}}} - """ work!(worker, sim::Simulation) @@ -32,18 +23,18 @@ Simplified interface, directly forwarded to `work!(worker, Simulation(db, a, aol 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::AnyMultiDomainSim) +function work!(worker, multisim::AbstractMultiDomainSim) _prepare_work!(multisim) for (name, sim) in multisim skip_this_domain(worker, name) && continue work_domain_sequential!(worker, sim) end end -function work!(worker, sim::AnySingleDomainSim) +function work!(worker, sim::AbstractSingleDomainSim) _prepare_work!(sim) work_domain_sequential!(worker, sim) end -function work!(worker, multisim::AnyMultiDomainThreadedSim) +function work!(worker, multisim::AbstractMultiDomainThreadedSim) _prepare_work!(multisim) if can_thread(worker) workers = TaskLocals(worker, num_tasks = get_num_tasks(multisim)) @@ -58,7 +49,7 @@ function work!(worker, multisim::AnyMultiDomainThreadedSim) end end end -function work!(worker, sim::AnySingleDomainThreadedSim) +function work!(worker, sim::AbstractSingleDomainThreadedSim) _prepare_work!(sim) if can_thread(worker) workers = TaskLocals(worker; num_tasks = get_num_tasks(sim)) @@ -68,7 +59,7 @@ function work!(worker, sim::AnySingleDomainThreadedSim) end end -function work_domain_sequential!(worker, sim::AnyDomainSim) +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) @@ -76,7 +67,7 @@ function work_domain_sequential!(worker, sim::AnyDomainSim) end end -function work_domain_threaded!(workers, sim::AnySingleDomainThreadedSim) +function work_domain_threaded!(workers, sim::AbstractSingleDomainThreadedSim) itembuffers = get_itembuffer(sim) #::TaskLocals scatter!(itembuffers) scatter!(workers) From 2995d4764436f1d16a353c0f449bca46aa11c627 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 15 Sep 2026 10:49:41 -0400 Subject: [PATCH 07/16] Remove custom accessors for CoupledSimulation --- src/Simulation.jl | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/src/Simulation.jl b/src/Simulation.jl index c2d7f8cb..e869cb70 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -80,13 +80,7 @@ struct CoupledSimulation{DB, S <: Simulation{DB}, P} <: AbstractSimulation{DB} 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 - -get_domainbuffer(sim::CoupledSimulation) = get_domainbuffer(getfield(sim, :sim)) +get_domainbuffer(sim::CoupledSimulation) = get_domainbuffer(sim.sim) replace_material(::CoupledSimulation, args...; kwargs...) = throw(ArgumentError( "replace_material on a CoupledSimulations member is not supported; use " * @@ -95,16 +89,16 @@ replace_material(::CoupledSimulation, args...; kwargs...) = throw(ArgumentError( # Per-domain iteration for a multi-domain member, mirroring `Simulation{<:DomainBuffers}`'s # own iteration but pairing each per-domain `Simulation` with its own slice of `partners`. function Base.iterate(csim::CoupledSimulation{<:DomainBuffers}) - it = iterate(getfield(csim, :sim)) + it = iterate(csim.sim) it === nothing && return nothing ((name, dsim), st) = it - return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st) + return ((name, CoupledSimulation(dsim, csim.partners[name])), st) end function Base.iterate(csim::CoupledSimulation{<:DomainBuffers}, st) - it = iterate(getfield(csim, :sim), st) + it = iterate(csim.sim, st) it === nothing && return nothing ((name, dsim), st2) = it - return ((name, CoupledSimulation(dsim, getfield(csim, :partners)[name])), st2) + return ((name, CoupledSimulation(dsim, csim.partners[name])), st2) end _scatter_partner_container!(c::TaskLocals) = scatter!(c) @@ -117,7 +111,7 @@ _flatten_partner_sims(partners_by_domain::Dict) = (psim for nt in values(partner # per-cell work, so a threaded reader always observes the partner's *current* state (e.g. its # time increment) even if the partner itself has not been `work!`ed since it last changed. function _scatter_all_partners!(csim::CoupledSimulation) - for psim in _flatten_partner_sims(getfield(csim, :partners)) + for psim in _flatten_partner_sims(csim.partners) _scatter_partner_container!(get_itembuffer(psim.db)) end return nothing From e5675deffd3d84197a12693556cfcecf378f0943 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 15 Sep 2026 10:55:09 -0400 Subject: [PATCH 08/16] Simplify iteration of CoupledSimulation (multidomain) --- src/Simulation.jl | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/src/Simulation.jl b/src/Simulation.jl index e869cb70..c2580a77 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -86,20 +86,13 @@ 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.")) -# Per-domain iteration for a multi-domain member, mirroring `Simulation{<:DomainBuffers}`'s -# own iteration but pairing each per-domain `Simulation` with its own slice of `partners`. -function Base.iterate(csim::CoupledSimulation{<:DomainBuffers}) - it = iterate(csim.sim) - it === nothing && return nothing - ((name, dsim), st) = it - return ((name, CoupledSimulation(dsim, csim.partners[name])), st) -end -function Base.iterate(csim::CoupledSimulation{<:DomainBuffers}, st) - it = iterate(csim.sim, st) - it === nothing && return nothing - ((name, dsim), st2) = it - return ((name, CoupledSimulation(dsim, csim.partners[name])), st2) +@inline function _iterate(csim::CoupledSimulation{<:DomainBuffers}, iter) + iter === nothing && return nothing + ((name, sim), state) = iter + return ((name, CoupledSimulation(sim, csim.partners[name])), state) end +Base.iterate(sim::CoupledSimulation{<:DomainBuffers}) = _iterate(sim, iterate(sim.sim)) +Base.iterate(sim::CoupledSimulation{<:DomainBuffers}, iter) = _iterate(sim, iterate(sim.sim, iter)) _scatter_partner_container!(c::TaskLocals) = scatter!(c) _scatter_partner_container!(::Any) = nothing From 55577bdba635be92626fa4d05e038307500aab99 Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 15 Sep 2026 11:07:17 -0400 Subject: [PATCH 09/16] Remove _prepare_work! and implement scatter(AbstractThreadedSim) --- src/Coupling.jl | 4 +++- src/ItemBuffers/CoupledCellBuffer.jl | 8 -------- src/Simulation.jl | 28 +++++++--------------------- src/work.jl | 11 +---------- 4 files changed, 11 insertions(+), 40 deletions(-) diff --git a/src/Coupling.jl b/src/Coupling.jl index fa3d0a5e..b9472cc3 100644 --- a/src/Coupling.jl +++ b/src/Coupling.jl @@ -15,7 +15,9 @@ further coupling initialization happens. """ function reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum::Int) reinit_buffer!(cb.primary, getfield(sim, :sim), cellnum) - reinit_partners!(cb.partner_buffers, getfield(sim, :partners), cellnum) + map(cb.partner_buffers, sim.partners) do (b, s) + reinit_buffer!(b, s, cellnum) + end return nothing end diff --git a/src/ItemBuffers/CoupledCellBuffer.jl b/src/ItemBuffers/CoupledCellBuffer.jl index 306b10f3..f0ee9fd2 100644 --- a/src/ItemBuffers/CoupledCellBuffer.jl +++ b/src/ItemBuffers/CoupledCellBuffer.jl @@ -36,14 +36,6 @@ for op = (:celldofs, :getcoordinates, :getfieldnames, :cellid) end Ferrite.dof_range(cb::CoupledCellBuffer, name::Symbol) = Ferrite.dof_range(cb.primary, name) -# `reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum::Int)` is defined in -# Coupling.jl, after `CoupledSimulation` exists (this file is included before Coupling.jl). - -function reinit_partners!(buffers::NamedTuple, sims::NamedTuple, cellnum::Int) - map((b, s) -> (reinit_buffer!(b, s, cellnum); nothing), buffers, sims) - return nothing -end - function _replace_material_with(cb::CoupledCellBuffer, new_material) new_primary = _replace_material_with(cb.primary, new_material) return CoupledCellBuffer(new_primary, cb.partner_buffers) diff --git a/src/Simulation.jl b/src/Simulation.jl index c2580a77..5ab17d5d 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -61,6 +61,8 @@ end Base.iterate(sim::Simulation{<:DomainBuffers}) = _iterate(sim, iterate(sim.db)) Base.iterate(sim::Simulation{<:DomainBuffers}, iter) = _iterate(sim, iterate(sim.db, iter)) +scatter!(sim::Simulation{<:ThreadedDomainBuffer}) = scatter!(get_itembuffer(sim)) + """ CoupledSimulation(sim, partners) @@ -82,6 +84,11 @@ end 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.")) @@ -93,24 +100,3 @@ replace_material(::CoupledSimulation, args...; kwargs...) = throw(ArgumentError( end Base.iterate(sim::CoupledSimulation{<:DomainBuffers}) = _iterate(sim, iterate(sim.sim)) Base.iterate(sim::CoupledSimulation{<:DomainBuffers}, iter) = _iterate(sim, iterate(sim.sim, iter)) - -_scatter_partner_container!(c::TaskLocals) = scatter!(c) -_scatter_partner_container!(::Any) = nothing - -_flatten_partner_sims(partners::NamedTuple) = values(partners) -_flatten_partner_sims(partners_by_domain::Dict) = (psim for nt in values(partners_by_domain) for psim in values(nt)) - -# Scatter every reachable partner's task-local buffers from its base once, before any -# per-cell work, so a threaded reader always observes the partner's *current* state (e.g. its -# time increment) even if the partner itself has not been `work!`ed since it last changed. -function _scatter_all_partners!(csim::CoupledSimulation) - for psim in _flatten_partner_sims(csim.partners) - _scatter_partner_container!(get_itembuffer(psim.db)) - end - return nothing -end - -# Hook hit once at the start of every top-level `work!` call (see `work.jl`); a no-op for a -# plain `Simulation`, overridden here so a `CoupledSimulation`'s partners are scattered before -# any per-cell work, without `work.jl` needing to know coupling exists. -_prepare_work!(csim::CoupledSimulation) = _scatter_all_partners!(csim) diff --git a/src/work.jl b/src/work.jl index f6babd0d..f58fbf37 100644 --- a/src/work.jl +++ b/src/work.jl @@ -2,11 +2,6 @@ function work!(worker, buffer::Union{AbstractDomainBuffer, DomainBuffers}; a = n return work!(worker, Simulation(buffer, a, aold)) end -# Hit once at the start of every top-level `work!` call, before any per-cell work. A no-op for -# a plain `Simulation`; `Coupling.jl` overrides this for `CoupledSimulation` to scatter its -# partners' task-local buffers, without this file needing to know coupling exists. -_prepare_work!(::Any) = nothing - """ work!(worker, sim::Simulation) @@ -24,18 +19,15 @@ The global degree of freedom vectors, `a` and `aold`, make their corresponding l available. If not passed, the local values are `NaN`s. """ function work!(worker, multisim::AbstractMultiDomainSim) - _prepare_work!(multisim) for (name, sim) in multisim skip_this_domain(worker, name) && continue work_domain_sequential!(worker, sim) end end function work!(worker, sim::AbstractSingleDomainSim) - _prepare_work!(sim) work_domain_sequential!(worker, sim) end function work!(worker, multisim::AbstractMultiDomainThreadedSim) - _prepare_work!(multisim) if can_thread(worker) workers = TaskLocals(worker, num_tasks = get_num_tasks(multisim)) for (name, sim) in multisim @@ -50,7 +42,6 @@ function work!(worker, multisim::AbstractMultiDomainThreadedSim) end end function work!(worker, sim::AbstractSingleDomainThreadedSim) - _prepare_work!(sim) if can_thread(worker) workers = TaskLocals(worker; num_tasks = get_num_tasks(sim)) work_domain_threaded!(workers, sim) @@ -68,8 +59,8 @@ function work_domain_sequential!(worker, sim::AbstractSimulation{<:AbstractDomai end function work_domain_threaded!(workers, sim::AbstractSingleDomainThreadedSim) + scatter!(sim) # Includes scatter of the `itembuffers` 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) From 3bf72a1e23fc3701f3b419ddc650e21c00f8731d Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Tue, 15 Sep 2026 11:08:49 -0400 Subject: [PATCH 10/16] Fix map in reinit buffer --- src/Coupling.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Coupling.jl b/src/Coupling.jl index b9472cc3..5724f52f 100644 --- a/src/Coupling.jl +++ b/src/Coupling.jl @@ -15,9 +15,7 @@ further coupling initialization happens. """ function reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum::Int) reinit_buffer!(cb.primary, getfield(sim, :sim), cellnum) - map(cb.partner_buffers, sim.partners) do (b, s) - reinit_buffer!(b, s, cellnum) - end + map((b, s) -> (reinit_buffer!(b, s, cellnum); nothing), cb.partner_buffers, sim.partners) return nothing end From 5e34dcf9763fcaa55893ba8d0b0695257bf88492 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Tue, 15 Sep 2026 11:25:25 -0400 Subject: [PATCH 11/16] Fix two regressions from the property-forwarding/scatter! simplification commits Reviewed 4 follow-up commits (2995d47, e5675de, 55577bd, 3bf72a1) that simplify CoupledSimulation. e5675de (restyle multi-domain iteration) and 3bf72a1 (fix a map(...) do (b, s) bug: that syntax binds one tuple-destructured argument, but map(f, coll1, coll2) calls the block with two separate positional args per iteration) are correct as-is. Found and fixed two real regressions in the other two, confirmed by direct REPL reproduction (Pkg.test() and a docs build both passed despite the bugs, since the only code path that exercises them - the fracture tutorial's `.a`/`.aold` dot access and its solve() call - sits inside a markdown-only comment block, never executed by `include()`): - 2995d47 removed CoupledSimulation's custom Base.getproperty, which forwarded unknown property names (`.a`, `.aold`, `.db`) to the wrapped `.sim`. `.sim`/ `.partners` still worked (real fields), but `g.a.a` now threw a FieldError - breaking the exact usage pattern the fracture tutorial documents and the CoupledSimulation docstring still (falsely) claimed to support. Restored the forwarding override. - 55577bd's scatter!(sim::CoupledSimulation{<:ThreadedDomainBuffer}) calls map(scatter!, sim.partners), but only defined scatter! for Simulation{<:ThreadedDomainBuffer} - no fallback for a non-threaded partner. validate_domain_pair explicitly allows a threaded reader with exactly 1 task to pair with a sequential partner (a sequential partner counts as 1 slot); work! on such a group threw MethodError. Added a generic scatter!(::AbstractSimulation) = nothing fallback (the ThreadedDomainBuffer- bound methods stay and take precedence). Also reworded Coupling.jl's now-stale "partner-scatter hook" comment (the _prepare_work! hook these commits removed) and unified reinit_buffer!'s getfield(sim, :sim)/sim.partners mix to plain dot-access, now that the restored getproperty override makes that the safe, consistent form throughout the file. Added regression tests for both bugs to test/coupled_simulations.jl. Reviewed with Codex (dual-review skill) at the plan stage (1 finding: a proposed test assertion needed to compare `.db` against the *rebuilt* coupled buffer, not the original source `Simulation`'s `.db` - fixed) and against the final diff (no findings). Pkg.test() passes in full (CoupledSimulations: 1571/1571, up from 1515 with the new regression tests). phasefield_fracture.jl and the full docs/make.jl build both run clean. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- src/Coupling.jl | 5 +++-- src/Simulation.jl | 7 +++++++ test/coupled_simulations.jl | 25 +++++++++++++++++++++++++ 3 files changed, 35 insertions(+), 2 deletions(-) diff --git a/src/Coupling.jl b/src/Coupling.jl index 5724f52f..ca185588 100644 --- a/src/Coupling.jl +++ b/src/Coupling.jl @@ -1,6 +1,7 @@ # CoupledSimulations group construction and the cell-buffer-dependent parts of coupling. # -# CoupledSimulation itself (struct, accessors, iteration, partner-scatter hook) lives in +# 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). @@ -14,7 +15,7 @@ Partner reinitialization does not recurse: partner buffers are plain `CellBuffer further coupling initialization happens. """ function reinit_buffer!(cb::CoupledCellBuffer, sim::CoupledSimulation, cellnum::Int) - reinit_buffer!(cb.primary, getfield(sim, :sim), cellnum) + reinit_buffer!(cb.primary, sim.sim, cellnum) map((b, s) -> (reinit_buffer!(b, s, cellnum); nothing), cb.partner_buffers, sim.partners) return nothing end diff --git a/src/Simulation.jl b/src/Simulation.jl index 5ab17d5d..58a1bb9e 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -61,6 +61,7 @@ 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)) """ @@ -82,6 +83,12 @@ struct CoupledSimulation{DB, S <: Simulation{DB}, P} <: AbstractSimulation{DB} 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 + get_domainbuffer(sim::CoupledSimulation) = get_domainbuffer(sim.sim) function scatter!(sim::CoupledSimulation{<:ThreadedDomainBuffer}) diff --git a/test/coupled_simulations.jl b/test/coupled_simulations.jl index 1b003afc..06994689 100644 --- a/test/coupled_simulations.jl +++ b/test/coupled_simulations.jl @@ -68,6 +68,14 @@ @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 + K = allocate_matrix(dh1) r = zeros(ndofs(dh1)) assembler = start_assemble(K, r) @@ -99,6 +107,23 @@ @test nalloc < 2_000_000 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 From d94cb7f1241f90376f4eee5f9973c894d842ec87 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Tue, 15 Sep 2026 13:35:11 -0400 Subject: [PATCH 12/16] Add tab-completion and constprop tests for CoupledSimulation getproperty Add Base.propertynames(csim::CoupledSimulation) so the forwarded Simulation properties (.a, .aold, .db) appear alongside the real fields (.sim, .partners) during tab-completion, mirroring what CoupledSimulations already had. Add test coverage: - is_concrete_inferred helper using Base.return_types on a closure with the literal dot-access baked in, since @inferred cannot check constant propagation through property-name dispatch (it only accepts call expressions and infers from argument runtime types, not the compile-time constant property symbol). - Constprop assertions for .a/.aold/.db/.sim/.partners on CoupledSimulation and .a on the CoupledSimulations group, confirming each resolves to a single concrete type rather than a Union across getproperty branches. - propertynames assertions confirming forwarded properties tab-complete. Regression tests for the two bugs fixed in 5e34dcf (lost property forwarding, missing scatter! fallback) already existed from that commit; no additional tests needed there. Validation: Pkg.test() 1659/1659 passing (up from 1571 before these test additions), phasefield_fracture.jl tutorial runs cleanly, docs/make.jl build succeeds. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- src/Simulation.jl | 3 +++ test/coupled_simulations.jl | 28 ++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+) diff --git a/src/Simulation.jl b/src/Simulation.jl index 58a1bb9e..13f42f44 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -89,6 +89,9 @@ function Base.getproperty(csim::CoupledSimulation, name::Symbol) 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}) diff --git a/test/coupled_simulations.jl b/test/coupled_simulations.jl index 06994689..cd3500fb 100644 --- a/test/coupled_simulations.jl +++ b/test/coupled_simulations.jl @@ -1,4 +1,14 @@ @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"))) @@ -76,6 +86,24 @@ @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. + pn = propertynames(g.a) + @test :a in pn + @test :aold in pn + @test :db in pn + @test :sim in pn + @test :partners in pn + K = allocate_matrix(dh1) r = zeros(ndofs(dh1)) assembler = start_assemble(K, r) From 04bf975a18e58a1754584fbb8dec551ae3259451 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Tue, 15 Sep 2026 13:40:15 -0400 Subject: [PATCH 13/16] Rename propertynames test variable to avoid spellcheck false-positive The Spell Check CI job flagged \`pn\` as a likely typo for \`on\`. Renamed to \`propnames\`. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- test/coupled_simulations.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/coupled_simulations.jl b/test/coupled_simulations.jl index cd3500fb..415a4495 100644 --- a/test/coupled_simulations.jl +++ b/test/coupled_simulations.jl @@ -97,12 +97,12 @@ @test is_concrete_inferred(grp -> grp.a, typeof(g)) # Forwarded properties must tab-complete, not just the two real struct fields. - pn = propertynames(g.a) - @test :a in pn - @test :aold in pn - @test :db in pn - @test :sim in pn - @test :partners in pn + propnames = propertynames(g.a) + @test :a in propnames + @test :aold in propnames + @test :db in propnames + @test :sim in propnames + @test :partners in propnames K = allocate_matrix(dh1) r = zeros(ndofs(dh1)) From bf17b0a8aa262c9b3210fd5bc6588b005ab60370 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Tue, 15 Sep 2026 14:35:56 -0400 Subject: [PATCH 14/16] Full-branch dual-review of PR92 coupling redesign; fix 3 test/validation gaps A complete review of the whole cb/coupled-simulations-redesign branch (main...HEAD), not just the latest commits, using the dual-review skill with Codex as independent reviewer (2 rounds: plan review, final-diff review). Fixes found: - src/Coupling.jl: unwrap_cb had no fallback method, so coupling an unsupported buffer kind (e.g. FacetBuffer) threw a raw MethodError from storage-identity validation (which runs before the actionable validate_domain_pair check, and even for a sole primary with no partners) instead of an ArgumentError. Added a fallback that throws an actionable ArgumentError. - test/coupled_simulations.jl: - The 3-member mutual-coupling test only checked construction, never called work!, so a positional mix-up between two partners' buffers/simulations could pass undetected. Now works g.a (sequentially and threaded) with distinguishable per-partner values for both partners. - The allocation-scaling smoke test only checked one 4-cell mesh against a fixed ceiling, which couldn't detect a small per-cell allocation. Added a dedicated testset comparing a 2x2 vs 20x20 grid against an uncoupled baseline, asserting exactly zero coupling-specific overhead (caught, in review, that the first version of this fix itself allocated a slice per cell, defeating its purpose). - Added facet-buffer coupling rejection tests (sole primary and partner). Also checked off BUG-017 in identified_bugs.md ("Coupling setup incompletely validates domain/task compatibility"): this branch's setup-time validation (domain key matching, task-count compatibility with actionable errors) is exactly the fix that bug called for. Test results: Pkg.test() all green with 1 and 4 threads (CoupledSimulations: 1724/1724, up from 1659). phasefield_fracture.jl tutorial and full docs/make.jl build (every tutorial/how-to) both run clean, no errors. --- src/Coupling.jl | 2 + test/coupled_simulations.jl | 139 ++++++++++++++++++++++++++++++++---- 2 files changed, 126 insertions(+), 15 deletions(-) diff --git a/src/Coupling.jl b/src/Coupling.jl index ca185588..29f6a96b 100644 --- a/src/Coupling.jl +++ b/src/Coupling.jl @@ -79,6 +79,8 @@ Base.propertynames(cs::CoupledSimulations) = (:primaries, :refs, :members, keys( 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 diff --git a/test/coupled_simulations.jl b/test/coupled_simulations.jl index 415a4495..f336a7df 100644 --- a/test/coupled_simulations.jl +++ b/test/coupled_simulations.jl @@ -44,6 +44,12 @@ @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 @@ -135,6 +141,82 @@ @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 @@ -161,22 +243,37 @@ ip3 = Lagrange{RefQuadrilateral,1}() dh3 = close!(add!(DofHandler(grid), :w, ip3)) cv3 = CellValues(qr, ip3, ip3) - a3 = zeros(ndofs(dh3)) + # 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)) - d1 = setup_domainbuffer(DomainSpec(dh1, CS_MA(), cvu); a = a1) - d2 = setup_domainbuffer(DomainSpec(dh2, CS_MB(), cvv); a = a2) - d3 = setup_domainbuffer(DomainSpec(dh3, CS_MC(), cv3); a = a3) - 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)}) + 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 @@ -310,5 +407,17 @@ # 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 From e0e7c4a0b47b947589391a6755445224f0f1dddd Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Tue, 15 Sep 2026 16:13:44 -0400 Subject: [PATCH 15/16] Add tests for PR92 coverage gaps identified by Codecov Codecov flagged 13 missed lines on this branch (project coverage 96.94% -> 96.43%, patch 94.77%), all in coupling code added/restructured by this PR: src/Coupling.jl (getproperty/propertynames on the CoupledSimulations group, the empty multi-domain reader branch of build_coupled_simulation, and the domain-selector path of replace_material), src/ItemBuffers/CoupledCellBuffer.jl (dof_range forwarding and _replace_material_with), and src/work.jl (the can_thread/skip_this_domain generic trait-default fallbacks, which every worker in the package overrides everywhere else). Added regression tests for each (test/coupled_simulations.jl, test/assemblers.jl), via the dual-review skill with Codex as independent reviewer: - CoupledSimulations group getproperty (unknown-member ArgumentError) and propertynames. - replace_material(group, member, f; domain=...) success and non-dict-member error paths. - CoupledCellBuffer's dof_range forwarding and _replace_material_with (calling replace_material directly on a member's `.db`, framed as an internal-plumbing check, not a supported user workflow, since the result isn't re-workable via work! without the partner Simulations that live on the owning CoupledSimulation). - A primary member with a genuinely empty multi-domain Dict (construction-only edge case). - can_thread/skip_this_domain defaults for a worker with no trait overrides. Left Coupling.jl:88 (_is_autodiff(::AutoDiffCellBuffer)) untested: reasoned, and Codex's plan review confirmed via a runtime probe, that it is already dispatched to by the existing sequential-autodiff coupling test - Codecov's miss there looks like a coverage-tool artifact for a trivial one-line method, not an untested path. Codex plan review caught that the initial domain-selector fixture mixed a multi-domain primary with a single-domain ref, which fails coupling's shape validation before reaching either intended branch; fixed by using two matching multi-domain members for the success path and a separate plain group for the error path. It also flagged that two comments overstated support for edge cases the tests don't actually exercise via work!; reworded both to state only what's verified. Final-diff review: no findings. Test results: Pkg.test() all green (CoupledSimulations 1784/1784, up from 1724; Assemblers 23/23, up from 22; full suite passes). docs/make.jl builds clean. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- test/assemblers.jl | 9 +++++ test/coupled_simulations.jl | 75 +++++++++++++++++++++++++++++++++++++ 2 files changed, 84 insertions(+) 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 index f336a7df..b7848600 100644 --- a/test/coupled_simulations.jl +++ b/test/coupled_simulations.jl @@ -110,6 +110,16 @@ @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) @@ -368,6 +378,71 @@ @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) From 367a73ca8f17e4ec6c47aa37b42f19caaa64a9ba Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Thu, 17 Sep 2026 14:10:32 +0200 Subject: [PATCH 16/16] Apply batched suggestions from code review Co-authored-by: Knut Andreas Meyer --- docs/src/literate_tutorials/phasefield_fracture.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/literate_tutorials/phasefield_fracture.jl b/docs/src/literate_tutorials/phasefield_fracture.jl index 2b0a146b..780041ec 100644 --- a/docs/src/literate_tutorials/phasefield_fracture.jl +++ b/docs/src/literate_tutorials/phasefield_fracture.jl @@ -158,9 +158,10 @@ db_d_uc, Kd, rd, ndofs_d = setup(PhaseFieldFracture{:d}(mbase), grid, :d; ip_quad = Lagrange{RefQuadrilateral, 2}() ) -sim_u0 = Simulation(db_u_uc, zeros(ndofs_u), zeros(ndofs_u)) -sim_d0 = Simulation(db_d_uc, zeros(ndofs_d), zeros(ndofs_d)) -g = CoupledSimulations((u = sim_u0, d = sim_d0)) # `:u` and `:d` mutually read each other +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 @@ -179,8 +180,7 @@ function get_reaction_dofs(dh) end; # ## Solving -# Write function to solve one simulation part. The coupling to the other part is already -# wired into `sim` (a `CoupledSimulations` member handle), so `work!` needs no extra input. +# 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))