diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md
index cd197add..a1def6f9 100644
--- a/docs/developer/subsystems/stress-transport.md
+++ b/docs/developer/subsystems/stress-transport.md
@@ -67,7 +67,8 @@ with the current velocity gradient, so over one step it stretches the conformati
 by about $(1 + \Delta t\,\dot\gamma)^2$ before relaxation acts. When
 $\Delta t\,\dot\gamma$ is of order one that update loses the positive-definiteness
 of the conformation $c = \sigma^*/G + I$ in the first step, and no arrangement of
-the split recovers it. On the confined cylinder at Courant one on the far-field
+the split recovers it. That is the default, linear step; the deformation step
+below does not have the limit. On the confined cylinder at Courant one on the far-field
 mesh the wall shear rate is ten times the far-field one: every history lost the
 conformation at the cylinder top in step one, the nodal history then ran away and
 the solve hung, and the integration-point history gave out at Wi 0.6. At a step
@@ -101,6 +102,67 @@ co-rotational part of the tangent grows with $|W||\sigma^*|/G$ and is a solver
 setting) from a solve that has lost its problem (nothing recovers it). Print this
 line every step on a new problem.
 
+## Past the conformation limit: the deformation step and the log-conformation history
+
+Two things lose the conformation at a re-entrant corner or a stagnation point at
+high Weissenberg number, and each has its own remedy.
+
+**The step.** Written in the conformation, the linear upper-convected BDF-1 step is
+
+$$c\,(1 + \Delta t/\lambda) = c^* + \Delta t\,(L c^* + c^* L^T) + (\Delta t/\lambda)\,I,$$
+
+the deformation $F c^* F^T$ with $F = I + \Delta t\,L$ less its second-order term
+$\Delta t^2 L c^* L^T$. Dropping that term is what makes the step indefinite once
+$\Delta t\,|L|$ is of order one. `convected_step="deformation"` keeps it: the
+step is $F c^* F^T + (\Delta t/\lambda) I$, positive-definite for any step and
+any velocity gradient, and still first order. With the exponential integrator the
+stretching of the relaxation target is completed to a product the same way. The
+relaxation itself needs nothing: it is linear in $c$. `max_elastic_timestep`
+returns no limit for this step.
+
+**The store.** A history stores the stress at its own points and hands it back by
+interpolation, projection or a per-cell fit. Near a corner singularity those
+undershoot, and a linear fit extrapolates to the cell edges; the stress they return
+can be indefinite although every stored value is not.
+`stress_history="log_conformation"` stores $\psi = \log c$ and the model reads
+$\sigma^* = G(e^{\psi^*} - I)$, which is a conformation whatever was done to
+$\psi$. It implies the deformation step. Every flavour transports the stored
+tensor by pure advection, so none of them changes; the step is still taken on $c$,
+so neither integrator changes.
+
+```python
+stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(
+    stokes.Unknowns, order=1, objective_rate="upper_convected",
+    stress_history="log_conformation")
+```
+
+The store then holds the dimensionless $\psi$. Read the carried stress through the
+model (`constitutive_model.stress_star`, in pascals in a units model) or
+`stokes.tau` (a projection of it); `DFDt.psi_star` is the raw record. The inflow
+datum is given as a stress and stored through the same encoding;
+`set_initial_history` takes stored values. The health check reports
+`fraction_floored`, where the logarithm's floor ($10^{-12}$, reached only by
+round-off after a positive-definite step) acted on the record.
+
+Measured on the cross-slot (creeping UCM, full resolution, forward history): the
+linear step with the stress stored stalls or hangs at every De from 0.5
+($\lambda U/H$, $H$ the half-width); with the log-conformation history the
+conformation stays above 0.33 everywhere, and a seeded run finds the purely
+elastic pitchfork, the symmetric state stable at De 0.65 and unstable at 0.8,
+onset near 0.71.
+
+Limits: first order only (the second-order schemes combine history levels with
+negative weights); the log-conformation history is 2-D only (closed-form 2x2
+logarithm and exponential); `objective_rate="upper_convected"` only. A
+geodynamic viscoelastic-plastic model with no objective rate carries stresses
+that are not conformations (compression beyond $G$ is legitimate there), so
+neither option applies to it. The deformation term makes the momentum equation
+quadratic in the velocity gradient: a solve at $\Delta t\,|L| \sim 1$ everywhere
+needs a starting velocity, as every step after the first has. The log costs about
+a quarter more error in the stress near a singular corner than storing the stress.
+`SNES_NavierStokes` (the Navier-Stokes solver that reads its history directly as
+a flux) and the multi-material model refuse the log-conformation history.
+
 ## The recommended configuration
 
 Integration-point history, the step set by the wall strain rate
diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py
index 0aa15984..a6bb2a52 100644
--- a/src/underworld3/constitutive_models.py
+++ b/src/underworld3/constitutive_models.py
@@ -1626,6 +1626,41 @@ def _object_viewer(self):
         return
 
 
