docs(tutorials): garnet-cpx thermometry + fix: garnet_growth radius arrays

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.

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."""

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.

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)