diff --git a/include/numsim-materials/core/derivatives.h b/include/numsim-materials/core/derivatives.h
new file mode 100644
index 0000000..1dc3b85
--- /dev/null
+++ b/include/numsim-materials/core/derivatives.h
@@ -0,0 +1,333 @@
+#ifndef NUMSIM_MATERIALS_CORE_DERIVATIVES_H
+#define NUMSIM_MATERIALS_CORE_DERIVATIVES_H
+
+#include <array>
+#include <cstddef>
+#include <functional>
+#include <numbers>
+#include <optional>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+#include <tmech/tmech.h>
+#include "numsim-materials/core/history_property.h"
+#include "numsim-materials/core/input_types.h"
+#include "numsim-materials/core/property.h"
+#include "numsim-materials/core/property_traits.h"
+
+/// Partial derivatives of a callback's outputs with respect to its reads (#217), packed into Mandel blocks.
+///
+/// Every index pair maps to one Mandel component with the weight sqrt(2) off the diagonal, the layout of
+/// solvers/unknown_layout.h. The weights make A : B = mandel(A) . mandel(B), so the chain rule over any
+/// middle quantity (scalar, tensor2, tensor4) is a plain matrix product of the packed blocks.
+namespace numsim::materials::derivatives {
+
+/// Rank of a value the graph differentiates: 0 for a scalar.
+template <typename V>
+constexpr std::size_t rank_of() {
+  if constexpr (std::is_arithmetic_v<V>)
+    return 0;
+  else
+    return V::rank();
+}
+
+/// The value behind an output or a read: T of history_property<T>, input_property<T> and input_history<T>.
+template <typename A>
+struct value_of {
+  using type = A;
+};
+template <typename T, typename Tr>
+struct value_of<numsim_core::history_property<T, Tr>> {
+  using type = T;
+};
+template <typename T, typename Tr>
+struct value_of<input_property<T, Tr>> {
+  using type = T;
+};
+template <typename T, typename Tr>
+struct value_of<input_history<T, Tr>> {
+  using type = T;
+};
+template <typename A>
+using value_of_t = typename value_of<A>::type;
+
+/// A value with a derivative: a scalar, or a symmetric tensor of rank 2 or 4 (minor symmetric).
+template <typename V, typename T, std::size_t Dim>
+concept node_value = std::is_same_v<V, T> || std::is_same_v<V, tmech::tensor<T, Dim, 2>> ||
+                     std::is_same_v<V, tmech::tensor<T, Dim, 4>>;
+
+/// The type a derivative is read as: node_value's scalar, or a tensor of rank 2 to 8 (d tensor4 / d tensor4), of the same T and Dim.
+template <typename R, typename T, std::size_t Dim>
+concept derivative_value = node_value<R, T, Dim> || std::is_same_v<R, tmech::tensor<T, Dim, 6>> ||
+                           std::is_same_v<R, tmech::tensor<T, Dim, 8>>;
+
+template <typename T, std::size_t Dim, std::size_t Rank>
+struct partial_value {
+  using type = tmech::tensor<T, Dim, Rank>;
+};
+template <typename T, std::size_t Dim>
+struct partial_value<T, Dim, 0> {
+  using type = T;
+};
+/// d out / d in: rank(out) + rank(in), a scalar for two scalars.
+template <typename Out, typename In, typename T, std::size_t Dim>
+using partial_t = typename partial_value<T, Dim, rank_of<Out>() + rank_of<In>()>::type;
+
+/// Mandel components of a symmetric index pair.
+template <std::size_t Dim>
+inline constexpr std::size_t pair_width{Dim * (Dim + 1) / 2};
+
+/// Mandel components of a value of rank `rank`, one Mandel index per index pair.
+template <std::size_t Dim>
+constexpr std::size_t width(std::size_t rank) {
+  std::size_t w{1};
+  for (std::size_t pair{0}; pair < rank / 2; ++pair)
+    w *= pair_width<Dim>;
+  return w;
+}
+
+namespace detail {
+/// The index pairs in tmech's Mandel order: 11, 22, 33, 23, 13, 12 (2D: 11, 22, 12).
+template <std::size_t Dim>
+constexpr auto pairs() {
+  if constexpr (Dim == 2)
+    return std::array<std::array<std::size_t, 2>, 3>{{{0, 0}, {1, 1}, {0, 1}}};
+  else
+    return std::array<std::array<std::size_t, 2>, 6>{{{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}}};
+}
+
+/// The Mandel component of the index pair (i, j).
+template <std::size_t Dim>
+constexpr std::size_t component(std::size_t i, std::size_t j) {
+  if (i == j)
+    return i;
+  if constexpr (Dim == 2)
+    return 2;
+  else
+    return 6 - i - j;
+}
+
+template <typename V, std::size_t... I>
+constexpr decltype(auto) entry(V& v, std::array<std::size_t, sizeof...(I)> const& idx, std::index_sequence<I...>) {
+  return v(idx[I]...);
+}
+
+constexpr std::size_t power(std::size_t base, std::size_t exponent) {
+  std::size_t p{1};
+  while (exponent-- > 0)
+    p *= base;
+  return p;
+}
+} // namespace detail
+
+/// Packs P = d out / d in (rank(out) = OutRank) into a width(OutRank) x width(rank(P) - OutRank) column-major block;
+/// each index pair is projected onto its symmetric part. Index loops: packing is the boundary to the matrices.
+template <std::size_t Dim, std::size_t OutRank, typename P, typename T>
+void pack(P const& p, T* block) {
+  if constexpr (std::is_arithmetic_v<P>) {
+    *block = p;
+  } else {
+    constexpr std::size_t rank{P::rank()}, n{rank / 2}, W{pair_width<Dim>};
+    constexpr std::size_t rows{width<Dim>(OutRank)}, cols{width<Dim>(rank - OutRank)};
+    constexpr auto pair{detail::pairs<Dim>()};
+    std::array<std::size_t, n> m{};
+    std::array<std::size_t, rank> idx{};
+    for (std::size_t flat{0}; flat < rows * cols; ++flat) {
+      for (std::size_t k{n}, rest{flat}; k-- > 0; rest /= W)
+        m[k] = rest % W; // the first pair varies slowest
+      T sum{0};
+      for (std::size_t swaps{0}; swaps < (std::size_t{1} << n); ++swaps) {
+        for (std::size_t k{0}; k < n; ++k) {
+          auto const flip{(swaps >> k) & 1};
+          idx[2 * k] = pair[m[k]][flip];
+          idx[2 * k + 1] = pair[m[k]][1 - flip];
+        }
+        sum += detail::entry(p, idx, std::make_index_sequence<rank>{});
+      }
+      T weight{T{1} / T(std::size_t{1} << n)};
+      for (auto const mk : m)
+        if (mk >= Dim)
+          weight *= std::numbers::sqrt2_v<T>;
+      block[flat / cols + (flat % cols) * rows] = weight * sum;
+    }
+  }
+}
+
+/// The inverse of pack(): a minor-symmetric tensor (or a scalar) from its column-major Mandel block.
+template <std::size_t Dim, std::size_t OutRank, typename P, typename T>
+void unpack(T const* block, P& p) {
+  if constexpr (std::is_arithmetic_v<P>) {
+    p = *block;
+  } else {
+    constexpr std::size_t rank{P::rank()}, W{pair_width<Dim>};
+    constexpr std::size_t rows{width<Dim>(OutRank)}, cols{width<Dim>(rank - OutRank)};
+    std::array<std::size_t, rank> idx{};
+    for (std::size_t flat{0}; flat < detail::power(Dim, rank); ++flat) {
+      for (std::size_t k{rank}, rest{flat}; k-- > 0; rest /= Dim)
+        idx[k] = rest % Dim;
+      std::size_t mandel{0};
+      T weight{1};
+      for (std::size_t k{0}; k < rank; k += 2) {
+        auto const c{detail::component<Dim>(idx[k], idx[k + 1])};
+        mandel = mandel * W + c;
+        if (c >= Dim)
+          weight *= std::numbers::sqrt2_v<T>;
+      }
+      detail::entry(p, idx, std::make_index_sequence<rank>{}) = block[mandel / cols + (mandel % cols) * rows] / weight;
+    }
+  }
+}
+
+/// The value of `source` as V: a property's value or a history property's new value.
+template <typename V>
+V& value_in(property_base& source) {
+  if (auto* plain{dynamic_cast<property<V>*>(&source)})
+    return plain->get();
+  if (auto* history{dynamic_cast<history_property<V, property_traits>*>(&source)})
+    return history->new_value();
+  throw std::logic_error("derivatives: '" + source.traits().id.owner + "::" + source.traits().id.name +
+                         "' holds another type than its reader");
+}
+
+/// The rank of the value `prop` holds when it has a derivative (node_value), else nothing.
+template <typename T, std::size_t Dim>
+std::optional<std::size_t> property_rank(property_base const& prop) {
+  auto const holds{[&]<typename V>(std::type_identity<V>) {
+    return dynamic_cast<property<V> const*>(&prop) || dynamic_cast<history_property<V, property_traits> const*>(&prop);
+  }};
+  if (holds(std::type_identity<T>{}))
+    return 0;
+  if (holds(std::type_identity<tmech::tensor<T, Dim, 2>>{}))
+    return 2;
+  if (holds(std::type_identity<tmech::tensor<T, Dim, 4>>{}))
+    return 4;
+  return std::nullopt;
+}
+
+/// A read the check perturbs one Mandel component at a time: its value is held by a property, given as `source`.
+template <typename T>
+class read_probe {
+public:
+  virtual ~read_probe() = default;
+  /// Mandel components of the read.
+  virtual std::size_t width() const noexcept = 0;
+  /// Mandel component k of the read's value, held by `source`.
+  virtual T read_component(property_base& source, std::size_t k) const = 0;
+  /// Adds h to Mandel component k of the read's value; the returned function restores the value.
+  virtual std::function<void()> perturb(property_base& source, std::size_t k, T h) const = 0;
+};
+
+template <typename T, std::size_t Dim, typename V>
+class typed_read : public read_probe<T> {
+public:
+  std::size_t width() const noexcept override { return derivatives::width<Dim>(rank_of<V>()); }
+
+  T read_component(property_base& source, std::size_t k) const override {
+    std::array<T, derivatives::width<Dim>(rank_of<V>())> packed{};
+    derivatives::pack<Dim, rank_of<V>()>(value_in<V>(source), packed.data());
+    return packed[k];
+  }
+
+  std::function<void()> perturb(property_base& source, std::size_t k, T h) const override {
+    auto& v{value_in<V>(source)};
+    std::function<void()> restore{[&v, old = v] { v = old; }};
+    if constexpr (std::is_arithmetic_v<V>) {
+      v += h;
+    } else {
+      std::array<T, derivatives::width<Dim>(rank_of<V>())> unit{};
+      unit[k] = h;
+      V step;
+      derivatives::unpack<Dim, rank_of<V>()>(unit.data(), step);
+      v = V{v + step};
+    }
+    return restore;
+  }
+};
+
+/// One declared partial d output / d read of a callback: the value the callback fills, packed on demand.
+template <typename T>
+class partial_block : public read_probe<T> {
+public:
+  partial_block(std::size_t rows, std::size_t cols) noexcept : m_rows(rows), m_cols(cols) {}
+  partial_block(partial_block const&) = delete;
+  partial_block& operator=(partial_block const&) = delete;
+
+  std::size_t rows() const noexcept { return m_rows; }
+  std::size_t cols() const noexcept { return m_cols; }
+  std::size_t width() const noexcept override { return m_cols; }
+  /// The partial as a rows() x cols() column-major Mandel block.
+  virtual void pack(T* block) const = 0;
+  /// The output's value, rows() Mandel components.
+  virtual void pack_output(T* value) const = 0;
+
+private:
+  std::size_t m_rows, m_cols;
+};
+
+template <typename T, std::size_t Dim, typename Out, typename In>
+class typed_partial final : public partial_block<T> {
+public:
+  using output_value = value_of_t<Out>;
+  using read_value = value_of_t<In>;
+  static constexpr std::size_t out_rank{rank_of<output_value>()};
+
+  explicit typed_partial(Out const& output)
+      : partial_block<T>{derivatives::width<Dim>(out_rank), derivatives::width<Dim>(rank_of<read_value>())}, m_output(output) {}
+
+  partial_t<output_value, read_value, T, Dim> value{};
+
+  void pack(T* block) const override { derivatives::pack<Dim, out_rank>(value, block); }
+
+  void pack_output(T* out) const override {
+    if constexpr (requires(Out const& o) { o.new_value(); })
+      derivatives::pack<Dim, out_rank>(m_output.new_value(), out);
+    else
+      derivatives::pack<Dim, out_rank>(m_output, out);
+  }
+
+  T read_component(property_base& source, std::size_t k) const override { return m_read.read_component(source, k); }
+
+  std::function<void()> perturb(property_base& source, std::size_t k, T h) const override { return m_read.perturb(source, k, h); }
+
+private:
+  Out const& m_output;
+  typed_read<T, Dim, read_value> m_read;
+};
+
+/// A requested derivative as the host reads it, unpacked from its accumulated block.
+template <typename T>
+class derivative_slot_base {
+public:
+  virtual ~derivative_slot_base() = default;
+  virtual void unpack(T const* block) = 0;
+};
+
+template <typename T, std::size_t Dim, typename R>
+class derivative_slot final : public derivative_slot_base<T> {
+public:
+  explicit derivative_slot(std::size_t out_rank) : m_out_rank(out_rank) {}
+
+  R value{};
+
+  void unpack(T const* block) override {
+    switch (m_out_rank) {
+    case 0: derivatives::unpack<Dim, 0>(block, value); break;
+    case 2:
+      if constexpr (rank_of<R>() >= 2)
+        derivatives::unpack<Dim, 2>(block, value);
+      break;
+    case 4:
+      if constexpr (rank_of<R>() >= 4)
+        derivatives::unpack<Dim, 4>(block, value);
+      break;
+    default: break;
+    }
+  }
+
+private:
+  std::size_t m_out_rank;
+};
+
+} // namespace numsim::materials::derivatives
+
+#endif // NUMSIM_MATERIALS_CORE_DERIVATIVES_H
diff --git a/include/numsim-materials/core/material_base.h b/include/numsim-materials/core/material_base.h
index 2ba1839..8b72894 100755
--- a/include/numsim-materials/core/material_base.h
+++ b/include/numsim-materials/core/material_base.h
@@ -169,8 +169,36 @@ class material_base : public material_interface<Traits> {
       return *this;
     }
 
