Skip to content

Partial derivatives: per-callback declarations and forward accumulation (#217) - #227

Open
petlenz wants to merge 4 commits into
mainfrom
feature/issue-217
Open

petlenz wants to merge 4 commits into
mainfrom
feature/issue-217

Conversation

@petlenz

@petlenz petlenz commented Oct 28, 2026 •

Copy link
Copy Markdown
Member

Closes #217 (epic #221). Merge after #223 (feature/graph-export): this branch is stacked on it and uses its typed edges and the upstream_of() cache.

What and why

Partial derivatives are declared per callback, next to the reads, and the engine assembles requested derivatives by forward accumulation over the graph. Hosts ask for a derivative the way they ask for an output. Opt-in: no material changes behaviour, and nothing is computed unless a derivative is requested.

base::add_callback(&linear_elasticity::update_stress, m_sig)
    .reads(m_C, m_eps, m_eigenstrain)
    .partials(m_sig, m_C, m_eps, m_eigenstrain);

void update_stress() {
  ...                                   // values unchanged
  if (base::want_partials()) {
    base::partial(m_sig, m_eps) = m_C;  // tensor<T, 3, 4>&: rank(out) + rank(read)
    ...
  }
}

auto const& ddsddt{point.derivative<tensor2>("elastic::stress", "temperature")};
point.evaluate();
  • Declarations (core/material_base.h, core/material_interface.h): callback_handle::partials(output, reads...) and independent_of(reads...) after reads(). Reads are inputs or outputs of the material, as references, pointers (null skipped) or pointer ranges. A read the callback does not declare in reads(), an output of another callback, a partial declared twice or a read both in partials() and independent_of() throws at construction. base::partial(output, read) returns the stored tmech value (scalar x scalar: value_type; otherwise tensor<T, Dim, rank(out) + rank(read)>). base::want_partials() is true only while a callback on a requested path runs.
  • Packing (core/derivatives.h, new): one Mandel index per index pair in tmech's order with the sqrt(2) weights of solvers/unknown_layout.h (width 1, 6 or 36), each pair projected onto its symmetric part. Every row of the contraction table becomes D[out] += P * D[in]. tmech's convert_tensor_to_mandel handles rank 2 and 4 only and does not project, so the packer here covers rank 0 to 8 itself.
  • Engine (core/property_engine.h): derivative<R>(target, seed) plans once. It walks the target's upstream (upstream_of()) from the seed in topological order. Every node a derivative reaches needs a partial (or an independent_of()) for each read that carries one. Otherwise it throws and names the whole path, e.g. strain::strain -> a::y -> b::y -> c::y. A node reached over an edge that is not a declared read (material_ref, drives, a callback without declared reads) throws the same way. differentiate() packs every partial on a plan once and multiplies the blocks with Eigen. The callbacks declaring partials are wrapped at build() so want_partials() follows the request. These are new members, and update() is untouched.
  • Context and host: material_context::update() calls differentiate() after the values (a no-op without requests); material_context::derivative<R>(material, property, seed_material, seed_property). material_point::derivative<T>("material::output", "input") resolves an input to its history output (strain::strain, temperature::state), or takes material::property.
  • Finite-difference check (core/partial_check.h, new): check_partials(ctx, material) compares every declared partial with central differences of that one callback. It perturbs one Mandel component of the read at a time, so a symmetric tensor moves by (e_ij + e_ji)/2. It also perturbs every read declared independent_of() and requires every output of the callback to stay unchanged. The step is sqrt(eps) times the largest component of the read block, and the roundoff allowance is bounded by the analytic block. It restores the state, also when a perturbed rerun throws, and returns one line per wrong declaration.
  • Materials: thermal_eigenstrain (d eps_th / d alpha, d eps_th / d T) and linear_elasticity (d sigma / d eps = C, d sigma / d eps_eig = -C, d sigma / d C = Isym (x) (eps - eps_eig)). The values and the tangent output are unchanged. This gives the first DDSDDT: derivative("elastic::stress", "temperature") = -C : alpha I.

Decisions and deviations from the issue

  • linear_elasticity declares the true d sigma / d C (rank 6), not independent_of(m_C) as in the issue's linear_stress sketch. independent_of() means "zero", and a tangent that depends on the state (damage) would make that wrong without an error. It also puts the last row of the table (rank 6, 6x36 blocks) in a real material. linear_stress is not migrated; one material plus the eigenstrain is enough to prove the end-to-end path, and Migrate materials to declared partials and residual systems #219 migrates the rest.
  • derivative<R> takes the value type or a tensor of rank 2 to 8 of the context's dimension (derivatives::derivative_value); on an evaluated context the callbacks of the new plan run once more, so the result is filled at once. has_path() tells independence from a typo.
  • independent_of() covers every output of the callback. A read cannot be in both partials() and independent_of().
  • Missing partials are found when the derivative is requested, not during the update. The graph is static, so the request throws before anything runs, and unrequested paths never throw.
  • Local reads (EdgeKind::Local) carry no derivative. They are solver internals, which Property engine: residual systems set up at run time (unknowns and equations from any nodes, Newton/Picard, implicit function theorem) #218 differentiates by the implicit function theorem. Old values (input_history::old_value()) are constants of the step. Partials are with respect to the new value.
  • Non-symmetric tensors (width 9) are out of scope. A tensor's symmetry cannot be read from its type (as in unknown_layout.h), so every rank-2/4 node is minor-symmetric here. A tag for the 9-wide layout can come with the first finite-strain derivative.
  • The NUMSIM_MATERIALS_DEBUG_CHECKS finite-difference fallback for a missing partial is not implemented. It would need a second execution path in the planner and in differentiate(). check_partials() already gives per-node FD, so the fallback can follow if prototyping asks for it.
  • A request reactivates the constant callbacks once (reactivate()), so a folded callback fills its partials too.

Tests

  • tests/test_partial_derivatives.cpp (new):
    • every row of the contraction table through hand-made two-node chains (scalar/scalar/tensor2, tensor2/scalar/tensor2, tensor2/tensor2/tensor2, scalar/tensor2/tensor2, tensor2/tensor2/scalar, tensor2/tensor4/tensor2), plus scalar/scalar, parallel paths, no path (zero) and a state-dependent result
    • a missing partial on a requested path throws and names the path; one off the requested paths does not throw; independent_of() gives zero; a result type of the wrong rank throws; every declaration mistake throws
    • no derivative requested: no partials filled (counted); after a request, only the callbacks on the path fill them
    • check_partials() passes correct partials of every kind and names a wrong one, restoring the state
    • thermal eigenstrain + linear elasticity through material_point: derivative("elastic::stress", "temperature") matches central differences of the whole evaluation and equals -3 K alpha I. derivative("elastic::stress", "strain") equals the tangent output.
  • tests/test_material_conformance.cpp: every header declaring .partials( or .independent_of( must be listed, and all their partials pass check_partials() in a model built from JSON.

Gate

GCC 14, reports on, -j4: build 0 warnings, ctest 100% tests passed, 0 tests failed out of 1524.

…on (#217)

Callbacks declare partials() / independent_of() next to reads(); the engine
plans a requested derivative over the target's upstream and accumulates the
Mandel-packed partials. Hosts request derivatives via material_context and
material_point. check_partials() verifies declared partials by central
differences. thermal_eigenstrain and linear_elasticity declare their partials,
giving the first DDSDDT.
With the partials requested, a rerun at a perturbed read refilled the blocks, so a later partial of the same callback was compared at the perturbed state. Found while building #218.
@petlenz

petlenz commented Oct 28, 2026

Copy link
Copy Markdown
Member Author

Pushed a fix to check_partials (found while building #218): with the partials requested, the perturbed reruns refilled the blocks before a later partial of the same callback was packed, so that partial was compared at a perturbed state. Regression test TheFiniteDifferenceCheckComparesWithThePartialsAtTheState (y = (a:a) b) fails without the fix. Gate: 1504/1504, 0 warnings.

…, fill at once

check_partials perturbs every read declared independent_of() and requires the
outputs of the callback to stay put; its step follows the read block, its
roundoff allowance is bounded by the analytic block, and the reads and outputs
are restored when a perturbed rerun throws.

derivative<R> accepts only the scalar or a tensor of the context's dimension,
refuses calls from a callback, builds its message on the error path only, and
on an evaluated context reruns the callbacks of the new plan so the result is
filled at once. has_path() tells independence from a typo.
@petlenz

petlenz commented Oct 29, 2026

Copy link
Copy Markdown
Member Author

Review findings, addressed in f7af4ae (gate: GCC 14, -O1, 0 warnings, 1524/1524 tests).

  1. independent_of() never verified (major). check_partials now perturbs every read declared independent_of() (each Mandel component, both signs, a step of 1e-3 of the read block) and reports the output that moved by more than the tolerance (relative). The declaration carries a type-erased read probe (independent_declaration::probe) and material_interface::watch_outputs(callback) gives the deviation of the outputs. The conformance guard also greps .independent_of(. Tests: a material with a wrong independent_of fails (red: passed), a true one passes.
  2. derivative<R> accepted any R (major). Constrained with derivatives::derivative_value<R, value_type, Dim>: the scalar or a tensor of rank 2, 4, 6 or 8 of the context's value_type and Dim. Rank 6 and 8 stay legal because the result of a tensor4 by a tensor2 or tensor4 has them, so node_value alone was too narrow. The same constraint sits on material_context::derivative and material_point::derivative. Test: a concept over derivative<R> rejects tensor<double,2,4>, tensor<float,3,2>, float, int, std::string.
  3. Reads not restored on a throwing rerun. Scope guards restore the read after each perturbed run, and the callback is rerun once at the original state when the scope unwinds. Test: a callback that throws at the perturbed state leaves read and output as they were (red: read off by 6e-6).
  4. Step and allowance. The step is h = sqrt(eps) * max(|x|, block_scale), the block scale being the inf-norm of the read with a floor of 1e-3; since it is a maximum over the block, the step is the same for every component. The roundoff allowance uses min(|output|, scale * (block_scale + h)), so a large output cannot hide a small partial. Test: a correct partial of y = c tr(x)^3 at strain 1e-3 (outputs about 6e8) passes (red: the old unit floor gave a step of 6e-6 and a truncation error of 7e-5 relative). A partial wrong by 2e-5 relative at the same outputs is detected. That test is green before and after: the old check was not blind in the linear regime either, and the offset regime (a large constant plus a small partial) cannot be resolved by any step. All materials of the conformance test still pass.
  5. derivative<R> hygiene. It throws logic_error when called from a callback (the m_running guard; an argument error does not mark the context failed). The what string is built on the error paths only. On an evaluated context (an update() ran to its end) the callbacks of the new plan rerun at the held values through guarded() and the plan is accumulated at once. A plan with a stateful callback is not filled at once, because rerunning it would repeat its side effect; it follows with the next update(). Test: a derivative requested after an update is correct immediately (red: zero). NoPartialsAreComputedWithoutARequest now expects the one fill of the request.
  6. Seed not upstream returns zero silently. New property_engine::has_path(target, seed), material_context::has_path(material, property, seed_material, seed_property) and material_point::has_path(output, input); the seed counts as upstream of itself. derivative() still gives zero without a path (NoPathGivesZero).

Design remark: independent_of() applies to all outputs of a callback (documented on the declaration); a read that only some outputs depend on is declared with partials() for those and cannot use it.

@petlenz

petlenz commented Oct 29, 2026

Copy link
Copy Markdown
Member Author

Merged main (#222, #223) into this branch. No conflicts.

Gate: g++-14 -O1, 1556 tests, 0 warnings.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Property engine: partial derivatives per callback, chain rule by forward accumulation

1 participant