From 631362bd485940e10dd200077e38906dd4fce51b Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 17 Jun 2024 15:40:26 +1000 Subject: [PATCH 1/7] feat: Add single-call API for multiphase texture updates Closes #198. Also rearranges definitions so that globals are above classes are above functions. This allows for mypy type hints of the multiphase function args. --- src/pydrex/minerals.py | 247 ++++++++++++++++++++++++----------------- 1 file changed, 143 insertions(+), 104 deletions(-) diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index 8f7ed4a3..9e3aa87f 100644 --- a/src/pydrex/minerals.py +++ b/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,121 @@ 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): + if i == len(minerals) - 1: # Collect new deformation gradient at last. + deformation_gradient = mineral.update_orientations( + params=params, + deformation_gradient=deformation_gradient, + get_velocity_gradient=get_velocity_gradient, + pathline=pathline, + get_regime=get_regime, + **kwargs, + ) + else: + mineral.update_orientations( + params=params, + deformation_gradient=deformation_gradient, + get_velocity_gradient=get_velocity_gradient, + pathline=pathline, + get_regime=get_regime, + **kwargs, + ) + return 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], +): + """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" + ) + + 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 + + +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") From f13779b553095348c791df54750ab98bf746f9ee Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 17 Jun 2024 15:42:42 +1000 Subject: [PATCH 2/7] fix: Update voigt averaging to use new StiffnessTensors API This fix was brought to you by mypy static type linting. --- src/pydrex/minerals.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index 9e3aa87f..d7c2a9a6 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -723,8 +723,8 @@ def voigt_averages( case _core.MineralPhase.enstatite: average_tensors[i] += _tensors.elastic_tensor_to_voigt( _tensors.rotate( - elastic_tensors["enstatite"], - minerals.orientations[i][n, ...].transpose(), + elastic_tensors.enstatite, + mineral.orientations[i][n, ...].transpose(), ) * mineral.fractions[i][n] * phase_fractions[phase_assemblage.index(mineral.phase)] From e1259e35bd5cb1ddda9b7a568585c0175a46b4ab Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 17 Jun 2024 15:53:36 +1000 Subject: [PATCH 3/7] fix: Update Numpy printoptions settings for Numpy 2.0 --- src/pydrex/logger.py | 8 +++++--- tests/test_doctests.py | 1 - 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/pydrex/logger.py b/src/pydrex/logger.py index 8c5b455f..08a4bf4e 100644 --- a/src/pydrex/logger.py +++ b/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 diff --git a/tests/test_doctests.py b/tests/test_doctests.py index ffb037b6..21c0c45c 100644 --- a/tests/test_doctests.py +++ b/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: From e731cdbf74ff748ed3c4d5be43aca663cca246c5 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 17 Jun 2024 15:54:33 +1000 Subject: [PATCH 4/7] fix: Export in top-level namespace --- src/pydrex/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pydrex/__init__.py b/src/pydrex/__init__.py index 18faa0dd..ade8550f 100644 --- a/src/pydrex/__init__.py +++ b/src/pydrex/__init__.py @@ -154,6 +154,7 @@ StiffnessTensors, peridotite_solidus, voigt_averages, + update_all, ) from pydrex.stats import ( misorientation_hist, From 51d7bf3080d99e2b8194225089c8b542da15a29f Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 17 Jun 2024 16:19:14 +1000 Subject: [PATCH 5/7] feat: Allow passing custom instance of StiffnessTensors to Voigt averaging --- src/pydrex/minerals.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index d7c2a9a6..4fbccb79 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -677,6 +677,7 @@ 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. @@ -702,8 +703,6 @@ def voigt_averages( "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)) From 618b5ce71db58354c2a2e2099518a2643b114b16 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 17 Jun 2024 16:19:31 +1000 Subject: [PATCH 6/7] docs: Update changelog with changes since 0.0.1 tag --- CHANGELOG.md | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4d54be74..8e11619f 100644 --- a/CHANGELOG.md +++ b/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. From 6a66c2b4423aec2e0e8a5af40884d93fca16cede Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jun 2024 10:31:26 +1000 Subject: [PATCH 7/7] fix: Address review comment --- src/pydrex/minerals.py | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index 4fbccb79..c74ab89e 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -651,25 +651,16 @@ def update_all( """ for i, mineral in enumerate(minerals): - if i == len(minerals) - 1: # Collect new deformation gradient at last. - deformation_gradient = mineral.update_orientations( - params=params, - deformation_gradient=deformation_gradient, - get_velocity_gradient=get_velocity_gradient, - pathline=pathline, - get_regime=get_regime, - **kwargs, - ) - else: - mineral.update_orientations( - params=params, - deformation_gradient=deformation_gradient, - get_velocity_gradient=get_velocity_gradient, - pathline=pathline, - get_regime=get_regime, - **kwargs, - ) - return deformation_gradient + # 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).