Size or age: a choice only a characteristics solver gets to make
plant @develop 53d0a20
PR #590 adds Control$node_density_in_birth_date (default FALSE). With it on, the SCM carries its size distribution as a density per unit birth date rather than per unit height.
That sounds like a bookkeeping preference. It is closer to a change of footing. The ∂g/∂h term in the transport equation stops being computed at all, a boundary special case disappears, a class of numerical blow-up becomes structurally impossible, the solver gets about twice as fast, and — for TF24 — one two-species run changes from competitive exclusion to coexistence.
This is Andrew O’Reilly-Nugent’s discovery. Chasing why the solver behaved oddly once TF24 gained a carbohydrate store, he worked out that the compression term the SCM computes stops being the quantity the transport equation asks for as soon as growth depends on any state other than size — then identified the birth-date coordinate as the fix, implemented it behind a flag, and validated it against an independent solver. The write-up is the primary source for everything below; this post is a gloss on it.
What follows is about why the choice exists, why it is a choice we get to make and a fixed-grid solver does not, and where that puts us relative to De Roos (1988).
Two coordinates, one population
The size-structured PDE is a conservation law. Expanding its flux:
\[ \frac{\partial N}{\partial a} + g\frac{\partial N}{\partial H} + N\frac{\partial g}{\partial H} = -d\,N \]
The first two terms move density along the height axis. The third is a compression term: where growth slows with size, neighbouring trajectories converge and density piles up. It is the term that requires knowing how the growth rate varies with size, and it is the whole difficulty.
Now change variables. Let \(\tau\) be a cohort’s introduction time — constant along its characteristic, by definition. Write \(J = \partial H/\partial\tau\) for the Jacobian of the map, and define \(\nu = N\,|J|\): abundance per unit birth date rather than per unit height. The first two terms are exactly the material derivative and collapse into \(\partial/\partial a|_\tau\); the third is precisely the rate of change of \(\log|J|\) and cancels against the Jacobian in \(\nu\). What is left is
\[ \frac{\partial \nu}{\partial a}\bigg|_\tau = -d\,\nu \]
The compression term has not been approximated or moved. It has been absorbed into the coordinate change.
Three ways to say the same thing:
| coordinate | transport velocity | compression term |
|---|---|---|
| size \(H\) | \(g(x, E_a)\) — state- and environment-dependent | \(\partial g/\partial H \neq 0\) |
| birth date \(\tau\) | \(0\) (a label; constant along characteristics) | \(0\) |
| age \(\alpha = a - \tau\) | \(1\), exactly and always | \(0\) |
The third row is the classical McKendrick–Von Foerster age-structured equation (McKendrick, 1926; Foerster, 1959). Age-structured models do not lack a compression term because they are simpler models — they lack it because their transport velocity is a constant, and the divergence of a constant field is zero. All the nonlinearity has moved out of the transport operator and into the rate functions, where nothing has to be differentiated.
And the biology stays size-structured throughout. Mortality and fecundity are still functions of the full state and the environment; the state is still obtained by integrating the size ODEs, at unchanged cost; the competition kernel is still \(A_{\rm l}(H)\,Q(z,H)\) evaluated at \(H(\tau,a)\). Only the coordinate in which abundance is booked changes.
Why the choice depends on how you solve
Here is the part that is easy to miss. The coordinate is not free to choose; it is downstream of the discretisation.
An Eulerian solver puts a fixed grid in state space and lets the population flow through it. Its unknowns are densities in fixed size bins. Because the grid does not move with the population, the flux divergence is unavoidable — the whole content of the scheme is bookkeeping how much mass crosses each bin boundary. That is a perfectly good engineering choice: bounded work per step, no schedule to tune, trivially parallel, no cohorts accumulating over time. It is what finite-volume schemes and Integral Projection Model transition matrices do. But such a solver has no birth-date coordinate available to it, because its grid is the coordinate. Asking it to count by birth date would mean discarding the grid.
A Lagrangian solver moves its unknowns with the population. Each unknown is a cohort travelling along a characteristic. Here the coordinate is a label, and the obvious label is the birth date. The transport term vanishes identically, not approximately.
plant is Lagrangian and always has been — the SCM is a characteristics method, following De Roos, Taljapurkar, & Caswell (1997) and Angulo & López-Marcos (2004). Which means the size-structured PDE was always the derivation, never the implementation. The ∂g/∂h term was the last vestige of the Eulerian formulation the model was derived from, sitting inside a scheme that had otherwise left it behind. #590 does not bolt a special case onto the solver; it removes a seam.
The general statement: a characteristics solver can choose its coordinate because its unknowns are labelled rather than located. That freedom is not available to a grid, and it is the reason this switch is even expressible.
What Andrew found
The prompt was NSC storage. TF24 carries a reserve pool that gates growth through \(G(r)\), \(r = S/S_{\max}\) — so height growth is no longer a function of size and environment alone.
The compression term the PDE needs is the total derivative of \(g\) along a cohort’s own trajectory:
\[ \frac{\rm d}{{\rm d}a}\log N = -d - \sum_k \frac{\partial g_H}{\partial x_k}\,\frac{{\rm d}x_k}{{\rm d}H}\bigg|_{\text{cohort curve}} \]
Node::growth_rate_gradient computed only the \(k = H\) term, by perturbing height and re-running the physiology. For FF16 and K93 those coincide exactly, because their non-height states are pure accumulators that never feed back into the rates — which is why this was correct for a decade. For TF24 they do not.
Andrew’s sharpest observation is about the mechanism rather than the size of the discrepancy. The probe moves height while holding absolute carbon fixed, so it shifts the reserve fraction; a real cohort grows with \(r\) roughly constant. The finite difference is therefore accurate about a physiological state the plant never occupies. In his phrasing, the estimate “stops computing the required quantity — not approximately, but at all — as soon as growth depends on any state other than size.” That is why the error is large rather than a few percent, and he derives the omitted term for TF24 in closed form.
What we gain
The boundary special case disappears. The height coordinate’s boundary condition divides by the growth rate, with a g > 0 branch and a requirement that \(S_{\rm G}/g \to 0\) as \(g \to 0\) (Falster, Brännström, Dieckmann, & Westoby, 2011). In birth date the birth density is birth_rate · pr_estab — finite always. The zero-density cohorts we had catalogued as expected SCM behaviour are ordinary cohorts occupying a compressed band of height space.
Density runaway becomes structurally impossible. plant carries a hard log_density ceiling of 50 because \(N \sim 1/|\partial h/\partial\tau|\) can diverge. Andrew measures the maximum reached across a run as \(+2.4168\) in height against \(-0.0015\) in birth date — the latter being the predicted bound, since \(\nu\) can only decrease. It is worth being precise about what that guard was catching: in a stationary patch the height-coordinate density genuinely diverges as \(N \sim (h_{\max}-h)^{d/c-1}\) when \(d < c\) (background mortality below the rate at which growth shuts off), and log_density_dt → c - d > 0, so it grows without bound. Any long run with a low-mortality strategy is guaranteed to trip that ceiling — and it is not an error, it is the correct solution viewed in a coordinate where it is singular.
Age is monotone by construction; size is monotone only by assumption. The size coordinate is valid only while \(\tau \mapsto H\) stays invertible, i.e. while cohorts never converge in height. Reserve-gated growth can violate that — under full-amplitude seasonal stress Andrew counts 66 crossings among 10,011 interior pairs, against zero under a constant environment. Age needs no such assumption: \({\rm d}\alpha/{\rm d}a = 1 > 0\) always, so cohorts can never overtake one another in age, and their ordering is fixed at birth. The fragility was inherited from the coordinate, not from the biology.
Accuracy, at the same cost. For FF16 and K93 the two coordinates must converge, and do — but the gap is almost entirely the height coordinate’s error. Over successive halvings of the node spacing, FF16’s offspring production:
| schedule | height | birth date | gap |
|---|---|---|---|
| default | 18.5071 | 18.6956 | 1.0e-2 |
| ×2 | 18.6211 | 18.6989 | 4.2e-3 |
| ×4 | 18.6777 | 18.6996 | 1.2e-3 |
| ×8 | 18.6941 | 18.6996 | 2.9e-4 |
The birth-date answer barely moves; the height answer climbs toward it. At the default schedule birth date is within 2e-4 of converged and height is 1.0e-2 away — roughly 50× more accurate for the same number of nodes. K93 behaves the same way (2.6e-3 → 3.5e-5), and so do the two-species FF16 and three-species K93 cases, per species, including a competitively marginal species ~90× below the dominant one.
Speed. growth_rate_gradient is not called, which removes a leaf solve per cohort per Runge–Kutta stage: Andrew measures 1.91–2.39×.
It is indifferent to state dimension. \(\partial\log\nu/\partial a = -d\) is unchanged for one state or twenty, because the label space has the dimension of the birth manifold, not of the state space. With acclimation and further NSC states coming to TF24, this is arguably the most valuable property of the lot.
And it can change the answer. On the two-species TF24 case, the height coordinate excludes the slower species (offspring production 1.4e-4 against 503) while birth date has the two coexisting (1004 against 2707), with refinement moving neither toward the other. The competitive-exclusion conclusion recorded for reserve-gated growth is therefore coordinate-dependent and needs re-deriving.
A partial return to the Escalator Boxcar Train
There is a history here. plant’s earliest incarnation — before this repository — was an Escalator Boxcar Train (De Roos, 1988; De Roos, Diekmann, & Metz, 1992). Daniel’s recollection of why it did not stay one is the instructive part: he still found himself needing to track the boundary cohort’s sizes, and ended up adding densities at the boundary, at which point the cohort moments looked like machinery being carried for no return. The model moved to densities along characteristics and stayed there for the next decade and a half.
So #590 is a return journey. But it is worth being exact about which leg of it we have travelled, because the two schemes differ in something more fundamental than their bookkeeping — and we have not become an EBT.
What is shared is the conservation principle. The characteristics method as conventionally presented (De Roos et al., 1997; Angulo & López-Marcos, 2004) — the form plant has implemented since — carries the density \(N\) along each characteristic, compression term and all. The EBT instead conserves each cohort’s abundance, precisely so that the growth rate never has to be differentiated. \(\partial\nu/\partial a = -d\,\nu\) contains no derivative of \(g\), and that is the EBT’s insight.
What is not shared is the representation, and that is not a detail. The EBT does not represent a density at all. It discretises the population as a sum of point masses: an internal cohort is a number \(N_i\) sitting at a representative state \(X_i\), and a population integral is the exact sum \(\sum_i N_i\,f(X_i)\). There is no quadrature rule, and no notion of the distance between neighbouring cohorts — spacing does not enter.
plant remains a density method. Each node carries \(\nu\), a density per unit birth date, and population integrals are trapezium quadratures \(\int f(x(\tau))\,\nu(\tau)\,{\rm d}\tau\) over the node grid. The spacing between adjacent nodes is load-bearing: it is the \({\rm d}\tau\). So \(\nu\) is not the EBT’s cohort abundance — \(\nu\,{\rm d}\tau\) is what corresponds to \(N_i\) — and the adaptive schedule that sets our spacing has no counterpart in an EBT at all.
The two therefore have different error structures, which is the clearest evidence they are not the same method. The EBT’s error comes from representing a cohort with real internal spread by a single point, which is what the first moment on its boundary cohort exists to mitigate. Ours is quadrature error in \(\tau\), controlled by refining the schedule. That is why the convergence tables above are the right way to validate this change — and why they would be meaningless for an EBT, which has no spacing to refine.
The boundary treatment differs for that reason rather than for a shared one. The EBT’s boundary cohort accumulates recruits arriving continuously across an interval, so a point mass cannot represent it: it carries a first moment until it is “internalised” into the interior. plant’s boundary node needs no moment — it carries a density value, birth_rate · pr_estab, and the exact Jacobian \(-g(H_0)\), with the spread represented by the interval to its neighbour. Both schemes need special handling at the boundary; they need different handling, and the original migration ran along precisely this seam.
Net: we have adopted the EBT’s conservation law while keeping the density-on-a-moving-grid discretisation of De Roos et al. (1997). That is one leg of a two-leg journey, and there is no plan to walk the second — the moments were dropped for a reason, and nothing here has changed it.
What it costs
The size distribution stops being a primitive of the model and becomes a derived output. Since the coordinate ought to be an implementation detail, plant converts at the reporting boundary: Species$log_densities, Patch$state, and therefore everything downstream in tidy_outputs.R and the plots, always receive a density in height, whichever coordinate was carried. The conversion is \(N = \nu/|{\rm d}h/{\rm d}\tau|\), with the Jacobian from central differences of adjacent node heights against their birth dates, and exactly \(-g(H_0)\) at the boundary node, which is recorded at birth. The typical node agrees with an otherwise identical height-coordinate run to 1.7e-3 in log space — 0.17% in density.
The honest caveat: the worst node does not improve with refinement. \(|{\rm d}h/{\rm d}\tau|\) is a ratio of two differences that both shrink, so cancellation sets a floor, and refining adds nodes in the near-empty tail where that floor is worst. Aggregate and plotting use is sound; individual node densities far out in the tail are not. Fixing it properly means evolving \(\partial h/\partial\tau\) as its own state, whose rate is the compression term this change exists to avoid — so it would be an opt-in extra carried purely for reporting fidelity.
Where the equivalence stops
Worth stating plainly, because it closes the loop with the grid-versus- characteristics theme. The coordinate change needs deterministic characteristics from a finite-dimensional birth manifold. plant satisfies that: birth state is a deterministic function of introduction time, so the label space is one-dimensional however many state variables a strategy carries.
Structured population models in general do not. If individual growth is stochastic, the governing equation is advection–diffusion rather than pure transport, characteristics do not exist, and age no longer determines state — individuals of the same age have a distribution of sizes. That is exactly the case for Integral Projection Models, whose kernel is a growth distribution. And a model with no characteristics has no choice but to be Eulerian in size.
So the ordering is: the biology decides whether characteristics exist; that decides whether a Lagrangian scheme is available; and only a Lagrangian scheme gets to choose its coordinate. plant sits at the end of that chain, which is why this switch is ours to make.
The flag is off by default and the solver is bit-identical with it off. Before it becomes the default: the competitive-exclusion result for reserve-gated growth needs re-deriving in the new coordinate, and the thermal-damage work of issue #566 was evaluated on the old one. The theory pages (the size-structured PDE, the characteristic method) still present the height-coordinate derivation, including the g > 0 boundary caveat that this change retires.