Repository navigation
Export the melt drainage quantities at the basic nodes - #66
Merged
Merged
Conversation
SolverOutput and the NetCDF snapshot now carry porosity_b, rho_solid_b, rho_melt_b and g_b, for consumers that compute melt drainage. They are read from the basic-node phase evaluator after get_state refreshes the state at the final entropy, so they describe the returned state. The three-regime permeability over porosity moves to a module-level mobility_function, and the density-derived porosity to EntropyPhaseEvaluator.porosity, so a consumer evaluates the same expressions the separation velocity uses. relative_velocity calls both and its result is unchanged. Under const_properties the porosity is 1 and both densities equal const_rho.
Below the solidus the soft-clipped porosity lies in [-2.5e-7, 5e-4], reaching 5e-4 at the solidus; the docstring and netcdf.md now say so, and relative_velocity's docstring says that where the melt is denser the porosity goes to the lower or the upper clip limit depending on the node density. mobility_function is exported from aragog.eos, and usage.md lists the four new SolverOutput fields. The solve-based export tests carry the smoke marker. New tests pin both blend centres of mobility_function, the porosity bound at the solidus, and g_b against an external gravity profile.
The exported porosity is (rho_s - rho)/(rho_s - rho_l), soft-clipped, not the melt fraction. A solid node whose density is below rho_s, which happens in tables where the solid density rises with entropy towards the solidus, gives values up to 1, and where the melt is denser the floored denominator gives a band about 1 kg/m^3 wide in which any value occurs. The labels, the porosity() and relative_velocity docstrings and netcdf.md now say this and point consumers to phi_basic for the presence of melt. Tests pin the clip value at the solidus, the porosity formula across the dense-melt band, and g_b against a non-linear external gravity profile. The export tests keep the unit marker, like the other short EOS-backed solves that PR CI runs.
The porosity denominator is floored at about 1 kg/m^3 wherever rho_s - rho_l is below that, a lighter melt included. There the porosity is (rho_s - rho)/(1 kg/m^3) clipped, so the intermediate band lies within 1 kg/m^3 below rho_s and a node at or above rho_s is about 0. The relative_velocity and mobility_function docstrings, netcdf.md and two_phase_flow.md now say this; two_phase_flow.md also notes that j_grav carries phi(1 - phi) of the melt fraction. The porosity formula test covers both sides of the band, the lighter-melt floor and the lower clip limit at 1e-12. Solve-based tests are smoke and formula tests unit, as build_tests.md defines the tiers.
relative_velocity is compared bit for bit with values computed before the extraction, for a normal and a dense-melt contrast in both separation modes, using a stub EOS so the unit tier runs it without tables. Unit tests also cover the order of phase_boundary_densities, the porosity ratio, the floor band for a melt lighter or denser by 0.5 kg/m^3, mobility_function at the negative clip values, and the const_properties copies; the NetCDF round trip checks the units of the four exports. Under const_properties phase_boundary_densities returns copies of the node density. The porosity() docstring and netcdf.md state that porosity_b is about 0 in a mush where rho_s <= rho_l, so drainage must test the sign of the contrast first, and that Phi_global_vol uses the unfloored ratio. The mobility_function docstring says the grain size enters as the radius of Stokes' law and notes the 4.4 percent overshoot of the upper blend; netcdf.md describes g_b.
There was a problem hiding this comment.
Copilot review overview
🟢 Approval recommended
The implementation preserves existing calculations while consistently exposing and testing the new diagnostics.
Review effort: Balanced
Findings: None
What changed in this PR
Exposes melt-drainage inputs at basic nodes without changing solver results.
Changes:
- Adds public porosity, phase-density, and mobility APIs.
- Exports drainage arrays through
SolverOutputand NetCDF. - Adds documentation and comprehensive regression coverage.
| File | Description |
|---|---|
src/aragog/eos/entropy_phase.py |
Extracts drainage calculations into public APIs. |
src/aragog/eos/__init__.py |
Exports mobility_function. |
src/aragog/solver/entropy_solver.py |
Adds final-state drainage arrays and serialization. |
tests/test_drainage_exports.py |
Tests formulas, exports, and solver modes. |
tests/test_solver_output_netcdf.py |
Tests NetCDF values, dimensions, and units. |
docs/Reference/api/aragog.eos.md |
Documents the mobility API. |
docs/How-to/usage.md |
Lists new output fields. |
docs/How-to/netcdf.md |
Explains drainage output semantics. |
docs/Explanations/two_phase_flow.md |
Clarifies density-crossing behavior. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
7 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
Exposes the quantities that the melt-solid separation velocity is built from, for consumers that compute melt drainage outside the solver.
mobility_function(porosity, grain_size)inaragog.eosis the three-regime permeability factor F = k(phi)/phi (Blake-Kozeny-Carman, Rumpf-Gupte, Stokes; Abe 1995, Bower et al. 2018).EntropyPhaseEvaluator.porosity()andphase_boundary_densities()return the density-derived porosity and the solid and melt densities at the phase boundary.relative_velocitynow calls these three; the expressions are unchanged.SolverOutputand the NetCDF file carry four new basic-node arrays:porosity_b,rho_solid_b,rho_melt_bandg_b. They are computed once from the final state, after the integration, on every backend.porosity_bis and is not: it is the porosity the separation velocity uses, not the melt fraction. Where the melt is denser than the solid it is about 0 although the node holds melt, so a drainage calculation tests the sign ofrho_solid_b - rho_melt_bfirst and usesphi_basicfor the presence of melt. The grain size enters as a radius (Stokes prefactor 2/9).No solver result changes.
Validation of changes
Tested on macOS 26 with Python 3.12.
Checklist