diff --git a/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb b/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb new file mode 100644 index 0000000..6f31633 --- /dev/null +++ b/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb @@ -0,0 +1,593 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Garnet–Clinopyroxene Fe²⁺–Mg Exchange Thermometry (S10 Eclogite)\n", + "\n", + "Garnet and clinopyroxene coexist over a wide range of high-pressure mafic rocks, from\n", + "granulites to eclogites, and the Fe²⁺–Mg distribution between them is a long-established\n", + "thermometer. The exchange reaction\n", + "\n", + "$$\n", + "\\tfrac{1}{3}\\mathrm{Mg_3Al_2Si_3O_{12}}\\ (\\text{pyrope}) + \\mathrm{CaFeSi_2O_6}\\ (\\text{hedenbergite})\n", + "\\;=\\;\n", + "\\tfrac{1}{3}\\mathrm{Fe_3Al_2Si_3O_{12}}\\ (\\text{almandine}) + \\mathrm{CaMgSi_2O_6}\\ (\\text{diopside})\n", + "$$\n", + "\n", + "has an equilibrium distribution coefficient\n", + "\n", + "$$\n", + "K_D = \\frac{(\\mathrm{Fe^{2+}/Mg})^{\\mathrm{garnet}}}{(\\mathrm{Fe^{2+}/Mg})^{\\mathrm{clinopyroxene}}}\n", + "$$\n", + "\n", + "whose temperature dependence is the basis of the thermometer: $K_D$ falls as temperature rises.\n", + "\n", + "This notebook:\n", + "\n", + "1. sets up the S10 eclogite bulk composition in the `mpe` database;\n", + "2. runs a P–T grid and extracts the garnet and clinopyroxene compositions;\n", + "3. maps $K_D$ across the grid using the traditional mixed Fe convention: garnet Fe is taken from its FeO-basis APFU (all Fe treated as Fe²⁺), while cpx Fe²⁺ comes from the excess-O heuristic split;\n", + "4. converts $K_D$ to temperature with Mysen & Heier (1972), Ganguly (1979), Krogh Ravna (2000), Ellis & Green (1979), and Räheim & Green (1974), and compares them against the true grid temperature.\n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## 1. Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from phasetools import MAGEMinPTGridCalculator" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## 2. Bulk composition and database\n", + "\n", + "The S10 eclogite bulk composition (in mol%) is taken from the `phasetools` test suite. The\n", + "`mpe` database combines the metapelite (`mp`) and metabasite (`mb`) solution models and is\n", + "the appropriate choice for a mafic bulk at eclogite-facies conditions. Pressures are in kbar\n", + "and temperatures in °C throughout." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "# Eclogite S10 bulk composition in mol% (from the phasetools test suite).\n", + "Xoxides = [\"H2O\", \"SiO2\", \"Al2O3\", \"CaO\", \"MgO\", \"FeO\", \"K2O\", \"Na2O\", \"TiO2\", \"MnO\", \"O\"]\n", + "X = [0.92, 54.57, 8.79, 11.20, 8.45, 12.89, 0.24, 2.24, 1.12, 0.22, 0.64]\n", + "\n", + "db = \"mpe\" # cross-lithology metapelite + metabasite solution models\n", + "dataset = 636\n", + "sys_in = \"mol\"\n", + "\n", + "calc = MAGEMinPTGridCalculator(db=db, dataset=dataset)\n", + "calc.setup_bulk_composition(Xoxides, X, sys_in=sys_in)" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## 3. P–T grid\n", + "\n", + "A grid spanning 15–40 kbar and 550–850 °C covers the eclogite-facies field where garnet and\n", + "clinopyroxene are both stable for this bulk composition." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "P = np.linspace(5.0, 20.0, 11) # kbar\n", + "T = np.linspace(350.0, 900.0, 13) # °C\n", + "Pgrid, Tgrid = np.meshgrid(P, T, indexing=\"xy\")\n", + "\n", + "out = calc.calculate_grid(Pgrid.ravel(), Tgrid.ravel())\n", + "print(f\"{len(out)} grid points; unique phases: {calc.get_all_unique_phases(out)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## 4. Resolve the clinopyroxene phase name\n", + "\n", + "MAGEMin reports clinopyroxene under different solution-model names depending on composition\n", + "(`dio` diopside, `omph` omphacite, `aug` augite). The name is resolved dynamically rather than\n", + "hard-coded, so the notebook follows whichever clinopyroxene solution is stable." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "CPX_CANDIDATES = [\"dio\", \"omph\", \"aug\", \"jac\"] # clinopyroxene solution names in mpe\n", + "\n", + "unique_phases = calc.get_all_unique_phases(out)\n", + "cpx = next((p for p in CPX_CANDIDATES if p in unique_phases), None)\n", + "if cpx is None:\n", + " raise RuntimeError(\"No clinopyroxene phase found in the grid.\")\n", + "print(\"Clinopyroxene phase resolved to:\", cpx)" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## 5. Extract Fe²⁺/Mg and compute $K_D$\n", + "\n", + "Only **divalent** iron takes part in the Fe²⁺–Mg exchange. For the traditional calibration convention\n", + "used here, Fe²⁺ would ideally be taken from the FeO basis for garnet, but for simplicity this\n", + "function uses the excess-O heuristic Fe²⁺/Fe³⁺ split consistently for all phases. Grid points\n", + "where a phase is absent — or where a solvus yields two coexisting limbs of the same phase, making\n", + "the pairing with the other phase ambiguous — are masked out." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "def extract_fe2_mg(calc, phase, grid_out):\n", + " \"\"\"Return (Fe2+/Mg, X_Ca, X_Mn, mask) arrays over the grid for ``phase``.\n", + "\n", + " For the traditional calibration convention used here, Fe2+ would ideally be\n", + " taken from the FeO basis for garnet, but for simplicity this function uses the\n", + " excess-O heuristic Fe2+/Fe3+ split consistently for all phases. Grid points where\n", + " the phase is absent, or where a solvus\n", + " gives two coexisting limbs of the same phase (an ambiguous pairing with the\n", + " other phase), are flagged True in ``mask`` and should be excluded before\n", + " forming a distribution coefficient.\n", + "\n", + " Parameters\n", + " ----------\n", + " calc : MAGEMinPTGridCalculator\n", + " Calculator that produced ``grid_out``.\n", + " phase : str\n", + " Phase key (e.g. ``\"g\"``).\n", + " grid_out : list\n", + " Output of :meth:`MAGEMinPTGridCalculator.calculate_grid`.\n", + "\n", + " Returns\n", + " -------\n", + " fe2_mg, x_ca, x_mn, mask : numpy.ndarray\n", + " Arrays over the grid. ``x_ca`` and ``x_mn`` are the garnet-style cation\n", + " fractions Ca/(Ca+Mn+Fe+Mg) and Mn/(Ca+Mn+Fe+Mg), with the Fe\n", + " convention above.\n", + " \"\"\"\n", + " bundles = calc.extract_from_grid(\n", + " phase, oxides=[\"MgO\", \"FeO\", \"CaO\", \"MnO\"], fe_split=True, grid_out=grid_out\n", + " )\n", + " n = len(grid_out)\n", + " if not bundles: # phase never appears in the grid\n", + " nan = np.full(n, np.nan)\n", + " return nan, nan, nan, np.ones(n, dtype=bool)\n", + "\n", + " b0 = bundles[0]\n", + " mg = b0[\"ox_apfu_MgO\"]\n", + " # Alternative: use the FeO-basis APFU for garnet to match the strict\n", + " # traditional convention (all garnet Fe treated as Fe2+).\n", + " # Kept here for reference; the active code uses the heuristic split.\n", + " fe2 = b0[\"Fe2\"]\n", + " ca = b0[\"ox_apfu_CaO\"]\n", + " mn = b0[\"ox_apfu_MnO\"]\n", + "\n", + " denom = ca + mn + fe2 + mg # Ca + Mn + Fe2+ + Mg, in atoms per formula unit\n", + " with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n", + " fe2_mg = fe2 / mg\n", + " x_ca = ca / denom\n", + " x_mn = mn / denom\n", + "\n", + " absent = np.isnan(fe2) | np.isnan(mg) | (mg <= 0.0)\n", + " if len(bundles) > 1:\n", + " ambiguous = ~np.isnan(bundles[1][\"Fe2\"]) # a second limb exists -> ambiguous pairing\n", + " else:\n", + " ambiguous = np.zeros(n, dtype=bool)\n", + " return fe2_mg, x_ca, x_mn, absent | ambiguous" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "fe2_mg_g, x_ca_g, x_mn_g, mask_g = extract_fe2_mg(calc, \"g\", out)\n", + "fe2_mg_c, x_ca_c, x_mn_c, mask_c = extract_fe2_mg(calc, cpx, out)\n", + "\n", + "valid = ~mask_g & ~mask_c\n", + "\n", + "Kd = np.full_like(fe2_mg_g, np.nan)\n", + "Kd[valid] = fe2_mg_g[valid] / fe2_mg_c[valid]\n", + "\n", + "print(f\"K_D defined at {int(valid.sum())} / {valid.size} grid points\")\n", + "print(f\"K_D range: {np.nanmin(Kd):.2f} – {np.nanmax(Kd):.2f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## 6. Map $K_D$ over the P–T grid\n", + "\n", + "$K_D$ decreases smoothly with rising temperature — the signature of the Fe²⁺–Mg exchange\n", + "equilibrium — and is only weakly dependent on pressure, which is why a pressure-independent\n", + "calibration such as Mysen & Heier (1972) can work at all." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6.5, 5))\n", + "Kd_map = Kd.reshape(Tgrid.shape)\n", + "\n", + "cs = ax.contourf(Tgrid, Pgrid, Kd_map, levels=14, cmap=\"viridis\")\n", + "ax.contour(Tgrid, Pgrid, Kd_map, levels=14, colors=\"k\", linewidths=0.4)\n", + "cbar = fig.colorbar(cs, ax=ax, label=r\"$K_D = (Fe^{2+}/Mg)^g \\,/\\, (Fe^{2+}/Mg)^{cpx}$\")\n", + "\n", + "ax.set_xlabel(\"Temperature (°C)\")\n", + "ax.set_ylabel(\"Pressure (kbar)\")\n", + "ax.set_title(\"Fe²⁺–Mg distribution coefficient, S10 eclogite (mpe)\")\n", + "ax.grid(ls=\":\", alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "## 7. Thermometer calibrations\n", + "\n", + "Five calibrations of the Fe²⁺–Mg exchange are implemented as pure-NumPy helpers. All return °C.\n", + "\n", + "**Mysen & Heier (1972)** — a pressure- and composition-independent calibration\n", + "(following Banno, 1970):\n", + "\n", + "$$\n", + "T\\,(\\mathrm{K}) = \\frac{2475}{\\ln K_D + 0.781}\n", + "$$\n", + "\n", + "**Ganguly (1979)** — includes pressure and garnet grossular content. This notebook retains the\n", + "piecewise coefficients tabulated by Yavuz & Yıldırım (2020), rather than silently attributing\n", + "that printed piecewise form to the original Ganguly equation.\n", + "in the piecewise form tabulated by Yavuz & Yıldırım (2020):\n", + "\n", + "$$\n", + "T\\,(°\\mathrm{C}) = \\frac{4100 + 1586\\,X_\\mathrm{Ca}^{g} + 11.07\\,P}{\\ln K_D + 2.40} - 273.15\n", + "\\quad (T \\geq 1060\\,°\\mathrm{C})\n", + "$$\n", + "\n", + "$$\n", + "T\\,(°\\mathrm{C}) = \\frac{4801 + 1586\\,X_\\mathrm{Ca}^{g} + 11.07\\,P}{\\ln K_D + 2.93} - 273.15\n", + "\\quad (T \\leq 1060\\,°\\mathrm{C})\n", + "$$\n", + "\n", + "where $X_\\mathrm{Ca}^{g} = \\mathrm{Ca}/(\\mathrm{Ca} + \\mathrm{Mn} + \\mathrm{Fe^{2+}} + \\mathrm{Mg})$\n", + "in garnet and $P$ is in kbar.\n", + "\n", + "**Krogh Ravna (2000)** uses $P_{GPa}=P_{kbar}/10$, $X_{Ca}=Ca/(Ca+Mn+Fe+Mg)$,\n", + "$X_{Mn}=Mn/(Ca+Mn+Fe+Mg)$, and $X_{Mg\\#}=Mg/(Mg+Fe)$.\n", + "\n", + "**Ellis & Green (1979)** uses $X_{Ca,binary}=Ca/(Ca+Fe+Mg)$, deliberately excluding Mn;\n", + "this differs from the garnet site $X_{Ca}$ used by Krogh Ravna and Ganguly.\n", + "\n", + "**Räheim & Green (1974)** is pressure-dependent but composition-independent." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "def T_mysen_heier_1972(Kd):\n", + " \"\"\"Mysen & Heier (1972) garnet–cpx Fe2+–Mg thermometer, in °C.\n", + "\n", + " T(K) = 2475 / (ln Kd + 0.781)\n", + " \"\"\"\n", + " lnKd = np.log(Kd)\n", + " return 2475.0 / (lnKd + 0.781) - 273.15\n", + "\n", + "\n", + "def T_ganguly_1979(Kd, P_kbar, X_Ca_grt):\n", + " \"\"\"Ganguly (1979) garnet–cpx Fe2+–Mg thermometer, in °C.\n", + "\n", + " Piecewise form tabulated by Yavuz & Yıldırım (2020, WinGrt):\n", + " T >= 1060 °C: T = (4100 + 1586·X_Ca + 11.07·P) / (ln Kd + 2.40) - 273.15\n", + " T <= 1060 °C: T = (4801 + 1586·X_Ca + 11.07·P) / (ln Kd + 2.93) - 273.15\n", + " \"\"\"\n", + " lnKd = np.log(Kd)\n", + " P = np.asarray(P_kbar, dtype=float)\n", + " XCa = np.asarray(X_Ca_grt, dtype=float)\n", + " T_hi = (4100.0 + 1586.0 * XCa + 11.07 * P) / (lnKd + 2.40) - 273.15\n", + " T_lo = (4801.0 + 1586.0 * XCa + 11.07 * P) / (lnKd + 2.93) - 273.15\n", + " return np.where(T_hi >= 1060.0, T_hi, T_lo)\n", + "\n", + "\n", + "def T_krogh_ravna_2000(Kd, P_kbar, X_Ca, X_Mn, X_Mg_number):\n", + " \"\"\"Krogh Ravna (2000) garnet-cpx thermometer, in °C.\"\"\"\n", + " lnKd = np.log(Kd)\n", + " P_GPa = np.asarray(P_kbar, dtype=float) / 10.0\n", + " numerator = (1939.9 + 3270.0 * X_Ca - 1396.0 * X_Ca**2\n", + " + 3319.0 * X_Mn + 3535.0 * X_Mn**2\n", + " + 1105.0 * X_Mg_number - 3561.0 * X_Mg_number**2\n", + " + 2324.0 * X_Mg_number**3 + 169.4 * P_GPa)\n", + " return numerator / (lnKd + 1.223) - 273.15\n", + "\n", + "\n", + "def T_ellis_green_1979(Kd, P_kbar, X_Ca_binary):\n", + " \"\"\"Ellis & Green (1979) garnet-cpx thermometer, in °C.\"\"\"\n", + " return (3104.0 * X_Ca_binary + 3030.0 + 10.86 * P_kbar) / (np.log(Kd) + 1.9034) - 273.15\n", + "\n", + "\n", + "def T_raheim_green_1974(Kd, P_kbar):\n", + " \"\"\"Räheim & Green (1974) garnet-cpx thermometer, in °C.\"\"\"\n", + " return (3686.0 + 28.35 * P_kbar) / (np.log(Kd) + 2.33) - 273.15" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "## 8. Apply the thermometers and compare with the true temperature\n", + "\n", + "Because every grid point has a known temperature, we can feed the modelled $K_D$ back into each\n", + "thermometer and see how well it recovers the input. Points lie above the 1:1 line when a\n", + "calibration overestimates the temperature." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "P_flat = Pgrid.ravel()\n", + "T_true = Tgrid.ravel()\n", + "\n", + "T_MH = T_mysen_heier_1972(Kd)\n", + "T_G = T_ganguly_1979(Kd, P_flat, x_ca_g)\n", + "X_mg_number_g = fe2_mg_g / (1.0 + fe2_mg_g)\n", + "X_ca_binary_g = x_ca_g / (1.0 - x_mn_g)\n", + "T_KR = T_krogh_ravna_2000(Kd, P_flat, x_ca_g, x_mn_g, X_mg_number_g)\n", + "T_EG = T_ellis_green_1979(Kd, P_flat, X_ca_binary_g)\n", + "T_RG = T_raheim_green_1974(Kd, P_flat)\n", + "thermometers = {\n", + " \"Mysen & Heier (1972)\": T_MH,\n", + " \"Ganguly (1979; Yavuz & Yıldırım piecewise)\": T_G,\n", + " \"Krogh Ravna (2000)\": T_KR,\n", + " \"Ellis & Green (1979)\": T_EG,\n", + " \"Räheim & Green (1974)\": T_RG,\n", + "}\n", + "\n", + "fig, ax = plt.subplots(figsize=(7, 5.5))\n", + "ax.plot([Tgrid.min(), Tgrid.max()], [Tgrid.min(), Tgrid.max()], \"k--\", lw=1, label=\"1:1 (true T)\")\n", + "markers = [\"o\", \"^\", \"s\", \"D\", \"v\"]\n", + "for (name, estimate), marker in zip(thermometers.items(), markers):\n", + " ax.scatter(T_true[valid], estimate[valid], s=16, alpha=0.65, marker=marker, label=name)\n", + "ax.set_xlabel(\"True grid temperature (°C)\")\n", + "ax.set_ylabel(\"Thermometer temperature (°C)\")\n", + "ax.legend()\n", + "ax.grid(ls=\":\", alpha=0.4)\n", + "plt.show()\n", + "\n", + "print(\"\\nThermometer offsets on valid masked points:\")\n", + "print(f\"{'Calibration':42s} {'Mean':>8s} {'Median':>8s} {'RMSE':>8s}\")\n", + "for name, estimate in thermometers.items():\n", + " offset = estimate[valid] - T_true[valid]\n", + " print(f\"{name:42s} {np.mean(offset):+8.1f} {np.median(offset):+8.1f} {np.sqrt(np.mean(offset**2)):8.1f}\")\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 4))\n", + "names = list(thermometers)\n", + "offsets = [np.mean(thermometers[name][valid] - T_true[valid]) for name in names]\n", + "ax.barh(names, offsets, color=plt.get_cmap(\"viridis\")(np.linspace(0.15, 0.85, len(names))))\n", + "ax.axvline(0.0, color=\"k\", lw=0.8)\n", + "ax.set_xlabel(\"Mean thermometer offset (°C)\")\n", + "ax.set_title(\"Calibration offsets on valid grid points\")\n", + "ax.grid(axis=\"x\", ls=\":\", alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6.5, 5))\n", + "Kd_map = Kd.reshape(Tgrid.shape)\n", + "T_RG_grid = T_RG.reshape(Tgrid.shape)\n", + "\n", + "cs = ax.contourf(Tgrid, Pgrid, T_RG_grid - Tgrid, levels=14, cmap=\"viridis\")\n", + "# ax.contour(Tgrid, Pgrid, Kd_map, levels=14, colors=\"k\", linewidths=0.4)\n", + "cbar = fig.colorbar(cs, ax=ax, label=r\"$\\Delta$T (T$_{RG}$ $-$ T$_{MAGEMin}$)\")\n", + "\n", + "ax.set_xlabel(\"Temperature (°C)\")\n", + "ax.set_ylabel(\"Pressure (kbar)\")\n", + "ax.set_title(\"Fe²⁺–Mg distribution coefficient, S10 eclogite (mpe)\")\n", + "ax.grid(ls=\":\", alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6.5, 5))\n", + "Kd_map = Kd.reshape(Tgrid.shape)\n", + "T_EG_grid = T_EG.reshape(Tgrid.shape)\n", + "\n", + "cs = ax.contourf(Tgrid, Pgrid, T_EG_grid - Tgrid, levels=14, cmap=\"viridis\")\n", + "# ax.contour(Tgrid, Pgrid, Kd_map, levels=14, colors=\"k\", linewidths=0.4)\n", + "cbar = fig.colorbar(cs, ax=ax, label=r\"$\\Delta$T (T$_{RG}$ $-$ T$_{MAGEMin}$)\")\n", + "\n", + "ax.set_xlabel(\"Temperature (°C)\")\n", + "ax.set_ylabel(\"Pressure (kbar)\")\n", + "ax.set_title(\"Fe²⁺–Mg distribution coefficient, S10 eclogite (mpe)\")\n", + "ax.grid(ls=\":\", alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "## 9. Why divalent iron — and why the calibrations disagree\n", + "\n", + "Clinopyroxene can carry substantial Fe³⁺ (in acmite/jadeite-type components), so its Fe²⁺\n", + "value is taken from the excess-O heuristic split. Garnet is intentionally different here: the\n", + "traditional thermometer convention uses $b_0[\\text{ox\\_apfu\\_FeO}]$, i.e. garnet Fe on an FeO\n", + "basis with all garnet Fe treated as Fe²⁺. The resulting garnet-FeO / cpx-heuristic-Fe²⁺ pairing\n", + "is a mixed convention adopted for calibration comparison, not a claim that both phases have\n", + "the same redox treatment." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "from phasetools import get_phase_fe_split\n", + "\n", + "i0 = int(np.flatnonzero(valid)[0]) # first valid grid point\n", + "o0 = out[i0]\n", + "cpx_split = get_phase_fe_split(o0, cpx)\n", + "grt_feo_basis = calc.extract_from_grid(\"g\", oxides=[\"FeO\"], grid_out=[o0])[0][\"ox_apfu_FeO\"][0]\n", + "\n", + "print(f\"At {T_true[i0]:.0f} °C, {P_flat[i0]:.0f} kbar:\")\n", + "print(f\" garnet : FeO-basis Fe = {grt_feo_basis:.3f} (all treated as Fe2+ by convention)\")\n", + "print(f\" {cpx:10s}: Fe2+ = {cpx_split['Fe2']:.3f}, Fe3+ = {cpx_split['Fe3']:.3f} \"\n", + " f\"(Fe3+/FeOt = {cpx_split['Fe3']/(cpx_split['Fe2']+cpx_split['Fe3']):.0%})\")" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "The calibrations disagree because they encode different pressure and composition corrections. Mysen & Heier\n", + "(1972) is pressure- and composition-independent; Ganguly (1979) is retained here only with its\n", + "Yavuz & Yıldırım (2020) piecewise coefficients explicitly labelled; Krogh Ravna (2000) includes\n", + "Mn, Ca, Mg#, and pressure; Ellis & Green (1979) uses binary garnet Ca and pressure; and Räheim &\n", + "Green (1974) uses pressure alone. The offsets therefore illustrate calibration uncertainty, compounded\n", + "by the intentional mixed garnet FeO / cpx heuristic Fe²⁺ convention. These equations should not be\n", + "extrapolated beyond their experimental compositional and P–T ranges." + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "---\n", + "## References\n", + "\n", + "- Mysen, B. O., & Heier, K. S. (1972). Petrogenesis of eclogites in high grade metamorphic\n", + " gneisses, exemplified by the Hareidland eclogite, western Norway. *Contributions to Mineralogy\n", + " and Petrology*, 36(1), 73–94. https://doi.org/10.1007/BF00372836\n", + "- Ganguly, J. (1979). Garnet and clinopyroxene solid solutions, and geothermometry based on\n", + " Fe–Mg distribution coefficient. *Geochimica et Cosmochimica Acta*, 43(7), 1021–1029.\n", + " https://doi.org/10.1016/0016-7037(79)90091-7\n", + "- Yavuz, F., & Yıldırım, D. K. (2020). WinGrt, a Windows program for garnet supergroup minerals.\n", + " *Journal of Geosciences*, 65(2), 71–95. https://doi.org/10.3190/jgeosci.303 (source of the\n", + " Ganguly 1979 piecewise coefficients)\n", + "- Johnson, C. A., Bohlen, S. R., & Essene, E. J. (1983). An evaluation of garnet-clinopyroxene\n", + " geothermometry in granulites. *Contributions to Mineralogy and Petrology*, 84(2–3), 191–198.\n", + " https://doi.org/10.1007/BF00371285\n", + "- Krogh Ravna, E. (2000). The garnet–clinopyroxene Fe²⁺–Mg geothermometer: an updated\n", + " calibration. *Journal of Metamorphic Geology*, 18(2), 211–219.\n", + " https://doi.org/10.1046/j.1525-1314.2000.00247.x\n", + "- Ellis, D. J., & Green, D. H. (1979). An experimental study of the effect of Ca upon garnet–\n", + " clinopyroxene Fe–Mg exchange equilibria. *Contributions to Mineralogy and Petrology*, 71, 13–22.\n", + " https://doi.org/10.1007/BF00371878\n", + "- Räheim, A., & Green, D. H. (1974). Experimental determination of the temperature and pressure\n", + " dependence of the Fe–Mg partition coefficient for coexisting garnet and clinopyroxene.\n", + " *Contributions to Mineralogy and Petrology*, 48, 179–203.\n", + " https://doi.org/10.1007/BF00392328\n", + "- Thomas, J. B., & Rana, S. (2024). Garnet–clinopyroxene thermometry (as discussed alongside\n", + " Yavuz & Yıldırım). Consult the source report for the exact calibration context; this notebook\n", + " does not implement an unverified Powell equation.\n", + " (The Ganguly piecewise coefficients above are attributed explicitly to Yavuz & Yıldırım.)" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "uw3-dev", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Tutorials/Garnet/README.md b/Tutorials/Garnet/README.md index 55fd92b..d8a391a 100644 --- a/Tutorials/Garnet/README.md +++ b/Tutorials/Garnet/README.md @@ -21,5 +21,14 @@ This directory contains tutorials focused on modelling garnet chemistry, growth, - Implementing geothermobarometry by minimising the misfit between measured chemistry and thermodynamic predictions. - Recovering P-T trajectories from zoned garnet crystals. +### [3. Garnet-Cpx Thermometry](./3-Garnet_cpx_thermometry.ipynb) +**Objective:** Mapping the garnet–clinopyroxene Fe²⁺–Mg distribution coefficient (Kd) and comparing thermometer calibrations. +- **Key Features:** + - Setting up the S10 eclogite bulk composition in the `mpe` database. + - Extracting garnet and clinopyroxene Fe²⁺/Mg from APFU with dynamic cpx phase resolution (`dio`/`omph`/`aug`). + - Masking absent and solvus-ambiguous grid cells. + - Mapping Kd over a P–T grid and applying Mysen & Heier (1972), Ganguly (1979; explicitly labelled Yavuz & Yıldırım piecewise form), Krogh Ravna (2000), Ellis & Green (1979), and Räheim & Green (1974) calibrations. + - Using the intentional traditional mixed Fe convention: garnet Fe from `ox_apfu_FeO` (all treated as Fe²⁺) and cpx Fe²⁺ from the excess-O heuristic split. + --- **Units Note:** All tutorials use `kbar` for pressure and `°C` for temperature. diff --git a/src/phasetools/calculators/garnet.py b/src/phasetools/calculators/garnet.py index 4565236..7e594d3 100644 --- a/src/phasetools/calculators/garnet.py +++ b/src/phasetools/calculators/garnet.py @@ -51,12 +51,18 @@ def _extract_garnet_elements_from_oxides(self, out, sys_in): return Mg, Mn, Fe, Ca def generate_2D_grid_gt_endmembers(self, P, T): - """Compute garnet endmember fractions over a P-T grid.""" + """Compute garnet endmember fractions over a P-T grid. + + Garnet is a single-instance phase, so its single per-instance + bundle is returned directly, preserving the historical key names + (``em_py``, ``em_alm``, ``mol_frac``, ...). If garnet ever + appeared as a solvus, the list of bundles would be returned + instead. + """ self.calculate_grid(P, T) # Automatic discovery of end-members res = self.extract_from_grid("g", end_members='auto') - - return res + return res[0] if len(res) == 1 else res def generate_2D_grid_gt_elements(self, P, T): """Compute garnet element fractions (Mg, Mn, Fe, Ca) over a P-T grid.""" diff --git a/src/phasetools/calculators/pt_grid.py b/src/phasetools/calculators/pt_grid.py index f5d1f94..7fee7ec 100644 --- a/src/phasetools/calculators/pt_grid.py +++ b/src/phasetools/calculators/pt_grid.py @@ -1,8 +1,9 @@ import numpy as np import sys +from juliacall import Main as jl, convert as jlconvert from ..core.base import MAGEMinBase -from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number -from ..utils.bulk_rock import atomic_mass_dict, convert_mol_percent_to_wt_percent +from ..core.phase_properties import extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, _phase_indices +from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac from phasetools import MAGEMin_C class MAGEMinPTGridCalculator(MAGEMinBase): @@ -20,17 +21,19 @@ def calculate_grid(self, P, T): P = np.atleast_1d(P) T = np.atleast_1d(T) - if P.shape != T.shape: - if P.ndim == 1 and T.ndim == 1: - P_orig, T_orig = P, T - P, T = np.meshgrid(P_orig, T_orig) - P = P.flatten() - T = T.flatten() - else: - raise ValueError(f"P and T must have the same shape or both be 1D. Got {P.shape} and {T.shape}") + if P.ndim == 1 and T.ndim == 1 and P.shape[0] != T.shape[0]: + P_orig, T_orig = P, T + P, T = np.meshgrid(P_orig, T_orig) + P = P.flatten() + T = T.flatten() + elif P.shape != T.shape: + raise ValueError(f"P and T must have the same shape or both be 1D. Got {P.shape} and {T.shape}") + + P_jl = jlconvert(jl.Vector[jl.Float64], P.astype(float)) + T_jl = jlconvert(jl.Vector[jl.Float64], T.astype(float)) out = MAGEMin_C.multi_point_minimization( - P, T, self.data, X=self.X, Xoxides=self.Xoxides, + P_jl, T_jl, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list ) sys.stdout.flush() @@ -77,38 +80,44 @@ def get_phase_endmembers(self, phase, grid_out=None): if out is None: return [] - # If it's a single output object (from single_point_calc) - if not isinstance(out, (list, np.ndarray)): + # If it's a single output object (from single_point_calc) — not a list/array/Vector + if not isinstance(out, (list, np.ndarray)) and hasattr(out, 'ph'): if phase in out.ph: ph_index = out.ph.index(phase) - return [str(n) for n in out.SS_vec[ph_index].emNames] + if ph_index < len(out.SS_vec): + return [str(n) for n in out.SS_vec[ph_index].emNames] return [] - # If it's a grid + # If it's a grid (list, numpy array, or Julia Vector) for o in out: if phase in o.ph: ph_index = o.ph.index(phase) - return [str(n) for n in o.SS_vec[ph_index].emNames] + if ph_index < len(o.SS_vec): + return [str(n) for n in o.SS_vec[ph_index].emNames] + return [] return [] def _extract_cations_from_apfu(self, out, phase, cations, sys_in): - """Internal: Extract cation ratios (e.g., XMg, XFe) for a specific phase.""" + """Internal: per-instance cation ratios (e.g., XMg, XFe). + + Returns ``cat_`` arrays, one entry per instance of the phase. + """ ox_to_query = ['MgO', 'MnO', 'CaO', 'FeO', 'Fe2O3', 'Fe', 'O'] - apfu = get_oxide_apfu(out, phase, ox_to_query) - - mg = apfu.get("MgO", 0.0) - mn = apfu.get("MnO", 0.0) - ca = apfu.get("CaO", 0.0) - - feo = apfu.get("FeO", 0.0) - fe2o3 = apfu.get("Fe2O3", 0.0) - fe_metal = apfu.get("Fe", 0.0) - atomic_o = apfu.get("O", 0.0) + apfu = get_oxide_apfu(out, phase, ox_to_query, instance='all') + + mg = np.asarray(apfu.get("MgO", 0.0), dtype=float) + mn = np.asarray(apfu.get("MnO", 0.0), dtype=float) + ca = np.asarray(apfu.get("CaO", 0.0), dtype=float) + + feo = np.asarray(apfu.get("FeO", 0.0), dtype=float) + fe2o3 = np.asarray(apfu.get("Fe2O3", 0.0), dtype=float) + fe_metal = np.asarray(apfu.get("Fe", 0.0), dtype=float) + atomic_o = np.asarray(apfu.get("O", 0.0), dtype=float) # Calculate total Fe as FeO equivalent (molar atoms) - if atomic_o > 0: + if np.any(atomic_o > 0): # MAGEMin O-basis (ig, mp) or sb24 basis - if fe_metal > 0: + if np.any(fe_metal > 0): fe = fe_metal else: fe = feo @@ -117,21 +126,45 @@ def _extract_cations_from_apfu(self, out, phase, cations, sys_in): fe = feo + 2.0 * fe2o3 total = mg + mn + fe + ca - if total <= 0: - return {c: 0.0 for c in cations} + total_safe = np.where(total > 0, total, np.nan) + with np.errstate(divide='ignore', invalid='ignore'): + vals = {"Mg": mg/total_safe, "Mn": mn/total_safe, + "Fe": fe/total_safe, "Ca": ca/total_safe} - vals = {"Mg": mg/total, "Mn": mn/total, "Fe": fe/total, "Ca": ca/total} - if sys_in.casefold() == 'wt': - keys = list(vals.keys()) - raw_vals = [vals[k] for k in keys] - wt_percents = convert_mol_percent_to_wt_percent(raw_vals, keys, atomic_mass_dict) - vals = {k: v/100.0 for k, v in zip(keys, wt_percents)} + vals = atomic_frac_to_wt_frac(vals, atomic_mass_dict) - return {f"cat_{c}": vals.get(c, 0.0) for c in cations} + return {f"cat_{c}": np.asarray(vals[c], dtype=float) for c in cations} def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False, grid_out=None): - """Extract phase properties from a previously calculated grid.""" + """Extract phase properties from a previously calculated grid. + + MAGEMin reports coexisting solvus limbs as repeated entries in + ``out.ph`` (e.g. two clinopyroxenes ``dio``). Returns a list + with one bundle per instance, indexed ``res[0]``, ``res[1]`` + ...: + + ``res[0]`` -- first (or only) instance; every key is an array + over the grid points (scalars for ``single_point_calc``). + + The key shape is identical for every phase -- a phase with a + single instance simply has a one-element list, always at + ``res[0]``. Sum the per-instance ``mol_frac`` arrays yourself + if you want the total. + + Each bundle key is NaN where the phase (or that instance) is + absent at a grid point. If the phase never appears in the grid, + the returned list is empty. + + Notes + ----- + - Bundle ``res[k]`` is the k-th occurrence of the phase in that + point's ``out.ph``; solvus branch ordering may swap between + grid points, so track both limbs when plotting isopleths. + - ``end_members='auto'`` discovers end-members from the first + occurrence of the phase; solvus limbs normally share the same + solution model and endmember set. + """ out = grid_out if grid_out is not None else self.last_grid_out if out is None: raise ValueError("No grid results found. Run calculate_grid first or provide grid_out.") @@ -141,50 +174,77 @@ def extract_from_grid(self, phase, end_members=None, oxides=None, chemistry=None end_members = self.get_phase_endmembers(phase, out) P_len = len(out) - results = { - "mol_frac": np.zeros(P_len), "wt_frac": np.zeros(P_len), "vol_frac": np.zeros(P_len), - } - + + # n_inst = max occurrences of the phase across the whole grid + n_inst = 0 + for o in out: + n_inst = max(n_inst, len(_phase_indices(o, phase, 'all'))) + + instances = [{} for _ in range(n_inst)] + + def _precreate(prefix, names): + for k in range(n_inst): + for n in names: + instances[k][f"{prefix}{n}"] = np.full(P_len, np.nan) + + for k in range(n_inst): + instances[k]["mol_frac"] = np.full(P_len, np.nan) + instances[k]["wt_frac"] = np.full(P_len, np.nan) + instances[k]["vol_frac"] = np.full(P_len, np.nan) if end_members: - for em in end_members: results[f"em_{em}"] = np.zeros(P_len) + _precreate("em_", end_members) if oxides: - for ox in oxides: results[f"ox_apfu_{ox}"] = np.zeros(P_len) + _precreate("ox_apfu_", oxides) if chemistry: - for ox in chemistry: results[f"chem_{ox}"] = np.zeros(P_len) + _precreate("chem_", chemistry) if cations: - for c in cations: results[f"cat_{c}"] = np.zeros(P_len) + _precreate("cat_", cations) if mg_number: - results["mg_number"] = np.zeros(P_len) + for k in range(n_inst): + instances[k]["Mg_number"] = np.full(P_len, np.nan) if fe_split: - results["fe2"] = np.zeros(P_len) - results["fe3"] = np.zeros(P_len) + for k in range(n_inst): + instances[k]["Fe2"] = np.full(P_len, np.nan) + instances[k]["Fe3"] = np.full(P_len, np.nan) for i in range(P_len): - if phase in out[i].ph: - results["mol_frac"][i] = phase_frac(phase, out[i], 'mol') - results["wt_frac"][i] = phase_frac(phase, out[i], 'wt') - results["vol_frac"][i] = phase_frac(phase, out[i], 'vol') - - if end_members: - for em in end_members: - results[f"em_{em}"][i] = extract_end_member(phase, out[i], em, self.sys_in) - if oxides: - apfu = get_oxide_apfu(out[i], phase, oxides) - for ox in oxides: results[f"ox_apfu_{ox}"][i] = apfu.get(ox, 0.0) - if chemistry: - chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in) - for ox in chemistry: results[f"chem_{ox}"][i] = chem.get(ox, 0.0) - if cations: - cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in) - for c in cations: results[f"cat_{c}"][i] = cat_vals[f"cat_{c}"] - if mg_number: - results["mg_number"][i] = get_phase_mg_number(out[i], phase) - if fe_split: - split = self._extract_fe_split_from_apfu(out[i], phase) - results["fe2"][i] = split["fe2"] - results["fe3"][i] = split["fe3"] - - return results + if phase not in out[i].ph: + continue + n_idx = _phase_indices(out[i], phase, 'all') + for k, j in enumerate(n_idx): + instances[k]["mol_frac"][i] = float(out[i].ph_frac[j]) + instances[k]["wt_frac"][i] = float(out[i].ph_frac_wt[j]) + instances[k]["vol_frac"][i] = float(out[i].ph_frac_vol[j]) + # instances k >= len(n_idx) stay NaN + + def _fill(key, value): + vals = np.atleast_1d(np.asarray(value, dtype=float)) + n = len(vals) + for k in range(n_inst): + if k < n: + instances[k][key][i] = vals[k] + # else stays NaN + + if end_members: + for em in end_members: + _fill(f"em_{em}", extract_end_member(phase, out[i], em, self.sys_in, instance='all')) + if oxides: + apfu = get_oxide_apfu(out[i], phase, oxides, instance='all') + for ox in oxides: _fill(f"ox_apfu_{ox}", apfu.get(ox, np.zeros(0))) + if chemistry: + chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in, instance='all') + for ox in chemistry: _fill(f"chem_{ox}", chem.get(ox, np.zeros(0))) + if cations: + cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in) + for c in cations: _fill(f"cat_{c}", cat_vals[f"cat_{c}"]) + if mg_number: + _fill("Mg_number", get_phase_mg_number(out[i], phase, instance='all')) + if fe_split: + split = self._extract_fe_split_from_apfu(out[i], phase, instance='all') + _fill("Fe2", split["Fe2"]) + _fill("Fe3", split["Fe3"]) + + return instances def generate_2D_grid(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False): """Convenience wrapper.""" @@ -192,31 +252,47 @@ def generate_2D_grid(self, P, T, phase, end_members=None, oxides=None, chemistry return self.extract_from_grid(phase, end_members, oxides, chemistry, cations, mg_number, fe_split) def single_point_calc(self, P, T, phase, end_members=None, oxides=None, chemistry=None, cations=None, mg_number=False, fe_split=False): - """Single-point calculation.""" + """Single-point calculation. + + Returns ``(bundles, out)`` where ``bundles`` is a list with one + bundle per instance of the phase (``bundles[0]``, ``bundles[1]``, + ...), keyed like the grid-level ``extract_from_grid`` with scalar + values. The list is empty when the phase is not present at this + P-T. + """ out = MAGEMin_C.single_point_minimization(P, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list) sys.stdout.flush() - results = {"mol_frac": 0.0, "wt_frac": 0.0, "vol_frac": 0.0, "present": False} + instances = [] if phase in out.ph: - results["present"] = True - results["mol_frac"] = phase_frac(phase, out, 'mol') - results["wt_frac"] = phase_frac(phase, out, 'wt') - results["vol_frac"] = phase_frac(phase, out, 'vol') + n_idx = _phase_indices(out, phase, 'all') + instances = [{"mol_frac": 0.0, "wt_frac": 0.0, "vol_frac": 0.0} + for _ in n_idx] + for k, j in enumerate(n_idx): + instances[k]["mol_frac"] = float(out.ph_frac[j]) + instances[k]["wt_frac"] = float(out.ph_frac_wt[j]) + instances[k]["vol_frac"] = float(out.ph_frac_vol[j]) + + def _store(key, value): + vals = np.atleast_1d(np.asarray(value, dtype=float)) + for k in range(len(vals)): + instances[k][key] = vals[k] if end_members: - for em in end_members: results[f"em_{em}"] = extract_end_member(phase, out, em, self.sys_in) + for em in end_members: + _store(f"em_{em}", extract_end_member(phase, out, em, self.sys_in, instance='all')) if oxides: - apfu = get_oxide_apfu(out, phase, oxides) - for ox in oxides: results[f"ox_apfu_{ox}"] = apfu.get(ox, 0.0) + apfu = get_oxide_apfu(out, phase, oxides, instance='all') + for ox in oxides: _store(f"ox_apfu_{ox}", apfu.get(ox, np.zeros(0))) if chemistry: - chem = get_phase_chemistry(out, phase, chemistry, self.sys_in) - for ox in chemistry: results[f"chem_{ox}"] = chem.get(ox, 0.0) + chem = get_phase_chemistry(out, phase, chemistry, self.sys_in, instance='all') + for ox in chemistry: _store(f"chem_{ox}", chem.get(ox, np.zeros(0))) if cations: cat_vals = self._extract_cations_from_apfu(out, phase, cations, self.sys_in) - for c in cations: results[f"cat_{c}"] = cat_vals[f"cat_{c}"] + for c in cations: _store(f"cat_{c}", cat_vals[f"cat_{c}"]) if fe_split: - split = self._extract_fe_split_from_apfu(out, phase) - results["fe2"] = split["fe2"] - results["fe3"] = split["fe3"] - - return results, out + split = self._extract_fe_split_from_apfu(out, phase, instance='all') + _store("Fe2", split["Fe2"]) + _store("Fe3", split["Fe3"]) + + return instances, out diff --git a/src/phasetools/core/README.md b/src/phasetools/core/README.md index c87f88e..2aa41f3 100644 --- a/src/phasetools/core/README.md +++ b/src/phasetools/core/README.md @@ -22,3 +22,43 @@ The `core` submodule provides the foundational classes and low-level bridging lo - `get_phase_mg2_number`: Phase-wide $Mg\#$ (Divalent Iron only). - `get_phase_fe_split`: Heuristic splitting of total iron into $\text{Fe}^{2+}$ and $\text{Fe}^{3+}$. - `calculate_kd_fe_mg`: Distribution coefficients between phases. + +### Solvus (multi-instance) handling +MAGEMin reports coexisting solvus limbs as repeated entries in `out.ph` +(e.g. two `dio` clinopyroxenes or two `amp` amphiboles). `phase_frac` +sums them; the low-level composition helpers in `phase_properties.py` +(`get_oxide_apfu`, `get_phase_chemistry`, `extract_end_member`, +`get_phase_mg_number`, `get_phase_mg2_number`, `get_phase_fe_split`) +take an `instance` argument: + +- integer index (default `0`) — that instance, with a `UserWarning` + noting how many other instances exist; out-of-range → zeros; +- `'all'` — one value per instance, returned as numpy arrays. + +`MAGEMinPTGridCalculator.extract_from_grid` / `single_point_calc` / +`generate_2D_grid` return a **list of per-instance bundles**, indexed +`res[0]`, `res[1]`, ... : + +```python +res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) +c0, c1 = res # limb 0, limb 1 +c0['mol_frac'] # array over grid points (limb 0) +c0['ox_apfu_Na2O'] # array over grid points (limb 0) +total_dio = c0['mol_frac'] + c1['mol_frac'] # sum the limbs yourself +``` + +- Bundle keys are unsuffixed (`mol_frac`, `wt_frac`, `vol_frac`, + `ox_apfu_*`, `chem_*`, `em_*`, `cat_*`, `Mg_number`, `Fe2`, `Fe3`). +- The shape is the **same for every phase** — a single-instance phase + (e.g. garnet) has a one-element list, always at `res[0]`. +- Each bundle value is an array over the grid points, **NaN** where the + phase (or that limb) is absent; if the phase never appears, the list + is empty. No totals are computed — sum the `mol_frac` bundles if + you want the whole-phase fraction. +- Bundle `res[k]` is the k-th occurrence of the phase in that point's + `out.ph`; solvus branch ordering may swap between grid points, so + track both limbs when plotting isopleths. `end_members='auto'` + discovers end-members from the first occurrence (solvus limbs share + the same solution model). `MAGEMinGarnetCalculator. + generate_2D_grid_gt_endmembers` returns the single-instance bundle + directly, preserving the historical `em_py`, `mol_frac`, ... keys. diff --git a/src/phasetools/core/base.py b/src/phasetools/core/base.py index 8a13bde..d81bd75 100644 --- a/src/phasetools/core/base.py +++ b/src/phasetools/core/base.py @@ -75,9 +75,9 @@ def setup_bulk_composition(self, Xoxides, X, sys_in, rm_list=None): except: pass - def _extract_fe_split_from_apfu(self, out, phase): + def _extract_fe_split_from_apfu(self, out, phase, instance=0): """ Internal: Calculate Fe2+ and Fe3+ amounts for a phase using an excess oxygen heuristic. """ from .phase_properties import get_phase_fe_split - return get_phase_fe_split(out, phase) + return get_phase_fe_split(out, phase, instance=instance) diff --git a/src/phasetools/core/phase_properties.py b/src/phasetools/core/phase_properties.py index f54efe8..4a78141 100644 --- a/src/phasetools/core/phase_properties.py +++ b/src/phasetools/core/phase_properties.py @@ -1,24 +1,93 @@ import numpy as np +import warnings -def get_oxide_apfu(out, ph, oxides): - """Extract oxide amounts from APFU output for a specific phase.""" - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - oxide_values = np.array(phase_obj.Comp_apfu, dtype=float) - oxide_names = [str(ox) for ox in out.oxides] - oxide_moles_dict = {ox: value for ox, value in zip(oxide_names, oxide_values)} +def _phase_indices(out: object, phase: str, instance: int | str = 0) -> list[int]: + """Return the index/indices of ``phase`` in ``out.ph``. + + Parameters + ---------- + out : object + MAGEMin output object with a ``ph`` attribute. + phase : str + Phase name. + instance : int or {'all'}, default=0 + Which instance to resolve. An integer index selects that instance + (negative indices count from the end). ``'all'`` returns every + occurrence. When the phase appears more than once (a solvus -- e.g. + two coexisting clinopyroxenes or amphiboles, reported by MAGEMin as + repeated entries in ``out.ph``), a warning notes how many other + instances exist. + + Returns + ------- + list[int] + Indices of ``phase`` in ``out.ph`` (empty list if the phase is + absent or the requested instance does not exist). + """ + idx = [i for i, p in enumerate(out.ph) if str(p) == phase] + if instance == 'all': + return idx + if not isinstance(instance, (int, np.integer)): + raise ValueError(f"instance must be an integer index or 'all', got {instance!r}") + if not idx: + return [] + if instance < 0: + instance = len(idx) + instance + if instance < 0 or instance >= len(idx): + warnings.warn( + f"phase '{phase}' has {len(idx)} instance(s); requested instance " + f"{instance} does not exist -- treating as absent.", + UserWarning, stacklevel=3) + return [] + if len(idx) > 1: + warnings.warn( + f"phase '{phase}' has {len(idx)} instances (solvus) in this " + f"output; using instance {instance} (there are {len(idx) - 1} " + f"others). Use instance='all' for per-instance arrays.", + UserWarning, stacklevel=3) + return [idx[instance]] + +def get_oxide_apfu(out, ph, oxides, instance=0): + """Extract oxide amounts from APFU output for a specific phase. - results = {} - for oxide in oxides: - results[oxide] = oxide_moles_dict.get(oxide, 0.0) - except (ValueError, IndexError): - results = {oxide: 0.0 for oxide in oxides} + Parameters + ---------- + out : object + MAGEMin output object. + ph : str + Phase name. + oxides : list[str] + Oxides to extract. + instance : int or {'all'}, default=0 + For a phase that appears multiple times (a solvus), an integer + index selects that instance (warning if more than one exists) and + ``'all'`` returns one value per instance as numpy arrays. - return results + Returns + ------- + dict + ``{oxide: value}`` for an integer index, or ``{oxide: ndarray}`` + (one entry per phase instance) for ``'all'``. + """ + idx = _phase_indices(out, ph, instance) + if not idx: + return {oxide: 0.0 for oxide in oxides} if instance != 'all' \ + else {oxide: np.zeros(0) for oxide in oxides} + oxide_names = [str(ox) for ox in out.oxides] + per = [] + for i in idx: + try: + values = np.array(out.SS_vec[i].Comp_apfu, dtype=float) + per.append({ox: float(values[oxide_names.index(ox)]) + if ox in oxide_names else 0.0 for ox in oxides}) + except (ValueError, IndexError, AttributeError): + per.append({ox: 0.0 for ox in oxides}) + if instance == 'all': + return {ox: np.array([d[ox] for d in per]) for ox in oxides} + return per[0] -def get_phase_chemistry(out, ph, oxides, sys_in): +def get_phase_chemistry(out, ph, oxides, sys_in, instance=0): """ Extract oxide concentrations (wt% or mol%) for a specific phase. @@ -32,46 +101,61 @@ def get_phase_chemistry(out, ph, oxides, sys_in): List of oxides to extract. sys_in : str Unit system ('wt' or 'mol'). + instance : int or {'all'}, default=0 + For a phase that appears multiple times (a solvus), an integer + index selects that instance (warning if more than one exists) and + ``'all'`` returns one value per instance as numpy arrays. Returns ------- dict - Oxide concentrations. + Oxide concentrations for an integer index, or per-instance arrays + for ``'all'``. """ - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - - if sys_in.casefold() == 'wt': - # Comp_wt is weight fraction (0-1) for oxides in the phase - values = np.array(phase_obj.Comp_wt, dtype=float) * 100.0 - else: - # Comp is molar fraction (0-1) for oxides in the phase - values = np.array(phase_obj.Comp, dtype=float) * 100.0 - - oxide_names = [str(ox) for ox in out.oxides] - oxide_dict = {ox: value for ox, value in zip(oxide_names, values)} - - results = {} - for oxide in oxides: - results[oxide] = oxide_dict.get(oxide, 0.0) - except (ValueError, IndexError): - results = {oxide: 0.0 for oxide in oxides} - - return results - -def extract_end_member(phase, MAGEMinOutput, end_member, sys_in): - """Extract specific end-member fraction from MAGEMin output.""" - try: - phase_ind = MAGEMinOutput.ph.index(phase) - em_index = MAGEMinOutput.SS_vec[phase_ind].emNames.index(end_member) - if sys_in.casefold() == 'wt': - data = MAGEMinOutput.SS_vec[phase_ind].emFrac_wt[em_index] - else: - data = MAGEMinOutput.SS_vec[phase_ind].emFrac[em_index] - except (ValueError, IndexError): - data = 0. - return data + idx = _phase_indices(out, ph, instance) + if not idx: + return {oxide: 0.0 for oxide in oxides} if instance != 'all' \ + else {oxide: np.zeros(0) for oxide in oxides} + oxide_names = [str(ox) for ox in out.oxides] + per = [] + for i in idx: + try: + phase_obj = out.SS_vec[i] + if sys_in.casefold() == 'wt': + # Comp_wt is weight fraction (0-1) for oxides in the phase + values = np.array(phase_obj.Comp_wt, dtype=float) * 100.0 + else: + # Comp is molar fraction (0-1) for oxides in the phase + values = np.array(phase_obj.Comp, dtype=float) * 100.0 + per.append({ox: float(values[oxide_names.index(ox)]) + if ox in oxide_names else 0.0 for ox in oxides}) + except (ValueError, IndexError, AttributeError): + per.append({ox: 0.0 for ox in oxides}) + if instance == 'all': + return {ox: np.array([d[ox] for d in per]) for ox in oxides} + return per[0] + +def extract_end_member(phase, MAGEMinOutput, end_member, sys_in, instance=0): + """Extract specific end-member fraction from MAGEMin output. + + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as a numpy array. + """ + idx = _phase_indices(MAGEMinOutput, phase, instance) + if not idx: + return 0.0 if instance != 'all' else np.zeros(0) + vals = [] + for i in idx: + try: + em_index = MAGEMinOutput.SS_vec[i].emNames.index(end_member) + if sys_in.casefold() == 'wt': + vals.append(float(MAGEMinOutput.SS_vec[i].emFrac_wt[em_index])) + else: + vals.append(float(MAGEMinOutput.SS_vec[i].emFrac[em_index])) + except (ValueError, IndexError, AttributeError): + vals.append(0.0) + return float(vals[0]) if instance != 'all' else np.array(vals) def phase_frac(phase, MAGEMinOutput, sys_in): """ @@ -97,7 +181,7 @@ def phase_frac(phase, MAGEMinOutput, sys_in): except: return 0.0 -def get_phase_mg_number(out, ph): +def get_phase_mg_number(out, ph, instance=0): """ Calculate Mg# (molar Mg / (Mg + Fe_total)) for a specific phase. @@ -107,110 +191,97 @@ def get_phase_mg_number(out, ph): Matches the logic used by MAGEMin's 'ss_MgNum' mode by pulling MgO and FeO directly from the phase's Comp_apfu array. Supports 'FeO', 'Fe' (sb24), and 'Fe2O3' fallback. + + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as a numpy array. """ - try: - ph_index = out.ph.index(ph) - phase_obj = out.SS_vec[ph_index] - - oxide_names = [str(ox) for ox in out.oxides] - - # Pull Mg + idx = _phase_indices(out, ph, instance) + if not idx: + return 0.0 if instance != 'all' else np.zeros(0) + oxide_names = [str(ox) for ox in out.oxides] + + def _mg(i): try: - mg_idx = oxide_names.index('MgO') - mg = float(phase_obj.Comp_apfu[mg_idx]) - except ValueError: - mg = 0.0 - - # Pull Fe (total iron atoms) - fe = 0.0 - if 'FeO' in oxide_names: - fe_idx = oxide_names.index('FeO') - fe += float(phase_obj.Comp_apfu[fe_idx]) - # If both are present, we sum them (though unlikely in standard MAGEMin output) - if 'Fe2O3' in oxide_names: - fe2o3_idx = oxide_names.index('Fe2O3') - fe += 2.0 * float(phase_obj.Comp_apfu[fe2o3_idx]) - elif 'Fe' in oxide_names: - fe_idx = oxide_names.index('Fe') - fe += float(phase_obj.Comp_apfu[fe_idx]) - elif 'Fe2O3' in oxide_names: - fe2o3_idx = oxide_names.index('Fe2O3') - fe += 2.0 * float(phase_obj.Comp_apfu[fe2o3_idx]) - else: - # No iron components found - if mg == 0: return 0.0 - return 1.0 # Pure Mg phase - - denominator = mg + fe - if denominator == 0: + phase_obj = out.SS_vec[i] + mg = float(phase_obj.Comp_apfu[oxide_names.index('MgO')]) if 'MgO' in oxide_names else 0.0 + fe = 0.0 + if 'FeO' in oxide_names: + fe += float(phase_obj.Comp_apfu[oxide_names.index('FeO')]) + if 'Fe2O3' in oxide_names: + fe += 2.0 * float(phase_obj.Comp_apfu[oxide_names.index('Fe2O3')]) + elif 'Fe' in oxide_names: + fe += float(phase_obj.Comp_apfu[oxide_names.index('Fe')]) + elif 'Fe2O3' in oxide_names: + fe += 2.0 * float(phase_obj.Comp_apfu[oxide_names.index('Fe2O3')]) + else: + return 1.0 if mg > 0 else 0.0 + denom = mg + fe + return mg / denom if denom else 0.0 + except (ValueError, IndexError, AttributeError): return 0.0 - return mg / denominator - except (ValueError, IndexError, AttributeError): - return 0.0 + vals = [_mg(i) for i in idx] + return float(vals[0]) if instance != 'all' else np.array(vals) -def get_phase_fe_split(out, ph): +def get_phase_fe_split(out, ph, instance=0): """ Calculate Fe2+ and Fe3+ amounts for a phase using an excess oxygen heuristic. Works for both traditional FeO-Fe2O3 bases and MAGEMin's O-basis (ig, mp). + + For a phase that appears multiple times (a solvus), ``instance=0`` + returns the first instance (warning if more than one exists) and + ``instance='all'`` returns one value per instance as numpy arrays. """ - try: - ox_to_query = ['FeO', 'Fe2O3', 'Fe', 'O'] - apfu = get_oxide_apfu(out, ph, ox_to_query) - - feo_val = apfu.get("FeO", 0.0) - fe2o3_val = apfu.get("Fe2O3", 0.0) - fe_metal_val = apfu.get("Fe", 0.0) - atomic_o = apfu.get("O", 0.0) - - # 1. Calculate Total Fe atoms (Atoms per formula unit) - if atomic_o > 0: - # MAGEMin O-basis (ig, mp) or sb24 basis - # If Fe component is present (sb24), use it; otherwise FeO is total iron. - if fe_metal_val > 0: - total_fe = fe_metal_val - else: - total_fe = feo_val - else: - # Traditional FeO/Fe2O3 basis - total_fe = feo_val + 2.0 * fe2o3_val - - # 2. Calculate Fe3+ atoms using excess oxygen heuristic - # excess_o identifies oxygen atoms added beyond the stoichiometric baseline. - # Works for both 'O as total oxygen' and 'O as excess oxygen' components. - excess_o = max(atomic_o - np.round(atomic_o), 0.0) - fe3 = 2.0 * fe2o3_val + 2.0 * excess_o - - # 3. Divalent iron is the remainder - fe2 = max(total_fe - fe3, 0.0) + apfu = get_oxide_apfu(out, ph, ['FeO', 'Fe2O3', 'Fe', 'O'], instance=instance) - return { - "fe2": fe2, - "fe3": fe3, - } - except: - return {"fe2": 0.0, "fe3": 0.0} + feo = np.asarray(apfu.get("FeO", 0.0), dtype=float) + fe2o3 = np.asarray(apfu.get("Fe2O3", 0.0), dtype=float) + fem = np.asarray(apfu.get("Fe", 0.0), dtype=float) + ato = np.asarray(apfu.get("O", 0.0), dtype=float) + + if ato.size == 0: + empty = {"Fe2": np.zeros(0), "Fe3": np.zeros(0)} + return {"Fe2": 0.0, "Fe3": 0.0} if instance != 'all' else empty + + # 1. Calculate Total Fe atoms (Atoms per formula unit) + # Use np.where so each instance follows its own basis (O-bearing vs traditional). + total_fe = np.where( + ato > 0, + np.where(fem > 0, fem, feo), + feo + 2.0 * fe2o3 + ) + + # 2. Calculate Fe3+ atoms using excess oxygen heuristic + excess_o = np.maximum(ato - np.floor(ato), 0.0) + fe3 = 2.0 * fe2o3 + 2.0 * excess_o + + # 3. Divalent iron is the remainder + fe2 = np.maximum(np.asarray(total_fe) - fe3, 0.0) -def get_phase_mg2_number(out, ph): + if instance != 'all': + return {"Fe2": float(fe2.item()), "Fe3": float(fe3.item())} + return {"Fe2": fe2, "Fe3": fe3} + +def get_phase_mg2_number(out: object, ph: str, instance: int | str = 0) -> float | np.ndarray: """ Calculate Mg# (molar Mg / (Mg + Fe2+)) for a specific phase. Uses an excess oxygen heuristic to split total iron into Fe2+ and Fe3+. """ try: - apfu = get_oxide_apfu(out, ph, ['MgO']) + apfu = get_oxide_apfu(out, ph, ['MgO'], instance=instance) mg = apfu.get('MgO', 0.0) - - split = get_phase_fe_split(out, ph) - fe2 = split['fe2'] - - denominator = mg + fe2 - if denominator == 0: - return 0.0 - - return mg / denominator - except: + + split = get_phase_fe_split(out, ph, instance=instance) + fe2 = split['Fe2'] + + denominator = np.asarray(mg + fe2, dtype=float) + if instance != 'all': + return float(mg / denominator) if denominator > 0 else 0.0 + return np.where(denominator > 0, mg / denominator, 0.0) + except Exception: return 0.0 def calculate_kd_fe_mg(out, phase1, phase2, use_fe2_only=False): diff --git a/src/phasetools/models/garnet_growth.py b/src/phasetools/models/garnet_growth.py index 5516cc3..bad72e9 100644 --- a/src/phasetools/models/garnet_growth.py +++ b/src/phasetools/models/garnet_growth.py @@ -3,6 +3,7 @@ from scipy.interpolate import interp1d from scipy.stats import norm from typing import Any +from juliacall import Main as jl, convert as jlconvert from ..calculators.garnet import MAGEMinGarnetCalculator def generate_distribution(n_classes, r_min, dr, fnr, Gn, tGn): @@ -42,8 +43,8 @@ class GarnetGenerator(MAGEMinGarnetCalculator): Generate synthetic garnet populations with compositional zoning along P-T-t paths. """ - def __init__(self, db="mpe", dataset=636, verbose=False): - super().__init__(db, dataset, verbose) + def __init__(self, db="mpe", dataset=636, verbose=False, fe_basis="FeOt"): + super().__init__(db, dataset, verbose, fe_basis=fe_basis) def setup_bulk_composition(self, Xoxides, X, sys_in, rm_list=None): super().setup_bulk_composition(Xoxides, X, sys_in, rm_list) @@ -79,6 +80,16 @@ def generate_garnet_data(self, normalise_start : bool, default=True If True, the initial garnet volume is set to 0 and only new growth is modeled. If False, the initial thermodynamic volume is used as the starting point. + + Notes + ----- + When both ``fractionate=True`` and ``normalise_start=False``, the + initial garnet fraction is removed from the reactive bulk at the first + P-T point (overstepped nucleation). Subsequent growth increments are + then fractionated from the depleted bulk. When + ``normalise_start=True``, only new growth beyond the initial fraction + is modelled and fractionated; the initial garnet is treated as a + non-reactive seed. """ @@ -288,13 +299,26 @@ def get_prograde_concentrations(self, new_t=None): def get_retrograde_concentrations(self, new_t=None): """Get the retrograde concentrations of garnet-forming elements. - Parameters: - new_t (array-like, optional): New time values to interpolate the data - and return the concentrations at these times. If None, the original - data is returned. Uses a linear interpolation between datapoints. - - Returns: - Concentrations (array): An array with the element concentrations and PTt data at each retrograde step. + Parameters + ---------- + new_t : array-like, optional + New time values to interpolate the data and return the + concentrations at these times. If ``None``, the original data is + returned. Uses linear interpolation between datapoints. + + Returns + ------- + numpy.ndarray + Array with element concentrations and PTt data at each retrograde + step. Rows are ``[t, T, P, Mn, Mg, Fe, Ca]``. + + Notes + ----- + Retrograde concentrations are recalculated using the bulk composition + at the **last growth** step, not the bulk at the end of the P-T path. + This means garnet resorption during retrograde is **not** accounted + for — the retrograde calculation assumes a fixed bulk equal to the + last-growth composition. """ GVi = np.array(self.gt_vol_frac) @@ -334,9 +358,10 @@ def get_retrograde_concentrations(self, new_t=None): Ca_eval = np.zeros_like(t_eval, dtype=float) for i in range(len(t_eval)): + x_jl = jlconvert(jl.Vector[jl.Float64], x_last_growth) (_gt_frac, _gt_wt, _gt_vol, - Mg_i, Mn_i, Fe_i, Ca_i, _out) = self.gt_single_point_calc_elements( - P_eval[i], T_eval[i], self.data, x_last_growth, self.Xoxides, self.sys_in, self.rm_list + Mg_i, Mn_i, Fe_i, Ca_i, _out) = self._gt_single_point_from_jl( + P_eval[i], T_eval[i], x_jl, self.Xoxides, self.sys_in, self.rm_list ) Mn_eval[i] = Mn_i Mg_eval[i] = Mg_i @@ -505,8 +530,8 @@ def plot_garnet_summary(self, size_dist='N', garnet_no=0, path=None, plot_fig=Tr axs[0, 0].set_title('Garnet size') axs[0, 0].set_xlabel('r') axs[0, 0].set_ylabel('f') - axs[0, 0].set_xlim([r_r.min(), r_r.max()]) - axs[0, 0].plot(r_r, finp, '-', label='Size Distribution') + axs[0, 0].set_xlim([r.min(), r.max()]) + axs[0, 0].plot(r, finp, '-', label='Size Distribution') for i in range(n_classes): axs[0, 0].plot([r_r[i], r_r[i]], [0, finp[i]], '-') axs[0, 0].legend() diff --git a/src/phasetools/utils/README.md b/src/phasetools/utils/README.md index 94aa1c2..ccd6c91 100644 --- a/src/phasetools/utils/README.md +++ b/src/phasetools/utils/README.md @@ -7,6 +7,8 @@ The `utils` submodule contains helper functions for chemistry, math, and general ### `bulk_rock.py` - Core stoichiometric engine for managing bulk compositions. - Provides molar mass lookups and unit conversion helpers (e.g., mol% to wt% and vice versa). +- `mol_fractions_to_wt_fractions` / `wt_fractions_to_mol_fractions`: non-normalising converters for single components or subsets (e.g., FeO/Fe2O3/FeOt iron redox conversions). +- `split_feot_to_feo_o` / `express_bulk_in_feo_o_basis`: split total iron (FeOt) into the MAGEMin **FeO + O** redox basis at a target Fe³⁺/FeOt fraction, **conserving FeOt** (FeO column = total Fe; O = f·FeOt/2 per 2FeO + O → Fe₂O₃). Useful for redox sweeps and for converting measured (FeOt-only) bulks into MAGEMin format. - Standardises oxide lists and manages cation-oxide mapping. ### `general.py` diff --git a/src/phasetools/utils/bulk_rock.py b/src/phasetools/utils/bulk_rock.py index 752d683..00da5f2 100644 --- a/src/phasetools/utils/bulk_rock.py +++ b/src/phasetools/utils/bulk_rock.py @@ -44,6 +44,25 @@ def convert_mol_percent_to_wt_percent(mol_percents, components, mass_dict): wt_percents = [(mol * mass_dict[comp] / total_mass) * 100 for comp, mol in zip(components, mol_percents)] return wt_percents +def atomic_frac_to_wt_frac(atomic_frac: dict[str, float], mass_dict: dict[str, float]) -> dict[str, float]: + """Convert atomic (molar) site fractions to weight-based site fractions. + + Parameters + ---------- + atomic_frac : dict + Mapping of component names to their atomic fractions (summing to 1.0). + mass_dict : dict + Mapping of component names to their atomic/molecular masses. + + Returns + ------- + dict + Weight-based fractions (summing to 1.0). + """ + total_mass = sum(atomic_frac[k] * mass_dict[k] for k in atomic_frac) + return {k: atomic_frac[k] * mass_dict[k] / total_mass for k in atomic_frac} + + def convert_wt_percent_to_mol_percent(wt_percents, components, mass_dict): """Generic conversion from weight (mass) percent to mole percent.""" total_moles = 0 @@ -52,6 +71,81 @@ def convert_wt_percent_to_mol_percent(wt_percents, components, mass_dict): mol_percents = [((wt / mass_dict[comp]) / total_moles) * 100 for comp, wt in zip(components, wt_percents)] return mol_percents +def mol_fractions_to_wt_fractions( + mol: float | list | np.ndarray, components: list[str], mass_dict: dict[str, float] +) -> float | list | np.ndarray: + """Convert mole fractions to weight fractions (no normalisation). + + Unlike :func:`convert_mol_percent_to_wt_percent`, the output is not + normalised to sum to 100. This means the function can be applied to a + single component (e.g. a single oxide for an iron redox conversion) or + a subset of a composition, as well as a full composition. + + Parameters + ---------- + mol : float or array_like + Mole fraction(s) of each component. + components : list of str + Component names corresponding to each input value. + mass_dict : dict + Mapping of component names to molecular masses. + + Returns + ------- + float or list or numpy.ndarray + Weight fraction(s) of each component. A scalar input returns a + scalar, a list input returns a list, and any other array-like + input returns a numpy array. + + Notes + ----- + Scalar input assumes a single component -- pass ``components=[comp]`` + with the oxide name. + """ + if np.isscalar(mol): + return float(mol) * mass_dict[components[0]] + mol_arr = np.asarray(mol, dtype=float) + if mol_arr.ndim == 0: + return float(mol_arr.item() * mass_dict[components[0]]) + masses = np.array([mass_dict[comp] for comp in components], dtype=float) + result = mol_arr * masses + return result.tolist() if isinstance(mol, list) else result + +def wt_fractions_to_mol_fractions( + wt: float | list | np.ndarray, components: list[str], mass_dict: dict[str, float] +) -> float | list | np.ndarray: + """Convert weight fractions to mole fractions (no normalisation). + + Unlike :func:`convert_wt_percent_to_mol_percent`, the output is not + normalised to sum to 100. This means the function can be applied to a + single component (e.g. a single oxide for an iron redox conversion) or + a subset of a composition, as well as a full composition. + + Parameters + ---------- + wt : float or array_like + Weight fraction(s) of each component. + components : list of str + Component names corresponding to each input value. + mass_dict : dict + Mapping of component names to molecular masses. + + Returns + ------- + float or list or numpy.ndarray + Mole fraction(s) of each component. A scalar input returns a + scalar, a list input returns a list, and any other array-like + input returns a numpy array. + """ + if np.isscalar(wt): + return float(wt) / mass_dict[components[0]] + wt_arr = np.asarray(wt, dtype=float) + if wt_arr.ndim == 0: + return float(wt_arr.item() / mass_dict[components[0]]) + masses = np.array([mass_dict[comp] for comp in components], dtype=float) + result = wt_arr / masses + return result.tolist() if isinstance(wt, list) else result + def convert_wt_percent_to_moles(wt_percents, components, mass_dict, total_weight): """Convert weight (mass) percentages to moles.""" moles = [] @@ -76,3 +170,119 @@ def convert_moles_to_mol_percent(moles, components): total = sum(moles_dict.values()) return {comp: (moles_dict[comp] / total) * 100 for comp in components} + +def split_feot_to_feo_o(feot_moles: float, fe3_frac: float) -> tuple[float, float]: + """ + Split total iron (FeOt) into the MAGEMin ``FeO + O`` redox pair at a + target Fe3+/FeOt fraction, conserving the total iron budget. + + Parameters + ---------- + feot_moles : float + Total iron in mole units (atoms of Fe, equivalently the amount of + FeO that would carry all iron as Fe2+). + fe3_frac : float + Target Fe3+/FeOt fraction in ``[0, 1]``. 0 = fully reduced (all + Fe2+), 1 = fully oxidised (all Fe3+). + + Returns + ------- + (feo, o) : tuple[float, float] + ``feo`` is the total-iron column (all Fe expressed as FeO, equal to + ``feot_moles``) and ``o`` is the excess oxygen required to oxidise + the target fraction, per ``2FeO + O -> Fe2O3``. + + Notes + ----- + Molar bookkeeping (each Fe2O3 carries 2 Fe atoms and needs 1 O): + Fe3+ atoms = 2 * O => O = fe3_frac * FeOt / 2 + Fe2+ atoms = FeOt - 2 * O (implicitly held by the FeO column) + """ + fe3_frac = float(fe3_frac) + if not 0.0 <= fe3_frac <= 1.0: + raise ValueError(f"fe3_frac must be in [0, 1], got {fe3_frac!r}") + feot_moles = float(feot_moles) + return feot_moles, fe3_frac * feot_moles / 2.0 + +def express_bulk_in_feo_o_basis( + X: list[float], + Xoxides: list[str], + fe3_frac: float, + feo_oxide: str = "FeO", + fe2o3_oxide: str = "Fe2O3", + o_oxide: str = "O", +) -> tuple[list[float], list[str]]: + """ + Express a bulk composition in the MAGEMin ``FeO + O`` redox basis at a + target Fe3+/FeOt fraction, conserving total iron. + + The ``FeO`` column is set to total iron (all Fe expressed as FeO, + ``FeOt``) and ``O`` is set to the excess oxygen giving the requested + Fe3+/FeOt partition. Any ``Fe2O3`` component is removed (set to 0) and + an existing ``O`` component is overwritten; missing components are + appended. The returned list preserves the input oxide order. + + Parameters + ---------- + X : array-like + Bulk composition values in MOLE units (mol fractions or mol%). Use + weight-based converters first if the input is in wt%. + Xoxides : list of str + Oxide names corresponding to ``X`` (may include 'FeO', 'Fe2O3', 'O' + in any combination). + fe3_frac : float + Target Fe3+/FeOt fraction in ``[0, 1]``. + feo_oxide, fe2o3_oxide, o_oxide : str + Component names used in ``Xoxides``. + + Returns + ------- + (X_new, Xoxides_new) : tuple[list, list] + Composition in the ``FeO + O`` basis (unnormalised -- pass to + ``convertBulk4MAGEMin`` or normalise afterwards). + + Examples + -------- + >>> X = [50.0, 8.0, 0.5] # SiO2, FeO, O (mol%) + >>> ox = ['SiO2', 'FeO', 'O'] + >>> X2, ox2 = express_bulk_in_feo_o_basis(X, ox, fe3_frac=0.1) + >>> ox2 + ['SiO2', 'FeO', 'O'] + >>> X2[2] == 0.1 * X[1] / 2.0 # O = fe3_frac * FeOt / 2 + True + """ + if len(Xoxides) != len(X): + raise ValueError( + f"X ({len(X)} items) and Xoxides ({len(Xoxides)} items) must have the same length" + ) + seen = set() + for ox in Xoxides: + if ox in seen: + raise ValueError(f"Duplicate oxide name in Xoxides: {ox!r}") + seen.add(ox) + + X = [float(v) for v in X] + Xoxides = list(Xoxides) + feo_i = Xoxides.index(feo_oxide) if feo_oxide in Xoxides else None + fe2o3_i = Xoxides.index(fe2o3_oxide) if fe2o3_oxide in Xoxides else None + o_i = Xoxides.index(o_oxide) if o_oxide in Xoxides else None + + feo_mol = X[feo_i] if feo_i is not None else 0.0 + fe2o3_mol = X[fe2o3_i] if fe2o3_i is not None else 0.0 + feot = feo_mol + 2.0 * fe2o3_mol # total Fe atoms (moles) + + feo_total, o_excess = split_feot_to_feo_o(feot, fe3_frac) + + if feo_i is not None: + X[feo_i] = feo_total + else: + X.append(feo_total) + Xoxides.append(feo_oxide) + if fe2o3_i is not None: + X[fe2o3_i] = 0.0 + if o_i is not None: + X[o_i] = o_excess + else: + X.append(o_excess) + Xoxides.append(o_oxide) + return X, Xoxides diff --git a/tests/test_solvus_instances.py b/tests/test_solvus_instances.py new file mode 100644 index 0000000..ac35947 --- /dev/null +++ b/tests/test_solvus_instances.py @@ -0,0 +1,365 @@ +"""Mock-based tests for solvus (multi-instance phase) handling. + +MAGEMin reports coexisting solvus limbs as repeated entries in ``out.ph`` +(e.g. two clinopyroxenes ``dio``, or two amphiboles ``amp``). + +Two layers are tested: + +* the low-level composition helpers (``get_oxide_apfu``, + ``extract_end_member``, ...) take ``instance`` (integer index or + ``'all'`` for per-instance numpy arrays); +* the grid-level ``extract_from_grid`` / ``single_point_calc`` return a + list of per-instance bundles indexed ``res[0]``, ``res[1]`` ... Bundle + keys are unsuffixed (``mol_frac``, ``ox_apfu_Na2O``, ``em_py``, ...), + and each value is an array over the grid points (NaN where the phase + or that limb is absent). The shape is the same for every phase; a + single-instance phase just has a one-element list. + +No live Julia runtime is needed -- ``out`` objects are mocked. +""" + +import unittest +import warnings +import numpy as np +from unittest.mock import MagicMock, patch + +from phasetools.calculators.pt_grid import MAGEMinPTGridCalculator +from phasetools.core.phase_properties import ( + get_oxide_apfu, get_phase_chemistry, extract_end_member, + get_phase_mg_number, get_phase_mg2_number, get_phase_fe_split, + _phase_indices, +) + +# mpe-style oxide ordering +OXIDES = ['H2O', 'SiO2', 'Al2O3', 'CaO', 'MgO', 'FeO', 'K2O', 'Na2O', + 'TiO2', 'MnO', 'O'] + + +def _make_ss(apfu, em_names, em_frac): + """Build a mocked SS_vec entry (one phase instance).""" + p = MagicMock() + p.Comp_apfu = apfu + # molar fractions (0-1) aligned with OXIDES + p.Comp = [0.0, 0.6, 0.03, 0.15, 0.13, 0.02, 0.0, 0.06, 0.0, 0.01, 0.0] + p.Comp_wt = [0.0, 0.55, 0.05, 0.12, 0.10, 0.02, 0.0, 0.05, 0.0, 0.01, 0.0] + p.emNames = em_names + p.emFrac = em_frac + p.emFrac_wt = [f * 0.9 for f in em_frac] + return p + + +def _dio_0(): + # diopside-rich limb: low Na + apfu = [0.0, 1.98, 0.05, 0.74, 0.64, 0.12, 0.0, 0.26, 0.0, 0.01, 6.0] + return _make_ss(apfu, ['di', 'hed', 'om', 'jac'], [0.60, 0.10, 0.25, 0.05]) + + +def _dio_1(): + # omphacite-rich limb: high Na + apfu = [0.0, 1.96, 0.09, 0.59, 0.51, 0.18, 0.0, 0.41, 0.0, 0.01, 6.0] + return _make_ss(apfu, ['di', 'hed', 'om', 'jac'], [0.40, 0.10, 0.40, 0.10]) + + +def _garnet(): + apfu = [0.0, 3.02, 2.0, 0.55, 0.70, 1.70, 0.0, 0.0, 0.0, 0.03, 12.0] + return _make_ss(apfu, ['py', 'alm', 'gr', 'spss'], [0.20, 0.60, 0.16, 0.04]) + + +def _make_out(ph_list, ss_vec, ph_frac=None): + out = MagicMock() + out.ph = ph_list + out.oxides = OXIDES + n = len(ph_list) + if ph_frac is None: + ph_frac = [0.1] * n + out.ph_frac = ph_frac + out.ph_frac_wt = list(ph_frac) + out.ph_frac_vol = list(ph_frac) + out.SS_vec = ss_vec + return out + + +class TestPhaseIndices(unittest.TestCase): + """_phase_indices resolution.""" + + def setUp(self): + self.out = _make_out(['dio', 'q', 'dio', 'g'], [_dio_0(), None, _dio_1(), _garnet()]) + + def test_all_returns_every_occurrence(self): + self.assertEqual(_phase_indices(self.out, 'dio', 'all'), [0, 2]) + + def test_default_zero_warns_with_others(self): + with self.assertWarnsRegex(UserWarning, r"2 instances.*1.*others"): + self.assertEqual(_phase_indices(self.out, 'dio', 0), [0]) + + def test_single_instance_does_not_warn(self): + with warnings.catch_warnings(): + warnings.simplefilter("error") + self.assertEqual(_phase_indices(self.out, 'g', 0), [3]) + + def test_out_of_range_warns_and_empty(self): + with self.assertWarnsRegex(UserWarning, r"instance 2 does not exist"): + self.assertEqual(_phase_indices(self.out, 'dio', 2), []) + + def test_absent_phase_empty(self): + self.assertEqual(_phase_indices(self.out, 'ep', 0), []) + + def test_invalid_instance_type_raises(self): + with self.assertRaises(ValueError): + _phase_indices(self.out, 'dio', 'first') + + +class TestSolvusHelpers(unittest.TestCase): + """Per-instance behaviour of the composition helpers.""" + + def setUp(self): + self.out = _make_out(['dio', 'q', 'dio', 'g'], [_dio_0(), None, _dio_1(), _garnet()]) + + def test_oxide_apfu_default_is_first(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO']) + self.assertAlmostEqual(apfu['Na2O'], 0.26) + self.assertAlmostEqual(apfu['MgO'], 0.64) + + def test_oxide_apfu_second_instance(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO'], instance=1) + self.assertAlmostEqual(apfu['Na2O'], 0.41) + self.assertAlmostEqual(apfu['MgO'], 0.51) + + def test_oxide_apfu_all_returns_arrays(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O', 'MgO'], instance='all') + self.assertTrue(np.allclose(apfu['Na2O'], [0.26, 0.41])) + self.assertTrue(np.allclose(apfu['MgO'], [0.64, 0.51])) + + def test_phase_chemistry_second_instance(self): + chem = get_phase_chemistry(self.out, 'dio', ['Na2O'], 'mol', instance=1) + self.assertAlmostEqual(chem['Na2O'], 0.06 * 100.0) # Comp[7] * 100 + + def test_extract_end_member_per_instance(self): + self.assertAlmostEqual(extract_end_member('dio', self.out, 'di', 'mol'), 0.60) + self.assertAlmostEqual(extract_end_member('dio', self.out, 'di', 'mol', instance=1), 0.40) + vals = extract_end_member('dio', self.out, 'di', 'mol', instance='all') + self.assertTrue(np.allclose(vals, [0.60, 0.40])) + + def test_mg_number_per_instance(self): + # Mg# = Mg/(Mg+Fe); first limb Mg-rich, second more Fe-rich + self.assertGreater(get_phase_mg_number(self.out, 'dio'), + get_phase_mg_number(self.out, 'dio', instance=1)) + vals = get_phase_mg_number(self.out, 'dio', instance='all') + self.assertEqual(vals.shape, (2,)) + self.assertAlmostEqual(vals[0], get_phase_mg_number(self.out, 'dio')) + + def test_mg2_number_per_instance(self): + m0 = get_phase_mg2_number(self.out, 'dio') + m1 = get_phase_mg2_number(self.out, 'dio', instance=1) + self.assertGreater(m0, m1) + vals = get_phase_mg2_number(self.out, 'dio', instance='all') + self.assertEqual(vals.shape, (2,)) + + def test_fe_split_per_instance(self): + # O-basis: excess O = 0, so Fe2 = total Fe, Fe3 = 0 for both limbs + s0 = get_phase_fe_split(self.out, 'dio') + s1 = get_phase_fe_split(self.out, 'dio', instance=1) + self.assertAlmostEqual(s0['Fe2'], 0.12) + self.assertAlmostEqual(s1['Fe2'], 0.18) + all_s = get_phase_fe_split(self.out, 'dio', instance='all') + self.assertTrue(np.allclose(all_s['Fe2'], [0.12, 0.18])) + + def test_out_of_range_returns_zeros(self): + apfu = get_oxide_apfu(self.out, 'dio', ['Na2O'], instance=5) + self.assertEqual(apfu['Na2O'], 0.0) + + +class TestExtractFromGridSolvus(unittest.TestCase): + """Grid extraction: list of per-instance bundles (res[0], res[1], ...).""" + + def _make_calc(self, grid_out): + calc = MAGEMinPTGridCalculator.__new__(MAGEMinPTGridCalculator) + calc.sys_in = 'mol' + calc.last_grid_out = grid_out + calc.rm_list = None + return calc + + def _dio_grid(self): + # point 0: dio solvus (2 limbs); point 1: single dio + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + g1 = _make_out(['q', 'g', 'dio'], + [None, _garnet(), _dio_0()], + ph_frac=[0.5, 0.2, 0.3]) + return [g0, g1] + + def test_uniform_shape_single_and_solvus(self): + """A single-instance phase and a solvus return the same bundle keys.""" + calc = self._make_calc(self._dio_grid()) + solvus = calc.extract_from_grid('dio', oxides=['Na2O', 'MgO'], + mg_number=True, fe_split=True) + single = calc.extract_from_grid('g', oxides=['MgO'], + mg_number=True, fe_split=True) + + # dio: two bundles; garnet: one bundle (always at res[0]) + self.assertEqual(len(solvus), 2) + self.assertEqual(len(single), 1) + + for bundle in solvus + single: + for key in ('mol_frac', 'wt_frac', 'vol_frac', + 'ox_apfu_Na2O' if 'ox_apfu_Na2O' in bundle else 'ox_apfu_MgO', + 'Mg_number', 'Fe2', 'Fe3'): + self.assertIn(key, bundle) + # unsuffixed keys inside each bundle (no _0/_1, no total_*) + self.assertIn('mol_frac', solvus[0]) + self.assertNotIn('mol_frac_0', solvus[0]) + self.assertNotIn('total_mol_frac', solvus[0]) + + def test_limb_bundles(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', oxides=['Na2O']) + c0, c1 = res + # limb fractions (0.3/0.1 at point 0; point 1 has only limb 0) + self.assertTrue(np.allclose(c0['mol_frac'], [0.3, 0.3])) + self.assertAlmostEqual(c0['wt_frac'][0], 0.3) + self.assertAlmostEqual(c0['vol_frac'][0], 0.3) + self.assertAlmostEqual(c1['mol_frac'][0], 0.1) + self.assertAlmostEqual(c0['ox_apfu_Na2O'][0], 0.26) + self.assertAlmostEqual(c1['ox_apfu_Na2O'][0], 0.41) + # user can sum the limbs themselves + total = c0['mol_frac'] + c1['mol_frac'] + self.assertAlmostEqual(total[0], 0.4) + + def test_nan_when_limb_absent(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) + c0, c1 = res + # point 1 has a single dio -> limb-1 values NaN + self.assertTrue(np.isnan(c1['mol_frac'][1])) + self.assertTrue(np.isnan(c1['wt_frac'][1])) + self.assertTrue(np.isnan(c1['ox_apfu_Na2O'][1])) + self.assertTrue(np.isnan(c1['Mg_number'][1])) + # limb 0 is filled + self.assertAlmostEqual(c0['mol_frac'][1], 0.3) + self.assertAlmostEqual(c0['ox_apfu_Na2O'][1], 0.26) + + def test_absent_phase_bundles_nan(self): + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + g2 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc([g0, g2]) + res = calc.extract_from_grid('dio', oxides=['Na2O']) + c0, c1 = res + # point 1: dio absent -> NaN in every bundle + self.assertTrue(np.isnan(c0['mol_frac'][1])) + self.assertTrue(np.isnan(c1['mol_frac'][1])) + self.assertTrue(np.isnan(c0['ox_apfu_Na2O'][1])) + + def test_absent_everywhere_empty(self): + """Phase never stable -> empty list.""" + g0 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + g1 = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.5, 0.5]) + calc = self._make_calc([g0, g1]) + res = calc.extract_from_grid('dio', oxides=['Na2O'], mg_number=True) + self.assertEqual(res, []) + + def test_single_instance_phase_filled(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('g', oxides=['MgO']) + bundle = res[0] + self.assertAlmostEqual(bundle['mol_frac'][0], 0.4) + self.assertAlmostEqual(bundle['ox_apfu_MgO'][0], 0.70) + + def test_cations_per_instance(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', cations=['Mg', 'Fe']) + c0, c1 = res + self.assertIn('cat_Mg', c0) + self.assertIn('cat_Mg', c1) + # first limb is Mg-richer + self.assertGreater(c0['cat_Mg'][0], c1['cat_Mg'][0]) + + def test_end_members_per_instance(self): + calc = self._make_calc(self._dio_grid()) + res = calc.extract_from_grid('dio', end_members=['di', 'om']) + c0, c1 = res + self.assertAlmostEqual(c0['em_di'][0], 0.60) + self.assertAlmostEqual(c1['em_di'][0], 0.40) + self.assertAlmostEqual(c0['em_om'][0], 0.25) + self.assertAlmostEqual(c1['em_om'][0], 0.40) + + +class TestSinglePointCalcSolvus(unittest.TestCase): + """single_point_calc returns the same nested-bundle schema.""" + + def _make_calc(self): + calc = MAGEMinPTGridCalculator.__new__(MAGEMinPTGridCalculator) + calc.sys_in = 'mol' + calc.data = None + calc.X = None + calc.Xoxides = None + calc.rm_list = None + return calc + + def test_solvus_returns_limb_bundles(self): + g0 = _make_out(['dio', 'q', 'dio', 'g'], + [_dio_0(), None, _dio_1(), _garnet()], + ph_frac=[0.3, 0.2, 0.1, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g0 + res, out = calc.single_point_calc(10.0, 600.0, 'dio', + oxides=['Na2O', 'MgO']) + self.assertEqual(len(res), 2) + c0, c1 = res + self.assertAlmostEqual(c0['mol_frac'], 0.3) + self.assertAlmostEqual(c1['mol_frac'], 0.1) + self.assertAlmostEqual(c0['wt_frac'], 0.3) + self.assertAlmostEqual(c1['wt_frac'], 0.1) + self.assertAlmostEqual(c0['ox_apfu_Na2O'], 0.26) + self.assertAlmostEqual(c1['ox_apfu_Na2O'], 0.41) + self.assertIs(out, g0) + + def test_single_instance_has_one_bundle(self): + g = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g + res, _ = calc.single_point_calc(10.0, 600.0, 'g', oxides=['MgO']) + self.assertEqual(len(res), 1) + bundle = res[0] + self.assertAlmostEqual(bundle['mol_frac'], 0.4) + self.assertAlmostEqual(bundle['ox_apfu_MgO'], 0.70) + + def test_absent_phase_no_bundles(self): + g = _make_out(['q', 'g'], [None, _garnet()], ph_frac=[0.6, 0.4]) + calc = self._make_calc() + with patch('phasetools.calculators.pt_grid.MAGEMin_C') as m: + m.single_point_minimization.return_value = g + res, _ = calc.single_point_calc(10.0, 600.0, 'dio') + self.assertEqual(res, []) + + +class TestGarnetEndmembersSuffix(unittest.TestCase): + """generate_2D_grid_gt_endmembers returns the single-instance bundle.""" + + def test_returns_bundle_with_historical_keys(self): + from phasetools.calculators.garnet import MAGEMinGarnetCalculator + calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator) + bundle = { + 'mol_frac': np.array([0.10, 0.15]), + 'wt_frac': np.array([0.11, 0.16]), + 'vol_frac': np.array([0.12, 0.17]), + 'em_py': np.array([0.20, 0.30]), + 'em_alm': np.array([0.60, 0.50]), + 'em_gr': np.array([0.16, 0.17]), + 'em_spss': np.array([0.04, 0.03]), + } + with patch.object(calc, 'calculate_grid', return_value=None), \ + patch.object(calc, 'extract_from_grid', return_value=[bundle]): + res = calc.generate_2D_grid_gt_endmembers([10.0], [600.0]) + self.assertIn('em_py', res) + self.assertIn('em_alm', res) + self.assertIn('mol_frac', res) + self.assertTrue(np.allclose(res['em_py'], [0.20, 0.30])) + + +if __name__ == '__main__': + unittest.main()