Repository navigation
Add a metal reservoir that partitions H, C, N and S into the core - #43
Merged
Merged
Conversation
The file holds one metal-silicate partitioning law per element (H, C, N, S, O, Si) with its basis, ranges, covariance and source. The tests load a valid entry and the packaged placeholder set, and expect PartitionTableError naming the element and field for each malformed field: a TOML boolean or string as a number, NaN, inf, an unknown form, quantity or basis, a mole basis without its silicate basis, a reversed or non-finite range, a pressure unit other than GPa, an interaction term for an unknown element, and a covariance of the wrong size, not square, not symmetric or with a negative eigenvalue.
src/calliope/metal.py reads a TOML file with schema_version 1, a meta table and one partition law per element (H, C, N, S, O, Si), and checks every field: the form, quantity and metal and silicate bases, finite numeric coefficients, ordered ranges (an empty T or dIW range is unknown), interaction terms of known elements, a covariance that is square, symmetric and positive semidefinite, and the source and DOI. The packaged placeholder file holds Tagawa et al. (2021) Eq. 2 for H, Fischer et al. (2020) Eq. 1 with its covariance for C, lines through the Suer et al. (2023) Table 1 redox pairs for N (20 GPa, 3000 K) and S (60 GPa, 4000 K), and the Tagawa et al. (2021) SI Table 3 exchange coefficients for O and Si.
The placeholder H law gives the Tagawa et al. (2021) Eq. 2 values at 3 pressures and temperatures, with the -dIW/4 term, and lies more than 0.3 from their Eq. 10. The C law gives the Fischer et al. (2020) single-stage D_C of 0.550 with the O and S terms and 0.447 without them. A point outside a range, or on an axis whose range is unknown, is flagged, and a non-finite condition or a metal fraction of 1 is refused.
log10_D gives a + b/T + c P/T + d dIW_factor dIW + e NBO/T + sum eps_j log10(1 - X_j) for one entry, with P in GPa, T in K and the metal mole fractions of the interaction terms, and refuses a non-finite condition, a temperature at or below zero and a metal fraction outside [0, 1). out_of_range is true outside any range of the entry and on an axis whose range is unknown.
core_light_elements turns molar D on the H- and C-free metal basis into core mass fractions with a KLB-1 melt. The 3 single-stage cases of Tsutsumi et al. (2025) SI Table 1 (710 ppm H2O and 120 ppm C in the BSE, D_H and D_C as listed, core O, Si and S fixed) give their core H and C at the 2 decimals of the table; D times the melt mass fraction without the conversion misses the 40 GPa case, and O and Si in the metal lower the H and C uptake.
core_light_elements solves x prime = D x for H and C with the melt on the single-cation basis (H as HO0.5) and the metal mole fraction over the elements other than H and C, the rest of the metal being Fe; the solve is closed form. The default melt is the KLB-1 peridotite of Takahashi (1986). The test of the effect of O and Si in the metal now states the right direction: they are lighter than Fe, so a kilogram of metal holds more H- and C-free moles and, at a fixed mole ratio, more H and C. constants.py gains the molar masses of the rock-forming elements.
With melt and metal alone the core holds D x / (1 + D x) of an element, x the metal to melt mass ratio, for D from 0.1 to 100, and the law without the melt back-reaction (f = D x) lies outside the tolerance. No metal gives zero core masses; metal with no melt and a negative D are refused.
metal_masses gives the element masses in the equilibrated metal, D_i times the metal mass times the melt mass fraction of the element, with the melt mass M_mantle times Phi_global. No metal gives zero; metal with no melt and a negative or non-finite D raise ValueError.
On the metal of the closed-form solve, the mass-basis D of an entry on the H- and C-free metal basis times the melt H gives the closed-form core H, and differs from the molar D by more than 6 %; a mass-basis entry is its law; a metal without iron is refused.
mass_D evaluates an entry on the metal composition, with the metal mole fractions for its interaction terms, and converts a molar D to the metal over melt mass ratio with the moles per kg of metal (without H and C for the H- and C-free basis) over the moles of melt units per kg; a mass-basis entry is returned as it is.
With the same mass-basis D for H, C, N and S, no O or Si in the metal and dIW as a parameter, the solve with metal equals the solve without metal and a melt mass larger by D times the metal mass, for D of 1 and 30: the partial pressures agree, the core mass of each element is D times the metal mass times its melt mass fraction, and gas, melt and core close each budget. The gas-free core fraction misses the result, and a metal mass of 0 gives the solve without metal bit for bit.
With metal_mass above zero, equilibrium_atmosphere iterates on the metal composition: each pass evaluates the mass-basis D of H, C, N and S on the current metal and melt, solves the element balance with the core term D M_metal w(melt), and takes the new metal from the core masses, until no metal mass fraction moves by more than tol_w_core. O and Si in the metal are fixed fractions, and their O is a reservoir of the O budget in the authoritative-O residual; dIW of the equilibrium is a parameter or 2 log10(x_FeO / x_Fe). The result gains the core mass, metal mass fraction, D and range flag of each element and checks the closure over gas, melt and core. Without metal_mass the residuals are unchanged. The solve test compares the two runs from the same random state and with an acceptance tight enough for the identity at 1e-6.
exchange_metal gives the metal mole fractions of O and Si from the exchange K_D of FeO = Fe + O and SiO2 + 2 Fe = Si + 2 FeO with ideal activities: x_O = K_D 10^(dIW/2) and x_Si = K_D x_SiO2 10^(-dIW) on the KLB-1 melt. A half-reaction reading of the Si law falls far short, one dIW unit more oxidised scales x_O by 10^0.5 and x_Si by 0.1, and a metal with no iron left is refused.
exchange_metal evaluates the O and Si exchange K_D of the partition file and returns the metal mole fractions x_O = K_D 10^(dIW/2) and x_Si = K_D x_SiO2 10^(-dIW), with x_SiO2 on the single-cation basis of the melt and activity coefficients of 1; O and Si that leave no iron raise ValueError.
With O_mode_core and Si_mode_core set to partition, the converged metal holds the O and Si mole fractions of the exchange laws at its dIW, and the core O mass is its mass fraction times the metal mass.
With O_mode_core or Si_mode_core set to partition, each pass of the metal fixed point puts the exchange mole fraction of the element at the dIW of the pass into the metal, converted to a mass fraction on the current metal. The docstring states that the closure check uses the root tolerance of the solve and that <e>_kg_total excludes the core.
Each field the metal reservoir reads (metal_mass, P_eq, T_eq, dIW_eq, silicate_nbo_t, silicate_oxides, w_O_core, w_Si_core, tol_w_core) raises ValueError naming it for a boolean, a string, NaN, inf or a value outside its range; the melt composition needs known oxides, finite and non-negative, summing to 95-101 wt%.
check_metal_fields refuses a boolean, a string, NaN, inf or an out-of-range value in metal_mass, P_eq, T_eq, dIW_eq, silicate_nbo_t, w_O_core, w_Si_core and tol_w_core, and a melt composition with unknown oxides, a non-finite or negative oxide, or an oxide sum outside 95-101 wt%, naming the field. The metal path calls it before its first solve.
The placeholder H law, Tagawa et al. (2021) Eq. 2, puts each of the 6 runs of their SI Table 1 within its 1 sigma, while their Eq. 10 misses every run by more than 2 sigma. A validation plot pair shows the measured against the fitted log10 K_D, and proteus-mpl joins the develop extra for the plot fixture.
…cases The core term lies on D x / (1 + D x) over D from 0.1 to 100, with f = D x more than 10 % above it from D = 0.3 on, and the 3 Tsutsumi et al. (2025) single-stage cases give their SI Table 1 core H and C within the rounding of the table; each case has a validation plot pair.
The solve validation page gains the metal reservoir: the gas-free core fraction, the equal-D identity with a larger melt, the Tagawa et al. (2021) H runs and the Tsutsumi et al. (2025) single-stage cores, each with its reference, tolerance, test id and command, and the statement that <e>_kg_total excludes the core.
With dIW_eq_mode from_FeO the equilibrium dIW is 2 log10(x_FeO / x_Fe) on the KLB-1 melt and the converged metal, near IW-2.3; an unknown O, Si or dIW mode raises ValueError before the first solve.
equilibrium_atmosphere_authoritative_O has no metal fixed point, so a metal_mass above zero would be left out of its balance without notice; it now raises ValueError and names equilibrium_atmosphere.
With a solve stand-in, a metal H that alternates between 2 values raises RuntimeError after 30 passes, and a state whose gas, melt and core miss a budget raises the closure error.
…gases The authoritative-O solve must close the O budget over gas, melt and metal O, reproduce the larger-melt identity at equal D with its O budget raised by the O of the extra melt, and refuse a metal O above the O budget before the solve. A metal solve with He must keep the He pressure in the warm start of each pass. The closure check of the fixed point must cover O when the solve has an O budget.
…metal solve The O budget of each authoritative-O metal test is the gas and melt O of the metal solve at IW+2, plus the metal O, so the derived fO2 must return IW+2. The closure check of the fixed point is tested at the boundary of the acceptance gate of each solve: the largest H, C, N, S budget times rtol plus atol for the buffered solve, each budget times rtol with an O budget.
equilibrium_atmosphere_authoritative_O runs the metal fixed point when metal_mass > 0, with the metal O as a reservoir of the O budget. A metal O at or above the O budget is refused before the solve. The warm start of each pass carries the noble gas pressures and, with an O budget, the derived fO2. The closure check uses the acceptance gate of the solve and covers O when the solve has an O budget; the authoritative output carries the core masses of H, C, N and S.
The passes of the metal fixed point read the buffer from ddict, so fO2_model='fischer' given as an argument ran the passes on the default buffer.
The identity holds to relative 1e-6 on both buffers; the plotted pressures move with the buffer, the lowest one from 10^-10.45 to 10^-10.71 bar. The gas and melt O budget of the metal O case is 8.76e19 kg.
The Validation tab gets a Metal reservoir page after Solve with the 7 metal cases and their 4 plot pairs. The Solve page links to it, and the overview rows and the partition file page point to the new anchors.
dIW_eq and the from_FeO dIW enter only the partition laws and are never converted to an absolute fO2, so they do not depend on fO2_model; a test pins dIW_eq, the D and the exchange O on 2 buffers. The test rules select a buffer through fO2_model, since a patched default does not reach a solve.
calculate_mantle_mass serves standalone runs. Inside PROTEUS the mantle mass comes from the interior structure, with the Zalmoxis core mass the mass enclosed at the core-mantle boundary, and PROTEUS passes no metal_mass. The structure validation page, the code architecture page and the partition file page say so.
…t by name The buffer reaches a solved dIW_eq and the D through the metal and melt composition, so the full-solve comparison of 2 buffers gave way to a check of the dIW itself on 3 buffers at a fixed composition. The argument test compares the derived shift in dex and refuses a contradicting buffer on the metal path of both entry points.
The solved from_FeO dIW and a composition-dependent D change with the buffer through the solved composition. A patch of calliope.solve.DEFAULT_FO2_MODEL changes every solve. Inside PROTEUS the structure solve gives the core mass and the first mantle mass, and Aragog or SPIDER report the mantle mass after each interior step; the code architecture page holds this and the other pages link to it. The review checklists name the buffer path of a solve.
A solve at IW-4 with little melt and metal moves the solved dIW_eq and D_C with the buffer while every element closes. The parameter-mode dIW check compares with pytest.approx, so the float literal rule holds.
PROTEUS passes hf_row['M_mantle'], which its interior structure and energetics steps set, with the core mass from the structure solve; the code architecture page says so and the structure validation, partition file, API and coupling pages link to it or use the same words. A patch of __init__.__defaults__ changes only the direct calls of that class, and the buffer moves the D of each law with a dIW or composition term.
In parameter mode the D of a law with only a dIW term is the same on every buffer, and the D of a law with a composition term or a mole basis moves; the from_FeO dIW moves the D of every law with a dIW term. The coupling diagram shows M_mantle under the interior step, the defaults rule names the calls that omit the buffer, and the full-solve buffer test names its 2 buffers.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #43 +/- ##
==========================================
+ Coverage 98.90% 99.07% +0.17%
==========================================
Files 7 8 +1
Lines 913 1297 +384
Branches 167 262 +95
==========================================
+ Hits 903 1285 +382
- Misses 4 5 +1
- Partials 6 7 +1
Flags with carried forward coverage won't be shown. Click here to find out more.
🚀 New features to boost your workflow:
|
The page opens with what it checks and why, and each case says why the check exists, what it tests, how the right answer is known and why it matters, with a reading guide for each plot. The reference, tolerance, test id and command of each case sit unchanged in a collapsed Test details box, and the case headings, and so their anchors, stay as they were.
Each page starts with a short section on the part of the model it covers, what validation means here, where the right answer comes from and why it matters. The mantle-mass sentences on the structure page and the code architecture page use shorter plain sentences.
The overview says what validation means here, how the tab is organised by process, how to read the case table and how the case pages are laid out; the table rows stay as they were.
The equal-D cases name a wrong factor or a left-out element in the metal term, and the Tagawa case names the wrong fit of the paper (Eq. 10 for Eq. 2), each a fault that fails its test.
7 of 8 tasks
out_of_range, silicate_units_per_kg, cation_fraction and exchange_metal get NumPy Parameters and Returns sections that state the unit of each argument (GPa, K, log10 units relative to IW, oxide wt%) and the basis of each result (single-cation melt, x' metal without H and C). No code change.
New Explanations page on the metal reservoir: the terms, the core term of the element balance and where solve.py adds it, the partition laws and their basis, the 2 dIW modes and how dIW_eq relates to the IW buffer of the solve, the O and Si exchange with the silicate FeO, a worked example with the whole Earth core at 50 GPa and 3500 K, and the limits. The model overview, mass balance, oxygen fugacity, code architecture, testing, usage, configuration and both coupling pages link to it where a reader needs the metal. The coupling theory page no longer lists the metal keys as hf_row keys: PROTEUS does not pass metal_mass and does not copy them. The metal reservoir tutorial runs the solve without and with metal, the 2 dIW modes, O and Si in fixed and partition mode, the out-of-range flags and warnings, a scan of the core share against dIW_eq with its plot, and the authoritative-O solve with metal O. The partition law tutorial puts the H law of Tsutsumi et al. (2025, Eq. 3) into a copy of the packaged file, maps each coefficient, shows 2 format errors, compares D_H with the placeholder law, and lists what a law in a new form or for a new element needs in code, docs and tests. The output in both pages is the output of their own code.
Define every symbol of the balance and the law on the Explanations page, name the mass ratio r, keep "metal" for the reservoir and say that the keys call it core, explain the K_D form of the H law, the dissolved H and C in the melt moles, and the stop rule of a partitioned O or Si. The configuration page states that metal_mass 0 still checks the other metal keys, and the code architecture page counts 8 modules plus the generated _version.py. The tutorials print the budget next to the bulk inventory, use "metal share", say that no P_eq and T_eq put the placeholder N and S laws in range together, describe the log axis and the range of each law in the plot caption, separate the letters of the paper from the file fields, define the test terms of the developer checklist, and name the kg_core loops a new element needs. Usage Pattern 3 shows its output.
The Explanations page says "metal term" and "metal share" throughout, defines ddict, the buffered solve, ppmw and the out-of-range flag before the example table, explains why a partitioned O or Si moves by only part of its step, splits the long K_D and partition-mode paragraphs, and says which non-converging case raises ValueError and which raises RuntimeError. Its -2.48 now points to the tutorial step that shows it. The tutorial gives 3.56 as the net exchange, the partition law tutorial takes the packaged name from the table and puts the note on paper letters after the mapping table, and the developer checklist names the tuples that gate a new value and the metal test files. The configuration and reference pages say that partition_table is the one metal key not checked before the solve, and the code architecture page lists the data folder.
load_partition_table names tomllib.TOMLDecodeError for a file that is not valid TOML, log10_D and mass_D get Raises sections for their non-finite inputs, the X_j range, a metal with no iron, a D beyond the float range and the form check, and metal_masses names a negative or non-finite metal_mass. No code change.
The Explanations page says that a partitioned O or Si takes the whole exchange step unless the secant slope is negative or the step would use more than half of the iron left, gives the 1e-15 floor its purpose, states the K_D form of the H law as the equation the file stores, and words the single-stage example as an effect at its own conditions. The partition law tutorial puts the note on paper letters before the mapping table, defines K_D and NBO/T on the page, separates the value tuples from COEFFICIENTS, and says that a new form for an existing element needs a call from _with_metal. The step 7 timing reads a minute or more.
test_metal_diw_takes_no_absolute_fo2 checked dIW_eq in parameter mode only on the O'Neill buffer, so a _metal_diw that shifted it on another buffer passed. The check now runs on Hirschmann, Fischer and O'Neill.
The Raises section of exchange_metal names the PartitionTableError for a law that is not exchange_KD, and the ValueError for P, T or dIW that is not finite, a T that is not positive, and a K_D or exchange term beyond the float range. No code change.
The Explanations page and the reference state the window from a scan on the worked-example input: O and Si both partitioned converge at 5500 K and stop with ValueError from 6000 K; O alone converges up to 8480 K, ends at the 30-pass cap from 8490 to 9700 K and stops with ValueError from 9720 K, where K_D(O) x_FeO has passed 1/4. The K_D bullet splits off the H example, the step rule defines the secant slope and halves a step while it takes at least half of the iron left, the metal reservoir tutorial says to run the steps in order and gives the measured time of step 7, and the partition law tutorial evaluates the new law with log10_D and gives a new form for an existing element its own checklist item.
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 a metal (core) reservoir to the CALLIOPE solve, so H, C, N and S can partition between gas, melt and a metal that equilibrates with the melt. It is the CALLIOPE part of core partitioning; the PROTEUS coupling follows in its own PR.
ddict['metal_mass'][kg] above 0 adds a core termequilibrium_atmosphereandequilibrium_atmosphere_authoritative_O.src/calliope/metal.py, schema indocs/Reference/partition_file.md), evaluated atP_eq,T_eq, the dIW of the equilibrium and the melt and metal composition. A fixed point on the metal composition updates'fixed'mass fractions or'partition'through their exchange reactions with FeO (ideal activities). dIW is a parameter or follows the melt FeO,O_kg_totalstays gas plus melt in both fO2 sources andO_kg_coreis reported apart. CALLIOPE does not hold the silicate FeO; the new outputsn_FeO_exchange,FeO_kg_exchange,n_FeO_silicateandFeO_exchange_out_of_rangegive the caller what it needs to apply the exchange.'partition'mode and an exchange above 10 % of the melt FeO each log a warning once per process.metal_massthe result equals main (all 954 output values of 6 bufferedequilibrium_atmospheresolves at IW-3, 0 and +2 and melt fraction 1 and 0.4, the same seeded inputs on main97ee527and on this branch); withmetal_mass = 0the outputs of main are unchanged and the result also carries<e>_kg_core = 0for H, C, N and S (tests/test_metal_solve.py::test_zero_metal_mass_is_the_solve_without_metalcompares the full result in both solves). Withmetal_massabove 0, the new outputs are listed indocs/Reference/partition_file.md; the validation plots below show them against their references.fO2_model, from main) reaches every pass of the metal fixed point in both solves; a buffer argument that contradictsddict['fO2_model']raises. In the authoritative-O solve the buffer is now set before the input checks, so a call with both an unknown buffer and a badtarget_draises the buffer error first, asequilibrium_atmospherealready did.Docs and agent files in the same PR: an Explanations page for the metal solve (
docs/Explanations/metal.md), linked from the model, mass balance, oxygen fugacity, coupling, code architecture, testing, usage and configuration pages; two tutorials, one on using the metal reservoir (docs/Tutorials/metal_reservoir.md) and one on adding a partition law (docs/Tutorials/partition_law.md); a Validation page of its own for the metal reservoir (docs/Validation/metal.md); NumPy docstrings with units, ranges and Raises for the public functions ofmetal.py. Validation is a docs tab with an overview of every pinned case (tests/test_validation_index.pykeeps the overview equal to thereference_pinnedtests), and the agent rules move from.github/copilot-instructions.mdand the rule files under.github/intoAGENTS.mdandtests/AGENTS.mdwith the shared ecosystem blocks, checked bytools/agents/check_agents_md.pyin pre-commit and CI, as in MORS (FormingWorlds/MORS#65).Limits, stated in the docs: Si uses ideal activities (no interaction terms for Si); with dIW from FeO, O and Si both partitioned stop at high temperature with no iron left, and O alone ends at the pass cap where$K_D(\mathrm{O}),x_\mathrm{FeO}$ approaches 1/4 (
docs/Explanations/metal.mdgives the measured windows); the melt FeO and dIW do not follow the exchange inside one solve.Validation of changes
macOS (Apple silicon), Python 3.12, tip
61c20271e3eacd56e457162eb19a82938d8f0333(contains main97ee527):pytest testsin one process: 677 passed.python -m coverage runovertests, then the added lines ofgit diff 97ee527 -- srcagainstcoverage json: 384 of 385 lines run (99.7 %; the missing one is the unknowndIW_eq_moderaise thatcheck_metal_fieldsmakes unreachable); total 99.3 %.zensical buildat the tip gives the same 157 issues as main, page by page; both tutorials run as written and print the values on their pages.proteus-mpl >= 1.5.0joins thedevelopextra for the validation plots, and the placeholder TOML is package data.Validation cases (
docs/Validation/metal.md), each areference_pinnedtest with its plot:tests/test_metal_solve.py::test_metal_with_equal_d_is_a_larger_melt).tests/test_metal.py::test_tagawa_runs_lie_within_one_sigma_of_eq2).tests/test_metal.py::TestBasis::test_tsutsumi_single_stage_core,tests/test_metal.py::test_gas_free_core_fraction_curve). The authoritative-O solve with metal O returns the fO2 of the buffered solve, so the metal O stays out of the O budget (tests/test_metal_solve.py::test_authoritative_o_budget_excludes_the_metal_o).Checklist