pydrex.minerals

  1"""> PyDRex: Computations of mineral texture and elasticity.
  2
  3**Acronyms:**
  4- CPO = Crystallographic preferred orientation,
  5    i.e. preferential clustering of polycrystal grain orientations in SO(3),
  6    leading to an overall anisotropic orientation distribution
  7
  8"""
  9
 10import io
 11from dataclasses import dataclass, field
 12from zipfile import ZipFile
 13
 14import numpy as np
 15from numpy import asarray as Ŋ
 16from scipy import linalg as la
 17from scipy.integrate import LSODA
 18from scipy.spatial.transform import Rotation
 19
 20from pydrex import core as _core
 21from pydrex import exceptions as _err
 22from pydrex import io as _io
 23from pydrex import logger as _log
 24from pydrex import tensors as _tensors
 25from pydrex import utils as _utils
 26
 27OLIVINE_PRIMARY_AXIS = {
 28    _core.MineralFabric.olivine_A: "a",
 29    _core.MineralFabric.olivine_B: "c",
 30    _core.MineralFabric.olivine_C: "c",
 31    _core.MineralFabric.olivine_D: "a",
 32    _core.MineralFabric.olivine_E: "a",
 33}
 34"""Primary slip axis name for for the given olivine `fabric`."""
 35
 36
 37OLIVINE_SLIP_SYSTEMS = (
 38    ([0, 1, 0], [1, 0, 0]),
 39    ([0, 0, 1], [1, 0, 0]),
 40    ([0, 1, 0], [0, 0, 1]),
 41    ([1, 0, 0], [0, 0, 1]),
 42)
 43"""Slip systems for olivine in conventional order.
 44
 45Tuples contain the slip plane normal and slip direction vectors.
 46The order of slip systems returned matches the order of critical shear stresses
 47returned by `pydrex.core.get_crss`.
 48
 49"""
 50
 51
 52@dataclass
 53class StiffnessTensors:
 54    olivine: np.ndarray = field(
 55        default_factory=lambda: np.array(
 56            [
 57                [320.71, 69.84, 71.22, 0.0, 0.0, 0.0],
 58                [69.84, 197.25, 74.8, 0.0, 0.0, 0.0],
 59                [71.22, 74.8, 234.32, 0.0, 0.0, 0.0],
 60                [0.0, 0.0, 0.0, 63.77, 0.0, 0.0],
 61                [0.0, 0.0, 0.0, 0.0, 77.67, 0.0],
 62                [0.0, 0.0, 0.0, 0.0, 0.0, 78.36],
 63            ]
 64        )
 65    )
 66    """Stiffness tensor for olivine (Voigt representation), with units of GPa.
 67
 68    The source of the values used here is unknown, but they are copied
 69    from the original DRex code: <http://www.ipgp.fr/~kaminski/web_doudoud/DRex.tar.gz>
 70    [88K download]
 71
 72    """
 73    enstatite: np.ndarray = field(
 74        default_factory=lambda: np.array(
 75            [
 76                [236.9, 79.6, 63.2, 0.0, 0.0, 0.0],
 77                [79.6, 180.5, 56.8, 0.0, 0.0, 0.0],
 78                [63.2, 56.8, 230.4, 0.0, 0.0, 0.0],
 79                [0.0, 0.0, 0.0, 84.3, 0.0, 0.0],
 80                [0.0, 0.0, 0.0, 0.0, 79.4, 0.0],
 81                [0.0, 0.0, 0.0, 0.0, 0.0, 80.1],
 82            ]
 83        )
 84    )
 85    """Stiffness tensor for enstatite (Voigt representation), with units of GPa.
 86
 87    The source of the values used here is unknown, but they are copied
 88    from the original DRex code: <http://www.ipgp.fr/~kaminski/web_doudoud/DRex.tar.gz>
 89    [88K download]
 90
 91    """
 92
 93    def __iter__(self):
 94        # So that [S for S in StiffnessTensors()] can be indexed with `MineralPhase`s.
 95        indexed = {
 96            _core.MineralPhase.olivine: self.olivine,
 97            _core.MineralPhase.enstatite: self.enstatite,
 98        }
 99        for _, v in sorted(indexed.items()):
