For Dt(g(u)) with a nonlinear g (Richards moisture content theta(h) is
the canonical case), PR firedrakeproject#204 already builds the conservative
two-evaluation form g(U_i) - g(u0) inside getFormStage, so the per-stage
equations are correct for any RK method. Mass conservation still failed
for non-stiffly-accurate methods like ImplicitMidpoint and
GaussLegendre(2) because _update_Ainv computes u_new as a linear
combination of u0 and the stage values. For linear Dt(u) that linear
combination is algebraically equivalent to the conservative time
integration, but for theta nonlinear it is not, and the conservation
property the stage equations establish is destroyed at the update step.

This routes the non-SA non-Bernstein path through a new conservative
update solve when the form has a composite Dt: solve
  replace(F_dtless, {u0: unew}) - replace(F_dtless, {u0: u0})
       + dt * sum_i b_i * F_remainder(stage_i) = 0
for u_new. For g = identity it collapses to the old
inner(unew - u0, v)*dx head, so the linear case is unchanged. For
nonlinear g it is a discrete statement of the global conservation
invariant and u_new satisfies it by construction. Linear Dt(u) keeps
_update_Ainv since that path also handles DAE structure correctly,
which the conservative variational head does not (the algebraic block
is left unconstrained, leading to a singular linear solve).

The composite-Dt detection lives in irksome/ufl/manipulation.py as the
public helper has_composite_time_derivative(F, u0). It walks
TimeDerivative nodes, exempting structural views (Indexed, ListTensor,
ComponentTensor), linear differential operators (Grad, Div, Curl), DG
facet operators (CellAvg, FacetAvg, restrictions), and arithmetic that
is linear in u0 (Sum, Product/Division with at most one operand
depending on u0). Pure-time forcing terms Dt(f(t,x)) are passed
straight through to the chain rule. The docstring carries an explicit
warning that detection is syntactic: pre-interpolating an expression
in u0 into a Function and writing Dt(that_function) loses the
dependence and silently produces a non-conservative discretisation.

The DIRK stage-derivative path folds the chain rule into the
K-substitution and has no clean conservative analogue without
restructuring. Rather than silently mis-discretising, DIRKTimeStepper
now raises NotImplementedError when the form has a composite Dt and
points users at stage_type=value. use_collocation_update=True
similarly refuses composite Dt because the collocation polynomials
terminal value is also a linear combination of stages.

New tests in tests/test_conservative_dt.py cover the seven implicit
tableaux Irksome ships (BackwardEuler, RadauIIA(2), LobattoIIIC(2),
Alexander, ImplicitMidpoint, GaussLegendre(2), QinZhang) on the
existing Richards problem with mass error below 1e-10 over twenty
steps; the DIRK and use_collocation_update refusals; the
linear-arithmetic classification matrix (Dt(2u), Dt(u/2),
Dt(u + sin(t)) classified linear; Dt(u*u), Dt(theta(u)), Dt(1/u)
classified composite); and identity-runs-and-decays sanity across
stage_value and DIRK. Targeted regression on test_dirk, test_dae,
test_explicit, test_galerkin, test_disc_galerkin, test_split,
test_time_form_splitting, test_mass_conservation, test_pc, test_bern,
test_butcher, test_differentiation, test_curl all pass: 243 passed,
2 skipped, 0 failed.

Per Pablo's review:

* Rename has_composite_time_derivative -> has_nonlinear_time_derivative.
  "Composite" was internal jargon; "nonlinear" is what we are actually
  detecting and reads more directly at the call sites.  Test file
  renamed in lockstep.

* Simplify the conservative update head from
  replace(F_dtless, {u0: unew}) - replace(F_dtless, {u0: u0})
  to replace(F_dtless, {u0: unew}) - F_dtless: the second replacement
  is an identity transformation, so the two forms are byte-identical.
  Pablo's suggested phrasing.

* Move the explanation of the two-evaluation conservative head from
  inline comments into the docstring of get_update_solver, and state
  the math in user-facing terms (inner(g(unew) - g(u0), v) * dx)
  rather than via the internal F_dtless symbol.  Drop the redundant
  "subtracting the two replaced forms keeps v unchanged" comment.

* Compress the warm-start comment in _update_general to a single line
  and drop the application-specific (Haverkamp / moisture capacity)
  language.  The failure mode is generic: g prime of u0 can vanish,
  in which case a zero-initialised iterate gives a singular initial
  Jacobian.

* Move the update_solver_parameters = solver_parameters default into
  the conservative-update branch that actually consumes it.  Pablo
  pointed out the previous placement implied the inheritance always
  applies; in fact the linear _update_Ainv path does not use it, so
  scoping the default makes the intent clearer.

Per Pablo's review on manipulation.py: the bespoke walker with its
_PASS_THROUGH_OPS exemption list and _is_pass_through / _depends_on
helpers is a maintenance liability -- every linear UFL operator added
upstream silently misclassifies until Irksome's exemption list is
updated to match.

ufl.derivative already knows the linearity of Grad, Div, Indexed,
ListTensor, ComponentTensor, restrictions, CellAvg, FacetAvg, and any
future linear node UFL adds.  Reformulating the classifier as

    D = expand_derivatives(derivative(f, u0, TrialFunction(V)))
    if u0 in extract_coefficients(D): nonlinear

inherits that classification for free and collapses the walker to its
extract_type(F, TimeDerivative) loop plus three lines.

Verified against tests/test_has_nonlinear_time_derivative.py (all
seven cases agree with the walker) and tests/test_mass_conservation.py
(all nine tableaux still conserve to machine precision).

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>