diff --git a/src/phasetools/models/magma_ocean.py b/src/phasetools/models/magma_ocean.py index 78f67c7..4106300 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 @@ -61,7 +62,7 @@ def depth_to_pressure(self, depth_km: float) -> float: return (self.rho_avg * self.g * depth_km * 1000.0) / 1e8 def radius_to_pressure(self, R_km: float) -> float: - """Convert radius from center (km) to pressure (kbar).""" + """Convert radius from centre (km) to pressure (kbar).""" depth = self.radius_body - R_km return self.depth_to_pressure(depth) @@ -89,7 +90,15 @@ 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}." + ) + raise RuntimeError( + f"find_temperature_at_vol_frac: bisection failed for P={P}, " + f"target_vol_frac={target_vol_frac}. Bracket T range: {bracket}" + ) 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 +130,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 +155,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 +205,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 +225,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/tests/test_fractionation.py b/tests/test_fractionation.py new file mode 100644 index 0000000..e601699 --- /dev/null +++ b/tests/test_fractionation.py @@ -0,0 +1,164 @@ +"""Mock-based tests for MagmaOcean fractional crystallisation fixes. + +These tests verify: +- 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. +""" + +import unittest +import numpy as np +from unittest.mock import MagicMock, patch, PropertyMock + + +# =========================================================================== +# 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)) + + +# =========================================================================== +# run_stage_0: no melt / no solid edge case +# =========================================================================== +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"], {}) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_solvus_instances.py b/tests/test_solvus_instances.py index ff1d058..ac35947 100644 --- a/tests/test_solvus_instances.py +++ b/tests/test_solvus_instances.py @@ -337,39 +337,6 @@ def test_absent_phase_no_bundles(self): self.assertEqual(res, []) -class TestFeSplitMixedBasis(unittest.TestCase): - """get_phase_fe_split uses the correct basis per instance.""" - - def test_mixed_o_and_traditional_basis(self): - """One O-basis instance and one traditional-basis instance are split correctly.""" - from phasetools.core.phase_properties import get_phase_fe_split - - # Two instances of 'sp': O-basis ignores Fe2O3 for total Fe; - # traditional basis counts Fe2O3 in total Fe. - mock_out = MagicMock() - mock_out.ph = ['sp', 'sp'] - mock_out.oxides = ['MgO', 'FeO', 'Fe2O3', 'O'] - - mock_sp_0 = MagicMock() # O-basis - mock_sp_0.Comp_apfu = [1.0, 2.0, 0.5, 6.0] - - mock_sp_1 = MagicMock() # Traditional basis - mock_sp_1.Comp_apfu = [1.0, 1.0, 1.0, 0.0] - - mock_out.SS_vec = [mock_sp_0, mock_sp_1] - - split = get_phase_fe_split(mock_out, 'sp', instance='all') - fe2 = split['Fe2'] - fe3 = split['Fe3'] - - # Instance 0: O-basis, total_fe = feo = 2.0; Fe3 = 2*0.5 + 0 = 1.0; Fe2 = 1.0 - self.assertAlmostEqual(fe2[0], 1.0) - self.assertAlmostEqual(fe3[0], 1.0) - # Instance 1: traditional, total_fe = 1.0 + 2*1.0 = 3.0; Fe3 = 2*1.0 = 2.0; Fe2 = 1.0 - self.assertAlmostEqual(fe2[1], 1.0) - self.assertAlmostEqual(fe3[1], 2.0) - - class TestGarnetEndmembersSuffix(unittest.TestCase): """generate_2D_grid_gt_endmembers returns the single-instance bundle."""