From 86e288d66ae92ce5c529fb08c0c59d72959479ec Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Mon, 16 Mar 2026 20:10:19 +0100 Subject: [PATCH 01/13] loads P-T and P-S melting curves --- src/proteus/interior/solidus_func.py | 641 +++++++++++++++++++++++++++ 1 file changed, 641 insertions(+) create mode 100644 src/proteus/interior/solidus_func.py diff --git a/src/proteus/interior/solidus_func.py b/src/proteus/interior/solidus_func.py new file mode 100644 index 000000000..64c2a0ce9 --- /dev/null +++ b/src/proteus/interior/solidus_func.py @@ -0,0 +1,641 @@ +from __future__ import annotations + +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +from scipy.interpolate import RegularGridInterpolator, interp1d + + +# ========================================================= +# Common pressure grid helper +# ========================================================= +def make_pressure_grid(Pmin=0.0, Pmax=1000.0, n=500): + """ + Create a pressure grid in GPa. + """ + return np.linspace(Pmin, Pmax, n) + + +# ========================================================= +# Generic helpers +# ========================================================= +def solidus_from_liquidus_stixrude(T_liq): + """ + Approximate solidus from liquidus using the Stixrude ratio. + """ + return T_liq / (1.0 - np.log(0.79)) + + +def liquidus_from_solidus_stixrude(T_sol): + """ + Approximate liquidus from solidus using the inverse Stixrude ratio. + """ + return T_sol * (1.0 - np.log(0.79)) + + +# ========================================================= +# Andrault et al. (2011) +# Formula expects P in Pa in the coefficients below +# ========================================================= +def andrault_2011(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): + P = make_pressure_grid(Pmin, Pmax, n) # GPa + P_pa = P * 1e9 # Pa + + 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 + + +# ========================================================= +# Fei et al. (2021) +# Expects P in GPa +# ========================================================= +def fei_2021(kind="liquidus", Pmin=1.0, Pmax=1000.0, n=500): + 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 + + +# ========================================================= +# Belonoshko et al. (2005) +# Expects P in GPa +# ========================================================= +def belonoshko_2005(kind="liquidus", Pmin=0.0, Pmax=1000.0, n=500): + 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 + + +# ========================================================= +# Fiquet et al. (2010) +# Expects P in GPa +# ========================================================= +def fiquet_2010(kind="liquidus", Pmin=0.0, Pmax=1000.0, n=500): + P = make_pressure_grid(Pmin, Pmax, n) + T_liq = np.zeros_like(P, dtype=float) + + low = P <= 20.0 + high = P > 20.0 + + 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 + + +# ========================================================= +# Monteux et al. (2016) +# Original coefficients are in Pa +# ========================================================= +def monteux_2016(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): + P = make_pressure_grid(Pmin, Pmax, n) # GPa + P_pa = P * 1e9 # Pa + + 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 + + +# ========================================================= +# Hirschmann (2000) +# Valid only at relatively low pressure +# ========================================================= +def hirschmann_2000(kind="solidus", Pmin=0.0, Pmax=10.0, n=500): + 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 + + +# ========================================================= +# Stixrude (2014) +# Expects P in GPa +# ========================================================= +def stixrude_2014(kind="liquidus", Pmin=1.0, Pmax=1000.0, n=500): + 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 + + +# ========================================================= +# Wolf & Bower (2018) / piecewise fits +# Expects P in GPa +# ========================================================= +def wolf_bower_2018(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): + 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 + + +# ========================================================= +# Katz (2003) +# Applies the same hydrous depression to both curves +# ========================================================= +def katz_2003(kind="solidus", X_h2o=30.0, Pmin=0.0, Pmax=1000.0, n=500): + 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 + + +# ========================================================= +# Lin et al. (2024) +# ========================================================= +def lin_2024(kind="solidus", fO2=-4.0, Pmin=0.0, Pmax=1000.0, n=500): + 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 + + +# ========================================================= +# Keep only the main physical interval where solidus < liquidus +# Works for Andrault and for high-P crossovers +# ========================================================= +def truncate_to_physical_interval(func): + 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] + elif kind == "liquidus": + return P_liq[main_block], T_liq[main_block] + else: + raise ValueError("kind must be 'solidus' or 'liquidus'") + + return wrapped + + +# Wrapped physical versions +andrault_2011_cut = truncate_to_physical_interval(andrault_2011) +monteux_2016_cut = truncate_to_physical_interval(monteux_2016) +wolf_bower_2018_cut = truncate_to_physical_interval(wolf_bower_2018) +katz_2003_cut = truncate_to_physical_interval(katz_2003) + + +# ========================================================= +# Dispatcher +# ========================================================= +def get_melting_curves(model_name, Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): + """ + Return physical solidus and liquidus curves for a given model name. + + Returns + ------- + P_sol, T_sol, P_liq, T_liq + """ + models = { + "andrault_2011": andrault_2011_cut, + "monteux_2016": monteux_2016_cut, + "wolf_bower_2018": wolf_bower_2018_cut, + "katz_2003": katz_2003_cut, + "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 + + +# ========================================================= +# Save helpers +# ========================================================= +def save_PT_table(path: Path, P_gpa, T_k): + """ + Save a P-T profile with header: + #pressure temperature + """ + data = np.column_stack([P_gpa, T_k]) + np.savetxt(path, data, fmt="%.18e %.18e", header="pressure temperature", comments="#") + + +# ================== EOS settings ================== +eos_solid_path = Path("temperature_solid.dat") +eos_liquid_path = Path("temperature_melt.dat") + +nP = 2020 +nS_solid = 125 +nS_liquid = 95 +skip_header = 5 + +SCALE_P_EOS = 1e9 +SCALE_T_EOS = 1.0 +SCALE_S_SOLID_EOS = 4.82426684604467e6 +SCALE_S_LIQUID_EOS = 4.805046659407042e6 + +SCALE_P_OUT = 1_000_000_000.0 +SCALE_S_OUT = 4_824_266.84604467 + +COMMON_HEADER = "\n".join([ + "# 5 400", + "# Pressure, Entropy, Quantity", + "# column * scaling factor should be SI units", + "# scaling factors (constant) for each column given on line below", + "# 1000000000.0 4824266.84604467", +]) + + +def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): + 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 invert_to_entropy_along_profile(P_gpa, T_k, S_axis, T_of_SP): + """ + Invert T(S,P) -> S along a P-T profile. + """ + 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] + + if T_i < T_sorted[0] or T_i > T_sorted[-1]: + continue + + try: + f = interp1d(T_sorted, S_sorted, kind="linear", assume_sorted=True) + S_out[i] = float(f(T_i)) + except Exception: + f = interp1d(T_sorted, S_sorted, kind="nearest", assume_sorted=True) + S_out[i] = float(f(T_i)) + + return S_out + + +def save_entropy_table_with_header(path: Path, P_gpa, S_jpk): + """ + Save entropy table in the scaled 2-column format. + """ + P_pa = P_gpa * 1e9 + data = np.column_stack([P_pa / SCALE_P_OUT, S_jpk / SCALE_S_OUT]) + np.savetxt(path, data, fmt="%.18e %.18e", header=COMMON_HEADER, comments="") + + +# Load EOS interpolators once +S_axis_solid, P_axis_solid, T_of_SP_solid = load_eos_T_of_SP( + eos_solid_path, nS_solid, SCALE_S_SOLID_EOS +) + +S_axis_liquid, P_axis_liquid, T_of_SP_liquid = load_eos_T_of_SP( + eos_liquid_path, nS_liquid, SCALE_S_LIQUID_EOS +) + + +# ========================================================= +# Main exporter: save P-T and P-S +# ========================================================= +def export_model_curves(model_name, out_root="outputs_entropy_curves", + Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): + """ + Generate melting curves for one model, save: + - solidus.dat + - liquidus.dat + - solidus_entropy.dat + - liquidus_entropy.dat + + Returns + ------- + dict with arrays + """ + 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 P-T + save_PT_table(out_dir / "solidus.dat", P_sol, T_sol) + save_PT_table(out_dir / "liquidus.dat", P_liq, T_liq) + + # Convert to entropy + 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 + ) + + # Remove NaNs before saving P-S files + mask_sol = np.isfinite(S_sol) + mask_liq = np.isfinite(S_liq) + + save_entropy_table_with_header( + out_dir / "solidus_entropy.dat", + P_sol[mask_sol], + S_sol[mask_sol], + ) + + save_entropy_table_with_header( + out_dir / "liquidus_entropy.dat", + P_liq[mask_liq], + S_liq[mask_liq], + ) + + 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, + } + + +# ========================================================= +# Convenience: export all models +# ========================================================= +def export_all_models(out_root="outputs_entropy_curves"): + models = [ + "andrault_2011", + "monteux_2016", + "wolf_bower_2018", + "katz_2003", + "fei_2021", + "belonoshko_2005", + "fiquet_2010", + "hirschmann_2000", + "stixrude_2014", + "lin_2024", + ] + + for model in models: + print(f"[INFO] Exporting {model}") + + if model == "katz_2003": + export_model_curves(model, out_root=out_root, X_h2o=30.0) + elif model == "lin_2024": + export_model_curves(model, out_root=out_root, fO2=-4.0) + elif model == "hirschmann_2000": + export_model_curves(model, out_root=out_root, Pmax=10.0) + elif model == "fei_2021": + export_model_curves(model, out_root=out_root, Pmin=1.0) + elif model == "stixrude_2014": + export_model_curves(model, out_root=out_root, Pmin=1.0) + else: + export_model_curves(model, out_root=out_root) + + print("[INFO] Done.") + + +# ========================================================= +# Optional quick plot +# ========================================================= +def show_melting_curves(model_name, Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): + P_sol, T_sol, P_liq, T_liq = get_melting_curves( + model_name, Pmin=Pmin, Pmax=Pmax, n=n, **kwargs + ) + + plt.figure(figsize=(6, 6)) + plt.plot(T_sol, P_sol, label="Solidus") + plt.plot(T_liq, P_liq, label="Liquidus") + plt.xlabel("Temperature (K)") + plt.ylabel("Pressure (GPa)") + plt.gca().invert_yaxis() + plt.title(model_name) + plt.legend() + plt.tight_layout() + plt.show() + + return P_sol, T_sol, P_liq, T_liq + +def show_entropy_curves(model_name, Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): + """ + Plot melting curves in P–S space. + """ + + P_sol, T_sol, P_liq, T_liq = get_melting_curves( + model_name, Pmin=Pmin, Pmax=Pmax, n=n, **kwargs + ) + + # convert to entropy + 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 + ) + + # remove NaNs before plotting + mask_sol = np.isfinite(S_sol) + mask_liq = np.isfinite(S_liq) + + plt.figure(figsize=(6,6)) + + plt.plot(S_sol[mask_sol], P_sol[mask_sol], label="Solidus") + plt.plot(S_liq[mask_liq], P_liq[mask_liq], label="Liquidus") + + plt.xlabel("Entropy (J kg$^{-1}$ K$^{-1}$)") + plt.ylabel("Pressure (GPa)") + plt.gca().invert_yaxis() + + plt.title(f"{model_name} (P–S)") + plt.legend() + plt.tight_layout() + plt.show() + + return ( + P_sol[mask_sol], S_sol[mask_sol], + P_liq[mask_liq], S_liq[mask_liq], + ) + +# ========================================================= +# Example usage +# ========================================================= +# One model: +# export_model_curves("andrault_2011") +# export_model_curves("katz_2003", X_h2o=30.0) +# export_model_curves("lin_2024", fO2=-4.0) + +# All models: +export_all_models() + +# Quick plot: +# show_melting_curves("wolf_bower_2018") +# show_melting_curves("andrault_2011") +# show_melting_curves("katz_2003", X_h2o=30.0) +# show_melting_curves("lin_2024", fO2=-4.0) +# show_melting_curves("fei_2021") +# show_melting_curves("belonoshko_2005") +# show_melting_curves("fiquet_2010") +# show_melting_curves("hirschmann_2000", Pmax=10.0) +# show_melting_curves("stixrude_2014", Pmin=1.0) +show_entropy_curves("wolf_bower_2018") +show_entropy_curves("andrault_2011") +show_entropy_curves("katz_2003", X_h2o=30.0) +show_entropy_curves("lin_2024", fO2=-4.0) +show_entropy_curves("fei_2021") +show_entropy_curves("belonoshko_2005") +show_entropy_curves("fiquet_2010") +show_entropy_curves("hirschmann_2000", Pmax=10.0) +show_entropy_curves("stixrude_2014", Pmin=1.0) From ade1a4aac0cb42e14c92583d70ee80f73ecc09ce Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Tue, 17 Mar 2026 17:46:38 +0100 Subject: [PATCH 02/13] version with docstrings --- src/proteus/interior/solidus_func.py | 1134 +++++++++++++++++++++----- 1 file changed, 920 insertions(+), 214 deletions(-) diff --git a/src/proteus/interior/solidus_func.py b/src/proteus/interior/solidus_func.py index 64c2a0ce9..31054240a 100644 --- a/src/proteus/interior/solidus_func.py +++ b/src/proteus/interior/solidus_func.py @@ -1,44 +1,184 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +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. + +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: 10.1016/j.epsl.2011.02.006 +- Monteux et al. (2016), DOI: 10.1016/j.epsl.2016.05.010 +- Lin et al. (2024), DOI: 10.1038/s41561-024-01495-1 +- Hirschmann (2000), DOI: 10.1029/2000GC000070 +- Katz et al. (2003), DOI: 10.1029/2002GC000433 +- Stixrude (2014), DOI: 10.1098/rsta.2013.0076 +- Fei et al. (2021), DOI: 10.1038/s41467-021-21170-y +- Belonoshko et al. (2005), DOI: 10.1103/PhysRevLett.94.195701 +- Fiquet et al. (2010), DOI: 10.1126/science.1192448 +""" from __future__ import annotations from pathlib import Path -import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import RegularGridInterpolator, interp1d +# ============================================================================= +# GENERAL HELPERS +# ============================================================================= -# ========================================================= -# Common pressure grid helper -# ========================================================= -def make_pressure_grid(Pmin=0.0, Pmax=1000.0, n=500): - """ - Create a pressure grid in GPa. +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. + + Notes + ----- + This helper is used by all melting parametrizations so that they share + a consistent pressure sampling unless otherwise specified. """ return np.linspace(Pmin, Pmax, n) -# ========================================================= -# Generic helpers -# ========================================================= -def solidus_from_liquidus_stixrude(T_liq): - """ - Approximate solidus from liquidus using the Stixrude ratio. +def solidus_from_liquidus_stixrude(T_liq: np.ndarray) -> np.ndarray: + r""" + Estimate the solidus from a liquidus using the Stixrude ratio. + + Parameters + ---------- + T_liq : np.ndarray + Liquidus temperature in K. + + Returns + ------- + np.ndarray + Estimated solidus temperature in K. + + Formula + ------- + .. math:: + + T_{\mathrm{sol}} = \frac{T_{\mathrm{liq}}}{1 - \ln(0.79)} + + Notes + ----- + This is used in cases where only a liquidus parametrization is available + and an approximate solidus is needed. + + Reference + --------- + Stixrude (2014), DOI: 10.1098/rsta.2013.0076 """ return T_liq / (1.0 - np.log(0.79)) -def liquidus_from_solidus_stixrude(T_sol): - """ - Approximate liquidus from solidus using the inverse Stixrude ratio. +def liquidus_from_solidus_stixrude(T_sol: np.ndarray) -> np.ndarray: + r""" + Estimate the liquidus from a solidus using the inverse Stixrude ratio. + + Parameters + ---------- + T_sol : np.ndarray + Solidus temperature in K. + + Returns + ------- + np.ndarray + Estimated liquidus temperature in K. + + Formula + ------- + .. math:: + + T_{\mathrm{liq}} = T_{\mathrm{sol}} \left(1 - \ln(0.79)\right) + + Notes + ----- + This is used in cases where only a solidus parametrization is available + and an approximate liquidus is needed. + + Reference + --------- + Stixrude (2014), DOI: 10.1098/rsta.2013.0076 """ return T_sol * (1.0 - np.log(0.79)) -# ========================================================= -# Andrault et al. (2011) -# Formula expects P in Pa in the coefficients below -# ========================================================= -def andrault_2011(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): +# ============================================================================= +# 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). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + The parametrization is written as + + .. math:: + + T(P) = T_0 \left( \frac{P}{a} + 1 \right)^{1/c} + + where the published coefficients are defined using pressure in Pa. + In this implementation the user-facing pressure grid is generated in GPa + and internally converted to Pa. + + Reference + --------- + Andrault et al. (2011), DOI: 10.1016/j.epsl.2011.02.006 + """ P = make_pressure_grid(Pmin, Pmax, n) # GPa P_pa = P * 1e9 # Pa @@ -53,11 +193,43 @@ def andrault_2011(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Fei et al. (2021) -# Expects P in GPa -# ========================================================= -def fei_2021(kind="liquidus", Pmin=1.0, Pmax=1000.0, n=500): +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). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + The liquidus is given by + + .. math:: + + T_{\mathrm{liq}}(P) = 6000 \left(\frac{P}{140}\right)^{0.26} + + with pressure in GPa. + + If ``kind="solidus"``, the solidus is estimated from the liquidus using + :func:`solidus_from_liquidus_stixrude`. + + Reference + --------- + Fei et al. (2021), DOI: 10.1038/s41467-021-21170-y + """ P = make_pressure_grid(Pmin, Pmax, n) T_liq = 6000.0 * (P / 140.0) ** 0.26 @@ -71,11 +243,43 @@ def fei_2021(kind="liquidus", Pmin=1.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Belonoshko et al. (2005) -# Expects P in GPa -# ========================================================= -def belonoshko_2005(kind="liquidus", Pmin=0.0, Pmax=1000.0, n=500): +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). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + The liquidus is given by + + .. math:: + + T_{\mathrm{liq}}(P) = 1831 \left(1 + \frac{P}{4.6}\right)^{0.33} + + with pressure in GPa. + + If ``kind="solidus"``, the solidus is estimated from the liquidus using + :func:`solidus_from_liquidus_stixrude`. + + Reference + --------- + Belonoshko et al. (2005), DOI: 10.1103/PhysRevLett.94.195701 + """ P = make_pressure_grid(Pmin, Pmax, n) T_liq = 1831.0 * (1.0 + P / 4.6) ** 0.33 @@ -89,11 +293,51 @@ def belonoshko_2005(kind="liquidus", Pmin=0.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Fiquet et al. (2010) -# Expects P in GPa -# ========================================================= -def fiquet_2010(kind="liquidus", Pmin=0.0, Pmax=1000.0, n=500): +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). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + The liquidus is implemented as a two-branch fit: + + for :math:`P \leq 20\ \mathrm{GPa}`: + + .. math:: + + T_{\mathrm{liq}} = 1982.1 \left(\frac{P}{6.594} + 1\right)^{1/5.374} + + and for :math:`P > 20\ \mathrm{GPa}`: + + .. math:: + + T_{\mathrm{liq}} = 78.74 \left(\frac{P}{0.004056} + 1\right)^{1/2.44} + + where pressure is in GPa. + + If ``kind="solidus"``, the solidus is estimated from the liquidus using + :func:`solidus_from_liquidus_stixrude`. + + Reference + --------- + Fiquet et al. (2010), DOI: 10.1126/science.1192448 + """ P = make_pressure_grid(Pmin, Pmax, n) T_liq = np.zeros_like(P, dtype=float) @@ -113,11 +357,42 @@ def fiquet_2010(kind="liquidus", Pmin=0.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Monteux et al. (2016) -# Original coefficients are in Pa -# ========================================================= -def monteux_2016(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): +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). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + Both solidus and liquidus are implemented as piecewise power-law fits: + + .. math:: + + T(P) = T_0 \left(\frac{P}{a} + 1\right)^{1/c} + + with one coefficient set below 20 GPa and another above 20 GPa. + The published coefficients use pressure in Pa, so pressure is converted + internally from GPa to Pa before evaluation. + + Reference + --------- + Monteux et al. (2016), DOI: 10.1016/j.epsl.2016.05.010 + """ P = make_pressure_grid(Pmin, Pmax, n) # GPa P_pa = P * 1e9 # Pa @@ -147,11 +422,46 @@ def monteux_2016(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Hirschmann (2000) -# Valid only at relatively low pressure -# ========================================================= -def hirschmann_2000(kind="solidus", Pmin=0.0, Pmax=10.0, n=500): +def hirschmann_2000(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 10.0, n: int = 500): + r""" + Melting curve from Hirschmann (2000). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + The fit is a quadratic polynomial in pressure: + + .. math:: + + T_{\mathrm{^\circ C}}(P) = A_1 + A_2 P + A_3 P^2 + + The result is then converted to Kelvin via + + .. math:: + + T_{\mathrm{K}} = T_{\mathrm{^\circ C}} + 273.15 + + This parametrization is only intended for relatively low pressures. + + Reference + --------- + Hirschmann (2000), DOI: 10.1029/2000GC000070 + """ P = make_pressure_grid(Pmin, Pmax, n) coeffs = { @@ -169,11 +479,43 @@ def hirschmann_2000(kind="solidus", Pmin=0.0, Pmax=10.0, n=500): return P, T -# ========================================================= -# Stixrude (2014) -# Expects P in GPa -# ========================================================= -def stixrude_2014(kind="liquidus", Pmin=1.0, Pmax=1000.0, n=500): +def stixrude_2014(kind: str = "liquidus", Pmin: float = 1.0, Pmax: float = 1000.0, n: int = 500): + r""" + Melting curve based on Stixrude (2014). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + The liquidus is given by + + .. math:: + + T_{\mathrm{liq}}(P) = 5400 \left(\frac{P}{140}\right)^{0.480} + + with pressure in GPa. + + If ``kind="solidus"``, the solidus is estimated from the liquidus using + :func:`solidus_from_liquidus_stixrude`. + + Reference + --------- + Stixrude (2014), DOI: 10.1098/rsta.2013.0076 + """ P = make_pressure_grid(Pmin, Pmax, n) T_liq = 5400.0 * (P / 140.0) ** 0.480 @@ -187,11 +529,58 @@ def stixrude_2014(kind="liquidus", Pmin=1.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Wolf & Bower (2018) / piecewise fits -# Expects P in GPa -# ========================================================= -def wolf_bower_2018(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): +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. + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + The curve is piecewise linear in pressure, with breakpoints + :math:`P_{\mathrm{cp1}}` and :math:`P_{\mathrm{cp2}}`: + + .. math:: + + T(P) = + \begin{cases} + s_1 P + c_1, & P < P_{\mathrm{cp1}} \\ + s_2 P + c_2, & P_{\mathrm{cp1}} \leq P < P_{\mathrm{cp2}} \\ + s_3 P + c_3, & P \geq P_{\mathrm{cp2}} + \end{cases} + + where continuity is enforced through + + .. math:: + + c_2 = c_1 + (s_1 - s_2) P_{\mathrm{cp1}} + + .. math:: + + c_3 = c_2 + (s_2 - s_3) P_{\mathrm{cp2}} + + Notes + ----- + This implementation uses coefficients already encoded in the script. + + Reference + --------- + Wolf & Bower (2018)-style parametrization used in the model workflow. + """ P = make_pressure_grid(Pmin, Pmax, n) params = { @@ -235,11 +624,51 @@ def wolf_bower_2018(kind="solidus", Pmin=0.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Katz (2003) -# Applies the same hydrous depression to both curves -# ========================================================= -def katz_2003(kind="solidus", X_h2o=30.0, Pmin=0.0, Pmax=1000.0, n=500): +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). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + X_h2o : float + Water content parameter used in the hydrous temperature depression. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + Starting from the anhydrous Wolf & Bower curve, a constant depression is + applied at all pressures: + + .. math:: + + \Delta T = K X_{\mathrm{H_2O}}^{\gamma} + + .. math:: + + T = T_{\mathrm{anhydrous}} - \Delta T + + where in this implementation + + .. math:: + + \gamma = 0.75,\quad K = 43 + + Reference + --------- + Katz et al. (2003), DOI: 10.1029/2002GC000433 + """ gamma = 0.75 K = 43.0 @@ -253,10 +682,50 @@ def katz_2003(kind="solidus", X_h2o=30.0, Pmin=0.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Lin et al. (2024) -# ========================================================= -def lin_2024(kind="solidus", fO2=-4.0, Pmin=0.0, Pmax=1000.0, n=500): +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). + + Parameters + ---------- + kind : {"solidus", "liquidus"} + Which branch to evaluate. + fO2 : float + Oxygen fugacity offset parameter used in the solidus shift. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of pressure samples. + + Returns + ------- + P : np.ndarray + Pressure in GPa. + T : np.ndarray + Temperature in K. + + Formula + ------- + The anhydrous solidus is first taken from :func:`wolf_bower_2018` and then + shifted by + + .. math:: + + \Delta T = \frac{340}{3.2} (2 - f\mathrm{O}_2) + + so that + + .. math:: + + T_{\mathrm{sol}} = T_{\mathrm{anhydrous}} + \Delta T + + If ``kind="liquidus"``, the liquidus is estimated from the shifted solidus + using :func:`liquidus_from_solidus_stixrude`. + + Reference + --------- + Lin et al. (2024), DOI: 10.1038/s41561-024-01495-1 + """ P, T_anhydrous = wolf_bower_2018(kind="solidus", Pmin=Pmin, Pmax=Pmax, n=n) delta_T = (340.0 / 3.2) * (2.0 - fO2) @@ -272,11 +741,39 @@ def lin_2024(kind="solidus", fO2=-4.0, Pmin=0.0, Pmax=1000.0, n=500): return P, T -# ========================================================= -# Keep only the main physical interval where solidus < liquidus -# Works for Andrault and for high-P crossovers -# ========================================================= +# ============================================================================= +# PHYSICAL-INTERVAL FILTER +# ============================================================================= + def truncate_to_physical_interval(func): + r""" + Wrap a melting-curve function so only the main interval with + :math:`T_{\mathrm{sol}} < T_{\mathrm{liq}}` is retained. + + Parameters + ---------- + func : callable + A function returning a melting curve in the form ``P, T`` for + ``kind="solidus"`` or ``kind="liquidus"``. + + Returns + ------- + callable + Wrapped function returning only the largest contiguous physically valid + pressure interval. + + Notes + ----- + Some parametrizations can cross or become unphysical at high pressure. + This wrapper evaluates both solidus and liquidus, identifies all points + where + + .. math:: + + T_{\mathrm{sol}} < T_{\mathrm{liq}} + + and keeps only the largest contiguous block of such points. + """ 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) @@ -301,23 +798,42 @@ def wrapped(kind="solidus", Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): return wrapped -# Wrapped physical versions +# Wrapped versions for models where physical truncation is needed. andrault_2011_cut = truncate_to_physical_interval(andrault_2011) monteux_2016_cut = truncate_to_physical_interval(monteux_2016) wolf_bower_2018_cut = truncate_to_physical_interval(wolf_bower_2018) katz_2003_cut = truncate_to_physical_interval(katz_2003) -# ========================================================= -# Dispatcher -# ========================================================= -def get_melting_curves(model_name, Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): - """ - Return physical solidus and liquidus curves for a given model name. +# ============================================================================= +# MODEL DISPATCHER +# ============================================================================= + +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. + + Parameters + ---------- + model_name : str + Identifier of the melting parametrization. + Pmin, Pmax : float + Pressure range in GPa. + n : int + Number of sampling points. + **kwargs + Additional keyword arguments passed to the selected model + (for example ``X_h2o`` or ``fO2``). Returns ------- - P_sol, T_sol, P_liq, T_liq + P_sol, T_sol, P_liq, T_liq : tuple of np.ndarray + Solidus and liquidus curves in P–T space. + + Raises + ------ + ValueError + If the requested model name is unknown. """ models = { "andrault_2011": andrault_2011_cut, @@ -342,35 +858,72 @@ def get_melting_curves(model_name, Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): return P_sol, T_sol, P_liq, T_liq -# ========================================================= -# Save helpers -# ========================================================= -def save_PT_table(path: Path, P_gpa, T_k): - """ - Save a P-T profile with header: - #pressure temperature +# ============================================================================= +# OUTPUT HELPERS +# ============================================================================= + +def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): + r""" + Save a pressure–temperature table to disk. + + Parameters + ---------- + path : pathlib.Path + Output filename. + P_gpa : np.ndarray + Pressure in GPa. + T_k : np.ndarray + Temperature in K. + + Notes + ----- + The file is written as two columns: + + .. code-block:: text + + pressure temperature + + using a simple header compatible with later inspection. """ data = np.column_stack([P_gpa, T_k]) np.savetxt(path, data, fmt="%.18e %.18e", header="pressure temperature", comments="#") -# ================== EOS settings ================== +# ============================================================================= +# EOS LOOKUP TABLE SETTINGS +# ============================================================================= +# +# The files temperature_solid.dat and temperature_melt.dat contain structured +# tables in (S, P) space. Each row stores: +# +# pressure, entropy, temperature +# +# in scaled units. The constants below decode those tables back to SI units. + eos_solid_path = Path("temperature_solid.dat") eos_liquid_path = Path("temperature_melt.dat") +# Number of pressure points in the EOS tables. nP = 2020 + +# Number of entropy points for the solid and liquid tables, respectively. nS_solid = 125 nS_liquid = 95 + +# Number of header lines to skip when reading the raw EOS files. skip_header = 5 +# Scaling factors used by the EOS tables. SCALE_P_EOS = 1e9 SCALE_T_EOS = 1.0 SCALE_S_SOLID_EOS = 4.82426684604467e6 SCALE_S_LIQUID_EOS = 4.805046659407042e6 +# Scaling factors used when exporting P–S tables back to SPIDER-like format. SCALE_P_OUT = 1_000_000_000.0 SCALE_S_OUT = 4_824_266.84604467 +# Header used for exported SPIDER-style entropy tables. COMMON_HEADER = "\n".join([ "# 5 400", "# Pressure, Entropy, Quantity", @@ -381,6 +934,37 @@ def save_PT_table(path: Path, P_gpa, T_k): def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): + r""" + Load an EOS lookup table and build an interpolator for :math:`T(S, P)`. + + Parameters + ---------- + eos_path : pathlib.Path + Path to the EOS file. + nS : int + Number of entropy grid points in the table. + scale_S_axis : float + Scaling factor needed to convert the stored entropy column to SI units. + + Returns + ------- + S_axis : np.ndarray + Entropy axis in :math:`\mathrm{J\,kg^{-1}\,K^{-1}}`. + P_axis_GPa : np.ndarray + Pressure axis in GPa. + T_interp : scipy.interpolate.RegularGridInterpolator + Interpolator returning temperature in K for points in ``(S, P)``. + + Notes + ----- + The EOS files are stored in a flattened 3-column format and reshaped into + + .. math:: + + (n_S, n_P, 3) + + where the three columns correspond to pressure, entropy, and temperature. + """ raw = np.genfromtxt(eos_path, skip_header=skip_header) resh = raw.reshape((nS, nP, 3)) @@ -397,9 +981,42 @@ def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): return S_axis, P_axis_GPa, T_interp -def invert_to_entropy_along_profile(P_gpa, T_k, S_axis, T_of_SP): - """ - Invert T(S,P) -> S along a P-T profile. +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 :math:`T(S, P)`. + + Parameters + ---------- + P_gpa : np.ndarray + Pressure profile in GPa. + T_k : np.ndarray + Temperature profile in K. + S_axis : np.ndarray + Entropy grid used by the EOS table, in + :math:`\mathrm{J\,kg^{-1}\,K^{-1}}`. + T_of_SP : callable + Interpolator returning temperature for input points ``(S, P)``. + + Returns + ------- + np.ndarray + Entropy values corresponding to the input P–T profile. Values are set + to NaN where inversion is not possible. + + Method + ------ + For each point :math:`(P_i, T_i)` along the melting curve: + + 1. Evaluate the EOS table along the full entropy axis at fixed pressure + :math:`P_i`. + 2. Build a 1D relation :math:`T(S)` at that pressure. + 3. Invert that relation numerically using interpolation to obtain + :math:`S(T_i)`. + + Notes + ----- + Repeated temperature values are removed before interpolation. If linear + inversion fails, a nearest-neighbor fallback is attempted. """ S_out = np.full_like(T_k, np.nan, dtype=float) @@ -418,29 +1035,145 @@ def invert_to_entropy_along_profile(P_gpa, T_k, S_axis, T_of_SP): T_sorted = Tv[order] S_sorted = Sv[order] - if T_i < T_sorted[0] or T_i > T_sorted[-1]: + # Remove repeated temperature values before inversion. + T_unique, idx_unique = np.unique(T_sorted, return_index=True) + S_unique = S_sorted[idx_unique] + + if len(T_unique) < 2: + continue + + # Skip points that lie outside the invertible temperature range. + if T_i < T_unique[0] or T_i > T_unique[-1]: continue try: - f = interp1d(T_sorted, S_sorted, kind="linear", assume_sorted=True) + f = interp1d(T_unique, S_unique, kind="linear", assume_sorted=True) S_out[i] = float(f(T_i)) except Exception: - f = interp1d(T_sorted, S_sorted, kind="nearest", assume_sorted=True) - S_out[i] = float(f(T_i)) + 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 save_entropy_table_with_header(path: Path, P_gpa, S_jpk): +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. + + Parameters + ---------- + P_sol, S_sol : np.ndarray + Solidus pressure and entropy arrays. + P_liq, S_liq : np.ndarray + Liquidus pressure and entropy arrays. + n_common : int or None, optional + Number of points in the common pressure grid. If omitted, the smaller + valid curve length is used. + + Returns + ------- + P_common, S_sol_common, S_liq_common : tuple of np.ndarray + Shared pressure grid in GPa, and the solidus/liquidus entropy values + interpolated onto that grid. + + Notes + ----- + The common pressure range is defined as the overlap between the valid + solidus and liquidus pressure intervals: + + .. math:: + + P_{\min,\mathrm{common}} = \max(P_{\min,\mathrm{sol}}, P_{\min,\mathrm{liq}}) + + .. math:: + + P_{\max,\mathrm{common}} = \min(P_{\max,\mathrm{sol}}, P_{\max,\mathrm{liq}}) + + After sorting and removing duplicate pressures, both entropy curves are + linearly interpolated onto the shared pressure array. """ - Save entropy table in the scaled 2-column format. + 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): + r""" + Save a pressure–entropy table in SPIDER-style scaled format. + + Parameters + ---------- + path : pathlib.Path + Output filename. + P_gpa : np.ndarray + Pressure in GPa. + S_jpk : np.ndarray + Entropy in :math:`\mathrm{J\,kg^{-1}\,K^{-1}}`. + + Notes + ----- + The data are rescaled before saving so that the output matches the expected + SPIDER-style two-column format: + + - pressure divided by ``SCALE_P_OUT`` + - entropy divided by ``SCALE_S_OUT`` + + The header also stores the scaling factors, allowing the file to be read + back into SI units later. """ P_pa = P_gpa * 1e9 data = np.column_stack([P_pa / SCALE_P_OUT, S_jpk / SCALE_S_OUT]) np.savetxt(path, data, fmt="%.18e %.18e", header=COMMON_HEADER, comments="") -# Load EOS interpolators once +# Load EOS interpolators once at import time so they can be reused efficiently. S_axis_solid, P_axis_solid, T_of_SP_solid = load_eos_T_of_SP( eos_solid_path, nS_solid, SCALE_S_SOLID_EOS ) @@ -450,21 +1183,46 @@ def save_entropy_table_with_header(path: Path, P_gpa, S_jpk): ) -# ========================================================= -# Main exporter: save P-T and P-S -# ========================================================= -def export_model_curves(model_name, out_root="outputs_entropy_curves", - Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): - """ - Generate melting curves for one model, save: - - solidus.dat - - liquidus.dat - - solidus_entropy.dat - - liquidus_entropy.dat +# ============================================================================= +# MAIN EXPORTER +# ============================================================================= + +def export_model_curves(model_name: str, out_root: str = "outputs_entropy_curves", + Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 2000, **kwargs): + r""" + Export one melting model in both P–T and P–S space. + + Parameters + ---------- + model_name : str + Name of the melting parametrization. + out_root : str, optional + Root output directory. + Pmin, Pmax : float, optional + Pressure range in GPa. + n : int, optional + Number of points along the pressure grid. + **kwargs + Additional keyword arguments passed to the selected model. Returns ------- - dict with arrays + dict + Dictionary containing the raw and converted curves. + + Files written + ------------- + The following files are created inside ``out_root/model_name``: + + - ``solidus_P-T.dat`` + - ``liquidus_P-T.dat`` + - ``solidus_P-S.dat`` + - ``liquidus_P-S.dat`` + + Notes + ----- + The exported P–S files are resampled onto a single common pressure grid so + that both branches have identical length and aligned pressure sampling. """ out_dir = Path(out_root) / model_name out_dir.mkdir(parents=True, exist_ok=True) @@ -473,11 +1231,11 @@ def export_model_curves(model_name, out_root="outputs_entropy_curves", model_name, Pmin=Pmin, Pmax=Pmax, n=n, **kwargs ) - # Save P-T - save_PT_table(out_dir / "solidus.dat", P_sol, T_sol) - save_PT_table(out_dir / "liquidus.dat", P_liq, T_liq) + # Save the direct pressure–temperature curves. + 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) - # Convert to entropy + # Convert both branches from P–T to P–S using the EOS tables. S_sol = invert_to_entropy_along_profile( P_sol, T_sol, S_axis_solid, T_of_SP_solid ) @@ -485,20 +1243,22 @@ def export_model_curves(model_name, out_root="outputs_entropy_curves", P_liq, T_liq, S_axis_liquid, T_of_SP_liquid ) - # Remove NaNs before saving P-S files - mask_sol = np.isfinite(S_sol) - mask_liq = np.isfinite(S_liq) + # Place both entropy curves on the same pressure grid. + P_common, S_sol_common, S_liq_common = build_common_entropy_grid( + P_sol, S_sol, P_liq, S_liq, n_common=n + ) + # Save the entropy-space curves. save_entropy_table_with_header( - out_dir / "solidus_entropy.dat", - P_sol[mask_sol], - S_sol[mask_sol], + out_dir / "solidus_P-S.dat", + P_common, + S_sol_common, ) save_entropy_table_with_header( - out_dir / "liquidus_entropy.dat", - P_liq[mask_liq], - S_liq[mask_liq], + out_dir / "liquidus_P-S.dat", + P_common, + S_liq_common, ) return { @@ -508,13 +1268,36 @@ def export_model_curves(model_name, out_root="outputs_entropy_curves", "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, } -# ========================================================= -# Convenience: export all models -# ========================================================= -def export_all_models(out_root="outputs_entropy_curves"): +# ============================================================================= +# BATCH EXPORTER +# ============================================================================= + +def export_all_models(out_root: str = "outputs_entropy_curves"): + r""" + Export all supported melting parametrizations. + + Parameters + ---------- + out_root : str, optional + Root output directory where subdirectories for each model will be made. + + Notes + ----- + Some models require custom keyword arguments or pressure limits: + + - ``katz_2003`` uses ``X_h2o=30.0`` + - ``lin_2024`` uses ``fO2=-4.0`` + - ``hirschmann_2000`` is truncated at ``Pmax=10.0`` + - ``fei_2021`` and ``stixrude_2014`` start from ``Pmin=1.0`` + + The function prints progress messages to standard output. + """ models = [ "andrault_2011", "monteux_2016", @@ -547,95 +1330,18 @@ def export_all_models(out_root="outputs_entropy_curves"): print("[INFO] Done.") -# ========================================================= -# Optional quick plot -# ========================================================= -def show_melting_curves(model_name, Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): - P_sol, T_sol, P_liq, T_liq = get_melting_curves( - model_name, Pmin=Pmin, Pmax=Pmax, n=n, **kwargs - ) - - plt.figure(figsize=(6, 6)) - plt.plot(T_sol, P_sol, label="Solidus") - plt.plot(T_liq, P_liq, label="Liquidus") - plt.xlabel("Temperature (K)") - plt.ylabel("Pressure (GPa)") - plt.gca().invert_yaxis() - plt.title(model_name) - plt.legend() - plt.tight_layout() - plt.show() - - return P_sol, T_sol, P_liq, T_liq - -def show_entropy_curves(model_name, Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): - """ - Plot melting curves in P–S space. - """ - - P_sol, T_sol, P_liq, T_liq = get_melting_curves( - model_name, Pmin=Pmin, Pmax=Pmax, n=n, **kwargs - ) - - # convert to entropy - 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 - ) - - # remove NaNs before plotting - mask_sol = np.isfinite(S_sol) - mask_liq = np.isfinite(S_liq) +# ============================================================================= +# EXAMPLE USAGE +# ============================================================================= +# +# To export all models: +# +# export_all_models() +# +# To export a single model: +# +# export_model_curves("wolf_bower_2018") +# +# This script currently runs the full export when executed directly. - plt.figure(figsize=(6,6)) - - plt.plot(S_sol[mask_sol], P_sol[mask_sol], label="Solidus") - plt.plot(S_liq[mask_liq], P_liq[mask_liq], label="Liquidus") - - plt.xlabel("Entropy (J kg$^{-1}$ K$^{-1}$)") - plt.ylabel("Pressure (GPa)") - plt.gca().invert_yaxis() - - plt.title(f"{model_name} (P–S)") - plt.legend() - plt.tight_layout() - plt.show() - - return ( - P_sol[mask_sol], S_sol[mask_sol], - P_liq[mask_liq], S_liq[mask_liq], - ) - -# ========================================================= -# Example usage -# ========================================================= -# One model: -# export_model_curves("andrault_2011") -# export_model_curves("katz_2003", X_h2o=30.0) -# export_model_curves("lin_2024", fO2=-4.0) - -# All models: export_all_models() - -# Quick plot: -# show_melting_curves("wolf_bower_2018") -# show_melting_curves("andrault_2011") -# show_melting_curves("katz_2003", X_h2o=30.0) -# show_melting_curves("lin_2024", fO2=-4.0) -# show_melting_curves("fei_2021") -# show_melting_curves("belonoshko_2005") -# show_melting_curves("fiquet_2010") -# show_melting_curves("hirschmann_2000", Pmax=10.0) -# show_melting_curves("stixrude_2014", Pmin=1.0) -show_entropy_curves("wolf_bower_2018") -show_entropy_curves("andrault_2011") -show_entropy_curves("katz_2003", X_h2o=30.0) -show_entropy_curves("lin_2024", fO2=-4.0) -show_entropy_curves("fei_2021") -show_entropy_curves("belonoshko_2005") -show_entropy_curves("fiquet_2010") -show_entropy_curves("hirschmann_2000", Pmax=10.0) -show_entropy_curves("stixrude_2014", Pmin=1.0) From a4fa57cd370fa671137e36512d99932d130e37f9 Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Fri, 20 Mar 2026 12:54:49 +0100 Subject: [PATCH 03/13] Final script executable --- input/ensembles/example.grid.toml | 6 +- .../{interior => utils}/solidus_func.py | 1119 ++++++----------- 2 files changed, 391 insertions(+), 734 deletions(-) rename src/proteus/{interior => utils}/solidus_func.py (50%) diff --git a/input/ensembles/example.grid.toml b/input/ensembles/example.grid.toml index 157b82a58..7f2e79105 100644 --- a/input/ensembles/example.grid.toml +++ b/input/ensembles/example.grid.toml @@ -10,11 +10,11 @@ symlink = "" ref_config = "input/demos/dummy.toml" # Use SLURM? -use_slurm = false +use_slurm = true # Execution limits -max_jobs = 10 # maximum number of concurrent tasks (e.g. 500 on Habrok) -max_days = 1 # maximum number of days to run (e.g. 1) +max_jobs = 500 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 3 # maximum number of days to run (e.g. 1) max_mem = 3 # maximum memory per CPU in GB (e.g. 3) # Now define grid axes... diff --git a/src/proteus/interior/solidus_func.py b/src/proteus/utils/solidus_func.py similarity index 50% rename from src/proteus/interior/solidus_func.py rename to src/proteus/utils/solidus_func.py index 31054240a..574903ae8 100644 --- a/src/proteus/interior/solidus_func.py +++ b/src/proteus/utils/solidus_func.py @@ -17,6 +17,21 @@ 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, @@ -38,13 +53,28 @@ - Belonoshko et al. (2005), DOI: 10.1103/PhysRevLett.94.195701 - Fiquet et al. (2010), DOI: 10.1126/science.1192448 """ + from __future__ import annotations +import argparse +import os from pathlib import Path import numpy as np +import platformdirs from scipy.interpolate import RegularGridInterpolator, interp1d +# ============================================================================= +# 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 # ============================================================================= @@ -66,11 +96,6 @@ def make_pressure_grid(Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500) -> ------- np.ndarray One-dimensional pressure array in GPa. - - Notes - ----- - This helper is used by all melting parametrizations so that they share - a consistent pressure sampling unless otherwise specified. """ return np.linspace(Pmin, Pmax, n) @@ -78,31 +103,6 @@ def make_pressure_grid(Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500) -> def solidus_from_liquidus_stixrude(T_liq: np.ndarray) -> np.ndarray: r""" Estimate the solidus from a liquidus using the Stixrude ratio. - - Parameters - ---------- - T_liq : np.ndarray - Liquidus temperature in K. - - Returns - ------- - np.ndarray - Estimated solidus temperature in K. - - Formula - ------- - .. math:: - - T_{\mathrm{sol}} = \frac{T_{\mathrm{liq}}}{1 - \ln(0.79)} - - Notes - ----- - This is used in cases where only a liquidus parametrization is available - and an approximate solidus is needed. - - Reference - --------- - Stixrude (2014), DOI: 10.1098/rsta.2013.0076 """ return T_liq / (1.0 - np.log(0.79)) @@ -110,33 +110,59 @@ def solidus_from_liquidus_stixrude(T_liq: np.ndarray) -> np.ndarray: 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)) - Parameters - ---------- - T_sol : np.ndarray - Solidus temperature in K. - - Returns - ------- - np.ndarray - Estimated liquidus temperature in K. - - Formula - ------- - .. math:: - - T_{\mathrm{liq}} = T_{\mathrm{sol}} \left(1 - \ln(0.79)\right) - - Notes - ----- - This is used in cases where only a solidus parametrization is available - and an approximate liquidus is needed. - Reference - --------- - Stixrude (2014), DOI: 10.1098/rsta.2013.0076 +def _fmt_range(arr: np.ndarray, fmt: str = ".3f") -> str: """ - return T_sol * (1.0 - np.log(0.79)) + 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, +): + """ + Print a compact summary of the exported P-T and P-S ranges for one model. + """ + label = DISPLAY_NAMES.get(model_name, model_name) + print(f"{label}:") + print(f" P-T solidus : P_range = {_fmt_range(P_sol)} GPa, T_range = {_fmt_range(T_sol)} K") + print(f" P-T liquidus : P_range = {_fmt_range(P_liq)} GPa, T_range = {_fmt_range(T_liq)} K") + print(f" P-S solidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_sol_common)} J kg^-1 K^-1") + print(f" P-S liquidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_liq_common)} J kg^-1 K^-1") + print() # ============================================================================= @@ -146,41 +172,9 @@ def liquidus_from_solidus_stixrude(T_sol: np.ndarray) -> np.ndarray: 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). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - The parametrization is written as - - .. math:: - - T(P) = T_0 \left( \frac{P}{a} + 1 \right)^{1/c} - - where the published coefficients are defined using pressure in Pa. - In this implementation the user-facing pressure grid is generated in GPa - and internally converted to Pa. - - Reference - --------- - Andrault et al. (2011), DOI: 10.1016/j.epsl.2011.02.006 """ - P = make_pressure_grid(Pmin, Pmax, n) # GPa - P_pa = P * 1e9 # Pa + P = make_pressure_grid(Pmin, Pmax, n) + P_pa = P * 1e9 if kind == "solidus": T0, a, c = 2045, 92e9, 1.3 @@ -196,39 +190,6 @@ def andrault_2011(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000.0 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). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - The liquidus is given by - - .. math:: - - T_{\mathrm{liq}}(P) = 6000 \left(\frac{P}{140}\right)^{0.26} - - with pressure in GPa. - - If ``kind="solidus"``, the solidus is estimated from the liquidus using - :func:`solidus_from_liquidus_stixrude`. - - Reference - --------- - Fei et al. (2021), DOI: 10.1038/s41467-021-21170-y """ P = make_pressure_grid(Pmin, Pmax, n) T_liq = 6000.0 * (P / 140.0) ** 0.26 @@ -246,39 +207,6 @@ def fei_2021(kind: str = "liquidus", Pmin: float = 1.0, Pmax: float = 1000.0, n: 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). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - The liquidus is given by - - .. math:: - - T_{\mathrm{liq}}(P) = 1831 \left(1 + \frac{P}{4.6}\right)^{0.33} - - with pressure in GPa. - - If ``kind="solidus"``, the solidus is estimated from the liquidus using - :func:`solidus_from_liquidus_stixrude`. - - Reference - --------- - Belonoshko et al. (2005), DOI: 10.1103/PhysRevLett.94.195701 """ P = make_pressure_grid(Pmin, Pmax, n) T_liq = 1831.0 * (1.0 + P / 4.6) ** 0.33 @@ -296,47 +224,6 @@ def belonoshko_2005(kind: str = "liquidus", Pmin: float = 0.0, Pmax: float = 100 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). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - The liquidus is implemented as a two-branch fit: - - for :math:`P \leq 20\ \mathrm{GPa}`: - - .. math:: - - T_{\mathrm{liq}} = 1982.1 \left(\frac{P}{6.594} + 1\right)^{1/5.374} - - and for :math:`P > 20\ \mathrm{GPa}`: - - .. math:: - - T_{\mathrm{liq}} = 78.74 \left(\frac{P}{0.004056} + 1\right)^{1/2.44} - - where pressure is in GPa. - - If ``kind="solidus"``, the solidus is estimated from the liquidus using - :func:`solidus_from_liquidus_stixrude`. - - Reference - --------- - Fiquet et al. (2010), DOI: 10.1126/science.1192448 """ P = make_pressure_grid(Pmin, Pmax, n) T_liq = np.zeros_like(P, dtype=float) @@ -344,6 +231,8 @@ def fiquet_2010(kind: str = "liquidus", Pmin: float = 0.0, Pmax: float = 1000.0, 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.004054 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) @@ -360,41 +249,9 @@ def fiquet_2010(kind: str = "liquidus", Pmin: float = 0.0, Pmax: float = 1000.0, 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). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - Both solidus and liquidus are implemented as piecewise power-law fits: - - .. math:: - - T(P) = T_0 \left(\frac{P}{a} + 1\right)^{1/c} - - with one coefficient set below 20 GPa and another above 20 GPa. - The published coefficients use pressure in Pa, so pressure is converted - internally from GPa to Pa before evaluation. - - Reference - --------- - Monteux et al. (2016), DOI: 10.1016/j.epsl.2016.05.010 """ - P = make_pressure_grid(Pmin, Pmax, n) # GPa - P_pa = P * 1e9 # Pa + P = make_pressure_grid(Pmin, Pmax, n) + P_pa = P * 1e9 params = { "solidus": { @@ -425,42 +282,6 @@ def monteux_2016(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000.0, def hirschmann_2000(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 10.0, n: int = 500): r""" Melting curve from Hirschmann (2000). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - The fit is a quadratic polynomial in pressure: - - .. math:: - - T_{\mathrm{^\circ C}}(P) = A_1 + A_2 P + A_3 P^2 - - The result is then converted to Kelvin via - - .. math:: - - T_{\mathrm{K}} = T_{\mathrm{^\circ C}} + 273.15 - - This parametrization is only intended for relatively low pressures. - - Reference - --------- - Hirschmann (2000), DOI: 10.1029/2000GC000070 """ P = make_pressure_grid(Pmin, Pmax, n) @@ -482,39 +303,6 @@ def hirschmann_2000(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 10.0 def stixrude_2014(kind: str = "liquidus", Pmin: float = 1.0, Pmax: float = 1000.0, n: int = 500): r""" Melting curve based on Stixrude (2014). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - The liquidus is given by - - .. math:: - - T_{\mathrm{liq}}(P) = 5400 \left(\frac{P}{140}\right)^{0.480} - - with pressure in GPa. - - If ``kind="solidus"``, the solidus is estimated from the liquidus using - :func:`solidus_from_liquidus_stixrude`. - - Reference - --------- - Stixrude (2014), DOI: 10.1098/rsta.2013.0076 """ P = make_pressure_grid(Pmin, Pmax, n) T_liq = 5400.0 * (P / 140.0) ** 0.480 @@ -532,54 +320,6 @@ def stixrude_2014(kind: str = "liquidus", Pmin: float = 1.0, Pmax: float = 1000. 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. - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - The curve is piecewise linear in pressure, with breakpoints - :math:`P_{\mathrm{cp1}}` and :math:`P_{\mathrm{cp2}}`: - - .. math:: - - T(P) = - \begin{cases} - s_1 P + c_1, & P < P_{\mathrm{cp1}} \\ - s_2 P + c_2, & P_{\mathrm{cp1}} \leq P < P_{\mathrm{cp2}} \\ - s_3 P + c_3, & P \geq P_{\mathrm{cp2}} - \end{cases} - - where continuity is enforced through - - .. math:: - - c_2 = c_1 + (s_1 - s_2) P_{\mathrm{cp1}} - - .. math:: - - c_3 = c_2 + (s_2 - s_3) P_{\mathrm{cp2}} - - Notes - ----- - This implementation uses coefficients already encoded in the script. - - Reference - --------- - Wolf & Bower (2018)-style parametrization used in the model workflow. """ P = make_pressure_grid(Pmin, Pmax, n) @@ -627,47 +367,6 @@ def wolf_bower_2018(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000 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). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - X_h2o : float - Water content parameter used in the hydrous temperature depression. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - Starting from the anhydrous Wolf & Bower curve, a constant depression is - applied at all pressures: - - .. math:: - - \Delta T = K X_{\mathrm{H_2O}}^{\gamma} - - .. math:: - - T = T_{\mathrm{anhydrous}} - \Delta T - - where in this implementation - - .. math:: - - \gamma = 0.75,\quad K = 43 - - Reference - --------- - Katz et al. (2003), DOI: 10.1029/2002GC000433 """ gamma = 0.75 K = 43.0 @@ -685,46 +384,6 @@ def katz_2003(kind: str = "solidus", X_h2o: float = 30.0, Pmin: float = 0.0, Pma 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). - - Parameters - ---------- - kind : {"solidus", "liquidus"} - Which branch to evaluate. - fO2 : float - Oxygen fugacity offset parameter used in the solidus shift. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of pressure samples. - - Returns - ------- - P : np.ndarray - Pressure in GPa. - T : np.ndarray - Temperature in K. - - Formula - ------- - The anhydrous solidus is first taken from :func:`wolf_bower_2018` and then - shifted by - - .. math:: - - \Delta T = \frac{340}{3.2} (2 - f\mathrm{O}_2) - - so that - - .. math:: - - T_{\mathrm{sol}} = T_{\mathrm{anhydrous}} + \Delta T - - If ``kind="liquidus"``, the liquidus is estimated from the shifted solidus - using :func:`liquidus_from_solidus_stixrude`. - - Reference - --------- - Lin et al. (2024), DOI: 10.1038/s41561-024-01495-1 """ P, T_anhydrous = wolf_bower_2018(kind="solidus", Pmin=Pmin, Pmax=Pmax, n=n) @@ -748,31 +407,7 @@ def lin_2024(kind: str = "solidus", fO2: float = -4.0, Pmin: float = 0.0, Pmax: def truncate_to_physical_interval(func): r""" Wrap a melting-curve function so only the main interval with - :math:`T_{\mathrm{sol}} < T_{\mathrm{liq}}` is retained. - - Parameters - ---------- - func : callable - A function returning a melting curve in the form ``P, T`` for - ``kind="solidus"`` or ``kind="liquidus"``. - - Returns - ------- - callable - Wrapped function returning only the largest contiguous physically valid - pressure interval. - - Notes - ----- - Some parametrizations can cross or become unphysical at high pressure. - This wrapper evaluates both solidus and liquidus, identifies all points - where - - .. math:: - - T_{\mathrm{sol}} < T_{\mathrm{liq}} - - and keeps only the largest contiguous block of such points. + 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) @@ -790,15 +425,13 @@ def wrapped(kind="solidus", Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): if kind == "solidus": return P_sol[main_block], T_sol[main_block] - elif kind == "liquidus": + if kind == "liquidus": return P_liq[main_block], T_liq[main_block] - else: - raise ValueError("kind must be 'solidus' or 'liquidus'") + raise ValueError("kind must be 'solidus' or 'liquidus'") return wrapped -# Wrapped versions for models where physical truncation is needed. andrault_2011_cut = truncate_to_physical_interval(andrault_2011) monteux_2016_cut = truncate_to_physical_interval(monteux_2016) wolf_bower_2018_cut = truncate_to_physical_interval(wolf_bower_2018) @@ -809,31 +442,23 @@ def wrapped(kind="solidus", Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): # MODEL DISPATCHER # ============================================================================= +SUPPORTED_MODELS = [ + "andrault_2011", + "monteux_2016", + "wolf_bower_2018", + "katz_2003", + "fei_2021", + "belonoshko_2005", + "fiquet_2010", + "hirschmann_2000", + "stixrude_2014", + "lin_2024", +] + + 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. - - Parameters - ---------- - model_name : str - Identifier of the melting parametrization. - Pmin, Pmax : float - Pressure range in GPa. - n : int - Number of sampling points. - **kwargs - Additional keyword arguments passed to the selected model - (for example ``X_h2o`` or ``fO2``). - - Returns - ------- - P_sol, T_sol, P_liq, T_liq : tuple of np.ndarray - Solidus and liquidus curves in P–T space. - - Raises - ------ - ValueError - If the requested model name is unknown. """ models = { "andrault_2011": andrault_2011_cut, @@ -864,26 +489,7 @@ def get_melting_curves(model_name: str, Pmin: float = 0.0, Pmax: float = 1000.0, def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): r""" - Save a pressure–temperature table to disk. - - Parameters - ---------- - path : pathlib.Path - Output filename. - P_gpa : np.ndarray - Pressure in GPa. - T_k : np.ndarray - Temperature in K. - - Notes - ----- - The file is written as two columns: - - .. code-block:: text - - pressure temperature - - using a simple header compatible with later inspection. + 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="#") @@ -892,38 +498,32 @@ def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): # ============================================================================= # EOS LOOKUP TABLE SETTINGS # ============================================================================= -# -# The files temperature_solid.dat and temperature_melt.dat contain structured -# tables in (S, P) space. Each row stores: -# -# pressure, entropy, temperature -# -# in scaled units. The constants below decode those tables back to SI units. - -eos_solid_path = Path("temperature_solid.dat") -eos_liquid_path = Path("temperature_melt.dat") - -# Number of pressure points in the EOS tables. -nP = 2020 -# Number of entropy points for the solid and liquid tables, respectively. +SCRIPT_DIR = Path(__file__).resolve().parent +spider_dir = (SCRIPT_DIR / "../../../SPIDER").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}") + +nP = 2020 nS_solid = 125 nS_liquid = 95 - -# Number of header lines to skip when reading the raw EOS files. skip_header = 5 -# Scaling factors used by the EOS tables. SCALE_P_EOS = 1e9 SCALE_T_EOS = 1.0 SCALE_S_SOLID_EOS = 4.82426684604467e6 SCALE_S_LIQUID_EOS = 4.805046659407042e6 -# Scaling factors used when exporting P–S tables back to SPIDER-like format. SCALE_P_OUT = 1_000_000_000.0 SCALE_S_OUT = 4_824_266.84604467 -# Header used for exported SPIDER-style entropy tables. COMMON_HEADER = "\n".join([ "# 5 400", "# Pressure, Entropy, Quantity", @@ -935,35 +535,7 @@ def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): r""" - Load an EOS lookup table and build an interpolator for :math:`T(S, P)`. - - Parameters - ---------- - eos_path : pathlib.Path - Path to the EOS file. - nS : int - Number of entropy grid points in the table. - scale_S_axis : float - Scaling factor needed to convert the stored entropy column to SI units. - - Returns - ------- - S_axis : np.ndarray - Entropy axis in :math:`\mathrm{J\,kg^{-1}\,K^{-1}}`. - P_axis_GPa : np.ndarray - Pressure axis in GPa. - T_interp : scipy.interpolate.RegularGridInterpolator - Interpolator returning temperature in K for points in ``(S, P)``. - - Notes - ----- - The EOS files are stored in a flattened 3-column format and reshaped into - - .. math:: - - (n_S, n_P, 3) - - where the three columns correspond to pressure, entropy, and temperature. + Load an EOS lookup table and build an interpolator for T(S, P). """ raw = np.genfromtxt(eos_path, skip_header=skip_header) resh = raw.reshape((nS, nP, 3)) @@ -983,40 +555,7 @@ def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): 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 :math:`T(S, P)`. - - Parameters - ---------- - P_gpa : np.ndarray - Pressure profile in GPa. - T_k : np.ndarray - Temperature profile in K. - S_axis : np.ndarray - Entropy grid used by the EOS table, in - :math:`\mathrm{J\,kg^{-1}\,K^{-1}}`. - T_of_SP : callable - Interpolator returning temperature for input points ``(S, P)``. - - Returns - ------- - np.ndarray - Entropy values corresponding to the input P–T profile. Values are set - to NaN where inversion is not possible. - - Method - ------ - For each point :math:`(P_i, T_i)` along the melting curve: - - 1. Evaluate the EOS table along the full entropy axis at fixed pressure - :math:`P_i`. - 2. Build a 1D relation :math:`T(S)` at that pressure. - 3. Invert that relation numerically using interpolation to obtain - :math:`S(T_i)`. - - Notes - ----- - Repeated temperature values are removed before interpolation. If linear - inversion fails, a nearest-neighbor fallback is attempted. + 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) @@ -1035,14 +574,12 @@ def invert_to_entropy_along_profile(P_gpa: np.ndarray, T_k: np.ndarray, S_axis: T_sorted = Tv[order] S_sorted = Sv[order] - # Remove repeated temperature values before inversion. T_unique, idx_unique = np.unique(T_sorted, return_index=True) S_unique = S_sorted[idx_unique] if len(T_unique) < 2: continue - # Skip points that lie outside the invertible temperature range. if T_i < T_unique[0] or T_i > T_unique[-1]: continue @@ -1059,41 +596,15 @@ def invert_to_entropy_along_profile(P_gpa: np.ndarray, T_k: np.ndarray, S_axis: 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): +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. - - Parameters - ---------- - P_sol, S_sol : np.ndarray - Solidus pressure and entropy arrays. - P_liq, S_liq : np.ndarray - Liquidus pressure and entropy arrays. - n_common : int or None, optional - Number of points in the common pressure grid. If omitted, the smaller - valid curve length is used. - - Returns - ------- - P_common, S_sol_common, S_liq_common : tuple of np.ndarray - Shared pressure grid in GPa, and the solidus/liquidus entropy values - interpolated onto that grid. - - Notes - ----- - The common pressure range is defined as the overlap between the valid - solidus and liquidus pressure intervals: - - .. math:: - - P_{\min,\mathrm{common}} = \max(P_{\min,\mathrm{sol}}, P_{\min,\mathrm{liq}}) - - .. math:: - - P_{\max,\mathrm{common}} = \min(P_{\max,\mathrm{sol}}, P_{\max,\mathrm{liq}}) - - After sorting and removing duplicate pressures, both entropy curves are - linearly interpolated onto the shared pressure array. """ mask_sol = np.isfinite(S_sol) mask_liq = np.isfinite(S_liq) @@ -1146,38 +657,16 @@ def build_common_entropy_grid(P_sol: np.ndarray, S_sol: np.ndarray, P_liq: np.nd def save_entropy_table_with_header(path: Path, P_gpa: np.ndarray, S_jpk: np.ndarray): r""" - Save a pressure–entropy table in SPIDER-style scaled format. - - Parameters - ---------- - path : pathlib.Path - Output filename. - P_gpa : np.ndarray - Pressure in GPa. - S_jpk : np.ndarray - Entropy in :math:`\mathrm{J\,kg^{-1}\,K^{-1}}`. - - Notes - ----- - The data are rescaled before saving so that the output matches the expected - SPIDER-style two-column format: - - - pressure divided by ``SCALE_P_OUT`` - - entropy divided by ``SCALE_S_OUT`` - - The header also stores the scaling factors, allowing the file to be read - back into SI units later. + Save a pressure-entropy table in SPIDER-style scaled format. """ P_pa = P_gpa * 1e9 data = np.column_stack([P_pa / SCALE_P_OUT, S_jpk / SCALE_S_OUT]) np.savetxt(path, data, fmt="%.18e %.18e", header=COMMON_HEADER, comments="") -# Load EOS interpolators once at import time so they can be reused efficiently. S_axis_solid, P_axis_solid, T_of_SP_solid = load_eos_T_of_SP( eos_solid_path, nS_solid, SCALE_S_SOLID_EOS ) - S_axis_liquid, P_axis_liquid, T_of_SP_liquid = load_eos_T_of_SP( eos_liquid_path, nS_liquid, SCALE_S_LIQUID_EOS ) @@ -1187,42 +676,16 @@ def save_entropy_table_with_header(path: Path, P_gpa: np.ndarray, S_jpk: np.ndar # MAIN EXPORTER # ============================================================================= -def export_model_curves(model_name: str, out_root: str = "outputs_entropy_curves", - Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 2000, **kwargs): +def export_model_curves( + model_name: str, + out_root: Path | str = MELTING_DIR, + Pmin: float = 0.0, + Pmax: float = 1000.0, + n: int = 2000, + **kwargs, +): r""" - Export one melting model in both P–T and P–S space. - - Parameters - ---------- - model_name : str - Name of the melting parametrization. - out_root : str, optional - Root output directory. - Pmin, Pmax : float, optional - Pressure range in GPa. - n : int, optional - Number of points along the pressure grid. - **kwargs - Additional keyword arguments passed to the selected model. - - Returns - ------- - dict - Dictionary containing the raw and converted curves. - - Files written - ------------- - The following files are created inside ``out_root/model_name``: - - - ``solidus_P-T.dat`` - - ``liquidus_P-T.dat`` - - ``solidus_P-S.dat`` - - ``liquidus_P-S.dat`` - - Notes - ----- - The exported P–S files are resampled onto a single common pressure grid so - that both branches have identical length and aligned pressure sampling. + 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) @@ -1231,11 +694,9 @@ def export_model_curves(model_name: str, out_root: str = "outputs_entropy_curves model_name, Pmin=Pmin, Pmax=Pmax, n=n, **kwargs ) - # Save the direct pressure–temperature curves. 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) - # Convert both branches from P–T to P–S using the EOS tables. S_sol = invert_to_entropy_along_profile( P_sol, T_sol, S_axis_solid, T_of_SP_solid ) @@ -1243,12 +704,10 @@ def export_model_curves(model_name: str, out_root: str = "outputs_entropy_curves P_liq, T_liq, S_axis_liquid, T_of_SP_liquid ) - # Place both entropy curves on the same pressure grid. P_common, S_sol_common, S_liq_common = build_common_entropy_grid( P_sol, S_sol, P_liq, S_liq, n_common=n ) - # Save the entropy-space curves. save_entropy_table_with_header( out_dir / "solidus_P-S.dat", P_common, @@ -1261,6 +720,16 @@ def export_model_curves(model_name: str, out_root: str = "outputs_entropy_curves S_liq_common, ) + print_model_summary( + model_name, + P_sol, T_sol, + P_liq, T_liq, + P_common, S_sol_common, S_liq_common, + ) + + print(f" Saved to : {out_dir.resolve()}") + print() + return { "P_sol": P_sol, "T_sol": T_sol, @@ -1278,70 +747,258 @@ def export_model_curves(model_name: str, out_root: str = "outputs_entropy_curves # BATCH EXPORTER # ============================================================================= -def export_all_models(out_root: str = "outputs_entropy_curves"): +def export_all_models(out_root: Path | str = MELTING_DIR, n: int = 2000): r""" Export all supported melting parametrizations. - - Parameters - ---------- - out_root : str, optional - Root output directory where subdirectories for each model will be made. - - Notes - ----- - Some models require custom keyword arguments or pressure limits: - - - ``katz_2003`` uses ``X_h2o=30.0`` - - ``lin_2024`` uses ``fO2=-4.0`` - - ``hirschmann_2000`` is truncated at ``Pmax=10.0`` - - ``fei_2021`` and ``stixrude_2014`` start from ``Pmin=1.0`` - - The function prints progress messages to standard output. """ - models = [ - "andrault_2011", - "monteux_2016", - "wolf_bower_2018", - "katz_2003", - "fei_2021", - "belonoshko_2005", - "fiquet_2010", - "hirschmann_2000", - "stixrude_2014", - "lin_2024", - ] - - for model in models: - print(f"[INFO] Exporting {model}") - + for model in SUPPORTED_MODELS: if model == "katz_2003": - export_model_curves(model, out_root=out_root, X_h2o=30.0) + _ = export_model_curves(model, out_root=out_root, n=n, X_h2o=30.0) elif model == "lin_2024": - export_model_curves(model, out_root=out_root, fO2=-4.0) + _ = export_model_curves(model, out_root=out_root, n=n, fO2=-4.0) elif model == "hirschmann_2000": - export_model_curves(model, out_root=out_root, Pmax=10.0) + _ = export_model_curves(model, out_root=out_root, n=n, Pmax=10.0) elif model == "fei_2021": - export_model_curves(model, out_root=out_root, Pmin=1.0) + _ = export_model_curves(model, out_root=out_root, n=n, Pmin=1.0) elif model == "stixrude_2014": - export_model_curves(model, out_root=out_root, Pmin=1.0) + _ = export_model_curves(model, out_root=out_root, n=n, Pmin=1.0) else: - export_model_curves(model, out_root=out_root) - - print("[INFO] Done.") + _ = export_model_curves(model, out_root=out_root, n=n) # ============================================================================= -# EXAMPLE USAGE +# COMMAND-LINE INTERFACE # ============================================================================= -# -# To export all models: -# -# export_all_models() -# -# To export a single model: -# -# export_model_curves("wolf_bower_2018") -# -# This script currently runs the full export when executed directly. - -export_all_models() + +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( + "--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, + **kwargs, + ) + + +def main(): + 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." + ) + export_all_models(out_root=args.out_root, n=args.n) + return + + selected_models = [m for m in [explicit_model, shortcut_model] if m is not None] + + if len(selected_models) == 0: + raise SystemExit( + "Error: no model selected. Use --all or choose one model with " + "--model or a shortcut like --katz2003." + ) + + if len(selected_models) > 1: + raise SystemExit( + "Error: please choose only one of --model or one shortcut flag." + ) + + export_one_model_from_cli(selected_models[0], args) + + +if __name__ == "__main__": + main() From 4b8ca603360848f1d8b1457368ca8f00bb38b2fb Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Fri, 20 Mar 2026 13:38:09 +0100 Subject: [PATCH 04/13] Modify the donwload routine --- input/all_options.toml | 4 +- src/proteus/utils/data.py | 70 ++++++++++++++++++++++--------- src/proteus/utils/solidus_func.py | 2 +- 3 files changed, 53 insertions(+), 23 deletions(-) diff --git a/input/all_options.toml b/input/all_options.toml index bcf054c36..d142b258e 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%) @@ -311,7 +311,7 @@ version = "2.0" tidal_heat = false # enable tidal heat production rheo_phi_loc = 0.6 # Centre of rheological transition rheo_phi_wid = 0.2 # Width of rheological transition - melting_dir = "Monteux-600" # Melting curve set used by all interior modules (Zalmoxis, Aragog, SPIDER) + melting_dir = "andrault_2011" # Melting curve set used by all interior modules (Zalmoxis, Aragog, SPIDER) eos_dir = "WolfBower2018_MgSiO3" # Dynamic EOS for SPIDER + Aragog (in FWL_DATA/interior_lookup_tables/EOS/dynamic/) diff --git a/src/proteus/utils/data.py b/src/proteus/utils/data.py index 60a715582..127c2129f 100644 --- a/src/proteus/utils/data.py +++ b/src/proteus/utils/data.py @@ -1189,41 +1189,71 @@ def download_interior_lookuptables(clean=False): desc=f'Interior lookup tables: {dir}', ) - -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 + log.debug("Download melting curve data") + rel_dir = Path("Melting_curves") / config.interior.melting_dir - data_dir = GetFWLData() / 'interior_lookup_tables' + 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}', + osf_id=source_info["osf_project"], + zenodo_id=source_info["zenodo_id"], + 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. - for stem in ('solidus', 'liquidus'): - legacy = folder_dir / f'{stem}.dat' - canonical = folder_dir / f'{stem}_P-T.dat' + # ------------------------------------------------------------------ + # 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" if legacy.is_file() and not canonical.is_file(): shutil.copy2(legacy, canonical) - log.debug('Copied %s -> %s', legacy.name, canonical.name) - + log.debug("Copied %s -> %s", legacy.name, canonical.name) def download_stellar_spectra(*, folders: tuple[str, ...] | None = None): """ diff --git a/src/proteus/utils/solidus_func.py b/src/proteus/utils/solidus_func.py index 574903ae8..9d70c0359 100644 --- a/src/proteus/utils/solidus_func.py +++ b/src/proteus/utils/solidus_func.py @@ -525,7 +525,7 @@ def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): SCALE_S_OUT = 4_824_266.84604467 COMMON_HEADER = "\n".join([ - "# 5 400", + "# 5 2000", "# Pressure, Entropy, Quantity", "# column * scaling factor should be SI units", "# scaling factors (constant) for each column given on line below", From a5e7fdeb777dcb5040faca263b14878e9ebe4817 Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Sun, 22 Mar 2026 18:43:40 +0100 Subject: [PATCH 05/13] Add description in usage --- docs/How-to/usage.md | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/docs/How-to/usage.md b/docs/How-to/usage.md index c8a79bb16..dfd41d78d 100644 --- a/docs/How-to/usage.md +++ b/docs/How-to/usage.md @@ -296,3 +296,36 @@ 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 define the temperatures at which materials begin to melt and become fully molten as a function of pressure. + +### Available parametrizations + +- Andrault et al. (2011) +- Monteux et al. (2016) +- Wolf & Bower (2018) +- Katz et al. (2003) +- Lin et al. (2024) +- Hirschmann (2000) +- Stixrude (2014) +- Fei et al. (2021) +- Belonoshko et al. (2005) +- Fiquet et al. (2010) + +### Generate melting curves + +Before running PROTEUS, generate the lookup tables: + +```console +python src/proteus/utils/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/ +``` From 1945aa5653da45d09dec92186b5ad477ddfe0557 Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Sun, 22 Mar 2026 19:18:28 +0100 Subject: [PATCH 06/13] Set original metling_dir --- input/all_options.toml | 2 +- input/ensembles/example.grid.toml | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/input/all_options.toml b/input/all_options.toml index d142b258e..9eabab3cd 100644 --- a/input/all_options.toml +++ b/input/all_options.toml @@ -311,7 +311,7 @@ version = "2.0" tidal_heat = false # enable tidal heat production rheo_phi_loc = 0.6 # Centre of rheological transition rheo_phi_wid = 0.2 # Width of rheological transition - melting_dir = "andrault_2011" # Melting curve set used by all interior modules (Zalmoxis, Aragog, SPIDER) + melting_dir = "Monteux-600" # Melting curve set used by all interior modules (Zalmoxis, Aragog, SPIDER) eos_dir = "WolfBower2018_MgSiO3" # Dynamic EOS for SPIDER + Aragog (in FWL_DATA/interior_lookup_tables/EOS/dynamic/) diff --git a/input/ensembles/example.grid.toml b/input/ensembles/example.grid.toml index 7f2e79105..157b82a58 100644 --- a/input/ensembles/example.grid.toml +++ b/input/ensembles/example.grid.toml @@ -10,11 +10,11 @@ symlink = "" ref_config = "input/demos/dummy.toml" # Use SLURM? -use_slurm = true +use_slurm = false # Execution limits -max_jobs = 500 # maximum number of concurrent tasks (e.g. 500 on Habrok) -max_days = 3 # maximum number of days to run (e.g. 1) +max_jobs = 10 # maximum number of concurrent tasks (e.g. 500 on Habrok) +max_days = 1 # maximum number of days to run (e.g. 1) max_mem = 3 # maximum memory per CPU in GB (e.g. 3) # Now define grid axes... From 480cd0ef6077a375a33711c38f5388e8a63b0c4f Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Sun, 22 Mar 2026 19:22:54 +0100 Subject: [PATCH 07/13] Check ruff --- docs/How-to/usage.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/How-to/usage.md b/docs/How-to/usage.md index dfd41d78d..c5eb399ed 100644 --- a/docs/How-to/usage.md +++ b/docs/How-to/usage.md @@ -322,7 +322,7 @@ Before running PROTEUS, generate the lookup tables: python src/proteus/utils/solidus_func.py --all ``` -Alternatively, you can generate a single parametrization using a specific flag (e.g. --katz2003, --lin2024). +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: From 01db1c40a3c5a14fb4e0318590ac9d81aa026458 Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Sun, 22 Mar 2026 19:39:16 +0100 Subject: [PATCH 08/13] Fix lint/formatting --- src/proteus/utils/data.py | 36 ++- src/proteus/utils/solidus_func.py | 508 +++++++++++++++++------------- 2 files changed, 302 insertions(+), 242 deletions(-) diff --git a/src/proteus/utils/data.py b/src/proteus/utils/data.py index 127c2129f..2f1ae72de 100644 --- a/src/proteus/utils/data.py +++ b/src/proteus/utils/data.py @@ -1189,6 +1189,7 @@ def download_interior_lookuptables(clean=False): desc=f'Interior lookup tables: {dir}', ) + def download_melting_curves(config: Config, clean: bool = False): """ Ensure melting curve data are available locally. @@ -1202,10 +1203,10 @@ def download_melting_curves(config: Config, clean: bool = False): solidus_P-S.dat liquidus_P-S.dat """ - log.debug("Download melting curve data") - rel_dir = Path("Melting_curves") / config.interior.melting_dir + log.debug('Download melting curve data') + rel_dir = Path('Melting_curves') / config.interior.melting_dir - data_dir = GetFWLData() / "interior_lookup_tables" + data_dir = GetFWLData() / 'interior_lookup_tables' data_dir.mkdir(parents=True, exist_ok=True) folder_dir = data_dir / rel_dir @@ -1216,13 +1217,13 @@ def download_melting_curves(config: Config, clean: bool = False): # ------------------------------------------------------------------ # 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" + 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) + log.debug('Melting curve data already present locally: %s', folder_dir) return # ------------------------------------------------------------------ @@ -1231,16 +1232,16 @@ def download_melting_curves(config: Config, clean: bool = False): source_info = get_data_source_info(rel_dir.as_posix()) if not source_info: raise ValueError( - f"No data source mapping found for folder: {rel_dir}. " - f"Also did not find local melting curve data in: {folder_dir}" + f'No data source mapping found for folder: {rel_dir}. ' + f'Also did not find local melting curve data in: {folder_dir}' ) download( 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: {rel_dir}", + osf_id=source_info['osf_project'], + zenodo_id=source_info['zenodo_id'], + desc=f'Melting curve data: {rel_dir}', ) # ------------------------------------------------------------------ @@ -1248,12 +1249,13 @@ def download_melting_curves(config: Config, clean: bool = False): # - 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" + for stem in ('solidus', 'liquidus'): + legacy = folder_dir / f'{stem}.dat' + canonical = folder_dir / f'{stem}_P-T.dat' if legacy.is_file() and not canonical.is_file(): shutil.copy2(legacy, canonical) - log.debug("Copied %s -> %s", legacy.name, canonical.name) + log.debug('Copied %s -> %s', legacy.name, canonical.name) + def download_stellar_spectra(*, folders: tuple[str, ...] | None = None): """ diff --git a/src/proteus/utils/solidus_func.py b/src/proteus/utils/solidus_func.py index 9d70c0359..bb5bfebd0 100644 --- a/src/proteus/utils/solidus_func.py +++ b/src/proteus/utils/solidus_func.py @@ -68,17 +68,16 @@ # DEFAULT OUTPUT LOCATION # ============================================================================= -FWL_DATA_DIR = Path( - os.environ.get("FWL_DATA", platformdirs.user_data_dir("fwl_data")) -) +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" +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. @@ -114,7 +113,7 @@ def liquidus_from_solidus_stixrude(T_sol: np.ndarray) -> np.ndarray: return T_sol * (1.0 - np.log(0.79)) -def _fmt_range(arr: np.ndarray, fmt: str = ".3f") -> str: +def _fmt_range(arr: np.ndarray, fmt: str = '.3f') -> str: """ Format the finite minimum and maximum of an array as a string. """ @@ -122,24 +121,24 @@ def _fmt_range(arr: np.ndarray, fmt: str = ".3f") -> str: mask = np.isfinite(arr) if not np.any(mask): - return "[nan, nan]" + return '[nan, nan]' amin = np.min(arr[mask]) amax = np.max(arr[mask]) - return f"[{amin:{fmt}}, {amax:{fmt}}]" + 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)", + '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)', } @@ -157,11 +156,19 @@ def print_model_summary( Print a compact summary of the exported P-T and P-S ranges for one model. """ label = DISPLAY_NAMES.get(model_name, model_name) - print(f"{label}:") - print(f" P-T solidus : P_range = {_fmt_range(P_sol)} GPa, T_range = {_fmt_range(T_sol)} K") - print(f" P-T liquidus : P_range = {_fmt_range(P_liq)} GPa, T_range = {_fmt_range(T_liq)} K") - print(f" P-S solidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_sol_common)} J kg^-1 K^-1") - print(f" P-S liquidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_liq_common)} J kg^-1 K^-1") + print(f'{label}:') + print( + f' P-T solidus : P_range = {_fmt_range(P_sol)} GPa, T_range = {_fmt_range(T_sol)} K' + ) + print( + f' P-T liquidus : P_range = {_fmt_range(P_liq)} GPa, T_range = {_fmt_range(T_liq)} K' + ) + print( + f' P-S solidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_sol_common)} J kg^-1 K^-1' + ) + print( + f' P-S liquidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_liq_common)} J kg^-1 K^-1' + ) print() @@ -169,16 +176,17 @@ def print_model_summary( # LITERATURE MELTING CURVES # ============================================================================= -def andrault_2011(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500): + +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": + if kind == 'solidus': T0, a, c = 2045, 92e9, 1.3 - elif kind == "liquidus": + elif kind == 'liquidus': T0, a, c = 1940, 29e9, 1.9 else: raise ValueError("kind must be 'solidus' or 'liquidus'") @@ -187,16 +195,16 @@ def andrault_2011(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000.0 return P, T -def fei_2021(kind: str = "liquidus", Pmin: float = 1.0, Pmax: float = 1000.0, n: int = 500): +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": + if kind == 'liquidus': T = T_liq - elif kind == "solidus": + elif kind == 'solidus': T = solidus_from_liquidus_stixrude(T_liq) else: raise ValueError("kind must be 'solidus' or 'liquidus'") @@ -204,16 +212,18 @@ def fei_2021(kind: str = "liquidus", Pmin: float = 1.0, Pmax: float = 1000.0, n: return P, T -def belonoshko_2005(kind: str = "liquidus", Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500): +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": + if kind == 'liquidus': T = T_liq - elif kind == "solidus": + elif kind == 'solidus': T = solidus_from_liquidus_stixrude(T_liq) else: raise ValueError("kind must be 'solidus' or 'liquidus'") @@ -221,7 +231,7 @@ def belonoshko_2005(kind: str = "liquidus", Pmin: float = 0.0, Pmax: float = 100 return P, T -def fiquet_2010(kind: str = "liquidus", Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500): +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). """ @@ -236,9 +246,9 @@ def fiquet_2010(kind: str = "liquidus", Pmin: float = 0.0, Pmax: float = 1000.0, 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": + if kind == 'liquidus': T = T_liq - elif kind == "solidus": + elif kind == 'solidus': T = solidus_from_liquidus_stixrude(T_liq) else: raise ValueError("kind must be 'solidus' or 'liquidus'") @@ -246,7 +256,7 @@ def fiquet_2010(kind: str = "liquidus", Pmin: float = 0.0, Pmax: float = 1000.0, return P, T -def monteux_2016(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500): +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). """ @@ -254,14 +264,14 @@ def monteux_2016(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000.0, 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}, + '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}, }, - "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: @@ -273,21 +283,25 @@ def monteux_2016(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000.0, 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"]) + 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): +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), + 'solidus': (1085.7, 132.9, -5.1), + 'liquidus': (1475.0, 80.0, -3.2), } if kind not in coeffs: @@ -300,16 +314,18 @@ def hirschmann_2000(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 10.0 return P, T -def stixrude_2014(kind: str = "liquidus", Pmin: float = 1.0, Pmax: float = 1000.0, n: int = 500): +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": + if kind == 'liquidus': T = T_liq - elif kind == "solidus": + elif kind == 'solidus': T = solidus_from_liquidus_stixrude(T_liq) else: raise ValueError("kind must be 'solidus' or 'liquidus'") @@ -317,29 +333,31 @@ def stixrude_2014(kind: str = "liquidus", Pmin: float = 1.0, Pmax: float = 1000. return P, T -def wolf_bower_2018(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 500): +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": ( + 'solidus': ( 7.696777581585296, 870.4767697319186, 101.52655163737373, 15.959022187236807, 3.090844734784906, - 1417.4258954709148 + 1417.4258954709148, ), - "liquidus": ( + 'liquidus': ( 8.864665249317456, 408.58442302949794, 46.288444869816615, 17.549174419770257, 3.679647802112376, - 2019.967799687511 - ) + 2019.967799687511, + ), } if kind not in params: @@ -364,35 +382,47 @@ def wolf_bower_2018(kind: str = "solidus", Pmin: float = 0.0, Pmax: float = 1000 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): +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"}: + 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 + 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): +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) + 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": + if kind == 'solidus': T = T_sol - elif kind == "liquidus": + elif kind == 'liquidus': T = liquidus_from_solidus_stixrude(T_sol) else: raise ValueError("kind must be 'solidus' or 'liquidus'") @@ -404,28 +434,30 @@ def lin_2024(kind: str = "solidus", fO2: float = -4.0, Pmin: float = 0.0, Pmax: # 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) + + 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") + 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": + if kind == 'solidus': return P_sol[main_block], T_sol[main_block] - if kind == "liquidus": + if kind == 'liquidus': return P_liq[main_block], T_liq[main_block] raise ValueError("kind must be 'solidus' or 'liquidus'") @@ -443,42 +475,44 @@ def wrapped(kind="solidus", Pmin=0.0, Pmax=1000.0, n=2000, **kwargs): # ============================================================================= SUPPORTED_MODELS = [ - "andrault_2011", - "monteux_2016", - "wolf_bower_2018", - "katz_2003", - "fei_2021", - "belonoshko_2005", - "fiquet_2010", - "hirschmann_2000", - "stixrude_2014", - "lin_2024", + 'andrault_2011', + 'monteux_2016', + 'wolf_bower_2018', + 'katz_2003', + 'fei_2021', + 'belonoshko_2005', + 'fiquet_2010', + 'hirschmann_2000', + 'stixrude_2014', + 'lin_2024', ] -def get_melting_curves(model_name: str, Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 2000, **kwargs): +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": andrault_2011_cut, - "monteux_2016": monteux_2016_cut, - "wolf_bower_2018": wolf_bower_2018_cut, - "katz_2003": katz_2003_cut, - "fei_2021": fei_2021, - "belonoshko_2005": belonoshko_2005, - "fiquet_2010": fiquet_2010, - "hirschmann_2000": hirschmann_2000, - "stixrude_2014": stixrude_2014, - "lin_2024": lin_2024, + 'andrault_2011': andrault_2011_cut, + 'monteux_2016': monteux_2016_cut, + 'wolf_bower_2018': wolf_bower_2018_cut, + 'katz_2003': katz_2003_cut, + '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}") + 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) + 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 @@ -487,12 +521,13 @@ def get_melting_curves(model_name: str, Pmin: float = 0.0, Pmax: float = 1000.0, # 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="#") + np.savetxt(path, data, fmt='%.18e %.18e', header='pressure temperature', comments='#') # ============================================================================= @@ -500,16 +535,16 @@ def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): # ============================================================================= SCRIPT_DIR = Path(__file__).resolve().parent -spider_dir = (SCRIPT_DIR / "../../../SPIDER").resolve() +spider_dir = (SCRIPT_DIR / '../../../SPIDER').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" +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}") + raise FileNotFoundError(f'Missing EOS file: {eos_solid_path}') if not eos_liquid_path.exists(): - raise FileNotFoundError(f"Missing EOS file: {eos_liquid_path}") + raise FileNotFoundError(f'Missing EOS file: {eos_liquid_path}') nP = 2020 nS_solid = 125 @@ -524,13 +559,15 @@ def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): SCALE_P_OUT = 1_000_000_000.0 SCALE_S_OUT = 4_824_266.84604467 -COMMON_HEADER = "\n".join([ - "# 5 2000", - "# Pressure, Entropy, Quantity", - "# column * scaling factor should be SI units", - "# scaling factors (constant) for each column given on line below", - "# 1000000000.0 4824266.84604467", -]) +COMMON_HEADER = '\n'.join( + [ + '# 5 2000', + '# Pressure, Entropy, Quantity', + '# column * scaling factor should be SI units', + '# scaling factors (constant) for each column given on line below', + '# 1000000000.0 4824266.84604467', + ] +) def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): @@ -553,7 +590,9 @@ def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): return S_axis, P_axis_GPa, T_interp -def invert_to_entropy_along_profile(P_gpa: np.ndarray, T_k: np.ndarray, S_axis: np.ndarray, T_of_SP): +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). """ @@ -584,11 +623,11 @@ def invert_to_entropy_along_profile(P_gpa: np.ndarray, T_k: np.ndarray, S_axis: continue try: - f = interp1d(T_unique, S_unique, kind="linear", assume_sorted=True) + 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) + f = interp1d(T_unique, S_unique, kind='nearest', assume_sorted=True) S_out[i] = float(f(T_i)) except Exception: continue @@ -620,7 +659,11 @@ def build_common_entropy_grid( 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: + 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: @@ -645,8 +688,22 @@ def build_common_entropy_grid( 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) + 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) @@ -661,7 +718,7 @@ def save_entropy_table_with_header(path: Path, P_gpa: np.ndarray, S_jpk: np.ndar """ P_pa = P_gpa * 1e9 data = np.column_stack([P_pa / SCALE_P_OUT, S_jpk / SCALE_S_OUT]) - np.savetxt(path, data, fmt="%.18e %.18e", header=COMMON_HEADER, comments="") + np.savetxt(path, data, fmt='%.18e %.18e', header=COMMON_HEADER, comments='') S_axis_solid, P_axis_solid, T_of_SP_solid = load_eos_T_of_SP( @@ -676,6 +733,7 @@ def save_entropy_table_with_header(path: Path, P_gpa: np.ndarray, S_jpk: np.ndar # MAIN EXPORTER # ============================================================================= + def export_model_curves( model_name: str, out_root: Path | str = MELTING_DIR, @@ -694,52 +752,52 @@ def export_model_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) + 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_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 - ) + 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 ) save_entropy_table_with_header( - out_dir / "solidus_P-S.dat", + out_dir / 'solidus_P-S.dat', P_common, S_sol_common, ) save_entropy_table_with_header( - out_dir / "liquidus_P-S.dat", + 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, + P_sol, + T_sol, + P_liq, + T_liq, + P_common, + S_sol_common, + S_liq_common, ) - print(f" Saved to : {out_dir.resolve()}") + print(f' Saved to : {out_dir.resolve()}') print() 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, + '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, } @@ -747,20 +805,21 @@ def export_model_curves( # BATCH EXPORTER # ============================================================================= + def export_all_models(out_root: Path | str = MELTING_DIR, n: int = 2000): r""" Export all supported melting parametrizations. """ for model in SUPPORTED_MODELS: - if model == "katz_2003": + if model == 'katz_2003': _ = export_model_curves(model, out_root=out_root, n=n, X_h2o=30.0) - elif model == "lin_2024": + elif model == 'lin_2024': _ = export_model_curves(model, out_root=out_root, n=n, fO2=-4.0) - elif model == "hirschmann_2000": + elif model == 'hirschmann_2000': _ = export_model_curves(model, out_root=out_root, n=n, Pmax=10.0) - elif model == "fei_2021": + elif model == 'fei_2021': _ = export_model_curves(model, out_root=out_root, n=n, Pmin=1.0) - elif model == "stixrude_2014": + elif model == 'stixrude_2014': _ = export_model_curves(model, out_root=out_root, n=n, Pmin=1.0) else: _ = export_model_curves(model, out_root=out_root, n=n) @@ -770,125 +829,126 @@ def export_all_models(out_root: Path | str = MELTING_DIR, n: int = 2000): # 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." + '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" + '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.", + '--all', + action='store_true', + help='Export all supported models.', ) parser.add_argument( - "--model", + '--model', type=str, default=None, choices=SUPPORTED_MODELS, - help="Export a single model by internal name.", + help='Export a single model by internal name.', ) parser.add_argument( - "--katz2003", - action="store_true", - help="Export Katz et al. (2003). Requires --X-h2o.", + '--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.", + '--lin2024', + action='store_true', + help='Export Lin et al. (2024). Requires --fO2.', ) parser.add_argument( - "--wolfbower2018", - action="store_true", - help="Export Wolf & Bower (2018).", + '--wolfbower2018', + action='store_true', + help='Export Wolf & Bower (2018).', ) parser.add_argument( - "--andrault2011", - action="store_true", - help="Export Andrault et al. (2011).", + '--andrault2011', + action='store_true', + help='Export Andrault et al. (2011).', ) parser.add_argument( - "--monteux2016", - action="store_true", - help="Export Monteux et al. (2016).", + '--monteux2016', + action='store_true', + help='Export Monteux et al. (2016).', ) parser.add_argument( - "--fei2021", - action="store_true", - help="Export Fei et al. (2021).", + '--fei2021', + action='store_true', + help='Export Fei et al. (2021).', ) parser.add_argument( - "--belonoshko2005", - action="store_true", - help="Export Belonoshko et al. (2005).", + '--belonoshko2005', + action='store_true', + help='Export Belonoshko et al. (2005).', ) parser.add_argument( - "--fiquet2010", - action="store_true", - help="Export Fiquet et al. (2010).", + '--fiquet2010', + action='store_true', + help='Export Fiquet et al. (2010).', ) parser.add_argument( - "--hirschmann2000", - action="store_true", - help="Export Hirschmann (2000).", + '--hirschmann2000', + action='store_true', + help='Export Hirschmann (2000).', ) parser.add_argument( - "--stixrude2014", - action="store_true", - help="Export Stixrude (2014).", + '--stixrude2014', + action='store_true', + help='Export Stixrude (2014).', ) parser.add_argument( - "--out-root", + '--out-root', type=str, default=str(MELTING_DIR), - help="Root directory where output folders will be created.", + help='Root directory where output folders will be created.', ) parser.add_argument( - "--Pmin", + '--Pmin', type=float, default=0.0, - help="Minimum pressure in GPa.", + help='Minimum pressure in GPa.', ) parser.add_argument( - "--Pmax", + '--Pmax', type=float, default=1000.0, - help="Maximum pressure in GPa.", + help='Maximum pressure in GPa.', ) parser.add_argument( - "-n", + '-n', type=int, default=2000, - help="Number of pressure samples.", + help='Number of pressure samples.', ) parser.add_argument( - "--X-h2o", - dest="X_h2o", + '--X-h2o', + dest='X_h2o', type=float, default=None, - help="Water content parameter for Katz (2003). Required for --katz2003.", + help='Water content parameter for Katz (2003). Required for --katz2003.', ) parser.add_argument( - "--fO2", + '--fO2', type=float, default=None, - help="Oxygen fugacity offset for Lin (2024). Required for --lin2024.", + help='Oxygen fugacity offset for Lin (2024). Required for --lin2024.', ) return parser.parse_args() @@ -899,22 +959,22 @@ 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", + '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.") + raise SystemExit('Error: please select only one model shortcut flag at a time.') if len(chosen) == 1: return chosen[0] @@ -932,31 +992,31 @@ def export_one_model_from_cli(model_name: str, args): Pmax = args.Pmax n = args.n - if model_name == "katz_2003": + 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" + '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 + kwargs['X_h2o'] = args.X_h2o - elif model_name == "lin_2024": + 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" + 'Error: --fO2 is required when using Lin (2024).\n' + 'Example: python solidus_func.py --lin2024 --fO2 -4' ) - kwargs["fO2"] = args.fO2 + kwargs['fO2'] = args.fO2 - elif model_name == "hirschmann_2000": + elif model_name == 'hirschmann_2000': if args.Pmax == 1000.0: Pmax = 10.0 - elif model_name == "fei_2021": + elif model_name == 'fei_2021': if args.Pmin == 0.0: Pmin = 1.0 - elif model_name == "stixrude_2014": + elif model_name == 'stixrude_2014': if args.Pmin == 0.0: Pmin = 1.0 @@ -979,7 +1039,7 @@ def main(): 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." + 'Error: please use either --all or a single model selection, not both.' ) export_all_models(out_root=args.out_root, n=args.n) return @@ -988,17 +1048,15 @@ def main(): if len(selected_models) == 0: raise SystemExit( - "Error: no model selected. Use --all or choose one model with " - "--model or a shortcut like --katz2003." + 'Error: no model selected. Use --all or choose one model with ' + '--model or a shortcut like --katz2003.' ) if len(selected_models) > 1: - raise SystemExit( - "Error: please choose only one of --model or one shortcut flag." - ) + raise SystemExit('Error: please choose only one of --model or one shortcut flag.') export_one_model_from_cli(selected_models[0], args) -if __name__ == "__main__": +if __name__ == '__main__': main() From 0cba4c5a52500d24bdad3d51f1e32e594abcaff2 Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Tue, 24 Mar 2026 16:28:36 +0100 Subject: [PATCH 09/13] double check testing --- docs/How-to/usage.md | 24 +- src/proteus/utils/solidus_func.py | 1311 ++++++++++++++++++++++++++++- 2 files changed, 1283 insertions(+), 52 deletions(-) diff --git a/docs/How-to/usage.md b/docs/How-to/usage.md index c5eb399ed..b577dc6db 100644 --- a/docs/How-to/usage.md +++ b/docs/How-to/usage.md @@ -303,16 +303,20 @@ PROTEUS uses precomputed solidus and liquidus curves from laboratory experiments ### Available parametrizations -- Andrault et al. (2011) -- Monteux et al. (2016) -- Wolf & Bower (2018) -- Katz et al. (2003) -- Lin et al. (2024) -- Hirschmann (2000) -- Stixrude (2014) -- Fei et al. (2021) -- Belonoshko et al. (2005) -- Fiquet et al. (2010) +Available melting_dir options +----------------------------- +The following directory names are supported and should be used exactly as written in the TOML configuration in the 'melting_dir' parameter: + +- andrault_2011 +- monteux_2016 +- wolf_bower_2018 +- katz_2003 +- fei_2021 +- belonoshko_2005 +- fiquet_2010 +- hirschmann_2000 +- stixrude_2014 +- lin_2024 ### Generate melting curves diff --git a/src/proteus/utils/solidus_func.py b/src/proteus/utils/solidus_func.py index bb5bfebd0..47b89f669 100644 --- a/src/proteus/utils/solidus_func.py +++ b/src/proteus/utils/solidus_func.py @@ -1,3 +1,1068 @@ +# #!/usr/bin/env python3 +# # -*- coding: utf-8 -*- +# """ +# 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: 10.1016/j.epsl.2011.02.006 +# - Monteux et al. (2016), DOI: 10.1016/j.epsl.2016.05.010 +# - Lin et al. (2024), DOI: 10.1038/s41561-024-01495-1 +# - Hirschmann (2000), DOI: 10.1029/2000GC000070 +# - Katz et al. (2003), DOI: 10.1029/2002GC000433 +# - Stixrude (2014), DOI: 10.1098/rsta.2013.0076 +# - Fei et al. (2021), DOI: 10.1038/s41467-021-21170-y +# - Belonoshko et al. (2005), DOI: 10.1103/PhysRevLett.94.195701 +# - Fiquet et al. (2010), DOI: 10.1126/science.1192448 +# """ + +# from __future__ import annotations + +# import argparse +# import os +# from pathlib import Path + +# import numpy as np +# import platformdirs +# from scipy.interpolate import RegularGridInterpolator, interp1d + +# # ============================================================================= +# # 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, +# ): +# """ +# Print a compact summary of the exported P-T and P-S ranges for one model. +# """ +# label = DISPLAY_NAMES.get(model_name, model_name) +# print(f'{label}:') +# print( +# f' P-T solidus : P_range = {_fmt_range(P_sol)} GPa, T_range = {_fmt_range(T_sol)} K' +# ) +# print( +# f' P-T liquidus : P_range = {_fmt_range(P_liq)} GPa, T_range = {_fmt_range(T_liq)} K' +# ) +# print( +# f' P-S solidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_sol_common)} J kg^-1 K^-1' +# ) +# print( +# f' P-S liquidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_liq_common)} J kg^-1 K^-1' +# ) +# print() + + +# # ============================================================================= +# # 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.004054 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 + + +# andrault_2011_cut = truncate_to_physical_interval(andrault_2011) +# monteux_2016_cut = truncate_to_physical_interval(monteux_2016) +# wolf_bower_2018_cut = truncate_to_physical_interval(wolf_bower_2018) +# katz_2003_cut = truncate_to_physical_interval(katz_2003) + + +# # ============================================================================= +# # MODEL DISPATCHER +# # ============================================================================= + +# SUPPORTED_MODELS = [ +# 'andrault_2011', +# 'monteux_2016', +# 'wolf_bower_2018', +# 'katz_2003', +# 'fei_2021', +# 'belonoshko_2005', +# 'fiquet_2010', +# 'hirschmann_2000', +# 'stixrude_2014', +# 'lin_2024', +# ] + + +# 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': andrault_2011_cut, +# 'monteux_2016': monteux_2016_cut, +# 'wolf_bower_2018': wolf_bower_2018_cut, +# 'katz_2003': katz_2003_cut, +# '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='#') + + +# # ============================================================================= +# # EOS LOOKUP TABLE SETTINGS +# # ============================================================================= + +# SCRIPT_DIR = Path(__file__).resolve().parent +# spider_dir = (SCRIPT_DIR / '../../../SPIDER').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}') + +# nP = 2020 +# nS_solid = 125 +# nS_liquid = 95 +# skip_header = 5 + +# SCALE_P_EOS = 1e9 +# SCALE_T_EOS = 1.0 +# SCALE_S_SOLID_EOS = 4.82426684604467e6 +# SCALE_S_LIQUID_EOS = 4.805046659407042e6 + +# SCALE_P_OUT = 1_000_000_000.0 +# SCALE_S_OUT = 4_824_266.84604467 + +# COMMON_HEADER = '\n'.join( +# [ +# '# 5 2000', +# '# Pressure, Entropy, Quantity', +# '# column * scaling factor should be SI units', +# '# scaling factors (constant) for each column given on line below', +# '# 1000000000.0 4824266.84604467', +# ] +# ) + + +# def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): +# r""" +# Load an EOS lookup table and build an interpolator for T(S, P). +# """ +# 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 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): +# r""" +# Save a pressure-entropy table in SPIDER-style scaled format. +# """ +# P_pa = P_gpa * 1e9 +# data = np.column_stack([P_pa / SCALE_P_OUT, S_jpk / SCALE_S_OUT]) +# np.savetxt(path, data, fmt='%.18e %.18e', header=COMMON_HEADER, comments='') + + +# S_axis_solid, P_axis_solid, T_of_SP_solid = load_eos_T_of_SP( +# eos_solid_path, nS_solid, SCALE_S_SOLID_EOS +# ) +# S_axis_liquid, P_axis_liquid, T_of_SP_liquid = load_eos_T_of_SP( +# eos_liquid_path, nS_liquid, SCALE_S_LIQUID_EOS +# ) + + +# # ============================================================================= +# # 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, +# **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_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 +# ) + +# 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, +# ) + +# print(f' Saved to : {out_dir.resolve()}') +# print() + +# 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): +# 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) +# elif model == 'lin_2024': +# _ = export_model_curves(model, out_root=out_root, n=n, fO2=-4.0) +# elif model == 'hirschmann_2000': +# _ = export_model_curves(model, out_root=out_root, n=n, Pmax=10.0) +# elif model == 'fei_2021': +# _ = export_model_curves(model, out_root=out_root, n=n, Pmin=1.0) +# elif model == 'stixrude_2014': +# _ = export_model_curves(model, out_root=out_root, n=n, Pmin=1.0) +# else: +# _ = export_model_curves(model, out_root=out_root, n=n) + + +# # ============================================================================= +# # 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( +# '--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, +# **kwargs, +# ) + + +# def main(): +# 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.' +# ) +# export_all_models(out_root=args.out_root, n=args.n) +# return + +# selected_models = [m for m in [explicit_model, shortcut_model] if m is not None] + +# if len(selected_models) == 0: +# raise SystemExit( +# 'Error: no model selected. Use --all or choose one model with ' +# '--model or a shortcut like --katz2003.' +# ) + +# if len(selected_models) > 1: +# raise SystemExit('Error: please choose only one of --model or one shortcut flag.') + +# export_one_model_from_cli(selected_models[0], args) + + +# if __name__ == '__main__': +# main() + + + #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ @@ -69,7 +1134,6 @@ # ============================================================================= 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' @@ -242,7 +1306,7 @@ def fiquet_2010(kind: str = 'liquidus', Pmin: float = 0.0, Pmax: float = 1000.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.004054 GPa. + # 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) @@ -534,18 +1598,6 @@ def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): # EOS LOOKUP TABLE SETTINGS # ============================================================================= -SCRIPT_DIR = Path(__file__).resolve().parent -spider_dir = (SCRIPT_DIR / '../../../SPIDER').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}') - nP = 2020 nS_solid = 125 nS_liquid = 95 @@ -559,15 +1611,125 @@ def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): SCALE_P_OUT = 1_000_000_000.0 SCALE_S_OUT = 4_824_266.84604467 -COMMON_HEADER = '\n'.join( - [ - '# 5 2000', - '# Pressure, Entropy, Quantity', - '# column * scaling factor should be SI units', - '# scaling factors (constant) for each column given on line below', - '# 1000000000.0 4824266.84604467', - ] -) +def make_entropy_header(n_rows: int) -> str: + """ + Build the SPIDER-style header for a pressure-entropy table. + + Parameters + ---------- + n_rows : int + Number of data rows in the file. + + 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', + '# 1000000000.0 4824266.84604467', + ] + ) + +def get_default_spider_dir() -> Path: + """ + Return the default SPIDER directory relative to this module. + """ + script_dir = Path(__file__).resolve().parent + return (script_dir / '../../../SPIDER').resolve() + +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): @@ -590,6 +1752,39 @@ def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): return S_axis, P_axis_GPa, T_interp +def load_default_eos_interpolators( + spider_dir: Path | str | None = None, +) -> 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. + + 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 + ) + S_axis_liquid, _, T_of_SP_liquid = load_eos_T_of_SP( + eos_liquid_path, nS_liquid, SCALE_S_LIQUID_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 ): @@ -718,16 +1913,8 @@ def save_entropy_table_with_header(path: Path, P_gpa: np.ndarray, S_jpk: np.ndar """ P_pa = P_gpa * 1e9 data = np.column_stack([P_pa / SCALE_P_OUT, S_jpk / SCALE_S_OUT]) - np.savetxt(path, data, fmt='%.18e %.18e', header=COMMON_HEADER, comments='') - - -S_axis_solid, P_axis_solid, T_of_SP_solid = load_eos_T_of_SP( - eos_solid_path, nS_solid, SCALE_S_SOLID_EOS -) -S_axis_liquid, P_axis_liquid, T_of_SP_liquid = load_eos_T_of_SP( - eos_liquid_path, nS_liquid, SCALE_S_LIQUID_EOS -) - + header = make_entropy_header(len(P_gpa)) + np.savetxt(path, data, fmt='%.18e %.18e', header=header, comments='') # ============================================================================= # MAIN EXPORTER @@ -740,6 +1927,7 @@ def export_model_curves( Pmin: float = 0.0, Pmax: float = 1000.0, n: int = 2000, + spider_dir: Path | str | None = None, **kwargs, ): r""" @@ -755,6 +1943,10 @@ def export_model_curves( 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) @@ -762,6 +1954,13 @@ def export_model_curves( 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, @@ -806,23 +2005,37 @@ def export_model_curves( # ============================================================================= -def export_all_models(out_root: Path | str = MELTING_DIR, n: int = 2000): +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) + _ = 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) + _ = 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) + _ = 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) + _ = 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) + _ = 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) + _ = export_model_curves(model, out_root=out_root, n=n, spider_dir=spider_dir) # ============================================================================= @@ -918,6 +2131,13 @@ def parse_args(): 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, @@ -1026,6 +2246,7 @@ def export_one_model_from_cli(model_name: str, args): Pmin=Pmin, Pmax=Pmax, n=n, + spider_dir=args.spider_dir, **kwargs, ) @@ -1041,7 +2262,10 @@ def main(): raise SystemExit( 'Error: please use either --all or a single model selection, not both.' ) - export_all_models(out_root=args.out_root, n=args.n) + 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] @@ -1055,7 +2279,10 @@ def main(): if len(selected_models) > 1: raise SystemExit('Error: please choose only one of --model or one shortcut flag.') - export_one_model_from_cli(selected_models[0], args) + try: + export_one_model_from_cli(selected_models[0], args) + except FileNotFoundError as exc: + raise SystemExit(f'Error: {exc}') from exc if __name__ == '__main__': From f394895e79e4372ed64c6860d1ebe3345b81b073 Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Tue, 24 Mar 2026 16:38:46 +0100 Subject: [PATCH 10/13] modify solifus script --- src/proteus/utils/solidus_func.py | 1067 ----------------------------- 1 file changed, 1067 deletions(-) diff --git a/src/proteus/utils/solidus_func.py b/src/proteus/utils/solidus_func.py index 47b89f669..2fe51cb65 100644 --- a/src/proteus/utils/solidus_func.py +++ b/src/proteus/utils/solidus_func.py @@ -1,1070 +1,3 @@ -# #!/usr/bin/env python3 -# # -*- coding: utf-8 -*- -# """ -# 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: 10.1016/j.epsl.2011.02.006 -# - Monteux et al. (2016), DOI: 10.1016/j.epsl.2016.05.010 -# - Lin et al. (2024), DOI: 10.1038/s41561-024-01495-1 -# - Hirschmann (2000), DOI: 10.1029/2000GC000070 -# - Katz et al. (2003), DOI: 10.1029/2002GC000433 -# - Stixrude (2014), DOI: 10.1098/rsta.2013.0076 -# - Fei et al. (2021), DOI: 10.1038/s41467-021-21170-y -# - Belonoshko et al. (2005), DOI: 10.1103/PhysRevLett.94.195701 -# - Fiquet et al. (2010), DOI: 10.1126/science.1192448 -# """ - -# from __future__ import annotations - -# import argparse -# import os -# from pathlib import Path - -# import numpy as np -# import platformdirs -# from scipy.interpolate import RegularGridInterpolator, interp1d - -# # ============================================================================= -# # 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, -# ): -# """ -# Print a compact summary of the exported P-T and P-S ranges for one model. -# """ -# label = DISPLAY_NAMES.get(model_name, model_name) -# print(f'{label}:') -# print( -# f' P-T solidus : P_range = {_fmt_range(P_sol)} GPa, T_range = {_fmt_range(T_sol)} K' -# ) -# print( -# f' P-T liquidus : P_range = {_fmt_range(P_liq)} GPa, T_range = {_fmt_range(T_liq)} K' -# ) -# print( -# f' P-S solidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_sol_common)} J kg^-1 K^-1' -# ) -# print( -# f' P-S liquidus : P_range = {_fmt_range(P_common)} GPa, S_range = {_fmt_range(S_liq_common)} J kg^-1 K^-1' -# ) -# print() - - -# # ============================================================================= -# # 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.004054 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 - - -# andrault_2011_cut = truncate_to_physical_interval(andrault_2011) -# monteux_2016_cut = truncate_to_physical_interval(monteux_2016) -# wolf_bower_2018_cut = truncate_to_physical_interval(wolf_bower_2018) -# katz_2003_cut = truncate_to_physical_interval(katz_2003) - - -# # ============================================================================= -# # MODEL DISPATCHER -# # ============================================================================= - -# SUPPORTED_MODELS = [ -# 'andrault_2011', -# 'monteux_2016', -# 'wolf_bower_2018', -# 'katz_2003', -# 'fei_2021', -# 'belonoshko_2005', -# 'fiquet_2010', -# 'hirschmann_2000', -# 'stixrude_2014', -# 'lin_2024', -# ] - - -# 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': andrault_2011_cut, -# 'monteux_2016': monteux_2016_cut, -# 'wolf_bower_2018': wolf_bower_2018_cut, -# 'katz_2003': katz_2003_cut, -# '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='#') - - -# # ============================================================================= -# # EOS LOOKUP TABLE SETTINGS -# # ============================================================================= - -# SCRIPT_DIR = Path(__file__).resolve().parent -# spider_dir = (SCRIPT_DIR / '../../../SPIDER').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}') - -# nP = 2020 -# nS_solid = 125 -# nS_liquid = 95 -# skip_header = 5 - -# SCALE_P_EOS = 1e9 -# SCALE_T_EOS = 1.0 -# SCALE_S_SOLID_EOS = 4.82426684604467e6 -# SCALE_S_LIQUID_EOS = 4.805046659407042e6 - -# SCALE_P_OUT = 1_000_000_000.0 -# SCALE_S_OUT = 4_824_266.84604467 - -# COMMON_HEADER = '\n'.join( -# [ -# '# 5 2000', -# '# Pressure, Entropy, Quantity', -# '# column * scaling factor should be SI units', -# '# scaling factors (constant) for each column given on line below', -# '# 1000000000.0 4824266.84604467', -# ] -# ) - - -# def load_eos_T_of_SP(eos_path: Path, nS: int, scale_S_axis: float): -# r""" -# Load an EOS lookup table and build an interpolator for T(S, P). -# """ -# 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 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): -# r""" -# Save a pressure-entropy table in SPIDER-style scaled format. -# """ -# P_pa = P_gpa * 1e9 -# data = np.column_stack([P_pa / SCALE_P_OUT, S_jpk / SCALE_S_OUT]) -# np.savetxt(path, data, fmt='%.18e %.18e', header=COMMON_HEADER, comments='') - - -# S_axis_solid, P_axis_solid, T_of_SP_solid = load_eos_T_of_SP( -# eos_solid_path, nS_solid, SCALE_S_SOLID_EOS -# ) -# S_axis_liquid, P_axis_liquid, T_of_SP_liquid = load_eos_T_of_SP( -# eos_liquid_path, nS_liquid, SCALE_S_LIQUID_EOS -# ) - - -# # ============================================================================= -# # 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, -# **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_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 -# ) - -# 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, -# ) - -# print(f' Saved to : {out_dir.resolve()}') -# print() - -# 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): -# 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) -# elif model == 'lin_2024': -# _ = export_model_curves(model, out_root=out_root, n=n, fO2=-4.0) -# elif model == 'hirschmann_2000': -# _ = export_model_curves(model, out_root=out_root, n=n, Pmax=10.0) -# elif model == 'fei_2021': -# _ = export_model_curves(model, out_root=out_root, n=n, Pmin=1.0) -# elif model == 'stixrude_2014': -# _ = export_model_curves(model, out_root=out_root, n=n, Pmin=1.0) -# else: -# _ = export_model_curves(model, out_root=out_root, n=n) - - -# # ============================================================================= -# # 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( -# '--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, -# **kwargs, -# ) - - -# def main(): -# 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.' -# ) -# export_all_models(out_root=args.out_root, n=args.n) -# return - -# selected_models = [m for m in [explicit_model, shortcut_model] if m is not None] - -# if len(selected_models) == 0: -# raise SystemExit( -# 'Error: no model selected. Use --all or choose one model with ' -# '--model or a shortcut like --katz2003.' -# ) - -# if len(selected_models) > 1: -# raise SystemExit('Error: please choose only one of --model or one shortcut flag.') - -# export_one_model_from_cli(selected_models[0], args) - - -# if __name__ == '__main__': -# main() - - - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Generate melting curves in pressure–temperature (P–T) and pressure–entropy (P–S) space for several literature parametrizations of peridotite / silicate melting. From ca37dcef75db5f3c36fd3485908bdf76d06c178e Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Tue, 24 Mar 2026 16:41:18 +0100 Subject: [PATCH 11/13] ruff check --- src/proteus/utils/solidus_func.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/proteus/utils/solidus_func.py b/src/proteus/utils/solidus_func.py index 2fe51cb65..ee14570fe 100644 --- a/src/proteus/utils/solidus_func.py +++ b/src/proteus/utils/solidus_func.py @@ -544,6 +544,7 @@ def save_PT_table(path: Path, P_gpa: np.ndarray, T_k: np.ndarray): SCALE_P_OUT = 1_000_000_000.0 SCALE_S_OUT = 4_824_266.84604467 + def make_entropy_header(n_rows: int) -> str: """ Build the SPIDER-style header for a pressure-entropy table. @@ -568,6 +569,7 @@ def make_entropy_header(n_rows: int) -> str: ] ) + def get_default_spider_dir() -> Path: """ Return the default SPIDER directory relative to this module. @@ -575,6 +577,7 @@ def get_default_spider_dir() -> Path: script_dir = Path(__file__).resolve().parent return (script_dir / '../../../SPIDER').resolve() + def validate_entropy_export_arrays( P_common: np.ndarray, S_sol_common: np.ndarray, @@ -620,9 +623,8 @@ def validate_entropy_export_arrays( 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.' - ) + 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]: """ @@ -849,6 +851,7 @@ def save_entropy_table_with_header(path: Path, P_gpa: np.ndarray, S_jpk: np.ndar header = make_entropy_header(len(P_gpa)) np.savetxt(path, data, fmt='%.18e %.18e', header=header, comments='') + # ============================================================================= # MAIN EXPORTER # ============================================================================= @@ -876,8 +879,8 @@ def export_model_curves( 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_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) From 662f86e217fcdc88583ffb4d253ea0ceae6f63ca Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Tue, 24 Mar 2026 17:23:42 +0100 Subject: [PATCH 12/13] add test for data modifications --- tests/utils/test_data.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/tests/utils/test_data.py b/tests/utils/test_data.py index aaadc4004..e2718dbed 100644 --- a/tests/utils/test_data.py +++ b/tests/utils/test_data.py @@ -2338,3 +2338,31 @@ 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() From a0e7d02470df8e28ff75abf2ddf9538b554ee6d3 Mon Sep 17 00:00:00 2001 From: Mariana Sastre Date: Tue, 24 Mar 2026 17:26:01 +0100 Subject: [PATCH 13/13] Format test with ruff --- tests/utils/test_data.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/utils/test_data.py b/tests/utils/test_data.py index e2718dbed..9de33918e 100644 --- a/tests/utils/test_data.py +++ b/tests/utils/test_data.py @@ -2340,7 +2340,6 @@ def test_download_eos_static_delegates(mock_seager): mock_seager.assert_called_once() - @pytest.mark.unit @patch('proteus.utils.data.download') @patch('proteus.utils.data.GetFWLData')