Repository navigation
Split bulk H, C, N and S between core and mantle with partition laws - #47
Merged
Merged
Conversation
Tests for the Suer et al. 2017 S law (Eq. 7, log10 in every term, a from Table 3) at a hand value and against the 9 runs of their Table 1, the optional f and log10_C_S terms of a log10_linear entry, the range distance per axis, named laws as TOML entries or callables, and split_core_mantle (closed-form core fraction, element ledger to 1e-12, unclamped out-of-range evaluation with flags, input domain). The stored outputs of main 75faea5 pin the packaged laws and 3 solves, so the solve path stays as it is.
metal.py gains named laws (get_law, register_law) that hold TOML entries and callables with a declared basis and range, the optional f log10 X_FeO and log10_C_S terms of a log10_linear entry, and range_distance, the distance outside each range per axis. The packaged file partition_laws_core.toml holds the S law of Suer et al. 2017 (Eq. 7 with log10 in every term, a from Table 3, log10 C_S -5.372). core_split.split_core_mantle splits a bulk H, C, N and S inventory between core and mantle at an equilibration point (default 60 GPa, 4000 K) by a fixed point on the metal composition, evaluates each law unclamped with a flag and its distance when it runs outside its range, refuses O and Si as bulk elements and asserts bulk = core + mantle to 1e-12 for each element. The legend_clear fixture checks that a plot legend covers no line, marker or band. The unknown-key example of the schema test is now g, since f is a schema key.
The split takes the core and mantle shares as fractions of order 1 times the bulk, so a small bulk does not underflow in bulk / (M_silicate + D M_metal), and the ledger tolerance has a floor of 2 ulp of 0 for the round-off of a subnormal bulk. The tests of the registry and the split assert each clause on its own and compare floats with a stated tolerance.
The Explanation of the metal reservoir gains the optional f and log10_C_S terms, a section on the split at an equilibration point (no clamp, flags with distances, laws by name) and its limits: O and Si are not split, the single-cation melt basis of the S law is an assumption, and one C_S is used for every melt. The Validation page gains the 9 runs of Suer et al. 2017 Table 1 (run 7 read with X_FeO 0.20 from its IW and Table 2, a stated miss at +0.495 dex) and the closed-form split. The Reference lists the new keys, the named laws and the inputs and outputs of split_core_mantle; a tutorial runs the split, its flags, a registered law and the ledger. core_split.py is a physics source of the test rules.
Run 7 of Table 1 uses X_FeO 0.20 in place of the printed 0.02: its IW of -1.06 implies 0.2007 and its Table 2 FeO gives 0.204. The test pins its residual at +0.495 dex, outside the 0.47 dex of the other 8 runs, and the Validation page states the reading with its pages. The source of the S law names the Table 3 d as the O interaction term, and register_law gets tests for an entry with declarations, an entry of another element and a law that is neither.
A registered callable must receive P, T, dIW, the metal and melt fractions and NBO/T, and return a finite real (a bool is refused, a numpy real accepted); a declared range that is an array, an incomplete entry and a taken name without replace are refused; get_law returns a copy; a mass-basis law does not read the melt oxides. The split refuses law keys outside its bulk, None, a bare callable and a law of another element, names an overflow of D times the metal, a bulk beyond the float range and a bulk the metal cannot hold, and warns once per entry with rounded distances. The exchange refusal covers both optional keys, and legend_clear gets its own tests, scatter points included.
…s flagged The S law of the split must take X_FeO = X_Fe 10^(dIW/2) from the metal and dIW, so D_S rises with dIW by 10^0.75 over 1.5 dex, and the melt oxides change D_S only through the silicate units; NBO/T must move D_C by 10^(0.118 dNBO/T). A metal richer in S than FeS must set beyond_FeS in a metal flag, with the ledger still closed.
A law with an FeO term now takes X_FeO = X_Fe 10^(dIW/2) of the metal of the pass (ideal mixing, as Suer et al. 2017 define IW), so dIW is the one redox input of the split and the melt oxides set only the silicate units. The flags gain a metal entry with the light-element and S mole fractions and beyond_FeS above the S of FeS, with no clamp. The split refuses law keys outside its bulk and laws that are not D laws of their element, names an overflow of D times the metal and a metal with no iron left, and warns once per law with rounded distances; the range flag comes from the distances. A registered callable also gets NBO/T, a returned bool is refused and a numpy real accepted, get_law returns a copy, register_law checks entry fields and takes replace, and mass_D reads the melt oxides only for a law that needs them. legend_clear checks scatter points and closes the figure.
The Explanation states that dIW is the one redox input of the split, with X_FeO = x_Fe 10^(dIW/2) for the S law and the S core share at dIW -3, -1.5 and -0.5, and describes the metal flag and the checks of the laws argument; the Reference lists the NBO/T argument of a law function, replace, the copies of the registry and the metal flag; the tutorial gains a step on the redox state and its numbers come from a run. The stored outputs of the solve and the packaged laws are reference_outputs.json, compared to the bit on the platform that stored them and to 1e-12 elsewhere. The bad-input tests share one call, a Hypothesis test checks the ledger with the default laws and metal feedback, and the Suer test constants carry no plan labels.
… S/Fe flag A dIW above the ideal-mixing limit -2 log10 x_Fe must raise a ValueError on dIW, also at dIW 700; a law function in the split must get the melt fractions of the oxides; beyond_FeS must compare S with Fe also with O and Si in the metal, with S_over_Fe and x_S checked against a hand count. Hand-made laws (a list, an entry without coefficients, function or with a string value) must raise by name; a law warns once per set of ranges, with rounded distances; a packaged law cannot be replaced; a registered entry is stored as a copy and a function is shared; a returned int past the float range or a complex is refused by name. The Hypothesis ledger test checks the D of the returned metal, and the S law guard compares at dIW -3.
mass_D now takes X_FeO = x_Fe 10^(dIW/2) of the metal for every law with an FeO term, in the metal reservoir and the split alike, computed in log space, and raises a ValueError on dIW above the ideal-mixing limit -2 log10 x_Fe; a law function gets the melt fractions of the oxides as given. beyond_FeS compares the S and Fe moles of the metal and the flags give S_over_Fe. check_law validates an entry for register_law and the split as a partition file does; get_law copies entry data and returns a law function itself; a packaged law cannot be replaced; a law warns once per set of ranges; a returned int past the float range raises by name; a failing law names its element. The docs state the dIW limit, the shared FeO rule, the S to Fe flag and the absolute floor of the fixed-point step; the tutorial checks its ledger against 1e-12.
… with Fe at round-off Each bad laws argument (an empty list, a list, a scalar or dict range, an entry without coefficients or function, a string coefficient) must raise its own error type naming the fault; dIW -700 must raise a ValueError naming dIW; a metal with Fe at round-off must give finite flags. Guards pin the dIW limit of the start metal, the S law FeO of mass_D against the from_FeO dIW of the metal reservoir, and the S core shares quoted on the Explanation page.
mass_D and the metal flags take the Fe moles from the Fe mass fraction, not as total minus light elements, so a metal with Fe at round-off gives finite values; mass_D checks dIW first and raises a ValueError naming dIW when X_FeO would underflow. check_law passes a range other than None to the range check, so a scalar or dict range raises by name, and a callable entry without a name is labelled entry. The split takes laws=None as the default and refuses any other value that is not a dict, and its warning key holds the full source. The docs give the dIW limit of the start metal, +0.06 for 2 wt% O and +0.11 with 3 wt% Si, and the tutorial names S_over_Fe.
The interval [0, 1) in the w_fixed row read as a link reference together with a later [e]; the row now says from 0 to below 1.
New tests pin the messages of the w_fixed rules and of an unknown law name, the 2 ulp floor of the ledger on a subnormal bulk, the Fe moles of a metal with Fe at 1e-12 by mass, the flag of a law with an unknown range, a ledger that falls short of the bulk, a point range and a range of 3 values, X_FeO of exactly 1 at a pure Fe metal and the underflow bound at X_FeO 1e-300.
…rance On Linux the nightly found the solve without metal at another root within the acceptance rtol 1e-5 (about 4e-6 relative) and different round-off residuals in the metal cases, and pytest.approx refused the nested law grid. Off the stamping platform the solve outputs are now compared at rel 1e-4 without the *_res columns, and the law grid is flattened and compared at rel 1e-12; on the stamping platform both stay exact.
The split must return the bisection root of an S-only metal at the 3 S-rich, metal-poor cases where the plain substitution oscillates or leaves the iron-bearing range, with the ledger closed and beyond_FeS set, and must return the bulk as mantle with no metal at any dIW. split_reference.json stores the outputs of 32 splits that converge today (the tutorial, a trace of S, a grid of bulk, dIW and metal mass), compared to the bit on the stamping platform. The S law ranges and their flags at dIW -3 and -0.5 are pinned, the O and Si exchange is compared on the stored grid, the off-platform solve compare uses abs 0, and a test comment states that the ln reading of Eq. 8 lands 5.4 dex higher.
Where the metal takes so much S that D_S falls fast with the metal S, the plain passes of the split swing about the root or leave the metal without Fe. The split then starts again with relaxed passes: each element moves 1/(1 - s) of its step, with s the negative secant slope of its last 2 passes, and a pass adds at most half of the Fe left; an error in the first pass still raises. The relaxed passes run only when the plain passes fail. With no metal the split evaluates no law and returns the bulk as mantle. The Explanation states the bases of X_FeO and C_S as Suer et al. 2017 give them and the domain checked for the relaxed passes, and the code architecture page names range_distance for the flags.
…plit per part The non-settling mock now runs 50 plain and 200 relaxed passes, and an S bulk above the metal mass settles at its bisection root past FeS. Off the stamping platform the stored split is compared per part (core, mantle, metal), since pytest.approx takes no nested mapping.
A rising map takes full steps, the first relaxed move takes at most half of the iron left (halved by 2), and the relaxed passes start again from the start metal; these kill the path mutants that the focused mutation run left (start point, the clamp of the secant slope, the default step, the cap and its halving).
…a coupled S-rich metal A C bulk that fills the metal and an N bulk beyond what the metal holds must raise the named no-iron ValueError; S roots with 1e-6 and 2e-5 of iron must match bisection although g falls by about 1e6 per unit w; a coupled S, C, H, N metal at 2.5 % iron must reach the root of a damped plain iteration.
… rule When plain substitution fails, the relaxed passes first damp each element by 1/(1+|s|) with its own secant slope and then, if that sweep does not settle, damp all elements by the smallest factor. Each pass takes at most half of the iron left. A metal with iron below 1e-8 by mass raises the named ValueError, as does a metal with no iron outside H and C in mass_D. A root closer than the float grid of w resolves is accepted when the slope is negative, the move is within 2 ulp and the step is within 1e-9 of g. The relaxed passes also require the melt H and C basis to settle.
…only metal guard A move within 2 ulp must not end the relaxed passes while the step is 2e-6 of g; the passes must wait until the melt H has settled; sweep case 1870 must settle on the common step after the per-element passes; a metal with w_H = 1 - 2**-53 must raise the named ValueError.
The steep-step test compares w_A with pytest.approx and the melt-settle test also checks that the passes return the start metal.
The sweep of 3000 splits gives D residuals up to 7.8e-4 for H and 5.5e-4 for S on the plain passes, all in splits that carry H.
The sweep gives 5.49e-4 for S.
Sweep splits 1532 (H) and 1721 (H and S) have a small silicate mass, so the melt H that sets the silicate units of D still moves when the metal has settled; at the returned split D must reproduce the core of every element to 1e-9.
The plain passes now end only when the metal and the melt H and C mass fractions settle, as the relaxed passes do; the melt H and C set the silicate units of the mole bases of D. Of 2555 random splits that converge on the plain passes, 896 change, by the distance to the fixed point (the one-pass residual times 1/(1 - rho) of the plain map), with a metal change of at most 1.8e-9. Four stored reference splits change and are restamped; the largest relative change per case is in the mantle mass: grid_0.01_-0.5_2.0 mantle S 1.1759556228481907e20 to 1.1759556228488423e20 (5.5e-13), grid_0.05_-1.5_0.48 mantle S 9.563373493908819e20 to 9.563373493912334e20 (3.7e-13), grid_0.05_-1.5_2.0 mantle S 1.1692686962694767e21 to 1.1692686962698916e21 (3.5e-13), grid_0.05_-0.5_2.0 mantle S 3.996533854793722e20 to 3.996533854795723e20 (5.0e-13). The other 28 stored cases are unchanged. The docs state the settle rule for the metal and the melt and that a pass leaving Fe below 1e-8 raises ValueError while a split that does not settle raises RuntimeError.
…s iron running out Sweep splits 359 (S), 2268 (C with trace S) and 849 (N and S) must raise RuntimeError after 2 x 1000 relaxed passes. 359 and 2268 have a root with Fe 1.5e-8 and 1.8e-5 by mass, but one float step of the metal composition there leaves a D residual of at least 8e-2 and 1e-5 at the returned metal, so no metal on the float grid meets the settle test; 849 has no root on a grid scan, which does not prove that none exists. The ValueError of the iron floor now says that a pass leaves iron below 1e-8, not that no root exists, and the docs state the RuntimeError limit for S- or C-rich, iron-poor inputs.
… move A plain pass that settles the metal on the start but meets a new melt H must take a second pass, and the grid rule must accept a move between 1 and 2 ulp after a 2-ulp secant step.
…e in the docstrings The docs name a damped return at the evaluated metal and give the share of round-off starts that reach the iron-poor root without the rule.
flags['metal']['root_rule'] is 1 when the damped passes replaced the first root (the S-C-H-N split with 2 roots) and 0 for a plain split, for the 32 stored splits and for an N-S split whose damped passes do not settle. The stored compare reads the stored keys of the metal flags.
flags['metal']['root_rule'] is 1 when the damped passes of a relaxed split past FeS replace the first root with one that holds more iron, else 0. The reference, explanation and tutorial pages name it.
Round-off decides which root of the S-C-H-N split the relaxed passes reach first; the flag is 1 exactly when that first root is the one with iron 6e-7.
… the evaluated metal The 2-root S-C-H-N split changed by up to 1e-6 settles on its iron-poor root and must take the root with iron 0.0497. A trace C in a return at the evaluated metal keeps g (small-silicate split 1544). The relaxed passes return core, mantle and w; the rule compares the damped root as another metal. Gap tests: iron out in the common passes, metal flags at the returned core, D residual with other oxides and NBO/T, the melt undamped in the damped passes, a metal without S.
… g for settled elements at the evaluated metal A relaxed split that settles on an iron-poor root past FeS now runs the damped passes too; their root replaces it only when it holds more iron and is another metal (an element apart by more than GRID_RTOL). In a return at the evaluated metal, elements that met the settle test keep g, so a trace element stays on its root. The relaxed passes return core, mantle and w. The 2 strata keep their outcomes; 14 converged splits move by at most 1.2e-5 relative, all C or H in a return at the evaluated metal, each with a smaller D residual. tools/core_split_two_roots.py gives the round-off starts of the 2-root split with and without the rule. The docs state the iron-out rule of the common passes, the cost and the residual bands.
… for bit A settled element of a return at the evaluated metal returns its core and mantle unchanged. Gap tests: one element equal and another apart, another metal with the same iron, a real split whose damped passes settle on the same metal (root_rule 0), the iron error of the common passes through pytest.raises.
…ated metal A settled element returns the core and mantle of the pass unchanged, and a bracketed element takes the rest of its bulk as mantle + core (g - w) / g, which does not cancel when the core dominates. The Reference states the iron error of the common passes without 'also'; the metal page gives the same-metal bound as the grid tolerance and the measured distance of the same-root damped settles.
…metal An element with g 0 that reaches a return at the evaluated metal has w 0 and is settled, so the settled set alone decides which elements keep the pass values.
…iron The common passes that settle after an iron-out send their return to the damped passes too. The same-metal test docstring gives split 406 its iron 0.09.
…lable laws report a point they cannot evaluate
The explanation of the core and mantle split now goes from the physics to the use to the numerics: what the split does, the mass balance with a worked example, why it iterates, the inputs, a table of the output flags, the solver steps in a collapsed box, and the sweep results as a table. The tutorial, the named-law reference and the 2 validation cases use shorter sentences and lists.
The How-to page walks through 2 complete examples: the Tsutsumi et al. (2025) H law as a partition file entry, checked against their Eq. 3 by hand, and the Suer et al. (2017) S law as a Python function, checked against the packaged entry. It lists the arguments of a function law and the errors a user meets. The tutorial, the reference and the explanation link to it.
… file The docstring says what each argument of a function law holds, with units, that Fe is the rest of x_metal, and that the function returns log10 D and raises ValueError where it cannot evaluate. The comments of partition_laws_core.toml say what each group of fields is for.
The code lines of the guide and the tutorial wrap before the edge of the code box, and the reference lists the 2 law files under the table of packaged names.
…s are optional The guide says in which passes a ValueError of a function law ends the split, that T_range_K and dIW_range are optional, that X_j is on the basis of the metal without H and C, and that laws also takes a loaded entry. The explanation says where nbo_t enters.
Example B of the guide imports what it uses and sets its own planet. The guide lists the error for a law of another element, and the comment of placeholder says what the flag does. The explanation says which quantities of the mass balance are inputs and names the planet of the 0.871 S share; the S validation case says a wrong term can move D_S by a factor of 10.
… of the melt The pages of the metal reservoir and the core split cite each paper with a footnote that gives the authors, title, journal, volume, pages, year, DOI and SciX link, as the other pages of the site do. The explanation says what the melt composition sets for a law on the mole basis, describes the solver without the names of its passes outside the collapsed steps, and names the iron-rich root rule where it points into those steps.
…-basis conversion The explanation calls one round of the split an update and uses that word before the solver steps; the mole-basis paragraph names the moles of metal and of melt per kg. The partition law tutorial cites Tsutsumi et al. at the sentence on their errors.
… a footnote The core and mantle split moves from the metal reservoir page to Explanations/core_split.md, with one subsection per input. The solver steps move to Reference/core_split_solver.md and the sweep of random inputs to a case on the metal validation page; the explanation keeps a short summary of each with a link. The long sections of the metal reservoir page get subsections. The validation index, the testing page and the core split pages cite each paper as Author (year) with a footnote, checked against the paper record, and the partition file reference names the data of the O law by its source table.
…ple notes to the tutorial The core split page starts with a short answer, defines the epsilon terms, the fixed point and NBO/T, uses pass for one round of the split, and lists what the code does outside each bound. The detail of the melt basis of the S law becomes a known limit of the S runs case. The validation page says which cases check the split and that the sweep is a recorded run; the index says the same. The notes on why the worked example uses this metal and these conditions, its O and Si case and the N and S flags move to the metal reservoir tutorial, where the case runs. Author and year are in the standard form on the touched pages.
6 of 8 tasks
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
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 adds
calliope.core_split.split_core_mantle, which splits the bulk H, C, N and S of a planet between the metal of the core and the silicate of the mantle at one equilibration point (default 60 GPa, 4000 K), with metal-silicate partition laws from the literature. PROTEUS does not call it yet; the coupling comes in a later PR.calliope.metal: packaged laws for S (Suer et al. 2017 Eq. 7, withdata/partition_laws_core.toml;register_lawadds a user law by name, as a file entry or a function. Each law carries its P, T andsplit_core_mantle(bulk_kg, M_metal, M_silicate, *, dIW, ...)returns(core_kg, mantle_kg, flags). D depends on the metal composition, so the split is a fixed point on the metal mass fractions and the melt H and C. Plain passes first; where the metal takes so much S that the passes swing, relaxed passes with a per-element step factor, an iron floor ofRuntimeErrorwhen it fails.flags['metal']['root_rule']). One constructed S-C-H-N input with very little silicate has 2 roots (Feflags[e]['D_residual']reports for each element how far the returned split is from its law;flags['metal']gives the S/Fe ratio andbeyond_FeS.tools/core_split_sweep.py(random splits in 2 strata) andtools/core_split_two_roots.pyreproduce the counts in the docs.No related issue.
Validation of changes
Both validation cases are on Validation/metal.md and run in the unit tier (
pytest tests/test_metal.py -k TestSuer2017andpytest tests/test_core_split.py -k "closed_form or ledger"):Physics impact: no existing output changes. Nothing on the solve path calls the new code, and PROTEUS does not call it;$\Delta\mathrm{IW}$ modes) and the branch compares them with
tests/data/reference_outputs.jsonholds outputs of main (solves without metal and with the placeholder metal, 2==.On the branch source at 248cc66 (macOS arm64, Python 3.12):
pytest testsat 27bb612, the last source change: 897 passed. 248cc66 adds 1 test case; that test passes, and it fails when the iron-out path skips the Fe-rich rule. The commits after it change docs, docstrings and comments only (metal.pyhas the same syntax tree without docstrings, and the law data file is unchanged), and the code examples of the how-to and the tutorial print what the pages show. Coverage withcoverage run --source=src/calliope: 99.4 % in total, 315 of 316 added source lines (the one miss ismetal.py:750).PYTHONPATH=src python tools/core_split_sweep.py run broadgives 2785 converged, 212ValueError, 3RuntimeError;run small-silicategives 2289, 204 and 7. With the Fe-rich rule switched off the counts are the same and no split changes (compare).PYTHONPATH=src python tools/core_split_two_roots.py: 30 of 30 round-off starts of the 2-root input return Fetime). The damped passes run on the 226 relaxed splits past FeS of the 5500 (110 broad, 116 small-silicate, counted from therunoutput) and replace none; a split whose damped passes do not settle costs up to 20000 more law evaluations. This matters once PROTEUS calls the split at every step.Not checked: a PROTEUS coupled run (PROTEUS does not call the split yet); the sweep counts, the 30 starts and the timings are from macOS arm64 only (the Linux nightly runs the 5-start test of the 2-root input).
Limits, stated in the docs: small-silicate split 1078 raises the iron$2.9 \times 10^{-2}$ ; 2 broad splits raise
ValueErroralthough a damped run reaches a root with FeRuntimeErroralthough a root with iron exists; for one moreRuntimeErrorinput no root was found, but a grid scan is not a proof that none exists. A registered function law must raiseValueErrorfor a point it cannot evaluate.Lines by category (
pr-breakdown 75faea5 3c586c0):The 4 plot files are the light and dark SVGs of the 2 validation figures above.
Checklist