Add multi-material mixing, H2 binodal suppression, and PALEOS unified EOS#56
Merged
timlichtenberg merged 115 commits intomainfrom Mar 25, 2026
Merged
Add multi-material mixing, H2 binodal suppression, and PALEOS unified EOS#56timlichtenberg merged 115 commits intomainfrom
timlichtenberg merged 115 commits intomainfrom
Conversation
The adiabatic temperature mode produced identical results to linear mode (issue #55) because the mass convergence break in the outer loop fired before the adiabat gate could activate on the next iteration. This commit fixes the convergence loop with a blending approach: when mass converges and temperature_mode='adiabatic', the solver continues iterating while gradually transitioning from linear T to adiabatic T (0% -> 50% -> 100% blend over 2 extra iterations). This prevents the solver from diverging when the temperature profile changes abruptly. Additionally, this adds PALEOS MgSiO3 as a new EOS option with phase-aware adiabatic gradient computation: - PALEOS tables (Zenodo 18924171) provide density and nabla_ad for both solid and liquid MgSiO3, covering 1 bar to 100 TPa - Phase-aware adiabat uses solid nabla_ad below solidus, liquid above liquidus, and melt-fraction-weighted average in the mushy zone - nabla_ad is converted to dT/dP = nabla_ad * T/P for integration Changes: - New PALEOS table loader with log-log RegularGridInterpolator - PALEOS density lookup via existing get_Tdep_density() phase routing - Phase-aware adiabatic gradient with solid/liquid/mixed nabla_ad - Convergence loop blend fix (applies to all T-dep EOS, not just PALEOS) - Material dictionaries extended to 5-tuple (index 4 = PALEOS) - PALEOS mass limit: 50 M_earth - Data provisioning from Zenodo 18924171 - Unit tests for PALEOS loading, density, nabla_ad, registration - Integration tests for PALEOS convergence and adiabatic vs linear
The PALEOS tables have NaN holes at the corners of the P-T domain where the thermodynamic solver did not converge. The liquid table is only 53% filled, with the low-P / high-T corner missing (e.g. T_max ~ 5200 K at P ~ 1 GPa). Previously, only global T bounds were enforced, so bilinear interpolation near the ragged domain boundary would return NaN (valid cell neighboring a NaN cell). This caused density lookup failures and solver divergence for planets above ~2.8 M_earth in adiabatic mode. Changes: - Precompute per-pressure-row valid T bounds during table loading - Clamp queries to the local valid T range before interpolation - Add NearestNDInterpolator fallback for remaining NaN results at the ragged boundary (bilinear interpolation failure) - Log a warning (once per table file) when per-cell clamping activates, informing the user that boundary values may be inaccurate With per-cell clamping, planets up to ~10 M_earth produce physical profiles in adiabatic mode instead of diverging, though convergence is not guaranteed above ~2.8 M_earth due to boundary density errors.
…plots Add analytic solidus/liquidus from Monteux et al. (2016, EPSL 448) Eqs. 10-13 with both A-chondritic and F-peridotitic composition models. Includes comparison plot against existing Zalmoxis tabulated melting curves. Also add diagnostic script for visualizing PALEOS EOS table coverage (density, nabla_ad) with melting curves and adiabatic profiles overlaid.
Introduces rock_solidus and rock_liquidus as config parameters in the
[EOS] section of the TOML file, replacing the hardcoded Monteux-600
tabulated melting curves.
Available options:
Solidus: Monteux16-solidus (analytic, default),
Monteux600-solidus-tabulated (legacy)
Liquidus: Monteux16-liquidus-A-chondritic (analytic, default),
Monteux16-liquidus-F-peridotitic (analytic),
Monteux600-liquidus-tabulated (legacy)
The analytic Monteux+2016 curves (Eqs. 10-13) are defined for all
pressures, eliminating the NaN fill_value issue of the tabulated curves
at pressures outside the table range. This ensures correct phase routing
(solid/mixed/liquid) in get_Tdep_density() at all conditions.
Changes:
- New src/zalmoxis/melting_curves.py module with analytic solidus/liquidus
functions and a dispatcher for config string lookup
- Updated get_solidus_liquidus_functions() to accept curve identifiers
- Config parsing in load_zalmoxis_config() reads rock_solidus/rock_liquidus
- Backward compatible: old TOML files without these fields use the
Monteux16 analytic defaults
The Brent pressure solver explores trial center pressures during bracket search. With the previous T(r) parameterization, the temperature at each radius was fixed regardless of the ODE's actual pressure. At trial P_center values far from the true solution, this created unphysical (low P, high T) combinations at interior radii, hitting NaN gaps in the PALEOS tables and biasing the Brent solution (~14% mass error at 3 M_earth). The fix: parameterize the adiabatic temperature as T(P) instead of T(r). The adiabat is naturally defined in P-space (dT/dP = nabla_ad * T/P), so T(P) always returns thermodynamically consistent temperatures. During the Brent search, the temperature now tracks the ODE's actual pressure via log-P interpolation of the previous iteration's adiabat. The temperature_function signature changes from f(r) to f(r, P): - Adiabatic mode: uses P argument (log-P interpolation of T(P) adiabat) - Blended (0.5): T = (1-blend)*T_linear(r) + blend*T_adiabat(P) - Linear/isothermal/prescribed: ignores P, uses r only Results with PALEOS:MgSiO3 adiabatic mode: - 1 M_earth: converged, R=1.099 R_E, T_c=5262 K (unchanged) - 3 M_earth: converged, R=1.482 R_E, T_c=6994 K (was 14% mass error) - 5 M_earth: converged, R=1.669 R_E, T_c=8348 K (was diverging) - 10 M_earth: converged, R=1.944 R_E, T_c=10551 K (was hanging)
Tests: - test_melting_curves.py: 16 unit tests for Monteux16 analytic solidus/liquidus, melting curve dispatcher, per-cell clamping, and NN fallback - test_convergence_PALEOS.py: expanded with parametrized 3/5/10 M_earth adiabatic convergence regression tests (T(P) parameterization) Documentation: - configuration.md: add PALEOS:MgSiO3 EOS option, melting curve selection (rock_solidus, rock_liquidus), PALEOS config example, mass limit warning - model.md: add PALEOS EOS physics section with nabla_ad description, phase-aware adiabat, grid coverage, T(P) parameterization - data.md: add PALEOS table files, fix mass-radius curve inventory (removed 3 non-existent Zeng curve entries) - Add API reference for zalmoxis.melting_curves module
Implements the silicate melting curves from Stixrude (2014, Phil. Trans.
R. Soc. A 372, 20130076):
- Eq. 1.9: liquidus T_rock = 5400 K * (P / 140 GPa)^0.480
Simon-like power law fitted to Lindemann melting theory.
- Eq. 1.10: solidus = liquidus * (1 - ln(x0))^{-1} with x0 = 0.79
Cryoscopic depression for Earth-like mantle composition.
The Stixrude curves have two advantages over Monteux+2016:
1. Valid at all pressures (single power law, no piecewise transitions)
2. Liquidus > solidus guaranteed at all P (constant multiplicative
factor), whereas the Monteux piecewise parameterization crosses at
~500 GPa
New config identifiers: Stixrude14-solidus, Stixrude14-liquidus
Default changed from Monteux16 to Stixrude14 in all function signatures,
config defaults, and default.toml.
…x16 melting curves
Fixes across 6 documentation files: - api/index.md, code_architecture.md: add melting_curves.py to module trees - code_architecture.md: add PALEOS-2phase to all T-dep EOS lists, update data flow diagram default to PALEOS, document T(P) parameterization and phase-aware nabla_ad in compute_adiabatic_temperature() - configuration.md: add Stixrude14 solidus/liquidus to melting curves table, fix surface_temperature 3500->3000 in PALEOS example, fix adaptive_radial_fraction scope (all T-dep EOS, not just WB2018), add PALEOS to mass limit upgrade recommendation - usage.md: add PALEOS-2phase to EOS table, T-dep EOS lists, and plots documentation; mention adiabatic mode availability - model.md: update default mantle EOS example to PALEOS-2phase, add PALEOS to higher-mass recommendation - data.md: add full melting curves inventory (tabulated + analytic) with Stixrude14 entries and cross-reference to config guide
- Fix 3 broken cross-references to #pressure-solver-brents-method (anchor lives in process_flow.md, not model.md or code_architecture.md) - Add melting_curves module to mkdocs.yml navigation - Fix naming convention: PALEOS -> PALEOS-2phase in configuration.md - Update model.md introduction: "three families" -> enumerate all families - Update usage.md quick-start: describe actual default config (PALEOS-2phase mantle, 3000 K surface temperature, data files required)
Physics reviewer findings: - Add NaN guard in _paleos_clamp_temperature for all-NaN pressure rows (returns unclamped, lets NN fallback handle it explicitly) - Add forward Euler comment in adiabat integration loop documenting first-order accuracy limitation near phase boundaries Config/melting curves reviewer findings: - Fix integration tests to pass rock_solidus/rock_liquidus from config (was silently using code defaults instead of default.toml values) - Fix tautological blend step test (now inspects actual source code) - Fix stixrude14_solidus docstring: 4370 K at 140 GPa, not 4100 K - Add F-peridotitic liquidus > solidus test (up to 140 GPa)
Refactor material_dictionaries from a fragile 5-tuple indexed by position to a flat dict (EOS_REGISTRY) keyed by EOS identifier strings. This eliminates hard-coded index access and makes adding new EOS trivial. Add support for unified PALEOS tables (Zenodo 19000316) that contain all stable phases for each material in a single file: - PALEOS:iron (5 phases: alpha-bcc, delta-bcc, gamma-fcc, epsilon-hcp, liquid) - PALEOS:MgSiO3 (6 phases: 3 pyroxene, bridgmanite, postperovskite, liquid) - PALEOS:H2O (7 EOS: ice Ih-X, liquid, vapor, superionic) Key changes: - eos_properties.py: new EOS_REGISTRY dict, keep legacy dicts for compat - constants.py: add PALEOS:iron, PALEOS:MgSiO3, PALEOS:H2O to TDEP_EOS_NAMES - eos_functions.py: add load_paleos_unified_table() for single-file tables, get_paleos_unified_density() with mushy zone support, refactor calculate_density() to use dict dispatch, update compute_adiabatic_temperature() for unified tables (nabla_ad directly from table, T-dependent core) - zalmoxis.py: load_material_dictionaries() returns dict, add mushy_zone_factor config parameter, update mass limits and melting curve loading - structure_model.py: update docstrings for dict type - setup_utils.py: add downloads for 3 new Zenodo tables - default.toml: change defaults to PALEOS:iron + PALEOS:MgSiO3, add mushy_zone_factor parameter - test_adiabatic.py: update for dict-based material_dicts
Add PALEOS:iron, PALEOS:MgSiO3, PALEOS:H2O to the configuration guide (EOS options table, temperature profiles, mass limits, examples) and to the data reference (new Zenodo 19000316 files). Document mushy_zone_factor parameter for configurable mushy zone width.
New test file test_paleos_unified.py covers: - Registration of PALEOS:iron, PALEOS:MgSiO3, PALEOS:H2O in lookup tables - Dict-based material_dictionaries (return type, backward compat) - Unified table loader (grid structure, liquidus extraction) - Unified density lookup (sharp boundary and mushy zone modes) - Unified nabla_ad lookup - Mass limit enforcement - mushy_zone_factor config parsing - calculate_density() dict dispatch for all EOS types New validation plot script plot_paleos_full_validation.py generates: - T-P profiles overlaid on PALEOS table phase maps - Density-pressure profiles (PALEOS vs Seager2007) - Mass-radius comparison (PALEOS adiabatic vs Zeng+2019) - Radial profiles for 1 M_earth (linear vs adiabatic) - mushy_zone_factor comparison (1.0 vs 0.8)
- model.md: add unified PALEOS physics section with phase boundary extraction, mushy zone, T-dependent iron core; update EOS identifier table and validity ranges - code_architecture.md: update module description, data flow diagram, key function docs for dict-based dispatch and unified PALEOS adiabat - usage.md: update default config description, EOS table, temperature mode docs - installation.md: update disk space estimate, convergence troubleshooting - zalmoxis.eos_properties.md: rewrite API reference to document EOS_REGISTRY structure with all EOS families
- Fix post_processing() crash with unified PALEOS: uses_Tdep_mantle triggered get_Tdep_material() which requires external melting curves, but load_solidus_liquidus_functions returns None for unified PALEOS. Changed to uses_phase_detection using _NEEDS_MELTING_CURVES set. - Fix mat_PALEOS_2ph unbound local: if PALEOS-2phase:MgSiO3 was used on any layer (not just mantle), the variable was never assigned. Now checks tdep_eos_used instead of just mantle_eos. - Strip whitespace from phase strings in load_paleos_unified_table() to prevent silent liquidus extraction failure from trailing whitespace. - Add bounds check for liquidus extrapolation: when query pressure is outside the liquidus coverage range, fall back to direct table lookup instead of using np.interp flat extrapolation. - Document known limitation: mushy_zone_factor weighting is not applied to nabla_ad in the adiabat computation (exact for default factor=1.0).
New validate_config() function in zalmoxis.py checks all input parameters for physical and logical consistency before the solver runs. Called automatically at the end of load_zalmoxis_config(). Validates: - Planet mass: must be positive - Mass fractions: ranges, sum <= 1, 3-layer requires mantle_mass_fraction > 0 - Temperature: positive, valid mode, adiabatic requires T-dependent EOS - Mushy zone factor: [0.7, 1.0] range, requires unified PALEOS EOS - Numerical: num_layers [10, 10000], positive tolerances and step sizes - Pressure solver: non-negative target, positive tolerance - EOS cross-checks: melting curves required for WB2018/RTPress/PALEOS-2phase All errors raise ValueError with clear messages explaining the problem and how to fix it. This prevents the solver from silently producing incorrect results with invalid parameter combinations. 26 new unit tests in test_config_validation.py cover all validation paths.
New benchmark_paleos_eos.py generates comparison plots and runtime documentation across all EOS configurations: Configurations tested: - PALEOS:iron + PALEOS:MgSiO3 (unified, adiabatic) - Seager2007:iron + PALEOS-2phase:MgSiO3 (2-phase, adiabatic) - Seager2007:iron + WolfBower2018:MgSiO3 (WB2018, adiabatic, <= 7 M_E) - Seager2007:iron + Seager2007:MgSiO3 (cold reference, isothermal) - Mushy zone factor variations (1.0, 0.9, 0.8, 0.7) - 3-layer model with PALEOS:H2O ice layer Plots generated: 1. T-P phase diagram with solidus/liquidus overlay 2. Density-pressure comparison across EOS 3. Radial profiles (T, rho, P, g) multi-panel 4. Mass-radius diagram vs Zeng+2019 5. Mushy zone factor effect on profiles 6. Phase regime visualization along adiabat 7. Runtime comparison bar chart Output: output_files/benchmark_paleos/ (plots + runtimes.txt)
…hmark Fix: skip adiabat integration step when pressure is below 1 bar (1e5 Pa). At very low P, dT/dP = nabla_ad * T/P diverges because T/P becomes huge. This caused runaway temperatures in 3-layer models where the surface pressure can be far below 1 bar after the terminal event pads with zeros. The fix holds T constant through the low-pressure shell, matching the physical expectation (the adiabat is meaningless in the near-vacuum). This enables 3-layer PALEOS:iron + PALEOS:MgSiO3 + PALEOS:H2O models in adiabatic mode for the first time. New benchmark_phase_regimes.py tests phase transitions across temperatures: - Temperature sweep 300-5000 K at 1 and 5 M_earth (2-layer) - Phase identification along the adiabat (iron + MgSiO3) - Mushy zone sweep at T_surf=1500 K (where adiabat crosses liquidus) - 3-layer PALEOS models at 300-1000 K - EOS comparison (PALEOS vs PALEOS-2phase vs WB2018 vs Seager) at 1500 K
New plot_verification.py generates 6 diagnostic plots: 1. T-P phase diagram: all 3 PALEOS tables (Fe, MgSiO3, H2O) with density background, extracted liquidus, Monteux16 melting curves, and adiabats at 300-5000 K overlaid. 3-layer H2O adiabat on H2O panel. 2. Radial profiles colored by phase: 5 rows (300, 1000, 2000, 3000, 5000 K) x 4 columns (T, rho, P, g) with each line segment colored by the thermodynamically stable phase (epsilon-hcp, bridgmanite, ppv, liquid). 3. Phase transition map: T_surf vs radius scatter for iron core and MgSiO3 mantle, showing where phase transitions occur as surface T increases. 4. 3-layer model: adiabatic and isothermal profiles for iron + MgSiO3 + H2O with CMB and mantle/ice boundary markers, colored by phase per layer. 5. EOS comparison: PALEOS unified vs PALEOS-2phase vs WB2018 vs Seager at 1000, 2000, 3000 K surface temperature (4 panels x 3 rows). 6. Summary: R/R_E, T_center, rho_center vs T_surf for 1 and 5 M_earth.
New module mixing.py with LayerMixture dataclass, parser, and mixing functions. Each layer can now contain multiple EOS materials mixed by volume additivity (harmonic mean density). Config format: "PALEOS:MgSiO3:0.85+PALEOS:H2O:0.15" (+ separated components with mass fractions). Single-material strings are backward compatible: "PALEOS:iron" parses to LayerMixture(["PALEOS:iron"], [1.0]). Key changes: - mixing.py: LayerMixture dataclass with runtime-updatable fractions (for PROTEUS/CALLIOPE coupling), parse_layer_components(), calculate_mixed_density(), get_mixed_nabla_ad() - structure_model.py: get_layer_eos() returns LayerMixture, coupled_odes() uses calculate_mixed_density(), mushy_zone_factor now threaded through (fixes latent bug where ODE used default 1.0) - eos_functions.py: compute_adiabatic_temperature() accepts layer_mixtures, delegates nabla_ad to get_mixed_nabla_ad() for multi-component layers - zalmoxis.py: main() accepts optional layer_mixtures, validation and melting curve checks are component-aware - 27 new unit tests in test_mixing.py (parsing, harmonic mean, runtime update, Tdep detection, backward compat) All 145 unit tests pass.
validate_config() now checks: - Mass fraction sums in multi-material layers - Warns about mixing T-dependent with T-independent EOS in same layer - Warns about high H2O fraction (>30%) at high T (>2000 K) - Warns if core uses non-iron EOS Documentation updated: - configuration.md: new "Multi-material mixing" section with syntax, rules, mixing physics, warnings, and hydrated mantle example - usage.md: mention multi-material mixing format - model.md: add mixing physics (harmonic mean density, weighted nabla_ad, PROTEUS/CALLIOPE runtime fraction update)
Welcome to Codecov 🎉Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests. Thanks for integrating Codecov - We've got you covered ☂️ |
timlichtenberg
requested review from
Copilot,
ema-jungova,
maraattia and
planetmariana
There was a problem hiding this comment.
Pull request overview
This PR adds multi-material layer mixing (with phase-aware suppression), hydrogen miscibility/binodal suppression in mixtures, and a unified PALEOS EOS registry/table format, alongside new grid-running/plotting tools and expanded docs/tests.
Changes:
- Replaced per-layer single-EOS handling with layer mixtures and mixed-density evaluation in the structure ODE.
- Introduced unified EOS registry entries for PALEOS (and Chabrier H) and expanded data download tooling.
- Added parameter-grid runner/plotter utilities and substantial documentation + test suite expansion.
Reviewed changes
Copilot reviewed 62 out of 82 changed files in this pull request and generated 9 comments.
Show a summary per file
| File | Description |
|---|---|
| src/zalmoxis/structure_model.py | Switches ODE density closure to mixture-based density; updates temperature handling in RHS wiring |
| src/zalmoxis/eos_properties.py | Introduces EOS_REGISTRY keyed by EOS identifier strings; adds PALEOS unified + Chabrier H entries |
| src/zalmoxis/constants.py | Expands T-dependent EOS set and adds defaults for phase-aware mixing suppression |
| src/tools/setup_utils.py | Adds Zenodo tarball download/extract helper for EOS datasets |
| src/tools/setup_tests.py | Adds a PALEOS run helper for tooling-based validation |
| src/tools/run_grid.py | New utility to sweep TOML parameters, run in parallel, emit CSV/JSON summaries |
| src/tools/plot_grid.py | New utility to plot grid_summary.csv outputs with grouping options |
| src/tools/plot_paleos_tables.py | New diagnostics for PALEOS table coverage and adiabats (tool script) |
| src/tools/plot_paleos_full_validation.py | New end-to-end PALEOS validation plotting script (tool script) |
| src/tools/plot_paleos_adiabat.py | New PALEOS adiabat comparison plotting script (tool script) |
| src/tools/melting_curves_monteux2016.py | Adds analytic Monteux+2016 melting curves utility + plotting |
| src/tests/test_plot_profiles.py | Adds smoke tests for profile plotting output paths and phase-legend branch |
| src/tests/test_paleos.py | Adds unit tests for PALEOS-2phase table loading + interpolation behavior |
| src/tests/test_melting_curves.py | Adds melting-curve dispatcher and PALEOS clamping/NN fallback tests |
| src/tests/test_convergence_PALEOS.py | Adds PALEOS convergence integration tests for linear/adiabatic modes |
| src/tests/test_adiabatic.py | Updates adiabatic tests for mixture-aware interface and adds PALEOS adiabat cases |
| pyproject.toml | Updates docs extras (Zensical) |
| mkdocs.yml | Adds new docs pages to nav (mixing/binodal/EOS physics/grids + API refs) |
| input/default.toml | Updates defaults and documents new mixing/binodal/mushy-zone parameters |
| input/grids/mass_radius.toml | Adds example grid spec for mass-radius sweep |
| input/grids/mass_temperature.toml | Adds example grid spec for mass-temperature sweep |
| input/grids/h2_mixing.toml | Adds example grid spec for H2 mixing sweep |
| input/grids/h2o_mixing.toml | Adds example grid spec for H2O mixing sweep |
| docs/index.md | Updates landing page (coverage badge + getting-started pointer) |
| docs/Reference/data.md | Adds PALEOS unified + Chabrier table documentation and melting-curve identifiers |
| docs/Reference/api/index.md | Updates API index for new modules |
| docs/Reference/api/zalmoxis.eos_properties.md | Documents EOS_REGISTRY and PALEOS/Chabrier registry entries |
| docs/Reference/api/zalmoxis.mixing.md | Adds mixing API doc page |
| docs/Reference/api/zalmoxis.melting_curves.md | Adds melting-curves API doc page |
| docs/Reference/api/zalmoxis.binodal.md | Adds binodal API doc page |
| docs/How-to/installation.md | Updates disk space and adds “first simulation” step |
| docs/How-to/usage.md | Refactors usage docs; adds explanation of outputs and points to grids |
| docs/How-to/testing.md | Updates testing docs for expanded suite + coverage badge |
| docs/How-to/grids.md | New how-to guide for running/plotting parameter grids |
| docs/Explanations/process_flow.md | Updates process-flow doc for mixing suppression + wording tweaks |
| docs/Explanations/eos_physics.md | New EOS physics explanation page (PALEOS, Chabrier H, legacy EOS) |
| docs/Explanations/mixing.md | New explanation page for mixing + suppression + binodal integration |
| docs/Explanations/binodal.md | New explanation page for binodal/miscibility suppression models |
| docs/Explanations/code_architecture.md | Updates architecture doc for mixing/binodal/melting-curves modules |
| README.md | Updates README with features, quick start, docs links |
| CONTRIBUTING.md | Updates docs build instructions to Zensical |
| CLAUDE.md | Adds repo-specific agent/developer guidelines |
| .github/workflows/CI.yml | Adds Codecov upload and renames workflow |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Member
Author
|
@maraattia @ema-jungova @planetmariana This is ready for review. I cannot join the PROTEUS dev meeting on Monday, so let's please iterate here. My suggestion is to mainly check out the new docs; which holds most explanations. I tested the main code changes extensively, so please just check if the main docs can be followed, and Ema will do her MSc structure parameter sweep with this then; if there are further issues another PR can solve it. I would like to merge by Tuesday next week. |
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 3b6eacef44
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
timlichtenberg
added
documentation
maraattia
previously approved these changes
Mar 24, 2026
maraattia
left a comment
maraattia
left a comment
There was a problem hiding this comment.
This is a great addition to PROTEUS, implementing state-of-the-art features. Besides the minor points I flagged in comments, this should be ready for merging.
A question about the usage though: the mass non-conservation could be a problem if the user injects pathological amounts of volatiles (if it's suppressed from a structure perspective, the mass is just lost). Shouldn't the percentage of added volatiles into a condensed layer be capped at input stage?
Docs:
- Fix LaTeX subscript formatting in binodal.md: T_b -> T_\mathrm{b},
T_c -> T_\mathrm{c}, x_c -> x_\mathrm{c}, P_c -> P_\mathrm{c}
- Document PALEOS analytic liquidus (Belonoshko+2005 / Fei+2021) in
eos_physics.md, clarify that both analytic and table-extracted paths exist
- Expand mush zone docs to cover both density (specific volume weighting)
and nabla_ad treatment in all three phase regimes
- Add nabla_ad mixing limitation note citing Kempton+2023 App. B,
noting the different thermodynamic context (condensed vs gas mixtures)
- Strengthen mass non-conservation warning for large volatile fractions
- Add H2O-MgSiO3 binodal to future extensions list
- Fix x_c and T_c formatting in mixing.md
Plots:
- Regenerate sigmoid_suppression.png with systematic annotation format:
component phase (P, T) on top, rho = value below
- Regenerate suppression_weights.png with corrected Gas-like/Condensed
label placement, rho_min and T_b arrows, sigma labels in upper left
- Add generation scripts for both doc figures
- Make profile plot phase labels material-specific (Liquid Fe, Liquid
MgSiO3, etc.) instead of generic "Liquid"
Phase panel now shows "Liquid Fe", "Bridgmanite", "Liquid MgSiO3" instead of generic "Liquid" for all materials.
Extend P range from 1 bar to 570 GPa so both transitions are visible: density sigmoid at ~0.01 GPa and binodal at ~5-10 GPa. Uses actual Zalmoxis profile data for the mantle (3 GPa+) with isothermal extrapolation at lower pressures.
Use a cooler adiabat (T_surf=1500K) to separate the density transition (~0.1 GPa) from the binodal transition (~10 GPa) by two decades in pressure. Shade three regimes: H2 excluded, dense but immiscible, and H2 included. Both sigma curves now clearly visible and well separated.
Member
Author
|
Thanks @maraattia for the comments, very useful! I think the comments are reasonably addressed for this time; of course we will continue to update Zalmoxis in the next weeks further. I will merge sometime tomorrow. |
Clamp solid-side and liquid-side evaluation temperatures against PALEOS's own melting curve before looking up thermodynamic properties in the unified table. This prevents cross-phase lookups when external solidus/liquidus parameterizations (e.g. from PROTEUS #658) disagree with PALEOS's internal phase boundary at high pressures. The guard is a 1 K offset that is effectively a no-op for the current code (T_liq shifts from T_melt to T_melt+1 K, T_sol is untouched), but activates defensively if future external melting curves cross the PALEOS boundary. A warning is emitted only when a genuine cross-phase incursion is detected. Suggested by @maraattia in PR #56 review.
github-project-automation
bot
moved this from In Progress
to Done
in PROTEUS Development Roadmap
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
Closes #24, closes #55. Partially addresses #47 (H2-silicate and H2-H2O binodal models implemented; Fe-H2 binodal model (Young+2025) and He EOS not yet added).
This PR adds three major capabilities to Zalmoxis: multi-material mixing in planetary layers, hydrogen miscibility physics (binodal suppression), and a unified PALEOS equation-of-state format. It also substantially expands the test suite, documentation, and tooling.
What changed and why
1. PALEOS unified EOS (iron, MgSiO3, H2O, Chabrier H)
Previously, each EOS material (iron, silicate, water) used a different code path with different table formats. Now all PALEOS tables use a single unified reader that handles density, phase state, and adiabatic gradient lookups from one 10-column table per material. This also adds Chabrier H2/He tables in the same format, enabling hydrogen envelopes for sub-Neptunes.
The profile plot below shows the solver output for a 1 Earth-mass planet with a PALEOS iron core and MgSiO3 mantle. The six panels show how density, pressure, temperature, gravity, phase state, and enclosed mass vary from the center to the surface:
2. Multi-material mixing
Planetary layers can now contain mixtures of materials (e.g., 90% MgSiO3 + 10% H2 by mass). The mixed density is computed using a volume-additive harmonic mean, weighted by mass fractions. A sigmoid suppression function removes non-condensed volatile components (like H2 gas above its critical point) from the density calculation so they do not artificially inflate the planet:
3. H2 binodal suppression
When hydrogen is dissolved in a silicate or water mantle, it only contributes to the bulk density below its miscibility limit (the binodal curve). Above the binodal temperature, H2 and silicate/water are fully miscible and form a single phase. Below it, they separate into H2-rich and H2-poor phases.
This PR implements two binodal models from the literature:
H2-MgSiO3 (Rogers et al. 2025):
H2-H2O (Gupta et al. 2025):
4. Performance improvements
math.expinstead ofnp.expin tight loops (binodal, sigmoid)5. Parameter grid runner and plotter
New tools for sweeping any combination of input parameters and visualizing results:
Mass-radius curve (single parameter sweep):
H2 mixing (mass x composition x temperature, 24 models). Each line is a unique combination of mantle composition and surface temperature, showing how dissolved hydrogen inflates the radius:
H2O mixing (mass x water fraction x temperature, 24 models):
Temperature effect (mass x surface temperature, 20 models). Higher surface temperatures lead to larger radii through thermal expansion:
6. Test coverage expansion
Unit tests: 174 -> 425. Total: 217 -> 468. New test files cover EOS functions, core solver helpers, binodal suppression, mixing batch paths, and profile plotting.
7. Documentation overhaul
Validation of changes
Checklist
Relevant people
@maraattia @ema-jungova @planetmariana