Mathematical analysis of the pressure-solver discontinuity (issue #77)#1
Open
maresb wants to merge 11 commits into
Open
Mathematical analysis of the pressure-solver discontinuity (issue #77)#1maresb wants to merge 11 commits into
maresb wants to merge 11 commits into
Conversation
…ford#77) Merged synthesis of two independent, mutually-audited analyses (discontinuity_analysis/ and reports/issue77/) of the period-2 convergence failure. The fixed-point map is discontinuous because cape() integrates signed buoyancy to the topmost positive grid level and clamps at zero; for the PR dgilford#77 sounding (Hurricane Francine's pre-landfall environment under a -1.6 K counterfactual SST) the only diagonal crossing falls inside a 0.84 hPa jump at x_d ~ 950.829 hPa, so no fixed point exists. Smooth toy models separate genuine variational discontinuities (argmax jumps with continuous value; barrier transitions, which are free-convection rather than PI physics) from this artificial one; refuted alternatives (quadrature, grid resolution, interpolation smoothness, floating point) are dispatched in one subsection. Repair: CAPE = max_q W(q) (value function, matching forced eyewall ascent semantics) plus residual-checked bracketed root finding with an explicit non-hypercane endpoint hypothesis. Exact corrected roots 950.3136889226 hPa (dgilford#77) and 960.7687 hPa (1995), verified by both implementations; 7,093-column benchmark: zero flag changes, 98.1% of VMAX within 0.1 m/s. The compiled PDF is gitignored (build with latexmk -pdf). Generated-with: Claude Fable 5
The analysis stands on its own as the definitive writeup: drop the synthesis framing, the companion citations and appendix, and keep a neutral note that all headline quantities were reproduced by a second independently written implementation. Generated-with: Claude Fable 5
The previous 'generic column' was the sample's single most extreme case (VMAX 97 m/s), whose saturated parcel is still +5 K buoyant at the last retained level, so the plotted 'LNB' was the ptop-truncation clip at 70 hPa and buoyancy never changed sign. Replace it with a representative column (January, 10N 108W, VMAX 72.5): buoyancy crosses zero between 100 and 70 hPa, interior interpolated LNB at 91.2 hPa, fixed point 944.6 hPa. Add a remark quantifying the clipping (23/7093 columns, all VMAX > 86, with competing TO-cold / CAPE-low biases, unquantifiable from the shipped sample whose data stops at 50 hPa) and extend implementation phase C to flag clipped LNBs. Generated-with: Claude Fable 5
At the point where E_top / E_max / E_reach are defined, show concrete profiles spanning every way they can relate: (a) single hump, all agree; (b) dominated upper cell without barrier, E_top undercounts; (c) marginal bump behind a barrier (the issue-77 shape), E_top clamps to 0; (d) dominant upper cell behind a barrier, E_reach stalls. Each panel marks the three terminal points and values. The toy-model sweeps are linked to the archetypes they animate. Generated-with: Claude Fable 5
- Author line: break manually (maketitle sets authors in an unbreakable
tabular, causing a 109pt overfull).
- Convention 2/3: breakable set notation; replace the over-wide prose
display for E_top with prose plus the compact formula.
- Eq. (2): stack g and CAT instead of side-by-side.
- Status table: ragged-right p-columns (array package) to kill
underfulls; final-interpretation causal chain restacked as one
five-row alignment.
- Appendix table: footnotesize plus row splits to fit the text width.
- Mute the spurious 'Command \showhyphens has changed' kernel warning
(reproducible with bare \usepackage{microtype} on this TeX Live 2025
install) around microtype's begindocument check only.
- Add \emergencystretch for residual tight lines.
Final log: zero warnings, zero over/underfull boxes.
Generated-with: Claude Fable 5
… bug Verified against matlab_scripts/pc_min.m (identical INB loop, clamp, and partial-area formulas; Emanuel's comment: 'Find maximum level of positive buoyancy, INB') and the git history (the Python CAPE module and pc_min.m enter together in the first pyPI commit, April 2020, with the INB block translated verbatim and never altered). Neither Emanuel (1994) nor Bister & Emanuel (2002) specifies a multi-crossing rule, so the selector is an undocumented implementation choice of the reference code. State this explicitly in Convention 3 and in the conclusions. Generated-with: Claude Fable 5
Panel (a) compared the legacy and repaired maps in overlapping solid
blue/teal with terse legend labels. Restyle using the document's
established convention colors: legacy E_top map in solid red, repaired
E_max map in dashed teal drawn on top, so coincidence (red through the
dash gaps) and divergence are both visible; spell out the intent in the
legend and title ("same map, one changed ingredient").
The restyle also exposed that the old caption was wrong: the two maps
coincide only on the right branch (where INB already selects the
dominant cell); on the left branch they differ by the full effect of
the 222 J/kg clamp. Caption corrected.
Generated-with: Claude Fable 5
Include discontinuity_analysis.pdf in the repository so the analysis is readable without a TeX toolchain; keep ignoring LaTeX intermediates. Generated-with: Claude Fable 5
New subsection (sec:jumpmag): the 0.838 hPa jump decomposes as |dg| ~ g * dE_M/(Rd*Tv) * (1 - k/2 * SST/T_O), with all factors verified against the real cape()/pi(): partial cancellation c1 = 0.342 (E_M enters CAT twice with opposite signs), thermal denominator Rd*Tv = 87.0 kJ/kg, level g ~ 951 hPa. Prediction 0.838 vs measured 0.8381; the 1995 column's own factors give 0.840 vs measured 0.839. Generic failure scale: ~200 J/kg of clamped CAPE x 3.7e-3 hPa/(J/kg) = 0.7-0.9 hPa, just above the 0.5 hPa tolerance - the silent-failure window - reached by both known cases from different ingredients. Correction: the clamp does NOT maximize (or cause) the jump. Disabling max(.,0) leaves the discontinuity intact and increases the measured jump to 1.0845 hPa (E_M flips -65.90 <-> 224.23); the clamp pins one branch at exactly 0 and attenuates the jump by a quarter. Convention 3 and the mechanism section no longer overstate the clamp's role, and the -60.9 (phase A) vs -65.90 (at x_d) raw integrals are now disambiguated. Generated-with: Claude Fable 5
The table conflated the trigger with the deciding feature's energy and presented the g-jump as a separate 'map discontinuity'. Restructured: the perturbation that flips the outcome is infinitesimal (a sign crossing - that is what discontinuity means); ~1e-2 J/kg is the feature's own energy content; and the 0.838 hPa jump of g is the same discontinuity as the 224 J/kg jump of E_M, attenuated through eq. (2), not a third object. Generated-with: Claude Fable 5
Generated-with: Claude Fable 5
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
This PR adds
discontinuity_analysis/— a self-contained mathematical account (19-page PDF, LaTeX source, 13 figures, and all analysis scripts) of the period-2 convergence failure reported in dgilford#77, written for a reader with no atmospheric-physics background. The compiled PDF is included:discontinuity_analysis.pdf.Key findings
gjumps downward by 0.837 hPa at x_d ≈ 950.8289947515 hPa; the residualg(x) − xchanges sign only across the jump. The bit-exact 2-cycle (950.653… ↔ 951.279…) is the canonical dynamics of a contraction with a diagonal-straddling gap.cape(): "integrate signed buoyancy to the topmost grid level with positive buoyancy, clamp at zero." A ±10⁻¹¹ K sign flicker at the 350 hPa level toggles 287 J/kg of negative area; eyewall CAPE flips between 224.23 J/kg and 0. The deciding feature is energetically negligible (~10⁻² J/kg) and dynamically unreachable (W = −61 J/kg there).pc_min.m("Find maximum level of positive buoyancy, INB"), and neither Emanuel (1994) nor Bister & Emanuel (2002) specifies a rule for multiple buoyant layers.g; the fatal diagonal-in-gap placement occupies a ~0.04 K SST window for this column — exactly what attribution-style SST sweeps find.max_q W(q)over the same quadrature and candidates (continuous by Berge's maximum theorem; matches forced-eyewall-ascent physics) and solve with residual-checked bracketing. Exact corrected roots: 950.3136889226 hPa (VMAX 83.9065 m/s, PMIN 920.7517 hPa) for Add failing test to illustrate failure of convergence dgilford/tcpyPI#77 and 960.7687 hPa (VMAX 74.4149 m/s) for the 1995 regression case. On the 7,093-column 2024 benchmark: zero flag changes, 98.1% of VMAX within 0.1 m/s, and a 1.9% tail confined to configurations where the legacy value was fragile by construction. The period-2 "rescue" becomes dead code (its midpoint was never a root).Contents
discontinuity_analysis.pdf/.tex— the document (rebuild withlatexmk -pdf; final log is warning-clean)figures/— all figures (PDF + PNG)scripts/— the analysis scripts behind every number and figure (run againstdgilford/tcpyPIbranchmodernize_eg, commitabc4b4e); every headline quantity was reproduced by a second, independently written implementationSuggested follow-ups (not in this PR)
Staged implementation plan is in §7.5 of the document: diagnostics first, then the max-work kernel behind an explicit option, outflow/tie semantics (including flagging ptop-clipped LNBs), and a bracketed solver. The convention change should be coordinated upstream before becoming a default.
🤖 Generated with Claude Code