Add full-DOF matrix-free action entry point#307
Conversation
Introduces assemble_matrix_free_action_full!(asm, func, U_full, v_full, p),
a sibling of the existing free-DOF assemble_matrix_free_action! that
accepts full-DOF U_full and v_full vectors instead of the free-DOF Uu
and Vu. Both vectors must have length equal to p.field.data.
The existing free-DOF entry point writes Vu into p.hvp_scratch_field's
free slots only, leaving the BC slots implicitly zero. This is the
right semantics for Krylov-style J·v calls, where v is a perturbation
ranging over unknowns and the BC component genuinely is zero. It is
the wrong semantics for residual-style evaluations of c_M·M·dU in
implicit Newmark with time-varying Dirichlet BCs: there, dU has a
nonzero BC slice equal to g(t_{n+1}) − U_pred_BC, and the silently
dropped M_{f,BC}·dU_BC cross-term biases the inertial residual at
boundary-adjacent free rows.
The full-DOF entry point gives the caller a direct way to inject the
BC contribution: copyto! into p.field.data and p.hvp_scratch_field.data
bypasses _update_for_assembly!'s "scatter free into unknowns, write BC
from bc_cache.vals" logic; the element kernel then sees v_el including
the BC component and produces the cross-term naturally. This mirrors
Norma.jl's full-DOF-state design without changing how Carina-style
free-DOF integrators store U/V/A elsewhere.
Carefully restores the contract for the free-DOF entry point: after
the action loop, BC slots of p.hvp_scratch_field are zeroed out. Were
this skipped, a subsequent free-DOF Krylov J·v call would silently
read leftover v_BC values and return a contaminated action.
Test (TestAssemblers.jl) validates three invariants against a sibling
no-BC assembler that produces the unrestricted element-scatter K_full:
(a) On free rows, the full-DOF action equals (K_full · v_full)[free]
for v_full with nonzero BC slots (rtol=1e-6 accommodates sparse
SpMV vs per-element scatter floating-point reordering).
(b) Zeroing v_full's BC slots collapses the full-DOF entry point
onto the free-DOF entry point bit-for-bit.
(c) The two answers in (a) and (b) genuinely differ — i.e., the
cross-term is actually carrying signal on this mesh, so the
validation is non-trivial.
The condensed-assembler approach normally used elsewhere does NOT
work as a reference: FEC's sparsity-pattern construction
(SparsityPatterns.jl:200-249) drops entries whose row OR column is a
Dirichlet DOF, even in condensed mode. Removing the BCs from the
assembler entirely is the only way to retain K_{f,BC}.
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #307 +/- ##
==========================================
+ Coverage 65.59% 65.77% +0.18%
==========================================
Files 56 56
Lines 4868 4894 +26
==========================================
+ Hits 3193 3219 +26
Misses 1675 1675 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
The two failure modes here are pre-existing on 1. "CI" matrix (GitHub-hosted runners, Julia 1.10/1.12 × macOS/Ubuntu/Windows) — the Julia test step itself passes on every runner ( Same failure on main commit The self-hosted runners (which presumably take a different path that doesn't hit the broken codecov-action step) pass on both this PR and main. 2.
None of these are touched by this PR — Happy to send a separate small PR adding the missing My new |
I'm fixing up the codecov in #308 and should get around to the docs in a subsequent PR. Thank for this PR @lxmota this is a subtlety I had not thought about carefully before. |
2 participants
Summary
Adds
assemble_matrix_free_action_full!(asm, func, U_full, v_full, p), a sibling of the existingassemble_matrix_free_action!that accepts full-DOFU_fullandv_fullvectors instead of free-DOFUu/Vu. Both vectors must have lengthlength(p.field.data).Why this primitive
The existing free-DOF entry point writes
Vuintop.hvp_scratch_field's free slots only, leaving the BC slots implicitly zero. That's correct for Krylov-styleJ·vcalls —vis a perturbation ranging over unknowns and the BC component genuinely is zero.It is wrong for residual-style evaluations like
c_M·M·dUin implicit Newmark with time-varying Dirichlet BCs: there,dUhas a nonzero BC slicedU_BC = g(t_{n+1}) − U_pred_BC, and the silently-droppedM_{f,BC}·dU_BCcross-term biases the inertial residual at boundary-adjacent free rows. (The downstream Carina.jl integrator currently has this exact bug; the fix is a small Carina-side change once this primitive lands.)The new entry point gives the caller a direct way to inject the BC contribution:
copyto!intop.field.dataandp.hvp_scratch_field.databypasses_update_for_assembly!'s "scatter free into unknowns, write BC frombc_cache.vals" logic; the element kernel then seesv_elincluding the BC component and produces the cross-term naturally. This mirrors Norma.jl's full-DOF-state design without changing how Carina-style free-DOF integrators store U/V/A elsewhere.Carefully restoring the contract
After the action loop, BC slots of
p.hvp_scratch_fieldare zeroed out. Were this skipped, a subsequent free-DOF KrylovJ·vcall would silently read leftoverv_BCvalues and return a contaminated action — turning what is a useful new entry point into a tripwire for anyone alternating the two paths.