Merge pull request #200 from seismic-anisotropy/add-update-multiphase

Add multiphase texture update API

CHANGELOG.md

@@ -5,6 +5,28 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+
+## [Unreleased]
+
+### Added
+- `pydrex.update_all` to update texture of multiphase aggregates simultaneously
+- `get_regime` argument to `update_all`/`update_orientations` to allow for
+  temporally variable deformation regimes
+- Texture evolution for diffusion creep and yielding regimes (experimental)
+- Functions to get peridotite solidus and second order tensor invariants
+- Terse SCSV schema parser (only for Python >= 3.12)
+- Ability to toggle verbose doctest output using `pytest -vv`
+
+### Changed
+- Call signature for steady flow 2D box visualisation function
+- Symbol names for default stiffness tensors (now members of
+  `minerals.StiffnessTensors` — use your own preferred stiffness tensors by
+  passing a custom instance when e.g. calculating Voigt averages)
+
+### Fixed
+- Handling of enstatite in Voigt averaging
+
+
 ## [0.0.1] - 2024-04-24
 
 Alpha release.

src/pydrex/__init__.py

@@ -154,6 +154,7 @@
     StiffnessTensors,
     peridotite_solidus,
     voigt_averages,
+    update_all,
 )
 from pydrex.stats import (
     misorientation_hist,

src/pydrex/logger.py

@@ -63,10 +63,12 @@
 # NOTE: Do NOT import any pydrex submodules here to avoid cyclical imports.
 
 np.set_printoptions(
-    formatter={"float_kind": np.format_float_scientific},
+    formatter={
+        "float_kind": np.format_float_scientific,
+        "object": ft.partial(np.array2string, separator=", "),
+    },
     linewidth=1000,
 )
-np.set_string_function(ft.partial(np.array2string, separator=", "), repr=False)
 
 
 class ConsoleFormatter(logging.Formatter):
@@ -118,7 +120,7 @@ def handle_exception(exec_type, exec_value, exec_traceback):
 
 
 @cl.contextmanager
-def handler_level(level: str, handler: logging.Handler=CONSOLE_LOGGER):
+def handler_level(level: str, handler: logging.Handler = CONSOLE_LOGGER):
     """Set logging handler level for current context.
 
     - `level` — logging level name e.g. "DEBUG", "ERROR", etc. See Python's logging

src/pydrex/minerals.py

@@ -24,6 +24,30 @@
 from pydrex import tensors as _tensors
 from pydrex import utils as _utils
 
+OLIVINE_PRIMARY_AXIS = {
+    _core.MineralFabric.olivine_A: "a",
+    _core.MineralFabric.olivine_B: "c",
+    _core.MineralFabric.olivine_C: "c",
+    _core.MineralFabric.olivine_D: "a",
+    _core.MineralFabric.olivine_E: "a",
+}
+"""Primary slip axis name for for the given olivine `fabric`."""
+
+
+OLIVINE_SLIP_SYSTEMS = (
+    ([0, 1, 0], [1, 0, 0]),
+    ([0, 0, 1], [1, 0, 0]),
+    ([0, 1, 0], [0, 0, 1]),
+    ([1, 0, 0], [0, 0, 1]),
+)
+"""Slip systems for olivine in conventional order.
+
+Tuples contain the slip plane normal and slip direction vectors.
+The order of slip systems returned matches the order of critical shear stresses
+returned by `pydrex.core.get_crss`.
+
+"""
+
 
 @dataclass
 class StiffnessTensors:
@@ -76,109 +100,6 @@ def __iter__(self):
             yield v
 
 
-OLIVINE_PRIMARY_AXIS = {
-    _core.MineralFabric.olivine_A: "a",
-    _core.MineralFabric.olivine_B: "c",
-    _core.MineralFabric.olivine_C: "c",
-    _core.MineralFabric.olivine_D: "a",
-    _core.MineralFabric.olivine_E: "a",
-}
-"""Primary slip axis name for for the given olivine `fabric`."""
-
-
-OLIVINE_SLIP_SYSTEMS = (
-    ([0, 1, 0], [1, 0, 0]),
-    ([0, 0, 1], [1, 0, 0]),
-    ([0, 1, 0], [0, 0, 1]),
-    ([1, 0, 0], [0, 0, 1]),
-)
-"""Slip systems for olivine in conventional order.
-
-Tuples contain the slip plane normal and slip direction vectors.
-The order of slip systems returned matches the order of critical shear stresses
-returned by `pydrex.core.get_crss`.
-
-"""
-
-
-def peridotite_solidus(pressure, fit="Hirschmann2000"):
-    """Get peridotite solidus (i.e. melting) temperature based on experimental fits.
-
-    Pressure is expected to be in GPa.
-
-    Supported fits:
-    - ["Hirschmann2000"](https://doi.org/10.1029/2000GC000070)
-    - ["Herzberg2000"](https://doi.org/10.1029/2000GC000089)
-
-    """
-    match fit:
-        case "Herzberg2000":
-            return 1086 - 5.7 * pressure + 390 * np.log(pressure)
-        case "Hirschmann2000":
-            return 5.104 * pressure**2 + 132.899 * pressure + 1120.661
-        case _:
-            raise ValueError("unsupported fit")
-
-
-# TODO: Compare to [Man & Huang, 2011](https://doi.org/10.1007/s10659-011-9312-y).
-def voigt_averages(minerals, phase_assemblage, phase_fractions):
-    """Calculate elastic tensors as the Voigt averages of a collection of `mineral`s.
-
-    - `minerals` — list of `pydrex.minerals.Mineral` instances storing orientations and
-      fractional volumes of the grains within each distinct mineral phase
-    - `phase_assemblage` — collection of `pydrex.core.MineralPhase`s
-    - `phase_fractions` — collection of volume fractions for each phase in
-      `phase_assemblage` (values should sum to 1).
-
-    Raises a ValueError if the minerals contain an unequal number of grains or stored
-    texture results.
-
-    """
-    n_grains = minerals[0].n_grains
-    if not np.all([m.n_grains == n_grains for m in minerals[1:]]):
-        raise ValueError("cannot average minerals with unequal grain counts")
-    n_steps = len(minerals[0].orientations)
-    if not np.all([len(m.orientations) == n_steps for m in minerals[1:]]):
-        raise ValueError(
-            "cannot average minerals with variable-length orientation arrays"
-        )
-    if not np.all([len(m.fractions) == n_steps for m in minerals]):
-        raise ValueError(
-            "cannot average minerals with variable-length grain volume arrays"
-        )
-
-    elastic_tensors = StiffnessTensors()
-
-    # TODO: Perform rotation directly on the 6x6 matrices, see Carcione 2007.
-    # This trick is implemented in cpo_elastic_tensor.cc in Aspect.
-    average_tensors = np.zeros((n_steps, 6, 6))
-    for i in range(n_steps):
-        for mineral in minerals:
-            for n in range(n_grains):
-                match mineral.phase:
-                    case _core.MineralPhase.olivine:
-                        average_tensors[i] += _tensors.elastic_tensor_to_voigt(
-                            _tensors.rotate(
-                                elastic_tensors.olivine,
-                                mineral.orientations[i][n, ...].transpose(),
-                            )
-                            * mineral.fractions[i][n]
-                            * phase_fractions[phase_assemblage.index(mineral.phase)]
-                        )
-                    case _core.MineralPhase.enstatite:
-                        average_tensors[i] += _tensors.elastic_tensor_to_voigt(
-                            _tensors.rotate(
-                                elastic_tensors["enstatite"],
-                                minerals.orientations[i][n, ...].transpose(),
-                            )
-                            * mineral.fractions[i][n]
-                            * phase_fractions[phase_assemblage.index(mineral.phase)]
-                        )
-                    case _:
-                        raise ValueError(f"unsupported mineral phase: {mineral.phase}")
-    return average_tensors
-
-
 @dataclass
 class Mineral:
     """Class for storing polycrystal texture for a single mineral phase.
@@ -375,7 +296,7 @@ def update_orientations(
         pathline: tuple,
         get_regime=None,
         **kwargs,
-    ):
+    ) -> np.ndarray:
         """Update orientations and volume distribution for the `Mineral`.
 
         Update crystalline orientations and grain volume distribution
@@ -712,3 +633,111 @@ def from_file(cls, filename, postfix=None):
         mineral.fractions = fractions
         mineral.orientations = orientations
         return mineral
+
+
+def update_all(
+    minerals: list[Mineral],
+    params: dict,
+    deformation_gradient: np.ndarray,
+    get_velocity_gradient,
+    pathline: tuple,
+    get_regime=None,
+    **kwargs,
+) -> np.ndarray:
+    """Update orientations and volume distributions for all mineral phases.
+
+    Returns the updated deformation gradient tensor which measures the accumulated
+    macroscopic strain.
+
+    """
+    for i, mineral in enumerate(minerals):
+        # Deformation gradient is independent of mineral phase.
+        new_deformation_gradient = mineral.update_orientations(
+            params=params,
+            deformation_gradient=deformation_gradient,
+            get_velocity_gradient=get_velocity_gradient,
+            pathline=pathline,
+            get_regime=get_regime,
+            **kwargs,
+        )
+    return new_deformation_gradient
+
+
+# TODO: Compare to [Man & Huang, 2011](https://doi.org/10.1007/s10659-011-9312-y).
+def voigt_averages(
+    minerals: list[Mineral],
+    phase_assemblage: list[_core.MineralPhase],
+    phase_fractions: list[float],
+    elastic_tensors: StiffnessTensors = StiffnessTensors(),
+):
+    """Calculate elastic tensors as the Voigt averages of a collection of `mineral`s.
+
+    - `minerals` — mineral phases storing orientations and fractional volumes of grains
+    - `phase_assemblage` — collection of unique mineral phases in the aggregate
+    - `phase_fractions` — collection of volume fractions for each phase in
+      `phase_assemblage` (values should sum to 1).
+
+    Raises a ValueError if the minerals contain an unequal number of grains or stored
+    texture results.
+
+    """
+    n_grains = minerals[0].n_grains
+    if not np.all([m.n_grains == n_grains for m in minerals[1:]]):
+        raise ValueError("cannot average minerals with unequal grain counts")
+    n_steps = len(minerals[0].orientations)
+    if not np.all([len(m.orientations) == n_steps for m in minerals[1:]]):
+        raise ValueError(
+            "cannot average minerals with variable-length orientation arrays"
+        )
+    if not np.all([len(m.fractions) == n_steps for m in minerals]):
+        raise ValueError(
+            "cannot average minerals with variable-length grain volume arrays"
+        )
+
+    # TODO: Perform rotation directly on the 6x6 matrices, see Carcione 2007.
+    # This trick is implemented in cpo_elastic_tensor.cc in Aspect.
+    average_tensors = np.zeros((n_steps, 6, 6))
+    for i in range(n_steps):
+        for mineral in minerals:
+            for n in range(n_grains):
+                match mineral.phase:
+                    case _core.MineralPhase.olivine:
+                        average_tensors[i] += _tensors.elastic_tensor_to_voigt(
+                            _tensors.rotate(
+                                elastic_tensors.olivine,
+                                mineral.orientations[i][n, ...].transpose(),
+                            )
+                            * mineral.fractions[i][n]
+                            * phase_fractions[phase_assemblage.index(mineral.phase)]
+                        )
+                    case _core.MineralPhase.enstatite:
+                        average_tensors[i] += _tensors.elastic_tensor_to_voigt(
+                            _tensors.rotate(
+                                elastic_tensors.enstatite,
+                                mineral.orientations[i][n, ...].transpose(),
+                            )
+                            * mineral.fractions[i][n]
+                            * phase_fractions[phase_assemblage.index(mineral.phase)]
+                        )
+                    case _:
+                        raise ValueError(f"unsupported mineral phase: {mineral.phase}")
+    return average_tensors
+
+
+def peridotite_solidus(pressure, fit="Hirschmann2000"):
+    """Get peridotite solidus (i.e. melting) temperature based on experimental fits.
+
+    Pressure is expected to be in GPa.
+
+    Supported fits:
+    - ["Hirschmann2000"](https://doi.org/10.1029/2000GC000070)
+    - ["Herzberg2000"](https://doi.org/10.1029/2000GC000089)
+
+    """
+    match fit:
+        case "Herzberg2000":
+            return 1086 - 5.7 * pressure + 390 * np.log(pressure)
+        case "Hirschmann2000":
+            return 5.104 * pressure**2 + 132.899 * pressure + 1120.661
+        case _:
+            raise ValueError("unsupported fit")

tests/test_doctests.py

@@ -16,7 +16,6 @@ def _get_submodule_list():
     # Reset NumPy print options because doctests are just string matches, and typing out
     # so many significant digits in doctests is annoying.
     np.set_printoptions()
-    np.set_string_function(None)
     modules = ["pydrex." + m.name for m in pkgutil.iter_modules(pydrex.__path__)]
     for module in modules:
         try: