tests.test_core

  1"""> PyDRex: Tests for core D-Rex routines."""
  2
  3import contextlib as cl
  4
  5import numpy as np
  6import pytest
  7from numpy import testing as nt
  8from pydrex import core as _core
  9from pydrex import io as _io
 10from pydrex import logger as _log
 11from pydrex import minerals as _minerals
 12from pydrex import utils as _utils
 13from pydrex import visualisation as _vis
 14from scipy.spatial.transform import Rotation
 15
 16# Subdirectory of `outdir` used to store outputs from these tests.
 17SUBDIR = "core"
 18
 19
 20class TestDislocationCreepOPX:
 21    """Single-grain orthopyroxene crystallographic rotation rate tests."""
 22
 23    class_id = "dislocation_creep_OPX"
 24
 25    @pytest.mark.skipif(_utils.in_ci("win32"), reason="Not equal to tolerance")
 26    def test_shear_dudz(self, outdir):
 27        test_id = "dudz"
 28        optional_logging = cl.nullcontext()
 29        if outdir is not None:
 30            optional_logging = _io.logfile_enable(
 31                f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.log"
 32            )
 33        with optional_logging:
 34            for θ in np.mgrid[0 : 2 * np.pi : 360j]:
 35                _log.debug("θ (°): %s", np.rad2deg(θ))
 36                orientations_diff, fractions_diff = _core.derivatives(
 37                    regime=_core.DeformationRegime.matrix_dislocation,
 38                    phase=_core.MineralPhase.enstatite,
 39                    fabric=_core.MineralFabric.enstatite_AB,
 40                    n_grains=1,
 41                    orientations=Rotation.from_rotvec([[0, θ, 0]]).as_matrix(),
 42                    fractions=np.array([1.0]),
 43                    strain_rate=np.array([[0, 0, 1], [0, 0, 0], [1, 0, 0]]),
 44                    velocity_gradient=np.array([[0, 0, 2], [0, 0, 0], [0, 0, 0]]),
 45                    deformation_gradient_spin=np.full((3, 3), np.nan),
 46                    stress_exponent=1.5,
 47                    deformation_exponent=3.5,
 48                    nucleation_efficiency=5,
 49                    gbm_mobility=0,
 50                    volume_fraction=1.0,
 51                )
 52                cosθ = np.cos(θ)
 53                cos2θ = np.cos(2 * θ)
 54                sinθ = np.sin(θ)
 55                target_orientations_diff = (1 + cos2θ) * np.array(
 56                    [[sinθ, 0, -cosθ], [0, 0, 0], [cosθ, 0, sinθ]]
 57                )
 58                np.testing.assert_allclose(
 59                    orientations_diff[0], target_orientations_diff
 60                )
 61                assert np.isclose(np.sum(fractions_diff), 0.0)
 62
 63    def test_shear_dvdx(self, outdir):
 64        test_id = "dvdx"
 65        optional_logging = cl.nullcontext()
 66        if outdir is not None:
 67            optional_logging = _io.logfile_enable(
 68                f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.log"
 69            )
 70        with optional_logging:
 71            for θ in np.linspace(0, 2 * np.pi, 361):
 72                _log.debug("θ (°): %s", np.rad2deg(θ))
 73                orientations = Rotation.from_rotvec([[0, 0, θ]]).as_matrix()
 74                deformation_rate = _core._get_deformation_rate(
 75                    _core.MineralPhase.enstatite,
 76                    orientations[0],
 77                    np.array([0, 0, 0, 0]),
 78                )
 79                np.testing.assert_allclose(deformation_rate.flatten(), np.zeros(9))
 80                orientations_diff, fractions_diff = _core.derivatives(
 81                    regime=_core.DeformationRegime.matrix_dislocation,
 82                    phase=_core.MineralPhase.enstatite,
 83                    fabric=_core.MineralFabric.enstatite_AB,
 84                    n_grains=1,
 85                    orientations=orientations,
 86                    fractions=np.array([1.0]),
 87                    strain_rate=np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]]),
 88                    velocity_gradient=np.array([[0, 0, 0], [2, 0, 0], [0, 0, 0]]),
 89                    deformation_gradient_spin=np.full((3, 3), np.nan),
 90                    stress_exponent=1.5,
 91                    deformation_exponent=3.5,
 92                    nucleation_efficiency=5,
 93                    gbm_mobility=0,
 94                    volume_fraction=1.0,
 95                )
 96                # Can't activate the (100)[001] slip system, no plastic deformation.
 97                # Only passive (rigid-body) rotation due to the velocity gradient.
 98                sinθ = np.sin(θ)
 99                cosθ = np.cos(θ)
