Skip to content

Bounded coupling design: after sequential thermal mechanics

Status: bounded 3D experimental Step provider implemented; smooth-mode space/time convergence and local installed-wheel serial/MPI acceptance pass. PR #90 is merged after its final-head CI passed. Version 0.4.1 publication remains gated by release acceptance; general coupled validation is not claimed. See evidence/coupling/2026-10-09-installed-wheel.md for the exact candidate identity.

Ordinary workflow (experimental)

Use studies.coupled_thermoelastic() with an SI Model, register displacement and a temperature-departure field, then call:

step = model.step(
    target=u, temperature_departure=theta, dt=0.1, steps=3,
    heat_loads=(heating,), mechanical_loads=(), output="fields.xdmf",
)
result = step.solve_result()

See examples/coupled_thermoelastic_3d.py for a complete ordinary workflow. procedures.staggered_implicit_euler() describes the numerical route; the Study describes shared-mesh 3D small-strain thermoelastic physics. No numerical loop was added to Model. Each registered natural load must be assigned to an explicit participant. Existing nested strong constraints from model.fix() are reused, including their component spaces. Regional materials, additional eigenstrains, weak/contact boundaries, changing natural loads and non-SI units are rejected in this slice. Output is the final accepted snapshot at its real physical time. progress=True uses the existing throttled reporter; default console cadence is every 100 accepted windows plus the final window, with the reporter's 30-second heartbeat between accepted-window notifications. It does not spawn a background thread or retain a second event list. Long individual linear solves do not yet emit within-solve heartbeats. checkpoint=checkpointing.every(...) publishes only after joint acceptance; explicit keep_last bounds retained files generated by that continuation segment. Explicit save/load remains available between windows. Custom history and transient field-series options are not accepted yet.

The Step option contract is shared by capability discovery and execution. Manual advance() inherits the same configured iteration policy as run(); an explicit per-window override does not modify subsequent windows. Restart identity includes actual initial field contents, mesh, individual boundary membership and values, histories, material and executable natural load forms. Changing a load therefore cannot reuse an old checkpoint merely because its display name stayed the same. The result remains experimental and computed; this interface does not promote scientific maturity.

Implemented experiment

Independent nonuniform decay convergence

tools/benchmark_coupled_decay.py exercises the ordinary Step, not the private cube oracle. On a unit box, all normal displacements vanish and all thermal faces are insulated. With H=lambda+2*mu, beta=(3*lambda+2*mu)*alpha and w=2*pi, the independently derived continuous solution is

theta = cos(w*x) * exp(-gamma*t)
u_x = beta/(H*w) * sin(w*x) * exp(-gamma*t), u_y=u_z=0
gamma = k*w^2 / (C + beta^2*T0/H)

Indeed H*du_x/dx=beta*theta satisfies axial equilibrium and the normal constraints; substituting into the linearized heat equation gives gamma. Nonzero initial fields are interpolated explicitly. The reference follows the linearized equations cited below; it is our analytical verification problem, not a reproduction of the cited perforated-plate geometry.

The 2026-10-09 local Q1-hexahedron run kept one cell in each transverse direction (the exact solution is transverse-uniform). At 128 time steps, 8/16/32 axial cells give temperature relative L2 errors 0.0712221 / 0.0177099 / 0.00404966, and displacement errors 0.117874 / 0.0299607 / 0.00713390. Observed spatial orders are 2.008, 2.129 and 1.976, 2.070. At 128 axial cells, 4/8/16 time steps give temperature errors 0.0168227 / 0.00849863 / 0.00417228 and displacement errors 0.0166185 / 0.00829613 / 0.00397068; temporal orders are 0.985, 1.026 and 1.002, 1.063. The fixed acceptance requires both spatial orders above 1.7 and both temporal orders above 0.8.

Each participant matrix assembles once in all six cases. Maximum quadratic residual is below 5.3e-14 J. This is a bounded smooth-mode verification and reuse check, not extreme-scale performance, multidimensional mesh convergence, or a general first-law certificate. Reproduce with:

PYTHONPATH=src python tools/benchmark_coupled_decay.py --output /tmp/coupled-decay.json

operators.thermoelastic_heat_source(...) contributes explicit reversible volumetric heat feedback for constant coefficients on a shared 3D mesh. It does not enable coupling automatically. A separate private benchmark driver compares staggered implicit Euler with independently written mixed-space UFL forms, both executed by the existing prepared PETSc solver. No second solver backend, global Model state or nested transient-step clock was introduced.

PYTHONPATH=src python tools/benchmark_thermoelastic_coupling.py --nonuniform --output outputs/coupling-prototype.json
PYTHONPATH=src agentfem mpi-run -n 2 --timeout 120 -- python tools/benchmark_thermoelastic_coupling.py --alpha 0.01 --relaxation 0.25

