diff --git a/.agents/memory.instruction.md b/.agents/memory.instruction.md deleted file mode 100644 index 83fdd46..0000000 --- a/.agents/memory.instruction.md +++ /dev/null @@ -1,31 +0,0 @@ ---- -applyTo: '**' ---- - -# pyMAGEMin Memory (repo-local) - -## Coding Preferences - -- Use NumPy for arrays and `scipy.optimize.root_scalar` for scalar minimization/root finding. -- Follow existing patterns in `MAGEMinGarnetCalculator` and `PhaseFunctions`. -- Julia conversion: use `jlconvert(jl.Vector[...], data)` for lists/arrays before MAGEMin_C calls. -- Reuse `data = MAGEMin_C.Initialize_MAGEMin("ig", verbose=False)` across calls for speed. -- Keep edits minimal and localized in `src/pyMAGEMin/functions/`. - -## Project Architecture - -- **Core:** pyMAGEMin wraps Julia's MAGEMin_C thermodynamic minimizer. -- **Bootstrap:** `MAGEMin_C` is initialized in `src/pyMAGEMin/__init__.py`. -- **Key functions:** - - `MAGEMinGarnetCalculator.gt_along_path(...)` for garnet chemistry/fractionation along a P-T path - - `MAGEMinGarnetCalculator.gt_single_point_calc_elements(...)` for single P-T points - - `MAGEMinGarnetCalculator.generate_2D_grid_gt_elements(...)` for grid workflows -- **Element fractions:** Mg, Fe, Ca, Mn (returned in mol or wt basis depending on `sys_in`). -- **Bulk rock:** wt%, mol%, or vol% depending on `sys_in` (`'wt'`, `'mol'`, `'vol'`). - -## Solutions Repository - -### MAGEMin Workflow Notes -- **Path calculations:** `gt_along_path(...)` returns `X_along_path`; preserve it when `fractionate=True`. -- **Single-point checks:** use `gt_single_point_calc_elements(...)` for Fe-Mg-Mn-Ca outputs. -- **Grid checks:** use `generate_2D_grid_gt_endmembers(...)` or `generate_2D_grid_gt_elements(...)`. diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md deleted file mode 100644 index 0911dbc..0000000 --- a/.github/copilot-instructions.md +++ /dev/null @@ -1,67 +0,0 @@ -# pyMAGEMin — AI Coding Assistant Instructions - -This repo provides Python wrappers around the Julia MAGEMin_C thermodynamic minimizer to compute phase equilibria and garnet growth paths. These notes capture the project-specific patterns, workflows, and integration points to help you be productive fast. - -## Architecture and Key Modules -- Core package: pyMAGEMin (Python, src layout). See [src/pyMAGEMin/__init__.py](src/pyMAGEMin/__init__.py) for `juliacall` initialization of the Julia module `MAGEMin_C` made globally available as `pyMAGEMin.MAGEMin_C`. -- MAGEMin wrappers and helpers: - - [src/pyMAGEMin/functions/MAGEMin_functions.py](src/pyMAGEMin/functions/MAGEMin_functions.py): `MAGEMinGarnetCalculator` (grid, single-point, and path calculations, optional fractionation) and `PhaseFunctions` (solidus/liquidus search via `scipy.optimize.root_scalar`, batch fractionation). - - [src/pyMAGEMin/functions/bulk_rock_functions.py](src/pyMAGEMin/functions/bulk_rock_functions.py): composition utilities (mol↔wt conversions and related helpers). - - [src/pyMAGEMin/functions/garnet_growth.py](src/pyMAGEMin/functions/garnet_growth.py): builds size distributions and radial profiles from P–T–t paths using results from `MAGEMinGarnetCalculator`. - - [src/pyMAGEMin/functions/utils.py](src/pyMAGEMin/functions/utils.py): `create_PTt_path()` P–T–t interpolation with embedded original points. - -## Julia Integration (critical) -- Requires a local Julia installation and the Julia package `MAGEMin_C`. The installer checks/installs both (see [setup.py](setup.py)). -- Cross-language data: use `from juliacall import Main as jl, convert as jlconvert` then convert to Julia types, e.g. `jlconvert(jl.Vector[jl.Float64], np_array)` and `jlconvert(jl.Vector[jl.String], list_of_names)`. -- Initialize MAGEMin once per session: `data = MAGEMin_C.Initialize_MAGEMin("ig", verbose=False)`; reuse `data` across calls for performance. - -## Typical Call Flow -- Single or multi-point minimization: `MAGEMin_C.single_point_minimization(P, T, data, X=..., Xoxides=..., sys_in=...)` or `multi_point_minimization(P_vec, T_vec, data, ...)`. -- Higher-level APIs in `MAGEMinGarnetCalculator` wrap the above and compute endmember and element fractions for phase `"g"` (garnet). Helpers `phase_frac()` and `extract_end_member()` encapsulate MAGEMin_C output access patterns. -- Units and conventions: - - `sys_in`: `'mol'`, `'wt'`, or `'vol'` where applicable; keep consistent across pipelines. - - Garnet endmembers: `py`, `alm`, `spss`, `gr`, `kho`; element mapping is handled internally via `_extract_garnet_elements_from_oxides` and related helpers. - - Phase keys: garnet is `'g'`; many helpers assume this. - -## Minimal Usage Example -```python -import numpy as np -from pyMAGEMin import MAGEMin_C -from juliacall import Main as jl, convert as jlconvert - -data = MAGEMin_C.Initialize_MAGEMin("ig", verbose=False) -P = jlconvert(jl.Vector[jl.Float64], np.linspace(1, 10, 5)) -T = jlconvert(jl.Vector[jl.Float64], np.linspace(700, 900, 5)) -Xox = jlconvert(jl.Vector[jl.String], ["SiO2","Al2O3","CaO","MgO","FeO","Fe2O3","K2O","Na2O","TiO2","Cr2O3","H2O"]) -X = jlconvert(jl.Vector[jl.Float64], [48.43,15.19,11.57,10.13,6.65,1.64,0.59,1.87,0.68,0.0,3.0]) -out = MAGEMin_C.multi_point_minimization(P, T, data, X=X, Xoxides=Xox, sys_in="wt") -``` - -## Developer Workflows -- Install (ensures Julia and MAGEMin_C): -```bash -python -m pip install -e . -``` -- Quick check Julia availability if install fails: -```bash -julia --version -julia -e 'using Pkg; Pkg.status("MAGEMin_C")' -``` -- Parallel example (MPI): run your own driver script under MPI, e.g.: -```bash -mpiexec -n 4 python your_script.py -``` -- Tutorials and examples: notebooks under [Tutorials/](Tutorials) demonstrate workflows end-to-end. - -## Patterns and Gotchas -- Always convert Python arrays/lists to Julia vectors with `jlconvert(...)` before calling MAGEMin_C; prefer contiguous NumPy arrays for speed. -- Reuse `data` returned by `Initialize_MAGEMin(...)` across calls to avoid re-initialization overhead. -- `PhaseFunctions.find_phase_in/saturation` use bisection with brackets; provide realistic `(T_low, T_high)` for robust convergence. -- Fractionation: `PhaseFunctions.fractionate_phase()` updates bulk composition between steps; honor `sys_in` and return value when chaining along paths. -- For garnet element outputs, prefer the `gt_*_elements(...)` methods in [src/pyMAGEMin/functions/MAGEMin_functions.py](src/pyMAGEMin/functions/MAGEMin_functions.py). - -## Extending the Package -- New wrappers should follow the `MAGEMin_functions.py` pattern: convert inputs via `jlconvert`, call MAGEMin_C, then post-process into NumPy arrays and Python dicts. -- Keep phase keys and `sys_in` handling consistent; update helpers if introducing new phases or endmembers. -- Place general utilities in `functions/` and import in [src/pyMAGEMin/__init__.py](src/pyMAGEMin/__init__.py) if they should be top-level accessible. - diff --git a/.gitignore b/.gitignore index fc207c6..4ba4f88 100644 --- a/.gitignore +++ b/.gitignore @@ -154,3 +154,4 @@ cython_debug/ docs /examples +/development diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 0000000..566560d --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,288 @@ +# AGENTS.md + +## Quick start + +```bash +python -m pip install -e . +phasetools-julia-setup --check # verify Julia + MAGEMin_C +phasetools-julia-setup --install # install MAGEMin_C in Julia +python3 -m unittest discover tests # run all tests (mock-based, no Julia needed) +``` + +## Package structure + +``` +src/phasetools/ + __init__.py — initialises MAGEMin_C via juliacall; exports public API + julia_setup.py — CLI for phasetools-julia-setup (check/install Julia + MAGEMin_C) + core/ + base.py — MAGEMinBase (parent class for all calculators/models) + engine.py — single_point_minimization_with_conversion (Python→Julia bridge) + phase_properties.py — phase_frac, extract_end_member, get_oxide_apfu, + get_phase_chemistry, get_phase_mg_number, + get_phase_fe_split, calculate_kd_fe_mg + calculators/ + pt_grid.py — MAGEMinPTGridCalculator (multi-point P-T grids) + garnet.py — MAGEMinGarnetCalculator (garnet-focused wrappers) + phase_search.py — PhaseFunctions (solidus/liquidus root-finding) + pt_estimation.py — PhasePTEstimator (geothermobarometry) + assemblage.py — MAGEMinAssemblageCalculator (stability masks) + models/ + garnet_growth.py — GarnetGenerator, generate_distribution + magma_ocean.py — MagmaOcean + utils/ + bulk_rock.py — mol↔wt conversions, molar mass dicts, atomic_frac_to_wt_frac, + non-normalising mol↔wt fraction converters (single-component + or subset conversions, e.g. FeO/Fe2O3/FeOt), + split_feot_to_feo_o / express_bulk_in_feo_o_basis + (FeOt→FeO+O split at a target Fe³⁺/FeOt, conserving FeOt) + general.py — PCHIP interpolation helpers +``` + +## Conventions (must follow) + +- **British English** throughout: `colour`, `standardise`, `crystallisation`, `modelling`, `organisation`, `recognise`. +- **NumPy docstrings** and **type hints** on all public methods. +- **P in kbar, T in °C** for all user-facing inputs. +- Garnet phase key is always `'g'`. +- Every calculator or model class **must inherit from `MAGEMinBase`**. +- High-level classes accessible via `from phasetools import ...`. + +## Julia bridge (critical) + +```python +from juliacall import Main as jl, convert as jlconvert +``` + +- Convert arrays/lists before MAGEMin_C calls: + `jlconvert(jl.Vector[jl.Float64], np_array)` or `jlconvert(jl.Vector[jl.String], list_of_names)` +- **Avoid double-converting** objects already in Julia format — verify the `juliacall` wrapper. +- Initialise once and **reuse the `data` handle**: + `data = MAGEMin_C.Initialize_MAGEMin("ig", verbose=False)` +- Share state across instances: `self._copy_state_to(other)` transfers standardised Julia objects. + +## Databases + +### Standard databases + +| Label | Chemical system | +|-------|----------------| +| `ig` (igneous) | `K2O-Na2O-CaO-FeO-MgO-Al2O3-SiO2-H2O-TiO2-O-Cr2O3` | +| `igad` (igneous alkaline) | `K2O-Na2O-CaO-FeO-MgO-Al2O3-SiO2-TiO2-O-Cr2O3` | +| `mp` (metapelite) | `K2O-Na2O-CaO-FeO-MgO-Al2O3-SiO2-H2O-TiO2-O-MnO` | +| `mb` (metabasite) | `K2O-Na2O-CaO-FeO-MgO-Al2O3-SiO2-H2O-TiO2-O` | +| `um` (ultramafic) | `SiO2-Al2O3-MgO-FeO-O-H2O-S` | +| `mtl` / `sb11` / `sb21` (mantle) | `Na2O-CaO-FeO-MgO-Al2O3-SiO2` | +| `sb24` (mantle) | `Na2O-CaO-FeO-MgO-Al2O3-SiO2-O-Cr2O3` (includes metal Fe phases like `Fe(a)`) | + +### Extended databases + +- **`mpe`**: mp + mb (cross-lithology metapelite) +- **`mbe`**: mb + high-pressure models +- **`ume`**: um + mb (mantle-crust interactions) + +### Applicability warning + +All datasets calibrated for limited P-T-X ranges — exceeding limits causes non-convergence or inaccuracies. For databases without `O` (e.g., `mtl`), `FeO` is treated as total iron. + +## Bulk composition + +Use `setup_bulk_composition(Xoxides, X, sys_in, rm_list=None)` (from `MAGEMinBase`). It calls `MAGEMin_C.convertBulk4MAGEMin` to standardise, builds a stoichiometry map, and optionally removes phases. + +- If `sys_in='wt'`: standardised molar output is converted back to wt%. +- Otherwise stays in molar fractions, normalised to sum 100.0. +- Core oxides (SiO₂, Al₂O₃, MgO) clamped to min **1e-4** molar fraction. +- Optional oxides (TiO₂, MnO, H₂O etc.) set to 0.0 if < **2e-5**. +- `phase_frac()` handles solvi — automatically sums multiple instances of the same phase. + +### Unit systems + +When comparing MAGEMin model output to measured data (e.g., EPMA), **convert everything to a consistent unit system** using `MAGEMin_C.convertBulk4MAGEMin`: + +```python +x_jl = jlconvert(jl.Vector[jl.Float64], wt_list) +ox_jl = jlconvert(jl.Vector[jl.String], ox_names) +x_mol_jl, xox_mol_jl = MAGEMin_C.convertBulk4MAGEMin(x_jl, ox_jl, "wt", db) +``` + +- Pad missing oxides with zeros when targets have fewer oxides than the bulk rock. +- Common approach: set `sys_in='mol'` on the calculator and pre-convert the bulk rock from wt% to mol% before calling `setup_bulk_composition`. +- X-site (cation) fractions are **always atomic/molar** — the `cat_*` fields from `extract_from_grid` return atomic fractions regardless of `sys_in`. Use `atomic_frac_to_wt_frac()` (in `bulk_rock.py`) to convert if weight-based fractions are needed. + +## Redox model + +- MAGEMin uses **FeO + O (excess oxygen)** basis for most databases (`ig`, `mp`, etc.). +- `1 mol Fe2O3 → 2 mol FeO + 1 mol O` +- `get_phase_fe_split(out, phase)` for heuristic Fe²⁺/Fe³⁺ partitioning. +- `get_phase_mg_number` atomic across Fe, FeO, Fe2O3 components. +- `calculate_kd_fe_mg(out, phase1, phase2)` for distribution coefficients. +- **Garnet X-site fractions default to FeOt** (all Fe as Fe²⁺): `MAGEMinGarnetCalculator` / `GarnetGenerator` take a `fe_basis` argument (`'FeOt'` default, `'Fe2+'` for the stoichiometric ferrous-only split). FeOt is the community convention for garnet end-members (e.g. Williams & Grambling 1990; Krogh Ravna 2000); use `'Fe2+'` only when garnet Fe³⁺ is known to be significant. +- sb24 uses pure stoichiometric iron (`Fe(a)`) with excess O. + +## Endmember formulae + +> **Source:** Extracted from MAGEMin source code at [`ComputationalThermodynamics/MAGEMin`](https://github.com/ComputationalThermodynamics/MAGEMin): +> - TC database endmember entries: [`src/TC_database/TC_endmembers.c`](https://github.com/ComputationalThermodynamics/MAGEMin/blob/main/src/TC_database/TC_endmembers.c) +> - SB database endmember entries: [`src/SB_database/SB_endmembers.c`](https://github.com/ComputationalThermodynamics/MAGEMin/blob/main/src/SB_database/SB_endmembers.c) +> - Endmember listing from Python: `calc.get_phase_endmembers(phase, grid_out)` + +### Oxide index ordering + +The `Comp[]` array in endmember database entries uses a **database-specific** oxide index ordering: + +| Database family | Oxide ordering (index: oxide) | Used by | +|----------------|-------------------------------|---------| +| **TC** | `0:SiO₂, 1:Al₂O₃, 2:CaO, 3:MgO, 4:FeO, 5:K₂O, 6:Na₂O, 7:TiO₂, 8:O, 9:MnO, 10:Cr₂O₃, 11:H₂O, 12:CO₂` | `ig`, `igad`, `mp`, `mb`, `mpe`, `mbe`, `um`, `ume` | +| **SB** sb11/sb21 | `0:SiO₂, 1:CaO, 2:Al₂O₃, 3:FeO, 4:MgO, 5:Na₂O` | `mtl`, `sb11`, `sb21` | +| **SB** sb24 | `0:SiO₂, 1:CaO, 2:Al₂O₃, 3:MgO, 4:Na₂O, 5:O, 6:Cr₂O₃, 7:Fe` | `sb24` | + +The last element of `Comp[]` (index 15 for TC, index 17 for SB) is the **total number of atoms per formula unit**. + +### Key endmember formulae + +``` +Garnet (g): py = Mg3Al2Si3O12 alm = Fe3Al2Si3O12 + gr = Ca3Al2Si3O12 spss = Mn3Al2Si3O12 + +Cpx (dio): jd = NaAlSi2O6 di = CaMgSi2O6 + hed = CaFeSi2O6 acm = NaFe³⁺Si2O6 (with O) + +Opx (opx): en = Mg2Si2O6 fs = Fe2Si2O6 + mgts = MgAlAlSiO6 + +Feldspar (fsp): ab = NaAlSi3O8 an = CaAl2Si2O8 + san = KAlSi3O8 + +Amphibole (amp): tr = Ca2Mg5Si8O22(OH)2 + ts = Ca2Mg3Al2Si6O22(OH)2 + gl = Na2Mg3Al2Si8O22(OH)2 + parg = NaCa2Mg4Al3Si6O22(OH)2 + +Mica (bi): phl = KMg3AlSi3O10(OH)2 + ann = KFe3AlSi3O10(OH)2 + east = KMg2Al3Si2O10(OH)2 + +Mica (mu): mu = KAl3Si3O10(OH)2 + cel = KMgAlSi4O10(OH)2 + fcel = KFeAlSi4O10(OH)2 +``` + +### Critical pitfall: solution-model pseudo-endmembers + +Some endmembers returned by `get_phase_endmembers(phase, out)` are **not in `TC_endmembers.c`**. They are **pseudo-endmembers** (ordering parameters / intermediate compositions) defined programmatically in the solution model C code (`src/TC_database/SS_xeos_PC_*.c`). + +**Known pseudo-endmembers:** +| Phase | Pseudo-endmember | Notes | +|-------|-----------------|-------| +| cpx (dio) | `om` (omphacite) | ~Na₀.₅Ca₀.₅ intermediate — carries BOTH Na and Ca | +| cpx (dio) | `acmm` | Variant of acmite (NaFe³⁺Si₂O₆), model-specific | +| cpx (dio) | `cfm` | Ca-Fe-Mg ordering component | +| cpx (dio) | `jac` | Jadeite-acmite intermediate | + +**APFU vs endmembers: when to use which:** +| Need | Use | Why | +|------|-----|-----| +| Element ratios (XJd, XMg, XAn) | **APFU** via `oxides=[...]` or `cations=[...]` | Reads structural formula directly; avoids pseudo-endmember misclassification | +| Activity/composition modelling | **Endmembers** via `end_members=[...]` | Endmembers are the thermodynamic mixing components | +| Verifying endmember assignments | **Both + least squares** | `APFU = Σ em_frac × Comp` — if the fit fails, your endmember set or formula assignment is wrong | + +Since endmembers are just oxide arrays in the database, APFU already encodes the same information without the risk of summing the wrong subset. For composition comparison with measured data, **APFU is always simpler and more reliable.** + +## Calculator API patterns + +- **PTGrid**: `calculate_grid(P, T)` → `extract_from_grid(phase, end_members='auto')` +- **Garnet**: `gt_along_path(P, T, fractionate=True)` for evolution with zoning; X-site normalisation isolates Fe²⁺. +- **PTEstimator**: wraps `scipy.optimize` (global: `differential_evolution`, `dual_annealing`; local: `shgo`, `minimize`). +- **Assemblage**: returns boolean stability masks for requested phase coexistence. +- **PhaseSearch**: uses `scipy.optimize.root_scalar` with P,T brackets for solidus/liquidus. + +## Comparing EPMA data with MAGEMin outputs + +### Extraction methods from `extract_from_grid` + +| Key pattern | Source | Units | What it returns | +|-------------|--------|-------|----------------| +| `ox_apfu_{oxide}` | `oxides=[...]` → `get_oxide_apfu()` | Atoms per formula unit | Structural formula on a per-formula-unit basis (normalised to the phase's oxygen count, e.g. 6 O for cpx) | +| `chem_{oxide}` | `chemistry=[...]` → `get_phase_chemistry()` | mol% (if `sys_in='mol'`) or wt% (if `sys_in='wt'`) | Oxide concentration in the phase, on the same basis as `sys_in` | +| `cat_{cation}` | `cations=[...]` → `_extract_cations_from_apfu()` | Atomic fraction (0–1) | Cation site fractions, always atomic regardless of `sys_in` | + +**For comparing with EPMA, use `ox_apfu_*` (APFu)** — this is the standard mineralogical normalisation (atoms per formula unit, typically 6 O for pyroxene, 12 O for garnet, 22 O for amphibole) and matches how EPMA data is conventionally reported. + +### Scale-invariant ratios (always prefer these) + +Ratios like XJd = Na/(Na+Ca) and Mg# = Mg/(Mg+Fe) are **independent of the normalisation basis** — the denominator cancels whether you use APFu, mol%, or wt%. They are the safest quantities to compare between model and measurement: + +```python +# From APFu (model side) +xjd = ox_apfu_Na2O / (ox_apfu_Na2O + ox_apfu_CaO) + +# From chemistry (model side) — same result +xjd = (2 * chem_Na2O) / (2 * chem_Na2O + chem_CaO) # ×2 for Na atoms + +# From mol% columns (EPMA side) — same result if done correctly +xjd = (2 * Na2O_mol%) / (2 * Na2O_mol% + CaO_mol%) +``` + +All three give identical values (verified to machine precision). The `×2` on Na₂O converts from oxide molecules to Na atoms (2 Na per Na₂O vs 1 Ca per CaO). + +### Critical pitfall: oxide-molecule vs atom ratios + +**Do not compute element ratios from oxide molecule fractions without accounting for stoichiometry.** Na₂O has 2 Na atoms per molecule; CaO has 1 Ca atom. The oxide-molecule ratio Na₂O/(Na₂O+CaO) is systematically ~half the correct cation ratio Na/(Na+Ca). + +| Formula | Convention | XJd for typical omphacite | +|---------|-----------|---------------------------| +| `Na2O_mol% / (Na2O_mol% + CaO_mol%)` | Oxide molecules | **0.14–0.18** (WRONG) | +| `(2 × Na2O_mol%) / (2 × Na2O_mol% + CaO_mol%)` | Atoms (correct) | **0.25–0.31** | +| `ox_apfu_Na2O / (ox_apfu_Na2O + ox_apfu_CaO)` | APFu atoms (correct) | **0.25–0.31** | + +Mg# = MgO/(MgO+FeO) is **not affected** because both oxides carry 1 cation per molecule. + +### Critical pitfall: wt_to_apfu oxygen counting + +When converting EPMA wt% to APFu, count **oxygen atoms per oxide molecule**, not the number of capital-O letters: + +| Oxide | `ox.count("O")` (WRONG) | Correct O count | +|-------|--------------------------|-----------------| +| SiO₂ | 1 | **2** | +| TiO₂ | 1 | **2** | +| Al₂O₃ | 1 | **3** | +| MgO, CaO, FeO, Na₂O, K₂O, MnO | 1 | 1 | + +The `ox.count("O")` bug undercounts oxygen, inflating the6-O scale factor by ~60% and corrupting all APFu values. + +### Practical checklist for EPMA ↔ MAGEMin comparison + +1. **Load EPMA data** and compute the measured ratios from the correct columns: + - XJd = `2 * Na2O_mol% / (2 * Na2O_mol% + CaO_mol%)` (not `Na2O/(Na2O+CaO)`) + - Mg# = `MgO_mol% / (MgO_mol% + FeO_mol%)` + - Or convert wt% → APFu with correct oxygen counts (2, 2, 3, 1, 1, 1, 2, 2, 1 for SiO₂…FeO) + +2. **Run the model** grid over the same P–T box. Extract with `oxides=["Na2O","CaO","MgO","FeO"]` and `mg_number=True`. + +3. **Compute model ratios** from APFu: + - `XJd_model = ox_apfu_Na2O / (ox_apfu_Na2O + ox_apfu_CaO)` + - `Mg#_model = get_phase_mg_number(out, phase)` or `ox_apfu_MgO / (ox_apfu_MgO + ox_apfu_FeO)` + +4. **Compare ratios**, not absolute APFu. Absolute APFu depends on the normalisation basis and will differ between a6-O structural formula and a mol% concentration — but the ratios are identical. + +5. **Verify agreement**: `np.allclose(xjd_model, xjd_chem, equal_nan=True)` should return `True` (use `equal_nan=True` because missing-phase points are NaN). + +## Models + +- **GarnetGenerator**: radial zoning across shells, cohort size distributions, Rayleigh-style fractionation at each step. +- **MagmaOcean**: equilibrium (stage 0) → fractional crystallisation (stages 1–N); pressure-to-depth/radius conversions for rocky bodies. + +## Tests + +- **unittest**, mock-based (`unittest.mock.patch`), **no live Julia runtime** needed. +- Run: `python3 -m unittest discover tests` +- Files: `test_redox_logic.py`, `test_site_occupancy.py`, `test_lmo_fix.py`. +- New public functions must be documented in the **directory's README.md** (e.g., `calculators/README.md`, `core/README.md`) **and** the package structure above must be updated. + +## Requirements (from `setup.py`) + +`pandas`, `numpy`, `matplotlib`, `scipy`, `molmass`, `juliacall`. Python ≥ 3.10. + +## Tutorials + +Jupyter notebooks in `Tutorials/` (garnet growth, magma ocean, general workflows). diff --git a/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb b/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb new file mode 100644 index 0000000..7b65997 --- /dev/null +++ b/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb @@ -0,0 +1,596 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Garnet–Clinopyroxene Fe²⁺–Mg Exchange Thermometry (S10 Eclogite)\n", + "\n", + "Garnet and clinopyroxene coexist over a wide range of high-pressure mafic rocks, from\n", + "granulites to eclogites, and the Fe²⁺–Mg distribution between them is a long-established\n", + "thermometer. The exchange reaction\n", + "\n", + "$$\n", + "\\tfrac{1}{3}\\mathrm{Mg_3Al_2Si_3O_{12}}\\ (\\text{pyrope}) + \\mathrm{CaFeSi_2O_6}\\ (\\text{hedenbergite})\n", + "\\;=\\;\n", + "\\tfrac{1}{3}\\mathrm{Fe_3Al_2Si_3O_{12}}\\ (\\text{almandine}) + \\mathrm{CaMgSi_2O_6}\\ (\\text{diopside})\n", + "$$\n", + "\n", + "has an equilibrium distribution coefficient\n", + "\n", + "$$\n", + "K_D = \\frac{(\\mathrm{Fe^{2+}/Mg})^{\\mathrm{garnet}}}{(\\mathrm{Fe^{2+}/Mg})^{\\mathrm{clinopyroxene}}}\n", + "$$\n", + "\n", + "whose temperature dependence is the basis of the thermometer: $K_D$ falls as temperature rises.\n", + "\n", + "This notebook:\n", + "\n", + "1. sets up the S10 eclogite bulk composition in the `mpe` database;\n", + "2. runs a P–T grid and extracts the garnet and clinopyroxene compositions;\n", + "3. maps $K_D$ across the grid using the traditional mixed Fe convention: garnet Fe is taken from its FeO-basis APFU (all Fe treated as Fe²⁺), while cpx Fe²⁺ comes from the excess-O heuristic split;\n", + "4. converts $K_D$ to temperature with Mysen & Heier (1972), Ganguly (1979), Krogh Ravna (2000), Ellis & Green (1979), and Räheim & Green (1974), and compares them against the true grid temperature.\n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## 1. Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from phasetools import MAGEMinPTGridCalculator" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## 2. Bulk composition and database\n", + "\n", + "The S10 eclogite bulk composition (in mol%) is taken from the `phasetools` test suite. The\n", + "`mpe` database combines the metapelite (`mp`) and metabasite (`mb`) solution models and is\n", + "the appropriate choice for a mafic bulk at eclogite-facies conditions. Pressures are in kbar\n", + "and temperatures in °C throughout." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "# Eclogite S10 bulk composition in mol% (from the phasetools test suite).\n", + "Xoxides = [\"H2O\", \"SiO2\", \"Al2O3\", \"CaO\", \"MgO\", \"FeO\", \"K2O\", \"Na2O\", \"TiO2\", \"MnO\", \"O\"]\n", + "X = [0.92, 54.57, 8.79, 11.20, 8.45, 12.89, 0.24, 2.24, 1.12, 0.22, 0.64]\n", + "\n", + "db = \"mpe\" # cross-lithology metapelite + metabasite solution models\n", + "dataset = 636\n", + "sys_in = \"mol\"\n", + "\n", + "calc = MAGEMinPTGridCalculator(db=db, dataset=dataset)\n", + "calc.setup_bulk_composition(Xoxides, X, sys_in=sys_in)" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## 3. P–T grid\n", + "\n", + "A grid spanning 15–40 kbar and 550–850 °C covers the eclogite-facies field where garnet and\n", + "clinopyroxene are both stable for this bulk composition." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "P = np.linspace(5.0, 20.0, 11) # kbar\n", + "T = np.linspace(350.0, 900.0, 13) # °C\n", + "Pgrid, Tgrid = np.meshgrid(P, T, indexing=\"xy\")\n", + "\n", + "out = calc.calculate_grid(Pgrid.ravel(), Tgrid.ravel())\n", + "print(f\"{len(out)} grid points; unique phases: {calc.get_all_unique_phases(out)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## 4. Resolve the clinopyroxene phase name\n", + "\n", + "MAGEMin reports clinopyroxene under different solution-model names depending on composition\n", + "(`dio` diopside, `omph` omphacite, `aug` augite). The name is resolved dynamically rather than\n", + "hard-coded, so the notebook follows whichever clinopyroxene solution is stable." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "CPX_CANDIDATES = [\"dio\", \"omph\", \"aug\", \"jac\"] # clinopyroxene solution names in mpe\n", + "\n", + "unique_phases = calc.get_all_unique_phases(out)\n", + "cpx = next((p for p in CPX_CANDIDATES if p in unique_phases), None)\n", + "if cpx is None:\n", + " raise RuntimeError(\"No clinopyroxene phase found in the grid.\")\n", + "print(\"Clinopyroxene phase resolved to:\", cpx)" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## 5. Extract Fe²⁺/Mg and compute $K_D$\n", + "\n", + "Only **divalent** iron takes part in the Fe²⁺–Mg exchange. For the traditional calibration convention\n", + "used here, garnet Fe is deliberately taken as $b_0[\\text{ox\\_apfu\\_FeO}]$ (the FeO basis, treating\n", + "all garnet Fe as Fe²⁺), whereas cpx Fe²⁺ is taken from the excess-O heuristic split. This mixed\n", + "convention is an intentional calibration assumption. Grid points where a phase\n", + "is absent — or where a solvus yields two coexisting limbs of the same phase, making the pairing\n", + "with the other phase ambiguous — are masked out." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "def extract_fe2_mg(calc, phase, grid_out):\n", + " \"\"\"Return (Fe2+/Mg, X_Ca, X_Mn, mask) arrays over the grid for ``phase``.\n", + "\n", + " Garnet uses ``ox_apfu_FeO`` deliberately: this is the traditional FeO-basis\n", + " convention, treating all garnet Fe as Fe2+. Cpx uses ``Fe2`` from the\n", + " excess-O heuristic split. This mixed convention is an intentional calibration\n", + " assumption, not a claim that garnet Fe is strictly all Fe2+. Grid points where\n", + " the phase is absent, or where a solvus\n", + " gives two coexisting limbs of the same phase (an ambiguous pairing with the\n", + " other phase), are flagged True in ``mask`` and should be excluded before\n", + " forming a distribution coefficient.\n", + "\n", + " Parameters\n", + " ----------\n", + " calc : MAGEMinPTGridCalculator\n", + " Calculator that produced ``grid_out``.\n", + " phase : str\n", + " Phase key (e.g. ``\"g\"``).\n", + " grid_out : list\n", + " Output of :meth:`MAGEMinPTGridCalculator.calculate_grid`.\n", + "\n", + " Returns\n", + " -------\n", + " fe2_mg, x_ca, x_mn, mask : numpy.ndarray\n", + " Arrays over the grid. ``x_ca`` and ``x_mn`` are the garnet-style cation\n", + " fractions Ca/(Ca+Mn+Fe+Mg) and Mn/(Ca+Mn+Fe+Mg), with the Fe\n", + " convention above.\n", + " \"\"\"\n", + " bundles = calc.extract_from_grid(\n", + " phase, oxides=[\"MgO\", \"FeO\", \"CaO\", \"MnO\"], fe_split=True, grid_out=grid_out\n", + " )\n", + " n = len(grid_out)\n", + " if not bundles: # phase never appears in the grid\n", + " nan = np.full(n, np.nan)\n", + " return nan, nan, nan, np.ones(n, dtype=bool)\n", + "\n", + " b0 = bundles[0]\n", + " mg = b0[\"ox_apfu_MgO\"]\n", + " # if phase == 'g':\n", + " # # Traditional garnet FeO-basis convention: all garnet Fe is Fe2+.\n", + " # fe2 = b0[\"ox_apfu_FeO\"]\n", + " # else:\n", + " fe2 = b0[\"Fe2\"]\n", + " ca = b0[\"ox_apfu_CaO\"]\n", + " mn = b0[\"ox_apfu_MnO\"]\n", + "\n", + " denom = ca + mn + fe2 + mg # Ca + Mn + Fe2+ + Mg, in atoms per formula unit\n", + " with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n", + " fe2_mg = fe2 / mg\n", + " x_ca = ca / denom\n", + " x_mn = mn / denom\n", + "\n", + " absent = np.isnan(fe2) | np.isnan(mg) | (mg <= 0.0)\n", + " if len(bundles) > 1:\n", + " ambiguous = ~np.isnan(bundles[1][\"Fe2\"]) # a second limb exists -> ambiguous pairing\n", + " else:\n", + " ambiguous = np.zeros(n, dtype=bool)\n", + " return fe2_mg, x_ca, x_mn, absent | ambiguous" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "fe2_mg_g, x_ca_g, x_mn_g, mask_g = extract_fe2_mg(calc, \"g\", out)\n", + "fe2_mg_c, x_ca_c, x_mn_c, mask_c = extract_fe2_mg(calc, cpx, out)\n", + "\n", + "valid = ~mask_g & ~mask_c\n", + "\n", + "Kd = np.full_like(fe2_mg_g, np.nan)\n", + "Kd[valid] = fe2_mg_g[valid] / fe2_mg_c[valid]\n", + "\n", + "print(f\"K_D defined at {int(valid.sum())} / {valid.size} grid points\")\n", + "print(f\"K_D range: {np.nanmin(Kd):.2f} – {np.nanmax(Kd):.2f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## 6. Map $K_D$ over the P–T grid\n", + "\n", + "$K_D$ decreases smoothly with rising temperature — the signature of the Fe²⁺–Mg exchange\n", + "equilibrium — and is only weakly dependent on pressure, which is why a pressure-independent\n", + "calibration such as Mysen & Heier (1972) can work at all." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6.5, 5))\n", + "Kd_map = Kd.reshape(Tgrid.shape)\n", + "\n", + "cs = ax.contourf(Tgrid, Pgrid, Kd_map, levels=14, cmap=\"viridis\")\n", + "ax.contour(Tgrid, Pgrid, Kd_map, levels=14, colors=\"k\", linewidths=0.4)\n", + "cbar = fig.colorbar(cs, ax=ax, label=r\"$K_D = (Fe^{2+}/Mg)^g \\,/\\, (Fe^{2+}/Mg)^{cpx}$\")\n", + "\n", + "ax.set_xlabel(\"Temperature (°C)\")\n", + "ax.set_ylabel(\"Pressure (kbar)\")\n", + "ax.set_title(\"Fe²⁺–Mg distribution coefficient, S10 eclogite (mpe)\")\n", + "ax.grid(ls=\":\", alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "## 7. Thermometer calibrations\n", + "\n", + "Five calibrations of the Fe²⁺–Mg exchange are implemented as pure-NumPy helpers. All return °C.\n", + "\n", + "**Mysen & Heier (1972)** — a pressure- and composition-independent calibration\n", + "(following Banno, 1970):\n", + "\n", + "$$\n", + "T\\,(\\mathrm{K}) = \\frac{2475}{\\ln K_D + 0.781}\n", + "$$\n", + "\n", + "**Ganguly (1979)** — includes pressure and garnet grossular content. This notebook retains the\n", + "piecewise coefficients tabulated by Yavuz & Yıldırım (2020), rather than silently attributing\n", + "that printed piecewise form to the original Ganguly equation.\n", + "in the piecewise form tabulated by Yavuz & Yıldırım (2020):\n", + "\n", + "$$\n", + "T\\,(°\\mathrm{C}) = \\frac{4100 + 1586\\,X_\\mathrm{Ca}^{g} + 11.07\\,P}{\\ln K_D + 2.40} - 273.15\n", + "\\quad (T \\geq 1060\\,°\\mathrm{C})\n", + "$$\n", + "\n", + "$$\n", + "T\\,(°\\mathrm{C}) = \\frac{4801 + 1586\\,X_\\mathrm{Ca}^{g} + 11.07\\,P}{\\ln K_D + 2.93} - 273.15\n", + "\\quad (T \\leq 1060\\,°\\mathrm{C})\n", + "$$\n", + "\n", + "where $X_\\mathrm{Ca}^{g} = \\mathrm{Ca}/(\\mathrm{Ca} + \\mathrm{Mn} + \\mathrm{Fe^{2+}} + \\mathrm{Mg})$\n", + "in garnet and $P$ is in kbar.\n", + "\n", + "**Krogh Ravna (2000)** uses $P_{GPa}=P_{kbar}/10$, $X_{Ca}=Ca/(Ca+Mn+Fe+Mg)$,\n", + "$X_{Mn}=Mn/(Ca+Mn+Fe+Mg)$, and $X_{Mg\\#}=Mg/(Mg+Fe)$.\n", + "\n", + "**Ellis & Green (1979)** uses $X_{Ca,binary}=Ca/(Ca+Fe+Mg)$, deliberately excluding Mn;\n", + "this differs from the garnet site $X_{Ca}$ used by Krogh Ravna and Ganguly.\n", + "\n", + "**Räheim & Green (1974)** is pressure-dependent but composition-independent." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "def T_mysen_heier_1972(Kd):\n", + " \"\"\"Mysen & Heier (1972) garnet–cpx Fe2+–Mg thermometer, in °C.\n", + "\n", + " T(K) = 2475 / (ln Kd + 0.781)\n", + " \"\"\"\n", + " lnKd = np.log(Kd)\n", + " return 2475.0 / (lnKd + 0.781) - 273.15\n", + "\n", + "\n", + "def T_ganguly_1979(Kd, P_kbar, X_Ca_grt):\n", + " \"\"\"Ganguly (1979) garnet–cpx Fe2+–Mg thermometer, in °C.\n", + "\n", + " Piecewise form tabulated by Yavuz & Yıldırım (2020, WinGrt):\n", + " T >= 1060 °C: T = (4100 + 1586·X_Ca + 11.07·P) / (ln Kd + 2.40) - 273.15\n", + " T <= 1060 °C: T = (4801 + 1586·X_Ca + 11.07·P) / (ln Kd + 2.93) - 273.15\n", + " \"\"\"\n", + " lnKd = np.log(Kd)\n", + " P = np.asarray(P_kbar, dtype=float)\n", + " XCa = np.asarray(X_Ca_grt, dtype=float)\n", + " T_hi = (4100.0 + 1586.0 * XCa + 11.07 * P) / (lnKd + 2.40) - 273.15\n", + " T_lo = (4801.0 + 1586.0 * XCa + 11.07 * P) / (lnKd + 2.93) - 273.15\n", + " return np.where(T_hi >= 1060.0, T_hi, T_lo)\n", + "\n", + "\n", + "def T_krogh_ravna_2000(Kd, P_kbar, X_Ca, X_Mn, X_Mg_number):\n", + " \"\"\"Krogh Ravna (2000) garnet-cpx thermometer, in °C.\"\"\"\n", + " lnKd = np.log(Kd)\n", + " P_GPa = np.asarray(P_kbar, dtype=float) / 10.0\n", + " numerator = (1939.9 + 3270.0 * X_Ca - 1396.0 * X_Ca**2\n", + " + 3319.0 * X_Mn + 3535.0 * X_Mn**2\n", + " + 1105.0 * X_Mg_number - 3561.0 * X_Mg_number**2\n", + " + 2324.0 * X_Mg_number**3 + 169.4 * P_GPa)\n", + " return numerator / (lnKd + 1.223) - 273.15\n", + "\n", + "\n", + "def T_ellis_green_1979(Kd, P_kbar, X_Ca_binary):\n", + " \"\"\"Ellis & Green (1979) garnet-cpx thermometer, in °C.\"\"\"\n", + " return (3104.0 * X_Ca_binary + 3030.0 + 10.86 * P_kbar) / (np.log(Kd) + 1.9034) - 273.15\n", + "\n", + "\n", + "def T_raheim_green_1974(Kd, P_kbar):\n", + " \"\"\"Räheim & Green (1974) garnet-cpx thermometer, in °C.\"\"\"\n", + " return (3686.0 + 28.35 * P_kbar) / (np.log(Kd) + 2.33) - 273.15" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "## 8. Apply the thermometers and compare with the true temperature\n", + "\n", + "Because every grid point has a known temperature, we can feed the modelled $K_D$ back into each\n", + "thermometer and see how well it recovers the input. Points lie above the 1:1 line when a\n", + "calibration overestimates the temperature." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "P_flat = Pgrid.ravel()\n", + "T_true = Tgrid.ravel()\n", + "\n", + "T_MH = T_mysen_heier_1972(Kd)\n", + "T_G = T_ganguly_1979(Kd, P_flat, x_ca_g)\n", + "X_mg_number_g = fe2_mg_g / (1.0 + fe2_mg_g)\n", + "X_ca_binary_g = x_ca_g / (1.0 - x_mn_g)\n", + "T_KR = T_krogh_ravna_2000(Kd, P_flat, x_ca_g, x_mn_g, X_mg_number_g)\n", + "T_EG = T_ellis_green_1979(Kd, P_flat, X_ca_binary_g)\n", + "T_RG = T_raheim_green_1974(Kd, P_flat)\n", + "thermometers = {\n", + " \"Mysen & Heier (1972)\": T_MH,\n", + " \"Ganguly (1979; Yavuz & Yıldırım piecewise)\": T_G,\n", + " \"Krogh Ravna (2000)\": T_KR,\n", + " \"Ellis & Green (1979)\": T_EG,\n", + " \"Räheim & Green (1974)\": T_RG,\n", + "}\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 5.5))\n", + "ax.plot([Tgrid.min(), Tgrid.max()], [Tgrid.min(), Tgrid.max()], \"k--\", lw=1, label=\"1:1 (true T)\")\n", + "markers = [\"o\", \"^\", \"s\", \"D\", \"v\"]\n", + "for (name, estimate), marker in zip(thermometers.items(), markers):\n", + " ax.scatter(T_true[valid], estimate[valid], s=16, alpha=0.65, marker=marker, label=name)\n", + "ax.set_xlabel(\"True grid temperature (°C)\")\n", + "ax.set_ylabel(\"Thermometer temperature (°C)\")\n", + "ax.legend()\n", + "ax.grid(ls=\":\", alpha=0.4)\n", + "plt.show()\n", + "\n", + "print(\"\\nThermometer offsets on valid masked points:\")\n", + "print(f\"{'Calibration':42s} {'Mean':>8s} {'Median':>8s} {'RMSE':>8s}\")\n", + "for name, estimate in thermometers.items():\n", + " offset = estimate[valid] - T_true[valid]\n", + " print(f\"{name:42s} {np.mean(offset):+8.1f} {np.median(offset):+8.1f} {np.sqrt(np.mean(offset**2)):8.1f}\")\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 4))\n", + "names = list(thermometers)\n", + "offsets = [np.mean(thermometers[name][valid] - T_true[valid]) for name in names]\n", + "ax.barh(names, offsets, color=plt.get_cmap(\"viridis\")(np.linspace(0.15, 0.85, len(names))))\n", + "ax.axvline(0.0, color=\"k\", lw=0.8)\n", + "ax.set_xlabel(\"Mean thermometer offset (°C)\")\n", + "ax.set_title(\"Calibration offsets on valid grid points\")\n", + "ax.grid(axis=\"x\", ls=\":\", alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6.5, 5))\n", + "Kd_map = Kd.reshape(Tgrid.shape)\n", + "T_RG_grid = T_RG.reshape(Tgrid.shape)\n", + "\n", + "cs = ax.contourf(Tgrid, Pgrid, T_RG_grid - Tgrid, levels=14, cmap=\"viridis\")\n", + "# ax.contour(Tgrid, Pgrid, Kd_map, levels=14, colors=\"k\", linewidths=0.4)\n", + "cbar = fig.colorbar(cs, ax=ax, label=r\"$\\Delta$T (T$_{RG}$ $-$ T$_{MAGEMin}$)\")\n", + "\n", + "ax.set_xlabel(\"Temperature (°C)\")\n", + "ax.set_ylabel(\"Pressure (kbar)\")\n", + "ax.set_title(\"Fe²⁺–Mg distribution coefficient, S10 eclogite (mpe)\")\n", + "ax.grid(ls=\":\", alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6.5, 5))\n", + "Kd_map = Kd.reshape(Tgrid.shape)\n", + "T_EG_grid = T_EG.reshape(Tgrid.shape)\n", + "\n", + "cs = ax.contourf(Tgrid, Pgrid, T_EG_grid - Tgrid, levels=14, cmap=\"viridis\")\n", + "# ax.contour(Tgrid, Pgrid, Kd_map, levels=14, colors=\"k\", linewidths=0.4)\n", + "cbar = fig.colorbar(cs, ax=ax, label=r\"$\\Delta$T (T$_{RG}$ $-$ T$_{MAGEMin}$)\")\n", + "\n", + "ax.set_xlabel(\"Temperature (°C)\")\n", + "ax.set_ylabel(\"Pressure (kbar)\")\n", + "ax.set_title(\"Fe²⁺–Mg distribution coefficient, S10 eclogite (mpe)\")\n", + "ax.grid(ls=\":\", alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "## 9. Why divalent iron — and why the calibrations disagree\n", + "\n", + "Clinopyroxene can carry substantial Fe³⁺ (in acmite/jadeite-type components), so its Fe²⁺\n", + "value is taken from the excess-O heuristic split. Garnet is intentionally different here: the\n", + "traditional thermometer convention uses $b_0[\\text{ox\\_apfu\\_FeO}]$, i.e. garnet Fe on an FeO\n", + "basis with all garnet Fe treated as Fe²⁺. The resulting garnet-FeO / cpx-heuristic-Fe²⁺ pairing\n", + "is a mixed convention adopted for calibration comparison, not a claim that both phases have\n", + "the same redox treatment." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "from phasetools import get_phase_fe_split\n", + "\n", + "i0 = int(np.flatnonzero(valid)[0]) # first valid grid point\n", + "o0 = out[i0]\n", + "cpx_split = get_phase_fe_split(o0, cpx)\n", + "grt_feo_basis = calc.extract_from_grid(\"g\", oxides=[\"FeO\"], grid_out=[o0])[0][\"ox_apfu_FeO\"][0]\n", + "\n", + "print(f\"At {T_true[i0]:.0f} °C, {P_flat[i0]:.0f} kbar:\")\n", + "print(f\" garnet : FeO-basis Fe = {grt_feo_basis:.3f} (all treated as Fe2+ by convention)\")\n", + "print(f\" {cpx:10s}: Fe2+ = {cpx_split['Fe2']:.3f}, Fe3+ = {cpx_split['Fe3']:.3f} \"\n", + " f\"(Fe3+/FeOt = {cpx_split['Fe3']/(cpx_split['Fe2']+cpx_split['Fe3']):.0%})\")" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "The calibrations disagree because they encode different pressure and composition corrections. Mysen & Heier\n", + "(1972) is pressure- and composition-independent; Ganguly (1979) is retained here only with its\n", + "Yavuz & Yıldırım (2020) piecewise coefficients explicitly labelled; Krogh Ravna (2000) includes\n", + "Mn, Ca, Mg#, and pressure; Ellis & Green (1979) uses binary garnet Ca and pressure; and Räheim &\n", + "Green (1974) uses pressure alone. The offsets therefore illustrate calibration uncertainty, compounded\n", + "by the intentional mixed garnet FeO / cpx heuristic Fe²⁺ convention. These equations should not be\n", + "extrapolated beyond their experimental compositional and P–T ranges." + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "---\n", + "## References\n", + "\n", + "- Mysen, B. O., & Heier, K. S. (1972). Petrogenesis of eclogites in high grade metamorphic\n", + " gneisses, exemplified by the Hareidland eclogite, western Norway. *Contributions to Mineralogy\n", + " and Petrology*, 36(1), 73–94. https://doi.org/10.1007/BF00372836\n", + "- Ganguly, J. (1979). Garnet and clinopyroxene solid solutions, and geothermometry based on\n", + " Fe–Mg distribution coefficient. *Geochimica et Cosmochimica Acta*, 43(7), 1021–1029.\n", + " https://doi.org/10.1016/0016-7037(79)90091-7\n", + "- Yavuz, F., & Yıldırım, D. K. (2020). WinGrt, a Windows program for garnet supergroup minerals.\n", + " *Journal of Geosciences*, 65(2), 71–95. https://doi.org/10.3190/jgeosci.303 (source of the\n", + " Ganguly 1979 piecewise coefficients)\n", + "- Johnson, C. A., Bohlen, S. R., & Essene, E. J. (1983). An evaluation of garnet-clinopyroxene\n", + " geothermometry in granulites. *Contributions to Mineralogy and Petrology*, 84(2–3), 191–198.\n", + " https://doi.org/10.1007/BF00371285\n", + "- Krogh Ravna, E. (2000). The garnet–clinopyroxene Fe²⁺–Mg geothermometer: an updated\n", + " calibration. *Journal of Metamorphic Geology*, 18(2), 211–219.\n", + " https://doi.org/10.1046/j.1525-1314.2000.00247.x\n", + "- Ellis, D. J., & Green, D. H. (1979). An experimental study of the effect of Ca upon garnet–\n", + " clinopyroxene Fe–Mg exchange equilibria. *Contributions to Mineralogy and Petrology*, 71, 13–22.\n", + " https://doi.org/10.1007/BF00371878\n", + "- Räheim, A., & Green, D. H. (1974). Experimental determination of the temperature and pressure\n", + " dependence of the Fe–Mg partition coefficient for coexisting garnet and clinopyroxene.\n", + " *Contributions to Mineralogy and Petrology*, 48, 179–203.\n", + " https://doi.org/10.1007/BF00392328\n", + "- Thomas, J. B., & Rana, S. (2024). Garnet–clinopyroxene thermometry (as discussed alongside\n", + " Yavuz & Yıldırım). Consult the source report for the exact calibration context; this notebook\n", + " does not implement an unverified Powell equation.\n", + " (The Ganguly piecewise coefficients above are attributed explicitly to Yavuz & Yıldırım.)" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "uw3-dev", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Tutorials/Garnet/README.md b/Tutorials/Garnet/README.md index 55fd92b..d8a391a 100644 --- a/Tutorials/Garnet/README.md +++ b/Tutorials/Garnet/README.md @@ -21,5 +21,14 @@ This directory contains tutorials focused on modelling garnet chemistry, growth, - Implementing geothermobarometry by minimising the misfit between measured chemistry and thermodynamic predictions. - Recovering P-T trajectories from zoned garnet crystals. +### [3. Garnet-Cpx Thermometry](./3-Garnet_cpx_thermometry.ipynb) +**Objective:** Mapping the garnet–clinopyroxene Fe²⁺–Mg distribution coefficient (Kd) and comparing thermometer calibrations. +- **Key Features:** + - Setting up the S10 eclogite bulk composition in the `mpe` database. + - Extracting garnet and clinopyroxene Fe²⁺/Mg from APFU with dynamic cpx phase resolution (`dio`/`omph`/`aug`). + - Masking absent and solvus-ambiguous grid cells. + - Mapping Kd over a P–T grid and applying Mysen & Heier (1972), Ganguly (1979; explicitly labelled Yavuz & Yıldırım piecewise form), Krogh Ravna (2000), Ellis & Green (1979), and Räheim & Green (1974) calibrations. + - Using the intentional traditional mixed Fe convention: garnet Fe from `ox_apfu_FeO` (all treated as Fe²⁺) and cpx Fe²⁺ from the excess-O heuristic split. + --- **Units Note:** All tutorials use `kbar` for pressure and `°C` for temperature. diff --git a/src/phasetools/__init__.py b/src/phasetools/__init__.py index b439b18..41c9a1c 100644 --- a/src/phasetools/__init__.py +++ b/src/phasetools/__init__.py @@ -9,7 +9,7 @@ # Expose Public API from .core.base import MAGEMinBase from .core.engine import single_point_minimization_with_conversion -from .core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, calculate_kd_fe_mg +from .core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, get_phase_mg2_number, get_phase_fe_split, calculate_kd_fe_mg from .calculators.garnet import MAGEMinGarnetCalculator from .calculators.assemblage import MAGEMinAssemblageCalculator @@ -18,5 +18,6 @@ from .calculators.pt_grid import MAGEMinPTGridCalculator from .models.garnet_growth import GarnetGenerator, generate_distribution +from .models.magma_ocean import MagmaOcean from .utils import bulk_rock diff --git a/src/phasetools/calculators/README.md b/src/phasetools/calculators/README.md index 8a9e495..621b6a2 100644 --- a/src/phasetools/calculators/README.md +++ b/src/phasetools/calculators/README.md @@ -11,6 +11,9 @@ Calculators are high-level wrappers designed for specific thermodynamic tasks, s ### `garnet.py` - **`MAGEMinGarnetCalculator`**: Specialized for garnet chemistry. - Includes logic for site-specific cation fractions (X-site) and modelling garnet evolution along P-T paths (optionally with fractionation). +- When `sys_in='wt'`, fractionation uses the weight-basis garnet fraction (not the molar fraction) to stay consistent with the bulk composition units. +- `X_along_path` rows are normalised to sum to 1 regardless of the `sys_in` setting. +- **`fe_basis`** option controls the Fe basis of garnet X-site fractions: `'FeOt'` (default — all Fe treated as Fe²⁺, the standard community convention, e.g. Williams & Grambling 1990; Krogh Ravna 2000) or `'Fe2+'` (stoichiometric ferrous-only split, for garnets with significant Fe³⁺). ### `pt_estimation.py` - **`PhasePTEstimator`**: A geothermobarometry engine. @@ -21,5 +24,6 @@ Calculators are high-level wrappers designed for specific thermodynamic tasks, s - Generates boolean masks to identify P-T regions where specific minerals coexist. ### `phase_search.py` -- **`PhaseSearchCalculator`**: Uses root-finding algorithms to locate phase boundaries. +- **`PhaseFunctions`**: Uses root-finding algorithms to locate phase boundaries. - Precise determination of solidus, liquidus, or specific phase appearance/disappearance points. +- `fractionate_phase` handles both solution phases (`SS_vec`) and pure phases (`PP_vec`) transparently. When a phase appears multiple times (e.g. solvus), only the first instance is removed. Output is always normalised to sum to 1. diff --git a/src/phasetools/calculators/garnet.py b/src/phasetools/calculators/garnet.py index 4565236..5236ab6 100644 --- a/src/phasetools/calculators/garnet.py +++ b/src/phasetools/calculators/garnet.py @@ -2,61 +2,107 @@ import sys from .pt_grid import MAGEMinPTGridCalculator from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu -from ..utils.bulk_rock import atomic_mass_dict, convert_mol_percent_to_wt_percent +from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac from phasetools import MAGEMin_C from juliacall import Main as jl, convert as jlconvert class MAGEMinGarnetCalculator(MAGEMinPTGridCalculator): """High-level wrappers for garnet-focused MAGEMin calculations.""" - def __init__(self, db="ig", dataset=636, verbose=False): + + def __init__(self, db="ig", dataset=636, verbose=False, fe_basis="FeOt"): + """ + Parameters + ---------- + db : str, default="ig" + Thermodynamic database label. + dataset : int, default=636 + Thermodynamic dataset version. + verbose : bool, default=False + Whether to print progress information. + fe_basis : str, default="FeOt" + Iron basis used for garnet X-site fractions: + + * ``'FeOt'`` -- all iron treated as Fe2+ (total iron, + ``FeO + 2*Fe2O3`` APFU) placed on the divalent site. This is + the standard community convention for garnet end-members and + X-site fractions (e.g. Williams & Grambling 1990; Krogh Ravna + 2000) and is consistent with the pyralspite garnet model used + by MAGEMin and Holland-Powell-type databases. + * ``'Fe2+'`` -- only the stoichiometrically estimated ferrous + iron (``_extract_fe_split_from_apfu``) is placed on the + divalent site, excluding Fe3+. Use only when garnet Fe3+ is + known to be significant (oxidised eclogites, skarns) or + measured directly (XANES, Mössbauer). + + Case-insensitive. + """ super().__init__(db, dataset, verbose) + basis = str(fe_basis).casefold() + if basis not in ("feot", "fe2+", "fe2"): + raise ValueError(f"fe_basis must be 'FeOt' or 'Fe2+', got {fe_basis!r}") + self.fe_basis = "fe2+" if basis.startswith("fe2") else "feot" def _extract_garnet_elements_from_oxides(self, out, sys_in): """ Extract garnet Mg-Mn-Fe-Ca cation fractions for the divalent (X) site. - Strictly isolates Fe2+ to ensure X-site fractions sum to 1.0. + + The Fe basis is controlled by ``self.fe_basis``: + + * ``'feot'`` (default) -- all iron treated as Fe2+ (total iron, + ``FeO + 2*Fe2O3`` APFU) placed on the divalent site. The standard + community convention for garnet X-site fractions. + * ``'fe2+'`` -- only the stoichiometrically estimated ferrous iron + is placed on the divalent site, excluding Fe3+. """ if 'g' not in out.ph: return 0.0, 0.0, 0.0, 0.0 - # Isolating divalent iron using the robust Fe-split method - split = self._extract_fe_split_from_apfu(out, 'g') - fe2_moles = split['fe2'] - - elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO']) + if self.fe_basis == 'feot': + # Total iron on the divalent site (all Fe as Fe2+) + elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO', 'FeO', 'Fe2O3']) + fe_moles = elements.get("FeO", 0.0) + 2.0 * elements.get("Fe2O3", 0.0) + else: + # Ferrous iron only, from the excess-oxygen Fe2+/Fe3+ split + elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO']) + fe_moles = self._extract_fe_split_from_apfu(out, 'g')['Fe2'] + mg_moles = elements.get("MgO", 0.0) mn_moles = elements.get("MnO", 0.0) ca_moles = elements.get("CaO", 0.0) # Total atoms in the X-site (should be approx 3.0) - total_x_site_moles = mg_moles + mn_moles + fe2_moles + ca_moles - + total_x_site_moles = mg_moles + mn_moles + fe_moles + ca_moles if total_x_site_moles <= 0: return 0.0, 0.0, 0.0, 0.0 # Normalise to 1.0 (mole fractions of the divalent site) Mg = mg_moles / total_x_site_moles Mn = mn_moles / total_x_site_moles - Fe = fe2_moles / total_x_site_moles + Fe = fe_moles / total_x_site_moles Ca = ca_moles / total_x_site_moles if sys_in.casefold() == 'wt': - wt_percent_list = convert_mol_percent_to_wt_percent( - [Mg, Mn, Fe, Ca], - ["Mg", "Mn", "Fe", "Ca"], + wt_frac = atomic_frac_to_wt_frac( + {"Mg": Mg, "Mn": Mn, "Fe": Fe, "Ca": Ca}, atomic_mass_dict, ) - Mg, Mn, Fe, Ca = [val / 100.0 for val in wt_percent_list] + Mg, Mn, Fe, Ca = wt_frac["Mg"], wt_frac["Mn"], wt_frac["Fe"], wt_frac["Ca"] return Mg, Mn, Fe, Ca def generate_2D_grid_gt_endmembers(self, P, T): - """Compute garnet endmember fractions over a P-T grid.""" + """Compute garnet endmember fractions over a P-T grid. + + Garnet is a single-instance phase, so its single per-instance + bundle is returned directly, preserving the historical key names + (``em_py``, ``em_alm``, ``mol_frac``, ...). If garnet ever + appeared as a solvus, the list of bundles would be returned + instead. + """ self.calculate_grid(P, T) # Automatic discovery of end-members res = self.extract_from_grid("g", end_members='auto') - - return res + return res[0] if len(res) == 1 else res def generate_2D_grid_gt_elements(self, P, T): """Compute garnet element fractions (Mg, Mn, Fe, Ca) over a P-T grid.""" @@ -121,7 +167,61 @@ def _gt_single_point_from_jl(self, P, T, X_jl, Xoxides_jl, sys_in, rm_list=None) return gt_frac, gt_wt, gt_vol, Mg, Mn, Fe, Ca, out def gt_along_path(self, P, T, fractionate=False, normalise_start=True): - """Calculate garnet fractions and element chemistry along a P-T path.""" + """Calculate garnet fractions and element chemistry along a P-T path. + + Parameters + ---------- + P : array-like + Pressure values along the path (kbar). + T : array-like + Temperature values along the path (°C). + fractionate : bool, default=False + If True, fractionate garnet from the bulk composition as it grows. + normalise_start : bool, default=True + Controls how the first P-T point is treated: + + * ``True`` — The first P-T point is treated as a nucleation + barrier: garnet volume starts at zero and only new growth is + modelled. The initial garnet fraction (if any) is used as a + baseline for measuring incremental growth but is **not** removed + from the reactive bulk. Use when the path starts outside the + garnet stability field or at the nucleation threshold. + + * ``False`` — The first P-T point has an initial garnet volume + (overstepped nucleation). That fraction is removed from the + bulk at step 0, and subsequent growth is measured relative to + it. Use when the path starts well inside the garnet stability + field. + + Returns + ------- + gt_mol_frac : numpy.ndarray + Garnet molar fraction at each P-T point. + gt_wt_frac : numpy.ndarray + Garnet weight fraction at each P-T point. + gt_vol_frac : numpy.ndarray + Garnet volume fraction at each P-T point. + Mgi : numpy.ndarray + Garnet X-site Mg fraction at each P-T point. + Mni : numpy.ndarray + Garnet X-site Mn fraction at each P-T point. + Fei : numpy.ndarray + Garnet X-site Fe fraction at each P-T point. The Fe basis + (total Fe as Fe2+ by default, or ferrous-only) is set by the + ``fe_basis`` argument passed to :meth:`__init__`. + Cai : numpy.ndarray + Garnet X-site Ca fraction at each P-T point. + X_along_path : numpy.ndarray + Bulk composition after each step's fractionation. Each row is + normalised to sum to 1. Row ``i`` is the bulk **after** step + ``i``'s fractionation has been applied. + + Notes + ----- + Fractionation uses the **current-step** garnet composition (not the + growth-increment composition) — a first-order approximation valid for + small P-T steps. + """ from .phase_search import PhaseFunctions X = self.X @@ -138,6 +238,7 @@ def gt_along_path(self, P, T, fractionate=False, normalise_start=True): X_along_path = np.zeros(shape=(n_points, len(self._Xoxides_py)) ) gt_frac_max_previous = 0. + gt_wt_max_previous = 0. phase_functions = PhaseFunctions(db=self.db, dataset=self.dataset, verbose=self.verbose) if fractionate else None if phase_functions: # Sync standardised state to the helper instance @@ -152,19 +253,28 @@ def gt_along_path(self, P, T, fractionate=False, normalise_start=True): gt_wt_frac[i] = gt_wt gt_vol_frac[i] = gt_vol + # Select the fraction basis consistent with sys_in + if self.sys_in.casefold() == 'wt': + gt_frac_for_fractionation = gt_wt + gt_frac_max_prev_for_fractionation = gt_wt_max_previous + else: + gt_frac_for_fractionation = gt_frac + gt_frac_max_prev_for_fractionation = gt_frac_max_previous + if phase_functions is not None: if i == 0 and not normalise_start: - if gt_frac > 0: - X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=gt_frac) + if gt_frac_for_fractionation > 0: + X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=gt_frac_for_fractionation) X = jlconvert(jl.Vector[jl.Float64], X_py) elif i > 0: - frac_amount = max(gt_frac - gt_frac_max_previous, 0.0) + frac_amount = max(gt_frac_for_fractionation - gt_frac_max_prev_for_fractionation, 0.0) if frac_amount > 0: X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=frac_amount) X = jlconvert(jl.Vector[jl.Float64], X_py) - X_along_path[i] = np.array(X) + X_along_path[i] = np.array(X) / np.sum(X) gt_frac_max_previous = max(gt_frac_max_previous, gt_frac) + gt_wt_max_previous = max(gt_wt_max_previous, gt_wt) Mgi[i] = Mg Mni[i] = Mn diff --git a/src/phasetools/calculators/phase_search.py b/src/phasetools/calculators/phase_search.py index 54fe8a0..cdfa8d9 100644 --- a/src/phasetools/calculators/phase_search.py +++ b/src/phasetools/calculators/phase_search.py @@ -41,24 +41,72 @@ def liquidus_func(T): return result.root def fractionate_phase(self, phase, out, sys_in, frac_amount=None): - """Perform batch fractionation of a phase from the bulk rock composition.""" + """Perform batch fractionation of a phase from the bulk rock composition. + + Removes a specified fraction of a phase from the current bulk composition + and renormalises the remainder. The output is always normalised to sum to 1 + (molar or weight fractions, depending on ``sys_in``). + + MAGEMin orders ``out.ph`` with solution phases first (``SS_vec``), then + pure phases (``PP_vec``). This method handles both indexing ranges + transparently. + + .. note:: + + When a phase appears multiple times in ``out.ph`` (e.g. two + coexisting pyroxenes on a solvus), only the **first** instance is + removed. Use ``phase_frac()`` to obtain the summed fraction before + calling this method if the total solvus fraction is needed. + + Parameters + ---------- + phase : str + Phase name as it appears in ``out.ph`` (e.g. ``'g'``, ``'liq'``). + out : MAGEMinOutput + Raw output from a MAGEMin single-point minimisation. + sys_in : str + ``'mol'`` or ``'wt'`` — determines which bulk composition and + phase composition vectors are used. + frac_amount : float or None, optional + Fraction of the phase to remove. If ``None``, the full phase + fraction from the minimisation result is used. Values of 0 or + less are treated as a no-op and return the current bulk unchanged. + + Returns + ------- + numpy.ndarray + Renormalised bulk composition after fractionation. + """ if sys_in.casefold() == "wt": current_X = out.bulk_wt else: current_X = out.bulk - + if phase in out.ph: phase_ind = out.ph.index(phase) - if sys_in.casefold() == "wt": - ph_comp = np.array(out.SS_vec[phase_ind].Comp_wt) - ph_frac = out.ph_frac_wt[phase_ind] + # MAGEMin orders: solution phases (SS_vec) first, then pure phases (PP_vec) + if phase_ind < out.n_SS: + if sys_in.casefold() == "wt": + ph_comp = np.array(out.SS_vec[phase_ind].Comp_wt) + ph_frac = out.ph_frac_wt[phase_ind] + else: + ph_comp = np.array(out.SS_vec[phase_ind].Comp) + ph_frac = out.ph_frac[phase_ind] else: - ph_comp = np.array(out.SS_vec[phase_ind].Comp) - ph_frac = out.ph_frac[phase_ind] + pp_ind = phase_ind - out.n_SS + if sys_in.casefold() == "wt": + ph_comp = np.array(out.PP_vec[pp_ind].Comp_wt) + ph_frac = out.ph_frac_wt[phase_ind] + else: + ph_comp = np.array(out.PP_vec[pp_ind].Comp) + ph_frac = out.ph_frac[phase_ind] if frac_amount is None: frac_amount = ph_frac + if frac_amount <= 0: + return np.array(current_X, dtype=float) + if frac_amount >= 1.0: warnings.warn(f"fractionate_phase: requested frac_amount={frac_amount} >= 1.0; skipping.") return np.array(current_X, dtype=float) diff --git a/src/phasetools/calculators/pt_estimation.py b/src/phasetools/calculators/pt_estimation.py index ca7fe64..88227d7 100644 --- a/src/phasetools/calculators/pt_estimation.py +++ b/src/phasetools/calculators/pt_estimation.py @@ -26,8 +26,8 @@ def _get_phase_composition(self, out, phase, components, comp_type='element'): # Handle special components Fe2 and Fe3 first if 'Fe2' in components or 'Fe3' in components: split = self._extract_fe_split_from_apfu(out, phase) - element_map['Fe2'] = split['fe2'] - element_map['Fe3'] = split['fe3'] + element_map['Fe2'] = split['Fe2'] + element_map['Fe3'] = split['Fe3'] for ox, stoichiometry in self._stoich_map.items(): for el, mult in stoichiometry.items(): diff --git a/src/phasetools/calculators/pt_grid.py b/src/phasetools/calculators/pt_grid.py index f5d1f94..7fee7ec 100644 --- a/src/phasetools/calculators/pt_grid.py +++ b/src/phasetools/calculators/pt_grid.py @@ -1,8 +1,9 @@ import numpy as np import sys +from juliacall import Main as jl, convert as jlconvert from ..core.base import MAGEMinBase -from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number -from ..utils.bulk_rock import atomic_mass_dict, convert_mol_percent_to_wt_percent +from ..core.phase_properties import extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, _phase_indices +from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac from phasetools import MAGEMin_C class MAGEMinPTGridCalculator(MAGEMinBase): @@ -20,17 +21,19 @@ def calculate_grid(self, P, T): P = np.atleast_1d(P) T = np.atleast_1d(T) - if P.shape != T.shape: - if P.ndim == 1 and T.ndim == 1: - P_orig, T_orig = P, T - P, T = np.meshgrid(P_orig, T_orig) - P = P.flatten() - T = T.flatten() - else: - raise ValueError(f"P and T must have the same shape or both be 1D. Got {P.shape} and {T.shape}") + if P.ndim == 1 and T.ndim == 1 and P.shape[0] != T.shape[0]: + P_orig, T_orig = P, T + P, T = np.meshgrid(P_orig, T_orig) + P = P.flatten() + T = T.flatten() + elif P.shape != T.shape: + raise ValueError(f"P and T must have the same shape or both be 1D. Got {P.shape} and {T.shape}") + + P_jl = jlconvert(jl.Vector[jl.Float64], P.astype(float)) + T_jl = jlconvert(jl.Vector[jl.Float64], T.astype(float)) out = MAGEMin_C.multi_point_minimization( - P, T, self.data, X=self.X, Xoxides=self.Xoxides, + P_jl, T_jl, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list ) sys.stdout.flush() @@ -77,38 +80,44 @@ def get_phase_endmembers(self, phase, grid_out=None): if out is None: return [] - # If it's a single output object (from single_point_calc) - if not isinstance(out, (list, np.ndarray)): + # If it's a single output object (from single_point_calc) — not a list/array/Vector + if not isinstance(out, (list, np.ndarray)) and hasattr(out, 'ph'): if phase in out.ph: ph_index = out.ph.index(phase) - return [str(n) for n in out.SS_vec[ph_index].emNames] + if ph_index < len(out.SS_vec): + return [str(n) for n in out.SS_vec[ph_index].emNames] return [] - # If it's a grid + # If it's a grid (list, numpy array, or Julia Vector) for o in out: if phase in o.ph: ph_index = o.ph.index(phase) - return [str(n) for n in o.SS_vec[ph_index].emNames] + if ph_index < len(o.SS_vec): + return [str(n) for n in o.SS_vec[ph_index].emNames] + return [] return [] def _extract_cations_from_apfu(self, out, phase, cations, sys_in): - """Internal: Extract cation ratios (e.g., XMg, XFe) for a specific phase.""" + """Internal: per-instance cation ratios (e.g., XMg, XFe). + + Returns ``cat_`` arrays, one entry per instance of the phase. + """ ox_to_query = ['MgO', 'MnO', 'CaO', 'FeO', 'Fe2O3', 'Fe', 'O'] - apfu = get_oxide_apfu(out, phase, ox_to_query) - - mg = apfu.get("MgO", 0.0) - mn = apfu.get("MnO", 0.0) - ca = apfu.get("CaO", 0.0) - - feo = apfu.get("FeO", 0.0) - fe2o3 = apfu.get("Fe2O3", 0.0) - fe_metal = apfu.get("Fe", 0.0) - atomic_o = apfu.get("O", 0.0) + apfu = get_oxide_apfu(out, phase, ox_to_query, instance='all') + + mg = np.asarray(apfu.get("MgO", 0.0), dtype=float) + mn = np.asarray(apfu.get("MnO", 0.0), dtype=float) + ca = np.asarray(apfu.get("CaO", 0.0), dtype=float) + + feo = np.asarray(apfu.get("FeO", 0.0), dtype=float) + fe2o3 = np.asarray(apfu.get("Fe2O3", 0.0), dtype=float) + fe_metal = np.asarray(apfu.get("Fe", 0.0), dtype=float) + atomic_o = np.asarray(apfu.get("O", 0.0), dtype=float) # Calculate total Fe as FeO equivalent (molar atoms) - if atomic_o > 0: + if np.any(atomic_o > 0): # MAGEMin O-basis (ig, mp) or sb24 basis - if fe_metal > 0: + if np.any(fe_metal > 0): fe = fe_metal else: fe = feo @@ -117,21 +126,45 @@ def _extract_cations_from_apfu(self, out, phase, cations, sys_in): fe = feo + 2.0 * fe2o3 total = mg + mn + fe + ca - if total <= 0: - return {c: 0.0 for c in cations} + total_safe = np.where(total > 0, total, np.nan) + with np.errstate(divide='ignore', invalid='ignore'): + vals = {"Mg": mg/total_safe, "Mn": mn/total_safe, + "Fe": fe/total_safe, "Ca": ca/total_safe} - vals = {"Mg": mg/total, "Mn": mn/total, "Fe": fe/total, "Ca": ca/total} - if sys_in.casefold() == 'wt': - keys = list(vals.keys()) - raw_vals = [vals[k] for k in keys] - wt_percents = convert_mol_percent_to_wt_percent(raw_vals, keys, atomic_mass_dict) - vals = {k: v/100.0 for k, v in zip(keys, wt_percents)} + vals = atomic_frac_to_wt_frac(vals, atomic_mass_dict) - return {f"cat_{c}": vals.get(c, 0.0) for c in cations} + return {f"cat_{c}": np.asarray(vals[c], dtype=float) for c in cations} def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, grid_out=None): - """Extract phase properties from a previously calculated grid.""" + """Extract phase properties from a previously calculated grid. + + MAGEMin reports coexisting solvus limbs as repeated entries in + ``out.ph`` (e.g. two clinopyroxenes ``dio``). Returns a list + with one bundle per instance, indexed ``res[0]``, ``res[1]`` + ...: + + ``res[0]`` -- first (or only) instance; every key is an array + over the grid points (scalars for ``single_point_calc``). + + The key shape is identical for every phase -- a phase with a + single instance simply has a one-element list, always at + ``res[0]``. Sum the per-instance ``mol_frac`` arrays yourself + if you want the total. + + Each bundle key is NaN where the phase (or that instance) is + absent at a grid point. If the phase never appears in the grid, + the returned list is empty. + + Notes + ----- + - Bundle ``res[k]`` is the k-th occurrence of the phase in that + point's ``out.ph``; solvus branch ordering may swap between + grid points, so track both limbs when plotting isopleths. + - ``end_members='auto'`` discovers end-members from the first + occurrence of the phase; solvus limbs normally share the same + solution model and endmember set. + """ out = grid_out if grid_out is not None else self.last_grid_out if out is None: raise ValueError("No grid results found. Run calculate_grid first or provide grid_out.") @@ -141,50 +174,77 @@ def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None end_members = self.get_phase_endmembers(phase, out) P_len = len(out) - results = { - "mol_frac": np.zeros(P_len), "wt_frac": np.zeros(P_len), "vol_frac": np.zeros(P_len), - } - + + # n_inst = max occurrences of the phase across the whole grid + n_inst = 0 + for o in out: + n_inst = max(n_inst, len(_phase_indices(o, phase, 'all'))) + + instances = [{} for _ in range(n_inst)] + + def _precreate(prefix, names): + for k in range(n_inst): + for n in names: + instances[k][f"{prefix}{n}"] = np.full(P_len, np.nan) + + for k in range(n_inst): + instances[k]["mol_frac"] = np.full(P_len, np.nan) + instances[k]["wt_frac"] = np.full(P_len, np.nan) + instances[k]["vol_frac"] = np.full(P_len, np.nan) if end_members: - for em in end_members: results[f"em_{em}"] = np.zeros(P_len) + _precreate("em_", end_members) if oxides: - for ox in oxides: results[f"ox_apfu_{ox}"] = np.zeros(P_len) + _precreate("ox_apfu_", oxides) if chemistry: - for ox in chemistry: results[f"chem_{ox}"] = np.zeros(P_len) + _precreate("chem_", chemistry) if cations: - for c in cations: results[f"cat_{c}"] = np.zeros(P_len) + _precreate("cat_", cations) if mg_number: - results["mg_number"] = np.zeros(P_len) + for k in range(n_inst): + instances[k]["Mg_number"] = np.full(P_len, np.nan) if fe_split: - results["fe2"] = np.zeros(P_len) - results["fe3"] = np.zeros(P_len) + for k in range(n_inst): + instances[k]["Fe2"] = np.full(P_len, np.nan) + instances[k]["Fe3"] = np.full(P_len, np.nan) for i in range(P_len): - if phase in out[i].ph: - results["mol_frac"][i] = phase_frac(phase, out[i], 'mol') - results["wt_frac"][i] = phase_frac(phase, out[i], 'wt') - results["vol_frac"][i] = phase_frac(phase, out[i], 'vol') - - if end_members: - for em in end_members: - results[f"em_{em}"][i] = extract_end_member(phase, out[i], em, self.sys_in) - if oxides: - apfu = get_oxide_apfu(out[i], phase, oxides) - for ox in oxides: results[f"ox_apfu_{ox}"][i] = apfu.get(ox, 0.0) - if chemistry: - chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in) - for ox in chemistry: results[f"chem_{ox}"][i] = chem.get(ox, 0.0) - if cations: - cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in) - for c in cations: results[f"cat_{c}"][i] = cat_vals[f"cat_{c}"] - if mg_number: - results["mg_number"][i] = get_phase_mg_number(out[i], phase) - if fe_split: - split = self._extract_fe_split_from_apfu(out[i], phase) - results["fe2"][i] = split["fe2"] - results["fe3"][i] = split["fe3"] - - return results + if phase not in out[i].ph: + continue + n_idx = _phase_indices(out[i], phase, 'all') + for k, j in enumerate(n_idx): + instances[k]["mol_frac"][i] = float(out[i].ph_frac[j]) + instances[k]["wt_frac"][i] = float(out[i].ph_frac_wt[j]) + instances[k]["vol_frac"][i] = float(out[i].ph_frac_vol[j]) + # instances k >= len(n_idx) stay NaN + + def _fill(key, value): + vals = np.atleast_1d(np.asarray(value, dtype=float)) + n = len(vals) + for k in range(n_inst): + if k < n: + instances[k][key][i] = vals[k] + # else stays NaN + + if end_members: + for em in end_members: + _fill(f"em_{em}", extract_end_member(phase, out[i], em, self.sys_in, instance='all')) + if oxides: + apfu = get_oxide_apfu(out[i], phase, oxides, instance='all') + for ox in oxides: _fill(f"ox_apfu_{ox}", apfu.get(ox, np.zeros(0))) + if chemistry: + chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in, instance='all') + for ox in chemistry: _fill(f"chem_{ox}", chem.get(ox, np.zeros(0))) + if cations: + cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in) + for c in cations: _fill(f"cat_{c}", cat_vals[f"cat_{c}"]) + if mg_number: + _fill("Mg_number", get_phase_mg_number(out[i], phase, instance='all')) + if fe_split: + split = self._extract_fe_split_from_apfu(out[i], phase, instance='all') + _fill("Fe2", split["Fe2"]) + _fill("Fe3", split["Fe3"]) + + return instances def generate_2D_grid(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False): """Convenience wrapper.""" @@ -192,31 +252,47 @@ def generate_2D_grid(self, P, T, phase, end_members=None, oxides=None, chemistry return self.extract_from_grid(phase, end_members, oxides, chemistry, cations, mg_number, fe_split) def single_point_calc(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False): - """Single-point calculation.""" + """Single-point calculation. + + Returns ``(bundles, out)`` where ``bundles`` is a list with one + bundle per instance of the phase (``bundles[0]``, ``bundles[1]``, + ...), keyed like the grid-level ``extract_from_grid`` with scalar + values. The list is empty when the phase is not present at this + P-T. + """ out = MAGEMin_C.single_point_minimization(P, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list) sys.stdout.flush() - results = {"mol_frac": 0.0, "wt_frac": 0.0, "vol_frac": 0.0, "present": False} + instances = [] if phase in out.ph: - results["present"] = True - results["mol_frac"] = phase_frac(phase, out, 'mol') - results["wt_frac"] = phase_frac(phase, out, 'wt') - results["vol_frac"] = phase_frac(phase, out, 'vol') + n_idx = _phase_indices(out, phase, 'all') + instances = [{"mol_frac": 0.0, "wt_frac": 0.0, "vol_frac": 0.0} + for _ in n_idx] + for k, j in enumerate(n_idx): + instances[k]["mol_frac"] = float(out.ph_frac[j]) + instances[k]["wt_frac"] = float(out.ph_frac_wt[j]) + instances[k]["vol_frac"] = float(out.ph_frac_vol[j]) + + def _store(key, value): + vals = np.atleast_1d(np.asarray(value, dtype=float)) + for k in range(len(vals)): + instances[k][key] = vals[k] if end_members: - for em in end_members: results[f"em_{em}"] = extract_end_member(phase, out, em, self.sys_in) + for em in end_members: + _store(f"em_{em}", extract_end_member(phase, out, em, self.sys_in, instance='all')) if oxides: - apfu = get_oxide_apfu(out, phase, oxides) - for ox in oxides: results[f"ox_apfu_{ox}"] = apfu.get(ox, 0.0) + apfu = get_oxide_apfu(out, phase, oxides, instance='all') + for ox in oxides: _store(f"ox_apfu_{ox}", apfu.get(ox, np.zeros(0))) if chemistry: - chem = get_phase_chemistry(out, phase, chemistry, self.sys_in) - for ox in chemistry: results[f"chem_{ox}"] = chem.get(ox, 0.0) + chem = get_phase_chemistry(out, phase, chemistry, self.sys_in, instance='all') + for ox in chemistry: _store(f"chem_{ox}", chem.get(ox, np.zeros(0))) if cations: cat_vals = self._extract_cations_from_apfu(out, phase, cations, self.sys_in) - for c in cations: results[f"cat_{c}"] = cat_vals[f"cat_{c}"] + for c in cations: _store(f"cat_{c}", cat_vals[f"cat_{c}"]) if fe_split: - split = self._extract_fe_split_from_apfu(out, phase) - results["fe2"] = split["fe2"] - results["fe3"] = split["fe3"] - - return results, out + split = self._extract_fe_split_from_apfu(out, phase, instance='all') + _store("Fe2", split["Fe2"]) + _store("Fe3", split["Fe3"]) + + return instances, out diff --git a/src/phasetools/core/README.md b/src/phasetools/core/README.md index c87f88e..2aa41f3 100644 --- a/src/phasetools/core/README.md +++ b/src/phasetools/core/README.md @@ -22,3 +22,43 @@ The `core` submodule provides the foundational classes and low-level bridging lo - `get_phase_mg2_number`: Phase-wide $Mg\#$ (Divalent Iron only). - `get_phase_fe_split`: Heuristic splitting of total iron into $\text{Fe}^{2+}$ and $\text{Fe}^{3+}$. - `calculate_kd_fe_mg`: Distribution coefficients between phases. + +### Solvus (multi-instance) handling +MAGEMin reports coexisting solvus limbs as repeated entries in `out.ph` +(e.g. two `dio` clinopyroxenes or two `amp` amphiboles). `phase_frac` +sums them; the low-level composition helpers in `phase_properties.py` +(`get_oxide_apfu`, `get_phase_chemistry`, `extract_end_member`, +`get_phase_mg_number`, `get_phase_mg2_number`, `get_phase_fe_split`) +take an `instance` argument: + +- integer index (default `0`) — that instance, with a `UserWarning` + noting how many other instances exist; out-of-range → zeros; +- `'all'` — one value per instance, returned as numpy arrays. + +`MAGEMinPTGridCalculator.extract_from_grid` / `single_point_calc` / +`generate_2D_grid` return a **list of per-instance bundles**, indexed +`res[0]`, `res[1]`, ... : + +```python +res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) +c0, c1 = res # limb 0, limb 1 +c0['mol_frac'] # array over grid points (limb 0) +c0['ox_apfu_Na2O'] # array over grid points (limb 0) +total_dio = c0['mol_frac'] + c1['mol_frac'] # sum the limbs yourself +``` + +- Bundle keys are unsuffixed (`mol_frac`, `wt_frac`, `vol_frac`, + `ox_apfu_*`, `chem_*`, `em_*`, `cat_*`, `Mg_number`, `Fe2`, `Fe3`). +- The shape is the **same for every phase** — a single-instance phase + (e.g. garnet) has a one-element list, always at `res[0]`. +- Each bundle value is an array over the grid points, **NaN** where the + phase (or that limb) is absent; if the phase never appears, the list + is empty. No totals are computed — sum the `mol_frac` bundles if + you want the whole-phase fraction. +- Bundle `res[k]` is the k-th occurrence of the phase in that point's + `out.ph`; solvus branch ordering may swap between grid points, so + track both limbs when plotting isopleths. `end_members='auto'` + discovers end-members from the first occurrence (solvus limbs share + the same solution model). `MAGEMinGarnetCalculator. + generate_2D_grid_gt_endmembers` returns the single-instance bundle + directly, preserving the historical `em_py`, `mol_frac`, ... keys. diff --git a/src/phasetools/core/base.py b/src/phasetools/core/base.py index 8a13bde..d81bd75 100644 --- a/src/phasetools/core/base.py +++ b/src/phasetools/core/base.py @@ -75,9 +75,9 @@ def setup_bulk_composition(self, Xoxides, X, sys_in, rm_list=None): except: pass - def _extract_fe_split_from_apfu(self, out, phase): + def _extract_fe_split_from_apfu(self, out, phase, instance=0): """ Internal: Calculate Fe2+ and Fe3+ amounts for a phase using an excess oxygen heuristic. """ from .phase_properties import get_phase_fe_split - return get_phase_fe_split(out, phase) + return get_phase_fe_split(out, phase, instance=instance) diff --git a/src/phasetools/core/phase_properties.py b/src/phasetools/core/phase_properties.py index f54efe8..02a0e1f 100644 --- a/src/phasetools/core/phase_properties.py +++ b/src/phasetools/core/phase_properties.py @@ -1,24 +1,93 @@ import numpy as np +import warnings -def get_oxide_apfu(out, ph, oxides): - """Extract oxide amounts from APFU output for a specific phase.""" - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - oxide_values = np.array(phase_obj.Comp_apfu, dtype=float) - oxide_names = [str(ox) for ox in out.oxides] - oxide_moles_dict = {ox: value for ox, value in zip(oxide_names, oxide_values)} +def _phase_indices(out, phase, instance=0): + """Return the index/indices of ``phase`` in ``out.ph``. - results = {} - for oxide in oxides: - results[oxide] = oxide_moles_dict.get(oxide, 0.0) - except (ValueError, IndexError): - results = {oxide: 0.0 for oxide in oxides} + Parameters + ---------- + out : object + MAGEMin output object with a ``ph`` attribute. + phase : str + Phase name. + instance : int or {'all'}, default=0 + Which instance to resolve. An integer index selects that instance + (negative indices count from the end). ``'all'`` returns every + occurrence. When the phase appears more than once (a solvus -- e.g. + two coexisting clinopyroxenes or amphiboles, reported by MAGEMin as + repeated entries in ``out.ph``), a warning notes how many other + instances exist. + + Returns + ------- + list[int] + Indices of ``phase`` in ``out.ph`` (empty list if the phase is + absent or the requested instance does not exist). + """ + idx = [i for i, p in enumerate(out.ph) if str(p) == phase] + if instance == 'all': + return idx + if not isinstance(instance, (int, np.integer)): + raise ValueError(f"instance must be an integer index or 'all', got {instance!r}") + if not idx: + return [] + if instance < 0: + instance = len(idx) + instance + if instance < 0 or instance >= len(idx): + warnings.warn( + f"phase '{phase}' has {len(idx)} instance(s); requested instance " + f"{instance} does not exist -- treating as absent.", + UserWarning, stacklevel=3) + return [] + if len(idx) > 1: + warnings.warn( + f"phase '{phase}' has {len(idx)} instances (solvus) in this " + f"output; using instance {instance} (there are {len(idx) - 1} " + f"others). Use instance='all' for per-instance arrays.", + UserWarning, stacklevel=3) + return [idx[instance]] - return results +def get_oxide_apfu(out, ph, oxides, instance=0): + """Extract oxide amounts from APFU output for a specific phase. -def get_phase_chemistry(out, ph, oxides, sys_in): + Parameters + ---------- + out : object + MAGEMin output object. + ph : str + Phase name. + oxides : list[str] + Oxides to extract. + instance : int or {'all'}, default=0 + For a phase that appears multiple times (a solvus), an integer + index selects that instance (warning if more than one exists) and + ``'all'`` returns one value per instance as numpy arrays. + + Returns + ------- + dict + ``{oxide: value}`` for an integer index, or ``{oxide: ndarray}`` + (one entry per phase instance) for ``'all'``. + """ + idx = _phase_indices(out, ph, instance) + if not idx: + return {oxide: 0.0 for oxide in oxides} if instance != 'all' \ + else {oxide: np.zeros(0) for oxide in oxides} + oxide_names = [str(ox) for ox in out.oxides] + per = [] + for i in idx: + try: + values = np.array(out.SS_vec[i].Comp_apfu, dtype=float) + per.append({ox: float(values[oxide_names.index(ox)]) + if ox in oxide_names else 0.0 for ox in oxides}) + except (ValueError, IndexError, AttributeError): + per.append({ox: 0.0 for ox in oxides}) + if instance == 'all': + return {ox: np.array([d[ox] for d in per]) for ox in oxides} + return per[0] + +def get_phase_chemistry(out, ph, oxides, sys_in, instance=0): """ Extract oxide concentrations (wt% or mol%) for a specific phase. @@ -32,46 +101,61 @@ def get_phase_chemistry(out, ph, oxides, sys_in): List of oxides to extract. sys_in : str Unit system ('wt' or 'mol'). + instance : int or {'all'}, default=0 + For a phase that appears multiple times (a solvus), an integer + index selects that instance (warning if more than one exists) and + ``'all'`` returns one value per instance as numpy arrays. Returns ------- dict - Oxide concentrations. + Oxide concentrations for an integer index, or per-instance arrays + for ``'all'``. """ - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - - if sys_in.casefold() == 'wt': - # Comp_wt is weight fraction (0-1) for oxides in the phase - values = np.array(phase_obj.Comp_wt, dtype=float) * 100.0 - else: - # Comp is molar fraction (0-1) for oxides in the phase - values = np.array(phase_obj.Comp, dtype=float) * 100.0 - - oxide_names = [str(ox) for ox in out.oxides] - oxide_dict = {ox: value for ox, value in zip(oxide_names, values)} - - results = {} - for oxide in oxides: - results[oxide] = oxide_dict.get(oxide, 0.0) - except (ValueError, IndexError): - results = {oxide: 0.0 for oxide in oxides} + idx = _phase_indices(out, ph, instance) + if not idx: + return {oxide: 0.0 for oxide in oxides} if instance != 'all' \ + else {oxide: np.zeros(0) for oxide in oxides} + oxide_names = [str(ox) for ox in out.oxides] + per = [] + for i in idx: + try: + phase_obj = out.SS_vec[i] + if sys_in.casefold() == 'wt': + # Comp_wt is weight fraction (0-1) for oxides in the phase + values = np.array(phase_obj.Comp_wt, dtype=float) * 100.0 + else: + # Comp is molar fraction (0-1) for oxides in the phase + values = np.array(phase_obj.Comp, dtype=float) * 100.0 + per.append({ox: float(values[oxide_names.index(ox)]) + if ox in oxide_names else 0.0 for ox in oxides}) + except (ValueError, IndexError, AttributeError): + per.append({ox: 0.0 for ox in oxides}) + if instance == 'all': + return {ox: np.array([d[ox] for d in per]) for ox in oxides} + return per[0] - return results +def extract_end_member(phase, MAGEMinOutput, end_member, sys_in, instance=0): + """Extract specific end-member fraction from MAGEMin output. -def extract_end_member(phase, MAGEMinOutput, end_member, sys_in): - """Extract specific end-member fraction from MAGEMin output.""" - try: - phase_ind = MAGEMinOutput.ph.index(phase) - em_index = MAGEMinOutput.SS_vec[phase_ind].emNames.index(end_member) - if sys_in.casefold() == 'wt': - data = MAGEMinOutput.SS_vec[phase_ind].emFrac_wt[em_index] - else: - data = MAGEMinOutput.SS_vec[phase_ind].emFrac[em_index] - except (ValueError, IndexError): - data = 0. - return data + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as a numpy array. + """ + idx = _phase_indices(MAGEMinOutput, phase, instance) + if not idx: + return 0.0 if instance != 'all' else np.zeros(0) + vals = [] + for i in idx: + try: + em_index = MAGEMinOutput.SS_vec[i].emNames.index(end_member) + if sys_in.casefold() == 'wt': + vals.append(float(MAGEMinOutput.SS_vec[i].emFrac_wt[em_index])) + else: + vals.append(float(MAGEMinOutput.SS_vec[i].emFrac[em_index])) + except (ValueError, IndexError, AttributeError): + vals.append(0.0) + return float(vals[0]) if instance != 'all' else np.array(vals) def phase_frac(phase, MAGEMinOutput, sys_in): """ @@ -97,7 +181,7 @@ def phase_frac(phase, MAGEMinOutput, sys_in): except: return 0.0 -def get_phase_mg_number(out, ph): +def get_phase_mg_number(out, ph, instance=0): """ Calculate Mg# (molar Mg / (Mg + Fe_total)) for a specific phase. @@ -107,109 +191,97 @@ def get_phase_mg_number(out, ph): Matches the logic used by MAGEMin's 'ss_MgNum' mode by pulling MgO and FeO directly from the phase's Comp_apfu array. Supports 'FeO', 'Fe' (sb24), and 'Fe2O3' fallback. + + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as a numpy array. """ - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - - oxide_names = [str(ox) for ox in out.oxides] - - # Pull Mg + idx = _phase_indices(out, ph, instance) + if not idx: + return 0.0 if instance != 'all' else np.zeros(0) + oxide_names = [str(ox) for ox in out.oxides] + + def _mg(i): try: - mg_idx = oxide_names.index('MgO') - mg = float(phase_obj.Comp_apfu[mg_idx]) - except ValueError: - mg = 0.0 - - # Pull Fe (total iron atoms) - fe = 0.0 - if 'FeO' in oxide_names: - fe_idx = oxide_names.index('FeO') - fe += float(phase_obj.Comp_apfu[fe_idx]) - # If both are present, we sum them (though unlikely in standard MAGEMin output) - if 'Fe2O3' in oxide_names: - fe2o3_idx = oxide_names.index('Fe2O3') - fe += 2.0 * float(phase_obj.Comp_apfu[fe2o3_idx]) - elif 'Fe' in oxide_names: - fe_idx = oxide_names.index('Fe') - fe += float(phase_obj.Comp_apfu[fe_idx]) - elif 'Fe2O3' in oxide_names: - fe2o3_idx = oxide_names.index('Fe2O3') - fe += 2.0 * float(phase_obj.Comp_apfu[fe2o3_idx]) - else: - # No iron components found - if mg == 0: return 0.0 - return 1.0 # Pure Mg phase - - denominator = mg + fe - if denominator == 0: + phase_obj = out.SS_vec[i] + mg = float(phase_obj.Comp_apfu[oxide_names.index('MgO')]) if 'MgO' in oxide_names else 0.0 + fe = 0.0 + if 'FeO' in oxide_names: + fe += float(phase_obj.Comp_apfu[oxide_names.index('FeO')]) + if 'Fe2O3' in oxide_names: + fe += 2.0 * float(phase_obj.Comp_apfu[oxide_names.index('Fe2O3')]) + elif 'Fe' in oxide_names: + fe += float(phase_obj.Comp_apfu[oxide_names.index('Fe')]) + elif 'Fe2O3' in oxide_names: + fe += 2.0 * float(phase_obj.Comp_apfu[oxide_names.index('Fe2O3')]) + else: + return 1.0 if mg > 0 else 0.0 + denom = mg + fe + return mg / denom if denom else 0.0 + except (ValueError, IndexError, AttributeError): return 0.0 - return mg / denominator - except (ValueError, IndexError, AttributeError): - return 0.0 + vals = [_mg(i) for i in idx] + return float(vals[0]) if instance != 'all' else np.array(vals) -def get_phase_fe_split(out, ph): +def get_phase_fe_split(out, ph, instance=0): """ Calculate Fe2+ and Fe3+ amounts for a phase using an excess oxygen heuristic. Works for both traditional FeO-Fe2O3 bases and MAGEMin's O-basis (ig, mp). + + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as numpy arrays. """ - try: - ox_to_query = ['FeO', 'Fe2O3', 'Fe', 'O'] - apfu = get_oxide_apfu(out, ph, ox_to_query) - - feo_val = apfu.get("FeO", 0.0) - fe2o3_val = apfu.get("Fe2O3", 0.0) - fe_metal_val = apfu.get("Fe", 0.0) - atomic_o = apfu.get("O", 0.0) - - # 1. Calculate Total Fe atoms (Atoms per formula unit) - if atomic_o > 0: - # MAGEMin O-basis (ig, mp) or sb24 basis - # If Fe component is present (sb24), use it; otherwise FeO is total iron. - if fe_metal_val > 0: - total_fe = fe_metal_val - else: - total_fe = feo_val - else: - # Traditional FeO/Fe2O3 basis - total_fe = feo_val + 2.0 * fe2o3_val - - # 2. Calculate Fe3+ atoms using excess oxygen heuristic - # excess_o identifies oxygen atoms added beyond the stoichiometric baseline. - # Works for both 'O as total oxygen' and 'O as excess oxygen' components. - excess_o = max(atomic_o - np.round(atomic_o), 0.0) - fe3 = 2.0 * fe2o3_val + 2.0 * excess_o - - # 3. Divalent iron is the remainder - fe2 = max(total_fe - fe3, 0.0) + apfu = get_oxide_apfu(out, ph, ['FeO', 'Fe2O3', 'Fe', 'O'], instance=instance) - return { - "fe2": fe2, - "fe3": fe3, - } - except: - return {"fe2": 0.0, "fe3": 0.0} + feo = np.asarray(apfu.get("FeO", 0.0), dtype=float) + fe2o3 = np.asarray(apfu.get("Fe2O3", 0.0), dtype=float) + fem = np.asarray(apfu.get("Fe", 0.0), dtype=float) + ato = np.asarray(apfu.get("O", 0.0), dtype=float) + + if ato.size == 0: + empty = {"Fe2": np.zeros(0), "Fe3": np.zeros(0)} + return {"Fe2": 0.0, "Fe3": 0.0} if instance != 'all' else empty -def get_phase_mg2_number(out, ph): + # 1. Calculate Total Fe atoms (Atoms per formula unit) + if np.any(ato > 0): + # MAGEMin O-basis (ig, mp) or sb24 basis + total_fe = fem if np.any(fem > 0) else feo + else: + # Traditional FeO/Fe2O3 basis + total_fe = feo + 2.0 * fe2o3 + + # 2. Calculate Fe3+ atoms using excess oxygen heuristic + excess_o = np.maximum(ato - np.floor(ato), 0.0) + fe3 = 2.0 * fe2o3 + 2.0 * excess_o + + # 3. Divalent iron is the remainder + fe2 = np.maximum(np.asarray(total_fe) - fe3, 0.0) + + if instance != 'all': + return {"Fe2": float(np.squeeze(fe2)), "Fe3": float(np.squeeze(fe3))} + return {"Fe2": fe2, "Fe3": fe3} + +def get_phase_mg2_number(out, ph, instance=0): """ Calculate Mg# (molar Mg / (Mg + Fe2+)) for a specific phase. Uses an excess oxygen heuristic to split total iron into Fe2+ and Fe3+. """ try: - apfu = get_oxide_apfu(out, ph, ['MgO']) + apfu = get_oxide_apfu(out, ph, ['MgO'], instance=instance) mg = apfu.get('MgO', 0.0) - split = get_phase_fe_split(out, ph) - fe2 = split['fe2'] + split = get_phase_fe_split(out, ph, instance=instance) + fe2 = split['Fe2'] denominator = mg + fe2 - if denominator == 0: - return 0.0 + if np.any(np.asarray(denominator) <= 0): + return 0.0 if instance != 'all' else np.zeros(len(mg)) - return mg / denominator + return (mg / denominator) if instance != 'all' else np.asarray(mg / denominator, dtype=float) except: return 0.0 diff --git a/src/phasetools/models/README.md b/src/phasetools/models/README.md index f6f2477..65191d7 100644 --- a/src/phasetools/models/README.md +++ b/src/phasetools/models/README.md @@ -7,8 +7,11 @@ Models are complex, multi-step simulations that chain multiple thermodynamic cal ### `garnet_growth.py` - **`GarnetGenerator`**: Simulates the fractional growth of garnet crystals. - Handles radial shell zoning, cohort formation times, and size-frequency distributions. +- When `fractionate=True`, garnet is removed from the reactive bulk at each P-T step. The fractionation amount is based on the appropriate unit basis (mol or wt) matching `sys_in`. +- `get_retrograde_concentrations` recalculates retrograde compositions using the bulk at last growth (fixed-bulk assumption — garnet resorption during retrograde is not modelled). ### `magma_ocean.py` -- **`MagmaOceanModel`**: Simulates the cooling and solidification of a planetary magma ocean. +- **`MagmaOcean`**: Simulates the cooling and solidification of a planetary magma ocean. - Transitions between equilibrium and fractional crystallisation stages. - Integrates planetary-scale geophysical parameters (radius, gravity) to calculate pressure-depth relationships. +- `run_fractional_stages` preserves and restores the instance's bulk composition (`self.X`) after execution, so repeated calls do not mutate state. diff --git a/src/phasetools/models/garnet_growth.py b/src/phasetools/models/garnet_growth.py index 5516cc3..b9dfa52 100644 --- a/src/phasetools/models/garnet_growth.py +++ b/src/phasetools/models/garnet_growth.py @@ -3,6 +3,7 @@ from scipy.interpolate import interp1d from scipy.stats import norm from typing import Any +from juliacall import Main as jl, convert as jlconvert from ..calculators.garnet import MAGEMinGarnetCalculator def generate_distribution(n_classes, r_min, dr, fnr, Gn, tGn): @@ -42,8 +43,8 @@ class GarnetGenerator(MAGEMinGarnetCalculator): Generate synthetic garnet populations with compositional zoning along P-T-t paths. """ - def __init__(self, db="mpe", dataset=636, verbose=False): - super().__init__(db, dataset, verbose) + def __init__(self, db="mpe", dataset=636, verbose=False, fe_basis="FeOt"): + super().__init__(db, dataset, verbose, fe_basis=fe_basis) def setup_bulk_composition(self, Xoxides, X, sys_in, rm_list=None): super().setup_bulk_composition(Xoxides, X, sys_in, rm_list) @@ -79,6 +80,16 @@ def generate_garnet_data(self, normalise_start : bool, default=True If True, the initial garnet volume is set to 0 and only new growth is modeled. If False, the initial thermodynamic volume is used as the starting point. + + Notes + ----- + When both ``fractionate=True`` and ``normalise_start=False``, the + initial garnet fraction is removed from the reactive bulk at the first + P-T point (overstepped nucleation). Subsequent growth increments are + then fractionated from the depleted bulk. When + ``normalise_start=True``, only new growth beyond the initial fraction + is modelled and fractionated; the initial garnet is treated as a + non-reactive seed. """ @@ -288,13 +299,26 @@ def get_prograde_concentrations(self, new_t=None): def get_retrograde_concentrations(self, new_t=None): """Get the retrograde concentrations of garnet-forming elements. - Parameters: - new_t (array-like, optional): New time values to interpolate the data - and return the concentrations at these times. If None, the original - data is returned. Uses a linear interpolation between datapoints. - - Returns: - Concentrations (array): An array with the element concentrations and PTt data at each retrograde step. + Parameters + ---------- + new_t : array-like, optional + New time values to interpolate the data and return the + concentrations at these times. If ``None``, the original data is + returned. Uses linear interpolation between datapoints. + + Returns + ------- + numpy.ndarray + Array with element concentrations and PTt data at each retrograde + step. Rows are ``[t, T, P, Mn, Mg, Fe, Ca]``. + + Notes + ----- + Retrograde concentrations are recalculated using the bulk composition + at the **last growth** step, not the bulk at the end of the P-T path. + This means garnet resorption during retrograde is **not** accounted + for — the retrograde calculation assumes a fixed bulk equal to the + last-growth composition. """ GVi = np.array(self.gt_vol_frac) @@ -334,9 +358,10 @@ def get_retrograde_concentrations(self, new_t=None): Ca_eval = np.zeros_like(t_eval, dtype=float) for i in range(len(t_eval)): + x_jl = jlconvert(jl.Vector[jl.Float64], x_last_growth) (_gt_frac, _gt_wt, _gt_vol, - Mg_i, Mn_i, Fe_i, Ca_i, _out) = self.gt_single_point_calc_elements( - P_eval[i], T_eval[i], self.data, x_last_growth, self.Xoxides, self.sys_in, self.rm_list + Mg_i, Mn_i, Fe_i, Ca_i, _out) = self._gt_single_point_from_jl( + P_eval[i], T_eval[i], x_jl, self.Xoxides, self.sys_in, self.rm_list ) Mn_eval[i] = Mn_i Mg_eval[i] = Mg_i @@ -506,9 +531,9 @@ def plot_garnet_summary(self, size_dist='N', garnet_no=0, path=None, plot_fig=Tr axs[0, 0].set_xlabel('r') axs[0, 0].set_ylabel('f') axs[0, 0].set_xlim([r_r.min(), r_r.max()]) - axs[0, 0].plot(r_r, finp, '-', label='Size Distribution') + axs[0, 0].plot(r, finp, '-', label='Size Distribution') for i in range(n_classes): - axs[0, 0].plot([r_r[i], r_r[i]], [0, finp[i]], '-') + axs[0, 0].plot([r[i], r[i]], [0, finp[i]], '-') axs[0, 0].legend() # Subplot 2: Classes' birth place @@ -588,34 +613,3 @@ def plot_garnet_summary(self, size_dist='N', garnet_no=0, path=None, plot_fig=Tr plt.show() else: plt.close() - GVi = np.array(self.gt_vol_frac) - GVn = self._compute_normalised_GVG(GVi) - first_one_idx = self._first_one_index(GVn) - try: - last_zero_idx = self._last_zero_before(GVn, first_one_idx, strict=True) - except IndexError: - last_zero_idx = -1 - ind = np.arange(last_zero_idx+1, first_one_idx+1) - tG, TG, PG, MnG, MgG, FeG, CaG = self._slice_arrays(ind) - n_classes, r, dr, finp, fnr = self._build_size_distribution(size_dist) - Gn = GVn[ind] / np.max(GVn[ind]) - G, t_arr, r_r, R = generate_distribution(n_classes, self.r_min, dr, fnr, Gn, tG) - PGrw, TGrw, Mnrw, Mgrw, Ferw = [self._interp(tG, arr, t_arr) for arr in [PG, TG, MnG, MgG, FeG]] - Carw = 1 - Mnrw - Mgrw - Ferw - - fig, axs = plt.subplots(3, 2, figsize=(10, 15)) - fig.suptitle('Garnet formation summary') - axs[0, 0].plot(r_r, finp, '-'); [axs[0, 0].plot([r_r[i], r_r[i]], [0, finp[i]], '-') for i in range(n_classes)] - axs[0, 1].plot(self.Ti, self.Pi, 'k-'); axs[0, 1].plot(TGrw, PGrw, 'r.') - for i in range(0, n_classes, 10): axs[1, 0].plot(t_arr[i:], R[i, i:], '.-') - axs[1, 1].plot(tG, Gn, 'k', drawstyle='steps-post'); axs[1, 1].plot(t_arr, G, 'mx') - for i in np.arange(0, n_classes, 10): - idx = np.arange(i, n_classes); rplt = R[i, idx] - axs[2, 0].plot(rplt, Mnrw[idx], '-b'); axs[2, 0].plot(rplt, Mgrw[idx], '-g') - axs[2, 0].plot(rplt, Ferw[idx], '-r'); axs[2, 0].plot(rplt, Carw[idx], '-', c='gold') - i = garnet_no; idx = np.arange(i, n_classes); rplt = R[i, idx] - axs[2, 1].plot(rplt, Mnrw[idx], 'b-'); axs[2, 1].plot(rplt, Mgrw[idx], 'g-') - axs[2, 1].plot(rplt, Ferw[idx], 'r-'); axs[2, 1].plot(rplt, Carw[idx], '-', c='gold') - plt.tight_layout(rect=(0, 0.03, 1, 0.95)) - if path: plt.savefig(path) - plt.show() if plot_fig else plt.close() diff --git a/src/phasetools/models/magma_ocean.py b/src/phasetools/models/magma_ocean.py index 78f67c7..555e7a9 100644 --- a/src/phasetools/models/magma_ocean.py +++ b/src/phasetools/models/magma_ocean.py @@ -1,5 +1,6 @@ import numpy as np import sys +import warnings from scipy import optimize from typing import List, Dict, Any, Tuple, Optional from ..core.base import MAGEMinBase @@ -89,7 +90,12 @@ def func(T): except ValueError: f_low = func(bracket[0]) f_high = func(bracket[1]) - return float(bracket[0] if abs(f_low) < abs(f_high) else bracket[1]) + endpoint = bracket[0] if abs(f_low) < abs(f_high) else bracket[1] + warnings.warn( + f"find_temperature_at_vol_frac: bisection failed for P={P}, " + f"target_vol_frac={target_vol_frac}; returning bracket endpoint T={endpoint}." + ) + return float(endpoint) def get_phase_chemistry_at_index(self, out, i: int) -> np.ndarray: """Extract the chemical composition vector of a phase at a specific index.""" @@ -121,6 +127,7 @@ def run_stage_0(self, p_start: float, p_end: float, solid_frac: float = 0.5, p_i } melt_sum = np.zeros(len(self._Xoxides_py)) + n_melt_samples = 0 layer_modes_sum = {} for P in pressures: @@ -145,15 +152,19 @@ def run_stage_0(self, p_start: float, p_end: float, solid_frac: float = 0.5, p_i if ph_str == 'liq': melt_comp = self.get_phase_chemistry_at_index(out, i) melt_sum += melt_comp + n_melt_samples += 1 else: layer_modes_sum[ph_str] = layer_modes_sum.get(ph_str, 0.0) + vfrac results["modes"].append(modes) results["densities"].append(densities) - avg_melt = melt_sum / p_intervals + avg_melt = melt_sum / max(n_melt_samples, 1) total_solid = sum(layer_modes_sum.values()) - results["layer_modes"] = {ph: val / total_solid for ph, val in layer_modes_sum.items()} + if total_solid > 0: + results["layer_modes"] = {ph: val / total_solid for ph, val in layer_modes_sum.items()} + else: + results["layer_modes"] = {} return results, avg_melt @@ -191,6 +202,12 @@ def run_fractional_stages(self, """ all_stage_results = [] current_melt_comp = starting_melt + + if len(starting_melt) != len(self._Xoxides_py): + raise ValueError( + f"starting_melt length ({len(starting_melt)}) does not match " + f"bulk composition oxides ({len(self._Xoxides_py)})." + ) # Calculate volume of total LMO based on the initial melt ocean bounds v_total_mo_init = self.get_volume_between_radii(self.pressure_to_radius(p_start), self.pressure_to_radius(p_end)) @@ -205,76 +222,80 @@ def run_fractional_stages(self, r_top = self.pressure_to_radius(p_end) current_liquid_vol_frac = starting_vol_frac - - for stage in range(1, n_stages + 1): - # Set composition - self.X = jlconvert(jl.Vector[jl.Float64], current_melt_comp) - - # Base pressure of the current liquid ocean - p_base = self.radius_to_pressure(r_bottom) - - # Per Johnson et al. 2021: Stage 1 concludes when 5 vol% solid is reached - # for the WHOLE melt ocean. Target solid frac = vol_step / current_liquid_vol_frac. - # Clamp to 1.0 to prevent minimization failure in the final stage. - target_solid_frac = min(vol_step / current_liquid_vol_frac, 1.0) - - # Find temperature at base pressure for target solid fraction - T = self.find_temperature_at_vol_frac(p_base, target_solid_frac) - - # Run minimization at the base pressure - out = MAGEMin_C.single_point_minimization(p_base, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list) - - stage_results = { - "stage": stage, - "p_base": p_base, - "p_top": self.radius_to_pressure(r_top), - "T": T, - "modes": {}, - "densities": {}, - "layer_modes": {} - } - - layer_modes_sum = {} - for i, ph_name in enumerate(out.ph): - ph_str = str(ph_name) - vfrac = float(out.ph_frac_vol[i]) - stage_results["modes"][ph_str] = vfrac + + saved_X = self.X + try: + for stage in range(1, n_stages + 1): + # Set composition + self.X = jlconvert(jl.Vector[jl.Float64], current_melt_comp) - if i < out.n_SS: - rho = float(out.SS_vec[i].rho) - else: - rho = float(out.PP_vec[i - out.n_SS].rho) - stage_results["densities"][ph_str] = rho + # Base pressure of the current liquid ocean + p_base = self.radius_to_pressure(r_bottom) - if ph_str == 'liq': - current_melt_comp = self.get_phase_chemistry_at_index(out, i) + # Per Johnson et al. 2021: Stage 1 concludes when 5 vol% solid is reached + # for the WHOLE melt ocean. Target solid frac = vol_step / current_liquid_vol_frac. + # Clamp to 1.0 to prevent minimization failure in the final stage. + target_solid_frac = min(vol_step / current_liquid_vol_frac, 1.0) + + # Find temperature at base pressure for target solid fraction + T = self.find_temperature_at_vol_frac(p_base, target_solid_frac) + + # Run minimization at the base pressure + out = MAGEMin_C.single_point_minimization(p_base, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list) + + stage_results = { + "stage": stage, + "p_base": p_base, + "p_top": self.radius_to_pressure(r_top), + "T": T, + "modes": {}, + "densities": {}, + "layer_modes": {} + } + + layer_modes_sum = {} + for i, ph_name in enumerate(out.ph): + ph_str = str(ph_name) + vfrac = float(out.ph_frac_vol[i]) + stage_results["modes"][ph_str] = vfrac + + if i < out.n_SS: + rho = float(out.SS_vec[i].rho) + else: + rho = float(out.PP_vec[i - out.n_SS].rho) + stage_results["densities"][ph_str] = rho + + if ph_str == 'liq': + current_melt_comp = self.get_phase_chemistry_at_index(out, i) + else: + layer_modes_sum[ph_str] = layer_modes_sum.get(ph_str, 0.0) + vfrac + + # Normalise solid modes for the layer + total_solid = sum(layer_modes_sum.values()) + if total_solid > 0: + stage_results["layer_modes"] = {ph: val / total_solid for ph, val in layer_modes_sum.items()} else: - layer_modes_sum[ph_str] = layer_modes_sum.get(ph_str, 0.0) + vfrac - - # Normalise solid modes for the layer - total_solid = sum(layer_modes_sum.values()) - if total_solid > 0: - stage_results["layer_modes"] = {ph: val / total_solid for ph, val in layer_modes_sum.items()} - else: - stage_results["layer_modes"] = {} - - # Update Geometry: Sinking minerals raise the bottom radius, floating ones lower the top radius. - pl_frac = sum(stage_results["layer_modes"].get(ph, 0.0) for ph in float_phases) - v_float = v_step * pl_frac - v_sink = v_step * (1.0 - pl_frac) - - # Ensure we don't exceed the available ocean volume (safety bound) - v_ocean = self.get_volume_between_radii(r_bottom, r_top) - v_float = min(v_float, v_ocean) - v_sink = min(v_sink, v_ocean - v_float) - - r_bottom = np.power(r_bottom**3 + (3 * v_sink) / (4 * np.pi), 1/3) - # Ensure r_top doesn't go below r_bottom - r_top = np.power(max(r_top**3 - (3 * v_float) / (4 * np.pi), r_bottom**3), 1/3) - - current_liquid_vol_frac -= vol_step - all_stage_results.append(stage_results) - - if current_liquid_vol_frac <= 0: break - + stage_results["layer_modes"] = {} + + # Update Geometry: Sinking minerals raise the bottom radius, floating ones lower the top radius. + pl_frac = sum(stage_results["layer_modes"].get(ph, 0.0) for ph in float_phases) + v_float = v_step * pl_frac + v_sink = v_step * (1.0 - pl_frac) + + # Ensure we don't exceed the available ocean volume (safety bound) + v_ocean = self.get_volume_between_radii(r_bottom, r_top) + v_float = min(v_float, v_ocean) + v_sink = min(v_sink, v_ocean - v_float) + + r_bottom = np.power(r_bottom**3 + (3 * v_sink) / (4 * np.pi), 1/3) + # Ensure r_top doesn't go below r_bottom + r_top = np.power(max(r_top**3 - (3 * v_float) / (4 * np.pi), r_bottom**3), 1/3) + + current_liquid_vol_frac -= vol_step + all_stage_results.append(stage_results) + + if current_liquid_vol_frac <= 0: break + finally: + self.X = saved_X + return all_stage_results diff --git a/src/phasetools/utils/README.md b/src/phasetools/utils/README.md index 94aa1c2..ccd6c91 100644 --- a/src/phasetools/utils/README.md +++ b/src/phasetools/utils/README.md @@ -7,6 +7,8 @@ The `utils` submodule contains helper functions for chemistry, math, and general ### `bulk_rock.py` - Core stoichiometric engine for managing bulk compositions. - Provides molar mass lookups and unit conversion helpers (e.g., mol% to wt% and vice versa). +- `mol_fractions_to_wt_fractions` / `wt_fractions_to_mol_fractions`: non-normalising converters for single components or subsets (e.g., FeO/Fe2O3/FeOt iron redox conversions). +- `split_feot_to_feo_o` / `express_bulk_in_feo_o_basis`: split total iron (FeOt) into the MAGEMin **FeO + O** redox basis at a target Fe³⁺/FeOt fraction, **conserving FeOt** (FeO column = total Fe; O = f·FeOt/2 per 2FeO + O → Fe₂O₃). Useful for redox sweeps and for converting measured (FeOt-only) bulks into MAGEMin format. - Standardises oxide lists and manages cation-oxide mapping. ### `general.py` diff --git a/src/phasetools/utils/bulk_rock.py b/src/phasetools/utils/bulk_rock.py index 752d683..ea799c0 100644 --- a/src/phasetools/utils/bulk_rock.py +++ b/src/phasetools/utils/bulk_rock.py @@ -44,6 +44,25 @@ def convert_mol_percent_to_wt_percent(mol_percents, components, mass_dict): wt_percents = [(mol * mass_dict[comp] / total_mass) * 100 for comp, mol in zip(components, mol_percents)] return wt_percents +def atomic_frac_to_wt_frac(atomic_frac, mass_dict): + """Convert atomic (molar) site fractions to weight-based site fractions. + + Parameters + ---------- + atomic_frac : dict + Mapping of component names to their atomic fractions (summing to 1.0). + mass_dict : dict + Mapping of component names to their atomic/molecular masses. + + Returns + ------- + dict + Weight-based fractions (summing to 1.0). + """ + total_mass = sum(atomic_frac[k] * mass_dict[k] for k in atomic_frac) + return {k: atomic_frac[k] * mass_dict[k] / total_mass for k in atomic_frac} + + def convert_wt_percent_to_mol_percent(wt_percents, components, mass_dict): """Generic conversion from weight (mass) percent to mole percent.""" total_moles = 0 @@ -52,6 +71,68 @@ def convert_wt_percent_to_mol_percent(wt_percents, components, mass_dict): mol_percents = [((wt / mass_dict[comp]) / total_moles) * 100 for comp, wt in zip(components, wt_percents)] return mol_percents +def mol_fractions_to_wt_fractions(mol, components, mass_dict): + """Convert mole fractions to weight fractions (no normalisation). + + Unlike :func:`convert_mol_percent_to_wt_percent`, the output is not + normalised to sum to 100. This means the function can be applied to a + single component (e.g. a single oxide for an iron redox conversion) or + a subset of a composition, as well as a full composition. + + Parameters + ---------- + mol : float or array_like + Mole fraction(s) of each component. + components : list of str + Component names corresponding to each input value. + mass_dict : dict + Mapping of component names to molecular masses. + + Returns + ------- + float or list or numpy.ndarray + Weight fraction(s) of each component. A scalar input returns a + scalar, a list input returns a list, and any other array-like + input returns a numpy array. + """ + if np.isscalar(mol): + return float(mol) * mass_dict[components[0]] + mol_arr = np.asarray(mol, dtype=float) + masses = np.array([mass_dict[comp] for comp in components], dtype=float) + result = mol_arr * masses + return result.tolist() if isinstance(mol, list) else result + +def wt_fractions_to_mol_fractions(wt, components, mass_dict): + """Convert weight fractions to mole fractions (no normalisation). + + Unlike :func:`convert_wt_percent_to_mol_percent`, the output is not + normalised to sum to 100. This means the function can be applied to a + single component (e.g. a single oxide for an iron redox conversion) or + a subset of a composition, as well as a full composition. + + Parameters + ---------- + wt : float or array_like + Weight fraction(s) of each component. + components : list of str + Component names corresponding to each input value. + mass_dict : dict + Mapping of component names to molecular masses. + + Returns + ------- + float or list or numpy.ndarray + Mole fraction(s) of each component. A scalar input returns a + scalar, a list input returns a list, and any other array-like + input returns a numpy array. + """ + if np.isscalar(wt): + return float(wt) / mass_dict[components[0]] + wt_arr = np.asarray(wt, dtype=float) + masses = np.array([mass_dict[comp] for comp in components], dtype=float) + result = wt_arr / masses + return result.tolist() if isinstance(wt, list) else result + def convert_wt_percent_to_moles(wt_percents, components, mass_dict, total_weight): """Convert weight (mass) percentages to moles.""" moles = [] @@ -76,3 +157,103 @@ def convert_moles_to_mol_percent(moles, components): total = sum(moles_dict.values()) return {comp: (moles_dict[comp] / total) * 100 for comp in components} + +def split_feot_to_feo_o(feot_moles, fe3_frac): + """ + Split total iron (FeOt) into the MAGEMin ``FeO + O`` redox pair at a + target Fe3+/FeOt fraction, conserving the total iron budget. + + Parameters + ---------- + feot_moles : float + Total iron in mole units (atoms of Fe, equivalently the amount of + FeO that would carry all iron as Fe2+). + fe3_frac : float + Target Fe3+/FeOt fraction in ``[0, 1]``. 0 = fully reduced (all + Fe2+), 1 = fully oxidised (all Fe3+). + + Returns + ------- + (feo, o) : tuple[float, float] + ``feo`` is the total-iron column (all Fe expressed as FeO, equal to + ``feot_moles``) and ``o`` is the excess oxygen required to oxidise + the target fraction, per ``2FeO + O -> Fe2O3``. + + Notes + ----- + Molar bookkeeping (each Fe2O3 carries 2 Fe atoms and needs 1 O): + Fe3+ atoms = 2 * O => O = fe3_frac * FeOt / 2 + Fe2+ atoms = FeOt - 2 * O (implicitly held by the FeO column) + """ + fe3_frac = float(fe3_frac) + if not 0.0 <= fe3_frac <= 1.0: + raise ValueError(f"fe3_frac must be in [0, 1], got {fe3_frac!r}") + feot_moles = float(feot_moles) + return feot_moles, fe3_frac * feot_moles / 2.0 + +def express_bulk_in_feo_o_basis(X, Xoxides, fe3_frac, feo_oxide="FeO", + fe2o3_oxide="Fe2O3", o_oxide="O"): + """ + Express a bulk composition in the MAGEMin ``FeO + O`` redox basis at a + target Fe3+/FeOt fraction, conserving total iron. + + The ``FeO`` column is set to total iron (all Fe expressed as FeO, + ``FeOt``) and ``O`` is set to the excess oxygen giving the requested + Fe3+/FeOt partition. Any ``Fe2O3`` component is removed (set to 0) and + an existing ``O`` component is overwritten; missing components are + appended. The returned list preserves the input oxide order. + + Parameters + ---------- + X : array-like + Bulk composition values in MOLE units (mol fractions or mol%). Use + weight-based converters first if the input is in wt%. + Xoxides : list of str + Oxide names corresponding to ``X`` (may include 'FeO', 'Fe2O3', 'O' + in any combination). + fe3_frac : float + Target Fe3+/FeOt fraction in ``[0, 1]``. + feo_oxide, fe2o3_oxide, o_oxide : str + Component names used in ``Xoxides``. + + Returns + ------- + (X_new, Xoxides_new) : tuple[list, list] + Composition in the ``FeO + O`` basis (unnormalised -- pass to + ``convertBulk4MAGEMin`` or normalise afterwards). + + Examples + -------- + >>> X = [50.0, 8.0, 0.5] # SiO2, FeO, O (mol%) + >>> ox = ['SiO2', 'FeO', 'O'] + >>> X2, ox2 = express_bulk_in_feo_o_basis(X, ox, fe3_frac=0.1) + >>> ox2 + ['SiO2', 'FeO', 'O'] + >>> X2[2] == 0.1 * X[1] / 2.0 # O = fe3_frac * FeOt / 2 + True + """ + X = [float(v) for v in X] + Xoxides = list(Xoxides) + feo_i = Xoxides.index(feo_oxide) if feo_oxide in Xoxides else None + fe2o3_i = Xoxides.index(fe2o3_oxide) if fe2o3_oxide in Xoxides else None + o_i = Xoxides.index(o_oxide) if o_oxide in Xoxides else None + + feo_mol = X[feo_i] if feo_i is not None else 0.0 + fe2o3_mol = X[fe2o3_i] if fe2o3_i is not None else 0.0 + feot = feo_mol + 2.0 * fe2o3_mol # total Fe atoms (moles) + + feo_total, o_excess = split_feot_to_feo_o(feot, fe3_frac) + + if feo_i is not None: + X[feo_i] = feo_total + else: + X.append(feo_total) + Xoxides.append(feo_oxide) + if fe2o3_i is not None: + X[fe2o3_i] = 0.0 + if o_i is not None: + X[o_i] = o_excess + else: + X.append(o_excess) + Xoxides.append(o_oxide) + return X, Xoxides diff --git a/tests/test_fe_basis.py b/tests/test_fe_basis.py new file mode 100644 index 0000000..1565e5c --- /dev/null +++ b/tests/test_fe_basis.py @@ -0,0 +1,100 @@ +"""Mock-based tests for the garnet ``fe_basis`` option. + +Verifies that ``MAGEMinGarnetCalculator._extract_garnet_elements_from_oxides`` +honours the ``fe_basis`` flag: + +* ``'FeOt'`` (default) -- total Fe (FeO + 2*Fe2O3 APFU) placed on the + divalent X-site. This is the standard community convention for garnet + end-members / X-site fractions. +* ``'Fe2+'`` -- only the stoichiometrically estimated ferrous iron is + placed on the divalent site, excluding Fe3+. + +No live Julia runtime is needed -- ``get_oxide_apfu`` and the Fe2+/Fe3+ +split are mocked. +""" + +import unittest +from unittest.mock import MagicMock, patch + +from phasetools.calculators.garnet import MAGEMinGarnetCalculator + + +class TestGarnetFeBasis(unittest.TestCase): + """Fe-basis flag on the garnet X-site extraction.""" + + def _make_calc(self, fe_basis): + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + calc.fe_basis = fe_basis # canonical lowercase form set by __init__ + return calc + + def test_feot_uses_total_iron(self): + """FeOt puts FeO + 2*Fe2O3 on the X-site and normalises to 1.""" + calc = self._make_calc('feot') + out = MagicMock() + out.ph = ['g', 'q'] + apfu = {'MgO': 0.6, 'MnO': 0.1, 'CaO': 0.8, 'FeO': 1.2, 'Fe2O3': 0.15} + with patch('phasetools.calculators.garnet.get_oxide_apfu', return_value=apfu): + Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol') + fe_total = 1.2 + 2.0 * 0.15 + total = 0.6 + 0.1 + 0.8 + fe_total + self.assertAlmostEqual(Fe, fe_total / total, places=6) + self.assertAlmostEqual(Mg, 0.6 / total, places=6) + self.assertAlmostEqual(Mn, 0.1 / total, places=6) + self.assertAlmostEqual(Ca, 0.8 / total, places=6) + self.assertAlmostEqual(Mg + Mn + Fe + Ca, 1.0, places=6) + + def test_fe2_uses_split(self): + """Fe2+ uses only the ferrous split, excluding Fe3+.""" + calc = self._make_calc('fe2+') + out = MagicMock() + out.ph = ['g', 'q'] + apfu = {'MgO': 0.6, 'MnO': 0.1, 'CaO': 0.8} + split = {'Fe2': 1.35, 'Fe3': 0.15} + with patch('phasetools.calculators.garnet.get_oxide_apfu', return_value=apfu), \ + patch.object(calc, '_extract_fe_split_from_apfu', return_value=split): + Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol') + total = 0.6 + 0.1 + 0.8 + 1.35 + self.assertAlmostEqual(Fe, 1.35 / total, places=6) + self.assertAlmostEqual(Mg + Mn + Fe + Ca, 1.0, places=6) + self.assertAlmostEqual(Mn, 0.1 / total, places=6) + + def test_absent_garnet_returns_zeros(self): + """No garnet in the assemblage -> all-zero X-site fractions.""" + calc = self._make_calc('feot') + out = MagicMock() + out.ph = ['q', 'dio'] + Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol') + self.assertEqual((Mg, Mn, Fe, Ca), (0.0, 0.0, 0.0, 0.0)) + + def test_invalid_basis_raises(self): + """Unsupported fe_basis values must raise ValueError at construction.""" + with patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__', + return_value=None): + with self.assertRaises(ValueError): + MAGEMinGarnetCalculator(db='ig', fe_basis='Fe3') + + def test_default_is_feot(self): + """The default fe_basis is 'FeOt' (community convention).""" + with patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__', + return_value=None): + calc = MAGEMinGarnetCalculator(db='ig') + self.assertEqual(calc.fe_basis, 'feot') + calc = MAGEMinGarnetCalculator(db='ig', fe_basis='Fe2+') + self.assertEqual(calc.fe_basis, 'fe2+') + calc = MAGEMinGarnetCalculator(db='ig', fe_basis='feot') + self.assertEqual(calc.fe_basis, 'feot') + + +class TestGarnetGeneratorFeBasis(unittest.TestCase): + """GarnetGenerator forwards fe_basis to the garnet calculator.""" + + @patch('phasetools.models.garnet_growth.MAGEMinGarnetCalculator.__init__') + def test_forwards_fe_basis(self, mock_init): + from phasetools.models.garnet_growth import GarnetGenerator + GarnetGenerator(db='mpe', fe_basis='Fe2+') + _, kwargs = mock_init.call_args + self.assertEqual(kwargs.get('fe_basis'), 'Fe2+') + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_fractionation.py b/tests/test_fractionation.py new file mode 100644 index 0000000..a67bc8d --- /dev/null +++ b/tests/test_fractionation.py @@ -0,0 +1,363 @@ +"""Mock-based tests for fractionation correctness fixes. + +These tests verify the fixes described in the implementation plan: +- Fix A: get_retrograde_concentrations crash (garnet_growth.py) +- Fix B: mol/wt unit mismatch in gt_along_path (garnet.py) +- Fix C: pure-phase IndexError in fractionate_phase (phase_search.py) +- Fix D: X_along_path normalise to 1 (garnet.py) +- Fix H1: self.X permanent mutation in run_fractional_stages (magma_ocean.py) + +No live Julia runtime is needed — all MAGEMin calls are mocked. +""" + +import unittest +import numpy as np +from unittest.mock import MagicMock, patch, PropertyMock + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _make_mock_out(phases, n_SS, bulk_mol, bulk_wt, ph_frac_mol, ph_frac_wt, comps_mol, comps_wt): + """Build a minimal mock MAGEMin output object.""" + out = MagicMock() + out.ph = phases + out.n_SS = n_SS + out.bulk = np.array(bulk_mol, dtype=float) + out.bulk_wt = np.array(bulk_wt, dtype=float) + out.ph_frac = list(ph_frac_mol) + out.ph_frac_wt = list(ph_frac_wt) + + ss_vec = [] + pp_vec = [] + for idx, ph in enumerate(phases): + obj = MagicMock() + obj.Comp = np.array(comps_mol[idx], dtype=float) + obj.Comp_wt = np.array(comps_wt[idx], dtype=float) + if idx < n_SS: + ss_vec.append(obj) + else: + pp_vec.append(obj) + + out.SS_vec = ss_vec + out.PP_vec = pp_vec + return out + + +# =========================================================================== +# Fix A: get_retrograde_concentrations crash +# =========================================================================== +class TestRetrogradeNoCrash(unittest.TestCase): + """Fix A: get_retrograde_concentrations must not raise TypeError.""" + + @patch('phasetools.models.garnet_growth.MAGEMinGarnetCalculator.__init__', return_value=None) + def test_retrograde_no_crash(self, mock_init): + """Verify the fixed call to _gt_single_point_from_jl works.""" + from phasetools.models.garnet_growth import GarnetGenerator + + gen = GarnetGenerator.__new__(GarnetGenerator) + # Set up minimal state + gen.Pi = np.array([10.0, 11.0, 12.0, 13.0]) + gen.Ti = np.array([800.0, 810.0, 820.0, 830.0]) + gen.ti = np.array([0.0, 1.0, 2.0, 3.0]) + gen.gt_vol_frac = np.array([0.01, 0.02, 0.05, 0.05]) + gen.Mgi = np.array([0.3, 0.3, 0.3, 0.3]) + gen.Mni = np.array([0.05, 0.05, 0.05, 0.05]) + gen.Fei = np.array([0.4, 0.4, 0.4, 0.4]) + gen.Cai = np.array([0.25, 0.25, 0.25, 0.25]) + # X_along_path: last_growth_index=2 has different bulk than [-1] + gen.X_along_path = np.array([ + [50.0, 50.0], + [48.0, 52.0], + [45.0, 55.0], + [40.0, 60.0], + ], dtype=float) + gen.last_growth_index = 2 + gen.Xoxides = ['SiO2', 'Al2O3'] + gen.sys_in = 'mol' + gen.rm_list = None + + # Mock _gt_single_point_from_jl to return controlled output + mock_return = (0.05, 0.05, 0.05, 0.3, 0.05, 0.4, 0.25, MagicMock()) + gen._gt_single_point_from_jl = MagicMock(return_value=mock_return) + + # Should not raise TypeError + result = gen.get_retrograde_concentrations() + self.assertEqual(result.shape[0], 7) # t, T, P, Mn, Mg, Fe, Ca + self.assertTrue(gen._gt_single_point_from_jl.called) + + +# =========================================================================== +# Fix B: wt fraction uses wt basis +# =========================================================================== +class TestWtFractionUsesWtBasis(unittest.TestCase): + """Fix B: when sys_in='wt', fractionate_phase should use gt_wt, not gt_frac.""" + + @patch('phasetools.calculators.garnet.MAGEMin_C') + @patch('phasetools.calculators.garnet.jlconvert') + @patch('phasetools.calculators.phase_search.PhaseFunctions.__init__', return_value=None) + @patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__', return_value=None) + def test_wt_fraction_uses_wt_basis(self, mock_grid_init, mock_pf_init, mock_jlconvert, mock_magemin): + from phasetools.calculators.garnet import MAGEMinGarnetCalculator + + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + calc.db = 'ig' + calc.dataset = 636 + calc.verbose = False + calc.sys_in = 'wt' + calc.X = np.array([50.0, 50.0]) + calc.Xoxides = MagicMock() + calc._Xoxides_py = ['SiO2', 'Al2O3'] + calc.rm_list = None + calc.data = MagicMock() + + # gt_frac (mol) = 0.10, gt_wt = 0.05 — they differ + # Step 0 returns mol=0.10, wt=0.05; step 1 returns mol=0.12, wt=0.07 + mock_return_0 = (0.10, 0.05, 0.08, 0.3, 0.05, 0.4, 0.25, MagicMock()) + mock_return_1 = (0.12, 0.07, 0.10, 0.3, 0.05, 0.4, 0.25, MagicMock()) + calc._gt_single_point_from_jl = MagicMock(side_effect=[mock_return_0, mock_return_1]) + + # Mock PhaseFunctions + mock_pf = MagicMock() + mock_pf.fractionate_phase = MagicMock(return_value=np.array([50.0, 50.0])) + + # Patch PhaseFunctions at its source module (local import in gt_along_path) + with patch('phasetools.calculators.phase_search.PhaseFunctions', return_value=mock_pf): + calc._copy_state_to = MagicMock() + # jlconvert should pass through the array so np.array(X) works + mock_jlconvert.side_effect = lambda t, v: np.array(v, dtype=float) + + P = np.array([10.0, 11.0]) + T = np.array([800.0, 810.0]) + calc.gt_along_path(P, T, fractionate=True, normalise_start=True) + + # Step 0: normalise_start=True, so no fractionation at i=0 + # Step 1: i>0, frac_amount = gt_wt[1] - gt_wt[0] = 0.07 - 0.05 = 0.02 + calls = mock_pf.fractionate_phase.call_args_list + self.assertEqual(len(calls), 1, f"Expected 1 fractionation call, got {len(calls)}") + _, kwargs = calls[0] + self.assertAlmostEqual(kwargs['frac_amount'], 0.02, places=10, + msg="frac_amount should be based on wt fraction (0.07-0.05), not mol (0.12-0.10)") + + +# =========================================================================== +# Fix D: X_along_path normalised to 1 +# =========================================================================== +class TestXAlongPathNormalised(unittest.TestCase): + """Fix D: each row of X_along_path must sum to 1.0.""" + + @patch('phasetools.calculators.garnet.MAGEMin_C') + @patch('phasetools.calculators.garnet.jlconvert') + @patch('phasetools.calculators.phase_search.PhaseFunctions.__init__', return_value=None) + @patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__', return_value=None) + def test_x_along_path_normalised_to_one(self, mock_grid_init, mock_pf_init, mock_jlconvert, mock_magemin): + from phasetools.calculators.garnet import MAGEMinGarnetCalculator + + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + calc.db = 'ig' + calc.dataset = 636 + calc.verbose = False + calc.sys_in = 'mol' + calc.X = np.array([50.0, 50.0]) # Julia vector — jlconvert will wrap + calc.Xoxides = MagicMock() + calc._Xoxides_py = ['SiO2', 'Al2O3'] + calc.rm_list = None + calc.data = MagicMock() + + mock_return = (0.05, 0.05, 0.05, 0.3, 0.05, 0.4, 0.25, MagicMock()) + calc._gt_single_point_from_jl = MagicMock(return_value=mock_return) + mock_jlconvert.return_value = calc.X + + P = np.array([10.0, 11.0]) + T = np.array([800.0, 810.0]) + + _, _, _, _, _, _, _, X_along_path = calc.gt_along_path(P, T, fractionate=False) + + for i in range(len(P)): + self.assertAlmostEqual(np.sum(X_along_path[i]), 1.0, places=10, + msg=f"Row {i} of X_along_path does not sum to 1.0") + + +# =========================================================================== +# Fix C: pure-phase IndexError +# =========================================================================== +class TestFractionatePurePhase(unittest.TestCase): + """Fix C: fractionate_phase must handle pure phases (PP_vec) without IndexError.""" + + def test_fractionate_pure_phase(self): + from phasetools.calculators.phase_search import PhaseFunctions + + pf = PhaseFunctions.__new__(PhaseFunctions) + + # ph=['q', 'liq'], n_SS=1 => 'q' is a pure phase at index 0 in PP_vec + out = _make_mock_out( + phases=['q', 'liq'], + n_SS=1, + bulk_mol=[60.0, 40.0], + bulk_wt=[62.0, 38.0], + ph_frac_mol=[0.15, 0.85], + ph_frac_wt=[0.16, 0.84], + comps_mol=[[100.0, 0.0], [50.0, 50.0]], + comps_wt=[[100.0, 0.0], [48.0, 52.0]], + ) + + # 'q' is at index 0 in out.ph, n_SS=1, so it's a pure phase (PP_vec[0]) + result = pf.fractionate_phase('q', out, 'mol', frac_amount=0.1) + self.assertIsNotNone(result) + self.assertTrue(np.all(np.isfinite(result))) + + def test_fractionate_solution_phase_still_works(self): + """Verify solution-phase path (SS_vec) still works after the fix.""" + from phasetools.calculators.phase_search import PhaseFunctions + + pf = PhaseFunctions.__new__(PhaseFunctions) + + out = _make_mock_out( + phases=['liq', 'g'], + n_SS=2, + bulk_mol=[60.0, 40.0], + bulk_wt=[62.0, 38.0], + ph_frac_mol=[0.85, 0.15], + ph_frac_wt=[0.84, 0.16], + comps_mol=[[50.0, 50.0], [40.0, 60.0]], + comps_wt=[[48.0, 52.0], [38.0, 62.0]], + ) + + result = pf.fractionate_phase('g', out, 'mol', frac_amount=0.1) + self.assertIsNotNone(result) + self.assertTrue(np.all(np.isfinite(result))) + + +# =========================================================================== +# Zero guard: frac_amount=0 is a no-op +# =========================================================================== +class TestFractionateZeroIsNoop(unittest.TestCase): + """frac_amount=0 must return the bulk unchanged (not normalised to sum=1).""" + + def test_fractionate_zero_is_noop(self): + from phasetools.calculators.phase_search import PhaseFunctions + + pf = PhaseFunctions.__new__(PhaseFunctions) + + out = _make_mock_out( + phases=['g', 'liq'], + n_SS=2, + bulk_mol=[6000.0, 4000.0], # sum=10000, not 1 + bulk_wt=[6200.0, 3800.0], + ph_frac_mol=[0.15, 0.85], + ph_frac_wt=[0.16, 0.84], + comps_mol=[[40.0, 60.0], [50.0, 50.0]], + comps_wt=[[38.0, 62.0], [48.0, 52.0]], + ) + + result = pf.fractionate_phase('g', out, 'mol', frac_amount=0.0) + np.testing.assert_array_equal(result, np.array(out.bulk, dtype=float)) + + +# =========================================================================== +# Fix H1: self.X restored after run_fractional_stages +# =========================================================================== +class TestMagmaOceanXRestored(unittest.TestCase): + """Fix H1: run_fractional_stages must restore self.X after execution.""" + + @patch('phasetools.models.magma_ocean.MAGEMinBase.__init__', return_value=None) + def test_magma_ocean_x_restored(self, mock_base_init): + from phasetools.models.magma_ocean import MagmaOcean + import phasetools.models.magma_ocean as mo_module + + mo = MagmaOcean.__new__(MagmaOcean) + mo._Xoxides_py = ['SiO2', 'Al2O3'] + mo.sys_in = 'mol' + mo.data = MagicMock() + mo.Xoxides = MagicMock() + mo.rm_list = None + mo.X = np.array([50.0, 50.0]) + mo.radius_body = 1737.1 + mo.radius_core = 330.0 + mo.g = 1.62 + mo.rho_avg = 3350.0 + + saved_X = mo.X.copy() + + # Mock find_temperature_at_vol_frac to return a fixed T + mo.find_temperature_at_vol_frac = MagicMock(return_value=1200.0) + + # Build mock MAGEMin output + mock_out = MagicMock() + mock_out.ph = ['ol', 'liq'] + mock_out.n_SS = 2 + mock_out.ph_frac_vol = [0.3, 0.7] + + mock_ol = MagicMock() + mock_ol.rho = 3300.0 + mock_ol.Comp = np.array([30.0, 10.0]) + mock_ol.Comp_wt = np.array([28.0, 12.0]) + + mock_liq = MagicMock() + mock_liq.rho = 2800.0 + mock_liq.Comp = np.array([45.0, 55.0]) + mock_liq.Comp_wt = np.array([43.0, 57.0]) + + mock_out.SS_vec = [mock_ol, mock_liq] + mock_out.PP_vec = [] + + mock_magemin_c = MagicMock() + mock_magemin_c.single_point_minimization = MagicMock(return_value=mock_out) + + # Patch MAGEMin_C and jlconvert at the module level + with patch.object(mo_module, 'MAGEMin_C', mock_magemin_c), \ + patch.object(mo_module, 'jlconvert', side_effect=lambda t, v: np.array(v, dtype=float)): + mo.get_phase_chemistry_at_index = MagicMock(return_value=np.array([45.0, 55.0])) + mo.get_volume_between_radii = MagicMock(return_value=1e12) + mo.pressure_to_radius = MagicMock(return_value=1400.0) + mo.radius_to_pressure = MagicMock(return_value=5.0) + + starting_melt = np.array([45.0, 55.0]) + mo.run_fractional_stages( + starting_melt=starting_melt, + p_start=5.0, + p_end=0.001, + vol_step=0.05, + starting_vol_frac=0.5, + n_stages=2, + ) + + np.testing.assert_array_equal( + mo.X, saved_X, + err_msg="self.X was not restored after run_fractional_stages" + ) + + +# =========================================================================== +# Fix H4: starting_melt length validation +# =========================================================================== +class TestStartingMeltValidation(unittest.TestCase): + """Fix H4: run_fractional_stages must reject mismatched starting_melt length.""" + + @patch('phasetools.models.magma_ocean.MAGEMinBase.__init__', return_value=None) + def test_starting_melt_length_mismatch(self, mock_base_init): + from phasetools.models.magma_ocean import MagmaOcean + + mo = MagmaOcean.__new__(MagmaOcean) + mo._Xoxides_py = ['SiO2', 'Al2O3', 'MgO'] + mo.sys_in = 'mol' + mo.data = MagicMock() + mo.X = np.array([33.0, 33.0, 34.0]) + mo.radius_body = 1737.1 + mo.radius_core = 330.0 + mo.g = 1.62 + mo.rho_avg = 3350.0 + + # starting_melt has 2 elements, but _Xoxides_py has 3 + with self.assertRaises(ValueError) as ctx: + mo.run_fractional_stages( + starting_melt=np.array([50.0, 50.0]), + p_start=5.0, + p_end=0.001, + ) + self.assertIn("does not match", str(ctx.exception)) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_iron_oxide_conversions.py b/tests/test_iron_oxide_conversions.py new file mode 100644 index 0000000..dec9399 --- /dev/null +++ b/tests/test_iron_oxide_conversions.py @@ -0,0 +1,176 @@ +import unittest +import numpy as np + +from phasetools.utils.bulk_rock import ( + get_molar_mass_dict, + get_atomic_mass_dict, + mol_fractions_to_wt_fractions, + wt_fractions_to_mol_fractions, +) + +class TestFractionConverters(unittest.TestCase): + """Tests for the non-normalising mole/weight fraction converters.""" + + def setUp(self): + self.mass_dict = get_molar_mass_dict() + + # ------------------------------------------------------------------ # + # Scalar behaviour + # ------------------------------------------------------------------ # + def test_scalar_mol_to_wt(self): + """A scalar input returns a scalar weight fraction (no normalisation).""" + result = mol_fractions_to_wt_fractions(0.5, ['FeO'], self.mass_dict) + self.assertIsInstance(result, float) + # Not normalised: must be mol * mass, not 100 mol% + self.assertAlmostEqual(result, 0.5 * self.mass_dict['FeO']) + + def test_scalar_wt_to_mol(self): + """A scalar input returns a scalar mole fraction (no normalisation).""" + result = wt_fractions_to_mol_fractions(71.844, ['FeO'], self.mass_dict) + self.assertIsInstance(result, float) + self.assertAlmostEqual(result, 71.844 / self.mass_dict['FeO']) + + # ------------------------------------------------------------------ # + # List behaviour + # ------------------------------------------------------------------ # + def test_list_mol_to_wt(self): + """A list input returns a list (no normalisation to sum to 100).""" + result = mol_fractions_to_wt_fractions([0.4, 0.6], ['FeO', 'Fe2O3'], self.mass_dict) + self.assertIsInstance(result, list) + self.assertAlmostEqual(result[0], 0.4 * self.mass_dict['FeO']) + self.assertAlmostEqual(result[1], 0.6 * self.mass_dict['Fe2O3']) + + def test_list_wt_to_mol(self): + """A list input returns a list (no normalisation to sum to 100).""" + result = wt_fractions_to_mol_fractions( + [self.mass_dict['FeO'], self.mass_dict['Fe2O3']], ['FeO', 'Fe2O3'], self.mass_dict + ) + self.assertIsInstance(result, list) + self.assertAlmostEqual(result[0], 1.0) + self.assertAlmostEqual(result[1], 1.0) + + def test_list_output_not_normalised(self): + """Two oxides given as fractions must NOT be normalised to 100 mol%.""" + wt = mol_fractions_to_wt_fractions([0.4, 0.6], ['FeO', 'Fe2O3'], self.mass_dict) + back = wt_fractions_to_mol_fractions(wt, ['FeO', 'Fe2O3'], self.mass_dict) + self.assertAlmostEqual(back[0], 0.4) + self.assertAlmostEqual(back[1], 0.6) + + # ------------------------------------------------------------------ # + # numpy array behaviour + # ------------------------------------------------------------------ # + def test_array_input_returns_array(self): + """A numpy array input returns a numpy array.""" + result = mol_fractions_to_wt_fractions( + np.array([0.4, 0.6]), ['FeO', 'Fe2O3'], self.mass_dict + ) + self.assertIsInstance(result, np.ndarray) + np.testing.assert_allclose( + result, + np.array([0.4, 0.6]) * np.array([self.mass_dict['FeO'], self.mass_dict['Fe2O3']]), + ) + + def test_single_component_broadcast(self): + """A single component can be applied to multiple values.""" + result = mol_fractions_to_wt_fractions([0.1, 0.2, 0.3], ['FeO'], self.mass_dict) + self.assertEqual(len(result), 3) + self.assertAlmostEqual(result[2], 0.3 * self.mass_dict['FeO']) + + # ------------------------------------------------------------------ # + # Round trips + # ------------------------------------------------------------------ # + def test_round_trip_scalar(self): + """wt -> mol -> wt recovers the input.""" + mol = 0.35 + wt = mol_fractions_to_wt_fractions(mol, ['FeO'], self.mass_dict) + back = wt_fractions_to_mol_fractions(wt, ['FeO'], self.mass_dict) + self.assertAlmostEqual(back, mol) + + def test_round_trip_multi_oxide(self): + """mol -> wt -> mol recovers the input for a pair of oxides.""" + mol = [0.7, 0.3] + wt = mol_fractions_to_wt_fractions(mol, ['FeO', 'Fe2O3'], self.mass_dict) + back = wt_fractions_to_mol_fractions(wt, ['FeO', 'Fe2O3'], self.mass_dict) + self.assertAlmostEqual(back[0], mol[0]) + self.assertAlmostEqual(back[1], mol[1]) + + # ------------------------------------------------------------------ # + # Consistency with the normalising converters + # ------------------------------------------------------------------ # + def test_consistency_with_percent_converters(self): + """The normalised converter is the non-normalised one scaled by a constant.""" + from phasetools.utils.bulk_rock import convert_mol_percent_to_wt_percent + mol_pct = [70.0, 30.0] # sums to 100 + norm = convert_mol_percent_to_wt_percent(mol_pct, ['FeO', 'Fe2O3'], self.mass_dict) + non_norm = mol_fractions_to_wt_fractions(mol_pct, ['FeO', 'Fe2O3'], self.mass_dict) + # Component ratios are identical; normalised output is a scaled version + self.assertAlmostEqual(norm[0] / norm[1], non_norm[0] / non_norm[1]) + scale = sum(norm) / sum(non_norm) + self.assertAlmostEqual(norm[0], non_norm[0] * scale) + self.assertAlmostEqual(norm[1], non_norm[1] * scale) + + # ------------------------------------------------------------------ # + # Iron oxide workflows built on the converters + # ------------------------------------------------------------------ # + def test_feot_from_wt_fractions(self): + """FeOt (wt) = FeO + Fe2O3 * (2 * M_FeO / M_Fe2O3).""" + feo, fe2o3 = 5.0, 2.0 + factor = 2 * self.mass_dict['FeO'] / self.mass_dict['Fe2O3'] + self.assertAlmostEqual(feo + fe2o3 * factor, feo + fe2o3 * factor) + + def test_feot_from_mol_fractions(self): + """FeOt (mol) = FeO + 2 * Fe2O3 (molar basis).""" + feo_mol, fe2o3_mol = 0.05, 0.01 + self.assertAlmostEqual(feo_mol + 2 * fe2o3_mol, 0.07) + + def test_garnet_weight_fractions_sum_to_one(self): + """atomic_frac_to_wt_frac keeps weight fractions summing to 1.0.""" + from phasetools.utils.bulk_rock import atomic_frac_to_wt_frac + atomic = {'Mg': 0.5, 'Mn': 0.1, 'Fe': 0.3, 'Ca': 0.1} + wt = atomic_frac_to_wt_frac(atomic, get_atomic_mass_dict()) + self.assertAlmostEqual(sum(wt.values()), 1.0) + + # ------------------------------------------------------------------ # + # FeOt -> FeO + O (MAGEMin redox basis) split + # ------------------------------------------------------------------ # + def test_split_feot_to_feo_o_conserves_feot(self): + """The FeO column always equals FeOt; only O changes.""" + from phasetools.utils.bulk_rock import split_feot_to_feo_o + feot = 8.33 + for f in (0.0, 0.05, 0.209, 0.5, 1.0): + feo, o = split_feot_to_feo_o(feot, f) + self.assertAlmostEqual(feo, feot) # FeOt conserved + self.assertAlmostEqual(o, f * feot / 2.0) # 2FeO + O -> Fe2O3 + + def test_split_feot_to_feo_o_invalid(self): + """fe3_frac outside [0, 1] must raise ValueError.""" + from phasetools.utils.bulk_rock import split_feot_to_feo_o + with self.assertRaises(ValueError): + split_feot_to_feo_o(8.33, 1.2) + with self.assertRaises(ValueError): + split_feot_to_feo_o(8.33, -0.1) + + def test_express_bulk_in_feo_o_basis(self): + """Bulk split conserves FeOt and sets O to f*FeOt/2 in place.""" + from phasetools.utils.bulk_rock import express_bulk_in_feo_o_basis + X = [50.0, 8.33, 0.0, 0.05] # SiO2, FeO, Fe2O3, H2O + ox = ['SiO2', 'FeO', 'Fe2O3', 'H2O'] + X2, ox2 = express_bulk_in_feo_o_basis(X, ox, fe3_frac=0.2) + # FeOt = FeO + 2*Fe2O3 = 8.33 + self.assertAlmostEqual(X2[1], 8.33) # FeO = FeOt + self.assertAlmostEqual(X2[2], 0.0) # Fe2O3 removed + self.assertAlmostEqual(X2[4], 0.2 * 8.33 / 2.0) # O appended + self.assertEqual(ox2, ['SiO2', 'FeO', 'Fe2O3', 'H2O', 'O']) + + def test_express_bulk_in_feo_o_basis_overwrites_o(self): + """An existing O component is overwritten by the split.""" + from phasetools.utils.bulk_rock import express_bulk_in_feo_o_basis + X = [50.0, 8.33, 0.87] + ox = ['SiO2', 'FeO', 'O'] + X2, ox2 = express_bulk_in_feo_o_basis(X, ox, fe3_frac=0.1) + self.assertAlmostEqual(X2[1], 8.33) + self.assertAlmostEqual(X2[2], 0.1 * 8.33 / 2.0) # O overwritten + self.assertEqual(ox2, ox) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_redox_logic.py b/tests/test_redox_logic.py index 07fb1d7..0d924ac 100644 --- a/tests/test_redox_logic.py +++ b/tests/test_redox_logic.py @@ -33,8 +33,8 @@ def test_garnet_redox_split(self, mock_get_apfu): # excess_o = 12.023 - 12.0 = 0.023 # fe3 = 2 * 0.023 = 0.046 # fe2 = 1.715 - 0.046 = 1.669 - self.assertAlmostEqual(result['fe3'], 0.046, places=3) - self.assertAlmostEqual(result['fe2'], 1.669, places=3) + self.assertAlmostEqual(result['Fe3'], 0.046, places=3) + self.assertAlmostEqual(result['Fe2'], 1.669, places=3) @patch('phasetools.core.phase_properties.get_oxide_apfu') def test_pyroxene_redox_split(self, mock_get_apfu): @@ -52,8 +52,8 @@ def test_pyroxene_redox_split(self, mock_get_apfu): # excess_o = 6.05 - 6.0 = 0.05 # fe3 = 0.1 # fe2 = 0.9 - self.assertAlmostEqual(result['fe3'], 0.1, places=2) - self.assertAlmostEqual(result['fe2'], 0.9, places=2) + self.assertAlmostEqual(result['Fe3'], 0.1, places=2) + self.assertAlmostEqual(result['Fe2'], 0.9, places=2) @patch('phasetools.core.phase_properties.get_oxide_apfu') def test_spinel_standard_split(self, mock_get_apfu): @@ -71,8 +71,8 @@ def test_spinel_standard_split(self, mock_get_apfu): # total_fe = 0.8 + 2*0.1 = 1.0 # fe3 = 2 * 0.1 = 0.2 # fe2 = 0.8 - self.assertAlmostEqual(result['fe3'], 0.2, places=2) - self.assertAlmostEqual(result['fe2'], 0.8, places=2) + self.assertAlmostEqual(result['Fe3'], 0.2, places=2) + self.assertAlmostEqual(result['Fe2'], 0.8, places=2) @patch('phasetools.core.phase_properties.get_oxide_apfu') def test_negative_clamp(self, mock_get_apfu): @@ -86,8 +86,8 @@ def test_negative_clamp(self, mock_get_apfu): result = self.base._extract_fe_split_from_apfu(None, 'dio') - self.assertEqual(result['fe2'], 0.0) - self.assertAlmostEqual(result['fe3'], 0.1, places=2) + self.assertEqual(result['Fe2'], 0.0) + self.assertAlmostEqual(result['Fe3'], 0.1, places=2) @patch('phasetools.core.phase_properties.get_oxide_apfu') def test_sb24_iron_handling(self, mock_get_apfu): @@ -104,8 +104,61 @@ def test_sb24_iron_handling(self, mock_get_apfu): # total_fe should be 1.0 (from 'Fe') # fe3 should be 0.1 (from 2 * 0.05) - self.assertAlmostEqual(result['fe2'], 0.9, places=2) - self.assertAlmostEqual(result['fe3'], 0.1, places=2) + self.assertAlmostEqual(result['Fe2'], 0.9, places=2) + self.assertAlmostEqual(result['Fe3'], 0.1, places=2) + + @patch('phasetools.core.phase_properties.get_oxide_apfu') + def test_bankers_rounding_half_integer(self, mock_get_apfu): + """O APFU at x.5 boundary — banker's rounding gives correct result.""" + mock_get_apfu.return_value = {'FeO': 1.0, 'O': 12.5, 'Fe2O3': 0.0} + result = self.base._extract_fe_split_from_apfu(None, 'g') + # int(12.5) = 12, excess_o = 0.5, fe3 = 1.0 + self.assertAlmostEqual(result['Fe3'], 1.0, places=4) + self.assertAlmostEqual(result['Fe2'], 0.0, places=4) + + @patch('phasetools.core.phase_properties.get_oxide_apfu') + def test_bankers_rounding_above_half(self, mock_get_apfu): + """O APFU above x.5 — int() truncates, giving correct excess.""" + mock_get_apfu.return_value = {'FeO': 1.0, 'O': 13.5, 'Fe2O3': 0.0} + result = self.base._extract_fe_split_from_apfu(None, 'g') + # int(13.5) = 13, excess_o = 0.5, fe3 = 1.0 + self.assertAlmostEqual(result['Fe3'], 1.0, places=4) + self.assertAlmostEqual(result['Fe2'], 0.0, places=4) + + @patch('phasetools.core.phase_properties.get_oxide_apfu') + def test_mixed_fe2o3_and_o_basis(self, mock_get_apfu): + """Both Fe2O3 and O present — no double-counting.""" + mock_get_apfu.return_value = {'FeO': 0.5, 'O': 6.1, 'Fe2O3': 0.1} + result = self.base._extract_fe_split_from_apfu(None, 'dio') + # excess_o = 0.1, fe3 = 2*0.1 + 2*0.1 = 0.4, total_fe = 0.5, fe2 = 0.1 + self.assertAlmostEqual(result['Fe3'], 0.4, places=4) + self.assertAlmostEqual(result['Fe2'], 0.1, places=4) + + @patch('phasetools.core.phase_properties.get_oxide_apfu') + def test_no_iron_phase(self, mock_get_apfu): + """Phase with no iron — returns zero split.""" + mock_get_apfu.return_value = {'FeO': 0.0, 'O': 12.0, 'Fe2O3': 0.0} + result = self.base._extract_fe_split_from_apfu(None, 'q') + self.assertEqual(result['Fe2'], 0.0) + self.assertEqual(result['Fe3'], 0.0) + + @patch('phasetools.core.phase_properties.get_oxide_apfu') + def test_fe3_exceeds_total_fe_clamp(self, mock_get_apfu): + """excess_o suggests more Fe3+ than total Fe — Fe2+ clamped to 0.""" + mock_get_apfu.return_value = {'FeO': 0.1, 'O': 6.3, 'Fe2O3': 0.0} + result = self.base._extract_fe_split_from_apfu(None, 'dio') + # excess_o = 0.3, fe3 = 0.6, total_fe = 0.1, fe2 = max(0.1-0.6, 0) = 0 + self.assertAlmostEqual(result['Fe3'], 0.6, places=4) + self.assertEqual(result['Fe2'], 0.0) + + @patch('phasetools.core.phase_properties.get_oxide_apfu') + def test_zero_o_basis(self, mock_get_apfu): + """Traditional FeO/Fe2O3 basis when O = 0.""" + mock_get_apfu.return_value = {'FeO': 0.8, 'O': 0.0, 'Fe2O3': 0.1} + result = self.base._extract_fe_split_from_apfu(None, 'sp') + # total_fe = 0.8 + 2*0.1 = 1.0, fe3 = 2*0.1 = 0.2, fe2 = 0.8 + self.assertAlmostEqual(result['Fe3'], 0.2, places=4) + self.assertAlmostEqual(result['Fe2'], 0.8, places=4) def test_get_phase_mg_number_robustness(self): """Test that get_phase_mg_number handles Fe and Fe2O3 components.""" diff --git a/tests/test_site_occupancy.py b/tests/test_site_occupancy.py index 400fc7a..ed8a4e4 100644 --- a/tests/test_site_occupancy.py +++ b/tests/test_site_occupancy.py @@ -1,8 +1,15 @@ import unittest import numpy as np from phasetools.calculators.pt_grid import MAGEMinPTGridCalculator -from phasetools.core.phase_properties import get_phase_mg_number, get_phase_mg2_number +from phasetools.core.phase_properties import get_phase_mg_number, get_phase_mg2_number, phase_frac +try: + from juliacall import Main as jl + HAS_JULIA = True +except ImportError: + HAS_JULIA = False + +@unittest.skipUnless(HAS_JULIA, "Julia + MAGEMin_C not available") class TestSiteOccupancy(unittest.TestCase): @classmethod def setUpClass(cls): @@ -17,121 +24,230 @@ def setUpClass(cls): # 30 kbar (3 GPa), 700 C - where g and dio are stable in mpe cls.out = cls.calc.calculate_grid(30.0, 700.0)[0] - def test_garnet_comparison(self): - """ - Compare Garnet Fe2+ and Mg# (X-site) from heuristic vs site occupancy. - - In the Garnet model (White et al. 2014), divalent cations (Mg, Fe2+, Ca, Mn) - are restricted to the 8-fold X-site. However, the Bulk Mg# (Total Iron) - and Site Mg# (Divalent) will still differ if Fe3+ is present. - - This is because Fe3+ sits on the octahedral (Y) site (e.g., in the - khoharite end-member, Mg3Fe2Si3O12). - - - Bulk Mg# (Total): Mg_total / (Mg_total + Fe2_total + Fe3_total) - - Site Mg# (X-site): Mg_X / (Mg_X + Fe2_X) + PT_POINTS = [ + (5.0, 450.0), + (10.0, 550.0), + (20.0, 650.0), + (30.0, 700.0), + (15.0, 800.0), + ] + + def test_garnet_heuristic_across_pt(self): + """Heuristic Fe2+/Fe3+ and Mg# match endmember stoichiometry across P-T range.""" + for P, T in self.PT_POINTS: + out = self.calc.calculate_grid(P, T)[0] + if 'g' not in out.ph: + continue + + split = self.calc._extract_fe_split_from_apfu(out, 'g') + fe2_h, fe3_h = split['Fe2'], split['Fe3'] + mg2_h = get_phase_mg2_number(out, 'g') + + ph_idx = out.ph.index('g') + em = {str(n): float(f) for n, f in zip(out.SS_vec[ph_idx].emNames, out.SS_vec[ph_idx].emFrac)} + fe2_s = 3.0 * em.get('alm', 0.0) + fe3_s = 2.0 * em.get('kho', 0.0) + mg_full = 3.0 * em.get('py', 0.0) + 3.0 * em.get('kho', 0.0) + mg_num_s = mg_full / (mg_full + fe2_s) if (mg_full + fe2_s) > 0 else 0.0 + + with self.subTest(P=P, T=T): + self.assertAlmostEqual(fe2_h, fe2_s, places=2) + self.assertAlmostEqual(fe3_h, fe3_s, places=2) + self.assertAlmostEqual(mg2_h, mg_num_s, places=2) + + def test_cpx_heuristic_across_pt(self): + """Heuristic Fe2+/Fe3+ and Mg# match endmember stoichiometry for cpx across P-T range.""" + for P, T in self.PT_POINTS: + out = self.calc.calculate_grid(P, T)[0] + ph_name = 'dio' if 'dio' in out.ph else ('omph' if 'omph' in out.ph else None) + if ph_name is None: + continue + + split = self.calc._extract_fe_split_from_apfu(out, ph_name) + fe2_h, fe3_h = split['Fe2'], split['Fe3'] + mg2_h = get_phase_mg2_number(out, ph_name) + + ph_idx = out.ph.index(ph_name) + em = {str(n): float(f) for n, f in zip(out.SS_vec[ph_idx].emNames, out.SS_vec[ph_idx].emFrac)} + fe2_s = 1.0 * em.get('hed', 0.0) + 0.5 * em.get('cfm', 0.0) + fe3_s = 1.0 * em.get('acmm', 0.0) + 0.5 * em.get('jac', 0.0) + mg_full = 1.0 * em.get('di', 0.0) + 0.5 * em.get('om', 0.0) + 0.5 * em.get('cfm', 0.0) + mg_num_s = mg_full / (mg_full + fe2_s) if (mg_full + fe2_s) > 0 else 0.0 + + with self.subTest(P=P, T=T): + self.assertAlmostEqual(fe2_h, fe2_s, places=2) + self.assertAlmostEqual(fe3_h, fe3_s, places=2) + self.assertAlmostEqual(mg2_h, mg_num_s, places=2) + + def test_epidote_heuristic_across_pt(self): + """Heuristic Fe3+ matches expected values for epidote across P-T range. - Because khoharite adds Magnesium to the X-site but Ferric iron to the - Y-site, the Bulk Mg# (Total Iron) will be lower than the Site Mg#. + Epidote carries Fe³⁺ (not Fe²⁺). The excess-O heuristic should give + Fe²⁺ ≈ 0 and Fe³⁺ ≈ total Fe at all P-T points where ep is stable. """ - print("\n>>> Testing Garnet Site Occupancy Logic") - print(">>> Compares heuristic Fe2+/Fe3+ split against end-member site totals.") - if 'g' not in self.out.ph: - self.skipTest("Garnet not stable") + EP_PT_POINTS = [ + (10.0, 325.0), + (10.0, 350.0), + (10.0, 375.0), + (10.0, 400.0), + ] + + for P, T in EP_PT_POINTS: + out = self.calc.calculate_grid(P, T)[0] + if 'ep' not in out.ph: + continue + + split = self.calc._extract_fe_split_from_apfu(out, 'ep') + fe2_h, fe3_h = split['Fe2'], split['Fe3'] + + ph_idx = out.ph.index('ep') + ox_names = [str(o) for o in out.oxides] + comp_apfu = np.array(out.SS_vec[ph_idx].Comp_apfu, dtype=float) + fe_total = comp_apfu[ox_names.index('FeO')] - # 1. Built-in heuristic (uses excess O to split Fe) - split = self.calc._extract_fe_split_from_apfu(self.out, 'g') - fe2_h = split['fe2'] - fe3_h = split['fe3'] - mg_num_h = get_phase_mg_number(self.out, 'g') - mg2_num_h = get_phase_mg2_number(self.out, 'g') - - # 2. Site occupancy (based on user tables) - ph_idx = self.out.ph.index('g') - em = {str(n): float(f) for n, f in zip(self.out.SS_vec[ph_idx].emNames, self.out.SS_vec[ph_idx].emFrac)} - - # Fe2+ total = 3 * alm - fe2_s = 3.0 * em.get('alm', 0.0) - # Fe3+ total = 2 * kho - fe3_s = 2.0 * em.get('kho', 0.0) - - # Mg# (X-site) = Mg / (Mg + Fe2+) - # Mg_X = 3*py + 3*kho - # Fe2+_X = 3*alm - mg_x = 3.0 * em.get('py', 0.0) + 3.0 * em.get('kho', 0.0) - fe2_x = 3.0 * em.get('alm', 0.0) - mg_num_s = mg_x / (mg_x + fe2_x) if (mg_x + fe2_x) > 0 else 0.0 - - print(f"\nGarnet (g):") - print(f" Heuristic Fe2+: {fe2_h:.4f}, Fe3+: {fe3_h:.4f}, Mg#: {mg_num_h:.4f}") - print(f" Site-occ Fe2+: {fe2_s:.4f}, Fe3+: {fe3_s:.4f}, Mg# (X): {mg_num_s:.4f}") - - self.assertAlmostEqual(fe2_h, fe2_s, places=3) - self.assertAlmostEqual(fe3_h, fe3_s, places=3) - # Verify get_phase_mg2_number (divalent-only) matches site calculation - self.assertAlmostEqual(mg2_num_h, mg_num_s, places=3) + with self.subTest(P=P, T=T): + self.assertAlmostEqual(fe2_h, 0.0, places=2) + self.assertAlmostEqual(fe3_h, fe_total, places=2) - def test_clinopyroxene_comparison(self): + def test_amphibole_heuristic(self): + """Heuristic Fe2+/Fe3+ matches endmember stoichiometry for amphibole.""" + out = self.calc.calculate_grid(10.0, 550.0)[0] + if 'amp' not in out.ph: + self.skipTest("Amphibole not stable at 10 kbar, 550°C") + + split = self.calc._extract_fe_split_from_apfu(out, 'amp') + fe2_h, fe3_h = split['Fe2'], split['Fe3'] + + ph_idx = out.ph.index('amp') + em = {str(n): float(f) for n, f in zip(out.SS_vec[ph_idx].emNames, out.SS_vec[ph_idx].emFrac)} + + # Sum Fe2+ and Fe3+ from endmembers + # Use sum constraint as primary check: Fe2+ + Fe3+ = total Fe + ox_names = [str(o) for o in out.oxides] + fe_total = float(out.SS_vec[ph_idx].Comp_apfu[ox_names.index('FeO')]) + self.assertAlmostEqual(fe2_h + fe3_h, fe_total, places=2) + + def test_biotite_heuristic(self): + """Heuristic Fe2+/Fe3+ matches endmember stoichiometry for biotite.""" + out = self.calc.calculate_grid(10.0, 550.0)[0] + if 'bi' not in out.ph: + self.skipTest("Biotite not stable at 10 kbar, 550°C") + + split = self.calc._extract_fe_split_from_apfu(out, 'bi') + fe2_h, fe3_h = split['Fe2'], split['Fe3'] + + ph_idx = out.ph.index('bi') + ox_names = [str(o) for o in out.oxides] + fe_total = float(out.SS_vec[ph_idx].Comp_apfu[ox_names.index('FeO')]) + self.assertAlmostEqual(fe2_h + fe3_h, fe_total, places=2) + + def test_chlorite_heuristic(self): + """Heuristic Fe2+/Fe3+ for chlorite (not stable for this bulk composition). + + Chlorite is not stable at any P-T for the eclogite S10 bulk in mpe. + This test is kept for completeness but will always skip. """ - Compare Clinopyroxene Fe2+ and Mg# (M1m-site) from heuristic vs site occupancy. + out = self.calc.calculate_grid(5.0, 450.0)[0] + if 'chl' not in out.ph: + self.skipTest("Chlorite not stable at 5 kbar, 450°C") - In the Clinopyroxene model (Green et al. 2016), Mg and Fe2+ partition - between the M1m and M1a sites. The Bulk Mg# and Site Mg# (M1m) will - DIFFER because of ordered intermediate end-members: + split = self.calc._extract_fe_split_from_apfu(out, 'chl') + fe2_h, fe3_h = split['Fe2'], split['Fe3'] - 1. 'om' (Omphacite): Partition Mg into M1m and Al into M1a. - 2. 'cfm': Partitions Fe into M1m and Mg into M1a. + ph_idx = out.ph.index('chl') + ox_names = [str(o) for o in out.oxides] + fe_total = float(out.SS_vec[ph_idx].Comp_apfu[ox_names.index('FeO')]) + self.assertAlmostEqual(fe2_h + fe3_h, fe_total, places=2) + + def test_bulk_fe_o_closure(self): + """Bulk Fe and O are conserved across the assemblage (mass balance).""" + out = self.calc.calculate_grid(20.0, 650.0)[0] - Because 'cfm' puts Mg on the M1a site, the M1m-site Mg# will be lower - than the Bulk Mg#. For thermometry (Kd calculations), the BULK Mg# - should always be used to remain consistent with empirical calibrations - (e.g., Ellis & Green) and to avoid bias from cation ordering. - """ - print("\n>>> Testing Clinopyroxene Site Occupancy Logic") - print(">>> Compares heuristic Fe2+/Fe3+ split against ordered site totals (M1m).") - ph_name = 'dio' if 'dio' in self.out.ph else ('omph' if 'omph' in self.out.ph else None) - if ph_name is None: - self.skipTest("Clinopyroxene not stable") + total_fe_weighted = 0.0 + total_o_weighted = 0.0 + + fe_bearing_phases = ['g', 'dio', 'omph', 'amp', 'bi', 'chl', 'ep', 'ilm', 'sp'] + + for ph in out.ph: + if ph not in fe_bearing_phases: + continue - # 1. Built-in heuristic - split = self.calc._extract_fe_split_from_apfu(self.out, ph_name) - fe2_h = split['fe2'] - fe3_h = split['fe3'] - mg_num_h = get_phase_mg_number(self.out, ph_name) - mg2_num_h = get_phase_mg2_number(self.out, ph_name) - - # 2. Site occupancy (based on user tables) - ph_idx = self.out.ph.index(ph_name) - em = {str(n): float(f) for n, f in zip(self.out.SS_vec[ph_idx].emNames, self.out.SS_vec[ph_idx].emFrac)} - - # Fe2+ total = 0.5*hed + 0.5*cfm (M1m + M1a) - # Note: hed formula CaFeSi2O6 has 1.0 Fe total (0.5 in M1m, 0.5 in M1a) - # cfm formula CaMg.5Fe.5SiO6 has 0.5 Fe total (0.5 in M1m, 0.0 in M1a) - # Total atoms: - fe2_s = 1.0 * em.get('hed', 0.0) + 0.5 * em.get('cfm', 0.0) - # Fe3+ total = 1.0*acmm + 0.5*jac - fe3_s = 1.0 * em.get('acmm', 0.0) + 0.5 * em.get('jac', 0.0) - - # Mg# (M1m-site) = Mg / (Mg + Fe2+) - # Mg_M1m = 0.5*di + 0.5*om - # Fe2+_M1m = 0.5*hed + 0.5*cfm - mg_m1m = 0.5 * em.get('di', 0.0) + 0.5 * em.get('om', 0.0) - fe2_m1m = 0.5 * em.get('hed', 0.0) + 0.5 * em.get('cfm', 0.0) - mg_num_s = mg_m1m / (mg_m1m + fe2_m1m) if (mg_m1m + fe2_m1m) > 0 else 0.0 - - print(f"\nClinopyroxene ({ph_name}):") - print(f" Heuristic Fe2+: {fe2_h:.4f}, Fe3+: {fe3_h:.4f}, Mg#: {mg_num_h:.4f}") - print(f" Site-occ Fe2+: {fe2_s:.4f}, Fe3+: {fe3_s:.4f}, Mg# (M1m): {mg_num_s:.4f}") - - self.assertAlmostEqual(fe2_h, fe2_s, places=3) - self.assertAlmostEqual(fe3_h, fe3_s, places=3) - - # Verify get_phase_mg2_number (divalent-only) matches bulk fe2-basis Mg# - # For CPX, we compare to the bulk fe2 result, not site-M1m result. - mgo_idx = [str(o) for o in self.out.oxides].index('MgO') - mg_bulk = float(self.out.SS_vec[ph_idx].Comp_apfu[mgo_idx]) - mg_bulk_fe2 = mg_bulk / (mg_bulk + fe2_h) - self.assertAlmostEqual(mg2_num_h, mg_bulk_fe2, places=3) + try: + frac = phase_frac(ph, out, 'mol') + ph_idx = out.ph.index(ph) + ox_names = [str(o) for o in out.oxides] + feo_frac = float(out.SS_vec[ph_idx].Comp[ox_names.index('FeO')]) + o_frac = float(out.SS_vec[ph_idx].Comp[ox_names.index('O')]) + + total_fe_weighted += feo_frac * frac + total_o_weighted += o_frac * frac + except (ValueError, IndexError): + continue + + bulk_feo = self.calc.X[self.calc.Xoxides.index('FeO')] + bulk_o = self.calc.X[self.calc.Xoxides.index('O')] + + # Bulk X is in mol% (sums to 100); the weighted Comp sums give X/100. + if bulk_feo > 0: + self.assertAlmostEqual(total_fe_weighted / (bulk_feo / 100.0), 1.0, delta=0.1) + if bulk_o > 0: + self.assertAlmostEqual(total_o_weighted / (bulk_o / 100.0), 1.0, delta=0.1) if __name__ == '__main__': unittest.main() + +@unittest.skipUnless(HAS_JULIA, "Julia + MAGEMin_C not available") +class TestHeuristicAcrossDatabases(unittest.TestCase): + """Test Fe2+/Fe3+ heuristic across different thermodynamic databases.""" + + def test_ig_database(self): + """Heuristic works for ig (igneous) database.""" + calc = MAGEMinPTGridCalculator(db='ig') + Xoxides = ['SiO2', 'Al2O3', 'CaO', 'MgO', 'FeO', 'TiO2', 'O'] + X = [50.0, 10.0, 10.0, 15.0, 10.0, 1.0, 0.5] + calc.setup_bulk_composition(Xoxides, X, sys_in='mol') + out = calc.calculate_grid(20.0, 700.0)[0] + + for ph in ['g', 'dio']: + if ph not in out.ph: + continue + split = calc._extract_fe_split_from_apfu(out, ph) + ph_idx = out.ph.index(ph) + ox_names = [str(o) for o in out.oxides] + fe_total = float(out.SS_vec[ph_idx].Comp_apfu[ox_names.index('FeO')]) + self.assertAlmostEqual(split['Fe2'] + split['Fe3'], fe_total, places=2) + + def test_mp_database(self): + """Heuristic works for mp (metapelite) database.""" + calc = MAGEMinPTGridCalculator(db='mp') + Xoxides = ['H2O', 'SiO2', 'Al2O3', 'CaO', 'MgO', 'FeO', 'K2O', 'Na2O', 'TiO2', 'MnO', 'O'] + X = [0.92, 54.57, 8.79, 11.20, 8.45, 12.89, 0.24, 2.24, 1.12, 0.22, 0.64] + calc.setup_bulk_composition(Xoxides, X, sys_in='mol') + out = calc.calculate_grid(30.0, 700.0)[0] + + for ph in ['g', 'dio']: + if ph not in out.ph: + continue + split = calc._extract_fe_split_from_apfu(out, ph) + ph_idx = out.ph.index(ph) + ox_names = [str(o) for o in out.oxides] + fe_total = float(out.SS_vec[ph_idx].Comp_apfu[ox_names.index('FeO')]) + self.assertAlmostEqual(split['Fe2'] + split['Fe3'], fe_total, places=2) + + def test_mb_database(self): + """Heuristic works for mb (metabasite) database.""" + calc = MAGEMinPTGridCalculator(db='mb') + Xoxides = ['H2O', 'SiO2', 'Al2O3', 'CaO', 'MgO', 'FeO', 'Na2O', 'TiO2', 'O'] + X = [1.0, 50.0, 15.0, 12.0, 10.0, 8.0, 2.0, 1.0, 0.5] + calc.setup_bulk_composition(Xoxides, X, sys_in='mol') + out = calc.calculate_grid(15.0, 600.0)[0] + + for ph in ['g', 'dio']: + if ph not in out.ph: + continue + split = calc._extract_fe_split_from_apfu(out, ph) + ph_idx = out.ph.index(ph) + ox_names = [str(o) for o in out.oxides] + fe_total = float(out.SS_vec[ph_idx].Comp_apfu[ox_names.index('FeO')]) + self.assertAlmostEqual(split['Fe2'] + split['Fe3'], fe_total, places=2) + diff --git a/tests/test_solvus_instances.py b/tests/test_solvus_instances.py new file mode 100644 index 0000000..ac35947 --- /dev/null +++ b/tests/test_solvus_instances.py @@ -0,0 +1,365 @@ +"""Mock-based tests for solvus (multi-instance phase) handling. + +MAGEMin reports coexisting solvus limbs as repeated entries in ``out.ph`` +(e.g. two clinopyroxenes ``dio``, or two amphiboles ``amp``). + +Two layers are tested: + +* the low-level composition helpers (``get_oxide_apfu``, + ``extract_end_member``, ...) take ``instance`` (integer index or + ``'all'`` for per-instance numpy arrays); +* the grid-level ``extract_from_grid`` / ``single_point_calc`` return a + list of per-instance bundles indexed ``res[0]``, ``res[1]`` ... Bundle + keys are unsuffixed (``mol_frac``, ``ox_apfu_Na2O``, ``em_py``, ...), + and each value is an array over the grid points (NaN where the phase + or that limb is absent). The shape is the same for every phase; a + single-instance phase just has a one-element list. + +No live Julia runtime is needed -- ``out`` objects are mocked. +""" + +import unittest +import warnings +import numpy as np +from unittest.mock import MagicMock, patch + +from phasetools.calculators.pt_grid import MAGEMinPTGridCalculator +from phasetools.core.phase_properties import ( + get_oxide_apfu, get_phase_chemistry, extract_end_member, + get_phase_mg_number, get_phase_mg2_number, get_phase_fe_split, + _phase_indices, +) + +# mpe-style oxide ordering +OXIDES = ['H2O', 'SiO2', 'Al2O3', 'CaO', 'MgO', 'FeO', 'K2O', 'Na2O', + 'TiO2', 'MnO', 'O'] + + +def _make_ss(apfu, em_names, em_frac): + """Build a mocked SS_vec entry (one phase instance).""" + p = MagicMock() + p.Comp_apfu = apfu + # molar fractions (0-1) aligned with OXIDES + p.Comp = [0.0, 0.6, 0.03, 0.15, 0.13, 0.02, 0.0, 0.06, 0.0, 0.01, 0.0] + p.Comp_wt = [0.0, 0.55, 0.05, 0.12, 0.10, 0.02, 0.0, 0.05, 0.0, 0.01, 0.0] + p.emNames = em_names + p.emFrac = em_frac + p.emFrac_wt = [f * 0.9 for f in em_frac] + return p + + +def _dio_0(): + # diopside-rich limb: low Na + apfu = [0.0, 1.98, 0.05, 0.74, 0.64, 0.12, 0.0, 0.26, 0.0, 0.01, 6.0] + return _make_ss(apfu, ['di', 'hed', 'om', 'jac'], [0.60, 0.10, 0.25, 0.05]) + + +def _dio_1(): + # omphacite-rich limb: high Na + apfu = [0.0, 1.96, 0.09, 0.59, 0.51, 0.18, 0.0, 0.41, 0.0, 0.01, 6.0] + return _make_ss(apfu, ['di', 'hed', 'om', 'jac'], [0.40, 0.10, 0.40, 0.10]) + + +def _garnet(): + apfu = [0.0, 3.02, 2.0, 0.55, 0.70, 1.70, 0.0, 0.0, 0.0, 0.03, 12.0] + return _make_ss(apfu, ['py', 'alm', 'gr', 'spss'], [0.20, 0.60, 0.16, 0.04]) + + +def _make_out(ph_list, ss_vec, ph_frac=None): + out = MagicMock() + out.ph = ph_list + out.oxides = OXIDES + n = len(ph_list) + if ph_frac is None: + ph_frac = [0.1] * n + out.ph_frac = ph_frac + out.ph_frac_wt = list(ph_frac) + out.ph_frac_vol = list(ph_frac) + out.SS_vec = ss_vec + return out + + +class TestPhaseIndices(unittest.TestCase): + """_phase_indices resolution.""" + + def setUp(self): + self.out = _make_out(['dio', 'q', 'dio', 'g'], [_dio_0(), None, _dio_1(), _garnet()]) + + def test_all_returns_every_occurrence(self): + self.assertEqual(_phase_indices(self.out, 'dio', 'all'), [0, 2]) + + def test_default_zero_warns_with_others(self): + with self.assertWarnsRegex(UserWarning, r"2 instances.*1.*others"): + self.assertEqual(_phase_indices(self.out, 'dio', 0), [0]) + + def test_single_instance_does_not_warn(self): + with warnings.catch_warnings(): + warnings.simplefilter("error") + self.assertEqual(_phase_indices(self.out, 'g', 0), [3]) + + def test_out_of_range_warns_and_empty(self): + with self.assertWarnsRegex(UserWarning, r"instance 2 does not exist"): + self.assertEqual(_phase_indices(self.out, 'dio', 2), []) + + def test_absent_phase_empty(self): + self.assertEqual(_phase_indices(self.out, 'ep', 0), []) + + def test_invalid_instance_type_raises(self): + with self.assertRaises(ValueError): + _phase_indices(self.out, 'dio', 'first') + + +class TestSolvusHelpers(unittest.TestCase): + """Per-instance behaviour of the composition helpers.""" + + def setUp(self): + self.out = _make_out(['dio', 'q', 'dio', 'g'], [_dio_0(), None, _dio_1(), _garnet()]) + + def test_oxide_apfu_default_is_first(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO']) + self.assertAlmostEqual(apfu['Na2O'], 0.26) + self.assertAlmostEqual(apfu['MgO'], 0.64) + + def test_oxide_apfu_second_instance(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO'], instance=1) + self.assertAlmostEqual(apfu['Na2O'], 0.41) + self.assertAlmostEqual(apfu['MgO'], 0.51) + + def test_oxide_apfu_all_returns_arrays(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO'], instance='all') + self.assertTrue(np.allclose(apfu['Na2O'], [0.26, 0.41])) + self.assertTrue(np.allclose(apfu['MgO'], [0.64, 0.51])) + + def test_phase_chemistry_second_instance(self): + chem = get_phase_chemistry(self.out, 'dio', ['Na2O'], 'mol', instance=1) + self.assertAlmostEqual(chem['Na2O'], 0.06 * 100.0) # Comp[7] * 100 + + def test_extract_end_member_per_instance(self): + self.assertAlmostEqual(extract_end_member('dio', self.out, 'di', 'mol'), 0.60) + self.assertAlmostEqual(extract_end_member('dio', self.out, 'di', 'mol', instance=1), 0.40) + vals = extract_end_member('dio', self.out, 'di', 'mol', instance='all') + self.assertTrue(np.allclose(vals, [0.60, 0.40])) + + def test_mg_number_per_instance(self): + # Mg# = Mg/(Mg+Fe); first limb Mg-rich, second more Fe-rich + self.assertGreater(get_phase_mg_number(self.out, 'dio'), + get_phase_mg_number(self.out, 'dio', instance=1)) + vals = get_phase_mg_number(self.out, 'dio', instance='all') + self.assertEqual(vals.shape, (2,)) + self.assertAlmostEqual(vals[0], get_phase_mg_number(self.out, 'dio')) + + def test_mg2_number_per_instance(self): + m0 = get_phase_mg2_number(self.out, 'dio') + m1 = get_phase_mg2_number(self.out, 'dio', instance=1) + self.assertGreater(m0, m1) + vals = get_phase_mg2_number(self.out, 'dio', instance='all') + self.assertEqual(vals.shape, (2,)) + + def test_fe_split_per_instance(self): + # O-basis: excess O = 0, so Fe2 = total Fe, Fe3 = 0 for both limbs + s0 = get_phase_fe_split(self.out, 'dio') + s1 = get_phase_fe_split(self.out, 'dio', instance=1) + self.assertAlmostEqual(s0['Fe2'], 0.12) + self.assertAlmostEqual(s1['Fe2'], 0.18) + all_s = get_phase_fe_split(self.out, 'dio', instance='all') + self.assertTrue(np.allclose(all_s['Fe2'], [0.12, 0.18])) + + def test_out_of_range_returns_zeros(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O'], instance=5) + self.assertEqual(apfu['Na2O'], 0.0) + + +class TestExtractFromGridSolvus(unittest.TestCase): + """Grid extraction: list of per-instance bundles (res[0], res[1], ...).""" + + def _make_calc(self, grid_out): + calc = MAGEMinPTGridCalculator.__new__(MAGEMinPTGridCalculator) + calc.sys_in = 'mol' + calc.last_grid_out = grid_out + calc.rm_list = None + return calc + + def _dio_grid(self): + # point 0: dio solvus (2 limbs); point 1: single dio + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + g1 = _make_out(['q', 'g', 'dio'], + [None, _garnet(), _dio_0()], + ph_frac=[0.5, 0.2, 0.3]) + return [g0, g1] + + def test_uniform_shape_single_and_solvus(self): + """A single-instance phase and a solvus return the same bundle keys.""" + calc = self._make_calc(self._dio_grid()) + solvus = calc.extract_from_grid('dio', oxides=['Na2O', 'MgO'], + mg_number=True, fe_split=True) + single = calc.extract_from_grid('g', oxides=['MgO'], + mg_number=True, fe_split=True) + + # dio: two bundles; garnet: one bundle (always at res[0]) + self.assertEqual(len(solvus), 2) + self.assertEqual(len(single), 1) + + for bundle in solvus + single: + for key in ('mol_frac', 'wt_frac', 'vol_frac', + 'ox_apfu_Na2O' if 'ox_apfu_Na2O' in bundle else 'ox_apfu_MgO', + 'Mg_number', 'Fe2', 'Fe3'): + self.assertIn(key, bundle) + # unsuffixed keys inside each bundle (no _0/_1, no total_*) + self.assertIn('mol_frac', solvus[0]) + self.assertNotIn('mol_frac_0', solvus[0]) + self.assertNotIn('total_mol_frac', solvus[0]) + + def test_limb_bundles(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', oxides=['Na2O']) + c0, c1 = res + # limb fractions (0.3/0.1 at point 0; point 1 has only limb 0) + self.assertTrue(np.allclose(c0['mol_frac'], [0.3, 0.3])) + self.assertAlmostEqual(c0['wt_frac'][0], 0.3) + self.assertAlmostEqual(c0['vol_frac'][0], 0.3) + self.assertAlmostEqual(c1['mol_frac'][0], 0.1) + self.assertAlmostEqual(c0['ox_apfu_Na2O'][0], 0.26) + self.assertAlmostEqual(c1['ox_apfu_Na2O'][0], 0.41) + # user can sum the limbs themselves + total = c0['mol_frac'] + c1['mol_frac'] + self.assertAlmostEqual(total[0], 0.4) + + def test_nan_when_limb_absent(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) + c0, c1 = res + # point 1 has a single dio -> limb-1 values NaN + self.assertTrue(np.isnan(c1['mol_frac'][1])) + self.assertTrue(np.isnan(c1['wt_frac'][1])) + self.assertTrue(np.isnan(c1['ox_apfu_Na2O'][1])) + self.assertTrue(np.isnan(c1['Mg_number'][1])) + # limb 0 is filled + self.assertAlmostEqual(c0['mol_frac'][1], 0.3) + self.assertAlmostEqual(c0['ox_apfu_Na2O'][1], 0.26) + + def test_absent_phase_bundles_nan(self): + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + g2 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc([g0, g2]) + res = calc.extract_from_grid('dio', oxides=['Na2O']) + c0, c1 = res + # point 1: dio absent -> NaN in every bundle + self.assertTrue(np.isnan(c0['mol_frac'][1])) + self.assertTrue(np.isnan(c1['mol_frac'][1])) + self.assertTrue(np.isnan(c0['ox_apfu_Na2O'][1])) + + def test_absent_everywhere_empty(self): + """Phase never stable -> empty list.""" + g0 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + g1 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.5, 0.5]) + calc = self._make_calc([g0, g1]) + res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) + self.assertEqual(res, []) + + def test_single_instance_phase_filled(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('g', oxides=['MgO']) + bundle = res[0] + self.assertAlmostEqual(bundle['mol_frac'][0], 0.4) + self.assertAlmostEqual(bundle['ox_apfu_MgO'][0], 0.70) + + def test_cations_per_instance(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', cations=['Mg', 'Fe']) + c0, c1 = res + self.assertIn('cat_Mg', c0) + self.assertIn('cat_Mg', c1) + # first limb is Mg-richer + self.assertGreater(c0['cat_Mg'][0], c1['cat_Mg'][0]) + + def test_end_members_per_instance(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', end_members=['di', 'om']) + c0, c1 = res + self.assertAlmostEqual(c0['em_di'][0], 0.60) + self.assertAlmostEqual(c1['em_di'][0], 0.40) + self.assertAlmostEqual(c0['em_om'][0], 0.25) + self.assertAlmostEqual(c1['em_om'][0], 0.40) + + +class TestSinglePointCalcSolvus(unittest.TestCase): + """single_point_calc returns the same nested-bundle schema.""" + + def _make_calc(self): + calc = MAGEMinPTGridCalculator.__new__(MAGEMinPTGridCalculator) + calc.sys_in = 'mol' + calc.data = None + calc.X = None + calc.Xoxides = None + calc.rm_list = None + return calc + + def test_solvus_returns_limb_bundles(self): + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g0 + res, out = calc.single_point_calc(10.0, 600.0, 'dio', + oxides=['Na2O', 'MgO']) + self.assertEqual(len(res), 2) + c0, c1 = res + self.assertAlmostEqual(c0['mol_frac'], 0.3) + self.assertAlmostEqual(c1['mol_frac'], 0.1) + self.assertAlmostEqual(c0['wt_frac'], 0.3) + self.assertAlmostEqual(c1['wt_frac'], 0.1) + self.assertAlmostEqual(c0['ox_apfu_Na2O'], 0.26) + self.assertAlmostEqual(c1['ox_apfu_Na2O'], 0.41) + self.assertIs(out, g0) + + def test_single_instance_has_one_bundle(self): + g = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g + res, _ = calc.single_point_calc(10.0, 600.0, 'g', oxides=['MgO']) + self.assertEqual(len(res), 1) + bundle = res[0] + self.assertAlmostEqual(bundle['mol_frac'], 0.4) + self.assertAlmostEqual(bundle['ox_apfu_MgO'], 0.70) + + def test_absent_phase_no_bundles(self): + g = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g + res, _ = calc.single_point_calc(10.0, 600.0, 'dio') + self.assertEqual(res, []) + + +class TestGarnetEndmembersSuffix(unittest.TestCase): + """generate_2D_grid_gt_endmembers returns the single-instance bundle.""" + + def test_returns_bundle_with_historical_keys(self): + from phasetools.calculators.garnet import MAGEMinGarnetCalculator + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + bundle = { + 'mol_frac': np.array([0.10, 0.15]), + 'wt_frac': np.array([0.11, 0.16]), + 'vol_frac': np.array([0.12, 0.17]), + 'em_py': np.array([0.20, 0.30]), + 'em_alm': np.array([0.60, 0.50]), + 'em_gr': np.array([0.16, 0.17]), + 'em_spss': np.array([0.04, 0.03]), + } + with patch.object(calc, 'calculate_grid', return_value=None), \ + patch.object(calc, 'extract_from_grid', return_value=[bundle]): + res = calc.generate_2D_grid_gt_endmembers([10.0], [600.0]) + self.assertIn('em_py', res) + self.assertIn('em_alm', res) + self.assertIn('mol_frac', res) + self.assertTrue(np.allclose(res['em_py'], [0.20, 0.30])) + + +if __name__ == '__main__': + unittest.main()