100                target_orientations_diff = np.array(
101                    [[sinθ, cosθ, 0], [-cosθ, sinθ, 0], [0, 0, 0]]
102                )
103                np.testing.assert_allclose(
104                    orientations_diff[0], target_orientations_diff, atol=1e-15
105                )
106                assert np.isclose(np.sum(fractions_diff), 0.0)
107
108
109class TestDislocationCreepOlivineA:
110    """Single-grain A-type olivine analytical rotation rate tests."""
111
112    class_id = "dislocation_creep_OlA"
113
114    @pytest.mark.skipif(_utils.in_ci("win32"), reason="Not equal to tolerance")
115    def test_shear_dvdx_slip_010_100(self, outdir):
116        r"""Single grain of A-type olivine, slip on (010)[100].
117
118        Velocity gradient:
119        $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$
120
121        """
122        test_id = "dvdx_010_100"
123        nondim_velocity_gradient = np.array([[0, 0, 0], [2, 0, 0], [0, 0, 0]])
124        # Strain rate is 0.5*(L + Lᵀ).
125        nondim_strain_rate = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])
126
127        optional_logging = cl.nullcontext()
128        if outdir is not None:
129            optional_logging = _io.logfile_enable(
130                f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.log"
131            )
132            initial_angles = []
133            rotation_rates = []
134            target_rotation_rates = []
135
136        with optional_logging:
137            for θ in np.mgrid[0 : 2 * np.pi : 3600j]:
138                _log.debug("θ (°): %s", np.rad2deg(θ))
139                # Initial grain rotations around Z (anti-clockwise).
140                initial_orientations = Rotation.from_rotvec([[0, 0, θ]])
141                orientation_matrix = initial_orientations[0].as_matrix()
142                slip_invariants = _core._get_slip_invariants(
143                    nondim_strain_rate, orientation_matrix
144                )
145                _log.debug("slip invariants: %s", slip_invariants)
146
147                crss = _core.get_crss(
148                    _core.MineralPhase.olivine,
149                    _core.MineralFabric.olivine_A,
150                )
151                slip_indices = np.argsort(np.abs(slip_invariants / crss))
152                slip_system = _minerals.OLIVINE_SLIP_SYSTEMS[slip_indices[-1]]
153                assert slip_system == ([0, 1, 0], [1, 0, 0])
154
155                slip_rates = _core._get_slip_rates_olivine(
156                    slip_invariants, slip_indices, crss, 3.5
157                )
158                _log.debug("slip rates: %s", slip_rates)
159
160                deformation_rate = _core._get_deformation_rate(
161                    _core.MineralPhase.olivine,
162                    orientation_matrix,
163                    slip_rates,
164                )
165                _log.debug("deformation rate:\n%s", deformation_rate)
166
167                orientations_diff, fractions_diff = _core.derivatives(
168                    regime=_core.DeformationRegime.matrix_dislocation,
169                    phase=_core.MineralPhase.olivine,
170                    fabric=_core.MineralFabric.olivine_A,
171                    n_grains=1,
172                    orientations=initial_orientations.as_matrix(),
173                    fractions=np.array([1.0]),
174                    strain_rate=nondim_strain_rate,
175                    velocity_gradient=nondim_velocity_gradient,
176                    deformation_gradient_spin=np.full((3, 3), np.nan),
177                    stress_exponent=1.5,
178                    deformation_exponent=3.5,
179                    nucleation_efficiency=5,
180                    gbm_mobility=0,
181                    volume_fraction=1.0,
182                )
183                cosθ = np.cos(θ)
184                cos2θ = np.cos(2 * θ)
185                sinθ = np.sin(θ)
186                target_orientations_diff = np.array(
187                    [
188                        [sinθ * (1 + cos2θ), cosθ * (1 + cos2θ), 0],
189                        [-cosθ * (1 + cos2θ), sinθ * (1 + cos2θ), 0],
190                        [0, 0, 0],
191                    ]
192                )
193                np.testing.assert_allclose(
194                    orientations_diff[0], target_orientations_diff
195                )
196                if outdir is not None:
197                    initial_angles.append(np.rad2deg(θ))
198                    rotation_rates.append(
199                        np.sqrt(
200                            orientations_diff[0][0, 0] ** 2
201                            + orientations_diff[0][0, 1] ** 2
202                        )
203                    )
204                    target_rotation_rates.append(1 + cos2θ)
205                assert np.isclose(np.sum(fractions_diff), 0.0)
206
207        if outdir is not None:
208            fig = _vis.figure(figscale=(1, 2 / 3))
209            axs = fig.subplot_mosaic(  # Put the legend on a special axes.
210                [["legend"], ["spin"]], gridspec_kw={"height_ratios": [0.1, 1]}
211            )
212            ax = axs["spin"]
213            fig, ax, colors = _vis.spin(
214                ax,
215                initial_angles[::25],
216                rotation_rates[::25],
217                initial_angles,
218                target_rotation_rates,
219                shear_axis=90,
220            )
221            _utils.redraw_legend(
222                ax,
223                legendax=axs["legend"],
224                loc="upper center",
225                ncols=3,
226                mode="expand",
227            )
228            fig.savefig(
229                _io.resolve_path(f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.pdf")
230            )
231
232    def test_shear_dudz_slip_001_100(self, outdir):
233        r"""Single grain of A-type olivine, slip on (001)[100].
234
235        Velocity gradient:
236        $$\bm{L} = \begin{bmatrix} 0 & 0 & 2 \cr 0 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$
237
238        """
239        test_id = "dudz_001_100"
240        nondim_velocity_gradient = np.array([[0, 0, 2], [0, 0, 0], [0, 0, 0]])
241        # Strain rate is 0.5*(L + Lᵀ).
242        nondim_strain_rate = np.array([[0, 0, 1], [0, 0, 0], [1, 0, 0]])
243
244        optional_logging = cl.nullcontext()
245        if outdir is not None:
246            optional_logging = _io.logfile_enable(
247                f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.log"
248            )
249            initial_angles = []
250            rotation_rates = []
251            target_rotation_rates = []
252
253        with optional_logging:
254            for θ in np.mgrid[0 : 2 * np.pi : 360j]:
255                _log.debug("θ (°): %s", np.rad2deg(θ))
256                # Initial grain rotations around Y (anti-clockwise).
257                initial_orientations = Rotation.from_rotvec([[0, θ, 0]])
258                orientation_matrix = initial_orientations[0].as_matrix()
259                slip_invariants = _core._get_slip_invariants(
260                    nondim_strain_rate, orientation_matrix
261                )
262                _log.debug("slip invariants: %s", slip_invariants)
263
264                crss = _core.get_crss(
265                    _core.MineralPhase.olivine,
266                    _core.MineralFabric.olivine_A,
267                )
268                slip_indices = np.argsort(np.abs(slip_invariants / crss))
269                slip_system = _minerals.OLIVINE_SLIP_SYSTEMS[slip_indices[-1]]
270                assert slip_system == ([0, 0, 1], [1, 0, 0])
271
272                slip_rates = _core._get_slip_rates_olivine(
273                    slip_invariants, slip_indices, crss, 3.5
274                )
275                _log.debug("slip rates: %s", slip_rates)
276
277                deformation_rate = _core._get_deformation_rate(
278                    _core.MineralPhase.olivine,
279                    orientation_matrix,
280                    slip_rates,
281                )
282                _log.debug("deformation rate:\n%s", deformation_rate)
283
284                orientations_diff, fractions_diff = _core.derivatives(
285                    regime=_core.DeformationRegime.matrix_dislocation,
286                    phase=_core.MineralPhase.olivine,
287                    fabric=_core.MineralFabric.olivine_A,
288                    n_grains=1,
289                    orientations=initial_orientations.as_matrix(),
290                    fractions=np.array([1.0]),
291                    strain_rate=nondim_strain_rate,
292                    velocity_gradient=nondim_velocity_gradient,
293                    deformation_gradient_spin=np.full((3, 3), np.nan),
294                    stress_exponent=1.5,
295                    deformation_exponent=3.5,
296                    nucleation_efficiency=5,
297                    gbm_mobility=0,
298                    volume_fraction=1.0,
299                )
300                cosθ = np.cos(θ)
301                cos2θ = np.cos(2 * θ)
302                sinθ = np.sin(θ)
303                target_orientations_diff = np.array(
304                    [
305                        [-sinθ * (cos2θ - 1), 0, cosθ * (cos2θ - 1)],
306                        [0, 0, 0],
307                        [cosθ * (1 - cos2θ), 0, sinθ * (1 - cos2θ)],
308                    ]
309                )
310                np.testing.assert_allclose(
311                    orientations_diff[0], target_orientations_diff
312                )
313                if outdir is not None:
314                    initial_angles.append(np.rad2deg(θ))
315                    rotation_rates.append(
316                        np.sqrt(
317                            orientations_diff[0][0, 0] ** 2
318                            + orientations_diff[0][0, 2] ** 2
319                        )
320                    )
321                    target_rotation_rates.append(1 - cos2θ)
322                assert np.isclose(np.sum(fractions_diff), 0.0)
323
324        if outdir is not None:
325            fig = _vis.figure(figscale=(1, 2 / 3))
326            axs = fig.subplot_mosaic(  # Put the legend on a special axes.
327                [["legend"], ["spin"]], gridspec_kw={"height_ratios": [0.1, 1]}
328            )
329            ax = axs["spin"]
330            fig, ax, colors = _vis.spin(
331                ax,
332                initial_angles[::5],
333                rotation_rates[::5],
334                initial_angles,
335                target_rotation_rates,
336                shear_axis=0,
337            )
338            _utils.redraw_legend(
339                ax,
340                legendax=axs["legend"],
341                loc="upper center",
342                ncols=3,
343                mode="expand",
344            )
345            fig.savefig(
346                _io.resolve_path(f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.pdf")
347            )
348
349    @pytest.mark.skipif(_utils.in_ci("win32"), reason="Not equal to tolerance")
350    def test_shear_dwdx_slip_001_100(self, outdir):
351        r"""Single grain of A-type olivine, slip on (001)[100].
352
353        Velocity gradient:
354        $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 0 & 0 & 0 \cr 2 & 0 & 0 \end{bmatrix}$$
355
356        """
357        test_id = "dwdx_001_100"
358        nondim_velocity_gradient = np.array([[0, 0, 0], [0, 0, 0], [2, 0, 0]])
359        # Strain rate is 0.5*(L + Lᵀ).
360        nondim_strain_rate = np.array([[0, 0, 1], [0, 0, 0], [1, 0, 0]])
361
362        optional_logging = cl.nullcontext()
363        if outdir is not None:
364            optional_logging = _io.logfile_enable(
365                f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.log"
366            )
367            initial_angles = []
368            rotation_rates = []
369            target_rotation_rates = []
370
371        with optional_logging:
372            for θ in np.mgrid[0 : 2 * np.pi : 360j]:
373                _log.debug("θ (°): %s", np.rad2deg(θ))
374                # Initial grain rotations around Y (anti-clockwise).
375                initial_orientations = Rotation.from_rotvec([[0, θ, 0]])
376                orientation_matrix = initial_orientations[0].as_matrix()
377                slip_invariants = _core._get_slip_invariants(
378                    nondim_strain_rate, orientation_matrix
379                )
380                _log.debug("slip invariants: %s", slip_invariants)
381
382                crss = _core.get_crss(
383                    _core.MineralPhase.olivine,
384                    _core.MineralFabric.olivine_A,
385                )
386                slip_indices = np.argsort(np.abs(slip_invariants / crss))
387                slip_system = _minerals.OLIVINE_SLIP_SYSTEMS[slip_indices[-1]]
388                assert slip_system == ([0, 0, 1], [1, 0, 0])
389
390                slip_rates = _core._get_slip_rates_olivine(
391                    slip_invariants, slip_indices, crss, 3.5
392                )
393                _log.debug("slip rates: %s", slip_rates)
394
395                deformation_rate = _core._get_deformation_rate(
396                    _core.MineralPhase.olivine,
397                    orientation_matrix,
398                    slip_rates,
399                )
400                _log.debug("deformation rate:\n%s", deformation_rate)
401
402                orientations_diff, fractions_diff = _core.derivatives(
403                    regime=_core.DeformationRegime.matrix_dislocation,
404                    phase=_core.MineralPhase.olivine,
405                    fabric=_core.MineralFabric.olivine_A,
406                    n_grains=1,
407                    orientations=initial_orientations.as_matrix(),
408                    fractions=np.array([1.0]),
409                    strain_rate=nondim_strain_rate,
410                    velocity_gradient=nondim_velocity_gradient,
411                    deformation_gradient_spin=np.full((3, 3), np.nan),
412                    stress_exponent=1.5,
413                    deformation_exponent=3.5,
414                    nucleation_efficiency=5,
415                    gbm_mobility=0,
416                    volume_fraction=1.0,
417                )
418                cosθ = np.cos(θ)
419                cos2θ = np.cos(2 * θ)
420                sinθ = np.sin(θ)
421                target_orientations_diff = np.array(
422                    [
423                        [-sinθ * (1 + cos2θ), 0, cosθ * (1 + cos2θ)],
424                        [0, 0, 0],
425                        [-cosθ * (1 + cos2θ), 0, -sinθ * (1 + cos2θ)],
426                    ]
427                )
428                np.testing.assert_allclose(
429                    orientations_diff[0], target_orientations_diff
430                )
431                if outdir is not None:
432                    initial_angles.append(np.rad2deg(θ))
433                    rotation_rates.append(
434                        np.sqrt(
435                            orientations_diff[0][0, 0] ** 2
436                            + orientations_diff[0][0, 2] ** 2
437                        )
438                    )
439                    target_rotation_rates.append(1 + cos2θ)
440                assert np.isclose(np.sum(fractions_diff), 0.0)
441
442        if outdir is not None:
443            fig = _vis.figure(figscale=(1, 2 / 3))
444            axs = fig.subplot_mosaic(  # Put the legend on a special axes.
445                [["legend"], ["spin"]], gridspec_kw={"height_ratios": [0.1, 1]}
446            )
447            ax = axs["spin"]
448            fig, ax, colors = _vis.spin(
449                ax,
450                initial_angles[::5],
451                rotation_rates[::5],
452                initial_angles,
453                target_rotation_rates,
454                shear_axis=90,
455            )
456            _utils.redraw_legend(
457                ax,
458                legendax=axs["legend"],
459                loc="upper center",
460                ncols=3,
461                mode="expand",
462            )
463            fig.savefig(
464                _io.resolve_path(f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.pdf")
465            )
466
467    @pytest.mark.skipif(_utils.in_ci("win32"), reason="Not equal to tolerance")
468    def test_shear_dvdz_slip_010_001(self, outdir):
469        r"""Single grain of A-type olivine, slip on (010)[001].
470
471        Velocity gradient:
472        $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 0 & 0 & 2 \cr 0 & 0 & 0 \end{bmatrix}$$
473
474        """
475        test_id = "dvdz_010_001"
476        nondim_velocity_gradient = np.array([[0, 0, 0], [0, 0, 2], [0, 0, 0]])
477        # Strain rate is 0.5*(L + Lᵀ).
478        nondim_strain_rate = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0]])
479
480        optional_logging = cl.nullcontext()
481        if outdir is not None:
482            optional_logging = _io.logfile_enable(
483                f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.log"
484            )
485            initial_angles = []
486            rotation_rates = []
487            target_rotation_rates = []
488
489        with optional_logging:
490            for θ in np.mgrid[0 : 2 * np.pi : 360j]:
491                _log.debug("θ (°): %s", np.rad2deg(θ))
492                # Initial grain rotations around X (anti-clockwise).
493                initial_orientations = Rotation.from_rotvec([[θ, 0, 0]])
494                orientation_matrix = initial_orientations[0].as_matrix()
495                slip_invariants = _core._get_slip_invariants(
496                    nondim_strain_rate, orientation_matrix
497                )
498                _log.debug("slip invariants: %s", slip_invariants)
499
500                crss = _core.get_crss(
501                    _core.MineralPhase.olivine,
502                    _core.MineralFabric.olivine_A,
503                )
504                slip_indices = np.argsort(np.abs(slip_invariants / crss))
505                slip_system = _minerals.OLIVINE_SLIP_SYSTEMS[slip_indices[-1]]
506                assert slip_system == ([0, 1, 0], [0, 0, 1])
507
508                slip_rates = _core._get_slip_rates_olivine(
509                    slip_invariants, slip_indices, crss, 3.5
510                )
511                _log.debug("slip rates: %s", slip_rates)
512
513                deformation_rate = _core._get_deformation_rate(
514                    _core.MineralPhase.olivine,
515                    orientation_matrix,
516                    slip_rates,
517                )
518                _log.debug("deformation rate:\n%s", deformation_rate)
519
520                orientations_diff, fractions_diff = _core.derivatives(
521                    regime=_core.DeformationRegime.matrix_dislocation,
522                    phase=_core.MineralPhase.olivine,
523                    fabric=_core.MineralFabric.olivine_A,
524                    n_grains=1,
525                    orientations=initial_orientations.as_matrix(),
526                    fractions=np.array([1.0]),
527                    strain_rate=nondim_strain_rate,
528                    velocity_gradient=nondim_velocity_gradient,
529                    deformation_gradient_spin=np.full((3, 3), np.nan),
530                    stress_exponent=1.5,
531                    deformation_exponent=3.5,
532                    nucleation_efficiency=5,
533                    gbm_mobility=0,
534                    volume_fraction=1.0,
535                )
536                cosθ = np.cos(θ)
537                cos2θ = np.cos(2 * θ)
538                sinθ = np.sin(θ)
539                target_orientations_diff = np.array(
540                    [
541                        [0, 0, 0],
542                        [0, -sinθ * (1 + cos2θ), -cosθ * (1 + cos2θ)],
543                        [0, cosθ * (1 + cos2θ), -sinθ * (1 + cos2θ)],
544                    ]
545                )
546                np.testing.assert_allclose(
547                    orientations_diff[0], target_orientations_diff
548                )
549                if outdir is not None:
550                    initial_angles.append(np.rad2deg(θ))
551                    rotation_rates.append(
552                        np.sqrt(
553                            orientations_diff[0][1, 1] ** 2
554                            + orientations_diff[0][1, 2] ** 2
555                        )
556                    )
557                    target_rotation_rates.append(1 + cos2θ)
558                assert np.isclose(np.sum(fractions_diff), 0.0)
559
560        if outdir is not None:
561            fig = _vis.figure(figscale=(1, 2 / 3))
562            axs = fig.subplot_mosaic(  # Put the legend on a special axes.
563                [["legend"], ["spin"]], gridspec_kw={"height_ratios": [0.1, 1]}
564            )
565            ax = axs["spin"]
566            fig, ax, colors = _vis.spin(
567                ax,
568                initial_angles[::5],
569                rotation_rates[::5],
570                initial_angles,
571                target_rotation_rates,
572                shear_axis=90,
573            )
574            _utils.redraw_legend(
575                ax,
576                legendax=axs["legend"],
577                loc="upper center",
578                ncols=3,
579                mode="expand",
580            )
581            fig.savefig(
582                _io.resolve_path(f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}.pdf")
583            )
584
585
586class TestRecrystallisation2D:
587    """Basic recrystallisation tests for 2D simple shear."""
588
589    class_id = "recrystallisation_2D"
590
591    def test_shear_dvdx_circle_inplane(self, outdir):
592        r"""360000 grains of A-type olivine with uniform spread of a-axes on a circle.
593
594        Grain growth rates are compared to analytical calculations.
595        The a-axes are distributed in the YX plane (i.e.\ rotated around Z).
596
597        Velocity gradient:
598        $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$
599
600        """
601        test_id = "dvdx_circle_inplane"
602
603        optional_logging = cl.nullcontext()
604        # Initial uniform distribution of orientations on a 2D circle.
605        initial_angles = np.mgrid[0 : 2 * np.pi : 360000j]
606        cos2θ = np.cos(2 * initial_angles)
607        if outdir is not None:
608            out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}"
609            optional_logging = _io.logfile_enable(f"{out_basepath}.log")
610
611        with optional_logging:
612            initial_orientations = Rotation.from_rotvec(
613                [[0, 0, θ] for θ in initial_angles]
614            )
615            orientations_diff, fractions_diff = _core.derivatives(
616                regime=_core.DeformationRegime.matrix_dislocation,
617                phase=_core.MineralPhase.olivine,
618                fabric=_core.MineralFabric.olivine_A,
619                n_grains=360000,
620                orientations=initial_orientations.as_matrix(),
621                fractions=np.full(360000, 1 / 360000),
622                strain_rate=np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]]),
623                velocity_gradient=np.array([[0, 0, 0], [2, 0, 0], [0, 0, 0]]),
624                deformation_gradient_spin=np.full((3, 3), np.nan),
625                stress_exponent=1.5,
626                deformation_exponent=3.5,
627                nucleation_efficiency=5,
628                gbm_mobility=125,
629                volume_fraction=1.0,
630            )
631            ρ = np.abs(cos2θ) ** (1.5 / 3.5)
632            # No need to sum over slip systems, only (010)[100] can activate in this
633            # 2D simple shear deformation geometry (ρₛ=ρ).
634            target_strain_energies = ρ * np.exp(-5 * ρ**2)
635            target_fractions_diff = np.array(
636                [  # df/dt* = - M* f (E - E_mean)
637                    -125 * 1 / 360000 * (E - np.mean(target_strain_energies))
638                    for E in target_strain_energies
639                ]
640            )
641
642        if outdir is not None:
643            fig = _vis.figure(figscale=(1, 4 / 3))
644            axs = fig.subplot_mosaic(  # Put the legend on a special axes.
645                [["legend"], ["spin"], ["growth"]],
646                gridspec_kw={"height_ratios": [0.2, 1, 1]},
647                sharex=True,
648            )
649            ax = axs["spin"]
650            xvals = np.rad2deg(initial_angles)
651            fig, ax, colors = _vis.spin(
652                ax,
653                xvals[::2500],
654                np.sqrt(
655                    [
656                        o[0, 0] ** 2 + o[0, 1] ** 2 + o[0, 2] ** 2
657                        for o in orientations_diff
658                    ][::2500]
659                ),
660                xvals,
661                1 + cos2θ,
662                shear_axis=90,
663            )
664            ax.get_legend().remove()
665            ax.label_outer()
666            ax2 = axs["growth"]
667            fig, ax2, colors = _vis.growth(
668                ax2,
669                xvals[::2500],
670                fractions_diff[::2500],
671                xvals,
672                target_fractions_diff,
673                shear_axis=90,
674            )
675            _utils.redraw_legend(
676                ax2,
677                legendax=axs["legend"],
678                loc="upper center",
679                ncols=3,
680                mode="expand",
681            )
682            _utils.add_subplot_labels(axs, labelmap={"spin": "a)", "growth": "b)"})
683            fig.savefig(f"{out_basepath}.pdf")
684
685        nt.assert_allclose(fractions_diff, target_fractions_diff, atol=1e-15, rtol=0)
686
687    def test_shear_dvdx_circle_shearplane(self, outdir):
688        r"""360000 grains of A-type olivine with uniform spread of a-axes on a circle.
689
690        Unlike `test_shear_dvdx_circle_inplane`, two slip systems are active here,
691        with cyclical variety in which one is dominant depending on grain orientation.
692        The a-axes are distributed in the YZ plane
693        (i.e.\ extrinsic rotation around Z by 90° and then around X).
694
695        Velocity gradient:
696        $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$
697
698        """
699        test_id = "dvdx_circle_shearplane"
700
701        optional_logging = cl.nullcontext()
702        # Initial uniform distribution of orientations on a 2D circle.
703        initial_angles = np.mgrid[0 : 2 * np.pi : 360000j]
704        if outdir is not None:
705            out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}"
706            optional_logging = _io.logfile_enable(f"{out_basepath}.log")
707
708        with optional_logging:
709            initial_orientations = Rotation.from_euler(
710                "zx", [[np.pi / 2, θ] for θ in initial_angles]
711            )
712            orientations_diff, fractions_diff = _core.derivatives(
713                regime=_core.DeformationRegime.matrix_dislocation,
714                phase=_core.MineralPhase.olivine,
715                fabric=_core.MineralFabric.olivine_A,
716                n_grains=360000,
717                orientations=initial_orientations.as_matrix(),
718                fractions=np.full(360000, 1 / 360000),
719                strain_rate=np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]]),
720                velocity_gradient=np.array([[0, 0, 0], [2, 0, 0], [0, 0, 0]]),
721                deformation_gradient_spin=np.full((3, 3), np.nan),
722                stress_exponent=1.5,
723                deformation_exponent=3.5,
724                nucleation_efficiency=5,
725                gbm_mobility=125,
726                volume_fraction=1.0,
727            )
728
729        if outdir is not None:
730            fig = _vis.figure(figscale=(1, 4 / 3))
731            axs = fig.subplot_mosaic([["spin"], ["growth"]], sharex=True)
732            ax = axs["spin"]
733            xvals = np.rad2deg(initial_angles)
734            fig, ax, colors = _vis.spin(
735                ax,
736                xvals[::2500],
737                np.sqrt(
738                    [
739                        o[0, 0] ** 2 + o[0, 1] ** 2 + o[0, 2] ** 2
740                        for o in orientations_diff
741                    ][::2500]
742                ),
743            )
744            ax.get_legend().remove()
745            ax.label_outer()
746            ax2 = axs["growth"]
747            fig, ax2, colors = _vis.growth(ax2, xvals[::2500], fractions_diff[::2500])
748            _utils.add_subplot_labels(axs, labelmap={"spin": "a)", "growth": "b)"})
749            fig.savefig(f"{out_basepath}.pdf")
750
751        # Check dominant slip system every 1°.
752        for θ in initial_angles[::1000]:
753            slip_invariants = _core._get_slip_invariants(
754                np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]]),
755                Rotation.from_euler("zx", [np.pi / 2, θ]).as_matrix(),
756            )
757            θ = np.rad2deg(θ)
758            crss = _core.get_crss(
759                _core.MineralPhase.olivine,
760                _core.MineralFabric.olivine_A,
761            )
762            slip_indices = np.argsort(np.abs(slip_invariants / crss))
763            slip_system = _minerals.OLIVINE_SLIP_SYSTEMS[slip_indices[-1]]
764
765            if 0 <= θ < 64:
766                assert slip_system == ([0, 1, 0], [1, 0, 0])
767            elif 64 <= θ < 117:
768                assert slip_system == ([0, 0, 1], [1, 0, 0])
769            elif 117 <= θ < 244:
770                assert slip_system == ([0, 1, 0], [1, 0, 0])
771            elif 244 <= θ < 297:
772                assert slip_system == ([0, 0, 1], [1, 0, 0])
773            else:
774                assert slip_system == ([0, 1, 0], [1, 0, 0])