+    /// The callback fills d output / d read for each of `reads` (as in reads()) through partial() while
+    /// want_partials() is true. After reads().
+    template <typename Output, typename... Reads>
+    callback_handle& partials(Output const& output, Reads const&... reads) {
+      (each_read(reads, [&](auto const& read) { m_owner.declare_partial(m_outputs.front(), output, read); }), ...);
+      return *this;
+    }
+
+    /// The outputs of the callback do not depend on `reads`: their partials are zero. After reads().
+    template <typename... Reads>
+    callback_handle& independent_of(Reads const&... reads) {
+      (each_read(reads, [&](auto const& read) { m_owner.declare_independent(m_outputs.front(), read); }), ...);
+      return *this;
+    }
+
   private:
     friend class material_base;
+    template <typename A, typename F>
+    static void each_read(A const& read, F const& f) {
+      if constexpr (std::is_pointer_v<A>) {
+        if (read)
+          f(*read);
+      } else if constexpr (detail::pointer_range<A>) {
+        for (auto const* p : read)
+          if (p)
+            f(*p);
+      } else {
+        f(read);
+      }
+    }
     template <typename A>
     static void sort_read(std::vector<input_wire_base const*>& inputs, std::vector<void const*>& outputs, A const& read) {
       if constexpr (std::derived_from<detail::output_t<A>, input_wire_base>)
diff --git a/include/numsim-materials/core/material_context.h b/include/numsim-materials/core/material_context.h
index 78b97f3..cbd0311 100644
--- a/include/numsim-materials/core/material_context.h
+++ b/include/numsim-materials/core/material_context.h
@@ -82,9 +82,39 @@ class material_context {
 
   // --- Execution (delegated to property_engine) ---
 
+  /// Runs the callbacks, then accumulates the requested derivatives (derivative()).
   void update() {
     check_finalized("update");
     m_engine.update();
+    m_engine.differentiate();
+  }
+
+  /// d material::property / d seed_material::seed_property as R (rank of both added, a scalar for two scalars),
+  /// filled by every update() after the values. A requested path through a callback without its partial throws.
+  /// R is the value type or a tensor of rank 2 to 8 of the context's dimension.
+  template <typename R>
+    requires derivatives::derivative_value<R, typename Traits::value_type, Traits::Dim>
+  const R& derivative(const std::string& material, const std::string& property, const std::string& seed_material,
+                      const std::string& seed_property) {
+    check_finalized("derivative");
+    auto const* target{find_property(material, property)};
+    auto const* seed{find_property(seed_material, seed_property)};
+    if (!target || !seed)
+      throw std::runtime_error("derivative(): property '" + (target ? seed_material + "::" + seed_property : material + "::" + property) +
+                               "' not found");
+    return m_engine.template derivative<R>(target, seed);
+  }
+
+  /// Whether the seed is upstream of the property: otherwise derivative() is zero. Throws on an unknown name.
+  bool has_path(const std::string& material, const std::string& property, const std::string& seed_material,
+                const std::string& seed_property) {
+    check_finalized("has_path");
+    auto const* target{find_property(material, property)};
+    auto const* seed{find_property(seed_material, seed_property)};
+    if (!target || !seed)
+      throw std::runtime_error("has_path(): property '" + (target ? seed_material + "::" + seed_property : material + "::" + property) +
+                               "' not found");
+    return m_engine.has_path(target, seed);
   }
 
   void update_property(const std::string& material, const std::string& property,
diff --git a/include/numsim-materials/core/material_interface.h b/include/numsim-materials/core/material_interface.h
index d6aed65..750dd1c 100644
--- a/include/numsim-materials/core/material_interface.h
+++ b/include/numsim-materials/core/material_interface.h
@@ -20,6 +20,7 @@
 #include <vector>
 #include <tmech/tmech.h>
 
+#include "numsim-materials/core/derivatives.h"
 #include "numsim-materials/core/property_traits.h"
 #include "numsim-materials/core/property.h"
 #include "numsim-materials/core/history_property.h"
@@ -61,6 +62,22 @@ bool same_value(T const& a, T const& b) {
   }
 }
 
+/// How far `b` is from `a` relative to the larger of the two: 0 when equal, infinite when they differ and one is not finite.
+template <verifiable T>
+double relative_deviation(T const& a, T const& b) {
+  double distance, size;
+  if constexpr (std::is_arithmetic_v<T>) {
+    distance = std::abs(static_cast<double>(a) - static_cast<double>(b));
+    size = std::max(std::abs(static_cast<double>(a)), std::abs(static_cast<double>(b)));
+  } else if constexpr (requires { tmech::norm(a - b); }) {
+    distance = static_cast<double>(tmech::norm(a - b));
+    size = static_cast<double>(std::max(tmech::norm(a), tmech::norm(b)));
+  } else {
+    return a == b ? 0 : 1;
+  }
+  return distance == 0 ? 0 : std::isfinite(distance) && std::isfinite(size) ? distance / size : std::numeric_limits<double>::infinity();
+}
+
 /// Whether a value is finite; types without a norm are.
 template <verifiable T>
 bool is_finite(T const& a) {
@@ -209,6 +226,7 @@ class material_interface {
   using property_base_type = numsim_core::property_base<property_traits>;
   /// Copies a value; the returned check compares it with the current one. Empty for types never constant.
   using snapshot = std::function<std::function<detail::verdict()>()>;
+  using deviation_probe = std::function<std::function<double()>()>;
 
   /// An output declared a function of the parameters, of `reads` (all inputs when unset) and of the outputs
   /// `own_reads` of this material (material_base::constant_output, callback_handle::reads).
@@ -226,6 +244,7 @@ class material_interface {
     std::vector<input_wire_base const*> inputs;
     std::vector<property_base_type const*> outputs;
     std::vector<input_wire_base const*> drives; ///< local inputs whose source callbacks it reruns
+    std::vector<property_base_type const*> own_reads; ///< outputs of this material it reads
   };
   const std::vector<read_declaration>& declared_reads() const noexcept { return m_declared_reads; }
 
@@ -267,7 +286,126 @@ class material_interface {
   /// The engine's revision counter, bumped by set_parameter(); bound at finalize().
   void bind_revision(std::size_t* revision) noexcept { m_revision = revision; }
 
+  /// d output / d read declared by a callback (callback_handle::partials()); the read is an input or an output of
+  /// this material, the block holds what partial() returns.
+  struct partial_declaration {
+    property_base_type const* callback;
+    property_base_type const* output;
+    input_wire_base const* input;  ///< the read when it is an input ...
+    property_base_type const* own; ///< ... else the output of this material it is
+    void const* output_value;      ///< the addresses partial() is called with
+    void const* read_value;
+    std::unique_ptr<derivatives::partial_block<value_type>> block;
+  };
+  const std::vector<partial_declaration>& partials() const noexcept { return m_partials; }
+
+  /// A read whose partials are declared zero (callback_handle::independent_of()), as in partial_declaration;
+  /// `probe` perturbs the read for check_partials(), empty for a read that is not a scalar or tensor.
+  struct independent_declaration {
+    property_base_type const* callback;
+    input_wire_base const* input;
+    property_base_type const* own;
+    std::unique_ptr<derivatives::read_probe<value_type>> probe;
+  };
+  const std::vector<independent_declaration>& independent_reads() const noexcept { return m_independent; }
+
+  /// The worst relative deviation of an output of a callback from a copy taken earlier.
+  struct output_deviation {
+    property_base_type const* output{nullptr};
+    double relative{0};
+  };
+  /// Copies the outputs of the callback filling `callback` (with declared reads); the returned function reports the
+  /// largest relative deviation of any from its copy. Outputs of types without a comparison are not watched.
+  std::function<output_deviation()> watch_outputs(property_base_type const* callback) const {
+    auto const it{std::ranges::find(m_declared_reads, callback, &read_declaration::callback)};
+    if (it == m_declared_reads.end())
+      throw std::logic_error("watch_outputs(): the callback has no declared reads");
+    std::vector<std::pair<property_base_type const*, std::function<double()>>> watched;
+    for (auto const* output : it->outputs)
+      for (auto const& [address, entry] : m_outputs)
+        if (entry.property == output && entry.deviation)
+          watched.emplace_back(output, entry.deviation());
+    return [watched = std::move(watched)] {
+      output_deviation worst;
+      for (auto const& [output, deviation] : watched)
+        if (auto const relative{deviation()}; relative > worst.relative || !worst.output)
+          worst = {output, relative};
+      return worst;
+    };
+  }
+
+  /// Whether `callback` fills its partials when it runs: the engine sets it for the callbacks on a requested derivative.
+  void request_partials(property_base_type const* callback, bool on) {
+    std::erase(m_partials_requested, callback);
+    if (on)
+      m_partials_requested.push_back(callback);
+  }
+  bool partials_requested(property_base_type const* callback) const noexcept {
+    return std::ranges::contains(m_partials_requested, callback);
+  }
+
+  /// Runs `run`, the update of `callback`, with want_partials() as requested for it; the engine wraps those callbacks.
+  template <typename Run>
+  void run_callback(property_base_type const* callback, Run const& run) {
+    struct restore {
+      bool& flag;
+      bool outer;
+      ~restore() { flag = outer; }
+    } const guard{m_want_partials, std::exchange(m_want_partials, partials_requested(callback))};
+    run();
+  }
+
 protected:
+  /// True while the running callback is to fill its declared partials: a requested derivative passes through it.
+  bool want_partials() const noexcept { return m_want_partials; }
+
+  /// The partial d output / d read the callback fills, declared by partials(): rank(output) + rank(read).
+  template <typename Out, typename In>
+  auto& partial(Out const& output, In const& read) {
+    using block_type = derivatives::typed_partial<value_type, Traits::Dim, Out, In>;
+    auto const it{std::ranges::find_if(m_partials, [&](auto const& p) {
+      return p.output_value == static_cast<void const*>(&output) && p.read_value == static_cast<void const*>(&read);
+    })};
+    auto* block{it != m_partials.end() ? dynamic_cast<block_type*>(it->block.get()) : nullptr};
+    if (!block)
+      throw std::logic_error("partial(): material '" + m_name + "' fills a partial its callback did not declare in partials()");
+    return block->value;
+  }
+
+  /// The callback filling `callback` (after its reads()) fills d output / d read through partial().
+  template <typename Out, typename In>
+  void declare_partial(void const* callback, Out const& output, In const& read) {
+    static_assert(derivatives::node_value<derivatives::value_of_t<Out>, value_type, Traits::Dim> &&
+                      derivatives::node_value<derivatives::value_of_t<In>, value_type, Traits::Dim>,
+                  "partials(): outputs and reads are scalars or symmetric tensors of rank 2 or 4");
+    auto const& [d, entry]{declared_read_of("partials", callback)};
+    auto const& out{output_entry_of("partials", &output)};
+    if (!std::ranges::contains(d->outputs, out.property))
+      throw std::logic_error("partials(): '" + key_of(out) + "' is not an output of the callback of '" + key_of(entry) + "'");
+    auto const [input, own]{read_of<In>("partials", *d, entry, read)};
+    if (std::ranges::any_of(m_partials, [&](auto const& p) { return p.output == out.property && p.input == input && p.own == own; }))
+      throw std::logic_error("partials(): '" + key_of(out) + "' declares a partial twice");
+    if (is_independent(d->callback, input, own))
+      throw std::logic_error("partials(): '" + key_of(entry) + "' declares a read both in partials() and independent_of()");
+    m_partials.push_back({d->callback, out.property, input, own, &output, &read,
+                          std::make_unique<derivatives::typed_partial<value_type, Traits::Dim, Out, In>>(output)});
+  }
+
+  /// The outputs of the callback filling `callback` do not depend on `read`: its partials are zero. It holds for every
+  /// output of the callback, and check_partials() perturbs the read to verify it.
+  template <typename In>
+  void declare_independent(void const* callback, In const& read) {
+    auto const& [d, entry]{declared_read_of("independent_of", callback)};
+    auto const [input, own]{read_of<In>("independent_of", *d, entry, read)};
+    if (is_independent(d->callback, input, own) ||
+        std::ranges::any_of(m_partials, [&](auto const& p) { return p.callback == d->callback && p.input == input && p.own == own; }))
+      throw std::logic_error("independent_of(): '" + key_of(entry) + "' declares a read twice or also in partials()");
+    std::unique_ptr<derivatives::read_probe<value_type>> probe;
+    if constexpr (derivatives::node_value<derivatives::value_of_t<In>, value_type, Traits::Dim>)
+      probe = std::make_unique<derivatives::typed_read<value_type, Traits::Dim, derivatives::value_of_t<In>>>();
+    m_independent.push_back({d->callback, input, own, std::move(probe)});
+  }
+
   template <typename T>
   T& add_property(std::string const& property_name, T&& data = T()) {
     T& value{m_property_handler.template add_property<T>(name(), property_name, std::forward<T>(data))};
@@ -279,14 +417,22 @@ class material_interface {
           return !detail::is_finite(value) ? not_finite : detail::same_value(copy, value) ? same : changed;
         }};
       };
-    m_outputs[&value] = {*m_property_handler.find(name(), property_name), std::move(take)};
+    deviation_probe deviation;
+    if constexpr (detail::verifiable<T>)
+      deviation = [&value] { return std::function<double()>{[&value, copy = value] { return detail::relative_deviation(copy, value); }}; };
+    m_outputs[&value] = {*m_property_handler.find(name(), property_name), std::move(take), std::move(deviation)};
     return value;
   }
 
   template <typename T>
   history_property<T>& add_history(std::string const& property_name, T&& data = T()) {
     auto& value{m_property_handler.template add_history<T>(name(), property_name, std::forward<T>(data))};
-    m_outputs[&value] = {*m_property_handler.find(name(), property_name), {}}; // never constant
+    deviation_probe deviation;
+    if constexpr (detail::verifiable<T>)
+      deviation = [&value] {
+        return std::function<double()>{[&value, copy = value.new_value()] { return detail::relative_deviation(copy, value.new_value()); }};
+      };
+    m_outputs[&value] = {*m_property_handler.find(name(), property_name), {}, std::move(deviation)}; // never constant
     return value;
   }
 
@@ -350,7 +496,8 @@ class material_interface {
     auto const entries{outputs | std::views::transform([this](void const* output) { return &output_entry_of("reads", output); }) |
                        std::ranges::to<std::vector>()};
     m_declared_reads.push_back({callback.property, inputs,
-                                entries | std::views::transform(&output_entry::property) | std::ranges::to<std::vector>(), {}});
+                                entries | std::views::transform(&output_entry::property) | std::ranges::to<std::vector>(), {},
+                                own_reads});
     if (verifiable && std::ranges::none_of(entries, [this](auto const* entry) {
           return entry->property->is_history() || std::ranges::contains(m_stateful, entry->property);
         }))
@@ -432,6 +579,7 @@ class material_interface {
   struct output_entry {
     property_base_type const* property;
     snapshot take;
+    deviation_probe deviation; ///< copies the value; the returned function gives its relative deviation from the copy
   };
   std::unordered_map<void const*, output_entry> m_outputs; ///< value address -> property
   std::vector<constant_declaration> m_constant_outputs;
@@ -466,6 +614,35 @@ class material_interface {
   bool has_writer(output_entry const& entry) const {
     return std::ranges::contains(m_written_by, entry.property, &decltype(m_written_by)::value_type::first);
   }
+  /// The read declaration of the callback filling `callback`, which partials() and independent_of() follow.
+  std::pair<read_declaration*, output_entry const&> declared_read_of(char const* what, void const* callback) {
+    auto const& entry{output_entry_of(what, callback)};
+    auto const it{std::ranges::find(m_declared_reads, entry.property, &read_declaration::callback)};
+    if (it == m_declared_reads.end())
+      throw std::logic_error(std::string{what} + "(): '" + key_of(entry) + "' declares partials before its reads");
+    return {&*it, entry};
+  }
+  /// A read of `d` as (input, own output); throws when the callback does not read it.
+  template <typename In>
+  std::pair<input_wire_base const*, property_base_type const*> read_of(char const* what, read_declaration const& d,
+                                                                      output_entry const& entry, In const& read) const {
+    if constexpr (std::derived_from<In, input_wire_base>) {
+      if (std::ranges::contains(d.inputs, static_cast<input_wire_base const*>(&read)))
+        return {&read, nullptr};
+    } else if (auto const it{m_outputs.find(&read)};
+               it != m_outputs.end() && std::ranges::contains(d.own_reads, it->second.property)) {
+      return {nullptr, it->second.property};
+    }
+    throw std::logic_error(std::string{what} + "(): '" + key_of(entry) + "' names a read its callback does not declare in reads()");
+  }
+  bool is_independent(property_base_type const* callback, input_wire_base const* input, property_base_type const* own) const {
+    return std::ranges::any_of(m_independent, [&](auto const& i) { return i.callback == callback && i.input == input && i.own == own; });
+  }
+
+  std::vector<partial_declaration> m_partials;
+  std::vector<independent_declaration> m_independent;
+  std::vector<property_base_type const*> m_partials_requested; ///< callbacks filling their partials
+  bool m_want_partials{false};
   std::size_t* m_revision{nullptr};
   std::vector<property_base_type const*> m_reads_all; ///< callbacks reading every input (declare_reads_all)
 };
diff --git a/include/numsim-materials/core/partial_check.h b/include/numsim-materials/core/partial_check.h
new file mode 100644
index 0000000..d90678a
--- /dev/null
+++ b/include/numsim-materials/core/partial_check.h
@@ -0,0 +1,159 @@
+#ifndef NUMSIM_MATERIALS_CORE_PARTIAL_CHECK_H
+#define NUMSIM_MATERIALS_CORE_PARTIAL_CHECK_H
+
+#include <algorithm>
+#include <cmath>
+#include <format>
+#include <limits>
+#include <ranges>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+#include "numsim-materials/core/material_context.h"
+
+namespace numsim::materials {
+
+namespace partial_check_detail {
+/// Runs `f` when the scope ends, also while an exception unwinds it; what `f` throws is dropped.
+template <typename F>
+struct on_exit {
+  F f;
+  ~on_exit() {
+    try {
+      f();
+    } catch (...) {
+    }
+  }
+};
+template <typename F>
+on_exit(F) -> on_exit<F>;
+} // namespace partial_check_detail
+
+/// Compares every partial `material` declares with central differences of its own callback, one Mandel component
+/// of the read at a time, at the state the context holds (after an update()), and perturbs every read declared
+/// independent_of() to see that no output of the callback moves. Returns one line per wrong declaration; the
+/// outputs and reads are restored, also when a perturbed callback throws.
+template <typename Traits>
+std::vector<std::string> check_partials(material_context<Traits>& ctx, std::string const& material,
+                                        typename Traits::value_type tolerance = 1e-6) {
+  using T = typename Traits::value_type;
+  using partial_check_detail::on_exit;
+  constexpr T eps{std::numeric_limits<T>::epsilon()};
+  constexpr T scale_floor{1e-3}; // the scale of a read that is zero, or nearly
+  constexpr T independence_step{1e-3};
+  auto* mat{ctx.find(material)};
+  if (!mat)
+    throw std::invalid_argument("check_partials(): no material '" + material + "'");
+
+  auto const run_of{[&](auto const* callback) {
+    auto* runs{ctx.find_property(material, callback->traits().id.name)};
+    return [runs] { runs->traits().update(); };
+  }};
+  auto const source_of{[&](auto const& d) {
+    return d.input ? ctx.find_property(d.input->source_owner(), d.input->source_name())
+                   : ctx.find_property(material, d.own->traits().id.name);
+  }};
+  auto const key_of{[](auto const* source) { return source->traits().id.owner + "::" + source->traits().id.name; }};
+  // the largest component of the read: one step for all of them
+  auto const read_scale{[&](auto const& read, auto& source) {
+    T scale{scale_floor};
+    for (std::size_t k{0}; k < read.width(); ++k)
+      scale = std::max(scale, std::abs(read.read_component(source, k)));
+    return scale;
+  }};
+
+  std::vector<std::string> mismatches;
+  auto const& partials{mat->partials()};
+  std::vector<property_base const*> callbacks;
+  for (auto const& p : partials)
+    if (!std::ranges::contains(callbacks, p.callback))
+      callbacks.push_back(p.callback);
+
+  for (auto const* callback : callbacks) {
+    auto const run{run_of(callback)};
+    on_exit rerun{run}; // the outputs at the state the check found
+    bool const requested{mat->partials_requested(callback)};
+    mat->request_partials(callback, true);
+    {
+      on_exit reset{[&] { mat->request_partials(callback, requested); }};
+      run();
+    }
+    auto of_callback{partials | std::views::filter([callback](auto const& q) { return q.callback == callback; })};
+    // packed before any perturbed rerun overwrites them
+    std::vector<std::vector<T>> packed;
+    for (auto const& p : of_callback)
+      p.block->pack(packed.emplace_back(p.block->rows() * p.block->cols()).data());
+    for (auto const& [p, analytic] : std::views::zip(of_callback, packed)) {
+      auto const& block{*p.block};
+      std::vector<T> plus(block.rows()), minus(block.rows());
+      auto* source{source_of(p)};
+      auto const read{key_of(source)};
+      T scale{0};
+      for (auto const a : analytic)
+        scale = std::max(scale, std::abs(a));
+      T const read_size{read_scale(block, *source)};
+      T const h{std::sqrt(eps) * read_size};
+      for (std::size_t k{0}; k < block.cols(); ++k) {
+        for (auto [step, values] : {std::pair{h, &plus}, std::pair{-h, &minus}}) {
+          on_exit restore{block.perturb(*source, k, step)};
+          run();
+          block.pack_output(values->data());
+        }
+        T error{0}, size{scale}, output{0};
+        for (std::size_t r{0}; r < block.rows(); ++r) {
+          T const fd{(plus[r] - minus[r]) / (2 * h)};
+          error = std::max(error, std::abs(fd - analytic[r + k * block.rows()]));
+          size = std::max(size, std::abs(fd));
+          output = std::max({output, std::abs(plus[r]), std::abs(minus[r])});
+        }
+        // Rounding of the differences, for an output no larger than the analytic block can produce: a large
+        // output must not hide a small partial.
+        T const allowance{10 * eps * std::min(output, scale * (read_size + h)) / h};
+        if (error > tolerance * size + allowance) {
+          mismatches.push_back(std::format("'{}::{}' wrt '{}': Mandel column {} differs from central differences by {:g} "
+                                           "(largest entry {:g})",
+                                           material, p.output->traits().id.name, read, k, error, size));
+          break;
+        }
+      }
+    }
+  }
+
+  std::vector<property_base const*> independents;
+  for (auto const& i : mat->independent_reads())
+    if (i.probe && !std::ranges::contains(independents, i.callback))
+      independents.push_back(i.callback);
+  for (auto const* callback : independents) {
+    auto const run{run_of(callback)};
+    on_exit rerun{run};
+    for (auto const& i : mat->independent_reads()) {
+      if (i.callback != callback || !i.probe)
+        continue;
+      auto* source{source_of(i)};
+      auto const read{key_of(source)};
+      T const h{independence_step * read_scale(*i.probe, *source)};
+      run(); // the outputs at the state the check found
+      auto const deviation{mat->watch_outputs(callback)};
+      bool moved{false};
+      for (std::size_t k{0}; k < i.probe->width() && !moved; ++k)
+        for (T const step : {h, -h}) {
+          on_exit restore{i.probe->perturb(*source, k, step)};
+          run();
+          auto const worst{deviation()};
+          if (worst.relative > tolerance) {
+            mismatches.push_back(std::format("'{}::{}' is declared independent_of '{}', but Mandel component {} of the read "
+                                             "changes it by {:g} (relative)",
+                                             material, worst.output->traits().id.name, read, k, worst.relative));
+            moved = true;
+            break;
+          }
+        }
+    }
+  }
+  return mismatches;
+}
+
+} // namespace numsim::materials
+
+#endif // NUMSIM_MATERIALS_CORE_PARTIAL_CHECK_H
diff --git a/include/numsim-materials/core/property_engine.h b/include/numsim-materials/core/property_engine.h
index 5d38990..1a9565d 100644
--- a/include/numsim-materials/core/property_engine.h
+++ b/include/numsim-materials/core/property_engine.h
@@ -3,6 +3,7 @@
 
 #include <algorithm>
 #include <cstddef>
+#include <format>
 #include <functional>
 #include <iterator>
 #include <memory>
@@ -19,6 +20,8 @@
 #include <unordered_set>
 #include <utility>
 #include <vector>
+#include <Eigen/Core>
+#include "numsim-materials/core/derivatives.h"
 #include "numsim-materials/core/property_traits.h"
 #include "numsim-materials/core/property.h"
 #include "numsim-materials/core/history_property.h"
@@ -88,6 +91,7 @@ class property_engine {
     if constexpr (NUMSIM_MATERIALS_DEBUG_CHECKS)
       track_reads(interfaces);
     route_to_writers(interfaces);
+    index_partials(interfaces);
   }
 
   /// Run the callbacks still to run, in topological order: a constant one only until it has run.
@@ -99,6 +103,7 @@ class property_engine {
       for (auto* prop : m_active)
         prop->traits().update();
       retire([](property_base*) { return true; });
+      m_evaluated = true;
     });
   }
 
@@ -234,6 +239,60 @@ class property_engine {
     return g;
   }
 
+  /// d target / d seed as R (rank(target) + rank(seed), a scalar for two scalars), filled by differentiate().
+  /// Resolved once: a path from the seed to the target through a callback without the partial it needs throws.
+  /// On an evaluated context the callbacks on the path run once more and fill the result at once.
+  template <typename R>
+    requires derivatives::derivative_value<R, typename Traits::value_type, Traits::Dim>
+  const R& derivative(const property_base* target, const property_base* seed) {
+    using slot_type = derivatives::derivative_slot<value_type, Traits::Dim, R>;
+    auto const what{[&] { return "derivative(): d '" + prop_key(target->traits().id) + "' / d '" + prop_key(seed->traits().id) + "'"; }};
+    if (m_running)
+      throw std::logic_error(what() + ": called from a callback");
+    if (auto const it{std::ranges::find_if(m_plans, [&](auto const& p) { return p.target == target && p.seed == seed; })};
+        it != m_plans.end()) {
+      if (auto* slot{dynamic_cast<slot_type*>(it->slot.get())})
+        return slot->value;
+      throw std::logic_error(what() + " is requested as another type already");
+    }
+    auto const target_rank{derivatives::property_rank<value_type, Traits::Dim>(*target)};
+    auto const seed_rank{derivatives::property_rank<value_type, Traits::Dim>(*seed)};
+    if (!target_rank || !seed_rank)
+      throw std::logic_error(what() + ": both must be scalars or symmetric tensors of rank 2 or 4");
+    if (*target_rank + *seed_rank != derivatives::rank_of<R>())
+      throw std::logic_error(std::format("{} has rank {}, the requested type rank {}", what(), *target_rank + *seed_rank,
+                                         derivatives::rank_of<R>()));
+    auto plan{plan_derivative(target, seed, *seed_rank)};
+    auto slot{std::make_unique<slot_type>(*target_rank)};
+    std::vector<value_type> const zero(derivatives::width<Traits::Dim>(*target_rank) * derivatives::width<Traits::Dim>(*seed_rank));
+    slot->unpack(zero.data());
+    auto& value{slot->value};
+    plan.slot = std::move(slot);
+    for (auto const& step : plan.steps)
+      step.material->request_partials(step.callback, true);
+    m_plans.push_back(std::move(plan));
+    reactivate(); // a constant callback fills its partials once more
+    if (m_evaluated)
+      fill_now(m_plans.back());
+    return value;
+  }
+
+  /// Whether `seed` is upstream of `target` (or the target itself): a derivative of one by the other can be non-zero.
+  /// Without a path derivative() stays zero.
+  bool has_path(const property_base* target, const property_base* seed) const {
+    return upstream_of(target).contains(const_cast<property_base*>(seed));
+  }
+
+  /// Accumulates the requested derivatives from the partials the last update filled: D[out] += P * D[in].
+  void differentiate() {
+    if (m_plans.empty())
+      return;
+    for (auto& [block, packed] : m_blocks)
+      block->pack(packed.data());
+    for (auto& plan : m_plans)
+      accumulate(plan);
+  }
+
   void dump() const {
     std::println(stderr, "=== Property execution order ({} properties) ===",
                  m_property_execution_order.size());
@@ -577,6 +636,7 @@ class property_engine {
     } catch (...) {
       m_running = false;
       m_failed = true;
+      m_evaluated = false;
       throw;
     }
     m_running = false;
@@ -702,6 +762,177 @@ class property_engine {
     return result;
   }
 
+  // --- derivatives (#217) ---
+  using value_type = typename Traits::value_type;
+  using matrix = Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic>;
+  using read_declaration = typename material_interface_type::read_declaration;
+
+  struct derivative_step {
+    std::size_t out, in; ///< nodes of the plan
+    std::size_t block;   ///< in m_blocks
+    material_interface_type* material;
+    const property_base* callback;
+  };
+  /// One requested derivative: the nodes between seed and target that carry one, in topological order.
+  struct derivative_plan {
+    const property_base* target;
+    const property_base* seed;
+    std::vector<matrix> D;              ///< d node / d seed; node 0 is the seed
+    std::vector<derivative_step> steps; ///< in topological order of their outputs
+    std::optional<std::size_t> target_node;
+    std::unique_ptr<derivatives::derivative_slot_base<value_type>> slot;
+  };
+
+  /// D[out] += P * D[in] over the steps of one plan, from the packed blocks.
+  void accumulate(derivative_plan& plan) {
+    plan.D.front().setIdentity();
+    for (auto& d : plan.D | std::views::drop(1))
+      d.setZero();
+    for (auto const& step : plan.steps)
+      plan.D[step.out].noalias() += m_blocks[step.block].second * plan.D[step.in];
+    if (plan.target_node)
+      plan.slot->unpack(plan.D[*plan.target_node].data());
+  }
+
+  /// Reruns the callbacks of a new plan at the values held and accumulates it. A stateful callback is not rerun: its
+  /// partials come with the next update().
+  void fill_now(derivative_plan& plan) {
+    auto const on_plan{[&](auto const* prop) {
+      return std::ranges::any_of(plan.steps, [&](auto const& step) { return step.callback == prop; });
+    }};
+    if (std::ranges::any_of(m_callbacks, [&](auto const* prop) { return on_plan(prop) && m_stateful.contains(prop); }))
+      return;
+    guarded("derivative", [&] {
+      for (auto* prop : m_callbacks | std::views::filter(on_plan))
+        prop->traits().update();
+    });
+    for (auto& [block, packed] : m_blocks)
+      block->pack(packed.data());
+    accumulate(plan);
+  }
+
+  /// Indexes the declared reads by output and runs every callback declaring partials with want_partials() as requested.
+  void index_partials(const std::vector<material_interface_type*>& interfaces) {
+    m_plans.clear();
+    m_blocks.clear();
+    m_read_declarations.clear();
+    for (auto* mat : interfaces) {
+      for (auto const& d : mat->declared_reads())
+        for (auto const* out : d.outputs)
+          m_read_declarations.emplace(out, std::pair{mat, &d});
+      for (auto const& p : mat->partials())
+        if (auto* prop{find_in_graph(mat->name(), p.callback->traits().id.name)}; m_partials_wrapped.insert(prop).second) {
+          auto& update{prop->traits().update};
+          update = [inner = std::move(update), mat, callback = p.callback] { mat->run_callback(callback, inner); };
+        }
+    }
+  }
+
+  /// The steps of d target / d seed over the target's upstream: every node a derivative reaches needs the partial
+  /// of each read that carries one, declared or declared zero (independent_of()).
+  derivative_plan plan_derivative(const property_base* target, const property_base* seed, std::size_t seed_rank) {
+    auto const seed_width{static_cast<Eigen::Index>(derivatives::width<Traits::Dim>(seed_rank))};
+    derivative_plan plan{target, seed, {}, {}, std::nullopt, nullptr};
+    plan.D.emplace_back(seed_width, seed_width);
+    std::unordered_map<const property_base*, std::size_t> node{{seed, 0}};
+    std::unordered_map<const property_base*, const property_base*> via; ///< the read each node is first reached over
+    auto const& up{upstream_of(target)};
+    if (!up.contains(const_cast<property_base*>(seed)))
+      return plan;
+    for (auto const* prop : m_property_execution_order | std::views::filter([&](auto* p) { return p != seed && up.contains(p); })) {
+      std::vector<const property_base*> reached;
+      if (auto const edges{m_reverse_edges.find(const_cast<property_base*>(prop))}; edges != m_reverse_edges.end())
+        std::ranges::copy_if(edges->second, std::back_inserter(reached), [&](auto const* src) { return node.contains(src); });
+      if (reached.empty())
+        continue;
+      auto const declared{m_read_declarations.find(prop)};
+      if (declared == m_read_declarations.end())
+        throw missing_partial(target, seed, reached.front(), prop, via, "its callback declares no reads");
+      auto* mat{declared->second.first};
+      auto const& d{*declared->second.second};
+      std::vector<derivative_step> steps;
+      std::vector<const property_base*> sources;
+      Eigen::Index width{0};
+      auto const take{[&](input_wire_base const* input, const property_base* own, const property_base* src) {
+        sources.push_back(src);
+        if (!node.contains(src))
+          return;
+        auto const same_read{[&](auto const& p) { return p.input == input && p.own == own; }};
+        auto const& partials{mat->partials()};
+        if (auto const p{std::ranges::find_if(partials, [&](auto const& q) { return q.output == prop && same_read(q); })};
+            p != partials.end()) {
+          steps.push_back({0, node.at(src), block_of(*p->block), mat, p->callback});
+          width = static_cast<Eigen::Index>(p->block->rows());
+          via.try_emplace(prop, src);
+        } else if (!std::ranges::any_of(mat->independent_reads(), [&](auto const& i) { return i.callback == d.callback && same_read(i); })) {
+          throw missing_partial(target, seed, src, prop, via, "its callback declares neither partials() nor independent_of() for that read");
+        }
+      }};
+      for (auto const* input : d.inputs)
+        if (input->edge_kind() == EdgeKind::Global)
+          if (auto const* src{find_in_graph(input->source_owner(), input->source_name())})
+            take(input, nullptr, src);
+      for (auto const* own : d.own_reads)
+        take(nullptr, own, own);
+      for (auto const* src : reached)
+        if (!std::ranges::contains(sources, src) && !std::ranges::contains(d.outputs, src))
+          throw missing_partial(target, seed, src, prop, via, "it is reached over a '" + std::string{to_string(kind_of(src, prop))} +
+                                                              "' edge, not a declared read, which carries no partials");
+      if (steps.empty())
+        continue;
+      auto const index{plan.D.size()};
+      plan.D.emplace_back(width, seed_width);
+      for (auto const& step : steps)
+        plan.steps.push_back({index, step.in, step.block, step.material, step.callback});
+      node.emplace(prop, index);
+    }
+    if (auto const it{node.find(target)}; it != node.end())
+      plan.target_node = it->second;
+    return plan;
+  }
+
+  /// The packed block of a partial, shared by every plan through it.
+  std::size_t block_of(const derivatives::partial_block<value_type>& block) {
+    if (auto const it{std::ranges::find(m_blocks, &block, &std::pair<const derivatives::partial_block<value_type>*, matrix>::first)};
+        it != m_blocks.end())
+      return static_cast<std::size_t>(it - m_blocks.begin());
+    m_blocks.emplace_back(&block, matrix(static_cast<Eigen::Index>(block.rows()), static_cast<Eigen::Index>(block.cols())));
+    return m_blocks.size() - 1;
+  }
+
+  edge_kind kind_of(const property_base* from, const property_base* to) const {
+    auto const it{std::ranges::find_if(m_edges, [&](auto const& e) { return e.from == from && e.to == to; })};
+    return it != m_edges.end() ? it->kind : edge_kind::dependency;
+  }
+
+  /// Names the path seed -> ... -> from -> at -> ... -> target in the error.
+  std::logic_error missing_partial(const property_base* target, const property_base* seed, const property_base* from,
+                                   const property_base* at, const std::unordered_map<const property_base*, const property_base*>& via,
+                                   const std::string& why) const {
+    std::vector<const property_base*> path{from, at};
+    for (auto it{via.find(from)}; it != via.end() && path.front() != seed; it = via.find(it->second))
+      path.insert(path.begin(), it->second);
+    auto const& up{upstream_of(target)};
+    for (auto const* p{at}; p != target;) { // any edge into the target's upstream leads on to the target
+      auto const next{std::ranges::find_if(m_edges, [&](auto const& e) {
+        return e.from == p && e.kind != edge_kind::local && up.contains(const_cast<property_base*>(e.to));
+      })};
+      if (next == m_edges.end())
+        break;
+      path.push_back(p = next->to);
+    }
+    std::string text;
+    for (auto const* p : path)
+      text += (text.empty() ? "" : " -> ") + prop_key(p->traits().id);
+    return std::logic_error("derivative(): d '" + prop_key(target->traits().id) + "' / d '" + prop_key(seed->traits().id) + "' passes '" +
+                            prop_key(at->traits().id) + "' from '" + prop_key(from->traits().id) + "', but " + why + "; path " + text);
+  }
+
+  std::vector<derivative_plan> m_plans;
+  std::vector<std::pair<const derivatives::partial_block<value_type>*, matrix>> m_blocks; ///< packed once per differentiate()
+  std::unordered_map<const property_base*, std::pair<material_interface_type*, const read_declaration*>> m_read_declarations;
+  std::unordered_set<const property_base*> m_partials_wrapped; ///< callbacks wrapped by index_partials()
+
   std::vector<material_interface_type*> m_execution_order;
   std::vector<property_base*> m_property_execution_order;
   std::unordered_map<std::string, property_base*> m_key_to_prop;  // for string-based find
@@ -723,6 +954,7 @@ class property_engine {
   bool m_verify{NUMSIM_MATERIALS_DEBUG_CHECKS != 0};
   bool m_fold{true};
   bool m_running{false}; ///< an update is in progress
+  bool m_evaluated{false}; ///< an update() ran to its end since the last failure
   bool m_failed{false};  ///< the last update threw
   std::unique_ptr<bool> m_check_reads{std::make_unique<bool>(true)}; ///< stable address for the wrapped callbacks
 };
diff --git a/include/numsim-materials/driver/material_point.h b/include/numsim-materials/driver/material_point.h
index 4eb4b29..ca301e5 100644
--- a/include/numsim-materials/driver/material_point.h
+++ b/include/numsim-materials/driver/material_point.h
@@ -131,6 +131,28 @@ class material_point {
     return output<T>(src.material, src.property);
   }
 
+  /// d output / d input, filled by evaluate(): output "material::property", input an input of this point (its history
+  /// output: the strain, the temperature) or "material::property". T is a scalar for two scalars, else a tensor of
+  /// the summed rank (stress by strain: rank 4).
+  template <typename T>
+    requires derivatives::derivative_value<T, typename Traits::value_type, Traits::Dim>
+  T const &derivative(std::string const &output, std::string const &input) const {
+    auto const out{connection_source::parse(output)};
+    find(out.material, out.property);
+    auto const seed{seed_of(input)};
+    find(seed.material, seed.property);
+    return m_ctx.template derivative<T>(out.material, out.property, seed.material, seed.property);
+  }
+
+  /// Whether the input (as in derivative()) is upstream of the output: otherwise derivative() is zero.
+  bool has_path(std::string const &output, std::string const &input) const {
+    auto const out{connection_source::parse(output)};
+    find(out.material, out.property);
+    auto const seed{seed_of(input)};
+    find(seed.material, seed.property);
+    return m_ctx.has_path(out.material, out.property, seed.material, seed.property);
+  }
+
   /// A history output, for both ends of the step.
   template <typename T>
   history_property<T> const &history(std::string const &material, std::string const &property) const {
@@ -198,6 +220,22 @@ class material_point {
     return parameters;
   }
 
+  /// An input's history output, or "material::property".
+  connection_source seed_of(std::string const &input) const {
+    if (input.contains("::"))
+      return connection_source::parse(input);
+    if (!std::ranges::contains(m_inputs, input))
+      throw umat::fatal_error("material_point: '" + input + "' is not an input of this material point");
+    std::vector<std::string> history;
+    for (auto const *prop : m_ctx.property_execution_order())
+      if (prop->is_history() && prop->traits().id.owner == input)
+        history.push_back(prop->traits().id.name);
+    if (history.size() != 1)
+      throw umat::fatal_error("material_point: input '" + input + "' has " + std::to_string(history.size()) +
+                              " history outputs; name the one to differentiate by as 'material::property'");
+    return {input, history.front()};
+  }
+
   property_base *find(std::string const &material, std::string const &property) const {
     if (auto *prop = m_ctx.find_property(material, property))
       return prop;
diff --git a/include/numsim-materials/materials/linear_elasticity.h b/include/numsim-materials/materials/linear_elasticity.h
index 8111a6f..82acd96 100755
--- a/include/numsim-materials/materials/linear_elasticity.h
+++ b/include/numsim-materials/materials/linear_elasticity.h
@@ -42,7 +42,9 @@ class linear_elasticity
   {
     update_tangent(); // readable before the first update
     base::add_callback(&linear_elasticity::update_tangent, m_C).reads();
-    base::add_callback(&linear_elasticity::update_stress, m_sig).reads(m_C, m_eps, m_eigenstrain);
+    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);
   }
 
   /**
@@ -69,6 +71,19 @@ class linear_elasticity
       m_sig = tmech::dcontract(m_C, tensor2{m_eps.get() - m_eigenstrain->get()});
     else
       m_sig = tmech::dcontract(m_C, m_eps.get());
+    if (base::want_partials())
+      update_partials();
+  }
+
+  /// d sigma / d eps = C, d sigma / d eps_eig = -C, d sigma / d C = Isym (x) (eps - eps_eig).
+  void update_partials() {
+    base::partial(m_sig, m_eps) = m_C;
+    if (m_eigenstrain) {
+      base::partial(m_sig, *m_eigenstrain) = value_type{-1} * m_C;
+      base::partial(m_sig, m_C) = tmech::otimes(sym_identity<value_type, Dim>(), tensor2{m_eps.get() - m_eigenstrain->get()});
+    } else {
+      base::partial(m_sig, m_C) = tmech::otimes(sym_identity<value_type, Dim>(), m_eps.get());
+    }
   }
 
   void update_tangent(){
diff --git a/include/numsim-materials/materials/thermal_eigenstrain.h b/include/numsim-materials/materials/thermal_eigenstrain.h
index edae4f6..060c1f5 100644
--- a/include/numsim-materials/materials/thermal_eigenstrain.h
+++ b/include/numsim-materials/materials/thermal_eigenstrain.h
@@ -30,7 +30,7 @@ class thermal_eigenstrain final : public material_base<thermal_eigenstrain<Trait
         m_T_ref(base::template get_parameter<value_type>("reference_temperature")),
         m_T(base::template add_input_history<value_type>(
             connection_source{base::template get_parameter<std::string>("temperature_source"), "state"}, EdgeKind::Global)) {
-    base::add_callback(&thermal_eigenstrain::compute, m_eigenstrain).reads(m_alpha, m_T);
+    base::add_callback(&thermal_eigenstrain::compute, m_eigenstrain).reads(m_alpha, m_T).partials(m_eigenstrain, m_alpha, m_T);
   }
 
   static input_parameter_controller parameters() {
@@ -41,7 +41,13 @@ class thermal_eigenstrain final : public material_base<thermal_eigenstrain<Trait
     return para;
   }
 
-  void compute() noexcept { m_eigenstrain = m_alpha.get() * (m_T.new_value() - m_T_ref) * tmech::eye<value_type, Dim, 2>(); }
+  void compute() {
+    m_eigenstrain = m_alpha.get() * (m_T.new_value() - m_T_ref) * tmech::eye<value_type, Dim, 2>();
+    if (base::want_partials()) {
+      base::partial(m_eigenstrain, m_alpha) = (m_T.new_value() - m_T_ref) * tmech::eye<value_type, Dim, 2>();
+      base::partial(m_eigenstrain, m_T) = m_alpha.get() * tmech::eye<value_type, Dim, 2>();
+    }
+  }
 
 private:
   tensor2 &m_eigenstrain;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 1844164..85f425b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -81,6 +81,7 @@ add_numsim_test(test_orientation test_orientation.cpp)
 add_numsim_test(test_short_fibre_composite test_short_fibre_composite.cpp)
 add_numsim_test(test_incremental_elasticity test_incremental_elasticity.cpp)
 add_numsim_test(test_thermal_eigenstrain test_thermal_eigenstrain.cpp)
+add_numsim_test(test_partial_derivatives test_partial_derivatives.cpp)
 add_numsim_test(test_kamal_sourour_rate test_kamal_sourour_rate.cpp)
 add_numsim_test(test_curing_eigenstrain test_curing_eigenstrain.cpp)
 add_numsim_test(test_damage_laws test_damage_laws.cpp)
diff --git a/tests/test_material_conformance.cpp b/tests/test_material_conformance.cpp
index 5d94152..fe2b539 100644
--- a/tests/test_material_conformance.cpp
+++ b/tests/test_material_conformance.cpp
@@ -5,8 +5,10 @@
 #include <regex>
 #include <set>
 #include <string>
+#include "numsim-materials/core/partial_check.h"
 #include "numsim-materials/default_materials.h"
 #include "numsim-materials/io/json_material_factory.h"
+#include "numsim-materials/umat/external_state_source.h"
 
 /// Generic conformance checks over every material, so a new one is covered
 /// without a bespoke test file.
@@ -153,4 +155,62 @@ TEST(MaterialConformance, MissingRequiredParametersThrow) {
   }
 }
 
+/// Materials whose callbacks declare partials or independence (.partials(), .independent_of()), checked by finite differences below.
+const std::set<std::string> kPartialsChecked = {"linear_elasticity", "thermal_eigenstrain"};
+
+/// Header stems under materials/ and solvers/ that declare partials or independent reads.
+std::set<std::string> headers_declaring_partials() {
+  namespace fs = std::filesystem;
+  std::set<std::string> out;
+  for (const char* directory : {NUMSIM_MATERIALS_INCLUDE_DIR, NUMSIM_SOLVERS_INCLUDE_DIR})
+    for (const auto& e : fs::directory_iterator(fs::path{directory})) {
+      if (e.path().extension() != ".h") continue;
+      std::ifstream in(e.path());
+      const std::string src{std::istreambuf_iterator<char>(in), std::istreambuf_iterator<char>()};
+      if (src.find(".partials(") != std::string::npos || src.find(".independent_of(") != std::string::npos)
+        out.insert(e.path().stem().string());
+    }
+  return out;
+}
+
+TEST(MaterialConformance, EveryMaterialDeclaringPartialsIsChecked) {
+  EXPECT_EQ(headers_declaring_partials(), kPartialsChecked)
+      << "a material declaring partials or independent reads joins kPartialsChecked and the model of DeclaredPartialsMatchFiniteDifferences";
+}
+
+/// Every declared partial against central differences of its own callback, at a state off every symmetry.
+TEST(MaterialConformance, DeclaredPartialsMatchFiniteDifferences) {
+  nm::material_context<policy> ctx;
+  for (const auto* source : {"strain", "temperature"}) {
+    policy::ParameterHandler p;
+    p.insert<std::string>("name", source);
+    if (std::string{source} == "strain")
+      ctx.create<nm::external_strain_source<policy>>(p);
+    else
+      ctx.create<nm::external_scalar_source<policy>>(p);
+  }
+  for (const auto* document : {
+           R"({"type": "constant_scalar", "name": "alpha", "value": 1.2e-5})",
+           R"({"type": "thermal_eigenstrain", "name": "thermal_eigenstrain", "alpha_source": "alpha::value", "reference_temperature": 20.0})",
+           R"({"type": "linear_elasticity", "name": "linear_elasticity", "strain_producer_name": "strain", "eigenstrain_source": "thermal_eigenstrain", "K": 150.0, "G": 70.0})"}) {
+    const nlohmann::json j = nlohmann::json::parse(document);
+    nm::create_from_json<policy>(ctx, j);
+  }
+  ctx.finalize();
+  tmech::tensor<double, 3, 2> strain;
+  strain(0, 0) = 1e-3, strain(1, 1) = -4e-4, strain(2, 2) = 2e-4;
+  strain(0, 1) = strain(1, 0) = 3e-4, strain(0, 2) = strain(2, 0) = -1e-4, strain(1, 2) = strain(2, 1) = 5e-4;
+  dynamic_cast<nm::external_strain_source<policy>*>(ctx.find("strain"))->bind({}, strain);
+  dynamic_cast<nm::external_scalar_source<policy>*>(ctx.find("temperature"))->bind(20, 75);
+  ctx.update();
+  std::set<std::string> checked;
+  for (const auto* material : ctx.materials())
+    if (!material->partials().empty()) {
+      checked.insert(material->name());
+      for (const auto& mismatch : nm::check_partials(ctx, material->name()))
+        ADD_FAILURE() << mismatch;
+    }
+  EXPECT_EQ(checked, kPartialsChecked);
+}
+
 }  // namespace
diff --git a/tests/test_partial_derivatives.cpp b/tests/test_partial_derivatives.cpp
new file mode 100644
index 0000000..ba05c0a
--- /dev/null
+++ b/tests/test_partial_derivatives.cpp
@@ -0,0 +1,780 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include <functional>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <tmech/tmech.h>
+#include "numsim-materials/core/material_base.h"
+#include "numsim-materials/core/material_context.h"
+#include "numsim-materials/core/parameter_helpers.h"
+#include "numsim-materials/core/partial_check.h"
+#include "numsim-materials/core/tensor_identities.h"
+#include "numsim-materials/driver/material_point.h"
+#include "numsim-materials/materials/constant_scalar.h"
+#include "numsim-materials/materials/linear_elasticity.h"
+#include "numsim-materials/materials/thermal_eigenstrain.h"
+#include "numsim-materials/umat/external_state_source.h"
+
+// Partial derivatives per callback and their forward accumulation (#217).
+
+namespace {
+
+namespace nm = numsim::materials;
+using policy = nm::material_policy_default;
+using T = policy::value_type;
+using tensor2 = tmech::tensor<T, 3, 2>;
+using tensor4 = tmech::tensor<T, 3, 4>;
+using tensor6 = tmech::tensor<T, 3, 6>;
+using ctx_type = nm::material_context<policy>;
+using param_type = policy::ParameterHandler;
+
+tensor2 const I{tmech::eye<T, 3, 2>()};
+tensor4 const Is{nm::sym_identity<T, 3>()};
+
+tensor2 sym(T a, T b, T c, T d, T e, T f) {
+  tensor2 x;
+  x(0, 0) = a, x(1, 1) = b, x(2, 2) = c;
+  x(1, 2) = x(2, 1) = d, x(0, 2) = x(2, 0) = e, x(0, 1) = x(1, 0) = f;
+  return x;
+}
+tensor2 const B2{sym(1, -2, 0.5, 0.3, -0.7, 1.1)};
+tensor2 const C2{sym(0.4, 0.2, -1, 0.6, 0.1, -0.3)};
+tensor2 const E2{sym(-0.5, 1.5, 0.25, -0.2, 0.9, 0.35)};
+tensor4 const A4{tensor4{2 * Is + 0.5 * tmech::otimes(I, I) + tmech::otimes(B2, B2)}};
+tensor4 const B4{tensor4{3 * Is - tmech::otimes(I, I) + 0.5 * tmech::otimes(C2, C2)}};
+
+// One law per edge of the contraction table: value and its partial in tmech.
+struct norm_trace { // scalar of tensor2
+  using in = tensor2;
+  using out = T;
+  static T value(tensor2 const &x) { return tmech::dcontract(x, x) + tmech::trace(x); }
+  static tensor2 partial(tensor2 const &x) { return 2 * x + I; }
+};
+struct square { // scalar of scalar
+  using in = T;
+  using out = T;
+  static T value(T s) { return s * s; }
+  static T partial(T s) { return 2 * s; }
+};
+struct scale_b { // tensor2 of scalar
+  using in = T;
+  using out = tensor2;
+  static tensor2 value(T s) { return s * B2; }
+  static tensor2 partial(T) { return B2; }
+};
+struct map_a { // tensor2 of tensor2
+  using in = tensor2;
+  using out = tensor2;
+  static tensor2 value(tensor2 const &x) { return tmech::dcontract(A4, x); }
+  static tensor4 partial(tensor2 const &) { return A4; }
+};
+struct map_b {
+  using in = tensor2;
+  using out = tensor2;
+  static tensor2 value(tensor2 const &x) { return tmech::dcontract(B4, x); }
+  static tensor4 partial(tensor2 const &) { return B4; }
+};
+struct project_c { // scalar of tensor2, linear
+  using in = tensor2;
+  using out = T;
+  static T value(tensor2 const &x) { return tmech::dcontract(C2, x); }
+  static tensor2 partial(tensor2 const &) { return C2; }
+};
+struct dyad { // tensor4 of tensor2: Y = tr(x) C (x) x
+  using in = tensor2;
+  using out = tensor4;
+  static tensor4 value(tensor2 const &x) { return tmech::trace(x) * tmech::otimes(C2, x); }
+  static tensor6 partial(tensor2 const &x) {
+    return tensor6{tmech::trace(x) * tmech::otimes(C2, Is) + tmech::otimes(tensor4{tmech::otimes(C2, x)}, I)};
+  }
+};
+struct contract_e { // tensor2 of tensor4: z = Y : E
+  using in = tensor4;
+  using out = tensor2;
+  static tensor2 value(tensor4 const &Y) { return tmech::dcontract(Y, E2); }
+  static tensor6 partial(tensor4 const &) { return tmech::otimes(Is, E2); }
+};
+struct constant_law { // tensor2 of tensor2, independent of its read
+  using in = tensor2;
+  using out = tensor2;
+  static tensor2 value(tensor2 const &) { return B2; }
+  static tensor4 partial(tensor2 const &) { return tensor4{}; }
+};
+// y = c tr(x)^3, about 1e8 at a strain of 1e-3; the partial is off by `error` (relative).
+template <int Error>
+struct cubic_strain {
+  using in = tensor2;
+  using out = T;
+  static constexpr T c{1e19};
+  static T value(tensor2 const &x) { return c * std::pow(tmech::trace(x), 3); }
+  static tensor2 partial(tensor2 const &x) {
+    return tensor2{(1 + Error * T{1e-6}) * 3 * c * std::pow(tmech::trace(x), 2) * I};
+  }
+};
+struct throwing_map { // throws once the read leaves the state it was armed at
+  using in = tensor2;
+  using out = tensor2;
+  static inline tensor2 armed{};
+  static inline bool active{false};
+  static tensor2 value(tensor2 const &x) {
+    if (active && tmech::norm(tensor2{x - armed}) > 0)
+      throw std::runtime_error("perturbed state");
+    return tmech::dcontract(A4, x);
+  }
+  static tensor4 partial(tensor2 const &) { return A4; }
+};
+struct wrong_map { // declares A4 for B4 : x
+  using in = tensor2;
+  using out = tensor2;
+  static tensor2 value(tensor2 const &x) { return tmech::dcontract(B4, x); }
+  static tensor4 partial(tensor2 const &) { return A4; }
+};
+
+/// y = Law::value(x) for x read from "x_source"; "partials" is "declared" (the default), "independent" or "none".
+template <typename Law, typename Traits>
+class law_node final : public nm::material_base<law_node<Law, Traits>, Traits> {
+public:
+  using base = nm::material_base<law_node<Law, Traits>, Traits>;
+  using input_parameter_controller = typename base::input_parameter_controller;
+
+  template <typename... Args>
+  explicit law_node(Args &&...args)
+      : base(std::forward<Args>(args)...), m_y(base::template add_output<typename Law::out>("y")),
+        m_source(nm::connection_source::parse(base::template get_parameter<std::string>("x_source"))),
+        m_x(base::template add_input<typename Law::in>(m_source.material, m_source.property, nm::EdgeKind::Global)),
+        m_mode(base::template get_parameter<std::string>("partials")) {
+    auto handle{base::add_callback(&law_node::compute, m_y)};
+    handle.reads(m_x);
+    if (m_mode == "declared")
+      handle.partials(m_y, m_x);
+    else if (m_mode == "independent")
+      handle.independent_of(m_x);
+  }
+
+  static input_parameter_controller parameters() {
+    input_parameter_controller para{base::parameters()};
+    nm::required<std::string>(para, "x_source");
+    nm::optional<std::string>(para, "partials", std::string{"declared"});
+    return para;
+  }
+
+  static inline std::function<void()> hook; ///< called by every run, when set
+
+  void compute() {
+    ++m_runs;
+    if (hook)
+      hook();
+    m_y = Law::value(m_x.get());
+    if (base::want_partials()) {
+      base::partial(m_y, m_x) = Law::partial(m_x.get());
+      ++m_fills;
+    }
+  }
+
+  int runs() const noexcept { return m_runs; }
+  int fills() const noexcept { return m_fills; }
+
+private:
+  typename Law::out &m_y;
+  nm::connection_source const m_source;
+  nm::input_property<typename Law::in, nm::property_traits> const &m_x;
+  std::string const &m_mode;
+  int m_runs{0};
+  int m_fills{0};
+};
+
+/// Declares its partials wrongly, as "mistake" says.
+template <typename Traits>
+class misdeclared final : public nm::material_base<misdeclared<Traits>, Traits> {
+public:
+  using base = nm::material_base<misdeclared<Traits>, Traits>;
+  using input_parameter_controller = typename base::input_parameter_controller;
+
+  template <typename... Args>
+  explicit misdeclared(Args &&...args)
+      : base(std::forward<Args>(args)...), m_y(base::template add_output<tensor2>("y")),
+        m_z(base::template add_output<tensor2>("z")),
+        m_x(base::template add_input<tensor2>("strain", "strain", nm::EdgeKind::Global)),
+        m_w(base::template add_input<T>("temperature", "state", nm::EdgeKind::Global)) {
+    auto const &mistake{base::template get_parameter<std::string>("mistake")};
+    base::add_callback(&misdeclared::compute_z, m_z).reads(m_x);
+    auto handle{base::add_callback(&misdeclared::compute, m_y)};
+    if (mistake == "before_reads")
+      handle.partials(m_y, m_x);
+    handle.reads(m_x);
+    if (mistake == "unread")
+      handle.partials(m_y, m_w);
+    else if (mistake == "twice")
+      handle.partials(m_y, m_x).partials(m_y, m_x);
+    else if (mistake == "both")
+      handle.partials(m_y, m_x).independent_of(m_x);
+    else if (mistake == "other_output")
+      handle.partials(m_z, m_x);
+  }
+
+  static input_parameter_controller parameters() {
+    input_parameter_controller para{base::parameters()};
+    nm::required<std::string>(para, "mistake");
+    return para;
+  }
+
+  void compute() { m_y = m_x.get(); }
+  void compute_z() { m_z = m_x.get(); }
+
+private:
+  tensor2 &m_y;
+  tensor2 &m_z;
+  nm::input_property<tensor2, nm::property_traits> const &m_x;
+  nm::input_property<T, nm::property_traits> const &m_w;
+};
+
+/// y = a::y + b::y, both read and differentiated.
+template <typename Traits>
+class sum_node final : public nm::material_base<sum_node<Traits>, Traits> {
+public:
+  using base = nm::material_base<sum_node<Traits>, Traits>;
+  using input_parameter_controller = typename base::input_parameter_controller;
+
+  template <typename... Args>
+  explicit sum_node(Args &&...args)
+      : base(std::forward<Args>(args)...), m_y(base::template add_output<tensor2>("y")),
+        m_a(base::template add_input<tensor2>("a", "y", nm::EdgeKind::Global)),
+        m_b(base::template add_input<tensor2>("b", "y", nm::EdgeKind::Global)) {
+    base::add_callback(&sum_node::compute, m_y).reads(m_a, m_b).partials(m_y, m_a, m_b);
+  }
+
+  static input_parameter_controller parameters() { return input_parameter_controller{base::parameters()}; }
+
+  void compute() {
+    m_y = m_a.get() + m_b.get();
+    if (base::want_partials()) {
+      base::partial(m_y, m_a) = Is;
+      base::partial(m_y, m_b) = Is;
+    }
+  }
+
+private:
+  tensor2 &m_y;
+  nm::input_property<tensor2, nm::property_traits> const &m_a;
+  nm::input_property<tensor2, nm::property_traits> const &m_b;
+};
+
+/// y = (a::y : a::y) b::y: the partial with respect to b depends on a.
+template <typename Traits>
+class product_node final : public nm::material_base<product_node<Traits>, Traits> {
+public:
+  using base = nm::material_base<product_node<Traits>, Traits>;
+  using input_parameter_controller = typename base::input_parameter_controller;
+
+  template <typename... Args>
+  explicit product_node(Args &&...args)
+      : base(std::forward<Args>(args)...), m_y(base::template add_output<tensor2>("y")),
+        m_a(base::template add_input<tensor2>("a", "y", nm::EdgeKind::Global)),
+        m_b(base::template add_input<tensor2>("b", "y", nm::EdgeKind::Global)) {
+    base::add_callback(&product_node::compute, m_y).reads(m_a, m_b).partials(m_y, m_a, m_b);
+  }
+
+  static input_parameter_controller parameters() { return input_parameter_controller{base::parameters()}; }
+
+  void compute() {
+    auto const aa{tmech::dcontract(m_a.get(), m_a.get())};
+    m_y = aa * m_b.get();
+    if (base::want_partials()) {
+      base::partial(m_y, m_a) = tensor4{T{2} * tmech::otimes(m_b.get(), m_a.get())};
+      base::partial(m_y, m_b) = tensor4{aa * Is};
+    }
+  }
+
+private:
+  tensor2 &m_y;
+  nm::input_property<tensor2, nm::property_traits> const &m_a;
+  nm::input_property<tensor2, nm::property_traits> const &m_b;
+};
+
+/// A context with the sources "strain" and "temperature" and law nodes added by node().
+struct chain {
+  std::unique_ptr<ctx_type> owner{std::make_unique<ctx_type>()};
+  ctx_type &ctx{*owner};
+
+  chain() {
+    param_type p;
+    p.insert<std::string>("name", "strain");
+    ctx.create<nm::external_strain_source<policy>>(p);
+    p.clear();
+    p.insert<std::string>("name", "temperature");
+    ctx.create<nm::external_scalar_source<policy>>(p);
+  }
+
+  template <typename Law>
+  law_node<Law, policy> &node(std::string const &name, std::string const &x_source, std::string const &partials = "declared") {
+    param_type p;
+    p.insert<std::string>("name", name);
+    p.insert<std::string>("x_source", x_source);
+    p.insert<std::string>("partials", partials);
+    return ctx.create<law_node<Law, policy>>(p);
+  }
+
+  void load(tensor2 const &strain, T temperature) {
+    dynamic_cast<nm::external_strain_source<policy> *>(ctx.find("strain"))->bind(tensor2{}, strain);
+    dynamic_cast<nm::external_scalar_source<policy> *>(ctx.find("temperature"))->bind(T{0}, temperature);
+    ctx.update();
+  }
+};
+
+tensor2 const x0{sym(0.3, -0.1, 0.2, 0.05, -0.15, 0.25)};
+constexpr T s0{1.7};
+
+template <typename A, typename B>
+T distance(A const &a, B const &b) {
+  if constexpr (std::is_arithmetic_v<A>)
+    return std::abs(a - b);
+  else
+    return tmech::norm(a - b);
+}
+
+// --- the contraction table, one hand-made chain per row ---
+
+TEST(PartialDerivatives, ScalarOfScalarOfTensor) {
+  chain c;
+  c.node<norm_trace>("y", "strain::strain");
+  c.node<square>("z", "y::y");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor2>("z", "y", "strain", "strain")};
+  c.load(x0, s0);
+  tensor2 const expected{2 * norm_trace::value(x0) * norm_trace::partial(x0)};
+  EXPECT_LT(distance(d, expected), 1e-12);
+}
+
+TEST(PartialDerivatives, TensorOfScalarOfTensor) {
+  chain c;
+  c.node<norm_trace>("y", "strain::strain");
+  c.node<scale_b>("z", "y::y");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor4>("z", "y", "strain", "strain")};
+  c.load(x0, s0);
+  EXPECT_LT(distance(d, tensor4{tmech::otimes(B2, norm_trace::partial(x0))}), 1e-12);
+}
+
+TEST(PartialDerivatives, TensorOfTensorOfTensor) {
+  chain c;
+  c.node<map_a>("y", "strain::strain");
+  c.node<map_b>("z", "y::y");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor4>("z", "y", "strain", "strain")};
+  c.load(x0, s0);
+  EXPECT_LT(distance(d, tensor4{tmech::dcontract(B4, A4)}), 1e-12);
+}
+
+TEST(PartialDerivatives, ScalarOfTensorOfTensor) {
+  chain c;
+  c.node<map_a>("y", "strain::strain");
+  c.node<project_c>("z", "y::y");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor2>("z", "y", "strain", "strain")};
+  c.load(x0, s0);
+  EXPECT_LT(distance(d, tensor2{tmech::dcontract(C2, A4)}), 1e-12);
+}
+
+TEST(PartialDerivatives, TensorOfTensorOfScalar) {
+  chain c;
+  c.node<scale_b>("y", "temperature::state");
+  c.node<map_a>("z", "y::y");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor2>("z", "y", "temperature", "state")};
+  c.load(x0, s0);
+  EXPECT_LT(distance(d, tensor2{tmech::dcontract(A4, B2)}), 1e-12);
+}
+
+TEST(PartialDerivatives, TensorOfFourthOrderTensorOfTensor) {
+  chain c;
+  c.node<dyad>("y", "strain::strain");
+  c.node<contract_e>("z", "y::y");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor4>("z", "y", "strain", "strain")};
+  c.load(x0, s0);
+  // z = tr(x) (x : E) C
+  tensor4 const expected{tmech::trace(x0) * tmech::otimes(C2, E2) + tmech::dcontract(x0, E2) * tmech::otimes(C2, I)};
+  EXPECT_LT(distance(d, expected), 1e-12);
+}
+
+TEST(PartialDerivatives, ScalarOfScalarIsAScalar) {
+  chain c;
+  c.node<square>("y", "temperature::state");
+  c.node<square>("z", "y::y");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<T>("z", "y", "temperature", "state")};
+  c.load(x0, s0);
+  EXPECT_NEAR(d, 4 * s0 * s0 * s0, 1e-12); // z = s^4
+}
+
+TEST(PartialDerivatives, ParallelPathsAreSummed) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.node<map_b>("b", "strain::strain");
+  param_type p;
+  p.insert<std::string>("name", "s");
+  c.ctx.create<sum_node<policy>>(p);
+  c.node<project_c>("z", "s::y");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor2>("z", "y", "strain", "strain")};
+  c.load(x0, s0);
+  EXPECT_LT(distance(d, tensor2{tmech::dcontract(C2, tensor4{A4 + B4})}), 1e-12);
+}
+
+TEST(PartialDerivatives, NoPathGivesZero) {
+  chain c;
+  c.node<map_a>("y", "strain::strain");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor2>("y", "y", "temperature", "state")};
+  c.load(x0, s0);
+  EXPECT_EQ(tmech::norm(d), 0);
+}
+
+TEST(PartialDerivatives, TheResultFollowsTheState) {
+  chain c;
+  c.node<norm_trace>("y", "strain::strain");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor2>("y", "y", "strain", "strain")};
+  for (T const f : {1.0, -2.0, 0.5}) {
+    c.load(tensor2{f * x0}, s0);
+    EXPECT_LT(distance(d, norm_trace::partial(tensor2{f * x0})), 1e-12);
+  }
+}
+
+// --- declarations and their errors ---
+
+TEST(PartialDerivatives, AMissingPartialOnARequestedPathThrowsAndNamesThePath) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.node<map_b>("b", "a::y", "none");
+  c.node<map_a>("c", "b::y");
+  c.ctx.finalize();
+  try {
+    c.ctx.derivative<tensor4>("c", "y", "strain", "strain");
+    FAIL() << "a requested path without a partial was accepted";
+  } catch (std::logic_error const &e) {
+    std::string const what{e.what()};
+    EXPECT_TRUE(what.contains("'b::y'")) << what;
+    EXPECT_TRUE(what.contains("strain::strain -> a::y -> b::y -> c::y")) << what;
+  }
+}
+
+TEST(PartialDerivatives, AMissingPartialOffTheRequestedPathsDoesNotThrow) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.node<map_b>("b", "a::y", "none");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor4>("a", "y", "strain", "strain")};
+  EXPECT_NO_THROW(c.load(x0, s0));
+  EXPECT_LT(distance(d, A4), 1e-12);
+}
+
+TEST(PartialDerivatives, IndependentOfDeclaresAZeroPartial) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.node<map_b>("b", "a::y", "independent");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor4>("b", "y", "strain", "strain")};
+  c.load(x0, s0);
+  EXPECT_EQ(tmech::norm(d), 0);
+}
+
+TEST(PartialDerivatives, TheResultTypeMustMatchTheRanks) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.ctx.finalize();
+  EXPECT_THROW(c.ctx.derivative<tensor2>("a", "y", "strain", "strain"), std::logic_error);
+  EXPECT_THROW(c.ctx.derivative<tensor4>("a", "y", "strain", "nothing"), std::runtime_error);
+}
+
+TEST(PartialDerivatives, PartialsAreDeclaredForReadsAndOutputsOfTheCallback) {
+  for (auto const *mistake : {"unread", "twice", "both", "other_output", "before_reads"}) {
+    chain c;
+    param_type p;
+    p.insert<std::string>("name", "m");
+    p.insert<std::string>("mistake", mistake);
+    EXPECT_THROW(c.ctx.create<misdeclared<policy>>(p), std::logic_error) << mistake;
+  }
+}
+
+// --- independent_of() is verified, not trusted (review of #227) ---
+
+TEST(PartialDerivatives, AWrongIndependentOfFailsTheFiniteDifferenceCheck) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.node<map_b>("liar", "a::y", "independent"); // B4 : x does depend on x
+  c.ctx.finalize();
+  c.load(x0, s0);
+  auto const mismatches{nm::check_partials(c.ctx, "liar")};
+  ASSERT_FALSE(mismatches.empty());
+  EXPECT_TRUE(mismatches.front().contains("liar::y")) << mismatches.front();
+  EXPECT_TRUE(mismatches.front().contains("a::y")) << mismatches.front();
+  EXPECT_TRUE(mismatches.front().contains("independent_of")) << mismatches.front();
+  EXPECT_LT(distance(c.ctx.get<tensor2>("liar", "y"), map_b::value(c.ctx.get<tensor2>("a", "y"))), 1e-15) << "the state is restored";
+}
+
+TEST(PartialDerivatives, ATrueIndependentOfPassesTheFiniteDifferenceCheck) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.node<constant_law>("k", "a::y", "independent");
+  c.ctx.finalize();
+  c.load(x0, s0);
+  EXPECT_TRUE(nm::check_partials(c.ctx, "k").empty());
+}
+
+// --- the step follows the read, the allowance the analytic block (review of #227) ---
+
+TEST(PartialDerivatives, ACorrectPartialPassesAtTheScaleOfAStrain) {
+  chain c;
+  c.node<cubic_strain<0>>("n", "strain::strain");
+  c.ctx.finalize();
+  c.load(tensor2{1e-3 * x0}, s0);
+  ASSERT_GT(std::abs(c.ctx.get<T>("n", "y")), 1e7);
+  auto const mismatches{nm::check_partials(c.ctx, "n")};
+  EXPECT_TRUE(mismatches.empty()) << mismatches.front();
+}
+
+TEST(PartialDerivatives, ASmallRelativeErrorIsDetectedAtOutputsOfOneHundredMillion) {
+  chain c;
+  c.node<cubic_strain<20>>("n", "strain::strain"); // 2e-5 relative
+  c.ctx.finalize();
+  c.load(tensor2{1e-3 * x0}, s0);
+  ASSERT_GT(std::abs(c.ctx.get<T>("n", "y")), 1e7);
+  EXPECT_FALSE(nm::check_partials(c.ctx, "n").empty());
+}
+
+// --- a throwing rerun leaves the context as it was (review of #227) ---
+
+TEST(PartialDerivatives, AThrowingPerturbedRerunRestoresReadsAndOutputs) {
+  chain c;
+  c.node<throwing_map>("t", "strain::strain");
+  c.ctx.finalize();
+  c.load(x0, s0);
+  throwing_map::armed = x0;
+  throwing_map::active = true;
+  EXPECT_THROW(nm::check_partials(c.ctx, "t"), std::runtime_error);
+  throwing_map::active = false;
+  EXPECT_EQ(tmech::norm(tensor2{dynamic_cast<nm::external_strain_source<policy> *>(c.ctx.find("strain"))->strain().new_value() - x0}), 0) << "the read";
+  EXPECT_LT(distance(c.ctx.get<tensor2>("t", "y"), throwing_map::value(x0)), 1e-15) << "the output";
+}
+
+// --- the result type is the one the context can fill (review of #227) ---
+
+template <typename R>
+concept derivable = requires(ctx_type &c, std::string const &name) { c.template derivative<R>(name, name, name, name); };
+
+TEST(PartialDerivatives, OnlyTheScalarAndTensorsOfTheContextAreDerivativeTypes) {
+  static_assert(derivable<T>);
+  static_assert(derivable<tensor2>);
+  static_assert(derivable<tensor4>);
+  static_assert(derivable<tensor6>);
+  static_assert(!derivable<tmech::tensor<T, 2, 4>>); // another dimension would write out of bounds
+  static_assert(!derivable<tmech::tensor<float, 3, 2>>);
+  static_assert(!derivable<float>);
+  static_assert(!derivable<int>);
+  static_assert(!derivable<std::string>);
+}
+
+
+TEST(PartialDerivatives, AFourthOrderResultOfATensorIsRankSix) {
+  chain c;
+  c.node<dyad>("d", "strain::strain");
+  c.ctx.finalize();
+  auto const &d{c.ctx.derivative<tensor6>("d", "y", "strain", "strain")};
+  c.load(x0, s0);
+  EXPECT_LT(distance(d, dyad::partial(x0)), 1e-12);
+}
+
+// --- derivative() on an evaluated context, and the path query (review of #227) ---
+
+TEST(PartialDerivatives, ADerivativeRequestedAfterAnUpdateIsFilledAtOnce) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.node<map_b>("b", "a::y");
+  c.ctx.finalize();
+  c.load(x0, s0);
+  auto const &d{c.ctx.derivative<tensor4>("b", "y", "strain", "strain")};
+  EXPECT_LT(distance(d, tensor4{tmech::dcontract(B4, A4)}), 1e-12);
+  c.load(tensor2{2 * x0}, s0);
+  EXPECT_LT(distance(d, tensor4{tmech::dcontract(B4, A4)}), 1e-12) << "and by the next update";
+}
+
+TEST(PartialDerivatives, ADerivativeFromACallbackIsRefused) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.ctx.finalize();
+  bool refused{false};
+  law_node<map_a, policy>::hook = [&] {
+    try {
+      c.ctx.derivative<tensor4>("a", "y", "strain", "strain");
+    } catch (std::logic_error const &) {
+      refused = true;
+    }
+  };
+  c.load(x0, s0);
+  law_node<map_a, policy>::hook = nullptr;
+  EXPECT_TRUE(refused);
+}
+
+TEST(PartialDerivatives, HasPathTellsIndependenceFromATypo) {
+  chain c;
+  c.node<map_a>("y", "strain::strain");
+  c.ctx.finalize();
+  EXPECT_TRUE(c.ctx.has_path("y", "y", "strain", "strain"));
+  EXPECT_TRUE(c.ctx.has_path("y", "y", "y", "y"));
+  EXPECT_FALSE(c.ctx.has_path("y", "y", "temperature", "state"));
+  EXPECT_FALSE(c.ctx.has_path("strain", "strain", "y", "y"));
+  EXPECT_THROW(c.ctx.has_path("y", "y", "temperature", "nothing"), std::runtime_error);
+}
+
+// --- zero cost when nothing is requested ---
+
+TEST(PartialDerivatives, NoPartialsAreComputedWithoutARequest) {
+  chain c;
+  auto &a{c.node<map_a>("a", "strain::strain")};
+  auto &b{c.node<map_b>("b", "a::y")};
+  auto &q{c.node<map_b>("q", "strain::strain")};
+  c.ctx.finalize();
+  for (int i{0}; i < 3; ++i)
+    c.load(x0, s0);
+  EXPECT_EQ(a.runs(), 3);
+  EXPECT_EQ(a.fills() + b.fills() + q.fills(), 0);
+
+  c.ctx.derivative<tensor4>("b", "y", "strain", "strain"); // fills once at the held state
+  EXPECT_EQ(a.fills(), 1);
+  EXPECT_EQ(b.fills(), 1);
+  for (int i{0}; i < 2; ++i)
+    c.load(x0, s0);
+  EXPECT_EQ(a.fills(), 3);
+  EXPECT_EQ(b.fills(), 3);
+  EXPECT_EQ(q.fills(), 0) << "q is off the requested path";
+}
+
+// --- the per-node finite-difference check ---
+
+TEST(PartialDerivatives, TheFiniteDifferenceCheckPassesCorrectPartials) {
+  chain c;
+  c.node<norm_trace>("n", "strain::strain");
+  c.node<square>("s", "n::y");
+  c.node<scale_b>("sb", "s::y");
+  c.node<map_a>("a", "sb::y");
+  c.node<project_c>("p", "a::y");
+  c.node<dyad>("d", "a::y");
+  c.node<contract_e>("e", "d::y");
+  c.ctx.finalize();
+  c.load(x0, s0);
+  for (auto const *name : {"n", "s", "sb", "a", "p", "d", "e"})
+    EXPECT_TRUE(nm::check_partials(c.ctx, name).empty()) << name;
+}
+
+TEST(PartialDerivatives, TheFiniteDifferenceCheckComparesWithThePartialsAtTheState) {
+  chain c;
+  c.node<map_a>("a", "strain::strain");
+  c.node<scale_b>("b", "temperature::state");
+  param_type p;
+  p.insert<std::string>("name", "q");
+  c.ctx.create<product_node<policy>>(p);
+  c.ctx.finalize();
+  c.ctx.derivative<tensor4>("q", "y", "strain", "strain"); // keeps the partials of q requested
+  c.load(x0, s0);
+  auto const mismatches{nm::check_partials(c.ctx, "q")};
+  EXPECT_TRUE(mismatches.empty()) << mismatches.front();
+}
+
+TEST(PartialDerivatives, TheFiniteDifferenceCheckNamesAWrongPartial) {
+  chain c;
+  c.node<wrong_map>("w", "strain::strain");
+  c.ctx.finalize();
+  c.load(x0, s0);
+  auto const mismatches{nm::check_partials(c.ctx, "w")};
+  ASSERT_FALSE(mismatches.empty());
+  EXPECT_TRUE(mismatches.front().contains("w::y")) << mismatches.front();
+  EXPECT_TRUE(mismatches.front().contains("strain::strain")) << mismatches.front();
+  EXPECT_LT(distance(c.ctx.get<tensor2>("w", "y"), wrong_map::value(x0)), 1e-15) << "the check restores the state";
+}
+
+// --- thermal eigenstrain with linear elasticity: the first DDSDDT ---
+
+constexpr T K{100}, G{40}, alpha{1e-5}, T_ref{20};
+
+struct thermoelastic {
+  std::unique_ptr<ctx_type> owner{std::make_unique<ctx_type>()};
+  ctx_type &ctx{*owner};
+  std::unique_ptr<nm::material_point<policy>> point;
+
+  thermoelastic() {
+    param_type p;
+    p.insert<std::string>("name", "strain");
+    ctx.create<nm::external_strain_source<policy>>(p);
+    p.clear();
+    p.insert<std::string>("name", "temperature");
+    ctx.create<nm::external_scalar_source<policy>>(p);
+    p.clear();
+    p.insert<std::string>("name", "alpha");
+    p.insert<T>("value", alpha);
+    ctx.create<nm::constant_scalar<policy>>(p);
+    p.clear();
+    p.insert<std::string>("name", "thermal");
+    p.insert<std::string>("alpha_source", "alpha::value");
+    p.insert<T>("reference_temperature", T_ref);
+    ctx.create<nm::thermal_eigenstrain<policy>>(p);
+    p.clear();
+    p.insert<std::string>("name", "elastic");
+    p.insert<std::string>("strain_producer_name", "strain");
+    p.insert<std::string>("eigenstrain_source", "thermal");
+    p.insert<T>("K", K);
+    p.insert<T>("G", G);
+    ctx.create<nm::linear_elasticity<policy>>(p);
+    ctx.finalize();
+    param_type q;
+    q.insert<std::vector<std::string>>("inputs", {"strain", "temperature"});
+    point = std::make_unique<nm::material_point<policy>>(ctx, q);
+  }
+
+  void evaluate(tensor2 const &strain, T temperature) {
+    point->reset();
+    point->input<nm::external_strain_source<policy>>("strain").bind(tensor2{}, strain);
+    point->input<nm::external_scalar_source<policy>>("temperature").bind(T_ref, temperature);
+    point->evaluate();
+  }
+};
+
+TEST(PartialDerivatives, StressWithRespectToTemperatureMatchesFiniteDifferences) {
+  thermoelastic m;
+  auto const &ddsddt{m.point->derivative<tensor2>("elastic::stress", "temperature")};
+  tensor2 const eps{sym(1e-3, -2e-4, 5e-4, 1e-4, -3e-4, 2e-4)};
+  T const theta{80}, h{1e-3};
+  m.evaluate(eps, theta + h);
+  tensor2 const plus{m.point->output<tensor2>("elastic::stress")};
+  m.evaluate(eps, theta - h);
+  tensor2 const minus{m.point->output<tensor2>("elastic::stress")};
+  m.evaluate(eps, theta);
+  tensor2 const fd{(plus - minus) / (2 * h)};
+  EXPECT_LT(tmech::norm(tensor2{ddsddt - fd}), 1e-8 * tmech::norm(fd));
+  EXPECT_LT(tmech::norm(tensor2{ddsddt + 3 * K * alpha * I}), 1e-12) << "-C : alpha I";
+}
+
+TEST(PartialDerivatives, StressWithRespectToStrainEqualsTheTangent) {
+  thermoelastic m;
+  auto const &ddsdde{m.point->derivative<tensor4>("elastic::stress", "strain")};
+  m.evaluate(sym(1e-3, -2e-4, 5e-4, 1e-4, -3e-4, 2e-4), 60);
+  auto const &tangent{m.point->output<tensor4>("elastic::tangent")};
+  EXPECT_LT(tmech::norm(tensor4{ddsdde - tangent}), 1e-12 * tmech::norm(tangent));
+}
+
+TEST(PartialDerivatives, TheMaterialPointNamesUnknownInputs) {
+  thermoelastic m;
+  EXPECT_THROW(m.point->derivative<tensor2>("elastic::stress", "pressure"), nm::umat::fatal_error);
+  EXPECT_THROW(m.point->derivative<tensor2>("elastic::nothing", "temperature"), nm::umat::fatal_error);
+}
+
+TEST(PartialDerivatives, ThermoelasticPartialsPassTheFiniteDifferenceCheck) {
+  thermoelastic m;
+  m.evaluate(sym(1e-3, -2e-4, 5e-4, 1e-4, -3e-4, 2e-4), 60);
+  for (auto const *name : {"thermal", "elastic"})
+    for (auto const &mismatch : nm::check_partials(m.ctx, name))
+      ADD_FAILURE() << mismatch;
+}
+
+} // namespace