These are repository development commands, not installed-user templates. The private class under benchmarks/thermoelastic_coupling.py is not a stable public Procedure. The result remains computed; successful execution does not claim general coupled validation.

The experiment uses an insulated unit cube with symmetry rollers, homogeneous isotropic elasticity, constant-strain volumetric capacity and no external mechanical work. The uniform-source oracle is independent:

theta(t) = Q*t / (rho*c_epsilon + 9*K*alpha^2*T0)
u(x,t) = alpha*theta(t)*x

A nonuniform source additionally activates conduction. The comparison tests two mesh/time resolutions, fixed relaxation, zero coupling, strong-coupling nonconvergence, single-rank acceptance rejection and retry. They test agreement with the discrete reference, not a full mesh/time convergence certificate. Residuals are unrelaxed volume-normalized FE RMS norms, with each field accepted against atol + rtol*norm(field). Both criteria must pass. Exactly zero coupling needs only one solve of each independent block.

Constant thermal and elastic matrices are each assembled once and reused through outer iterations and accepted windows. Solver execution counters retain failed attempts; accepted time/history and fields do not. No output is published from rejected windows. Non-finite global diagnostics and rank-inconsistent iteration controls are rejected explicitly.

Energy evidence: deliberately separate identities

For insulated boundaries, the volume integral of the linearized heat equation checks rho*c_epsilon*delta_theta + beta*T0*div(delta_u) - dt*Q. This is not a nonlinear total-energy closure claim.

With homogeneous mechanical boundary work, the prototype also verifies the discrete quadratic identity. Write a(u,v) for elastic virtual work and let

E_n = a(u_n,u_n)/2 + integral(rho*c_epsilon*theta_n^2/(2*T0))
D_dt = a(delta_u,delta_u)/2 + integral(rho*c_epsilon*delta_theta^2/(2*T0))
E_n - E_old + D_dt + dt*integral(k*|grad(theta_n)|^2/T0)
    = dt*integral(Q*theta_n/T0)

The time-discretization term is not material dissipation. This positive quadratic balance follows by testing the mechanical equation with delta_u and the thermal equation with dt*theta_n/T0; opposite coupling terms cancel. Keep it separate from physical total internal energy and boundary heating.

One local two-cell-per-axis nonuniform run (three dt=0.1 windows, alpha=0.002) used 9 outer iterations per window and one assembly per participant; maximum temperature/displacement reference errors were approximately 2.2e-12 K and 4.4e-15 m. These are RMS differences between two discrete solves, not physical prediction errors or a platform speedup. Reproduce with the command above.

Current acceptance and remaining gates

Ordinary registered boundaries, physical-only joint restart, the bounded Step option contract, and smooth-mode space/time convergence are implemented. tests/portable_coupled_checkpoint_driver.py additionally passes public-provider 1 -> 2 -> 1 continuation with moving strong boundaries. At every station it compares displacement, temperature, both force stations, and numeric accepted history against uninterrupted execution on the current partition. This driver is included in the installed-wheel release workflow. Local acceptance used 24 targeted serial tests and nine public-provider tests on each of two ranks.

PR #90 passed remote FEniCSx, macOS installed-wheel, lint, and documentation checks at 008d419d; subsequent fixes require their own check status before integration. The original installed-wheel record remains bound to 9cd7e644, not silently relabelled as evidence for later code. Release integration and any broader multidimensional/external scientific promotion remain separate gates. Neither the smooth mode nor the quadratic balance establishes a general first-law certificate.

Boundary and recovery increment (2026-10-08)

PR #89 was merged after remote FEniCSx and macOS installed-wheel acceptance. The follow-on private experiment now accepts uniform inward boundary heat flux (positive heating), prescribed boundary temperature rate, or insulated faces; temperature and flux cannot be prescribed on the same complete boundary. It also accepts a prescribed affine dilation rate on the upper cube faces. These options are benchmark inputs, not new application-specific public APIs.

Prescribed-temperature heat is extracted from the unconstrained thermal residual at owned boundary degrees of freedom, including capacity and coupling terms. Ghost contributions are accumulated before reduction; they are not counted twice. Flux heat is separately integrated on the boundary. The temperature-weighted residual supply belongs to the quadratic identity, not to the physical incoming heat record.

The virtual displacement v=x gives the generalized dilation reaction from the mechanical weak residual. Accepted old/new reactions and the prescribed dilation increment give trapezoidal path work. Backward Euler's quadratic identity instead uses endpoint work; both are recorded with distinct names. The prescribed-dilation oracle independently checks theta=(Q-3*beta*T0*dilation_rate)*t/C and reaction=3*(3*K*dilation_rate*t-beta*theta) on the unit cube. Zero-coupling uniform heating also checks prescribed-temperature reactions. These references do not establish arbitrary-boundary first-law closure.

