Explicit setup-time coupling for CoupledSimulations (replaces couple_buffers)

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
 ```

docs/src/literate_tutorials/phasefield_fracture.jl

@@ -158,8 +158,11 @@ db_d_uc, Kd, rd, ndofs_d = setup(PhaseFieldFracture{:d}(mbase), grid, :d;
     ip_quad = Lagrange{RefQuadrilateral, 2}()
     )
 
-sim_u = Simulation(couple_buffers(db_u_uc; d = db_d_uc), zeros(ndofs_u), zeros(ndofs_u))
-sim_d = Simulation(couple_buffers(db_d_uc; u = db_u_uc), zeros(ndofs_d), zeros(ndofs_d));
+g = CoupledSimulations((
+    u = Simulation(db_u_uc, zeros(ndofs_u), zeros(ndofs_u)), 
+    d = Simulation(db_d_uc, zeros(ndofs_d), zeros(ndofs_d)),
+    ))
+sim_u, sim_d = g.u, g.d;
 
 # Setup loading and boundary conditions
 load_function(t) = 1e-4 * t
@@ -177,16 +180,16 @@ function get_reaction_dofs(dh)
 end;
 
 # ## Solving
-# Write function to solve one simulation part, given the other as input.
-function solve_single_part(sim, coupled, K, r, ch; firsttol = 1e-5, tol = 1e-6, maxiter = 100)
+# Common function to solve each simulation
+function solve_single_part(sim, K, r, ch; firsttol = 1e-5, tol = 1e-6, maxiter = 100)
     if ch !== nothing # Displacement part
         reaction_dofs = get_reaction_dofs(FerriteAssembly.get_dofhandler(sim))
     else
         reaction_dofs = Int[]
     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

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,18 +57,14 @@ 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
         return AutoDiffCellBuffer(cb)
     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