Add an escape regime dispatcher with per-species fractionation - #34
Draft
maraattia wants to merge 92 commits into
Draft
Add an escape regime dispatcher with per-species fractionation#34maraattia wants to merge 92 commits into
maraattia wants to merge 92 commits into
Conversation
…assertions Pin the tidally corrected critical exobase temperature away from the x = 1 point where a dropped 1/x factor is invisible, pin the nitrogen recombination coefficient with a swapped-coefficient guard, pin both Tang timescales and both verdict directions, prove the diagnostics container is boxed by sabotaging its producers instead of re-running, pin the scaled kinetic diameters in the assembled table, and assert the shallow-column CO2 trapping branch. The Volkov flat-factor test now exercises the source function beside the published-coefficient oracle. Document the wind-temperature sensitivity of the EL/RR crossover and the dense-base clamp behavior of the thermostat in their modules.
maraattia
added
enhancement
The physical base pressure the base method targets is a property of the level, not a warning, so it moves onto the level dict as p_physical and into the dispatch diagnostics beside the level it describes. The flags a dispatch returns are now warnings only: a call with nothing to report returns an empty flags dict instead of one that always carried a value.
Runs the framework over synthetic atmospheres: one verdict read field by field, flux sweeps that cross two boundaries on two compositions, the boil-off and overflow labels, the diagnostics container, the four knobs that move a boundary, the hysteresis window, and one frozen atmosphere dispatched along a stellar XUV history. Every step is a function that returns its results, so the pieces can be imported one at a time, and the figure follows the framework figure conventions with the escape domain color and one marker per regime.
Three pages, one per purpose. The tutorial drives the framework on synthetic atmospheres: one verdict read field by field, two boundary crossings on one planet with the criterion band measured around one of them, the boil-off and overflow labels, the diagnostics container in the order the questions occur, the four knobs that move a boundary, and one atmosphere dispatched along a stellar history. The results reference tabulates every field, flag, and diagnostics group, with the two reporting conventions a caller needs: a rate floor for numerical noise and the depletion timescale for relevance. The triage guide goes from a symptom to the entry that explains it.
The example tests exercise the shipped script's own entry points: the regime sequence of the documented flux sweep, the per-species closure on both the fractionating and the hydrostatic path, the zero-flux limit, the error contract on a malformed state, the bisection that produces the boundary fluxes the tutorial quotes, and the label change along a stellar history with the stellar lookup mocked. The dispatcher test now pins the unclamped base pressure the diagnostics report, and both rule files record the naming exception this file takes.
The new test carried its own tier marker beside the module-level one, which the structure validator rejects: exactly one tier per test, and the module already declares smoke.
The tutorial's opening snippet imports the worked example by its path from the repository root, which is what a reader with a clone would type, and that resolves only with the root on sys.path. Whether it is there depends on how pytest was started: running it as a module puts the working directory on the path and the console script does not, so the test passed locally and failed in CI. It now prepends the root itself.
The snippet test executed the whole page, and the page closes on a track dispatched along a real MORS stellar history, which reads the Spada grid under FWL_DATA. The unit and smoke tier does not download that grid, so the test could not pass on a fresh runner however it was invoked. That block now runs in the integration tier beside the other test that uses the real tracks, with the same character-for-character comparison, and the fast tier covers the other eighteen output blocks. The block is excluded by name rather than by whether the data happens to be present, so both tiers check the same thing on every machine, and the fast tier is verified against an empty data root.
Twelve accuracy corrections, each checked against the code or the arithmetic. The wind-base clamp does not depend on a single profile top: the Lopez level is tens of nanobars on an Earth-mass carbon dioxide planet and below a nanobar on a low-gravity hydrogen envelope, where a profile stopping at a nanobar misses it, so both pages that promised one number now say to read the flag. The tidal setting divides the interior-luminosity cap as well as the energy-limited rate. The diagnostics API page was missing the function behind a documented group. The floor illustration said 66 grams a year where two nanograms a second gives 63. The README still described the tidal domain from before the reference radius followed the scaling. The bolometric instellation is carried and validated but no branch reads it, which the parameter table now discloses. The hysteresis flag says the window was available, not that it changed anything. The self-consistency span is 65 to 767 Myr, not 40 to 800. The boil-off label does not test whether the atmosphere is inflated past its sonic radius. And the evaluation order listed the Roche screen before the bolometric residual in both the page and the module comments, while the code runs the residual first, which matters because the radius the screen tests belongs to whichever branch actually won. The hydrostatic section also called an effusion velocity a flux, used an undefined flux symbol in the harmonic mean, and left the base density and its mean-mass convention unstated; all four are now written out.
The accent is restored on Chassefière, in the module, its tests, and two pages. The Caldiroli threshold potential is quoted on the untidal binding energy, which is the convention the screen uses and the one Salz share; the tidal factor enters the rate and not the potential the threshold is measured against. The interior-luminosity cap is their Eq. 8, the update equation that carries the minimum against the cooling luminosity, not Eq. 9, which is the period distribution of their population model. The upper edge of the Knudsen band is not Johnson's number: their simulations give 0.1 and 1, and the 3 comes from Chatterjee and Pierrehumbert arguing the energy limit may survive that far, so the band is asymmetric in what supports it and the page now says which end rests on what. And the activation parameter is written with the composition mean molecular mass where Fossati write the atomic hydrogen mass, a factor mu over m_H between the two, which is now stated where the equation is given rather than left for a reader to discover by comparing numbers.
The metaphor family was the project's own coinage for the coefficient fallback order, not a code identifier and not the field's standard name, so the prose rules replace it with what it describes: 33 occurrences across four modules, three test files, and three pages. Two test names carried it and are renamed with it. The same pass takes the comma after e.g. and i.e. in three shipped docstrings, settles the prefix hyphenation on the unhyphenated form for nonthermal and deexcitation, which the branch had split down the middle, writes numbers above ten as numerals in three places, and drops a borrowed noun from a dispatcher docstring. The worked example said it writes one figure and writes two, as a PDF and a PNG each. The energy-limited API entry was rendered on two pages and is now rendered once, on the released import path, with a pointer from the module page. And the validation page for that function pointed at a docstring and an overview this branch moved.
Three corrections of the same kind as the documentation pass. The figures named Spline Sans Mono unconditionally while the colours that come with it sat behind an optional import, so without proteus-mpl matplotlib substituted a face with a warning and a reader's figures differed from the committed ones with nothing saying why; the family now falls back to a generic monospace where the package is absent. The module notes on the recombination chain still said the flux scaling separates genuine saturation from barometric suppression, which the barometric factor does and the flux exponent does not. And the Jupiter constants had been inserted into the middle of the Earth block, leaving six Earth quantities under the Jupiter heading.
The tidally driven transfer through L1 from Jackson et al. (2017), ApJ 835, 145: their Eq. 3 rate with the Eq. 10 curvature fit, the Eggleton (1983) lobe radius, and the Eq. 14 volume averaged potential, closed form throughout. At and beyond lobe contact the exponent clamps at the lobe filling boundary value, which is their Figure 5 case. Validated against their Table 2 planet rates, their Figure 5 binary curve, and the exact corotating Roche potential at a numerically solved L1 point.
The nozzle competes in the final comparison on both sides of the activation gate, wherever the overflow description applies: the isothermal sonic radius must reach the L1 distance (Jackson et al. 2017, their Figure 9 crossover), because inward of that a spherical wind chokes at its own sonic surface first, and the nozzle area grows as the cube of the separation, which otherwise hands a loosely bound envelope an unbounded rate far from its lobe. A candidate below the one proton per year floor does not compete either: this label boundary is a rate crossing, and a crossing between two numerically empty numbers would rename the deeply bound corner on no content. A nozzle win labels roche_overflow with the transfer rate itself, continuous at the crossing; the geometric rename keeps its bound flow meaning and its rate. The subflag still reads the geometry, with a third value, nozzle, for a win whose structure sits inside the Hill sphere. The flow carries no energy cap, faithful to the primary, and the lift power is reported beside the interior and intercepted luminosities. nozzle_temperature selects the photospheric level or the thermostat wind state; eccentric wins are periapsis evaluated and flagged.
The regimes page derives the nozzle candidate, its applicability criterion, and the conventions it carries; the results reference separates the rate crossing reading from the geometric rename and adds the nozzle diagnostics group and the two new flags; the parameters page gains nozzle_temperature; the limitations and troubleshooting pages state what the transfer rate assumes and how to read a roche_overflow verdict, including the torque balance regulator the module deliberately does not compute.
The nozzle group brought the container to nineteen entries on a call that also carries the species closure, and to eighteen without it, so the reference page and the tutorial both understated it by one. The tutorial also sent the reader to three cases in the troubleshooting guide, which lists four. The tutorial's own printed key list already showed nineteen, so the pages disagreed with output quoted on the same page.
A secular caller integrates over many orbital periods, so the dispatched transfer is now the time average of Jackson et al. (2017) Eq. 3 over the orbit, Kepler-weighted through the eccentric anomaly and duty-cycled over the arc where the overflow description applies. That arc surrounds periapsis, because the L1 distance grows with separation while the sonic radius does not, and the criterion now uses the L1 distance itself rather than the Hill radius that stood in for it. The lift power is rebuilt from the barrier the rate applied plus the acceleration to the sonic speed, which is the heat an isothermal flow demands, so it stays finite at lobe contact instead of reading a potential the clamp had already discarded. A circular orbit is unchanged to machine precision.
The subflag claimed to say whether the atmosphere reaches its lobe while testing it against the Hill radius, which sits about 43 percent further out, so a structure past its own lobe could still be reported as staying inside. It now tests the lobe radius, which is the critical surface, and the third value is geometric rather than named after the candidate that won: the mechanism is already carried by the rate branch. The overflow label and its subflag are assigned in one block instead of two, so the precedence rule is written once. The winning branch's flow radius is no longer replaced by the lobe radius, which had pinned the reported ratio at a constant and left the near miss warning suppressed only by clause order.
Under the wind temperature setting the launch density came from the photosphere while the sound speed evaluating the barrier came from the wind, so the two belonged to different structures and the Bernoulli cancellation behind the launch-level convention did not hold: the rate tracked whichever level was nominated, moving by nearly three decades across the plausible range. The launch level now sits on the wind's own isothermal column, anchored at the wind base with the ideal-gas density there for the wind's temperature and mean mass, since pressure is continuous across the temperature transition and density is not. Across a factor of four in launch radius the rate now holds to half a percent, against a factor of eleven on the profile's colder column, and a test pins both halves. The column places the launch level and makes no claim about structure below the anchor.
Contributor
|
@maraattia I can review the docs part of this PR now, should I or are you still making changes? |
Author
|
@stuitje No not yet please, when I open the PR I'll explicitly tag you. Thanks! |
Each candidate's flags were merged into the result as they were computed, so a dispatched rate arrived carrying warnings about branches that lost: across three thousand random states, a third of the overflow results still advertised the luminosity cap of a bolometric candidate they had beaten, and hydrostatic extension warnings survived onto rates the hydrostatic branch did not produce. Each candidate now holds its own flags and only the winner's are merged, so the flag set always describes the dispatched rate. Input hygiene stays separate, since a stale input or a clamped base is a property of the state rather than of a candidate. This replaces two hand-maintained lists of flags to discard, which had already drifted once.
The rate branch reported alongside an overflow verdict took a value that no exported tuple contained, so a consumer validating it against the regime labels failed silently and nothing enumerated the legal set. The transfer now reports the label itself, which puts every value back inside the regime labels: reading roche_overflow there means the L1 transfer was dispatched, and any other value under that label means the geometric screen renamed a state whose rate that branch produced. The geometric case needs no name of its own, since it ships whatever rate won before the screen and the screen is already visible in the flags.
The overflow module is documented public API with its own reference page, but only the mass ratio was checked, so a direct caller got a negative mass-loss rate from a negative density, a complex cube root from a negative mass ratio, and a bare division by zero from a launch level at the origin or an eccentricity of one. Every exported function now raises with a message naming the argument, in the style the mass-ratio check already used. Nothing reachable through the dispatcher changes, since the state and profile validation upstream already excluded these cases.
The transfer rate is invariant to the launch level only along a column whose sound speed matches the density it carries, and every dispatcher fixture was isothermal at the equilibrium temperature, so no test could see that cancellation fail or tell a profile temperature from an equilibrium one. A fixture with a temperature rising through the column now measures the residual: a factor 2.4 across three decades of launch level, against 1.14 on the isothermal control. Three gaps close with it: the flow radius reported on an overflow win is pinned as the wind's rather than the lobe's, the bulk split on that branch is pinned against the Jeans split its rival would have produced, and the boxedness check gains an overflow state. The unreachable second clause of the saturation test is gone, since the potential rises monotonically to the lobe at every mass ratio and the exponent can only reach zero where the geometry already has.
The validation page described Jackson et al.'s Table 2 as six rows with two near-lobe hot Jupiters left unpinned, and quoted a residual of 1.4 to 1.7 that only exists after a radius conversion this implementation deliberately omits. The table lists twenty-one objects, seventeen of which this implementation does not reproduce, by factors running to 3.9e6 and rising monotonically with launch depth. The pin now covers every row whose launch level sits below 0.4 lobe radii, four instead of two, at the depth where the potential approximation holds, and the page states the two causes of the shortfall: the omitted radius conversion, which needs 2.2 to 5.2 percent in radius against an exponential sensitivity of 9.7 to 276, and the assumed masses the candidate rows carry. Three quoted figures are corrected: the agreement against the total potential is dominated by a constant that cancels, the Figure 5 anchor asserts at 30 percent and exercises no barrier, and the curvature row no longer cites the paper's own printed asymptotic, which carries half the coefficient the fit uses.
Four statements did not survive a check against Jackson et al. (2017). The isothermal flow was described as one the radiation field maintains, where their Section 3 names the neglected heating and cooling balance and calls isothermality an important limitation. The bound-flow rate under a geometric rename was called a lower limit on the tidal transfer without qualification, which has the wrong sign wherever the nozzle candidate was applicable and lost, since the module's own estimate of that transfer is then the smaller number. The fractionation page still said the overflow label never changes the rate. And the rate floor was quoted under three different names for one constant. The regimes page also splits into six paragraphs from one of nearly eight hundred words, and states the floor's reach: it marks what is distinguishable from zero in floating point, so about one overflow win in twenty is dispatched below 1e-20 kg per second, and the depletion screen rather than the floor is what says whether a rate matters.
The page recorded its last comparison against the source one day early.
Past the activation gate the luminosity-capped bolometric candidate is still computed and reported, but it now competes for the dispatched rate only under the new `residual` setting, whose default `off` leaves the XUV or hydrostatic verdict standing. In the window just past the gate the closed-form Parker and Bondi rates have not yet shut off, so the candidate there is the interior-luminosity cap itself, the core-powered rate whose persistence Tang et al. (2024) dispute, and it outran the XUV rate at interior fluxes near a watt per square meter on a three Earth-mass hydrogen envelope, enough to strip a one percent envelope in a few hundred million years. The bolometric diagnostics gain `residual_mode` and `competes` so a reported candidate can be told apart from a losing one. The tests whose families were built on a residual win admit it explicitly, the overflow flag case moves to a nozzle win under the defaults, and the tutorial's overflow example admits the residual and states what the default returns instead.
The regimes page states the dispute the residual sits on and the side the default takes: Gupta & Schlichting (2019) keep a bolometric wind fed by the cooling interior alive for gigayears, Tang et al. (2024) find it removes at most a tenth of a percent of the envelope once boil-off is initialized self-consistently, and the luminosity cap does not make the difference small, since just past the gate the residual is the cap itself and exceeds the XUV rate by more than a decade at interior fluxes of a few watts per square meter. The evaluation order, the flowchart, and the paragraph that previously claimed the two readings differ by little now say what the setting does instead. The results and parameter references gain the `residual` knob, the two new bolometric diagnostics, and the qualification that the residual flags are reachable only when it is admitted; the model overview and the troubleshooting guide follow.
The Bondi-limited cap is evaluated at the same wind temperature as the Parker rate, T_eq / 2^(1/4), where Gupta & Schlichting (2020) write their Eq. (10) at T_eq; Misener et al. (2025) recommend the cooler value for that same isothermal form, so the cap follows the correction rather than the original. The module header and the regimes page now say so, give the factor between the two conventions, which grows exponentially with the launch-level Jeans parameter and reaches about 30 at the shutoff, and state why it seldom reaches the dispatched rate: the cap and the Parker rate describe one isothermal wind, sitting a constant factor of about e^(3/2) apart at a launch level of unit optical depth, and on the default launch level, which is optically thick to the supplied opacity on an inflated hydrogen envelope, the cap does not bind at all. The reported launch-level optical depth is what tells the two cases apart.
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
2 participants
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.
Description
This PR adds an escape regime dispatcher to ZEPHYRUS. One call,
zephyrus.dispatch, takes a planetary state (scalars plus an atmosphere profile) and returns one regime label, one bulk mass-loss rate, per-species rates that sum to it, flags, and a boxed diagnostics container. Every physically posed input returns one and only one labeled result; exceptions are reserved for malformed input.Every state receives one of five regime labels (boiloff, hydrodynamic:EL, hydrodynamic:RR, hydrostatic, roche_overflow), decided in a fixed evaluation order:
roche_overflowand keeps the rate and the split its own branch produced, because that boundary is a comparison between rate candidates and the two sides of it hold the same branch; substituting another branch's formula would make the dispatched rate jump at a line the physics puts nowhere in particular. The label therefore asserts a geometry, and the rate beside it is a bound-flow estimate, a lower limit on what tides would do. A subflag separates an atmosphere whose own extent reaches the lobe from one where only the would-be sonic surface would sit outside it (Owen & Jackson 2012). A hysteresis window on the switch threshold is available for evolutionary use so tracks cannot chatter between branches.Confirmed hydrodynamic outflows are partitioned over species by a simultaneous N-species fractionation closure (the constant-composition closure of the Zahnle et al. 1990 multispecies wind system with active-set dropout, generalizing Hunten et al. 1987), on a binary-diffusion-coefficient library in which every pair carries a provenance class. The split follows the branch that produced the rate rather than the label, so a wind that reaches the Hill sphere keeps its fractionation.
New modules (flat layout):
composition,profiles,knudsen,diffusion,atomic_data,thermostat,boiloff,hydrodynamic,hydrostatic,fractionation,diagnostics,dispatcher. The energy-limited physics now lives inhydrodynamic.pybeside the other rate prescriptions, andescape.pyis a thin compatibility re-export, sofrom zephyrus.escape import EL_escapeand the PROTEUS version pin keep working unchanged; the released function's behavior is intact and additive apart from the tidal convention fix below. The dispatcher is a standalone API until the PROTEUS-side wiring lands in a follow-up PR (config schema, wrapper dispatch, profile transport, per-element reservoir debits).EL_escapeitself receives three documentation and convention fixes: the R_XUV docstring now names which level the 20 mbar value is and distinguishes it from the nanobar level where the wind is launched (Lopez 2017); the tidal factor's xi is now measured from the radius thescalingargument selects (R_p for the Erkaev form, restoring that paper's own convention; R_XUV for the single-radius form), which changes standalonetidal_contribution=True, scaling=2results; and the efficiency docstring carries the Caldiroli et al. (2022) threshold-potential caveat.constants.pymoves to the exact SI defined values and gains the proton mass, atomic mass constant, Planck constant, and eV conversion;planets_parameters.pygains the IAU nominal Jupiter values.The documentation is restructured around this model: a refreshed landing page, a hub page presenting the two loss channels (continuous escape and giant impacts) and the five regimes with a simplified decision flowchart, and topic pages covering the energy-limited default, the regime framework, fractionation (citing Attia & Lichtenberg 2026, in prep.), and the impact law, each written with every quantity defined and a bibliography. Alongside them, for the dispatcher specifically: a tutorial that crosses the boundaries on a worked planet, a reference page documenting every flag and diagnostics group the result carries, and a troubleshooting guide keyed to the symptoms a first user hits. Every printed number in those pages is the verbatim output of a snippet the reader can run, the two figures come from
examples/demo_dispatcher/, which is itself covered by tests, and the mermaid flowcharts now render on the built site instead of arriving as diagram source.Three boundaries are worth calling out for reviewers. Every regime criterion carries a stated physical band (the Knudsen threshold band 0.1 to 3, the boil-off activation band 15 to 35), and the diagnostics container prints the counterfactual labels at the band edges beside every verdict; the container is reporting only and nothing in the control flow reads it. Two known limitations are flagged on every affected result rather than hidden: hydrostatic heavy-element rates are lower limits (nonthermal channels are absent), and grid points where the neutral and plasma escape-temperature conventions disagree are marked contested with both branch rates recorded. And because a strongly bound heavy atmosphere returns rates dozens of decades below anything physical, every verdict reports whether it clears one proton per Julian year, which is the smallest rate with content; the module never applies that floor, since what counts as negligible belongs to the caller.
Validation of changes
ruff checkclean oversrc/,tests/, andtools/, and the documentation site builds without issues. (ruff checkover the whole tree reports two pre-existing findings inexamples/demo_earth/, untouched by this PR and not run by CI.)docs/Validation/. Highlights: the Erkaev et al. (2007) Table 1 enhancement factors within 1%, on both consumers of that factor; the Caldiroli et al. (2022) efficiency fit at six spot points within 5% and a clean rejection of its complex-valued region; the Murray-Clay et al. (2009) fiducial hot Jupiter (base Jeans parameter, sonic radius, and printed sonic-point Knudsen numbers); the Laricchiuta et al. (2009) collision integrals reproducing four measured viscosities within 7%; the Yelle (2024) Figure 1 Mars hydrogen flux through the Jeans-to-diffusion transition; the Volkov et al. (2011) coefficient table; the Owen & Wu (2016) sonic-point limit and shutoff; the Salz et al. (2016) wind-versus-thermosphere thresholds, reported with their hydrogen-dominated scope rather than applied; and the Hunten et al. (1987) Earth, Mars, and Venus anchors.pytest.approxkeeps its defaultabs=1e-12alongside anyrel=, and takes the larger, so a pin on a cross section, a recombination coefficient, or a physical constant accepted any value at all: inflating the proton mass, atomic mass constant, and Planck constant by a tenth, or a collision integral by 30%, left the suite green. Both are now caught, and the boxedness check runs on one state per branch at a rate of order unity or above rather than one near the denormal floor.EL_escapeat 1e-9 relative, and the worked example is driven end to end by its own test module.Checklist
Relevant people
@timlichtenberg @nichollsh @EmmaPostolec