Add external mesh input and EOS out-of-range safety#5
Conversation
Enable SPIDER to read P(r), rho(r), g(r) from an external file at initialization, replacing the Adams-Williamson EOS for mesh construction. This allows Zalmoxis (or any hydrostatic solver) to provide the static structure for super-Earths to sub-Neptunes. New options: -MESH_SOURCE 1 (0 = AW default, 1 = external file) -mesh_external_filename <path> The AW pathway is unchanged and remains the default. When MESH_SOURCE=1, the external file provides radius, pressure, density, and gravity at both basic and staggered nodes. SPIDER nondimensionalizes the input, computes shell masses, and derives dxi/dr from the density profile. Includes test scripts: - generate_aw_mesh.py: AW round-trip verification - generate_super_earth_mesh.py: Non-AW profile smoke test - validate_mesh_fields.py: Mesh field consistency checker
The external mesh pathway uses midpoint-rule quadrature for shell masses, while AW uses analytical integrals. For 50 nodes this gives ~2e-5 discretization error in dxidr that propagates through the time integration to ~1.5e-3 by t=1200 years. Set rtol=2e-3 accordingly. Also rename super-Earth test to non-AW mesh test using Earth radius to stay within EOS table parameter ranges.
The script assumed SPIDER JSON had a "mesh" section with structured
sub-objects, but mesh fields are actually in the flat "data" section
with names like "radius_b", "pressure_b", each as {"scaling": X,
"values": [...]}. The fallback path compared rho*g against itself,
passing trivially regardless of actual correctness.
Rewrite load_spider_json() to extract fields from the correct JSON
paths: data["data"] for mesh fields and data["atmosphere"] for
mantle_mass. All six validation checks now test real values.
At high CMB pressures (>430 GPa), the solid-phase P-S EOS tables can be evaluated beyond their valid entropy range. The clamped lookups produce negative thermal expansion coefficients, which cause NaN in the mixing-length convective transport (sqrt of negative) and crash CVode. Changes: - interp.c: emit first-occurrence warnings when P or S is clamped to table edges (one per table per direction, avoids log spam) - interp.h: add warning counters to Interp2d struct, remove const from SetInterp2dValue signature to allow counter mutation - matprop.c: clamp alpha < 0 to zero in GetEddyDiffusivity() with a one-time warning; guard the crit/alpha division against zero
Add section 4 to README documenting -MESH_SOURCE and -mesh_external_filename options, file format, and EOS out-of-range warnings. Add plot_validation_F.py for generating PR plots from Block F test outputs.
![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: a138d39eeb
ℹ️ 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".
There was a problem hiding this comment.
Pull request overview
Adds support for running SPIDER with a pre-computed radial mesh provided via an external file (instead of internal Adams–Williamson mesh generation), and introduces safety mechanisms to prevent silent NaN propagation when EOS lookup tables are queried out-of-range.
Changes:
- Add
-MESH_SOURCE/-mesh_external_filenameoptions and route mesh initialization through a new external-mesh reader when enabled. - Add EOS-table out-of-range one-time warnings and clamp negative thermal expansion coefficient (
alpha) to prevent NaNs in mixing-length convection logic. - Add/extend test coverage and utilities (mesh generators, validation scripts) plus README documentation for the external mesh format.
Reviewed changes
Copilot reviewed 12 out of 14 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
parameters.h |
Adds new parameters for choosing mesh source and external filename. |
parameters.c |
Parses new mesh options and conditionally creates/destroys the AW EOS used for mesh mapping. |
mesh.c |
Implements SetMeshFromExternalFile() and integrates it into mesh setup flow. |
interp.h |
Extends Interp2d with warning counters; updates SetInterp2dValue signature to allow mutation. |
interp.c |
Initializes per-table warning counters and emits first-occurrence out-of-range warnings with clamping. |
matprop.c |
Clamps negative alpha and guards divisions/sqrt usage to avoid NaN propagation. |
README.md |
Documents external mesh options, file format, and new EOS out-of-range behavior. |
tests/tests.yml |
Adds new SciATH test cases for external mesh round-trip and non-AW mesh runs. |
tests/generate_aw_mesh.py |
Adds a Python AW mesh generator used by tests (currently requires SciPy). |
tests/generate_super_earth_mesh.py |
Adds a synthetic non-AW mesh generator for external-mesh pathway testing. |
tests/validate_mesh_fields.py |
Adds a field-level validator comparing JSON mesh fields to external mesh input. |
tests/plot_validation_F.py |
Adds a plotting helper for PR validation outputs (not part of CI test execution). |
tests/data/.gitkeep |
Keeps tests/data/ present for generated mesh files. |
.gitignore |
Ignores generated mesh .dat files under tests/data/ and verification outputs. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
SPIDER uses short field names (temp_s, S_s, phi_s) not long names (temperature_s, entropy_s, melt_fraction_s). Previous version produced empty plots.
Reject unsupported MESH_SOURCE values (only 0 and 1 are valid) with a clear error instead of falling through to a null-pointer crash. Also check that -mesh_external_filename is actually provided when MESH_SOURCE=1, catching the common typo of forgetting the filename flag. Addresses bot review comment on PR #5.
The external_mesh_roundtrip CI test fails because generate_aw_mesh.py imports scipy.optimize.fsolve, but scipy is not in py/requirements.txt. Replace fsolve with a scalar Newton iteration using only numpy. The xi-to-r mapping is monotonic, so convergence is guaranteed from the evenly-spaced initial guess. Results agree with scipy to ~1e-11 relative error.
|
Does this PROTEUS PR FormingWorlds/PROTEUS#648 depend on the changes in this SPIDER branch? If it does, would it make sense to tackle the changes here first before looking into the PROTEUS PR? |
|
Yes, but it is not yet ready. I will tag you when it is. It is not quite as easy as this unfortunately, because for these coupled PRs the choices in one branch strongly affect the outcome in PROTEUS; they are all very much connected to each other. Zalmoxis T-P profile impacts strongly the SPIDER calculation and the PROTEUS time evolution and vice versa. |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 344919a364
ℹ️ 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".
|
Now it's ready @nichollsh |
After reading an external mesh file, compare the file's surface radius and CMB/surface ratio against the -radius and -coresize options. Emit a warning if they differ by more than 1%. This catches silent misconfiguration when running SPIDER standalone with a mesh file whose geometry doesn't match the command-line parameters (mixing length and other physics use -radius/-coresize directly). In the PROTEUS coupling pathway, _coresize_from_mesh() already ensures consistency, so this warning is a belt-and-suspenders safeguard for standalone use.
Response to automated review (round 2, 2 March 2026)Two new comments from the overnight Codex review:
No changes to existing tests or behaviour. The geometry warning is purely diagnostic. |
There was a problem hiding this comment.
Overview
Thanks, Tim. Will be fantastic to have SPIDER working for higher mass planets.
Creating these Python scripts to emulate the behaviour is a good idea. The code changes all make sense to me, and SPIDER compiles and runs fine on my computer. However, there are a couple of small things below.
1) Radius dimensionalisation
When creating a mesh using the generate_aw_mesh.py script and running SPIDER with the new options, I get this error:
WARNING: external mesh surface radius (6.371000e+06 m) differs from -radius (6.371000e-02 m) by >1%. Ensure -radius matches the file.
I believe that this is because a dimensionalised variable is being compared to a non-dimensionalised one, near line 151 of mesh.c. This might also be silently causing problems with the coresize calculation on the following lines.
2) Version numbering
These changes are the first to spider in a long time. We should iterate the version numbering scheme to mark this. I suggest switching to calver: 26.03.02.
3) Precision warning
I wonder if it would be worthwhile adding a warning (or even hard-exiting) if the float128 precision happens to be enabled? This is checked on line 92 of main.c.
The geometry validation added in 597eea9 compared a dimensionalised surface radius (metres) against the non-dimensional P->radius, so the warning fired every time the external mesh was used. Now both sides of the comparison use non-dimensional values; the warning message still prints metres for readability. Also adds a compile-time error if PETSC_USE_REAL___FLOAT128 is defined, since SetMeshFromExternalFile uses fscanf(%lf) which writes double into PetscScalar — undefined behaviour under quad precision. Reported by @nichollsh in PR review.
|
Thanks for the review @nichollsh — addressed all three points. 1) Radius dimensionalisationGood catch. The validation was comparing 2) Version numberingAgreed. I'll add a calver bump to 3) Float128 precisionAdded a compile-time |
Adopts calver as suggested by @nichollsh, marking the first SPIDER release with external mesh support.
nichollsh
left a comment
nichollsh
left a comment
There was a problem hiding this comment.
Changes look good! The new tools work for me. I also checked that PROTEUS all_options.toml still functions as expected.
Happy for this to be merged.
Are you planning on making a release?
|
Thanks! Yes, right now after merge. |
github-project-automation
bot
moved this from In Progress
to Done
in PROTEUS Development Roadmap
Overview
This PR adds two capabilities to SPIDER that were not previously available:
External mesh input: SPIDER can now read its radial mesh (the grid of radius, pressure, density, and gravity values that defines the planet interior) from an external file, rather than computing it internally from the Adams-Williamson analytical equation of state. This is the foundation for coupling SPIDER with external structure solvers that provide more physically accurate density profiles.
Equation of state safety: When SPIDER encounters pressure or entropy values outside the range of its lookup tables — which happens for massive planets with extreme compositions — it now warns the user and clamps values gracefully, instead of silently producing NaN values that crash the simulation.
How this connects to the other two PRs
This is one of three coordinated PRs that together enable coupled interior structure and thermal evolution for rocky planets:
The merge order is: this PR + Zalmoxis PR first (independent of each other), then the PROTEUS PR (which depends on both).
Why external mesh input matters
Previously, SPIDER computed its own radial mesh using the Adams-Williamson equation of state — an analytical approximation that assumes the interior is a single self-compressing material. This works reasonably well for Earth-like planets, but becomes a poor approximation for planets with very different iron core fractions (e.g. 10% vs Earth's 32.5%) or masses above about 2 Earth masses. By reading the mesh from a file, SPIDER can now use density profiles computed by a full structure solver (Zalmoxis) that accounts for distinct core and mantle layers, temperature-dependent material properties, and phase transitions.
Why equation of state safety matters
For massive planets (3+ Earth masses) with small iron cores (10% of total mass), the pressure at the core-mantle boundary exceeds 400–600 GPa. At these extreme pressures, the equation of state lookup tables for silicate materials are being extrapolated beyond their validated range. Before this PR, such extrapolation could produce a negative thermal expansion coefficient, which caused the mixing-length convection model to take the square root of a negative number, producing NaN values that immediately crashed the ODE time integrator. The new safety clamp prevents this crash, and the warning system helps users understand when their simulation conditions exceed the equation of state data coverage.
Changes
External mesh input (9 files changed, ~1,300 lines added):
parameters.h,parameters.c: Two new command-line options (-MESH_SOURCEto select internal vs external mesh, and-mesh_external_filenameto specify the file path), with input validationmesh.c: New functionSetMeshFromExternalFile()that reads radius, pressure, density, and gravitational acceleration for both cell-boundary and cell-center nodes from a plain text file. The file uses SI units with surface-to-core ordering and negative gravity (SPIDER's sign convention for inward-pointing acceleration)Equation of state safety (3 files changed):
interp.h,interp.c: Per-table warning counters that emit a message on the first out-of-range lookup for each table and direction (e.g. "entropy above table maximum in solid MgSiO3"), then silently clamp subsequent lookupsmatprop.c: Clamp negative thermal expansion coefficient to zero with a one-time warning; guard a division by the thermal expansion coefficient against the zero caseDocumentation and tests (5 files added):
README.md: External mesh file format specification and usage examplesValidation
Test 1: Round-trip consistency — does reading an external mesh reproduce the built-in result?
To verify that the external mesh pathway does not introduce errors, I exported SPIDER's internal Adams-Williamson mesh to a file and then ran the same simulation reading that file back. The two runs should produce identical results.
Result: Temperature, entropy, and melt fraction profiles match within 0.15% relative error at all radial positions and time steps. The small residual comes from the difference between SPIDER's analytical mass integral and the discrete midpoint-rule integral used to generate the file.
Test 2: Field-by-field validation of external mesh
To confirm that SPIDER correctly reads and uses every field from the external mesh file (not just radius), I compared all output fields from the external-mesh run against the internal-mesh reference.
Why this test: A mesh file contains four quantities (radius, pressure, density, gravity) and SPIDER computes three derived quantities from the thermal evolution (temperature, entropy, melt fraction). All seven must match between the internal and external mesh pathways — if any field is silently dropped or misread, the simulation would evolve on an incorrect state.
Result: All seven physical quantities pass validation:
The ~0.1% errors are the expected numerical floor from the discretisation difference between SPIDER's analytical mass integral and the midpoint-rule integral used to generate the external mesh file. The derived quantities (temperature, entropy, melt fraction) have smaller errors than the input quantities (radius, pressure, density, gravity) because they are determined primarily by the equation of state, which is evaluated identically in both runs.
Test 3: Non-analytical mesh — does SPIDER accept arbitrary density profiles?
The round-trip test uses the same density profile that SPIDER would have computed internally. To verify that SPIDER works with genuinely different density profiles, I generated a mesh with a linear density gradient (not physical, but a useful stress test) and ran SPIDER with it.
Why this test: The external mesh feature is only useful if SPIDER can handle density profiles that differ substantially from its built-in Adams-Williamson model. A linear gradient is maximally different from the self-compressing analytical profile.
Result: SPIDER accepts and evolves on the non-analytical mesh without errors. The surface temperature cools monotonically as expected for a blackbody cooling problem.
Test 4: Thermal expansion safety clamp — does SPIDER survive extreme pressures?
At pressures above ~400 GPa (reached at the core-mantle boundary for 3+ Earth mass planets with small iron cores), the equation of state lookup tables for solid-phase silicate are being extrapolated beyond their validated range. Before this PR, such extrapolation could produce a negative thermal expansion coefficient α, which entered the mixing-length convection model through the expression √(α · g · ΔT). Taking the square root of a negative number produced NaN values that immediately crashed the ODE time integrator at the first step.
Why this test: Without this fix, SPIDER cannot simulate any planet where the core-mantle boundary pressure exceeds the equation of state table range (~430 GPa for the WolfBower2018 solid-phase MgSiO₃ tables). This rules out all planets above ~2.5 Earth masses with iron core fractions below ~20%.
The fix: When the thermal expansion coefficient is negative (indicating extrapolation beyond the table), it is clamped to zero with a one-time warning. This effectively disables convective mixing at that node — a physically reasonable fallback, because the equation of state is unreliable there anyway. The guard also prevents a division-by-zero in the critical Rayleigh number computation (which divides by 4α).
Result: The most extreme test case (3 Earth masses, 10% iron core, core-mantle boundary pressure ~670 GPa) now completes 136 time steps and evolves to 1 million years of simulated time. Before this fix, it crashed at the first step. The clamp fires for the deepest ~5 nodes (those closest to the core-mantle boundary where pressure exceeds the table range) and does not affect the well-constrained upper mantle.
Test 5: Equation of state out-of-range warning system
When SPIDER requests a pressure or entropy value outside the range covered by the equation of state lookup tables, the code now emits a diagnostic warning on the first occurrence for each table and direction (e.g. "entropy above table maximum in solid MgSiO₃ at S = 2395.7 J/(kg·K), table max = 2395.3 J/(kg·K)"). Subsequent out-of-range lookups for the same table and direction are silently clamped to the nearest table edge, avoiding log spam while still alerting the user to the extrapolation.
Why this test: Before this PR, out-of-range lookups were silently clamped with no indication to the user. A simulation could complete but produce physically incorrect results because the equation of state was being extrapolated without warning. The per-table, per-direction tracking (implemented as counters in the
Interp2dstruct) ensures that users see exactly which tables are being exceeded and in which direction, making it straightforward to diagnose whether higher-coverage equation of state tables (such as RTPress100TPa) are needed for a given planet configuration.Result: Warnings correctly fire for the 3 Earth mass, 10% iron core case (where solid-phase entropy exceeds the WolfBower2018 table maximum) and do not fire for well-constrained cases (1 Earth mass, 32.5% iron core).
Challenges during development
What was not tackled / future work
Dependencies
Companion PRs: