tests.test_corner_flow_2d

PyDRex: 2D corner flow tests.

  1"""> PyDRex: 2D corner flow tests."""
  2
  3import contextlib as cl
  4import functools as ft
  5import pathlib as pl
  6from time import perf_counter
  7
  8import numpy as np
  9import pytest
 10
 11from pydrex import core as _core
 12from pydrex import diagnostics as _diagnostics
 13from pydrex import geometry as _geo
 14from pydrex import io as _io
 15from pydrex import logger as _log
 16from pydrex import minerals as _minerals
 17from pydrex import pathlines as _path
 18from pydrex import stats as _stats
 19from pydrex import utils as _utils
 20from pydrex import velocity as _velocity
 21from pydrex import visualisation as _vis
 22
 23Pool, HAS_RAY = _utils.import_proc_pool()
 24
 25# Subdirectory of `outdir` used to store outputs from these tests.
 26SUBDIR = "2d_cornerflow"
 27
 28
 29class TestOlivineA:
 30    """Tests for pure A-type olivine polycrystals in 2D corner flows."""
 31
 32    class_id = "corner_olivineA"
 33
 34    @classmethod
 35    def run(
 36        cls,
 37        params,
 38        seed,
 39        get_velocity,
 40        get_velocity_gradient,
 41        min_coords,
 42        max_coords,
 43        max_strain,
 44        n_timesteps,
 45        final_location,
 46    ):
 47        """Run 2D corner flow A-type olivine simulation."""
 48        mineral = _minerals.Mineral(
 49            _core.MineralPhase.olivine,
 50            _core.MineralFabric.olivine_A,
 51            _core.DeformationRegime.matrix_dislocation,
 52            n_grains=params["number_of_grains"],
 53            seed=seed,
 54        )
 55        deformation_gradient = np.eye(3)
 56        timestamps, get_position = _path.get_pathline(
 57            final_location,
 58            get_velocity,
 59            get_velocity_gradient,
 60            min_coords,
 61            max_coords,
 62            max_strain,
 63            regular_steps=n_timesteps,
 64        )
 65        positions = [get_position(t) for t in timestamps]
 66        velocity_gradients = [
 67            get_velocity_gradient(np.nan, np.asarray(x)) for x in positions
 68        ]
 69        strains = np.empty_like(timestamps)
 70        strains[0] = 0
 71        for t, time in enumerate(timestamps[:-1], start=1):
 72            strains[t] = strains[t - 1] + (
 73                _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
 74            )
 75            _log.info(
 76                "final location = %s; step %d/%d (ε = %.2f)",
 77                final_location.ravel(),
 78                t,
 79                len(timestamps) - 1,
 80                strains[t],
 81            )
 82
 83            deformation_gradient = mineral.update_orientations(
 84                params,
 85                deformation_gradient,
 86                get_velocity_gradient,
 87                pathline=(time, timestamps[t], get_position),
 88            )
 89        return timestamps, positions, strains, mineral, deformation_gradient
 90
 91    @pytest.mark.slow
 92    def test_steady4(self, outdir, seed, ncpus):
 93        """Test CPO evolution in steady 2D corner flow along 4 pathlines.
 94
 95        Initial condition: random orientations and uniform volumes in all `Mineral`s.
 96
 97        Plate velocity: 2 cm/yr
 98
 99        .. note::
100            This example takes about 11 CPU hours to run and uses around 60GB of RAM.
101            It is recommended to only use `ncpus=4` which matches the number of
102            pathlines, because higher numbers can lead to redundant cross-core
103            communication.
104
105        """
106        # Plate speed (half spreading rate), convert cm/yr to m/s.
107        plate_speed = 2.0 / (100.0 * 365.0 * 86400.0)
108        domain_height = 2.0e5  # Represents the depth of olivine-spinel transition.
109        domain_width = 1.0e6
110        params = _core.DefaultParams().as_dict()
111        params["number_of_grains"] = 5000
112        n_timesteps = 50  # Number of places along the pathline to compute CPO.
113        get_velocity, get_velocity_gradient = _velocity.corner_2d("X", "Z", plate_speed)
114        # Z-values at the end of each pathline.
115        z_ends = list(map(lambda x: x * domain_height, (-0.1, -0.3, -0.54, -0.78)))
116
117        # Optional plotting and logging setup.
118        optional_logging = cl.nullcontext()
119        if outdir is not None:
120            out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_prescribed"
121            optional_logging = _io.logfile_enable(f"{out_basepath}.log")
122            npzpath = pl.Path(f"{out_basepath}.npz")
123            labels = []
124            angles = []
125            indices = []
126            paths = []
127            path_strains = []
128            directions = []
129            max_grainsizes = []
130            timestamps = []
131
132        _run = ft.partial(
133            self.run,
134            params,
135            seed,
136            get_velocity,
137            get_velocity_gradient,
138            np.array([0.0, 0.0, -domain_height]),
139            np.array([domain_width, 0.0, 0.0]),
140            10,
141            n_timesteps,
142        )
143
144        final_locations = [np.array([domain_width, 0.0, z_exit]) for z_exit in z_ends]
145        with optional_logging:
146            _begin = perf_counter()
147            with Pool(processes=ncpus) as pool:
148                for i, out in enumerate(pool.imap(_run, final_locations)):
149                    times, positions, strains, mineral, deformation_gradient = out
150                    _log.info("final deformation gradient:\n%s", deformation_gradient)
151                    z_exit = z_ends[i]
152                    if outdir is not None:
153                        mineral.save(npzpath, _io.stringify(z_exit))
154
155                    misorient_angles = np.zeros(n_timesteps)
156                    bingham_vectors = np.zeros((n_timesteps, 3))
157                    orientations_resampled, f_resampled = _stats.resample_orientations(
158                        mineral.orientations, mineral.fractions, seed=seed
159                    )
160                    misorient_indices = _diagnostics.misorientation_indices(
161                        orientations_resampled,
162                        _geo.LatticeSystem.orthorhombic,
163                        pool=pool,
164                    )
165                    for idx, matrices in enumerate(orientations_resampled):
166                        direction_mean = _diagnostics.bingham_average(
167                            matrices,
168                            axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric],
169                        )
170                        misorient_angles[idx] = _diagnostics.smallest_angle(
171                            direction_mean,
172                            np.asarray([1, 0, 0], dtype=np.float64),
173                        )
174                        bingham_vectors[idx] = direction_mean
175
176                    _log.debug("Total walltime: %s", perf_counter() - _begin)
177
178                    if outdir is not None:
179                        labels.append(rf"$z_{{f}}$ = {z_exit/1e3:.1f} km")
180                        angles.append(misorient_angles)
181                        indices.append(misorient_indices)
182                        paths.append(positions)
183                        max_grainsizes.append(
184                            np.asarray([np.max(f) for f in f_resampled])
185                        )
186                        path_strains.append(strains)
187                        timestamps.append(times)
188                        directions.append(bingham_vectors)
189
190        if outdir is not None:
191            np.savez(
192                f"{out_basepath}_path.npz",
193                strains_1=path_strains[0],
194                strains_2=path_strains[1],
195                strains_3=path_strains[2],
196                strains_4=path_strains[3],
197                paths_1=np.asarray(paths[0]),
198                paths_2=np.asarray(paths[1]),
199                paths_3=np.asarray(paths[2]),
200                paths_4=np.asarray(paths[3]),
201            )
202            np.savez(
203                f"{out_basepath}_strength.npz",
204                strength_1=indices[0],
205                strength_2=indices[1],
206                strength_3=indices[2],
207                strength_4=indices[3],
208            )
209            np.savez(
210                f"{out_basepath}_angles.npz",
211                angles_1=angles[0],
212                angles_2=angles[1],
213                angles_3=angles[2],
214                angles_4=angles[3],
215            )
216            markers = ("s", "o", "v", "*")
217            cmaps = ["cmc.batlow_r"] * len(markers)
218            fig_domain = _vis.figure(figsize=(10, 3))
219            ax_domain = fig_domain.add_subplot()
220            for i, z_exit in enumerate(z_ends):
221                if i == 0:
222                    resolution = [25, 5]
223                else:
224                    resolution = None
225                _vis.pathline_box2d(
226                    ax_domain,
227                    get_velocity,
228                    "xz",
229                    path_strains[i],
230                    paths[i],
231                    markers[i],
232                    [0, -domain_height],
233                    [domain_width + 1e5, 0],
234                    resolution,
235                    cmap=cmaps[i],
236                    cpo_vectors=directions[i],
237                    cpo_strengths=indices[i],
238                    scale=25,
239                    scale_units="width",
240                )
241
242            fig_strength, ax_strength, colors = _vis.strengths(
243                None,
244                path_strains,
245                indices,
246                "CPO strength (M-index)",
247                markers,
248                labels,
249                colors=path_strains,
250                cmaps=cmaps,
251            )
252            fig_alignment, ax_alignment, colors = _vis.alignment(
253                None,
254                path_strains,
255                angles,
256                markers,
257                labels,
258                colors=path_strains,
259                cmaps=cmaps,
260            )
261
262            fig_domain.savefig(_io.resolve_path(f"{out_basepath}_path.pdf"))
263            fig_strength.savefig(_io.resolve_path(f"{out_basepath}_strength.pdf"))
264            fig_alignment.savefig(_io.resolve_path(f"{out_basepath}_angles.pdf"))
SUBDIR = '2d_cornerflow'
class TestOlivineA: 
 30class TestOlivineA:
 31    """Tests for pure A-type olivine polycrystals in 2D corner flows."""
 32
 33    class_id = "corner_olivineA"
 34
 35    @classmethod
 36    def run(
 37        cls,
 38        params,
 39        seed,
 40        get_velocity,
 41        get_velocity_gradient,
 42        min_coords,
 43        max_coords,
 44        max_strain,
 45        n_timesteps,
 46        final_location,
 47    ):
 48        """Run 2D corner flow A-type olivine simulation."""
 49        mineral = _minerals.Mineral(
 50            _core.MineralPhase.olivine,
 51            _core.MineralFabric.olivine_A,
 52            _core.DeformationRegime.matrix_dislocation,
 53            n_grains=params["number_of_grains"],
 54            seed=seed,
 55        )
 56        deformation_gradient = np.eye(3)
 57        timestamps, get_position = _path.get_pathline(
 58            final_location,
 59            get_velocity,
 60            get_velocity_gradient,
 61            min_coords,
 62            max_coords,
 63            max_strain,
 64            regular_steps=n_timesteps,
 65        )
 66        positions = [get_position(t) for t in timestamps]
 67        velocity_gradients = [
 68            get_velocity_gradient(np.nan, np.asarray(x)) for x in positions
 69        ]
 70        strains = np.empty_like(timestamps)
 71        strains[0] = 0
 72        for t, time in enumerate(timestamps[:-1], start=1):
 73            strains[t] = strains[t - 1] + (
 74                _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
 75            )
 76            _log.info(
 77                "final location = %s; step %d/%d (ε = %.2f)",
 78                final_location.ravel(),
 79                t,
 80                len(timestamps) - 1,
 81                strains[t],
 82            )
 83
 84            deformation_gradient = mineral.update_orientations(
 85                params,
 86                deformation_gradient,
 87                get_velocity_gradient,
 88                pathline=(time, timestamps[t], get_position),
 89            )
 90        return timestamps, positions, strains, mineral, deformation_gradient
 91
 92    @pytest.mark.slow
 93    def test_steady4(self, outdir, seed, ncpus):
 94        """Test CPO evolution in steady 2D corner flow along 4 pathlines.
 95
 96        Initial condition: random orientations and uniform volumes in all `Mineral`s.
 97
 98        Plate velocity: 2 cm/yr
 99
100        .. note::
101            This example takes about 11 CPU hours to run and uses around 60GB of RAM.
102            It is recommended to only use `ncpus=4` which matches the number of
103            pathlines, because higher numbers can lead to redundant cross-core
104            communication.
105
106        """
107        # Plate speed (half spreading rate), convert cm/yr to m/s.
108        plate_speed = 2.0 / (100.0 * 365.0 * 86400.0)
109        domain_height = 2.0e5  # Represents the depth of olivine-spinel transition.
110        domain_width = 1.0e6
111        params = _core.DefaultParams().as_dict()
112        params["number_of_grains"] = 5000
113        n_timesteps = 50  # Number of places along the pathline to compute CPO.
114        get_velocity, get_velocity_gradient = _velocity.corner_2d("X", "Z", plate_speed)
115        # Z-values at the end of each pathline.
116        z_ends = list(map(lambda x: x * domain_height, (-0.1, -0.3, -0.54, -0.78)))
117
118        # Optional plotting and logging setup.
119        optional_logging = cl.nullcontext()
120        if outdir is not None:
121            out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_prescribed"
122            optional_logging = _io.logfile_enable(f"{out_basepath}.log")
123            npzpath = pl.Path(f"{out_basepath}.npz")
124            labels = []
125            angles = []
126            indices = []
127            paths = []
128            path_strains = []
129            directions = []
130            max_grainsizes = []
131            timestamps = []
132
133        _run = ft.partial(
134            self.run,
135            params,
136            seed,
137            get_velocity,
138            get_velocity_gradient,
139            np.array([0.0, 0.0, -domain_height]),
140            np.array([domain_width, 0.0, 0.0]),
141            10,
142            n_timesteps,
143        )
144
145        final_locations = [np.array([domain_width, 0.0, z_exit]) for z_exit in z_ends]
146        with optional_logging:
147            _begin = perf_counter()
148            with Pool(processes=ncpus) as pool:
149                for i, out in enumerate(pool.imap(_run, final_locations)):
150                    times, positions, strains, mineral, deformation_gradient = out
151                    _log.info("final deformation gradient:\n%s", deformation_gradient)
152                    z_exit = z_ends[i]
153                    if outdir is not None:
154                        mineral.save(npzpath, _io.stringify(z_exit))
155
156                    misorient_angles = np.zeros(n_timesteps)
157                    bingham_vectors = np.zeros((n_timesteps, 3))
158                    orientations_resampled, f_resampled = _stats.resample_orientations(
159                        mineral.orientations, mineral.fractions, seed=seed
160                    )
161                    misorient_indices = _diagnostics.misorientation_indices(
162                        orientations_resampled,
163                        _geo.LatticeSystem.orthorhombic,
164                        pool=pool,
165                    )
166                    for idx, matrices in enumerate(orientations_resampled):
167                        direction_mean = _diagnostics.bingham_average(
168                            matrices,
169                            axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric],
170                        )
171                        misorient_angles[idx] = _diagnostics.smallest_angle(
172                            direction_mean,
173                            np.asarray([1, 0, 0], dtype=np.float64),
174                        )
175                        bingham_vectors[idx] = direction_mean
176
177                    _log.debug("Total walltime: %s", perf_counter() - _begin)
178
179                    if outdir is not None:
180                        labels.append(rf"$z_{{f}}$ = {z_exit/1e3:.1f} km")
181                        angles.append(misorient_angles)
182                        indices.append(misorient_indices)
183                        paths.append(positions)
184                        max_grainsizes.append(
185                            np.asarray([np.max(f) for f in f_resampled])
186                        )
187                        path_strains.append(strains)
188                        timestamps.append(times)
189                        directions.append(bingham_vectors)
190
191        if outdir is not None:
192            np.savez(
193                f"{out_basepath}_path.npz",
194                strains_1=path_strains[0],
195                strains_2=path_strains[1],
196                strains_3=path_strains[2],
197                strains_4=path_strains[3],
198                paths_1=np.asarray(paths[0]),
199                paths_2=np.asarray(paths[1]),
200                paths_3=np.asarray(paths[2]),
201                paths_4=np.asarray(paths[3]),
202            )
203            np.savez(
204                f"{out_basepath}_strength.npz",
205                strength_1=indices[0],
206                strength_2=indices[1],
207                strength_3=indices[2],
208                strength_4=indices[3],
209            )
210            np.savez(
211                f"{out_basepath}_angles.npz",
212                angles_1=angles[0],
213                angles_2=angles[1],
214                angles_3=angles[2],
215                angles_4=angles[3],
216            )
217            markers = ("s", "o", "v", "*")
218            cmaps = ["cmc.batlow_r"] * len(markers)
219            fig_domain = _vis.figure(figsize=(10, 3))
220            ax_domain = fig_domain.add_subplot()
221            for i, z_exit in enumerate(z_ends):
222                if i == 0:
223                    resolution = [25, 5]
224                else:
225                    resolution = None
226                _vis.pathline_box2d(
227                    ax_domain,
228                    get_velocity,
229                    "xz",
230                    path_strains[i],
231                    paths[i],
232                    markers[i],
233                    [0, -domain_height],
234                    [domain_width + 1e5, 0],
235                    resolution,
236                    cmap=cmaps[i],
237                    cpo_vectors=directions[i],
238                    cpo_strengths=indices[i],
239                    scale=25,
240                    scale_units="width",
241                )
242
243            fig_strength, ax_strength, colors = _vis.strengths(
244                None,
245                path_strains,
246                indices,
247                "CPO strength (M-index)",
248                markers,
249                labels,
250                colors=path_strains,
251                cmaps=cmaps,
252            )
253            fig_alignment, ax_alignment, colors = _vis.alignment(
254                None,
255                path_strains,
256                angles,
257                markers,
258                labels,
259                colors=path_strains,
260                cmaps=cmaps,
261            )
262
263            fig_domain.savefig(_io.resolve_path(f"{out_basepath}_path.pdf"))
264            fig_strength.savefig(_io.resolve_path(f"{out_basepath}_strength.pdf"))
265            fig_alignment.savefig(_io.resolve_path(f"{out_basepath}_angles.pdf"))

Tests for pure A-type olivine polycrystals in 2D corner flows.

class_id = 'corner_olivineA'
@classmethod
def run( cls, params, seed, get_velocity, get_velocity_gradient, min_coords, max_coords, max_strain, n_timesteps, final_location): 
35    @classmethod
36    def run(
37        cls,
38        params,
39        seed,
40        get_velocity,
41        get_velocity_gradient,
42        min_coords,
43        max_coords,
44        max_strain,
45        n_timesteps,
46        final_location,
47    ):
48        """Run 2D corner flow A-type olivine simulation."""
49        mineral = _minerals.Mineral(
50            _core.MineralPhase.olivine,
51            _core.MineralFabric.olivine_A,
52            _core.DeformationRegime.matrix_dislocation,
53            n_grains=params["number_of_grains"],
54            seed=seed,
55        )
56        deformation_gradient = np.eye(3)
57        timestamps, get_position = _path.get_pathline(
58            final_location,
59            get_velocity,
60            get_velocity_gradient,
61            min_coords,
62            max_coords,
63            max_strain,
64            regular_steps=n_timesteps,
65        )
66        positions = [get_position(t) for t in timestamps]
67        velocity_gradients = [
68            get_velocity_gradient(np.nan, np.asarray(x)) for x in positions
69        ]
70        strains = np.empty_like(timestamps)
71        strains[0] = 0
72        for t, time in enumerate(timestamps[:-1], start=1):
73            strains[t] = strains[t - 1] + (
74                _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
75            )
76            _log.info(
77                "final location = %s; step %d/%d (ε = %.2f)",
78                final_location.ravel(),
79                t,
80                len(timestamps) - 1,
81                strains[t],
82            )
83
84            deformation_gradient = mineral.update_orientations(
85                params,
86                deformation_gradient,
87                get_velocity_gradient,
88                pathline=(time, timestamps[t], get_position),
89            )
90        return timestamps, positions, strains, mineral, deformation_gradient

Run 2D corner flow A-type olivine simulation.

@pytest.mark.slow
def test_steady4(self, outdir, seed, ncpus): 
 92    @pytest.mark.slow
 93    def test_steady4(self, outdir, seed, ncpus):
 94        """Test CPO evolution in steady 2D corner flow along 4 pathlines.
 95
 96        Initial condition: random orientations and uniform volumes in all `Mineral`s.
 97
 98        Plate velocity: 2 cm/yr
 99
100        .. note::
101            This example takes about 11 CPU hours to run and uses around 60GB of RAM.
102            It is recommended to only use `ncpus=4` which matches the number of
103            pathlines, because higher numbers can lead to redundant cross-core
104            communication.
105
106        """
107        # Plate speed (half spreading rate), convert cm/yr to m/s.
108        plate_speed = 2.0 / (100.0 * 365.0 * 86400.0)
109        domain_height = 2.0e5  # Represents the depth of olivine-spinel transition.
110        domain_width = 1.0e6
111        params = _core.DefaultParams().as_dict()
112        params["number_of_grains"] = 5000
113        n_timesteps = 50  # Number of places along the pathline to compute CPO.
114        get_velocity, get_velocity_gradient = _velocity.corner_2d("X", "Z", plate_speed)
115        # Z-values at the end of each pathline.
116        z_ends = list(map(lambda x: x * domain_height, (-0.1, -0.3, -0.54, -0.78)))
117
118        # Optional plotting and logging setup.
119        optional_logging = cl.nullcontext()
120        if outdir is not None:
121            out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_prescribed"
122            optional_logging = _io.logfile_enable(f"{out_basepath}.log")
123            npzpath = pl.Path(f"{out_basepath}.npz")
124            labels = []
125            angles = []
126            indices = []
127            paths = []
128            path_strains = []
129            directions = []
130            max_grainsizes = []
131            timestamps = []
132
133        _run = ft.partial(
134            self.run,
135            params,
136            seed,
137            get_velocity,
138            get_velocity_gradient,
139            np.array([0.0, 0.0, -domain_height]),
140            np.array([domain_width, 0.0, 0.0]),
141            10,
142            n_timesteps,
143        )
144
145        final_locations = [np.array([domain_width, 0.0, z_exit]) for z_exit in z_ends]
146        with optional_logging:
147            _begin = perf_counter()
148            with Pool(processes=ncpus) as pool:
149                for i, out in enumerate(pool.imap(_run, final_locations)):
150                    times, positions, strains, mineral, deformation_gradient = out
151                    _log.info("final deformation gradient:\n%s", deformation_gradient)
152                    z_exit = z_ends[i]
153                    if outdir is not None:
154                        mineral.save(npzpath, _io.stringify(z_exit))
155
156                    misorient_angles = np.zeros(n_timesteps)
157                    bingham_vectors = np.zeros((n_timesteps, 3))
158                    orientations_resampled, f_resampled = _stats.resample_orientations(
159                        mineral.orientations, mineral.fractions, seed=seed
160                    )
161                    misorient_indices = _diagnostics.misorientation_indices(
162                        orientations_resampled,
163                        _geo.LatticeSystem.orthorhombic,
164                        pool=pool,
165                    )
166                    for idx, matrices in enumerate(orientations_resampled):
167                        direction_mean = _diagnostics.bingham_average(
168                            matrices,
169                            axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric],
170                        )
171                        misorient_angles[idx] = _diagnostics.smallest_angle(
172                            direction_mean,
173                            np.asarray([1, 0, 0], dtype=np.float64),
174                        )
175                        bingham_vectors[idx] = direction_mean
176
177                    _log.debug("Total walltime: %s", perf_counter() - _begin)
178
179                    if outdir is not None:
180                        labels.append(rf"$z_{{f}}$ = {z_exit/1e3:.1f} km")
181                        angles.append(misorient_angles)
182                        indices.append(misorient_indices)
183                        paths.append(positions)
184                        max_grainsizes.append(
185                            np.asarray([np.max(f) for f in f_resampled])
186                        )
187                        path_strains.append(strains)
188                        timestamps.append(times)
189                        directions.append(bingham_vectors)
190
191        if outdir is not None:
192            np.savez(
193                f"{out_basepath}_path.npz",
194                strains_1=path_strains[0],
195                strains_2=path_strains[1],
196                strains_3=path_strains[2],
197                strains_4=path_strains[3],
198                paths_1=np.asarray(paths[0]),
199                paths_2=np.asarray(paths[1]),
200                paths_3=np.asarray(paths[2]),
201                paths_4=np.asarray(paths[3]),
202            )
203            np.savez(
204                f"{out_basepath}_strength.npz",
205                strength_1=indices[0],
206                strength_2=indices[1],
207                strength_3=indices[2],
208                strength_4=indices[3],
209            )
210            np.savez(
211                f"{out_basepath}_angles.npz",
212                angles_1=angles[0],
213                angles_2=angles[1],
214                angles_3=angles[2],
215                angles_4=angles[3],
216            )
217            markers = ("s", "o", "v", "*")
218            cmaps = ["cmc.batlow_r"] * len(markers)
219            fig_domain = _vis.figure(figsize=(10, 3))
220            ax_domain = fig_domain.add_subplot()
221            for i, z_exit in enumerate(z_ends):
222                if i == 0:
223                    resolution = [25, 5]
224                else:
225                    resolution = None
226                _vis.pathline_box2d(
227                    ax_domain,
228                    get_velocity,
229                    "xz",
230                    path_strains[i],
231                    paths[i],
232                    markers[i],
233                    [0, -domain_height],
234                    [domain_width + 1e5, 0],
235                    resolution,
236                    cmap=cmaps[i],
237                    cpo_vectors=directions[i],
238                    cpo_strengths=indices[i],
239                    scale=25,
240                    scale_units="width",
241                )
242
243            fig_strength, ax_strength, colors = _vis.strengths(
244                None,
245                path_strains,
246                indices,
247                "CPO strength (M-index)",
248                markers,
249                labels,
250                colors=path_strains,
251                cmaps=cmaps,
252            )
253            fig_alignment, ax_alignment, colors = _vis.alignment(
254                None,
255                path_strains,
256                angles,
257                markers,
258                labels,
259                colors=path_strains,
260                cmaps=cmaps,
261            )
262
263            fig_domain.savefig(_io.resolve_path(f"{out_basepath}_path.pdf"))
264            fig_strength.savefig(_io.resolve_path(f"{out_basepath}_strength.pdf"))
265            fig_alignment.savefig(_io.resolve_path(f"{out_basepath}_angles.pdf"))

Test CPO evolution in steady 2D corner flow along 4 pathlines.

Initial condition: random orientations and uniform volumes in all Minerals.

Plate velocity: 2 cm/yr

This example takes about 11 CPU hours to run and uses around 60GB of RAM. It is recommended to only use ncpus=4 which matches the number of pathlines, because higher numbers can lead to redundant cross-core communication.