+#: Smallest conformation eigenvalue the log-conformation history takes the
+#: logarithm of. The positive-definite step never produces a smaller one except
+#: by round-off, so this only bounds psi; the health check counts where it acts.
+_CONFORMATION_FLOOR = 1.0e-12
+
+
+def _sym2_parts(m):
+    """Mean and half-spread of the eigenvalues of a symmetric 2x2 matrix; the
+    spread is lifted by 1e-8 so the closed forms below stay finite at an
+    isotropic point (the lift is far below the solver tolerance)."""
+    a, b, c = m[0, 0], m[1, 1], m[0, 1]
+    mean = (a + b) / 2
+    d = sympy.sqrt(((a - b) / 2) ** 2 + c ** 2 + 1.0e-16)
+    return a, b, c, mean, d
+
+
+def _expm_sym2(m):
+    r"""Closed-form :math:`e^{M}` of a symmetric 2x2 matrix: positive-definite for any M."""
+    a, b, c, mean, d = _sym2_parts(m)
+    ch, sh = sympy.cosh(d), sympy.sinh(d) / d
+    return sympy.exp(mean) * sympy.Matrix([[ch + sh * (a - mean), sh * c],
+                                           [sh * c, ch + sh * (b - mean)]])
+
+
+def _logm_sym2(m):
+    r"""Closed-form :math:`\log M` of a symmetric 2x2 matrix, eigenvalues floored
+    at :data:`_CONFORMATION_FLOOR`."""
+    a, b, c, mean, d = _sym2_parts(m)
+    l1 = sympy.log(sympy.Max(mean + d, _CONFORMATION_FLOOR))
+    l2 = sympy.log(sympy.Max(mean - d, _CONFORMATION_FLOOR))
+    k = (l1 - l2) / (2 * d)
+    return sympy.Matrix([[(l1 + l2) / 2 + k * (a - mean), k * c],
+                         [k * c, (l1 + l2) / 2 + k * (b - mean)]])
+
+
 class ViscoElasticPlasticFlowModel(ViscousFlowModel):
     r"""
     Viscoelastic-plastic flow constitutive model.
@@ -1649,7 +1684,8 @@ class ViscoElasticPlasticFlowModel(ViscousFlowModel):
     """
 
     def __init__(self, unknowns, order=1, integrator: str = "bdf",
-                 material_name: str = None, objective_rate: str = "none"):
+                 material_name: str = None, objective_rate: str = "none",
+                 stress_history: str = "stress", convected_step: str = None):
         """Construct a viscoelastic-plastic flow model.
 
         Parameters
@@ -1684,6 +1720,28 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf",
               See ``docs/developer/design/EXPONENTIAL_VE_INTEGRATOR.md``.
         material_name : str, optional
             Name identifier for this material.
+        objective_rate : {"none", "upper_convected", "jaumann"}, default "none"
+            The objective stress rate the carried stress obeys (see
+            :meth:`_objective_term`).
+        stress_history : {"stress", "log_conformation"}, default "stress"
+            What the stress history stores. ``"log_conformation"`` stores
+            :math:`\psi = \log(\sigma/G + I)` and the model reads
+            :math:`\sigma^* = G(e^{\psi^*} - I)`, a conformation whatever the
+            history's interpolation, projection or fit did to :math:`\psi`.
+            The step is still taken on the conformation, where the relaxation
+            is linear, so neither integrator changes. Upper-convected, first
+            order, 2-D; it implies ``convected_step="deformation"``.
+        convected_step : {"linear", "deformation"}, optional
+            How the upper-convected stretching is taken over one step. In the
+            conformation :math:`c = \sigma/G + I`, ``"linear"`` advances the
+            carried state by :math:`c^* + \Delta t\,(Lc^* + c^*L^T)`, which
+            loses positive-definiteness once :math:`\Delta t\,|L|` is of
+            order one; ``"deformation"`` uses :math:`F c^* F^T` with
+            :math:`F = I + \Delta t\,L`, positive-definite for any step and
+            equally first order (with the exponential integrator the
+            stretching of the relaxation target is completed the same way).
+            Default ``"linear"``, or ``"deformation"`` with the
+            log-conformation history, which requires it.
         """
         if integrator not in ("bdf", "etd"):
             raise ValueError(
@@ -1707,6 +1765,24 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf",
         # the CURRENT velocity gradient, so it is linear in the unknown and
         # first order in time.
         self._objective_rate = objective_rate
+        if stress_history not in ("stress", "log_conformation"):
+            raise ValueError(f"stress_history must be 'stress' or 'log_conformation', got {stress_history!r}")
+        if convected_step is None:
+            convected_step = "deformation" if stress_history == "log_conformation" else "linear"
+        if convected_step not in ("linear", "deformation"):
+            raise ValueError(f"convected_step must be 'linear' or 'deformation', got {convected_step!r}")
+        if stress_history == "log_conformation" and convected_step != "deformation":
+            raise ValueError("the log-conformation history needs convected_step='deformation': the "
+                             "linear step can hand it an indefinite conformation, which has no logarithm")
+        if convected_step == "deformation":
+            if objective_rate != "upper_convected":
+                raise ValueError("convected_step='deformation' and the log-conformation history describe a "
+                                 "conformation tensor: they need objective_rate='upper_convected'")
+            if stress_history == "log_conformation" and unknowns.u.mesh.dim != 2:
+                raise NotImplementedError("the log-conformation history uses the closed-form 2x2 matrix "
+                                          "logarithm and exponential; 3-D is not implemented")
+        self._stress_history = stress_history
+        self._convected_step = convected_step
 
         # Store material_name before creating expressions (needed by create_unique_symbol)
         self._material_name = material_name
@@ -1739,6 +1815,7 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf",
             "Equivalent value of strain rate 2nd invariant (accounting for stress history)",
         )
 
