Follow-up to #71 / #77, which bounded both root vulnerability curves past their grid. The two halves got different qualities of answer, and the reason is worth stating rather than leaving as an inconsistency.
Where it stands after #77
| curve |
limit past the grid |
how #77 handles it |
cumulative integral G |
finite and non-zero, G(inf) = (b/c)*Gamma(1/c) = 3.46578 |
capped at the exact closed form |
conductivity f_r |
exactly zero |
clamped to the last knot's ~0.0113 |
The integral's treatment is exact and free. The conductivity's is an approximation, and a fairly coarse one: at 10 MPa the clamp returns 0.0113 where the true value is 3.77e-06, roughly 3000x too conductive.
Why the obvious fix does not work
"Evaluate exp(-(psi/root_b)^root_c) past the last knot" is exact and costs one pow+exp on a path that is essentially never taken. It also underflows to exactly 0.0 in double precision above 45.98 MPa at the root defaults — and f_r is a divisor (r_R_H = r_R_H_min / f_ri, roots.hpp:672). So the exact form trips the f_ri <= 0.0 guard at roots.hpp:664 and calls util::stop().
Adding a small positive floor to keep the division safe works, but it only relocates the arbitrary constant that the 1% grid end already represents. That is not progress.
The actual question
What should a fully embolised layer do? Not "what number should we substitute for zero conductivity", but "what flux does a layer with no conductivity contribute". The physically right answer is zero, and expressing it means a branch in uptake_impl rather than a value from an accessor — the layer drops out of the sum instead of dividing by a floor.
Urgency: low, and here is why
root_vuln_at is read in exactly one place, roots.hpp:663, inside the equal-potentials branch, which requires psi_soil[i] == T_collar. The collar is bounded by psi_crit (5.870283 MPa at the defaults — which is itself exactly the 5%-conductivity point), so a layer drier than that can never reach this branch. Deep drought goes down the general branch, which reads only the integral, and #77 caps that exactly.
Measured over the golden grid, the conductivity curve is never read past 4.0 MPa. So this is near-dead code, not a live wrong answer — but it is the one place where the choice of the 1% grid end still reaches a returned number, and it is worth closing deliberately rather than by accident.
Not the fix: widening the grid
For the record, since it is the intuitive move. The grid end is psi_max = b * log(1/floor)^(1/c); resolution is a knot count, so step = psi_max/resolution and extending the domain coarsens the knots everywhere. Measured over [0, 5.87], the range the collar solve actually reaches:
| floor |
psi_max |
step |
max spline error in G |
in f_r |
integral still missing |
| 0.01 (current) |
6.892 |
0.0689 |
4.15e-08 |
8.92e-07 |
1.44e-03 |
| 0.001 |
8.017 |
0.0802 |
7.24e-08 |
1.34e-06 |
1.16e-04 |
| 0.0001 |
8.926 |
0.0893 |
1.08e-07 |
1.79e-06 |
9.82e-06 |
Going to 0.01% pays 2.6x the interpolation error where the model actually lives to recover 1.4e-3 of an integral that #77's closed-form cap already supplies exactly. It also raises vulnerability_x_max() = log(1/floor) from 4.605 to 9.21, which bounds the incomplete-gamma series in cumulative_vulnerability_integral_derivatives_at (#81) — more terms, less overflow headroom, and it would need re-measuring.
And it does not fix anything: plant's soil potential ceiling is 1e3 MPa, so the out-of-grid case survives a 0.01% floor by two orders of magnitude regardless. The bound is required either way.
Follow-up to #71 / #77, which bounded both root vulnerability curves past their grid. The two halves got different qualities of answer, and the reason is worth stating rather than leaving as an inconsistency.
Where it stands after #77
GG(inf) = (b/c)*Gamma(1/c)= 3.46578f_rThe integral's treatment is exact and free. The conductivity's is an approximation, and a fairly coarse one: at 10 MPa the clamp returns 0.0113 where the true value is 3.77e-06, roughly 3000x too conductive.
Why the obvious fix does not work
"Evaluate
exp(-(psi/root_b)^root_c)past the last knot" is exact and costs onepow+expon a path that is essentially never taken. It also underflows to exactly 0.0 in double precision above 45.98 MPa at the root defaults — andf_ris a divisor (r_R_H = r_R_H_min / f_ri,roots.hpp:672). So the exact form trips thef_ri <= 0.0guard atroots.hpp:664and callsutil::stop().Adding a small positive floor to keep the division safe works, but it only relocates the arbitrary constant that the 1% grid end already represents. That is not progress.
The actual question
What should a fully embolised layer do? Not "what number should we substitute for zero conductivity", but "what flux does a layer with no conductivity contribute". The physically right answer is zero, and expressing it means a branch in
uptake_implrather than a value from an accessor — the layer drops out of the sum instead of dividing by a floor.Urgency: low, and here is why
root_vuln_atis read in exactly one place,roots.hpp:663, inside the equal-potentials branch, which requirespsi_soil[i] == T_collar. The collar is bounded bypsi_crit(5.870283 MPa at the defaults — which is itself exactly the 5%-conductivity point), so a layer drier than that can never reach this branch. Deep drought goes down the general branch, which reads only the integral, and #77 caps that exactly.Measured over the golden grid, the conductivity curve is never read past 4.0 MPa. So this is near-dead code, not a live wrong answer — but it is the one place where the choice of the 1% grid end still reaches a returned number, and it is worth closing deliberately rather than by accident.
Not the fix: widening the grid
For the record, since it is the intuitive move. The grid end is
psi_max = b * log(1/floor)^(1/c);resolutionis a knot count, sostep = psi_max/resolutionand extending the domain coarsens the knots everywhere. Measured over [0, 5.87], the range the collar solve actually reaches:Going to 0.01% pays 2.6x the interpolation error where the model actually lives to recover 1.4e-3 of an integral that #77's closed-form cap already supplies exactly. It also raises
vulnerability_x_max() = log(1/floor)from 4.605 to 9.21, which bounds the incomplete-gamma series incumulative_vulnerability_integral_derivatives_at(#81) — more terms, less overflow headroom, and it would need re-measuring.And it does not fix anything: plant's soil potential ceiling is 1e3 MPa, so the out-of-grid case survives a 0.01% floor by two orders of magnitude regardless. The bound is required either way.