diff --git a/docs/How-to/usage.md b/docs/How-to/usage.md index c8a79bb16..f1bd17af4 100644 --- a/docs/How-to/usage.md +++ b/docs/How-to/usage.md @@ -296,3 +296,67 @@ Environment variables FWL_DATA: ok RAD_DIR: ok ``` + +## Melting Curves + +PROTEUS uses precomputed solidus and liquidus curves from laboratory experiments +and theoretical parametrizations of silicate melting. These curves define the +temperatures at which a silicate material begins to melt (solidus) and becomes +fully molten (liquidus) as a function of pressure. + +The melting-curve exporter generates lookup tables in both pressure–temperature +(P–T) and pressure–entropy (P–S) space for several literature parametrizations +of peridotite / silicate melting. + +### What the exporter does + +The script `tools/solidus_func.py` is designed to work with the legacy EOS lookup +tables: + +- `temperature_solid.dat` +- `temperature_melt.dat` + +These tables provide temperature as a function of entropy and pressure, + on structured grids. The exporter therefore performs the following +steps: + +1. Build solidus and liquidus curves in P–T space from literature fits. +2. Convert those curves into P–S space by inverting the EOS relation \(T(S, P)\). +3. Resample the solidus and liquidus entropy curves onto a common pressure grid. +4. Save both the P–T and P–S versions to disk for later use by PROTEUS. + +### Available parametrizations + +The following directory names are supported and should be used exactly as written +in the TOML configuration in the `melting_dir` parameter: + + +| Directory name | Reference | DOI | +|--------------------|------------------------------|-----| +| `andrault_2011` | Andrault et al. (2011) | [10.1016/j.epsl.2011.02.006](https://doi.org/10.1016/j.epsl.2011.02.006) | +| `monteux_2016` | Monteux et al. (2016) | [10.1016/j.epsl.2016.05.010](https://doi.org/10.1016/j.epsl.2016.05.010) | +| `wolf_bower_2018` | Wolf & Bower (2018) | [10.1016/j.pepi.2017.11.004](https://doi.org/10.1016/j.pepi.2017.11.004)
[10.1051/0004-6361/201935710](https://doi.org/10.1051/0004-6361/201935710) | +| `katz_2003` | Katz et al. (2003) | [10.1029/2002GC000433](https://doi.org/10.1029/2002GC000433) | +| `fei_2021` | Fei et al. (2021) | [10.1038/s41467-021-21170-y](https://doi.org/10.1038/s41467-021-21170-y) | +| `belonoshko_2005` | Belonoshko et al. (2005) | [10.1103/PhysRevLett.94.195701](https://doi.org/10.1103/PhysRevLett.94.195701) | +| `fiquet_2010` | Fiquet et al. (2010) | [10.1126/science.1192448](https://doi.org/10.1126/science.1192448) | +| `hirschmann_2000` | Hirschmann (2000) | [10.1029/2000GC000070](https://doi.org/10.1029/2000GC000070) | +| `stixrude_2014` | Stixrude (2014) | [10.1098/rsta.2013.0076](https://doi.org/10.1098/rsta.2013.0076) | +| `lin_2024` | Lin et al. (2024) | [10.1038/s41561-024-01495-1](https://doi.org/10.1038/s41561-024-01495-1) | + + +### Generate melting curves + +Before running PROTEUS, generate the lookup tables: + +```console +python tools/solidus_func.py --all +``` + +Alternatively, you can generate a single parametrization using a specific flag (e.g.--katz2003, --lin2024). + +This will compute all parametrizations, convert them to P–T and P–S space, and store them in: + +```console +$FWL_DATA/interior_lookup_tables/Melting_curves/ +``` diff --git a/input/all_options.toml b/input/all_options.toml index 2ba5c3014..1f412c9e3 100644 --- a/input/all_options.toml +++ b/input/all_options.toml @@ -167,7 +167,7 @@ version = "2.0" core_density = 10738.33 # Core density [kg m-3] core_heatcap = 880.0 # Core specific heat capacity [J K-1 kg-1] - module = "zalmoxis" # self | zalmoxis + module = "self" # self | zalmoxis update_interval = 0 # max interval (ceiling) between structure updates [yr]; 0 = only at init update_min_interval = 0 # min interval (floor) between updates [yr]; prevents thrashing update_dtmagma_frac = 0.03 # trigger on relative T_magma change (3%) diff --git a/src/proteus/utils/data.py b/src/proteus/utils/data.py index feadd6e22..31fda1846 100644 --- a/src/proteus/utils/data.py +++ b/src/proteus/utils/data.py @@ -1210,33 +1210,65 @@ def download_interior_lookuptables(clean=False): ) -def download_melting_curves(config: Config, clean=False): +def download_melting_curves(config: Config, clean: bool = False): """ - Download melting curve data + Ensure melting curve data are available locally. + + Expected layout: + interior_lookup_tables/ + Melting_curves/ + / + solidus_P-T.dat + liquidus_P-T.dat + solidus_P-S.dat + liquidus_P-S.dat """ log.debug('Download melting curve data') - dir = 'Melting_curves/' + config.interior.melting_dir + rel_dir = Path('Melting_curves') / config.interior.melting_dir data_dir = GetFWLData() / 'interior_lookup_tables' data_dir.mkdir(parents=True, exist_ok=True) - folder_dir = data_dir / dir + folder_dir = data_dir / rel_dir + if clean: safe_rm(folder_dir.as_posix()) - source_info = get_data_source_info(dir) + + # ------------------------------------------------------------------ + # Canonical flat layout: if files already exist locally, do not download. + # ------------------------------------------------------------------ + solidus_pt = folder_dir / 'solidus_P-T.dat' + liquidus_pt = folder_dir / 'liquidus_P-T.dat' + solidus_ps = folder_dir / 'solidus_P-S.dat' + liquidus_ps = folder_dir / 'liquidus_P-S.dat' + + if all(p.is_file() for p in (solidus_pt, liquidus_pt, solidus_ps, liquidus_ps)): + log.debug('Melting curve data already present locally: %s', folder_dir) + return + + # ------------------------------------------------------------------ + # Fallback: try remote source mapping. + # ------------------------------------------------------------------ + source_info = get_data_source_info(rel_dir.as_posix()) if not source_info: - raise ValueError(f'No data source mapping found for folder: {dir}') + raise ValueError( + f'No data source mapping found for folder: {rel_dir}. ' + f'Also did not find local melting curve data in: {folder_dir}' + ) download( - folder=dir, + folder=rel_dir.as_posix(), target=data_dir, osf_id=source_info['osf_project'], zenodo_id=source_info['zenodo_id'], - desc=f'Melting curve data: {dir}', + desc=f'Melting curve data: {rel_dir}', ) - # Create canonical _P-T copies from Zenodo legacy names (solidus.dat -> solidus_P-T.dat). - # Keep originals so md5sum checks don't trigger unnecessary re-downloads. + # ------------------------------------------------------------------ + # Legacy compatibility: + # - if download contains solidus.dat / liquidus.dat, treat them as P-T + # and create canonical *_P-T.dat copies + # ------------------------------------------------------------------ for stem in ('solidus', 'liquidus'): legacy = folder_dir / f'{stem}.dat' canonical = folder_dir / f'{stem}_P-T.dat' diff --git a/tests/tools/test_solidus_func.py b/tests/tools/test_solidus_func.py new file mode 100644 index 000000000..dd23dc1e4 --- /dev/null +++ b/tests/tools/test_solidus_func.py @@ -0,0 +1,1007 @@ +""" +Comprehensive unit tests for tools/solidus_func.py. + +This module tests the melting curve parametrizations and EOS inversion pipeline. +Covers: +1. All 10 literature melting curve models (parametric functions) +2. Helper functions for solidus/liquidus relationships +3. EOS inversion and resampling pipeline (with mocked interpolators) +4. Physical validation (solidus < liquidus, monotonicity, etc.) +5. Edge cases and error handling + +References: +- docs/How-to/test_infrastructure.md +- docs/How-to/test_categorization.md +- docs/How-to/test_building.md +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from tools.solidus_func import ( + _fmt_range, + andrault_2011, + belonoshko_2005, + build_common_entropy_grid, + fei_2021, + fiquet_2010, + get_melting_curves, + hirschmann_2000, + invert_to_entropy_along_profile, + katz_2003, + lin_2024, + liquidus_from_solidus_stixrude, + make_entropy_header, + make_pressure_grid, + monteux_2016, + solidus_from_liquidus_stixrude, + stixrude_2014, + truncate_to_physical_interval, + validate_entropy_export_arrays, + wolf_bower_2018, +) + +# ============================================================================= +# FIXTURES & HELPERS +# ============================================================================= + + +@pytest.fixture +def pressure_grid(): + """Standard pressure grid for testing (0–1000 GPa, 100 points).""" + return np.linspace(0.0, 1000.0, 100) + + +@pytest.fixture +def small_pressure_grid(): + """Small pressure grid for quick smoke tests (0–100 GPa, 10 points).""" + return np.linspace(0.0, 100.0, 10) + + +@pytest.fixture +def mock_eos_interpolator(): + """Mock EOS interpolator for T(S, P). Returns physically plausible T values.""" + + def mock_T_of_SP(points): + # points is shape (N, 2) with columns [S, P] + # Simple model: T = 1000 + 10*S + 100*P (linear, increases with entropy and pressure) + # Valid bounds: S in [100, 1000] J/kg/K, P in [0, 100] GPa + # Returns nans outside bounds (mimics RegularGridInterpolator behavior) + S = points[:, 0] + P = points[:, 1] + + T = 1000.0 + 10.0 * S + 100.0 * P + result = T.copy() + + # Set out-of-bounds to NaN + mask_invalid = (S < 100.0) | (S > 1000.0) | (P < 0.0) | (P > 100.0) + result[mask_invalid] = np.nan + + return result + + return mock_T_of_SP + + +@pytest.fixture +def mock_S_axes(): + """Mock entropy axes for solid and liquid EOS tables.""" + return ( + np.linspace(100.0, 1000.0, 125), # solid S axis + np.linspace(150.0, 1100.0, 95), # liquid S axis + ) + + +# ============================================================================= +# TEST: make_pressure_grid +# ============================================================================= + + +@pytest.mark.unit +class TestMakePressureGrid: + """Test pressure grid generation.""" + + def test_default_grid(self): + """Creates a 500-point grid from 0 to 1000 GPa (default).""" + P = make_pressure_grid() + assert len(P) == 500 + assert pytest.approx(P[0], rel=1e-10) == 0.0 + assert pytest.approx(P[-1], rel=1e-10) == 1000.0 + assert np.all(np.diff(P) > 0), 'Pressure grid must be monotonically increasing' + + def test_custom_range_and_count(self): + """Custom Pmin, Pmax, n parameters.""" + P = make_pressure_grid(Pmin=10.0, Pmax=100.0, n=20) + assert len(P) == 20 + assert pytest.approx(P[0], rel=1e-10) == 10.0 + assert pytest.approx(P[-1], rel=1e-10) == 100.0 + + def test_uniform_spacing(self): + """Pressure grid is uniformly spaced.""" + P = make_pressure_grid(Pmin=0.0, Pmax=100.0, n=101) + dP = np.diff(P) + assert np.allclose(dP, dP[0]), 'Spacing must be uniform' + + def test_single_point(self): + """Single-point grid returns array of length 1.""" + P = make_pressure_grid(Pmin=50.0, Pmax=50.0, n=1) + assert len(P) == 1 + assert pytest.approx(P[0], rel=1e-10) == 50.0 + + +# ============================================================================= +# TEST: Stixrude Ratio Functions +# ============================================================================= + + +@pytest.mark.unit +class TestStixrudeRatios: + """Test solidus/liquidus conversion via Stixrude ratios.""" + + @pytest.mark.parametrize( + 'T_liq_k', + [ + np.array([1500.0]), + np.array([1500.0, 2000.0, 2500.0]), + np.linspace(1200.0, 3000.0, 20), + ], + ) + def test_solidus_from_liquidus_stixrude(self, T_liq_k): + """Solidus < liquidus via Stixrude ratio.""" + T_sol = solidus_from_liquidus_stixrude(T_liq_k) + + assert T_sol.shape == T_liq_k.shape + assert np.all(T_sol > 0.0), 'Solidus must be positive' + assert np.all(T_sol < T_liq_k), 'Solidus < liquidus (physical constraint)' + + @pytest.mark.parametrize( + 'T_sol_k', + [ + np.array([1200.0]), + np.array([1200.0, 1500.0, 1800.0]), + np.linspace(1000.0, 2500.0, 20), + ], + ) + def test_liquidus_from_solidus_stixrude(self, T_sol_k): + """Liquidus from solidus via Stixrude ratio.""" + T_liq = liquidus_from_solidus_stixrude(T_sol_k) + + assert T_liq.shape == T_sol_k.shape + assert np.all(T_liq > 0.0), 'Liquidus must be positive' + assert np.all(T_liq > T_sol_k), 'Liquidus > solidus (physical constraint)' + + def test_stixrude_round_trip(self): + """Round-trip conversion: T_liq -> T_sol -> T_liq recovers original (approximately).""" + T_liq_orig = np.linspace(1500.0, 2500.0, 10) + T_sol = solidus_from_liquidus_stixrude(T_liq_orig) + T_liq_recovered = liquidus_from_solidus_stixrude(T_sol) + + assert np.allclose(T_liq_recovered, T_liq_orig, rtol=1e-10) + + +# ============================================================================= +# TEST: _fmt_range Helper +# ============================================================================= + + +@pytest.mark.unit +class TestFmtRange: + """Test formatting of array min/max ranges.""" + + def test_normal_finite_array(self): + """Format array with finite values.""" + arr = np.array([10.5, 50.0, 100.0]) + result = _fmt_range(arr, fmt='.1f') + assert '[10.5, 100.0]' in result + + def test_array_with_nans(self): + """Ignores NaN values, uses only finite elements.""" + arr = np.array([10.0, np.nan, 50.0, np.nan, 100.0]) + result = _fmt_range(arr, fmt='.1f') + assert '[10.0, 100.0]' in result + + def test_all_nans(self): + """Returns [nan, nan] for all-NaN array.""" + arr = np.array([np.nan, np.nan, np.nan]) + result = _fmt_range(arr) + assert '[nan, nan]' in result + + +# ============================================================================= +# TEST: Melting Curve Functions - Andrault (2011) +# ============================================================================= + + +@pytest.mark.unit +class TestAndrault2011: + """Test Andrault et al. (2011) melting curve.""" + + def test_solidus_basic_shape(self): + """Solidus returns correct shape for given pressure grid.""" + P, T = andrault_2011(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0), 'Temperature must be positive' + + def test_liquidus_basic_shape(self): + """Liquidus returns correct shape for given pressure grid.""" + P, T = andrault_2011(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0), 'Temperature must be positive' + + def test_solidus_less_than_liquidus_when_both_defined(self): + """Solidus < liquidus in regions where both are defined. + + Note: Andrault (2011) raw model may have unphysical regions where + the separation isn't maintained. Real usage with truncate_to_physical_interval + enforces solidus < liquidus. This tests the raw parametrization. + """ + P_sol, T_sol = andrault_2011(kind='solidus', Pmin=50.0, Pmax=100.0, n=50) + P_liq, T_liq = andrault_2011(kind='liquidus', Pmin=50.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq), 'Pressure grids must match' + # At higher pressures where both curves are well-defined + assert np.all(T_sol < T_liq), 'Solidus < liquidus in well-defined region' + + def test_temperature_increases_with_pressure(self): + """Temperature is monotonically increasing with pressure.""" + _, T_sol = andrault_2011(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + _, T_liq = andrault_2011(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + + assert np.all(np.diff(T_sol) > 0), 'Solidus must increase with pressure' + assert np.all(np.diff(T_liq) > 0), 'Liquidus must increase with pressure' + + def test_invalid_kind(self): + """Invalid kind parameter raises ValueError.""" + with pytest.raises(ValueError, match='kind must be'): + andrault_2011(kind='invalid') + + +# ============================================================================= +# TEST: Melting Curve Functions - Fei (2021) +# ============================================================================= + + +@pytest.mark.unit +class TestFei2021: + """Test Fei et al. (2021) melting curve.""" + + def test_solidus_basic_shape(self): + """Solidus returns correct shape.""" + P, T = fei_2021(kind='solidus', Pmin=1.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_liquidus_basic_shape(self): + """Liquidus returns correct shape.""" + P, T = fei_2021(kind='liquidus', Pmin=1.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_solidus_less_than_liquidus_at_valid_pressures(self): + """Solidus < liquidus in physically valid regions (fei_2021 uses stixrude ratio).""" + P_sol, T_sol = fei_2021(kind='solidus', Pmin=10.0, Pmax=100.0, n=50) + P_liq, T_liq = fei_2021(kind='liquidus', Pmin=10.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq) + # Stixrude ratio ensures solidus < liquidus by construction + assert np.all(T_sol < T_liq) + + def test_temperature_increases_with_pressure(self): + """Temperature increases monotonically with pressure.""" + _, T_sol = fei_2021(kind='solidus', Pmin=1.0, Pmax=100.0, n=50) + _, T_liq = fei_2021(kind='liquidus', Pmin=1.0, Pmax=100.0, n=50) + + assert np.all(np.diff(T_sol) > 0) + assert np.all(np.diff(T_liq) > 0) + + +# ============================================================================= +# TEST: Melting Curve Functions - Belonoshko (2005) +# ============================================================================= + + +@pytest.mark.unit +class TestBelonoshko2005: + """Test Belonoshko et al. (2005) melting curve.""" + + def test_solidus_basic_shape(self): + """Solidus returns correct shape.""" + P, T = belonoshko_2005(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_liquidus_basic_shape(self): + """Liquidus returns correct shape.""" + P, T = belonoshko_2005(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_solidus_less_than_liquidus_at_valid_pressures(self): + """Solidus < liquidus in physically valid regions (uses stixrude ratio).""" + P_sol, T_sol = belonoshko_2005(kind='solidus', Pmin=10.0, Pmax=100.0, n=50) + P_liq, T_liq = belonoshko_2005(kind='liquidus', Pmin=10.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq) + # Stixrude ratio ensures solidus < liquidus by construction + assert np.all(T_sol < T_liq) + + def test_temperature_increases_with_pressure(self): + """Temperature increases monotonically with pressure.""" + _, T_sol = belonoshko_2005(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + _, T_liq = belonoshko_2005(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + + assert np.all(np.diff(T_sol) > 0) + assert np.all(np.diff(T_liq) > 0) + + +# ============================================================================= +# TEST: Melting Curve Functions - Fiquet (2010, Piecewise) +# ============================================================================= + + +@pytest.mark.unit +class TestFiquet2010: + """Test Fiquet et al. (2010) piecewise melting curve.""" + + def test_solidus_basic_shape(self): + """Solidus returns correct shape.""" + P, T = fiquet_2010(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_liquidus_basic_shape(self): + """Liquidus returns correct shape.""" + P, T = fiquet_2010(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_solidus_less_than_liquidus_at_valid_pressures(self): + """Solidus < liquidus in physically valid regions (uses stixrude ratio).""" + P_sol, T_sol = fiquet_2010(kind='solidus', Pmin=10.0, Pmax=100.0, n=50) + P_liq, T_liq = fiquet_2010(kind='liquidus', Pmin=10.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq) + # Stixrude ratio ensures solidus < liquidus by construction + assert np.all(T_sol < T_liq) + + def test_piecewise_transition(self): + """Piecewise branches are smooth and the known jump at 20 GPa is explicit.""" + P, T = fiquet_2010(kind='liquidus', Pmin=18.0, Pmax=22.0, n=100) + + # Branch-wise smoothness: verify monotonic increase away from the 20 GPa split. + low = P < 20.0 + high = P > 20.0 + assert np.all(np.diff(T[low]) > 0.0), ( + 'Low-pressure branch should increase with pressure' + ) + assert np.all(np.diff(T[high]) > 0.0), ( + 'High-pressure branch should increase with pressure' + ) + + # The implementation uses two branches that are not exactly continuous at 20 GPa. + # Measure the one-sided jump from function output (left/right of the split): + # jump = T(20+)-T(20-) should be negative with magnitude ~0.56 K. + _, T_left = fiquet_2010(kind='liquidus', Pmin=20.0, Pmax=20.0, n=1) + _, T_right = fiquet_2010(kind='liquidus', Pmin=20.0 + 1e-6, Pmax=20.0 + 1e-6, n=1) + jump_at_20 = float(T_right[0] - T_left[0]) + + assert jump_at_20 < 0.0, 'Expected a downward jump at 20 GPa (T_right < T_left)' + assert abs(jump_at_20) == pytest.approx(0.55908, rel=1e-2) + + +# ============================================================================= +# TEST: Melting Curve Functions - Monteux (2016, Piecewise) +# ============================================================================= + + +@pytest.mark.unit +class TestMonteux2016: + """Test Monteux et al. (2016) piecewise melting curve.""" + + def test_solidus_basic_shape(self): + """Solidus returns correct shape.""" + P, T = monteux_2016(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_liquidus_basic_shape(self): + """Liquidus returns correct shape.""" + P, T = monteux_2016(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_solidus_less_than_liquidus(self): + """Solidus < liquidus at all pressures.""" + P_sol, T_sol = monteux_2016(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + P_liq, T_liq = monteux_2016(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq) + assert np.all(T_sol < T_liq) + + def test_piecewise_transition(self): + """Piecewise function has well-defined behavior at transition.""" + P, T = monteux_2016(kind='solidus', Pmin=15.0, Pmax=25.0, n=100) + + # Check that temperature changes appropriately (no large jumps) + dT = np.diff(T) + finite_dT = dT[np.isfinite(dT)] + if len(finite_dT) > 0: + # Most gradients should be small and positive for smooth parametrization + assert np.median(finite_dT) > 0, 'Median gradient should be positive' + assert np.max(np.abs(finite_dT)) < 100, 'No unreasonable jumps in T' + + +# ============================================================================= +# TEST: Melting Curve Functions - Hirschmann (2000) +# ============================================================================= + + +@pytest.mark.unit +class TestHirschmann2000: + """Test Hirschmann (2000) melting curve (low pressure, ~0–10 GPa).""" + + def test_solidus_basic_shape(self): + """Solidus returns correct shape.""" + P, T = hirschmann_2000(kind='solidus', Pmin=0.0, Pmax=10.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_liquidus_basic_shape(self): + """Liquidus returns correct shape.""" + P, T = hirschmann_2000(kind='liquidus', Pmin=0.0, Pmax=10.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_solidus_less_than_liquidus(self): + """Solidus < liquidus at all pressures.""" + P_sol, T_sol = hirschmann_2000(kind='solidus', Pmin=0.0, Pmax=10.0, n=50) + P_liq, T_liq = hirschmann_2000(kind='liquidus', Pmin=0.0, Pmax=10.0, n=50) + + assert np.allclose(P_sol, P_liq) + assert np.all(T_sol < T_liq) + + def test_temperature_increases_with_pressure(self): + """Temperature increases monotonically with pressure.""" + _, T_sol = hirschmann_2000(kind='solidus', Pmin=0.0, Pmax=10.0, n=50) + _, T_liq = hirschmann_2000(kind='liquidus', Pmin=0.0, Pmax=10.0, n=50) + + assert np.all(np.diff(T_sol) > 0) + assert np.all(np.diff(T_liq) > 0) + + +# ============================================================================= +# TEST: Melting Curve Functions - Stixrude (2014) +# ============================================================================= + + +@pytest.mark.unit +class TestStixrude2014: + """Test Stixrude (2014) melting curve.""" + + def test_solidus_basic_shape(self): + """Solidus returns correct shape.""" + P, T = stixrude_2014(kind='solidus', Pmin=1.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_liquidus_basic_shape(self): + """Liquidus returns correct shape.""" + P, T = stixrude_2014(kind='liquidus', Pmin=1.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_solidus_less_than_liquidus_at_valid_pressures(self): + """Solidus < liquidus in physically valid regions (uses stixrude ratio).""" + P_sol, T_sol = stixrude_2014(kind='solidus', Pmin=10.0, Pmax=100.0, n=50) + P_liq, T_liq = stixrude_2014(kind='liquidus', Pmin=10.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq) + # Stixrude ratio ensures solidus < liquidus by construction + assert np.all(T_sol < T_liq) + + +# ============================================================================= +# TEST: Melting Curve Functions - Wolf & Bower (2018, Piecewise) +# ============================================================================= + + +@pytest.mark.unit +class TestWolfBower2018: + """Test Wolf & Bower (2018) piecewise melting curve.""" + + def test_solidus_basic_shape(self): + """Solidus returns correct shape.""" + P, T = wolf_bower_2018(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_liquidus_basic_shape(self): + """Liquidus returns correct shape.""" + P, T = wolf_bower_2018(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_solidus_less_than_liquidus(self): + """Solidus < liquidus at all pressures.""" + P_sol, T_sol = wolf_bower_2018(kind='solidus', Pmin=0.0, Pmax=100.0, n=50) + P_liq, T_liq = wolf_bower_2018(kind='liquidus', Pmin=0.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq) + assert np.all(T_sol < T_liq) + + def test_piecewise_continuity(self): + """Piecewise function is continuous across transitions.""" + P, T = wolf_bower_2018(kind='solidus', Pmin=0.0, Pmax=100.0, n=200) + + # Check temperature is generally monotonically increasing (with small numerical tolerance) + dT = np.diff(T) + assert np.count_nonzero(dT > -1e-6) > len(dT) * 0.99, 'Temperature should increase' + + +# ============================================================================= +# TEST: Melting Curve Functions - Katz (2003, Hydrous) +# ============================================================================= + + +@pytest.mark.unit +class TestKatz2003: + """Test Katz et al. (2003) hydrous melting curve.""" + + def test_solidus_dry(self): + """Solidus with X_h2o = 0 matches dry Wolf & Bower (2018) at valid pressures.""" + P_sol_katz, T_sol_katz = katz_2003( + kind='solidus', X_h2o=0.0, Pmin=50.0, Pmax=100.0, n=50 + ) + P_wb, T_wb = wolf_bower_2018(kind='solidus', Pmin=50.0, Pmax=100.0, n=50) + + assert np.allclose(T_sol_katz, T_wb, rtol=1e-10) + + def test_hydrous_effect(self): + """Increasing water content decreases melting temperature (physical).""" + P, T_dry = katz_2003(kind='solidus', X_h2o=0.0, Pmin=50.0, Pmax=100.0, n=50) + _, T_wet = katz_2003(kind='solidus', X_h2o=30.0, Pmin=50.0, Pmax=100.0, n=50) + + assert np.all(T_wet < T_dry), 'Water lowers melting temperature' + + def test_default_water_content(self): + """Default water content is X_h2o = 30 ppm.""" + P1, T1 = katz_2003(kind='solidus', Pmin=50.0, Pmax=100.0, n=50) + P2, T2 = katz_2003(kind='solidus', X_h2o=30.0, Pmin=50.0, Pmax=100.0, n=50) + + assert np.allclose(T1, T2) + + def test_physical_constraint_solidus_less_than_liquidus(self): + """Solidus < liquidus with hydration at valid pressures.""" + P_sol, T_sol = katz_2003(kind='solidus', X_h2o=30.0, Pmin=50.0, Pmax=100.0, n=50) + P_liq, T_liq = katz_2003(kind='liquidus', X_h2o=30.0, Pmin=50.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq) + assert np.all(T_sol < T_liq) + + +# ============================================================================= +# TEST: Melting Curve Functions - Lin (2024, Oxygen Fugacity) +# ============================================================================= + + +@pytest.mark.unit +class TestLin2024: + """Test Lin et al. (2024) oxygen-fugacity-dependent solidus.""" + + def test_solidus_default_fo2(self): + """Solidus with default fO2 = -4.""" + P, T = lin_2024(kind='solidus', fO2=-4.0, Pmin=50.0, Pmax=100.0, n=50) + assert len(P) == 50 + assert len(T) == 50 + assert np.all(T > 0.0) + + def test_fo2_effect(self): + """fO2 parameter affects melting temperature (Lin et al. 2024 parametrization). + + Lower (more reducing) fO2 values increase the solidus temperature; + higher (more oxidizing) fO2 values decrease it. + """ + P, T_reducing = lin_2024(kind='solidus', fO2=-5.0, Pmin=50.0, Pmax=100.0, n=50) + _, T_oxidizing = lin_2024(kind='solidus', fO2=-3.0, Pmin=50.0, Pmax=100.0, n=50) + + # Lower fO2 (more reducing) -> higher T + assert np.all(T_reducing > T_oxidizing), 'Lower fO2 should increase solidus T' + + def test_physical_constraint_solidus_less_than_liquidus(self): + """Solidus < liquidus with varying fO2 at valid pressures.""" + P_sol, T_sol = lin_2024(kind='solidus', fO2=-4.0, Pmin=50.0, Pmax=100.0, n=50) + P_liq, T_liq = lin_2024(kind='liquidus', fO2=-4.0, Pmin=50.0, Pmax=100.0, n=50) + + assert np.allclose(P_sol, P_liq) + assert np.all(T_sol < T_liq) + + +# ============================================================================= +# TEST: get_melting_curves Dispatcher +# ============================================================================= + + +@pytest.mark.unit +class TestGetMeltingCurves: + """Test melting curve dispatcher.""" + + @pytest.mark.parametrize( + 'model_name,Pmin,Pmax', + [ + ('andrault_2011', 0.0, 100.0), + ('monteux_2016', 0.0, 100.0), + ('wolf_bower_2018', 0.0, 100.0), + ('katz_2003', 0.0, 100.0), + ('fei_2021', 1.0, 100.0), # Requires Pmin >= 1 + ('belonoshko_2005', 0.0, 100.0), + ('fiquet_2010', 0.0, 100.0), + ('hirschmann_2000', 0.0, 5.0), # Low pressure only + ('stixrude_2014', 1.0, 100.0), # Requires Pmin >= 1 + ('lin_2024', 0.0, 100.0), + ], + ) + def test_all_supported_models(self, model_name, Pmin, Pmax): + """Dispatcher returns non-empty, shape-consistent outputs for all models.""" + P_sol, T_sol, P_liq, T_liq = get_melting_curves(model_name, Pmin=Pmin, Pmax=Pmax, n=50) + + assert len(P_sol) > 0, f'{model_name}: no pressure values' + assert len(T_sol) > 0, f'{model_name}: no temperature values' + assert len(P_liq) > 0, f'{model_name}: no pressure values for liquidus' + assert len(T_liq) > 0, f'{model_name}: no liquidus temperature values' + assert len(P_sol) == len(T_sol), f'{model_name}: solidus P/T length mismatch' + assert len(P_liq) == len(T_liq), f'{model_name}: liquidus P/T length mismatch' + assert np.all(np.isfinite(P_sol)), f'{model_name}: non-finite values in P_sol' + assert np.all(np.isfinite(T_sol)), f'{model_name}: non-finite values in T_sol' + assert np.all(np.isfinite(P_liq)), f'{model_name}: non-finite values in P_liq' + assert np.all(np.isfinite(T_liq)), f'{model_name}: non-finite values in T_liq' + + def test_unknown_model(self): + """Unknown model raises ValueError.""" + with pytest.raises(ValueError, match='Unknown model'): + get_melting_curves('unknown_model') + + def test_custom_pressure_range(self): + """Custom pressure range is applied correctly.""" + P_sol, T_sol, _, _ = get_melting_curves('wolf_bower_2018', Pmin=10.0, Pmax=50.0, n=50) + + # After truncation, pressure should be within requested range + if len(P_sol) > 0: + assert P_sol[0] >= 10.0 - 1.0, 'Lowest pressure should be near requested Pmin' + assert P_sol[-1] <= 50.0 + 1.0, 'Highest pressure should be near requested Pmax' + + +# ============================================================================= +# TEST: truncate_to_physical_interval Decorator +# ============================================================================= + + +@pytest.mark.unit +class TestTruncateToPhysicalInterval: + """Test physical interval filtering decorator.""" + + def test_decorator_filters_unphysical_regions(self): + """Truncation removes regions where solidus >= liquidus.""" + # Create a wrapped function + wrapped_wolf = truncate_to_physical_interval(wolf_bower_2018) + + P_sol, T_sol = wrapped_wolf(kind='solidus', Pmin=0.0, Pmax=100.0, n=100) + P_liq, T_liq = wrapped_wolf(kind='liquidus', Pmin=0.0, Pmax=100.0, n=100) + + # In the physical interval, solidus < liquidus + assert np.all(T_sol < T_liq) + + def test_decorator_preserves_main_block(self): + """Main physical block is largest contiguous block.""" + wrapped_wolf = truncate_to_physical_interval(wolf_bower_2018) + + # Should not raise + P_sol, T_sol = wrapped_wolf(kind='solidus', Pmin=0.0, Pmax=1000.0, n=200) + assert len(P_sol) > 0 + + +# ============================================================================= +# TEST: EOS Inversion Pipeline - invert_to_entropy_along_profile +# ============================================================================= + + +@pytest.mark.unit +class TestInvertToEntropyAlongProfile: + """Test EOS inversion to entropy curves with mocked T(S,P) interpolator.""" + + def test_basic_inversion(self, mock_eos_interpolator, mock_S_axes): + """Basic entropy inversion with mocked EOS. + + Note: With the simplified mock interpolator, inversion may not always + find valid S values. This test verifies the function runs without error + and returns the correct shape. + """ + S_axis_solid, _ = mock_S_axes + # Use mid-range pressure/temperature values + P = np.array([50.0, 60.0, 70.0]) # GPa + # For mock T(S,P) = 1000 + 10*S + 100*P, with S in [100,1000] and P in [0,100]: + # T ranges from 1000+1000+0 = 2000 to 1000+10000+10000 = 21000 + # Use realistic T values for the given P, S ranges + T = np.array([6500.0, 7000.0, 7500.0]) # K (plausible for mid-range S,P) + + S = invert_to_entropy_along_profile(P, T, S_axis_solid, mock_eos_interpolator) + + assert S.shape == P.shape + # May or may not find valid values depending on mock bounds, but should not crash + assert isinstance(S, np.ndarray) + + def test_inversion_returns_entropy_values(self, mock_eos_interpolator, mock_S_axes): + """Entropy values are within expected bounds when valid.""" + S_axis_solid, _ = mock_S_axes + P = np.array([50.0]) # Mid-range pressure + T = np.array([2000.0]) # Mid-range temperature + + S = invert_to_entropy_along_profile(P, T, S_axis_solid, mock_eos_interpolator) + + # When valid, S should be in the S_axis range (approximately) + valid = np.isfinite(S) + if np.any(valid): + assert np.all(S[valid] >= S_axis_solid.min() - 100) + assert np.all(S[valid] <= S_axis_solid.max() + 100) + + def test_out_of_bounds_returns_nan(self, mock_eos_interpolator, mock_S_axes): + """Out-of-bounds P,T return NaN.""" + S_axis_solid, _ = mock_S_axes + P = np.array([1000.0]) # Outside mock bounds [0, 100] + T = np.array([10000.0]) # Unphysical temperature + + S = invert_to_entropy_along_profile(P, T, S_axis_solid, mock_eos_interpolator) + + assert np.isnan(S[0]), 'Out-of-bounds should return NaN' + + def test_empty_pressure_array(self, mock_eos_interpolator, mock_S_axes): + """Empty pressure array returns empty entropy array.""" + S_axis_solid, _ = mock_S_axes + P = np.array([]) + T = np.array([]) + + S = invert_to_entropy_along_profile(P, T, S_axis_solid, mock_eos_interpolator) + + assert len(S) == 0 + + +# ============================================================================= +# TEST: build_common_entropy_grid +# ============================================================================= + + +@pytest.mark.unit +class TestBuildCommonEntropyGrid: + """Test resampling solidus/liquidus onto common pressure grid.""" + + def test_basic_resampling(self): + """Basic resampling onto common pressure grid.""" + P_sol = np.array([10.0, 20.0, 30.0, 40.0, 50.0]) + S_sol = np.array([100.0, 150.0, 200.0, 250.0, 300.0]) + + P_liq = np.array([10.0, 20.0, 30.0, 40.0, 50.0]) + S_liq = np.array([150.0, 200.0, 250.0, 300.0, 350.0]) + + P_common, S_sol_common, S_liq_common = build_common_entropy_grid( + P_sol, S_sol, P_liq, S_liq, n_common=5 + ) + + assert len(P_common) > 0 + assert len(S_sol_common) == len(P_common) + assert len(S_liq_common) == len(P_common) + + def test_common_grid_within_overlap(self): + """Common pressure grid is within overlap of solidus/liquidus ranges.""" + P_sol = np.array([5.0, 10.0, 15.0, 20.0, 25.0]) + S_sol = np.array([100.0, 150.0, 200.0, 250.0, 300.0]) + + P_liq = np.array([8.0, 12.0, 16.0, 20.0, 24.0]) + S_liq = np.array([150.0, 200.0, 250.0, 300.0, 350.0]) + + P_common, _, _ = build_common_entropy_grid(P_sol, S_sol, P_liq, S_liq) + + if len(P_common) > 0: + P_min_overlap = max(P_sol.min(), P_liq.min()) + P_max_overlap = min(P_sol.max(), P_liq.max()) + assert np.all(P_common >= P_min_overlap - 1.0) + assert np.all(P_common <= P_max_overlap + 1.0) + + def test_nans_filtered_out(self): + """NaN values are filtered out before resampling.""" + P_sol = np.array([10.0, 20.0, 30.0, 40.0, 50.0]) + S_sol = np.array([100.0, np.nan, 200.0, 250.0, 300.0]) + + P_liq = np.array([10.0, 20.0, 30.0, 40.0, 50.0]) + S_liq = np.array([150.0, 200.0, 250.0, np.nan, 350.0]) + + P_common, S_sol_common, S_liq_common = build_common_entropy_grid( + P_sol, S_sol, P_liq, S_liq + ) + + assert np.all(np.isfinite(P_common)) + assert np.all(np.isfinite(S_sol_common)) + assert np.all(np.isfinite(S_liq_common)) + + def test_no_overlap_returns_empty(self): + """Non-overlapping pressure ranges return empty arrays.""" + P_sol = np.array([10.0, 20.0, 30.0]) + S_sol = np.array([100.0, 150.0, 200.0]) + + P_liq = np.array([50.0, 60.0, 70.0]) + S_liq = np.array([150.0, 200.0, 250.0]) + + P_common, S_sol_common, S_liq_common = build_common_entropy_grid( + P_sol, S_sol, P_liq, S_liq + ) + + assert len(P_common) == 0 + assert len(S_sol_common) == 0 + assert len(S_liq_common) == 0 + + def test_custom_n_common(self): + """Custom n_common parameter sets output grid size.""" + P_sol = np.linspace(10.0, 50.0, 100) + S_sol = np.linspace(100.0, 500.0, 100) + + P_liq = np.linspace(10.0, 50.0, 100) + S_liq = np.linspace(150.0, 550.0, 100) + + P_common, _, _ = build_common_entropy_grid(P_sol, S_sol, P_liq, S_liq, n_common=25) + + if len(P_common) > 0: + assert len(P_common) <= 25 + + +# ============================================================================= +# TEST: validate_entropy_export_arrays +# ============================================================================= + + +@pytest.mark.unit +class TestValidateEntropyExportArrays: + """Test validation of entropy export arrays.""" + + def test_valid_arrays_passes(self): + """Valid arrays pass validation silently.""" + P = np.array([10.0, 20.0, 30.0]) + S_sol = np.array([100.0, 150.0, 200.0]) + S_liq = np.array([150.0, 200.0, 250.0]) + + # Should not raise + validate_entropy_export_arrays(P, S_sol, S_liq, 'test_model') + + def test_empty_array_raises(self): + """Empty array raises ValueError.""" + P = np.array([]) + S_sol = np.array([]) + S_liq = np.array([]) + + with pytest.raises(ValueError, match='could not build|empty'): + validate_entropy_export_arrays(P, S_sol, S_liq, 'test_model') + + def test_mismatched_lengths_raises(self): + """Mismatched array lengths raise ValueError.""" + P = np.array([10.0, 20.0, 30.0]) + S_sol = np.array([100.0, 150.0]) # Wrong length + S_liq = np.array([150.0, 200.0, 250.0]) + + with pytest.raises(ValueError, match='inconsistent'): + validate_entropy_export_arrays(P, S_sol, S_liq, 'test_model') + + def test_nans_in_arrays_raises(self): + """NaN values raise ValueError.""" + P = np.array([10.0, 20.0, np.nan]) + S_sol = np.array([100.0, 150.0, 200.0]) + S_liq = np.array([150.0, 200.0, 250.0]) + + with pytest.raises(ValueError, match='non-finite'): + validate_entropy_export_arrays(P, S_sol, S_liq, 'test_model') + + def test_infs_in_arrays_raises(self): + """Infinite values raise ValueError.""" + P = np.array([10.0, 20.0, 30.0]) + S_sol = np.array([100.0, np.inf, 200.0]) + S_liq = np.array([150.0, 200.0, 250.0]) + + with pytest.raises(ValueError, match='non-finite'): + validate_entropy_export_arrays(P, S_sol, S_liq, 'test_model') + + +# ============================================================================= +# TEST: make_entropy_header +# ============================================================================= + + +@pytest.mark.unit +class TestMakeEntropyHeader: + """Test entropy table header generation.""" + + def test_default_header(self): + """Default header includes n_rows and scaling factors.""" + header = make_entropy_header(n_rows=100) + + assert '100' in header + assert '1e+09' in header or '1000000000' in header + assert 'Pressure' in header + assert 'Entropy' in header + assert 'scaling factor' in header + + def test_custom_scaling_factors(self): + """Custom scaling factors appear in header.""" + scale_p = 5e8 + scale_s = 1e6 + + header = make_entropy_header(n_rows=50, scale_p_out=scale_p, scale_s_out=scale_s) + + assert '50' in header + assert str(scale_p) in header or '5e+08' in header + assert str(scale_s) in header or '1e+06' in header + + def test_multiline_header(self): + """Header is multiline with proper structure.""" + header = make_entropy_header(n_rows=100) + lines = header.split('\n') + + assert len(lines) >= 5 + assert all(line.startswith('#') for line in lines if line.strip()) + + +# ============================================================================= +# INTEGRATION SMOKE TESTS (Basic End-to-End) +# ============================================================================= + + +@pytest.mark.unit +class TestMeltingCurveSmoke: + """Smoke tests: basic physical sanity checks across all models.""" + + @pytest.mark.parametrize( + 'model_name,Pmin,Pmax', + [ + ('andrault_2011', 0.0, 100.0), + ('monteux_2016', 0.0, 100.0), + ('wolf_bower_2018', 0.0, 100.0), + ('katz_2003', 0.0, 100.0), + ('fei_2021', 1.0, 100.0), + ('belonoshko_2005', 0.0, 100.0), + ('fiquet_2010', 0.0, 100.0), + ('hirschmann_2000', 0.0, 10.0), + ('stixrude_2014', 1.0, 100.0), + ('lin_2024', 0.0, 100.0), + ], + ) + def test_model_physical_constraints(self, model_name, Pmin, Pmax): + """Each model satisfies basic physical constraints.""" + P_sol, T_sol, P_liq, T_liq = get_melting_curves(model_name, Pmin=Pmin, Pmax=Pmax, n=50) + + # Soundness checks + assert len(P_sol) > 0, f'{model_name}: no output' + assert np.all(np.isfinite(P_sol)), f'{model_name}: P_sol contains NaNs/infs' + assert np.all(np.isfinite(T_sol)), f'{model_name}: T_sol contains NaNs/infs' + assert np.all(np.isfinite(P_liq)), f'{model_name}: P_liq contains NaNs/infs' + assert np.all(np.isfinite(T_liq)), f'{model_name}: T_liq contains NaNs/infs' + + # Physics constraints + assert np.all(T_sol > 0.0), f'{model_name}: T_sol <= 0' + assert np.all(T_liq > 0.0), f'{model_name}: T_liq <= 0' + assert np.all(T_sol < T_liq), f'{model_name}: solidus >= liquidus' + + # Monotonicity + assert np.all(np.diff(T_sol) >= -1e-10), f'{model_name}: T_sol not monotonic' + assert np.all(np.diff(T_liq) >= -1e-10), f'{model_name}: T_liq not monotonic' diff --git a/tests/utils/test_data.py b/tests/utils/test_data.py index ee72c82f0..5fcbb5610 100644 --- a/tests/utils/test_data.py +++ b/tests/utils/test_data.py @@ -2355,3 +2355,30 @@ def test_download_eos_static_delegates(mock_seager): download_eos_static() mock_seager.assert_called_once() + + +@pytest.mark.unit +@patch('proteus.utils.data.download') +@patch('proteus.utils.data.GetFWLData') +def test_download_melting_curves_skips_when_canonical_files_exist( + mock_getfwl, mock_download, tmp_path +): + """download_melting_curves should return early when all canonical files already exist.""" + from unittest.mock import MagicMock + + from proteus.utils.data import download_melting_curves + + mock_getfwl.return_value = tmp_path + + mc_dir = tmp_path / 'interior_lookup_tables' / 'Melting_curves' / 'Wolf_Bower+2018' + mc_dir.mkdir(parents=True, exist_ok=True) + + for name in ['solidus_P-T.dat', 'liquidus_P-T.dat', 'solidus_P-S.dat', 'liquidus_P-S.dat']: + (mc_dir / name).write_text('dummy\n') + + mock_config = MagicMock() + mock_config.interior.melting_dir = 'Wolf_Bower+2018' + + download_melting_curves(mock_config, clean=False) + + mock_download.assert_not_called() diff --git a/tools/solidus_func.py b/tools/solidus_func.py new file mode 100644 index 000000000..46dbf4c0b --- /dev/null +++ b/tools/solidus_func.py @@ -0,0 +1,1302 @@ +""" +Generate melting curves in pressure–temperature (P–T) and pressure–entropy (P–S) +space for several literature parametrizations of peridotite / silicate melting. + +This version is designed to work with the legacy EOS lookup tables: + + - temperature_solid.dat + - temperature_melt.dat + +These tables provide temperature as a function of entropy and pressure, +T(S, P), on structured grids. The script therefore: + +1. Builds solidus and liquidus curves in P–T space from literature fits. +2. Converts those curves into P–S space by inverting the EOS relation T(S, P). +3. Resamples the solidus and liquidus entropy curves onto a common pressure grid. +4. Saves both the P–T and P–S versions to disk. + +Usage +----- +Export one model with a dedicated shortcut: + + python solidus_func.py --katz2003 --X-h2o 30 + python solidus_func.py --lin2024 --fO2 -4 + +Export one model by explicit internal name: + + python solidus_func.py --model wolf_bower_2018 + +Export all supported models: + + python solidus_func.py --all + +Notes +----- +- Pressure is handled internally in GPa for the melting-curve parametrizations, + unless a given published fit is explicitly defined in Pa and converted. +- The EOS tables are assumed to have the SPIDER-style format with scaling factors + in the header. +- The exported entropy files use the same scaled SPIDER-like header so they can + be re-used by downstream tools expecting that format. + +References included here +------------------------ +- Andrault et al. (2011), DOI: https://doi.org/10.1016/j.epsl.2011.02.006 +- Monteux et al. (2016), DOI: https://doi.org/10.1016/j.epsl.2016.05.010 +- Lin et al. (2024), DOI: https://doi.org/10.1038/s41561-024-01495-1 +- Hirschmann (2000), DOI: https://doi.org/10.1029/2000GC000070 +- Katz et al. (2003), DOI: https://doi.org/10.1029/2002GC000433 +- Stixrude (2014), DOI: https://doi.org/10.1098/rsta.2013.0076 +- Fei et al. (2021), DOI: https://doi.org/10.1038/s41467-021-21170-y +- Belonoshko et al. (2005), DOI: https://doi.org/10.1103/PhysRevLett.94.195701 +- Fiquet et al. (2010), DOI: https://doi.org/10.1126/science.1192448 +""" + +from __future__ import annotations + +import argparse +import logging +import os +from pathlib import Path + +import numpy as np +import platformdirs +from scipy.interpolate import RegularGridInterpolator, interp1d + +from proteus.utils.coupler import get_proteus_directories + +log = logging.getLogger(__name__) +# ============================================================================= +# DEFAULT OUTPUT LOCATION +# ============================================================================= + +FWL_DATA_DIR = Path(os.environ.get('FWL_DATA', platformdirs.user_data_dir('fwl_data'))) +MELTING_DIR = FWL_DATA_DIR / 'interior_lookup_tables' / 'Melting_curves' + + +# ============================================================================= +# GENERAL HELPERS +# ============================================================================= + + +def make_pressure_grid(Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500) -> np.ndarray: + r""" + Create a uniformly sampled pressure grid in GPa. + + Parameters + ---------- + Pmin : float, optional + Minimum pressure in GPa. + Pmax : float, optional + Maximum pressure in GPa. + n : int, optional + Number of grid points. + + Returns + ------- + np.ndarray + One-dimensional pressure array in GPa. + """ + return np.linspace(Pmin, Pmax, n) + + +def solidus_from_liquidus_stixrude(T_liq: np.ndarray) -> np.ndarray: + r""" + Estimate the solidus from a liquidus using the Stixrude ratio. + """ + return T_liq / (1.0 - np.log(0.79)) + + +def liquidus_from_solidus_stixrude(T_sol: np.ndarray) -> np.ndarray: + r""" + Estimate the liquidus from a solidus using the inverse Stixrude ratio. + """ + return T_sol * (1.0 - np.log(0.79)) + + +def _fmt_range(arr: np.ndarray, fmt: str = '.3f') -> str: + """ + Format the finite minimum and maximum of an array as a string. + """ + arr = np.asarray(arr, dtype=float) + mask = np.isfinite(arr) + + if not np.any(mask): + return '[nan, nan]' + + amin = np.min(arr[mask]) + amax = np.max(arr[mask]) + return f'[{amin:{fmt}}, {amax:{fmt}}]' + + +DISPLAY_NAMES = { + 'andrault_2011': 'Andrault et al. (2011)', + 'monteux_2016': 'Monteux et al. (2016)', + 'wolf_bower_2018': 'Wolf & Bower (2018)', + 'katz_2003': 'Katz et al. (2003)', + 'fei_2021': 'Fei et al. (2021)', + 'belonoshko_2005': 'Belonoshko et al. (2005)', + 'fiquet_2010': 'Fiquet et al. (2010)', + 'hirschmann_2000': 'Hirschmann (2000)', + 'stixrude_2014': 'Stixrude (2014)', + 'lin_2024': 'Lin et al. (2024)', +} + + +def print_model_summary( + model_name: str, + P_sol: np.ndarray, + T_sol: np.ndarray, + P_liq: np.ndarray, + T_liq: np.ndarray, + P_common: np.ndarray, + S_sol_common: np.ndarray, + S_liq_common: np.ndarray, +): + """ + Log a compact summary of the exported P-T and P-S ranges for one model. + """ + label = DISPLAY_NAMES.get(model_name, model_name) + + log.info(f"{label}:") + log.info( + f" P-T solidus : P_range = {_fmt_range(P_sol)} GPa, " + f"T_range = {_fmt_range(T_sol)} K" + ) + log.info( + f" P-T liquidus : P_range = {_fmt_range(P_liq)} GPa, " + f"T_range = {_fmt_range(T_liq)} K" + ) + log.info( + f" P-S solidus : P_range = {_fmt_range(P_common)} GPa, " + f"S_range = {_fmt_range(S_sol_common)} J kg^-1 K^-1" + ) + log.info( + f" P-S liquidus : P_range = {_fmt_range(P_common)} GPa, " + f"S_range = {_fmt_range(S_liq_common)} J kg^-1 K^-1" + ) + log.info("") # blank line for spacing + + +# ============================================================================= +# LITERATURE MELTING CURVES +# ============================================================================= + + +def andrault_2011(kind: str = 'solidus', Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500): + r""" + Melting curve from Andrault et al. (2011). + """ + P = make_pressure_grid(Pmin, Pmax, n) + P_pa = P * 1e9 + + if kind == 'solidus': + T0, a, c = 2045, 92e9, 1.3 + elif kind == 'liquidus': + T0, a, c = 1940, 29e9, 1.9 + else: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + T = T0 * ((P_pa / a) + 1.0) ** (1.0 / c) + return P, T + + +def fei_2021(kind: str = 'liquidus', Pmin: float = 1.0, Pmax: float = 1000.0, n: int = 500): + r""" + Melting curve based on Fei et al. (2021). + """ + P = make_pressure_grid(Pmin, Pmax, n) + T_liq = 6000.0 * (P / 140.0) ** 0.26 + + if kind == 'liquidus': + T = T_liq + elif kind == 'solidus': + T = solidus_from_liquidus_stixrude(T_liq) + else: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + return P, T + + +def belonoshko_2005( + kind: str = 'liquidus', Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500 +): + r""" + Melting curve based on Belonoshko et al. (2005). + """ + P = make_pressure_grid(Pmin, Pmax, n) + T_liq = 1831.0 * (1.0 + P / 4.6) ** 0.33 + + if kind == 'liquidus': + T = T_liq + elif kind == 'solidus': + T = solidus_from_liquidus_stixrude(T_liq) + else: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + return P, T + + +def fiquet_2010(kind: str = 'liquidus', Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500): + r""" + Melting curve based on Fiquet et al. (2010). + """ + P = make_pressure_grid(Pmin, Pmax, n) + T_liq = np.zeros_like(P, dtype=float) + + low = P <= 20.0 + high = P > 20.0 + + # Original high-pressure constant is given in Pa in the paper. + # Here pressure is in GPa, so 4.054e6 Pa -> 0.004056 GPa. + T_liq[low] = 1982.1 * ((P[low] / 6.594) + 1.0) ** (1.0 / 5.374) + T_liq[high] = 78.74 * ((P[high] / 0.004056) + 1.0) ** (1.0 / 2.44) + + if kind == 'liquidus': + T = T_liq + elif kind == 'solidus': + T = solidus_from_liquidus_stixrude(T_liq) + else: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + return P, T + + +def monteux_2016(kind: str = 'solidus', Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500): + r""" + Melting curve from Monteux et al. (2016). + """ + P = make_pressure_grid(Pmin, Pmax, n) + P_pa = P * 1e9 + + params = { + 'solidus': { + 'low': {'T0': 1661.2, 'a': 1.336e9, 'c': 7.437}, + 'high': {'T0': 2081.8, 'a': 101.69e9, 'c': 1.226}, + }, + 'liquidus': { + 'low': {'T0': 1982.1, 'a': 6.594e9, 'c': 5.374}, + 'high': {'T0': 2006.8, 'a': 34.65e9, 'c': 1.844}, + }, + } + + if kind not in params: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + p = params[kind] + T = np.zeros_like(P_pa, dtype=float) + + mask_low = P_pa <= 20e9 + mask_high = P_pa > 20e9 + + T[mask_low] = p['low']['T0'] * ((P_pa[mask_low] / p['low']['a']) + 1.0) ** ( + 1.0 / p['low']['c'] + ) + T[mask_high] = p['high']['T0'] * ((P_pa[mask_high] / p['high']['a']) + 1.0) ** ( + 1.0 / p['high']['c'] + ) + + return P, T + + +def hirschmann_2000(kind: str = 'solidus', Pmin: float = 0.0, Pmax: float = 10.0, n: int = 500): + r""" + Melting curve from Hirschmann (2000). + """ + P = make_pressure_grid(Pmin, Pmax, n) + + coeffs = { + 'solidus': (1085.7, 132.9, -5.1), + 'liquidus': (1475.0, 80.0, -3.2), + } + + if kind not in coeffs: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + A1, A2, A3 = coeffs[kind] + T_c = A1 + A2 * P + A3 * P**2 + T = T_c + 273.15 + + return P, T + + +def stixrude_2014( + kind: str = 'liquidus', Pmin: float = 1.0, Pmax: float = 1000.0, n: int = 500 +): + r""" + Melting curve based on Stixrude (2014). + """ + P = make_pressure_grid(Pmin, Pmax, n) + T_liq = 5400.0 * (P / 140.0) ** 0.480 + + if kind == 'liquidus': + T = T_liq + elif kind == 'solidus': + T = solidus_from_liquidus_stixrude(T_liq) + else: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + return P, T + + +def wolf_bower_2018( + kind: str = 'solidus', Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500 +): + r""" + Piecewise melting curve based on Wolf & Bower (2018) style fits. + """ + P = make_pressure_grid(Pmin, Pmax, n) + + params = { + 'solidus': ( + 7.696777581585296, + 870.4767697319186, + 101.52655163737373, + 15.959022187236807, + 3.090844734784906, + 1417.4258954709148, + ), + 'liquidus': ( + 8.864665249317456, + 408.58442302949794, + 46.288444869816615, + 17.549174419770257, + 3.679647802112376, + 2019.967799687511, + ), + } + + if kind not in params: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + cp1, cp2, s1, s2, s3, intercept = params[kind] + + c1 = intercept + c2 = c1 + (s1 - s2) * cp1 + c3 = c2 + (s2 - s3) * cp2 + + T = np.zeros_like(P, dtype=float) + + m1 = P < cp1 + m2 = (P >= cp1) & (P < cp2) + m3 = P >= cp2 + + T[m1] = s1 * P[m1] + c1 + T[m2] = s2 * P[m2] + c2 + T[m3] = s3 * P[m3] + c3 + + return P, T + + +def katz_2003( + kind: str = 'solidus', + X_h2o: float = 30.0, + Pmin: float = 0.0, + Pmax: float = 1000.0, + n: int = 500, +): + r""" + Hydrous melting-curve correction following Katz et al. (2003). + """ + gamma = 0.75 + K = 43.0 + + if kind not in {'solidus', 'liquidus'}: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + P, T_anhydrous = wolf_bower_2018(kind=kind, Pmin=Pmin, Pmax=Pmax, n=n) + delta_T = K * X_h2o**gamma + T = T_anhydrous - delta_T + + return P, T + + +def lin_2024( + kind: str = 'solidus', + fO2: float = -4.0, + Pmin: float = 0.0, + Pmax: float = 1000.0, + n: int = 500, +): + r""" + Oxygen-fugacity-dependent solidus following Lin et al. (2024). + """ + P, T_anhydrous = wolf_bower_2018(kind='solidus', Pmin=Pmin, Pmax=Pmax, n=n) + + delta_T = (340.0 / 3.2) * (2.0 - fO2) + T_sol = T_anhydrous + delta_T + + if kind == 'solidus': + T = T_sol + elif kind == 'liquidus': + T = liquidus_from_solidus_stixrude(T_sol) + else: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + return P, T + + +# ============================================================================= +# PHYSICAL-INTERVAL FILTER +# ============================================================================= + + +def truncate_to_physical_interval(func): + r""" + Wrap a melting-curve function so only the main interval with + T_sol < T_liq is retained. + """ + + def wrapped(kind='solidus', Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): + P_sol, T_sol = func(kind='solidus', Pmin=Pmin, Pmax=Pmax, n=n, **kwargs) + P_liq, T_liq = func(kind='liquidus', Pmin=Pmin, Pmax=Pmax, n=n, **kwargs) + + good = T_sol < T_liq + idx = np.where(good)[0] + + if len(idx) == 0: + raise ValueError(f'{func.__name__}: no physical interval where solidus < liquidus') + + splits = np.where(np.diff(idx) > 1)[0] + 1 + blocks = np.split(idx, splits) + main_block = max(blocks, key=len) + + if kind == 'solidus': + return P_sol[main_block], T_sol[main_block] + if kind == 'liquidus': + return P_liq[main_block], T_liq[main_block] + raise ValueError("kind must be 'solidus' or 'liquidus'") + + return wrapped + + + + +# ============================================================================= +# MODEL DISPATCHER +# ============================================================================= + + +SUPPORTED_MODELS = list(DISPLAY_NAMES) + +def get_melting_curves( + model_name: str, Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 2000, **kwargs +): + r""" + Return solidus and liquidus curves for a given model. + """ + models = { + 'andrault_2011': truncate_to_physical_interval(andrault_2011), + 'monteux_2016': truncate_to_physical_interval(monteux_2016), + 'wolf_bower_2018': truncate_to_physical_interval(wolf_bower_2018), + 'katz_2003': truncate_to_physical_interval(katz_2003), + 'fei_2021': fei_2021, + 'belonoshko_2005': belonoshko_2005, + 'fiquet_2010': fiquet_2010, + 'hirschmann_2000': hirschmann_2000, + 'stixrude_2014': stixrude_2014, + 'lin_2024': lin_2024, + } + + if model_name not in models: + raise ValueError(f'Unknown model: {model_name}') + + func = models[model_name] + P_sol, T_sol = func(kind='solidus', Pmin=Pmin, Pmax=Pmax, n=n, **kwargs) + P_liq, T_liq = func(kind='liquidus', Pmin=Pmin, Pmax=Pmax, n=n, **kwargs) + + return P_sol, T_sol, P_liq, T_liq + + +# ============================================================================= +# OUTPUT HELPERS +# ============================================================================= + + +def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): + r""" + Save a pressure-temperature table to disk. + """ + data = np.column_stack([P_gpa, T_k]) + np.savetxt(path, data, fmt='%.18e %.18e', header='pressure temperature', comments='#') + + + +def make_entropy_header( + n_rows: int, + scale_p_out: float = 1_000_000_000.0, + scale_s_out: float = 4_824_266.84604467, +) -> str: + """ + Build the SPIDER-style header for a pressure-entropy table. + + Parameters + ---------- + n_rows : int + Number of data rows in the file. + scale_p_out : float, optional + Output pressure scaling factor written to the header. + scale_s_out : float, optional + Output entropy scaling factor written to the header. + + Returns + ------- + str + Multiline header string. + """ + return '\n'.join( + [ + f'# 5 {n_rows}', + '# Pressure, Entropy', + '# column * scaling factor should be SI units', + '# scaling factors (constant) for each column given on line below', + f'# {scale_p_out} {scale_s_out}', + ] + ) + + +def get_default_spider_dir() -> Path: + dirs = get_proteus_directories() + return Path(dirs["spider"]) + + +def validate_entropy_export_arrays( + P_common: np.ndarray, + S_sol_common: np.ndarray, + S_liq_common: np.ndarray, + model_name: str, +): + """ + Validate the common P-S arrays before writing them to disk. + + Parameters + ---------- + P_common : np.ndarray + Common pressure grid in GPa. + S_sol_common : np.ndarray + Solidus entropy values on the common pressure grid. + S_liq_common : np.ndarray + Liquidus entropy values on the common pressure grid. + model_name : str + Name of the melting model, used in error messages. + + Raises + ------ + ValueError + If the entropy export arrays are empty, mismatched, or non-finite. + """ + if len(P_common) == 0 or len(S_sol_common) == 0 or len(S_liq_common) == 0: + raise ValueError( + f'{model_name}: could not build a valid common P-S grid. ' + 'EOS inversion may have failed or the solidus/liquidus entropy ' + 'ranges may not overlap.' + ) + + if not (len(P_common) == len(S_sol_common) == len(S_liq_common)): + raise ValueError( + f'{model_name}: inconsistent P-S array lengths: ' + f'len(P_common)={len(P_common)}, ' + f'len(S_sol_common)={len(S_sol_common)}, ' + f'len(S_liq_common)={len(S_liq_common)}' + ) + + if not ( + np.all(np.isfinite(P_common)) + and np.all(np.isfinite(S_sol_common)) + and np.all(np.isfinite(S_liq_common)) + ): + raise ValueError(f'{model_name}: non-finite values found in common P-S export arrays.') + + +def resolve_eos_paths(spider_dir: Path | str | None = None) -> tuple[Path, Path]: + """ + Resolve and validate the default solid and liquid EOS lookup table paths. + + Parameters + ---------- + spider_dir : Path | str | None, optional + Root directory of the SPIDER repository. If None, a path relative to + this module is used. + + Returns + ------- + tuple[Path, Path] + Paths to the solid and liquid EOS tables. + + Raises + ------ + FileNotFoundError + If one or both EOS files are missing. + """ + if spider_dir is None: + spider_dir = get_default_spider_dir() + + spider_dir = Path(spider_dir).resolve() + + eos_solid_path = ( + spider_dir / 'lookup_data' / '1TPa-dK09-elec-free' / 'temperature_solid.dat' + ) + eos_liquid_path = ( + spider_dir / 'lookup_data' / '1TPa-dK09-elec-free' / 'temperature_melt.dat' + ) + + if not eos_solid_path.exists(): + raise FileNotFoundError(f'Missing EOS file: {eos_solid_path}') + + if not eos_liquid_path.exists(): + raise FileNotFoundError(f'Missing EOS file: {eos_liquid_path}') + + return eos_solid_path, eos_liquid_path + + +def load_eos_T_of_SP( + eos_path: Path, + nS: int, + scale_S_axis: float, + *, + nP: int = 2020, + skip_header: int = 5, + scale_p_eos: float = 1e9, + scale_t_eos: float = 1.0, +): + r""" + Load an EOS lookup table and build an interpolator for T(S, P). + + Parameters + ---------- + eos_path : Path + Path to the EOS lookup table. + nS : int + Number of entropy points in the EOS table. + scale_S_axis : float + Entropy scaling factor applied to the EOS entropy axis. + nP : int, optional + Number of pressure points in the EOS table. + skip_header : int, optional + Number of header lines to skip when reading the table. + scale_p_eos : float, optional + Pressure scaling factor used in the EOS table. + scale_t_eos : float, optional + Temperature scaling factor used in the EOS table. + """ + raw = np.genfromtxt(eos_path, skip_header=skip_header) + resh = raw.reshape((nS, nP, 3)) + + P_axis_GPa = resh[0, :, 0] * scale_p_eos / 1e9 + S_axis = resh[:, 0, 1] * scale_S_axis + T_grid = resh[:, :, 2] * scale_t_eos + + T_interp = RegularGridInterpolator( + (S_axis, P_axis_GPa), + T_grid, + bounds_error=False, + fill_value=np.nan, + ) + return S_axis, P_axis_GPa, T_interp + +def load_default_eos_interpolators( + spider_dir: Path | str | None = None, + *, + nP: int = 2020, + nS_solid: int = 125, + nS_liquid: int = 95, + skip_header: int = 5, + scale_p_eos: float = 1e9, + scale_t_eos: float = 1.0, + scale_s_solid_eos: float = 4.82426684604467e6, + scale_s_liquid_eos: float = 4.805046659407042e6, +) -> tuple[np.ndarray, RegularGridInterpolator, np.ndarray, RegularGridInterpolator]: + """ + Load the default solid and liquid EOS interpolators. + + Parameters + ---------- + spider_dir : Path | str | None, optional + Root directory of the SPIDER repository. If None, a path relative to + this module is used. + nP : int, optional + Number of pressure points in each EOS table. + nS_solid : int, optional + Number of entropy points in the solid EOS table. + nS_liquid : int, optional + Number of entropy points in the liquid EOS table. + skip_header : int, optional + Number of header lines to skip in the EOS tables. + scale_p_eos : float, optional + Pressure scaling factor used in the EOS tables. + scale_t_eos : float, optional + Temperature scaling factor used in the EOS tables. + scale_s_solid_eos : float, optional + Entropy scaling factor for the solid EOS table. + scale_s_liquid_eos : float, optional + Entropy scaling factor for the liquid EOS table. + + Returns + ------- + tuple + A tuple containing: + - solid entropy axis + - solid T(S, P) interpolator + - liquid entropy axis + - liquid T(S, P) interpolator + """ + eos_solid_path, eos_liquid_path = resolve_eos_paths(spider_dir=spider_dir) + + S_axis_solid, _, T_of_SP_solid = load_eos_T_of_SP( + eos_solid_path, + nS_solid, + scale_s_solid_eos, + nP=nP, + skip_header=skip_header, + scale_p_eos=scale_p_eos, + scale_t_eos=scale_t_eos, + ) + S_axis_liquid, _, T_of_SP_liquid = load_eos_T_of_SP( + eos_liquid_path, + nS_liquid, + scale_s_liquid_eos, + nP=nP, + skip_header=skip_header, + scale_p_eos=scale_p_eos, + scale_t_eos=scale_t_eos, + ) + + return S_axis_solid, T_of_SP_solid, S_axis_liquid, T_of_SP_liquid + +def invert_to_entropy_along_profile( + P_gpa: np.ndarray, T_k: np.ndarray, S_axis: np.ndarray, T_of_SP +): + r""" + Convert a P-T curve into a P-S curve by inverting T(S, P). + """ + S_out = np.full_like(T_k, np.nan, dtype=float) + + for i, (P_i_gpa, T_i) in enumerate(zip(P_gpa, T_k)): + P_col = np.full_like(S_axis, P_i_gpa) + T_vals = T_of_SP(np.column_stack([S_axis, P_col])) + + valid = np.isfinite(T_vals) + if np.count_nonzero(valid) < 2: + continue + + Tv = T_vals[valid] + Sv = S_axis[valid] + + order = np.argsort(Tv) + T_sorted = Tv[order] + S_sorted = Sv[order] + + T_unique, idx_unique = np.unique(T_sorted, return_index=True) + S_unique = S_sorted[idx_unique] + + if len(T_unique) < 2: + continue + + if T_i < T_unique[0] or T_i > T_unique[-1]: + continue + + try: + f = interp1d(T_unique, S_unique, kind='linear', assume_sorted=True) + S_out[i] = float(f(T_i)) + except Exception: + try: + f = interp1d(T_unique, S_unique, kind='nearest', assume_sorted=True) + S_out[i] = float(f(T_i)) + except Exception: + continue + + return S_out + + +def build_common_entropy_grid( + P_sol: np.ndarray, + S_sol: np.ndarray, + P_liq: np.ndarray, + S_liq: np.ndarray, + n_common: int | None = None, +): + r""" + Resample solidus and liquidus entropy curves onto a shared pressure grid. + """ + mask_sol = np.isfinite(S_sol) + mask_liq = np.isfinite(S_liq) + + P_sol_v = np.asarray(P_sol[mask_sol], dtype=float) + S_sol_v = np.asarray(S_sol[mask_sol], dtype=float) + P_liq_v = np.asarray(P_liq[mask_liq], dtype=float) + S_liq_v = np.asarray(S_liq[mask_liq], dtype=float) + + if len(P_sol_v) < 2 or len(P_liq_v) < 2: + return np.array([]), np.array([]), np.array([]) + + Pmin_common = max(np.min(P_sol_v), np.min(P_liq_v)) + Pmax_common = min(np.max(P_sol_v), np.max(P_liq_v)) + + if ( + not np.isfinite(Pmin_common) + or not np.isfinite(Pmax_common) + or Pmax_common <= Pmin_common + ): + return np.array([]), np.array([]), np.array([]) + + if n_common is None: + n_common = min(len(P_sol_v), len(P_liq_v)) + + order_sol = np.argsort(P_sol_v) + order_liq = np.argsort(P_liq_v) + + P_sol_s = P_sol_v[order_sol] + S_sol_s = S_sol_v[order_sol] + P_liq_s = P_liq_v[order_liq] + S_liq_s = S_liq_v[order_liq] + + P_sol_u, idx_sol = np.unique(P_sol_s, return_index=True) + S_sol_u = S_sol_s[idx_sol] + + P_liq_u, idx_liq = np.unique(P_liq_s, return_index=True) + S_liq_u = S_liq_s[idx_liq] + + if len(P_sol_u) < 2 or len(P_liq_u) < 2: + return np.array([]), np.array([]), np.array([]) + + P_common = np.linspace(Pmin_common, Pmax_common, n_common) + + f_sol = interp1d( + P_sol_u, + S_sol_u, + kind='linear', + bounds_error=False, + fill_value=np.nan, + assume_sorted=True, + ) + f_liq = interp1d( + P_liq_u, + S_liq_u, + kind='linear', + bounds_error=False, + fill_value=np.nan, + assume_sorted=True, + ) + + S_sol_common = f_sol(P_common) + S_liq_common = f_liq(P_common) + + mask = np.isfinite(S_sol_common) & np.isfinite(S_liq_common) + return P_common[mask], S_sol_common[mask], S_liq_common[mask] + + +def save_entropy_table_with_header( + path: Path, + P_gpa: np.ndarray, + S_jpk: np.ndarray, + *, + scale_p_out: float = 1_000_000_000.0, + scale_s_out: float = 4_824_266.84604467, +): + r""" + Save a pressure-entropy table in SPIDER-style scaled format. + + Parameters + ---------- + path : Path + Output file path. + P_gpa : np.ndarray + Pressure array in GPa. + S_jpk : np.ndarray + Entropy array in J kg^-1 K^-1. + scale_p_out : float, optional + Pressure scaling factor used in the written file. + scale_s_out : float, optional + Entropy scaling factor used in the written file. + """ + P_pa = P_gpa * 1e9 + data = np.column_stack([P_pa / scale_p_out, S_jpk / scale_s_out]) + header = make_entropy_header( + len(P_gpa), + scale_p_out=scale_p_out, + scale_s_out=scale_s_out, + ) + np.savetxt(path, data, fmt='%.18e %.18e', header=header, comments='') + +# ============================================================================= +# MAIN EXPORTER +# ============================================================================= + + +def export_model_curves( + model_name: str, + out_root: Path | str = MELTING_DIR, + Pmin: float = 0.0, + Pmax: float = 1000.0, + n: int = 2000, + spider_dir: Path | str | None = None, + **kwargs, +): + r""" + Export one melting model in both P-T and P-S space. + """ + out_dir = Path(out_root) / model_name + out_dir.mkdir(parents=True, exist_ok=True) + + P_sol, T_sol, P_liq, T_liq = get_melting_curves( + model_name, Pmin=Pmin, Pmax=Pmax, n=n, **kwargs + ) + + save_PT_table(out_dir / 'solidus_P-T.dat', P_sol, T_sol) + save_PT_table(out_dir / 'liquidus_P-T.dat', P_liq, T_liq) + + S_axis_solid, T_of_SP_solid, S_axis_liquid, T_of_SP_liquid = load_default_eos_interpolators( + spider_dir=spider_dir + ) + + S_sol = invert_to_entropy_along_profile(P_sol, T_sol, S_axis_solid, T_of_SP_solid) + S_liq = invert_to_entropy_along_profile(P_liq, T_liq, S_axis_liquid, T_of_SP_liquid) + + P_common, S_sol_common, S_liq_common = build_common_entropy_grid( + P_sol, S_sol, P_liq, S_liq, n_common=n + ) + + validate_entropy_export_arrays( + P_common, + S_sol_common, + S_liq_common, + model_name=model_name, + ) + + save_entropy_table_with_header( + out_dir / 'solidus_P-S.dat', + P_common, + S_sol_common, + ) + + save_entropy_table_with_header( + out_dir / 'liquidus_P-S.dat', + P_common, + S_liq_common, + ) + + print_model_summary( + model_name, + P_sol, + T_sol, + P_liq, + T_liq, + P_common, + S_sol_common, + S_liq_common, + ) + log.info(f"Saved to: {out_dir.resolve()}") + + return { + 'P_sol': P_sol, + 'T_sol': T_sol, + 'P_liq': P_liq, + 'T_liq': T_liq, + 'S_sol': S_sol, + 'S_liq': S_liq, + 'P_entropy_common': P_common, + 'S_sol_common': S_sol_common, + 'S_liq_common': S_liq_common, + } + + +# ============================================================================= +# BATCH EXPORTER +# ============================================================================= + + +def export_all_models( + out_root: Path | str = MELTING_DIR, + n: int = 2000, + spider_dir: Path | str | None = None, +): + r""" + Export all supported melting parametrizations. + """ + for model in SUPPORTED_MODELS: + if model == 'katz_2003': + _ = export_model_curves( + model, out_root=out_root, n=n, X_h2o=30.0, spider_dir=spider_dir + ) + elif model == 'lin_2024': + _ = export_model_curves( + model, out_root=out_root, n=n, fO2=-4.0, spider_dir=spider_dir + ) + elif model == 'hirschmann_2000': + _ = export_model_curves( + model, out_root=out_root, n=n, Pmax=10.0, spider_dir=spider_dir + ) + elif model == 'fei_2021': + _ = export_model_curves( + model, out_root=out_root, n=n, Pmin=1.0, spider_dir=spider_dir + ) + elif model == 'stixrude_2014': + _ = export_model_curves( + model, out_root=out_root, n=n, Pmin=1.0, spider_dir=spider_dir + ) + else: + _ = export_model_curves(model, out_root=out_root, n=n, spider_dir=spider_dir) + + +# ============================================================================= +# COMMAND-LINE INTERFACE +# ============================================================================= + + +def parse_args(): + parser = argparse.ArgumentParser( + description=( + 'Export solidus and liquidus melting curves in P-T and P-S space ' + 'for one or more literature parametrizations.' + ), + epilog=( + 'Examples:\n' + ' python solidus_func.py --all\n' + ' python solidus_func.py --katz2003 --X-h2o 30\n' + ' python solidus_func.py --lin2024 --fO2 -4\n' + ' python solidus_func.py --model wolf_bower_2018\n' + ), + formatter_class=argparse.RawTextHelpFormatter, + ) + + parser.add_argument( + '--all', + action='store_true', + help='Export all supported models.', + ) + + parser.add_argument( + '--model', + type=str, + default=None, + choices=SUPPORTED_MODELS, + help='Export a single model by internal name.', + ) + + parser.add_argument( + '--katz2003', + action='store_true', + help='Export Katz et al. (2003). Requires --X-h2o.', + ) + parser.add_argument( + '--lin2024', + action='store_true', + help='Export Lin et al. (2024). Requires --fO2.', + ) + parser.add_argument( + '--wolfbower2018', + action='store_true', + help='Export Wolf & Bower (2018).', + ) + parser.add_argument( + '--andrault2011', + action='store_true', + help='Export Andrault et al. (2011).', + ) + parser.add_argument( + '--monteux2016', + action='store_true', + help='Export Monteux et al. (2016).', + ) + parser.add_argument( + '--fei2021', + action='store_true', + help='Export Fei et al. (2021).', + ) + parser.add_argument( + '--belonoshko2005', + action='store_true', + help='Export Belonoshko et al. (2005).', + ) + parser.add_argument( + '--fiquet2010', + action='store_true', + help='Export Fiquet et al. (2010).', + ) + parser.add_argument( + '--hirschmann2000', + action='store_true', + help='Export Hirschmann (2000).', + ) + parser.add_argument( + '--stixrude2014', + action='store_true', + help='Export Stixrude (2014).', + ) + + parser.add_argument( + '--out-root', + type=str, + default=str(MELTING_DIR), + help='Root directory where output folders will be created.', + ) + + parser.add_argument( + '--spider-dir', + type=str, + default=None, + help='Path to the SPIDER root directory containing lookup_data/.', + ) + + parser.add_argument( + '--Pmin', + type=float, + default=0.0, + help='Minimum pressure in GPa.', + ) + parser.add_argument( + '--Pmax', + type=float, + default=1000.0, + help='Maximum pressure in GPa.', + ) + parser.add_argument( + '-n', + type=int, + default=2000, + help='Number of pressure samples.', + ) + + parser.add_argument( + '--X-h2o', + dest='X_h2o', + type=float, + default=None, + help='Water content parameter for Katz (2003). Required for --katz2003.', + ) + parser.add_argument( + '--fO2', + type=float, + default=None, + help='Oxygen fugacity offset for Lin (2024). Required for --lin2024.', + ) + + return parser.parse_args() + + +def resolve_requested_model(args) -> str | None: + """ + Resolve which single-model shortcut flag was requested. + """ + shortcut_map = { + 'katz2003': 'katz_2003', + 'lin2024': 'lin_2024', + 'wolfbower2018': 'wolf_bower_2018', + 'andrault2011': 'andrault_2011', + 'monteux2016': 'monteux_2016', + 'fei2021': 'fei_2021', + 'belonoshko2005': 'belonoshko_2005', + 'fiquet2010': 'fiquet_2010', + 'hirschmann2000': 'hirschmann_2000', + 'stixrude2014': 'stixrude_2014', + } + + chosen = [model for flag, model in shortcut_map.items() if getattr(args, flag)] + + if len(chosen) > 1: + raise SystemExit('Error: please select only one model shortcut flag at a time.') + + if len(chosen) == 1: + return chosen[0] + + return None + + +def export_one_model_from_cli(model_name: str, args): + """ + Export a single model, applying model-specific defaults and enforcing + required parameters. + """ + kwargs = {} + Pmin = args.Pmin + Pmax = args.Pmax + n = args.n + + if model_name == 'katz_2003': + if args.X_h2o is None: + raise SystemExit( + 'Error: --X-h2o is required when using Katz (2003).\n' + 'Example: python solidus_func.py --katz2003 --X-h2o 30' + ) + kwargs['X_h2o'] = args.X_h2o + + elif model_name == 'lin_2024': + if args.fO2 is None: + raise SystemExit( + 'Error: --fO2 is required when using Lin (2024).\n' + 'Example: python solidus_func.py --lin2024 --fO2 -4' + ) + kwargs['fO2'] = args.fO2 + + elif model_name == 'hirschmann_2000': + if args.Pmax == 1000.0: + Pmax = 10.0 + + elif model_name == 'fei_2021': + if args.Pmin == 0.0: + Pmin = 1.0 + + elif model_name == 'stixrude_2014': + if args.Pmin == 0.0: + Pmin = 1.0 + + _ = export_model_curves( + model_name, + out_root=args.out_root, + Pmin=Pmin, + Pmax=Pmax, + n=n, + spider_dir=args.spider_dir, + **kwargs, + ) + + +def main(): + + logging.basicConfig( + level=logging.INFO, + format="%(message)s", + ) + + args = parse_args() + + shortcut_model = resolve_requested_model(args) + explicit_model = args.model + + if args.all: + if explicit_model is not None or shortcut_model is not None: + raise SystemExit( + 'Error: please use either --all or a single model selection, not both.' + ) + try: + export_all_models(out_root=args.out_root, n=args.n, spider_dir=args.spider_dir) + except FileNotFoundError as exc: + raise SystemExit(f'Error: {exc}') from exc + return + + selected_models = [m for m in [explicit_model, shortcut_model] if m is not None] + + if len(selected_models) == 0: + models_list = "\n".join( + f" - --{key:<20} ({name})" + for key, name in sorted(DISPLAY_NAMES.items()) + ) + + raise SystemExit( + "Error: no model selected.\n\n" + "Use --all or choose one model with --model or a shortcut like --katz_2003.\n\n" + "Available models:\n" + f"{models_list}" + ) + + if len(selected_models) > 1: + raise SystemExit('Error: please choose only one of --model or one shortcut flag.') + + try: + export_one_model_from_cli(selected_models[0], args) + except FileNotFoundError as exc: + raise SystemExit(f'Error: {exc}') from exc + + +if __name__ == '__main__': + main()