+        self._check_order_supported(order)
         self._order = order
         self._yield_mode = "softmin"  # "min", "harmonic", "smooth", or "softmin"
         self._yield_softness = 0.1  # δ parameter for "softmin" mode
@@ -1938,6 +2015,7 @@ def order(self, value):
         warn if the DFDt was created with a lower order (since it can't be
         changed after creation — the DFDt allocates history buffers at init).
         """
+        self._check_order_supported(value)
         self._order = value
         self._reset()
 
@@ -2087,12 +2165,36 @@ def stress_history_ddt_kwargs(self):
             return {"with_forcing_history": True}
         return {}
 
+    def _check_order_supported(self, order):
+        """The positive-definite step (and so the log-conformation history) is
+        first order: the second-order schemes combine history levels with
+        negative weights, which does not preserve positivity."""
+        if order != 1 and self._convected_step == "deformation":
+            raise NotImplementedError("convected_step='deformation' and the log-conformation history "
+                                      "are first order only")
+
+    def _carried_stress_sym(self, level=0):
+        r"""The carried stress :math:`\sigma^*` of history level ``level``: the
+        stored level, or :math:`G(e^{\psi^*} - I)` for the log-conformation
+        history (the modulus as its expression, so a read of it carries units)."""
+        stored = self.Unknowns.DFDt.psi_star[level].sym
+        if self._stress_history == "stress":
+            return stored
+        return (_expm_sym2(sympy.Matrix(stored)) - sympy.eye(2)) * self.Parameters.shear_modulus
+
+    def encode_history(self, stress):
+        r"""What the history stores for a stress: the stress, or
+        :math:`\log(\sigma/G + I)` for the log-conformation history."""
+        if self._stress_history == "stress":
+            return stress
+        return _logm_sym2(sympy.Matrix(stress) / self.Parameters.shear_modulus + sympy.eye(2))
+
     # The following should have no setters
     @property
     def stress_star(self):
         r"""Previous timestep stress :math:`\boldsymbol{\sigma}^*` from history."""
         if self.Unknowns.DFDt is not None:
-            self._stress_star.sym = self.Unknowns.DFDt.psi_star[0].sym
+            self._stress_star.sym = self._carried_stress_sym(0)
 
         return self._stress_star
 
@@ -2104,7 +2206,7 @@ def stress_2star(self):
 
         if self.Unknowns.DFDt is not None:
             if self.Unknowns.DFDt.order >= 2:
-                self._stress_2star.sym = self.Unknowns.DFDt.psi_star[1].sym
+                self._stress_2star.sym = self._carried_stress_sym(1)
             else:
                 self._stress_2star.sym = sympy.sympify(0)
 
@@ -2150,32 +2252,47 @@ def E_eff(self):
             # expression tree, no separate code path needed.
             alpha = DDt._exp_alpha
             phi = DDt._exp_phi
-            sigma_star = DDt.psi_star[0].sym
+            sigma_star = self._carried_stress_sym(0)
             if DDt.forcing_star is not None:
                 edot_star = DDt.forcing_star.sym
             else:
                 edot_star = sympy.zeros(*E.shape)
             eta_raw = self.Parameters.shear_viscosity_0
-            self._E_eff.sym = (
+            E_eff = (
                 (1 - phi) * E
                 + (alpha / (2 * eta_raw)) * sigma_star
                 + (phi - alpha) * edot_star
                 # objective rate: the flux gains alpha * dt * S(sigma*)
                 + alpha * self.Parameters.dt_elastic * self._objective_term(sigma_star) / (2 * eta_raw)
             )
+            if self._convected_step == "deformation":
+                # The relaxation target is stretched during the step too. The
+                # exponential step carries that as (1-a) I + b (L + L^T) with
+                # b = lam (1 - (1+x) e^-x), x = dt/lam, exact to first order in L;
+                # adding (b^2/(1-a)) L L^T makes it the product
+                # (sqrt(1-a) I + b L/sqrt(1-a))(...)^T, positive semi-definite,
+                # without touching the first-order part. Written in x, not in the
+                # integrator's alpha, which is clamped to one in the elastic limit.
+                L = sympy.Matrix(self.Unknowns.u.sym).jacobian(self.Unknowns.u.mesh.X)
+                lam = eta_raw / self.Parameters.shear_modulus
+                x = self.Parameters.dt_elastic / lam
+                decay = sympy.exp(-x)
+                b = lam * (1 - (1 + x) * decay)
+                E_eff = E_eff + (b ** 2 / (1 - decay)) * L * L.T / (2 * lam)
+            self._E_eff.sym = E_eff
             return self._E_eff
 
         # BDF default
         mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus
         bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3]
         for i in range(DDt.order):
-            E += -bdf_cs[i] * DDt.psi_star[i].sym / (2 * mu_dt)
+            E += -bdf_cs[i] * self._carried_stress_sym(i) / (2 * mu_dt)
         # objective rate on the newest level, explicit and first order: the
         # constitutive law is sigma + (eta/mu)(Dsigma/Dt - S) = 2 eta E, so the
         # flux gains (eta/mu) S(sigma*) and E_eff gains S/(2 mu) with weight
         # ONE whatever the BDF order (the BDF weights belong to the time
         # derivative, not to the source; -c1 = 2 at order 2 would double it).
-        E += self._objective_term(DDt.psi_star[0].sym) / (2 * self.Parameters.shear_modulus)
+        E += self._objective_term(self._carried_stress_sym(0)) / (2 * self.Parameters.shear_modulus)
         self._E_eff.sym = E
         return self._E_eff
 
@@ -2192,6 +2309,13 @@ def _objective_term(self, sigma):
         sigma = sympy.Matrix(sigma)
         L = sympy.Matrix(self.Unknowns.u.sym).jacobian(self.Unknowns.u.mesh.X)
         if self._objective_rate == "upper_convected":
+            if self._convected_step == "deformation":
+                # F c* F^T - c* = dt (L c* + c* L^T) + dt^2 L c* L^T with
+                # c = sigma/G + I: the second-order term, in stress, is
+                # dt (L sigma L^T + G L L^T).
+                dt = self.Parameters.dt_elastic
+                G = self.Parameters.shear_modulus
+                return L * sigma + sigma * L.T + dt * (L * sigma * L.T + G * L * L.T)
             return L * sigma + sigma * L.T
         W = (L - L.T) / 2
         return W * sigma - sigma * W
@@ -2209,7 +2333,15 @@ def _carried_stress(self):
             raise NotImplementedError(f"{type(DDt).__name__} does not expose the carried stress as "
                                       "one tensor per point; the conformation check reads a "
                                       "per-point tensor history (the trace-back flavours)")
-        return DDt.carried_tensors()
+        tau, points = DDt.carried_tensors()
+        if self._stress_history == "log_conformation":
+            from underworld3.systems.ddt import _to_nondim_ndarray
+            G = np.asarray(_to_nondim_ndarray(
+                uw.function.evaluate(self.Parameters.shear_modulus.sym, points))).reshape(-1)
+            w, v = np.linalg.eigh(tau)
+            c = v @ (np.exp(w)[:, :, None] * np.transpose(v, (0, 2, 1)))
+            tau = G[:, None, None] * (c - np.eye(tau.shape[-1])[None])
+        return tau, points
 
     def max_elastic_timestep(self, safety: float = 0.3) -> float:
         r"""The largest step the explicit stretching term tolerates:
