Add non-elastic shear heating

[tan2]

Adding non-elastic shear heating to DES's rheology model: the elastic strain implied by a stress state is recovered by inverting isotropic Hooke's law at the start and end of each step, and what strain increment isn't accounted for elastically is dissipated as heat, evaluated against the trapezoidal-average stress. This works uniformly across the elastic/maxwell/plastic/viscoplastic branches since they share the same elastic moduli; the memoryless rh_viscous rheology has no elastic modulus tied to its stress update, so it gets the closed-form viscous dissipation formula instead. 2D plane-strain plastic branches fold in stressyy, matching geoFLAC's third principal stress.

Gated behind control.has_shear_heating (default off); the resulting per-element shear_heat field feeds update_temperature as a volumetric source term alongside radiogenic_source.

This feature is based on the discussion in Issue #61.

[echoi]

@tan2 This is a great improvement. Do you have an example cfg files that you used during the development? If so, please add them. A benchmark would be even better but at least an example would do if it clearly shows the effects of viscous and inelastic shear heating on the temperature field.

[tan2]

A simple example input file is added. It is a 20 km x 20 km elato-plastic box with uniform initial temperature, no thermal expansion, and no strain weakening.

The shear heating increases the mean temperature by 5.1 K after compressing at a constant rate 1e-9 m/s for 20,000 years.

[chaseshyu]

This is a valuable implementation to closes the energy budget between mechanical work and heat, which DynEarthSol has been missing. There are a few bugs to sort out first, though. I ran the branch through Claude Code, and the review below should help with them.

Following reviews from Claude Code might help with those fixes.

Reference fix branch: feat/shear-heating-fixes (https://github.com/GeoFLAC/DynEarthSol/tree/feat/shear-heating-fixes).

---

The design is right — inverting Hooke's law so one path covers all the branches, and the
one-step lag against update_temperature handled via the checkpoint. The formula itself
checks out: unclamped, its sum matches the analytic s':s'/(2η)·dt to seven digits.

Two things layered on top of it make the heat wrong. Quickest way in, @tan2: the +5.1 K
is what the current code gives; the correct value for that cfg is +4.1 K, and the
hottest element is +29.2 K rather than +38.1 K.

BLOCKER

1. rheology.cxx:380 — the per-component clamp destroys a cancellation that only holds in the sum. de_ne is deviatoric, so Σ sᵢᵢ·de_neᵢᵢ = Σ s'ᵢᵢ·de_neᵢᵢ; clamping
   each component against the total stress removes exactly the terms that cancel the
   pressure. Error scales as |p|/|s′|. Measured on 2-D elasto-visco-plastic at η=1e18
   (0 of 119 elements plastic, so purely the viscous winner — the production case):
   2.40e−5 vs 5.63e−7 W/m³, 43× too much; max ΔT +11.30 vs +0.09 K. 3-D maxwell
   gives 7.5×. Fix: drop the three clamps, add if (q < 0.) q = 0.; before return q.

2. rheology.cxx:368 — the 2-D Hooke inversion assumes σ_yy = 0, but DES's 2-D state is plane strain. 3λ+2μ needs the full 3-D trace, available only when
   mat.is_plane_strain is set — and it defaults to false. Leaves a
   −dv·λ²/(2μ(3λ+2μ)) residue per diagonal on any volume change, which the clamp then
   rectifies. Run your cfg with rheology_type = elastic: elasticity dissipates nothing,
   so the exact answer is max(shear_heat) == 0. As merged it is nonzero everywhere and
   climbs 500.000 → 500.404 K over 20 kyr with no bound; 3-D gives 1e−17, so it is
   2-D only. Fix: use 2λ+2μ when σ_yy is not stored — that run then returns 1.7e−15 and
   holds 500.00000 K.

MAJOR

3. 2vtk.py:307 / Dynearthsol.py:185 — the field never reaches ParaView.
   'shear heating' is not registered as an element scalar, so on the legacy binary
   format (the default) read_field raises NameError and the new except swallows it.
   Only the HDF5 path works. One string fixes it.

4. fields.cxx:171 — shear heat is wiped at every remesh. reallocate_variables
   re-news it to zero, nothing remaps it, and update_temperature runs before
   update_stress. Measured with a frame on each side of 270 remeshes: post-remesh
   shear_heat is identically zero in 270 of 271 pairs. radiogenic_source is the
   same kind of field and is already carried through nn_interpolate_elem_fields.

5. rheology.cxx:1123 — with hydraulic diffusion on, pore pressure is booked as heat.
   elastic_effective / elasto_plastic2d add dpp to the diagonal with no strain
   attached, so the inversion reads the whole jump as non-elastic strain. 2-D elastic +
   has_hydraulic_diffusion = yes gives 4.46e−6 W/m³ and +1.04 K where the exact answer is zero. Currently unguarded.

MINOR

6. rheology.cxx:369 — the denom guard can't fire: 3λ+2μ ≡ 3·bulkm, so it needs
   K < 3e−31 Pa. The 1e-30 has no stated provenance. Drop it or make it an assert.
7. rheology.cxx:1135 — has_syy is keyed on the rheology type, but evp/evp_rsf
   only write stressyy when the plastic stress wins; on a maxwell-winning step a fossil
   value enters both traces. Should key on the branch that actually ran.
8. Comments assert something that doesn't hold — the clamp "keeps every genuinely
   dissipative component untouched" (see 1), in both rheology.cxx and the cfg header.
   The header also gives the mean as 2.0e-5 in one place and 1.7e-5 in another.
9. Coverage — the cfg exercises 2-D elasto-plastic only. Nothing covers elastic or
   maxwell (where 2 lives), viscous, 3-D, or is_plane_strain = yes. The elastic
   variant in 2 is the cheapest addition: exact answer, two seconds, no baseline needed.
10. The cfg extends the box, it doesn't compress it — 20.0 → 21.26 km, mean
    ε_xx = +5.5e−2. Not incidental: the rectification in 2 keeps the extensional half, so
    this test's sign is the one that exposes it.

Reference commits, if useful

There are six commits sitting on a branch off yours — about five lines of actual logic. They're there purely as a reference.

One commit per finding, so you can take only the ones you agree with:

commit	finding
clamp on the sum, not per component	1
invert the Hooke's law matching the trace on hand	2, 6
read shear heating from the branch that set the stress	5, 7
register shear heating in the reader	3
remap shear_heat across a remesh	4
re-measure the cfg header, add the elastic zero case	8, 9

Gates they pass, in case it saves you rerunning them: flag-off bitwise identical to this
PR on four cases including 2-D and 3-D remeshing; restart determinism holds; 2-D, 3-D and opt=0 build clean; make openacc=1 builds on tierra4 (nvc++ 25.7, RTX A5000) with the on-device elastic case at 1.49e−17, maxwell matching CPU to 1e−4 relative, and GPU flag-off bitwise identical on the 3-D remeshing case.

Not checked: the geoFLAC cross-check of the trapezoidal-average form, and ep_rsf/evp_rsf.

---

Reviewed with Claude Code (Opus 5). Every number above comes from an A/B of this branch against its base — locally for the CPU runs, tierra4 for the GPU ones — not from
inspection alone.