diff --git a/src/phasetools/calculators/garnet.py b/src/phasetools/calculators/garnet.py index 7e594d3..8d2a30c 100644 --- a/src/phasetools/calculators/garnet.py +++ b/src/phasetools/calculators/garnet.py @@ -2,51 +2,94 @@ 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: str = "ig", dataset: int = 636, verbose: bool = False, fe_basis: str = "FeOt" + ) -> None: + """ + 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). + * ``'Fe2'`` or ``'fe2'`` -- alias for ``'Fe2+'``. + + 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 @@ -127,7 +170,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 @@ -144,6 +241,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 @@ -158,19 +256,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/tests/test_fe_basis.py b/tests/test_fe_basis.py new file mode 100644 index 0000000..125a2a5 --- /dev/null +++ b/tests/test_fe_basis.py @@ -0,0 +1,89 @@ +"""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') + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_fractionation.py b/tests/test_fractionation.py index e601699..297f9ff 100644 --- a/tests/test_fractionation.py +++ b/tests/test_fractionation.py @@ -1,10 +1,11 @@ -"""Mock-based tests for MagmaOcean fractional crystallisation fixes. +"""Mock-based tests for fractionation correctness fixes. -These tests verify: +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) -- Fix H4: starting_melt length validation -- run_stage_0: no melt / no solid edge case -- bisection failure raises RuntimeError No live Julia runtime is needed — all MAGEMin calls are mocked. """ @@ -14,150 +15,201 @@ 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 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 H1: self.X restored after run_fractional_stages +# Fix D: X_along_path normalised to 1 # =========================================================================== -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" - ) +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 H4: starting_melt length validation +# Fix C: pure-phase IndexError # =========================================================================== -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)) +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))) # =========================================================================== -# run_stage_0: no melt / no solid edge case +# Zero guard: frac_amount=0 is a no-op # =========================================================================== -class TestStageZeroNoMeltNoSolid(unittest.TestCase): - """run_stage_0 must return zeros for avg_melt and empty layer_modes when liq is absent.""" - - @patch('phasetools.models.magma_ocean.MAGEMinBase.__init__', return_value=None) - def test_no_melt_no_solid(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', 'MgO'] - mo.sys_in = 'mol' - mo.data = MagicMock() - mo.X = np.array([33.0, 33.0, 34.0]) - mo.Xoxides = MagicMock() - mo.rm_list = None - - # find_temperature_at_vol_frac is called but its return value is irrelevant - # because the mocked output has no phases. - mo.find_temperature_at_vol_frac = MagicMock(return_value=1200.0) - - mock_out = MagicMock() - mock_out.ph = [] - mock_out.n_SS = 0 - mock_out.SS_vec = [] - mock_out.PP_vec = [] - - mock_magemin_c = MagicMock() - mock_magemin_c.single_point_minimization = MagicMock(return_value=mock_out) - - with patch.object(mo_module, 'MAGEMin_C', mock_magemin_c): - results, avg_melt = mo.run_stage_0( - p_start=5.0, p_end=0.001, solid_frac=0.5, p_intervals=3 - ) - - self.assertTrue(np.all(np.isfinite(avg_melt))) - self.assertTrue(np.allclose(avg_melt, np.zeros(3))) - self.assertEqual(results["layer_modes"], {}) +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)) if __name__ == '__main__':