Joint checkpoints reuse agentfem.transient-checkpoint.v5: both accepted fields, independent reference fields, common time, history digest, and all benchmark physics inputs are bound together. Trial-window writes are refused. Restore stages fields before checking accepted history and before touching live fields or time. Altered physics or damaged history is rejected. The independent reference fields are verification overhead and will not belong to the production Procedure. Serial-to-two-rank-to-serial recovery was exercised; repeat after archive-contract changes using:

PYTHONPATH=src python tools/benchmark_thermoelastic_coupling.py --dilation-rate .0001 --inward-heat-flux 2 --stop-after 1 --checkpoint /tmp/coupled-one
PYTHONPATH=src agentfem mpi-run -n 2 --timeout 120 -- python tools/benchmark_thermoelastic_coupling.py --dilation-rate .0001 --inward-heat-flux 2 --resume /tmp/coupled-one --stop-after 2 --checkpoint /tmp/coupled-two
PYTHONPATH=src python tools/benchmark_thermoelastic_coupling.py --dilation-rate .0001 --inward-heat-flux 2 --resume /tmp/coupled-two

Next: separate the reusable physical participant/Procedure from this oracle, retain ordinary region/constraint/load semantics, add its Step option contract, then perform space/time convergence and installed-use acceptance. Do not promote the benchmark class itself or duplicate its cube geometry in core.

The first extraction is private time/_staggered.py: ordered prepared solves, MPI-global unrelaxed field criteria and fixed relaxation, with residual forms compiled once. It owns scratch iterates only; accepted time, rollback and checkpoint publication remain with the calling Procedure. It has no cube, material, thermal boundary or independent reference solve. This is an internal implementation seam, not yet a public coupled Step. A regression prohibits creating new residual forms during iteration; no end-to-end speedup is claimed.

The physical blocks are separately lowered by private operators.thermoelastic._thermoelastic_blocks: constant homogeneous material, current/old displacement and temperature departure, time increment and explicit integrated loads. They have no geometry factory, boundary selection, solver or reference problem. The monolithic oracle continues to write its equations independently. The seam rejects incompatible history spaces and loads on the wrong test space; departure from T0 is explicit, never guessed from a field name. The remaining public Step must own boundary/time-input handling and the joint accepted lifecycle rather than subclassing the benchmark.

Private time/thermoelastic.py now exercises that separation on rectangular 3D tetrahedral and hexahedral meshes: two physical prepared systems, one accepted time, atomic rejection/retry and immutable returned records. It does not construct a reference solution. Its fixed-input joint checkpoint contains only physical accepted fields and history; the caller must supply a reviewed scientific input identity, in addition to the material and mesh identity. History is integrity-checked before live state changes. This is not a built-in Step provider yet; public execution/output policy remains a gate.

The private Procedure now consumes existing TimeDependentDirichlet assets, updates them at the proposed physical time, and restores their constants on rejection. Restore re-evaluates the same registered histories at accepted time. No arbitrary operator-changing callback is accepted. Generic evidence uses uneliminated thermal/mechanical residuals on owned strong-boundary DOFs, with reverse ghost accumulation, instead of geometry-specific reaction formulas. Natural-load and prescribed-motion path work retain both accepted force stations; their force history is part of the joint checkpoint. Energy forms are compiled once. The result reports separate cumulative heat and path work, linearized heat residual and quadratic residual; trust remains computed. Numerical inputs/results in this private slice use SI. Non-SI Model metadata, ordinary output policies and explicit capability rejection need review before public registration.

The initial-state preflight now rejects unsupported constraints, boundaries attached to the wrong participant space, non-finite initial values, and initial fields inconsistent with the prescribed values at physical time zero. It checks each boundary independently, so overlapping conflicting values cannot be hidden by assignment order. This is not an equilibrium certificate: users still own the physical consistency of the initial state. No automatic repair or hidden initial jump is introduced. A failed preflight restores registered boundary constants and runs before prepared solver construction. The targeted Procedure suite passes 13 tests in serial and on each of two MPI ranks, including a non-finite initial value injected on only one rank.

Public lowering must keep both participants explicit: displacement as the Step target and a named temperature-departure field, with no field-name inference. Its option contract must distinguish integrated mechanical and heat loads, participant-specific strong boundaries, and outer iteration controls. Before registering this route, bind registered model inputs to the restart identity and reuse the common result/output/progress policy. Reject unsupported regional materials, weak/contact constraints, operator-changing inputs, and incompatible unit systems rather than ignoring them. Public registration remains pending these checks and convergence acceptance.

