pydrex.core

  1"""> PyDRex: Core D-Rex functions and enums.
  2
  3The function `derivatives` implements the core D-Rex solver, which computes the
  4crystallographic rotation rate and changes in fractional grain volumes.
  5
  6**Acronyms:**
  7- CRSS = Critical Resolved Shear Stress,
  8    i.e. threshold stress required to initiate slip on a slip system,
  9    normalised to the stress required to initiate slip on the softest slip system
 10
 11"""
 12
 13from dataclasses import asdict, dataclass
 14from enum import IntEnum, unique
 15
 16import numba as nb
 17import numpy as np
 18
 19# NOTE: Do NOT import any pydrex submodules here to avoid cyclical imports.
 20
 21
 22_USE_ORIGINAL_DREX = False  # Switch to use D-Rex 2004 version without corrections.
 23
 24
 25PERMUTATION_SYMBOL = np.array(
 26    [
 27        [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, -1.0, 0.0]],
 28        [[0.0, 0.0, -1.0], [0.0, 0.0, 0.0], [1.0, 0.0, 0.0]],
 29        [[0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
 30    ]
 31)
 32"""Sometimes called the Levi-Civita symbol."""
 33
 34
 35@unique
 36class MineralPhase(IntEnum):
 37    """Supported mineral phases.
 38
 39    Forsterite and fayalite are grouped into “olivine”, because we treat them as
 40    rheologically equivalent.
 41
 42    """
 43
 44    olivine = 0
 45    """(Mg,Fe)₂SiO₄"""
 46    enstatite = 1
 47    """MgSiO₃"""
 48
 49
 50@unique
 51class MineralFabric(IntEnum):
 52    """Supported mineral fabrics.
 53
 54    The following fabric types are supported:
 55    - olivine A-E type fabrics according to e.g.
 56      [Karato et al. (2008)](https://doi.org/10.1146%2Fannurev.earth.36.031207.124120).
 57    - enstatite AB fabric, see
 58      [Bernard et al. (2021)](https://doi.org/10.1016/j.tecto.2021.228954).
 59
 60    """
 61
 62    olivine_A = 0
 63    olivine_B = 1
 64    olivine_C = 2
 65    olivine_D = 3
 66    olivine_E = 4
 67    enstatite_AB = 5
 68
 69
 70@unique
 71class DeformationRegime(IntEnum):
 72    r"""Ordinals to track distinct regimes of dominant deformation mechanisms.
 73
 74    The mechanism of deformation that dominates accommodation of plastic deformation
 75    depends in general both on material properties such as grain size and mineral phase
 76    content as well as on thermodynamic properties such as temperature, pressure and
 77    water fugacity.
 78
 79    The activity of diffusive mechanisms depends more strongly on grain size, whereas
 80    that of dislocation mechanisms depends more strongly on temperature. High
 81    temperatures enable more frequent recovery of dislocation density equilibrium via
 82    e.g. dislocation climb. Dislocation mechanisms are often accompanied by dynamic
 83    recrystallisation, which acts as an additional recovery mechanism.
 84
 85    Diffusive mechanisms are expected, however, to become dominant at depth, as these
 86    mechanisms are primarily facilitated by Si diffusion in olivine, which is
 87    temperature-dependent.
 88
 89    Rheology in the intra-granular dislocation regime was classically described by
 90    separate flow laws depending on temperature; a power-law at high temperature
 91    [$\dot{ε} ∝ σⁿ$] and an exponential-law at low temperature [$\dot{ε} ∝ \exp(σ)$].
 92    More recent work has suggested unified dislocation creep flow laws,
 93    e.g. [Gouriet et al. 2019](http://dx.doi.org/10.1016/j.epsl.2018.10.049),
 94    [Garel et al. 2020](http://dx.doi.org/10.1016/j.epsl.2020.116243) and
 95    [Demouchy et al. 2023](http://dx.doi.org/10.2138/gselements.19.3.151).
 96
 97    .. note:: Although a draft texture evolution behaviour is implemented in the
 98        `frictional_yielding` regime, it is experimental and not yet configurable via
 99        the parameter interface.
100
101    """
102
103    min_viscosity = 0
104    """Arbitrary lower-bound viscosity regime."""
105    matrix_diffusion = 1
106    """Intra-granular Nabarro-Herring creep, i.e. grains diffuse through the matrix."""
107    boundary_diffusion = 2
108    """Inter-granular Coble creep, i.e. grains diffuse along grain boundaries."""
109    sliding_diffusion = 3
110    """Inter-granular diffusion-assisted grain-boundary sliding (diffGBS)."""
111    matrix_dislocation = 4
112    """Intra-granular dislocation creep (glide + climb) and dynamic recrystallisation."""
113    sliding_dislocation = 5
114    """Inter-granular dislocation-assisted grain-boundary sliding (disGBS)."""
115    frictional_yielding = 6
116    """Frictional sliding along micro-fractures (Byerlee's law for yield strength)."""
117    max_viscosity = 7
118    """Arbitrary upper-bound viscosity regime."""
119
120
121@dataclass(frozen=True)
122class DefaultParams:
123    phase_assemblage: tuple = (MineralPhase.olivine,)
124    """Mineral phases present in the aggregate."""
125    phase_fractions: tuple = (1.0,)
126    """Volume fractions of each mineral phase present in the aggregate."""
127    stress_exponent: float = 1.5
128    """The value for $p$ in $ρ ∝ τᵖ$ where $ρ$ is the dislocation density and $τ$ the shear stress.
129
130    Default value taken from [Kaminski et al. 2004](https://doi.org/10.1111/j.1365-246X.2004.02308.x)
131    based on studies referenced therein.
132
133    """
134    deformation_exponent: float = 3.5
135    """The value for $n$ in $τ ∝ |D|^{1/n}$ where $τ$ is the shear stress and D the deformation rate.
136
137    Default value taken from [Kaminski et al. 2004](https://doi.org/10.1111/j.1365-246X.2004.02308.x)
138    based on studies referenced therein.
139
140    """
141    gbm_mobility: int = 125
142    """Dimensionless grain boundary mobility parameter (M*).
143
144    This controls the rate of all dynamic recrystallisation processes in
145    `DeformationRegime.matrix_dislocation`.
146
147    .. note:: This parameter is not easy to constrain. Reasonable values may lie in the
148        range [10, 150]. Comparing outputs for multiple M* values is recommended.
149
150    """
151    gbs_threshold: float = 0.3
152    """Grain boundary sliding threshold.
153
154    In `DeformationRegime.matrix_dislocation` or `DeformationRegime.sliding_dislocation`
155    this controls the smallest size of a grain, relative to its original size,
156    at which it will still deform via dislocation creep.
157    Smaller grains will instead deform via grain boundary sliding,
158    therefore not contributing to any further texture evolution.
159
160    .. note:: Values for this parameter do NOT correspond to a physical grain size
161        threshold. Grains in PyDRex are a surrogate numerical discretisation. There are
162        likely to be feedbacks between this parameter and `number_of_grains`. Values of
163        ~0.3 have been used in numerous studies which employed ~3000 grains.
164
165    """
166    nucleation_efficiency: float = 5.0
167    """Dimensionless nucleation efficiency (λ*).
168
169    This controls the nucleation of subgrains in `DeformationRegime.matrix_dislocation`.
170
171    The default value comes from [Kaminski & Ribe 2001](https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9).
172
173    """
174    number_of_grains: int = 3500
175    """Number of surrogate grains for numerical discretisation of the aggregate.
176
177    This is a numerical discretisation, so generally more grains is better. However,
178    there is a sharp tradeoff in performance, especially for M-index calculation.
179
180    """
181    initial_olivine_fabric: MineralFabric = MineralFabric.olivine_A
182    """Olivine fabric (CRSS distribution) at the beginning of the simulation.
183
184    In general, the fabric should be allowed to change during the simulation,
185    but that is not yet implemented.
186
187    """
188    disl_Peierls_stress: float = 2
189    """Stress barrier in GPa for activation of dislocation motion at low temperatures.
190
191    - 2GPa suggested by [Demouchy et al. 2023](http://dx.doi.org/10.2138/gselements.19.3.151)
192
193    .. note:: Not relevant if the unified dislocation creep flow law is used.
194
195    """
196    # NOTE: For now, we just roll this into the prefactors.
197    # water_fugacity: float = 1.5
198    # """Approximate constant water fugacity of the aggregate (GPa)."""
199    # TODO: Check units of Garel 2020 pre-exp value below.
200    # TODO: Find and add references for the other value.
201    disl_prefactors: tuple = (1e-16, 1e-17)
202    """Prefactors for dislocation creep exponential & power laws (s⁻¹).
203
204    - B = 1.7e-16 s⁻¹ for the exponential law suggested by [Demouchy et al. 2023](http://dx.doi.org/10.2138/gselements.19.3.151)
205    - pre-exponential factor of 4.4e-17 Pa⁻ⁿs⁻¹ used for the exponential law by [Garel et al. 2020](http://dx.doi.org/10.1016/j.epsl.2020.116243)
206
207    """
208    # TODO: Add references, tweak default values if necessary.
209    diff_prefactors: tuple = (1e-10, 1e-10)
210    r"""Prefactors for (matrix and boundary) diffusion creep power laws (s⁻¹).
211
212    Dependence on molar volume and physical grain size are suppressed because these
213    diagnostics are not readily available. The prefactor is roughly equal to
214    $(Vₘ D⁰)/d²$ for $Vₘ$ the molar volume, $D⁰$ the reference diffusion coefficient and
215    $d²$ an average grain size assumed to be constant.
216
217    """
218    disl_lowtemp_switch: float = 0.7
219    """Threshold homologous temperature below which to use the exponential flow law.
220
221    The default value suggested here comes from
222    [Demouchy et al. 2023](http://dx.doi.org/10.2138/gselements.19.3.151).
223
224    .. note:: Not relevant if the unified dislocation creep flow law is used instead.
225
226    """
227    # TODO: Add more references from experimental studies/reviews.
228    disl_activation_energy: float = 460.0
229    """Activation energy for dislocation creep power law (kJ/mol).
230
231    - 443 kJ/mol used by [Gouriet et al. 2019](http://dx.doi.org/10.1016/j.epsl.2018.10.049),
232        but for an exponential law
233    - 540 kJ/mol used by [Garel et al. 2020](http://dx.doi.org/10.1016/j.epsl.2020.116243)
234
235    """
236    # TODO: Add more references, tweak default if necessary.
237    disl_activation_volume: float = 12.0
238    """Activation volume for dislocation creep power law (cm³/mol).
239
240    - 0, 6 and 12 cm³/mol used by [Garel et al. 2020](http://dx.doi.org/10.1016/j.epsl.2020.116243)
241
242    """
243    # TODO: Add more references, justify choice of lower value than Garel 2020.
244    diff_activation_energies: tuple = (430.0, 330)
245    """Activation energies for (matrix and boundary) diffusion creep power laws (kJ/mol).
246
247    - 530 kJ/mol reported for Si self-diffusion in olivine¹
248    - 410 kJ/mol used by [Garel et al. 2020](http://dx.doi.org/10.1016/j.epsl.2020.116243)
249
250    ¹[Dohmen et al. 2002](http://dx.doi.org/10.1029/2002GL015480)
251
252    """
253    # TODO: Add more references, tweak default values if necessary.
254    diff_activation_volumes: tuple = (4.0, 4.0)
255    """Activation volumes for (matrix and boundary) diffusion creep power laws (cm³/mol).
256
257    - 4 kJ/mol used by [Garel et al. 2020](http://dx.doi.org/10.1016/j.epsl.2020.116243)
258
259    """
260    disl_coefficients: tuple = (
261        4.4e8,  # a₀
262        -5.26e4,  # b₀
263        2.11e-2,  # a₁
264        1.74e-4,  # b₁
265        -41.8,  # a₂
266        4.21e-2,  # b₂
267        -1.14e-5,  # c₂
268    )
269    """Coefficients for polynomials used in the alternative dislocation creep flow law.
270
271    The defaults are for the TANH variant of the unified creep law, proposed in
272    [Garel et al. 2020](http://dx.doi.org/10.1016/j.epsl.2020.116243).
273    By contrast, [Gouriet et al. 2019](http://dx.doi.org/10.1016/j.epsl.2018.10.049)
274    used the ERF variant with: [4.4e8, -2.2e4, 3e-2, 1.3e-4, -42, 4.2e-2, -1.1e-5].
275
276    """
277
278    def as_dict(self):
279        """Return mutable copy of default arguments as a dictionary."""
280        return asdict(self)
281
282
283@nb.njit
284def get_crss(phase: MineralPhase, fabric: MineralFabric) -> np.ndarray:
285    """Get Critical Resolved Shear Stress for the mineral `phase` and `fabric`.
286
287    Returns an array of the normalised threshold stresses required to activate slip on
288    each slip system. Olivine slip systems are ordered according to the convention used
289    for `pydrex.minerals.OLIVINE_SLIP_SYSTEMS`.
290
291    """
292    if phase == MineralPhase.olivine:
293        match fabric:
294            case MineralFabric.olivine_A:
295                return np.array([1, 2, 3, np.inf])
296            case MineralFabric.olivine_B:
297                return np.array([3, 2, 1, np.inf])
298            case MineralFabric.olivine_C:
299                return np.array([3, 2, np.inf, 1])
300            case MineralFabric.olivine_D:
301                return np.array([1, 1, 3, np.inf])
302            case MineralFabric.olivine_E:
303                return np.array([3, 1, 2, np.inf])
304            case _:
305                raise ValueError(f"unsupported olivine fabric: {fabric}")
306    elif phase == MineralPhase.enstatite:
307        if fabric == MineralFabric.enstatite_AB:
308            return np.array([np.inf, np.inf, np.inf, 1])
309        raise ValueError(f"unsupported enstatite fabric: {fabric}")
310    raise ValueError(f"phase must be a valid `MineralPhase`, not {phase}")
311
312
313# 12+ args is a lot, but this way we can use numba
314# (only primitives and numpy containers allowed).
315@nb.njit(fastmath=True)
316def derivatives(
317    regime: DeformationRegime,
318    phase: MineralPhase,
319    fabric: MineralFabric,
320    n_grains: int,
321    orientations: np.ndarray,
322    fractions: np.ndarray,
323    strain_rate: np.ndarray,
324    velocity_gradient: np.ndarray,
325    deformation_gradient_spin: np.ndarray,
326    stress_exponent: float,
327    deformation_exponent: float,
328    nucleation_efficiency: float,
329    gbm_mobility: float,
330    volume_fraction: float,
331):
332    """Get derivatives of orientation and volume distribution.
333
334    - `regime` — ordinal number of the local deformation mechanism
335    - `phase` — ordinal number of the mineral phase
336    - `fabric` — ordinal number of the fabric type
337    - `n_grains` — number of "grains" i.e. discrete volume segments
338    - `orientations` — `n_grains`x3x3 orientations (direction cosines)
339    - `fractions` — volume fractions of the grains relative to aggregate volume
340    - `strain_rate` — 3x3 dimensionless macroscopic strain-rate tensor
341    - `velocity_gradient` — 3x3 dimensionless macroscopic velocity gradient
342    - `deformation_gradient_spin` — 3x3 spin tensor defining the rate of rotation of the
343      finite strain ellipse
344    - `stress_exponent` — `p` in `dislocation_density ∝ shear_stress^p`
345    - `deformation_exponent` — `n` in `shear_stress ∝ |deformation_rate|^(1/n)`
346    - `nucleation_efficiency` — parameter controlling grain nucleation
347    - `gmb_mobility` — grain boundary mobility parameter
348    - `volume_fraction` — volume fraction of the mineral phase relative to other phases
349      (multiphase simulations)
350
351    Returns a tuple with the rotation rates and grain volume fraction changes.
352
353    """
354    if regime == DeformationRegime.min_viscosity:
355        # Do absolutely nothing, all derivatives are zero.
356        return (
357            np.repeat(np.eye(3), n_grains).reshape(3, 3, n_grains).transpose(),
358            np.zeros(n_grains),
359        )
360    elif regime == DeformationRegime.matrix_diffusion:
361        # Passive rotation based on macroscopic vorticity for diffusion creep?
362        # vorticity = 0.5 * (velocity_gradient - velocity_gradient.transpose())
363        # Passive rotation based on spin of F for diffusion creep.
364        vorticity = deformation_gradient_spin
365        # Or just don't change at all?
366        # vorticity = np.zeros((3, 3))
367        # This 💃 is because numba doesn't let us use np.tile or even np.array([a] * n).
368        return (
369            np.repeat(vorticity.transpose(), n_grains)
370            .reshape(3, 3, n_grains)
371            .transpose(),
372            np.zeros(n_grains),
373        )
374    elif regime == DeformationRegime.boundary_diffusion:
375        raise ValueError("this deformation mechanism is not yet supported.")
376    elif regime == DeformationRegime.sliding_diffusion:
377        # Miyazaki et al. 2013 wrote controversial Nature article proposing that CPO can
378        # develop in the diffGBS regime. However, Karato 2024 gives some convincing
379        # arguments in the Journal of Geodynamics for why their results are likely
380        # inapplicable to Earth's upper mantle.
381        raise ValueError("this deformation mechanism is not yet supported.")
382    elif regime == DeformationRegime.matrix_dislocation:
383        # Based on subroutine DERIV in original Fortran.
384        strain_energies = np.empty(n_grains)
385        orientations_diff = np.empty((n_grains, 3, 3))
386        for grain_index in range(n_grains):
387            orientation_change, strain_energy = _get_rotation_and_strain(
388                phase,
389                fabric,
390                orientations[grain_index],
391                strain_rate,
392                velocity_gradient,
393                stress_exponent,
394                deformation_exponent,
395                nucleation_efficiency,
396            )
397            orientations_diff[grain_index] = orientation_change
398            strain_energies[grain_index] = strain_energy
399        # Volume average mean strain energy.
400        mean_energy = np.sum(fractions * strain_energies)
401        # Strain energy residual.
402        strain_residuals = mean_energy - strain_energies
403        fractions_diff = volume_fraction * gbm_mobility * fractions * strain_residuals
404        return orientations_diff, fractions_diff
405    elif regime == DeformationRegime.sliding_dislocation:
406        raise ValueError("this deformation mechanism is not yet supported.")
407    elif regime == DeformationRegime.frictional_yielding:
408        # For now, the same as matrix_dislocation, but we smooth the strain energy
409        # distribution and the orientation changes, since some energy is lost to
410        # micro-fracturing. Also increase the GBS threshold, dislocations will tend to
411        # pile up at grain boundaries.
412        # TODO: Maybe modify the stress/deformation exponents?
413        # TODO: Reduce nucleation efficiency?
414        strain_energies = np.empty(n_grains)
415        orientations_diff = np.empty((n_grains, 3, 3))
416        for grain_index in range(n_grains):
417            orientation_change, strain_energy = _get_rotation_and_strain(
418                phase,
419                fabric,
420                orientations[grain_index],
421                strain_rate,
422                velocity_gradient,
423                stress_exponent,
424                deformation_exponent,
425                nucleation_efficiency,
426            )
427            orientations_diff[grain_index] = 0.3 * orientation_change
428            strain_energies[grain_index] = strain_energy
429        # Volume average mean strain energy.
430        mean_energy = np.sum(fractions * strain_energies)
431        # Strain energy residuals, minus the energy lost to micro-fracturing.
432        strain_residuals = 0.3 * (mean_energy - strain_energies)
433        fractions_diff = volume_fraction * gbm_mobility * fractions * strain_residuals
434        return orientations_diff, fractions_diff
435    elif regime == DeformationRegime.max_viscosity:
436        # Do absolutely nothing, all derivatives are zero.
437        return (
438            np.repeat(np.eye(3), n_grains).reshape(3, 3, n_grains).transpose(),
439            np.zeros(n_grains),
440        )
441    else:
442        raise ValueError(f"regime must be a valid `DeformationRegime`, not {regime}")
443
444
445@nb.njit(fastmath=True)
446def _get_deformation_rate(
447    phase: MineralPhase, orientation: np.ndarray, slip_rates: np.ndarray
448):
449    """Calculate deformation rate tensor for olivine or enstatite.
450
451    Calculate the deformation rate with respect to the local coordinate frame,
452    defined by the principal strain axes (finite strain ellipsoid).
453
454    - `phase` — ordinal number of the mineral phase
455    - `orientation` — 3x3 orientation matrix (direction cosines)
456    - `slip_rates` — slip rates relative to slip rate on softest slip system
457
458    """
459    # This is called the 'G'/'slip' tensor in the Fortran code, aka the Schmid tensor.
460    deformation_rate = np.empty((3, 3))
461    for i in range(3):
462        for j in range(3):
463            deformation_rate[i, j] = 2 * (
464                slip_rates[0] * orientation[0, i] * orientation[1, j]
465                + slip_rates[1] * orientation[0, i] * orientation[2, j]
466                + slip_rates[2] * orientation[2, i] * orientation[1, j]
467                + slip_rates[3] * orientation[2, i] * orientation[0, j]
468            )
469    return deformation_rate
470
471
472@nb.njit(fastmath=True)
473def _get_slip_rate_softest(deformation_rate: np.ndarray, velocity_gradient: np.ndarray):
474    """Calculate dimensionless strain rate on the softest slip system.
475
476    - `deformation_rate` — 3x3 dimensionless deformation rate matrix
477    - `velocity_gradient` — 3x3 dimensionless velocity gradient matrix
478
479    """
480    # See eq. 4 in Fraters 2021.
481    enumerator = 0
482    denominator = 0
483
484    for j in range(3):
485        # NOTE: Mistake in original DRex code (j + 2), see Fraters & Billen 2021 S1.
486        if _USE_ORIGINAL_DREX:
487            k = (j + 2) % 3
488        else:
489            k = (j + 1) % 3
490        enumerator -= (velocity_gradient[j, k] - velocity_gradient[k, j]) * (
491            deformation_rate[j, k] - deformation_rate[k, j]
492        )
493        # NOTE: Mistake in Kaminski 2001 eq. 7: kl+1 instead of kk+1
494        # See Fraters & Billen 2021 supplementary informaton S1.
495        denominator -= (deformation_rate[j, k] - deformation_rate[k, j]) ** 2
496
497        for L in range(3):
498            enumerator += 2 * deformation_rate[j, L] * velocity_gradient[j, L]
499            denominator += 2 * deformation_rate[j, L] ** 2
500
501    # Avoid zero divisions:
502    # Tiny denominator means that relevant deformation_rate entries are zero.
503    # No deformation rate means no slip rate.
504    if -1e-15 < denominator < 1e-15:
505        return 0.0
506    return enumerator / denominator
507
508
509@nb.njit(fastmath=True)
510def _get_slip_rates_olivine(
511    invariants: np.ndarray,
512    slip_indices: np.ndarray,
513    crss: np.ndarray,
514    deformation_exponent: float,
515):
516    """Calculate relative slip rates of the active slip systems for olivine.
517
518    - `invariants` — strain rate invariants for the four slip systems
519    - `slip_indices` — indices that sort the CRSS by increasing slip-rate activity
520    - `crss` — reference resolved shear stresses (CRSS), see `pydrex.fabric`
521    - `deformation_exponent` — `n` in `shear_stress ∝ |deformation_rate|^(1/n)`
522
523    """
524    i_inac, i_min, i_int, i_max = slip_indices
525    # Ratio of slip rates on each slip system to slip rate on softest slip system.
526    # Softest slip system has max. slip rate (aka activity).
527    # See eq. 5, Kaminski 2001.
528    prefactor = crss[i_max] / invariants[i_max]
529    ratio_min = prefactor * invariants[i_min] / crss[i_min]
530    ratio_int = prefactor * invariants[i_int] / crss[i_int]
531    slip_rates = np.empty(4)
532    slip_rates[i_inac] = 0  # Hardest system is completely inactive in olivine.
533    slip_rates[i_min] = ratio_min * np.abs(ratio_min) ** (deformation_exponent - 1)
534    slip_rates[i_int] = ratio_int * np.abs(ratio_int) ** (deformation_exponent - 1)
535    slip_rates[i_max] = 1
536    return slip_rates
537
538
539@nb.njit(fastmath=True)
540def _get_slip_invariants(strain_rate: np.ndarray, orientation: np.ndarray):
541    r"""Calculate strain rate invariants for minerals with four slip systems.
542
543    Calculates $I = ∑_{ij} l_{i} n_{j} \dot{ε}_{ij}$ for each slip sytem of:
544    - (010)[100]
545    - (001)[100]
546    - (010)[001]
547    - (100)[001]
548    Only the last return value is relevant for enstatite.
549    These are not configurable for now.
550
551    - `strain_rate` — 3x3 dimensionless strain rate matrix
552    - `orientation` — 3x3 orientation matrix (direction cosines)
553
554    """
555    invariants = np.zeros(4)
556    for i in range(3):
557        for j in range(3):
558            # (010)[100]
559            invariants[0] += strain_rate[i, j] * orientation[0, i] * orientation[1, j]
560            # (001)[100]
561            invariants[1] += strain_rate[i, j] * orientation[0, i] * orientation[2, j]
562            # (010)[001]
563            invariants[2] += strain_rate[i, j] * orientation[2, i] * orientation[1, j]
564            # (100)[001]
565            invariants[3] += strain_rate[i, j] * orientation[2, i] * orientation[0, j]
566    return invariants
567
568
569@nb.njit(fastmath=True)
570def _get_orientation_change(
571    orientation: np.ndarray,
572    velocity_gradient: np.ndarray,
573    deformation_rate: np.ndarray,
574    slip_rate_softest: float,
575):
576    """Calculate the rotation rate for a grain undergoing dislocation creep.
577
578    - `orientation` — 3x3 orientation matrix (direction cosines)
579    - `velocity_gradient` — 3x3 dimensionless velocity gradient matrix
580    - `deformation_rate` — 3x3 dimensionless deformation rate matrix
581    - `slip_rate_softest` — slip rate on the softest (most active) slip system
582
583    """
584    orientation_change = np.zeros((3, 3))
585    # Spin vector for the grain, see eq. 3 in Fraters 2021.
586    # Includes rigid body rotation as well as the plastic deformation contribution.
587    # This means that even with no active slip systems, the rotation will be nonzero
588    # if there is a rotational (i.e. antisymmetric) component of the velocity_gradient.
589    spin_vector = np.empty(3)
590    for j in range(3):
591        r = (j + 1) % 3
592        s = (j + 2) % 3
593        spin_vector[j] = (
594            (velocity_gradient[s, r] - velocity_gradient[r, s])
595            - (deformation_rate[s, r] - deformation_rate[r, s]) * slip_rate_softest
596        ) / 2
597
598    # Calculate rotation rate, see eq. 9 Kaminski & Ribe (2001).
599    # Equivalent to:
600    #   spin_matrix = np.einsum("ikj,k->ij", PERMUTATION_SYMBOL, spin_vector)
601    #   orientation_change = spin_matrix.transpose() @ orientation
602    for p in range(3):
603        for q in range(3):
604            for r in range(3):
605                for s in range(3):
606                    orientation_change[p, q] += (
607                        PERMUTATION_SYMBOL[q, r, s] * orientation[p, s] * spin_vector[r]
608                    )
609
610    return orientation_change
611
612
613@nb.njit(fastmath=True)
614def _get_strain_energy(
615    crss: np.ndarray,
616    slip_rates: np.ndarray,
617    slip_indices: np.ndarray,
618    slip_rate_softest: float,
619    stress_exponent: float,
620    deformation_exponent: float,
621    nucleation_efficiency: float,
622):
623    """Calculate strain energy due to dislocations for an olivine grain.
624
625    - `crss` — reference resolved shear stresses (CRSS), see `pydrex.fabric`
626    - `slip_rates` — slip rates relative to slip rate on softest slip system
627    - `slip_indices` — indices that sort the CRSS by increasing slip-rate activity
628    - `slip_rate_softest` — slip rate on the softest (most active) slip system
629    - `stress_exponent` — `p` in `dislocation_density ∝ shear_stress^p`
630    - `deformation_exponent` — `n` in `shear_stress ∝ |deformation_rate|^(1/n)`
631    - `nucleation_efficiency` — parameter controlling grain nucleation
632
633    Note that "new" grains are assumed to rotate with their parent.
634
635    """
636    strain_energy = 0.0
637    # Dimensionless dislocation density for each slip system.
638    # See eq. 16 Fraters 2021.
639    # NOTE: Mistake in eq. 11, Kaminski 2004: spurious division by strain rate scale.
640    # NOTE: Here we call 'p' the 'stress_exponent' and 'n' the 'deformation_exponent',
641    # but in the original they use the variable 'stress_exponent' for 'n' (3.5).
642    for i in range(3):
643        dislocation_density = (1 / crss[i]) ** (
644            deformation_exponent - stress_exponent
645        ) * np.abs(slip_rates[i] * slip_rate_softest) ** (
646            stress_exponent / deformation_exponent
647        )
648        # Dimensionless strain energy for this grain, see eq. 14, Fraters 2021.
649        strain_energy += dislocation_density * np.exp(
650            -nucleation_efficiency * dislocation_density**2
651        )
652    return strain_energy
653
654
655@nb.njit(fastmath=True)
656def _get_rotation_and_strain(
657    phase: MineralPhase,
658    fabric: MineralFabric,
659    orientation: np.ndarray,
660    strain_rate: np.ndarray,
661    velocity_gradient: np.ndarray,
662    stress_exponent: float,
663    deformation_exponent: float,
664    nucleation_efficiency: float,
665):
666    """Get the crystal axes rotation rate and strain energy of individual grain.
667
668    - `phase` — ordinal number of the mineral phase
669    - `fabric` — ordinal number of the fabric type
670    - `orientation` — 3x3 orientation matrix (direction cosines)
671    - `strain_rate` — 3x3 dimensionless strain rate matrix
672    - `velocity_gradient` — 3x3 dimensionless velocity gradient matrix
673    - `stress_exponent` — `p` in `dislocation_density ∝ shear_stress^p`
674    - `deformation_exponent` — `n` in `shear_stress ∝ |deformation_rate|^(1/n)`
675    - `nucleation_efficiency — parameter controlling grain nucleation
676
677    Note that "new" grains are assumed to rotate with their parent.
678
679    Returns a tuple with the rotation rate of the crystalline axes
680    with respect to the principal strain axes and strain energy of the grain.
681
682    """
683    crss = get_crss(phase, fabric)
684    slip_invariants = _get_slip_invariants(strain_rate, orientation)
685    if phase == MineralPhase.olivine:
686        slip_indices = np.argsort(np.abs(slip_invariants / crss))
687        slip_rates = _get_slip_rates_olivine(
688            slip_invariants,
689            slip_indices,
690            crss,
691            deformation_exponent,
692        )
693    elif phase == MineralPhase.enstatite:
694        slip_indices = np.argsort(1 / crss)
695        slip_rates = np.zeros(4)
696        if _USE_ORIGINAL_DREX:
697            slip_rates[-1] = 1  # Original had an arbitrary slip always active (L1410).
698        else:
699            # Assumes exclusively (100)[001] slip for enstatite.
700            if np.abs(slip_invariants[-1]) > 1e-15:
701                slip_rates[-1] = 1
702    else:
703        assert False  # Should never happen.
704
705    deformation_rate = _get_deformation_rate(phase, orientation, slip_rates)
706    slip_rate_softest = _get_slip_rate_softest(deformation_rate, velocity_gradient)
707    orientation_change = _get_orientation_change(
708        orientation,
709        velocity_gradient,
710        deformation_rate,
711        slip_rate_softest,
712    )
713    strain_energy = _get_strain_energy(
714        crss,
715        slip_rates,
716        slip_indices,
717        slip_rate_softest,
718        stress_exponent,
719        deformation_exponent,
720        nucleation_efficiency,
721    )
722    return orientation_change, strain_energy