@@ -2233,11 +2365,12 @@ def max_elastic_timestep(self, safety: float = 0.3) -> float:
         the strain rate (what ``evaluate`` returns for a gradient) at the points
         of the carried stress. That projection sits a little below the per-cell
         gradient at a wall, which the safety factor covers. Reduced over
-        ranks. ``inf`` when there is no objective rate (nothing stretches) or
-        no flow; a velocity that is not finite raises rather than returning
-        a silent ``nan``.
+        ranks. ``inf`` when there is no objective rate (nothing stretches),
+        with ``convected_step="deformation"`` (the step is positive-definite
+        for any timestep, so this limit does not apply) or with no flow; a
+        velocity that is not finite raises rather than returning a silent ``nan``.
         """
-        if self._objective_rate == "none":
+        if self._objective_rate == "none" or self._convected_step == "deformation":
             return float("inf")
         E = sympy.Matrix(self.Unknowns.E)
         rate = sympy.sqrt(2 * (E.T * E).trace())
@@ -2296,7 +2429,16 @@ def conformation_min_eigenvalue(self):
         n_all = int(comm.allreduce(n_local, op=uw.MPI.SUM))
         n_neg_all = int(comm.allreduce(n_neg, op=uw.MPI.SUM))
         where = tuple(float(x) for x in points[int(lo.argmin())]) if (n_local and local_min == gmin) else None
-        return {"min": gmin, "max": gmax, "fraction_negative": (n_neg_all / n_all) if n_all else 0.0, "where": where}
+        health = {"min": gmin, "max": gmax, "fraction_negative": (n_neg_all / n_all) if n_all else 0.0, "where": where}
+        if self._stress_history == "log_conformation":
+            # A decoded conformation is positive by construction, so the check
+            # that matters is where the logarithm's floor acted on the record.
+            psi, _ = self.Unknowns.DFDt.carried_tensors()
+            n_floor = int((np.linalg.eigvalsh(psi)[:, 0] <= np.log(_CONFORMATION_FLOOR) + 1.0e-6).sum()) \
+                if psi.shape[0] else 0
+            n_floor_all = int(comm.allreduce(n_floor, op=uw.MPI.SUM))
+            health["fraction_floored"] = (n_floor_all / n_all) if n_all else 0.0
+        return health
 
     @property
     def E_eff_inv_II(self):
@@ -2499,13 +2641,11 @@ def stress_projection(self):
         stress = 2 * self.Parameters.ve_effective_viscosity * edot
 
         if self.Unknowns.DFDt is not None:
-            stress_star = self.Unknowns.DFDt.psi_star[0]
-
             if self.is_elastic:
                 # 1st order
                 stress += (
                     self.Parameters.ve_effective_viscosity
-                    * stress_star.sym
+                    * self._carried_stress_sym(0)
                     / (self.Parameters.dt_elastic * self.Parameters.shear_modulus)
                 )
 
@@ -4687,6 +4827,10 @@ def __init__(
         # regression in tests/test_0103_jit_rampable_constants.py.
         # Validate compatibility before initialization
         self._validate_model_compatibility(constitutive_models)
+        if any(getattr(m, "_stress_history", "stress") != "stress" for m in constitutive_models):
+            raise NotImplementedError(
+                "the multi-material model stores its averaged flux as a stress; a constituent "
+                "with stress_history='log_conformation' cannot share that history")
 
         self._material_var = material_swarmVariable
         self._constitutive_models = constitutive_models
diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py
index 7739c992..b21f6a02 100644
--- a/src/underworld3/systems/ddt.py
+++ b/src/underworld3/systems/ddt.py
@@ -1110,6 +1110,17 @@ def inflow_value(self, value):
     #: Whether this flavour compiles :attr:`inflow_value` into its transport.
     applies_inflow_value = False
 
+    #: The owner's map from the carried quantity to what the history stores
+    #: (a solver sets its constitutive model's ``encode_history``); ``None``
+    #: stores the quantity itself.
+    _encode = None
+
+    def _inflow_record(self):
+        """:attr:`inflow_value` in the form the history stores."""
+        if self._encode is None:
+            return self._inflow_value
+        return sympy.Matrix(self._encode(self._inflow_value))
+
     def _nondim_timestep(self, dt):
         """The timestep as a non-dimensional model time (:func:`_as_float`); a
         symbolic timestep passes through unchanged."""
@@ -1123,7 +1134,7 @@ def _write_inflow(self, var, coords, rows):
         only ``rows`` are written). Storage is non-dimensional, so a value
         that evaluates with units is reduced. A flavour that calls this sets
         ``_components`` (its stored columns)."""