Private solve_result() now delegates final-field publication and checkpoint contract attachment to the shared Result lifecycle. It advances only remaining windows, records the executed segment, and labels field output explicitly as a final accepted snapshot, not a transient series. Performance uses the shared rank-reduced timing schema with solve and result/output stages; constructor/JIT preparation is outside that call's timing scope. A rejected window publishes no completion dataset. Finite-strain OutputPlan is rejected before advancing; automatic progress/checkpoint policies now follow accepted windows; field-series cadence remains separate work. The initial public provider exposes this bounded completion; earlier private-stage notes above describe the extraction history rather than a separate user-facing solver.

Product decision

Finish the shared-mesh sequential route first. Keep Model as the problem definition, Operator as the physical residual contribution, Procedure as the iteration owner, State as the accepted/trial owner, and Result as evidence. Do not introduce a new multiphysics Model or a thermal-wall-specific solver.

The first two-way verification slice is small-strain, homogeneous linear thermoelasticity on one mesh with a reversible volumetric-strain-rate heat term. It is selected for an independent linear-system reference and explicit energy identity, not because it covers industrial forming. No plastic-work heat fraction, contact heating, arbitrary mesh transfer, or nonlinear material history belongs in this first slice.

For temperature departure theta = T - T0 and constant isotropic coefficients:

sigma = C:epsilon(u) - beta*theta*I
rho*c*theta_dot - div(k*grad(theta)) = r - beta*T0*div(u_dot)
beta = 3*K*alpha, K = bulk modulus

Here c is the heat capacity consistent with the chosen strain-based linear thermoelastic free energy. Derive units, signs, boundary work and entropy/ energy identities before coding the thermal feedback Operator. Do not insert this term into existing one-way heat models automatically. The 2D assumption must be explicitly derived; start with a 3D manufactured reference to avoid silently using a plane-stress reduction for a plane-strain model.

Numerical route and ownership

  1. Save the accepted time-window state of both participants.
  2. Solve thermal then mechanical subproblems at the same proposed physical time, with a bounded outer fixed-point iteration and optional fixed under-relaxation. Keep the accepted old-time state fixed during every retry.
  3. Measure the unrelaxed update/residual of each exchanged field, using MPI-global owned coefficients or physical FE norms. Report dimensional absolute and scaled relative residuals separately. A tiny relaxation factor must not create false convergence.
  4. Require every declared field criterion and the inner solver criteria. An inner Newton convergence flag does not imply outer coupling convergence. Near-zero reference fields require an absolute criterion, not division by an arbitrary tiny denominator. Tolerances and field scales are user-visible.
  5. Commit both states/time only after outer acceptance; on failure restore both, discard trial output, and report an explicit exhausted-iteration status. Time cutback, if enabled, restarts from the same accepted boundary.

Sequential EngineeringStep configuration remains declarative. An eventual coupled Procedure owns this loop; state.field_transaction alone cannot roll back a participant's completed time index or constitutive history. Before promoting the loop, audit participant snapshot/restore for all time, history, callbacks and operator invalidation, not just displacement/temperature.

With constant coefficients and dt, coupling updates change the RHS, allowing prepared subproblem reuse. A changed conductivity, elastic tangent, dt or constraint membership follows the existing operator invalidation contract. Timing and assembly counters must demonstrate reuse under unchanged physics.

Acceptance and stop conditions

  • Zero coupling reduces to the existing independently verified heat/solid solutions; no needless outer iteration is required.
  • Compare to an independently assembled monolithic block system for a tiny linear manufactured problem, including sign and unit checks.
  • Verify first-law accounting with imposed motion and thermal boundary input; do not equate mechanical strain energy alone to a thermodynamic balance.
  • Compare dt and mesh refinement, two rank counts, iteration tolerances and acceleration disabled/enabled. Under-relaxation must change convergence speed, not the converged solution.
  • Force one-rank failure and outer nonconvergence: neither participant may advance accepted time; restart/retry must match an uninterrupted solution.
  • Record outer residual history, inner solve counts, accepted time windows, rollback/cutback counts and immutable input references in existing results.

Timebox the first prototype to these gates. If participant rollback requires a broad rewrite, stop and publish the gap rather than adding a callback loop that merely appears coupled. Aitken/quasi-Newton acceleration follows a correct fixed-point baseline; a second backend and nonmatching meshes remain deferred.

References and design rationale

These references inform the contract. This plan neither adds a preCICE dependency nor claims that borrowing its iteration pattern verifies physics.

The zero-coupling regression also exposed two general edge cases: UFL erasing the argument of an exactly zero eigenstrain/projection form, and stale PETSc solution norms after direct NumPy edits. Both now have targeted regressions. The prepared-solve boundary increments vector state before KSP, following PETSc object-state semantics; it does not invalidate the matrix or its factorization.