100            yield v
101
102
103def peridotite_solidus(pressure, fit="Hirschmann2000"):
104    """Get peridotite solidus (i.e. melting) temperature based on experimental fits.
105
106    Pressure is expected to be in GPa.
107
108    Supported fits:
109    - ["Hirschmann2000"](https://doi.org/10.1029/2000GC000070)
110    - ["Herzberg2000"](https://doi.org/10.1029/2000GC000089)
111    - ["Duvernay2024"](https://doi.org/10.1029/2023GC011288)
112
113    """
114    match fit:
115        case "Herzberg2000":
116            return 1086 - 5.7 * pressure + 390 * np.log(pressure)
117        case "Hirschmann2000":
118            return -5.104 * pressure**2 + 132.899 * pressure + 1120.661
119        case "Duvernay2024":
120            return -6.8 * pressure**2 + 141.4 * pressure + 1101.3
121        case _:
122            raise ValueError(f"unsupported fit '{fit}'")
123
124
125@dataclass
126class Mineral:
127    """Class for storing polycrystal texture for a single mineral phase.
128
129    A `Mineral` stores texture data for an aggregate of grains¹.
130    Additionally, mineral fabric type and deformation regime are also tracked.
131    To provide an initial texture for the mineral, use the constructor arguments
132    `fractions_init` and `orientations_init`. By default,
133    a uniform volume distribution of random orientations is generated.
134
135    The `update_orientations` method computes new orientations and grain volumes
136    for a given velocity gradient. These results are stored in the `.orientations` and
137    `.fractions` attributes of the `Mineral` instance. The method also returns the
138    updated macroscopic deformation gradient based on the provided initial deformation
139    gradient.
140
141    ¹Note that the "number of grains" is a static integer value that
142    does not track the actual number of physical grains in the deforming polycrystal.
143    Instead, this number acts as a "number of bins" for the statistical resolution of
144    the crystallographic orientation distribution. The value is roughly equivalent to
145    (a multiple of) the number of initial, un-recrystallised grains in the polycrystal.
146    It is assumed that recrystallised grains do not grow large enough to require
147    rotation tracking.
148
149    **Examples:**
150
151    Mineral with isotropic initial texture:
152
153    >>> import pydrex
154    >>> olA = pydrex.Mineral(
155    ...     phase=pydrex.MineralPhase.olivine,
156    ...     fabric=pydrex.MineralFabric.olivine_A,
157    ...     regime=pydrex.DeformationRegime.matrix_dislocation,
158    ...     n_grains=2000
159    ... )
160    >>> olA.phase
161    <MineralPhase.olivine: 0>
162    >>> olA.fabric
163    <MineralFabric.olivine_A: 0>
164    >>> olA.regime
165    <DeformationRegime.matrix_dislocation: 4>
166    >>> olA.n_grains
167    2000
168
169    Mineral with specified initial texture and default phase, fabric and regime settings
170    which are for an olivine A-type mineral in the dislocation creep regime.
171    The initial grain volume fractions should be normalised.
172
173    >>> import numpy as np
174    >>> from scipy.spatial.transform import Rotation
175    >>> import pydrex
176    >>> rng = np.random.default_rng()
177    >>> n_grains = 2000
178    >>> olA = pydrex.Mineral(
179    ...     n_grains=n_grains,
180    ...     fractions_init=np.full(n_grains, 1 / n_grains),
181    ...     orientations_init=Rotation.from_euler(
182    ...         "zxz", [
183    ...             [x * np.pi / 2, np.pi / 2, np.pi / 2] for x in rng.random(n_grains)
184    ...         ]
185    ...     ).inv().as_matrix(),
186    ... )
187    >>> len(olA.orientations)
188    1
189    >>> type(olA.orientations)
190    <class 'list'>
191    >>> olA.orientations[0].shape
192    (2000, 3, 3)
193    >>> olA.fractions[0].shape
194    (2000,)
195
196    Note that minerals can also be constructed from serialized data,
197    see `Mineral.load` and `Mineral.from_file`.
198
199    **Attributes:**
200    - `phase` (`pydrex.core.MineralPhase`) — ordinal number of the mineral phase
201    - `fabric` (`pydrex.core.MineralFabric`) — ordinal number of the fabric type
202    - `regime` (`pydrex.core.DeformationRegime`) — ordinal number of the deformation
203      regime
204    - `n_grains` (int) — number of grains in the aggregate
205    - `fractions` (list of arrays) — grain volume fractions for each texture snapshot
206    - `orientations` (list of arrays) — grain orientation matrices for each texture
207      snapshot
208    - `seed` (int) — seed used by the random number generator to set up the isotropic
209      initial condition when `fractions_init` or `orientations_init` are not provided
210    - `lband` (int) — passed to the `scipy.integrate.LSODA` solver
211    - `uband` (int) — passed to the `scipy.integrate.LSODA` solver
212
213    """
214
215    phase: int = _core.MineralPhase.olivine
216    fabric: int = _core.MineralFabric.olivine_A
217    regime: int = _core.DeformationRegime.matrix_dislocation
218    n_grains: int = _core.DefaultParams().number_of_grains
219    # Initial condition, randomised if not given.
220    fractions_init: np.ndarray | None = None
221    orientations_init: np.ndarray | None = None
222    fractions: list = field(default_factory=list)
223    orientations: list = field(default_factory=list)
224    seed: int | None = None
225    lband: int | None = None
226    uband: int | None = None
227
228    def __str__(self):
229        # String output, used for str(self) and f"{self}", etc.
230        if hasattr(self.fractions[0], "shape"):
231            shape_of_fractions = str(self.fractions[0].shape)
232        else:
233            shape_of_fractions = "(?)"
234
235        if hasattr(self.orientations[0], "shape"):
236            shape_of_orientations = str(self.orientations[0].shape)
237        else:
238            shape_of_orientations = "(?)"
239
240        obj = self.__class__.__qualname__
241        phase = f"(phase={self.phase!r}, "
242        fabric = f"fabric={self.fabric!r}, "
243        regime = f"regime={self.regime!r}, "
244        n_grains = f"n_grains={self.n_grains}, "
245        _fclass = self.fractions.__class__.__qualname__
246        _f0class = self.fractions[0].__class__.__qualname__
247        frac = f"fractions=<{_fclass} of {_f0class} {shape_of_fractions}>, "
248        _oclass = self.orientations.__class__.__qualname__
249        _o0class = self.orientations[0].__class__.__qualname__
250        orient = f"orientations=<{_oclass} of {_o0class} {shape_of_orientations}>)"
251        return f"{obj}{phase}{fabric}{regime}{n_grains}{frac}{orient}"
252
253    def _repr_pretty_(self, p, cycle):
254        # Format to use when printing to IPython or other interactive console.
255        p.text(self.__str__() if not cycle else self.__class__.__qualname__ + "(...)")
256
257    def __post_init__(self):
258        """Initialise random orientations and grain volume fractions."""
259        if self.fractions_init is None:
260            self.fractions_init = np.full(self.n_grains, 1.0 / self.n_grains)
261        if self.orientations_init is None:
262            self.orientations_init = Rotation.random(
263                self.n_grains, random_state=self.seed
264            ).as_matrix()
265        # For large numbers of grains, the number of ODE's exceeds what LSODA can
266        # handle. Therefore, we specify the Jacobian matrix as banded.
267        # By default, we use a bandwidth o f 12000 with lband = uband = 6000.
268        # This should work for up to 10000 grains.
269        if self.lband is None and self.uband is None and self.n_grains > 4632:
270            _log.warning(
271                "using a banded Jacobian because of the large number of grains."
272                + " To manually control the bandwidth, set `lband` and/or `uband`"
273                + f" in calls to `{self.__class__.__qualname__}.update_orientations`."
274            )
275            self.lband = 6000
276            self.uband = 6000
277
278        # Copy the initial values to the storage lists.
279        self.fractions.append(self.fractions_init)
280        self.orientations.append(self.orientations_init)
281
282        # Delete the initial value duplicates to avoid confusion.
283        del self.fractions_init
284        del self.orientations_init
285
286        _log.info("created %s", self)
287
288    def __eq__(self, other):
289        if other.__class__ is self.__class__:
290            return (
291                self.phase == other.phase
292                and self.fabric == other.fabric
293                and self.regime == other.regime
294                and self.n_grains == other.n_grains
295                and len(self.fractions) == len(other.fractions)
296                and np.all(
297                    Ŋ([f.shape for f in self.fractions])
298                    == Ŋ([f.shape for f in other.fractions])
299                )
300                and np.all(Ŋ(self.fractions) == Ŋ(other.fractions))
301                and len(self.orientations) == len(other.orientations)
302                and np.all(
303                    Ŋ([f.shape for f in self.orientations])
304                    == Ŋ([f.shape for f in other.orientations])
305                )
306                and np.all(Ŋ(self.orientations) == Ŋ(other.orientations))
307                and self.seed == other.seed
308                and self.lband == other.lband
309                and self.uband == other.uband
310            )
311        return False
312
313    def update_orientations(
314        self,
315        params: dict,
316        deformation_gradient: np.ndarray,
317        get_velocity_gradient,
318        pathline: tuple,
319        get_regime=None,
320        **kwargs,
321    ) -> np.ndarray:
322        """Update orientations and volume distribution for the `Mineral`.
323
324        Update crystalline orientations and grain volume distribution
325        for minerals undergoing plastic deformation. Return the updated deformation
326        gradient measuring the corresponding macroscopic deformation.
327
328        .. note:: For multi-phase simulations, the return value should only be collected
329            for the last phase, i.e. the last call to `update_orientations`. This value
330            should be passed on to subsequent calls to `update_orientations` as the
331            `deformation_gradient` argument.
332
333        - `params` — dictionary with the same keys as the return value of
334          `pydrex.core.DefaultParams.as_dict()`
335        - `deformation_gradient` — 3x3 deformation gradient tensor
336        - `get_velocity_gradient` — callable with signature `f(t, x)` that returns a 3x3
337          velocity gradient matrix at time `t` and position `x` (3D vector)
338        - `pathline` — tuple consisting of:
339            1. the time at which to start the CPO integration
340            2. the time at which to stop the CPO integration
341            3. a callable with signature `f(t)` that returns the position of the mineral
342            at time `t` where $t ∈ [t_{start}, t_{end}]$
343        - `get_regime` (optional) — callable with signature f(t, x) that returns a
344          `pydrex.core.DeformationRegime` declaring the dominant deformation mechanism
345          at time `t` and position `x` (3D vector)
346
347        Any additional (optional) keyword arguments are passed to
348        `scipy.integrate.LSODA`.
349
350        Array values must provide a NumPy-compatible interface:
351        <https://numpy.org/doc/stable/user/whatisnumpy.html>
352
353        Examples:
354
355        >>> from itertools import pairwise
356        >>> import pydrex
357        >>> olA = pydrex.Mineral(
358        ...           phase=pydrex.MineralPhase.olivine,
359        ...           fabric=pydrex.MineralFabric.olivine_A,
360        ...           regime=pydrex.DeformationRegime.matrix_dislocation,
361        ...           n_grains=pydrex.DefaultParams().number_of_grains,
362        ... )
363        >>> def get_velocity_gradient(t, x):  # Simple L for simple shear.
364        ...     L = np.zeros((3, 3))
365        ...     L[0, 2] = 1
366        ...     return L
367        >>> def get_position(t):
368        ...     return np.zeros(3)  # Stationary polycrystal for this example.
369        >>> timestamps = range(2)  # Just 1 timestep for demonstration.
370        >>> deformation_gradient = np.eye(3)  # Start with an undeformed polycrystal.
371        >>> for t_start, t_end in pairwise(timestamps):
372        ...     # Update deformation_gradient, olA.orientations and olA.fractions.
373        ...     deformation_gradient = olA.update_orientations(
374        ...         pydrex.DefaultParams().as_dict(),
375        ...         deformation_gradient,
376        ...         get_velocity_gradient,
377        ...         (t_start, t_end, get_position),
378        ...     )
379        >>> from numpy import testing as nt
380        >>> nt.assert_allclose(deformation_gradient,
381        ...     np.eye(3) + get_velocity_gradient(t_end, np.nan),
382        ...     atol=1e-15, rtol=0
383        ... )
384
385        """
386
387        # ===== Set up callables for the ODE solver and internal processing =====
388
389        def eval_rhs(t, y):
390            """Evaluate right hand side of the D-Rex PDE."""
391            # assert not np.any(np.isnan(y)), y[np.isnan(y)].shape
392            position = get_position(t)
393            velocity_gradient = get_velocity_gradient(t, position)
394            # _log.debug(
395            #     "calculating CPO at %s (t=%e) with velocity gradient %s",
396            #     position,
397            #     t,
398            #     velocity_gradient.ravel(),
399            # )
400
401            if self.phase not in params["phase_assemblage"]:
402                # Warning rather than failure, so that we tolerate loops like:
403                # [m.update_orientations(...) for m in minerals]
404                # Where some minerals in the list might be (temporarily) superfluous.
405                _log.warning(
406                    "skipping %s mineral, phase omitted from configuration", self.phase
407                )
408                return
409
410            if get_regime is not None:
411                self.regime = get_regime(t, position)
412
413            try:
414                volume_fraction = params["phase_fractions"][
415                    params["phase_assemblage"].index(self.phase)
416                ]
417            except IndexError:
418                raise ValueError(
419                    f"phase must be a valid `MineralPhase`, not {self.phase}"
420                )
421
422            strain_rate = (velocity_gradient + velocity_gradient.transpose()) / 2
423            strain_rate_max = np.abs(la.eigvalsh(strain_rate)).max()
424            deformation_gradient, orientations, fractions = _utils.extract_vars(
425                y, self.n_grains
426            )
427            deformation_gradient_diff = velocity_gradient @ deformation_gradient
428            deformation_gradient_spin = _tensors.polar_decompose(
429                deformation_gradient_diff
430            )[1]
431            # Uses nondimensional values of strain rate and velocity gradient.
432            orientations_diff, fractions_diff = _core.derivatives(
433                regime=self.regime,
434                phase=self.phase,
435                fabric=self.fabric,
436                n_grains=self.n_grains,
437                orientations=orientations,
438                fractions=fractions,
439                strain_rate=strain_rate / strain_rate_max,
440                velocity_gradient=velocity_gradient / strain_rate_max,
441                deformation_gradient_spin=deformation_gradient_spin,
442                stress_exponent=params["stress_exponent"],
443                deformation_exponent=params["deformation_exponent"],
444                nucleation_efficiency=params["nucleation_efficiency"],
445                gbm_mobility=params["gbm_mobility"],
446                volume_fraction=volume_fraction,
447            )
448            return np.hstack(
449                (
450                    deformation_gradient_diff.flatten(),
451                    orientations_diff.flatten() * strain_rate_max,
452                    fractions_diff * strain_rate_max,
453                )
454            )
455
456        def perform_step(solver):
457            """Perform SciPy solver step and appropriate processing."""
458            message = solver.step()
459            if message is not None and solver.status == "failed":
460                raise _err.IterationError(message)
461            # _log.debug(
462            #     "%s step_size=%e", solver.__class__.__qualname__, solver.step_size
463            # )
464
465            deformation_gradient, orientations, fractions = _utils.extract_vars(
466                solver.y, self.n_grains
467            )
468            orientations, fractions = _utils.apply_gbs(
469                orientations,
470                fractions,
471                params["gbs_threshold"],
472                self.orientations[-1],
473                self.n_grains,
474            )
475            solver.y[9:] = np.hstack((orientations.flatten(), fractions))
476
477        # ===== Initialise and run the solver using the above callables =====
478
479        time_start, time_end, get_position = pathline
480        if not callable(get_velocity_gradient):
481            raise ValueError(
482                "unable to evaluate velocity gradient callable."
483                + " You must provide a callable with signature f(t, x)"
484                + " that returns a 3x3 matrix."
485            )
486        if not callable(get_position):
487            raise ValueError(
488                "unable to evaluate position callable."
489                + " You must provide a callable with signature f(t)"
490                + " that returns a 3-component array."
491            )
492        _log.debug(
493            "calculating CPO from %s (t=%s) to %s (t=%s)",
494            get_position(time_start),
495            time_start,
496            get_position(time_end),
497            time_end,
498        )
499        _log.debug("    with deformation gradient %s", deformation_gradient.ravel())
500        _log.debug(
501            "    with velocity gradient interpolated between %s and %s",
502            get_velocity_gradient(time_start, get_position(time_start)).ravel(),
503            get_velocity_gradient(time_end, get_position(time_end)).ravel(),
504        )
505        _log.debug(
506            "    intermediate velocity gradient = %s",
507            get_velocity_gradient(
508                (time_start + time_end) / 2, get_position((time_start + time_end) / 2)
509            ).ravel(),
510        )
511
512        y_start = np.hstack(
513            (
514                deformation_gradient.flatten(),
515                self.orientations[-1].flatten(),
516                self.fractions[-1],
517            )
518        )
519        solver = LSODA(
520            eval_rhs,
521            time_start,
522            y_start,
523            time_end,
524            atol=kwargs.pop("atol", np.abs(y_start * 1e-6) + 1e-12),
525            rtol=kwargs.pop("rtol", 1e-6),
526            first_step=kwargs.pop("first_step", np.abs(time_end - time_start) * 1e-1),
527            # max_step=kwargs.pop("max_step", np.abs(time_end - time_start)),
528            lband=self.lband,
529            uband=self.uband,
530            **kwargs,
531        )
532        perform_step(solver)
533        while solver.status == "running":
534            perform_step(solver)
535
536        # Extract final values for this simulation step, append to storage.
537        deformation_gradient, orientations, fractions = _utils.extract_vars(
538            solver.y.squeeze(), self.n_grains
539        )
540        self.orientations.append(orientations)
541        self.fractions.append(fractions)
542        return deformation_gradient
543
544    def save(self, filename, postfix=None):
545        """Save CPO data for all stored timesteps to a `numpy` NPZ file.
546
547        If `postfix` is not `None`, the data is appended to the NPZ file
548        in fields ending with "`_postfix`".
549
550        Raises a `ValueError` if the data shapes are not compatible.
551
552        See also: `numpy.savez`, `Mineral.load`, `Mineral.from_file`.
553
554        """
555        if len(self.fractions) != len(self.orientations):
556            raise ValueError(
557                "Length of stored results must match."
558                + " You've supplied currupted data with:\n"
559                + f"- {len(self.fractions)} grain size results, and\n"
560                + f"- {len(self.orientations)} orientation results."
561            )
562        if self.fractions[0].shape[0] == self.orientations[0].shape[0] == self.n_grains:
563            data = {
564                "meta": np.array(
565                    [self.phase, self.fabric, self.regime], dtype=np.uint8
566                ),
567                "fractions": np.stack(self.fractions),
568                "orientations": np.stack(self.orientations),
569            }
570            # Create parent directories, resolve relative paths.
571            _io.resolve_path(filename)
572            # Append to file, requires postfix (unique name).
573            if postfix is not None:
574                _log.info("saving Mineral to file %s (postfix: %s)", filename, postfix)
575                archive = ZipFile(filename, mode="a", allowZip64=True)
576                for key in data.keys():
577                    with archive.open(
578                        f"{key}_{postfix}", "w", force_zip64=True
579                    ) as file:
580                        buffer = io.BytesIO()
581                        np.save(buffer, data[key])
582                        file.write(buffer.getvalue())
583                        buffer.close()
584            else:
585                _log.info("saving Mineral to file %s", filename)
586                np.savez(filename, **data)
587        else:
588            raise ValueError(
589                "Size of CPO data arrays must match number of grains."
590                + " You've supplied corrupted data with:\n"
591                + f"- `n_grains = {self.n_grains}`,\n"
592                + f"- `fractions[0].shape = {self.fractions[0].shape}`, and\n"
593                + f"- `orientations[0].shape = {self.orientations[0].shape}`."
594            )
595
596    def load(self, filename, postfix=None):
597        """Load CPO data from a `numpy` NPZ file.
598
599        If `postfix` is not `None`, data is read from fields ending with "`_postfix`".
600
601        See also: `Mineral.save`, `Mineral.from_file`.
602
603        """
604        if not filename.endswith(".npz"):
605            raise ValueError(
606                f"Must only load from numpy NPZ format. Cannot load from {filename}."
607            )
608        data = np.load(filename)
609        if postfix is not None:
610            phase, fabric, regime = data[f"meta_{postfix}"]
611            self.fractions = list(data[f"fractions_{postfix}"])
612            self.orientations = list(data[f"orientations_{postfix}"])
613        else:
614            phase, fabric, regime = data["meta"]
615            self.fractions = list(data["fractions"])
616            self.orientations = list(data["orientations"])
617
618        self.phase = phase
619        self.fabric = fabric
620        self.regime = regime
621        self.orientations_init = self.orientations[0]
622        self.fractions_init = self.fractions[0]
623
624    @classmethod
625    def from_file(cls, filename, postfix=None):
626        """Construct a `Mineral` instance using data from a `numpy` NPZ file.
627
628        If `postfix` is not `None`, data is read from fields ending with “`_postfix`”.
629
630        See also: `Mineral.save`, `Mineral.load`.
631
632        """
633        if not filename.endswith(".npz"):
634            raise ValueError(
635                f"Must only load from numpy NPZ format. Cannot load from {filename}."
636            )
637        data = np.load(filename)
638        if postfix is not None:
639            phase, fabric, regime = data[f"meta_{postfix}"]
640            fractions = list(data[f"fractions_{postfix}"])
641            orientations = list(data[f"orientations_{postfix}"])
642        else:
643            phase, fabric, regime = data["meta"]
644            fractions = list(data["fractions"])
645            orientations = list(data["orientations"])
646
647        mineral = cls(
648            phase,
649            fabric,
650            regime,
651            n_grains=len(fractions[0]),
652            fractions_init=fractions[0],
653            orientations_init=orientations[0],
654        )
655        mineral.fractions = fractions
656        mineral.orientations = orientations
657        return mineral
658
659
660def update_all(
661    minerals: list[Mineral],
662    params: dict,
663    deformation_gradient: np.ndarray,
664    get_velocity_gradient,
665    pathline: tuple,
666    get_regime=None,
667    **kwargs,
668) -> np.ndarray:
669    """Update orientations and volume distributions for all mineral phases.
670
671    Returns the updated deformation gradient tensor which measures the accumulated
672    macroscopic strain.
673
674    """
675    for i, mineral in enumerate(minerals):
676        # Deformation gradient is independent of mineral phase.
677        new_deformation_gradient = mineral.update_orientations(
678            params=params,
679            deformation_gradient=deformation_gradient,
680            get_velocity_gradient=get_velocity_gradient,
681            pathline=pathline,
682            get_regime=get_regime,
683            **kwargs,
684        )
685    return new_deformation_gradient
686
687
688# TODO: Compare to [Man & Huang, 2011](https://doi.org/10.1007/s10659-011-9312-y).
689def voigt_averages(
690    minerals: list[Mineral],
691    phase_assemblage: list[_core.MineralPhase],
692    phase_fractions: list[float],
693    elastic_tensors: StiffnessTensors = StiffnessTensors(),
694):
695    """Calculate elastic tensors as the Voigt averages of a collection of `mineral`s.
696
697    - `minerals` — mineral phases storing orientations and fractional volumes of grains
698    - `phase_assemblage` — collection of unique mineral phases in the aggregate
699    - `phase_fractions` — collection of volume fractions for each phase in
700      `phase_assemblage` (values should sum to 1).
701
702    Raises a ValueError if the minerals contain an unequal number of grains or stored
703    texture results.
704
705    """
706    n_grains = minerals[0].n_grains
707    if not np.all([m.n_grains == n_grains for m in minerals[1:]]):
708        raise ValueError("cannot average minerals with unequal grain counts")
709    n_steps = len(minerals[0].orientations)
710    if not np.all([len(m.orientations) == n_steps for m in minerals[1:]]):
711        raise ValueError(
712            "cannot average minerals with variable-length orientation arrays"
713        )
714    if not np.all([len(m.fractions) == n_steps for m in minerals]):
715        raise ValueError(
716            "cannot average minerals with variable-length grain volume arrays"
717        )
718
719    # TODO: Perform rotation directly on the 6x6 matrices, see Carcione 2007.
720    # This trick is implemented in cpo_elastic_tensor.cc in Aspect.
721    average_tensors = np.zeros((n_steps, 6, 6))
722    for i in range(n_steps):
723        for mineral in minerals:
724            for n in range(n_grains):
725                match mineral.phase:
726                    case _core.MineralPhase.olivine:
727                        average_tensors[i] += _tensors.elastic_tensor_to_voigt(
728                            _tensors.rotate(
729                                elastic_tensors.olivine,
730                                mineral.orientations[i][n, ...].transpose(),
731                            )
732                            * mineral.fractions[i][n]
733                            * phase_fractions[phase_assemblage.index(mineral.phase)]
734                        )
735                    case _core.MineralPhase.enstatite:
736                        average_tensors[i] += _tensors.elastic_tensor_to_voigt(
737                            _tensors.rotate(
738                                elastic_tensors.enstatite,
739                                mineral.orientations[i][n, ...].transpose(),
740                            )
741                            * mineral.fractions[i][n]
742                            * phase_fractions[phase_assemblage.index(mineral.phase)]
743                        )
744                    case _:
745                        raise ValueError(f"unsupported mineral phase: {mineral.phase}")
746    return average_tensors