-        expr = self._inflow_value
+        expr = self._inflow_record()
         for column, (i, j) in enumerate(self._components):
             vals = uw.function.evaluate(expr[i, j], coords)
             var.data[rows, column] = np.asarray(
@@ -1744,6 +1755,9 @@ def initialise_history(self):
     def set_initial_history(self, values, dt=None):
         r"""Plant history values for BDF restart or analytical IC.
 
+        The values are what the history stores: for a stress history whose
+        model stores the log-conformation, ``log(sigma/G + I)``.
+
         Bypasses the automatic ``effective_order`` ramp so the very
         first solve runs at the full BDF order rather than starting at
         BDF-1. Use this when you have known values at :math:`t` and
@@ -2245,7 +2259,7 @@ class _HistoryTransport(uw.systems.SNES_MultiComponent):
             normal_flow = sum(a[0, i] * self.mesh.Gamma[i] for i in range(self.mesh.dim))
             entering = sympy.Min(normal_flow, 0)
             indices = self._transport_components()
-            incoming = self._inflow_value
+            incoming = self._inflow_record()
             # Sign: `entering` is non-positive, so -entering is |u.n| on the
             # inflow and zero elsewhere, and the term is dissipative in
             # (sigma - sigma_in). With the sign the other way it amplifies:
@@ -3162,6 +3176,9 @@ def initialise_history(self):
     def set_initial_history(self, values, dt=None):
         r"""Plant history values for BDF restart or analytical IC.
 
+        The values are what the history stores: for a stress history whose
+        model stores the log-conformation, ``log(sigma/G + I)``.
+
         Bypasses the automatic ``effective_order`` ramp so the very
         first solve runs at the full BDF order rather than starting
         at BDF-1. Use this when you have known values at :math:`t`
@@ -5667,7 +5684,7 @@ def _fit_arrivals(self, X, values, cell=None, inflow=None):
             # boundary-cell dof, whether or not any of its cells is short
             bcells = np.flatnonzero(np.isin(np.arange(ncell), self._bface_cell))
             filled = np.column_stack([
-                _to_nondim_ndarray(uw.function.evaluate(self._inflow_value[i, j],
+                _to_nondim_ndarray(uw.function.evaluate(self._inflow_record()[i, j],
                                                         dofs[bcells].reshape(-1, mesh.cdim))
                                    ).reshape(-1)
                 for (i, j) in self._components]).reshape(bcells.size, ndof, self.num_components)
diff --git a/src/underworld3/systems/navier_stokes_eulerian.py b/src/underworld3/systems/navier_stokes_eulerian.py
index f15656ea..c45a6dc7 100644
--- a/src/underworld3/systems/navier_stokes_eulerian.py
+++ b/src/underworld3/systems/navier_stokes_eulerian.py
@@ -431,7 +431,10 @@ def _viscous_flux(self):
                 # is rebuilt from the stored velocity, as the new level has it
                 eta_s = getattr(self.constitutive_model.Parameters, "solvent_viscosity", 0)
                 solvent = 2 * eta_s * sympy.Matrix(self.mesh.vector.strain_tensor(u_k))
-                total = total + w * (sympy.Matrix(stress_history.psi_star[level].sym) + solvent)
+                carried = self.constitutive_model._carried_stress_sym(level) \
+                    if hasattr(self.constitutive_model, "_carried_stress_sym") \
+                    else stress_history.psi_star[level].sym
+                total = total + w * (sympy.Matrix(carried) + solvent)
             else:
                 raise ValueError(
                     f"the time scheme weights the flux at level {level + 1}, but the "
diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index 8021763d..1c6614f3 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -85,13 +85,20 @@ def expression(*args, **kwargs):
     return public_expression(*args, _unique_name_generation=True, **kwargs)
 
 
+def _history_record(constitutive_model):
+    """What a stress history stores: the model's memory part of the flux when
+    it separates one out (a solvent viscosity is rebuilt each step), else the
+    whole flux, in the model's encoding of it (:meth:`encode_history`)."""
+    flux = getattr(constitutive_model, "history_flux", None)
+    if flux is None:
+        flux = constitutive_model.flux
+    encode = getattr(constitutive_model, "encode_history", None)
+    return flux if encode is None else encode(flux)
+
+
 def _history_psi_fn(constitutive_model):
-    """The flux a stress history carries: the model's memory part when it
-    separates one out (a solvent viscosity is rebuilt each step), else the
-    whole flux. Transposed to the history's row layout."""
-    if hasattr(constitutive_model, "history_flux"):
-        return constitutive_model.history_flux.T
-    return constitutive_model.flux.T
+    """:func:`_history_record`, transposed to the history's row layout."""
+    return _history_record(constitutive_model).T
 
 
 def _as_scalar(value):
@@ -1732,9 +1739,7 @@ def _stress_history_post_solve(self, timestep, verbose=False, evalf=False):
         # levels. A particle-carried history does that itself in its post-solve,
         # by evaluating the new stress at its own particles.
         if not self.DFDt.commits_flux_in_post_solve:
-            self.DFDt.commit_flux_to_history(
-                getattr(self.constitutive_model, "history_flux", self.constitutive_model.flux),
-                verbose=verbose)
+            self.DFDt.commit_flux_to_history(_history_record(self.constitutive_model), verbose=verbose)
 
         self.DFDt.update_post_solve(timestep, verbose=verbose, evalf=evalf)
         self._devss_refresh(verbose=verbose)
@@ -1780,9 +1785,10 @@ def _create_stress_history_ddt(self, order=2):
             bcs=None,
             order=order,
             smoothing=0.0001,
-            # the history carries a stress; the flux handed over at construction
-            # is a zero placeholder, so it cannot say so itself
-            units=uw.units.Pa,
+            # the history carries a stress (the flux handed over at construction
+            # is a zero placeholder, so it cannot say so itself), or the
+            # dimensionless log-conformation of one
+            units=uw.units.Pa if getattr(cm, "_stress_history", "stress") == "stress" else None,
         )
         if self.stress_transport == "integration_point":
             unsupported = set(ddt_kwargs) - {"with_forcing_history"}
@@ -1871,6 +1877,9 @@ def _create_stress_history_ddt(self, order=2):
         # flux→psi_star[0] becomes implicit in psi_star[0] and Min-mode at
         # yield admits the wrong fixed point under timestep change.
         self.Unknowns.DFDt.enable_source_snapshot()
+        # the history stores the model's encoding of a stress; an inflow datum
+        # is given as a stress and stored through the same encoding
+        self.Unknowns.DFDt._encode = getattr(cm, "encode_history", None)
 
     @timing.routine_timer_decorator
     @memprobe.instrument("Stokes.solve")
@@ -2078,11 +2087,24 @@ def tau(self):
         r"""Deviatoric stress from the most recent solve.
 
         When stress history is active (VEP), returns ``psi_star[0]`` which
-        contains the actual projected stress. Otherwise falls through to the
-        base class lazy projection.
+        contains the actual projected stress; a log-conformation history is
+        decoded into a projected stress variable. Otherwise falls through to
+        the base class lazy projection.
         """
         if self.Unknowns.DFDt is not None:
-            return self.DFDt.psi_star[0]
+            cm = self.constitutive_model
+            if getattr(cm, "_stress_history", "stress") == "stress":
+                return self.DFDt.psi_star[0]
+            if getattr(self, "_tau_decoded", None) is None:
+                self._tau_decoded = uw.discretisation.MeshVariable(
+                    "tau_decoded", self.mesh, (self.mesh.dim, self.mesh.dim),
+                    vtype=uw.VarType.SYM_TENSOR, degree=self.DFDt.psi_star[0].degree,
+                    continuous=True, units=uw.units.Pa if uw.get_default_model().has_units() else None)
+                self._tau_decode = uw.systems.Tensor_Projection(self.mesh, self._tau_decoded)
+                self._tau_decode.smoothing = 0.0
+            self._tau_decode.uw_function = cm._carried_stress_sym(0)
+            self._tau_decode.solve()
+            return self._tau_decoded
         return super().tau
 
     # =========================================================================
@@ -5502,6 +5524,10 @@ def F1(self):
 
         if DFDt is not None:
             # We can flag to only do this if the constitutive model has been updated
+            if getattr(self._constitutive_model, "_stress_history", "stress") != "stress":
+                raise NotImplementedError(
+                    "SNES_NavierStokes reads its history as a flux (Adams-Moulton), so it cannot "
+                    "decode a log-conformation history; use uw.systems.NavierStokes")
             DFDt.psi_fn = getattr(self._constitutive_model, 'history_flux', self._constitutive_model.flux).T
 
             F1 = expression(
diff --git a/tests/test_1064_stress_history_units.py b/tests/test_1064_stress_history_units.py
index 3a8f2db1..6d72401f 100644
--- a/tests/test_1064_stress_history_units.py
+++ b/tests/test_1064_stress_history_units.py
@@ -28,7 +28,7 @@
     + [(t, 2) for t in ("semi_lagrangian", "integration_point", "eulerian")]
 
 
-def _shear_box(transport, order, with_units):
+def _shear_box(transport, order, with_units, **model_options):
     uw.reset_default_model()
     if with_units:
         uw.get_default_model().set_reference_quantities(
@@ -40,7 +40,8 @@ def _shear_box(transport, order, with_units):
     p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1, units="Pa" if with_units else None)
     stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
     stokes.stress_transport = transport
-    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(
+        stokes.Unknowns, order=order, **model_options)
     parameters = stokes.constitutive_model.Parameters
     if with_units:
         parameters.shear_viscosity_0 = uw.quantity(1.0e21, "Pa*s")
@@ -83,3 +84,23 @@ def test_units_model_history_is_the_nondimensional_history(transport, order):
     value_pa = float(uw.units.Quantity(float(np.asarray(carried).reshape(-1)[0]), carried.units).to("Pa").magnitude)
     expected_pa = plain_xy * STRESS_SCALE_PA
     assert abs(value_pa - expected_pa) < 1.0e-6 * abs(expected_pa), (transport, order, value_pa, expected_pa)
+
+
+@pytest.mark.parametrize("transport", ["semi_lagrangian", "forward"])
+def test_units_model_log_conformation_history(transport):
+    """The log-conformation record is dimensionless: the same in both runs; the
+    stress the model reads from it, G (exp(psi) - I), is in Pa."""
+    options = dict(objective_rate="upper_convected", stress_history="log_conformation")
+    plain = _shear_box(transport, 1, False, **options)
+    plain_store = np.array(plain.DFDt.psi_star[0].data)
+    plain_xy = float(np.asarray(uw.function.evaluate(
+        plain.constitutive_model.stress_star.sym[0, 1], ORIGIN)).reshape(-1)[0])
+
+    stokes = _shear_box(transport, 1, True, **options)
+    store = np.asarray(stokes.DFDt.psi_star[0].data)
+    assert np.abs(store - plain_store).max() < 1.0e-6 * np.abs(plain_store).max(), transport
+
+    carried = uw.function.evaluate(stokes.constitutive_model.stress_star.sym[0, 1], ORIGIN)
+    value_pa = float(uw.units.Quantity(float(np.asarray(carried).reshape(-1)[0]), carried.units).to("Pa").magnitude)
+    expected_pa = plain_xy * STRESS_SCALE_PA
+    assert abs(value_pa - expected_pa) < 1.0e-6 * abs(expected_pa), (transport, value_pa, expected_pa)
diff --git a/tests/test_1065_log_conformation.py b/tests/test_1065_log_conformation.py
new file mode 100644
index 00000000..ba14f502
--- /dev/null
+++ b/tests/test_1065_log_conformation.py
@@ -0,0 +1,148 @@
+"""The positive-definite convected step and the log-conformation stress history.
+
+In the conformation c = sigma/G + I, with F = I + dt L, x = dt/lam, the
+deformation step is
+    BDF-1:  c (1 + x) = F c* F^T + x I
+    ETD-1:  c = a F c* F^T + ((1-a) I + b L)((1-a) I + b L^T)/(1-a),
+            a = exp(-x), b = lam (1 - (1 + x) exp(-x))
+and the linear step drops the dt^2 L c* L^T of F c* F^T.
+
+One step of uniform planar extension from rest, u = (e x, -e y), UCM with
+eta = G = 1, dt e = 1: the linear BDF step leaves c_yy = -0.818, not a
+conformation; the deformation step gives the closed forms above, whether the
+history stores the stress or its log-conformation.
+
+The start-up of simple shear: every flavour, both integrators, the
+log-conformation history against the same recurrence, step by step, in both
+the shear stress and the first normal-stress difference (every deformation term
+lives in sigma_xx here).
+"""
+
+import numpy as np
+import pytest
+import sympy
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_2, pytest.mark.tier_b]
+
+DT, RATE = 0.1, 10.0            # dt * rate = 1, lam = 1
+
+
+def _deformation_step(c_star, L, dt, lam, integrator):
+    I = np.eye(2)
+    F = I + dt * L
+    x = dt / lam
+    if integrator == "bdf":
+        return (F @ c_star @ F.T + x * I) / (1 + x)
+    a = np.exp(-x)
+    b = lam * (1 - (1 + x) * np.exp(-x))
+    return a * F @ c_star @ F.T + ((1 - a) * I + b * L) @ ((1 - a) * I + b * L.T) / (1 - a)
+
+
+def _stored_conformation(stokes, stress_history):
+    """Smallest and largest conformation eigenvalue over the history's stored
+    nodes (2-D symmetric storage: xx, yy, xy), decoded from the log if need be."""
+    stored = np.asarray(stokes.DFDt.psi_star[0].data).reshape(-1, 3)
+    m = np.empty((stored.shape[0], 2, 2))
+    m[:, 0, 0], m[:, 1, 1], m[:, 0, 1], m[:, 1, 0] = stored[:, 0], stored[:, 1], stored[:, 2], stored[:, 2]
+    if stress_history == "log_conformation":
+        w, v = np.linalg.eigh(m)
+        c = v @ (np.exp(w)[:, :, None] * np.transpose(v, (0, 2, 1)))
+    else:
+        c = m + np.eye(2)                          # G = 1
+    ev = np.linalg.eigvalsh(c)
+    return ev[:, 0].min(), ev[:, 1].max()
+
+
+def _extension_step(transport, convected_step, stress_history, integrator="bdf"):
+    uw.reset_default_model()
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4), minCoords=(-1, -1), maxCoords=(1, 1))
+    x, y = mesh.X
+    tag = f"{transport[:3]}{convected_step[:3]}{stress_history[:3]}{integrator}"
+    v = uw.discretisation.MeshVariable(f"U_{tag}", mesh, 2, degree=2)
+    p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1)
+    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    stokes.stress_transport = transport
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(
+        stokes.Unknowns, order=1, integrator=integrator, objective_rate="upper_convected",
+        convected_step=convected_step, stress_history=stress_history)
+    parameters = stokes.constitutive_model.Parameters
+    parameters.shear_viscosity_0 = 1.0
+    parameters.shear_modulus = 1.0
+    parameters.dt_elastic = DT
+    extension = (RATE * x, -RATE * y)
+    for boundary in ("Left", "Top", "Bottom"):     # the right side is free: sigma_xx - p = 0 there
+        stokes.add_dirichlet_bc(extension, boundary)
+    stokes.tolerance = 1.0e-8
+    stokes.petsc_options["snes_type"] = "newtonls"
+    # The deformation step makes the momentum equation quadratic in grad u; at
+    # dt |L| = 1 Newton does not converge from rest, so start from the answer,
+    # and from rest in the stress (set_initial_history takes stored values; zero
+    # is rest in both representations).
+    X = np.asarray(v.coords)
+    v.data[:, 0], v.data[:, 1] = RATE * X[:, 0], -RATE * X[:, 1]
+    stokes.DFDt.set_initial_history([0.0])
+    stokes.solve(timestep=DT, zero_init_guess=False)
+    return _stored_conformation(stokes, stress_history)
+
+
+@pytest.mark.parametrize("integrator", ["bdf", "etd"])
+@pytest.mark.parametrize("transport, stress_history", [
+    ("semi_lagrangian", "stress"), ("semi_lagrangian", "log_conformation"),
+    ("eulerian", "stress"), ("eulerian", "log_conformation")])
+def test_deformation_step_is_the_closed_form(transport, stress_history, integrator):
+    c = _deformation_step(np.eye(2), np.diag([RATE, -RATE]), DT, 1.0, integrator)
+    c_min, c_max = _extension_step(transport, "deformation", stress_history, integrator)
+    assert abs(c_min - c[1, 1]) < 1.0e-6, (c_min, c[1, 1])
+    assert abs(c_max - c[0, 0]) < 1.0e-6, (c_max, c[0, 0])
+
+
+def test_linear_step_loses_the_conformation():
+    """The defect the deformation step removes, measured against its closed form."""
+    c_min, _ = _extension_step("semi_lagrangian", "linear", "stress")
+    assert abs(c_min - (1 - 2 * DT * RATE + DT) / (1 + DT)) < 1.0e-6, c_min
+
+
+def _shear_startup(transport, integrator, steps=10, dt=0.1):
+    uw.reset_default_model()
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(16, 8), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5))
+    tag = f"s{transport[:3]}{integrator}"
+    v = uw.discretisation.MeshVariable(f"U_{tag}", mesh, 2, degree=2)
+    p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1)
+    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    stokes.stress_transport = transport
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(
+        stokes.Unknowns, order=1, integrator=integrator, objective_rate="upper_convected",
+        stress_history="log_conformation")
+    parameters = stokes.constitutive_model.Parameters
+    parameters.shear_viscosity_0 = 1.0
+    parameters.shear_modulus = 1.0
+    parameters.dt_elastic = dt
+    stokes.add_dirichlet_bc((0.5, 0.0), "Top")
+    stokes.add_dirichlet_bc((-0.5, 0.0), "Bottom")
+    stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left")
+    stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right")
+    stokes.tolerance = 1.0e-8
+    stokes.petsc_options["snes_type"] = "newtonls"
+    stokes.petsc_options["ksp_type"] = "fgmres"
+    for _ in range(steps):
+        stokes.solve(timestep=dt, zero_init_guess=False)
+    sigma = stokes.constitutive_model.stress_star.sym
+    origin = np.array([[0.0, 0.0]])
+    read = lambda e: float(np.asarray(uw.function.evaluate(e, origin)).reshape(-1)[0])
+    c = np.eye(2)
+    L = np.array([[0.0, 1.0], [0.0, 0.0]])              # shear rate one
+    for _ in range(steps):
+        c = _deformation_step(c, L, dt, 1.0, integrator)
+    return (read(sigma[0, 1]), c[0, 1]), (read(sigma[0, 0] - sigma[1, 1]), c[0, 0] - c[1, 1])
+
+
+@pytest.mark.parametrize("integrator", ["bdf", "etd"])
+@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point", "forward", "eulerian"])
+def test_log_conformation_history_follows_the_recurrence(transport, integrator):
+    """Uniform stress, so transport is a no-op: this checks the encoding, the
+    decoding and the step together, against the discrete recurrence."""
+    (xy, xy_exact), (n1, n1_exact) = _shear_startup(transport, integrator)
+    assert abs(xy - xy_exact) < 1.0e-5 * abs(xy_exact), (transport, integrator, xy, xy_exact)
+    assert abs(n1 - n1_exact) < 1.0e-5 * abs(n1_exact), (transport, integrator, n1, n1_exact)