tests.test_simple_shear_3d
PyDRex: Simple shear 3D tests.
1"""> PyDRex: Simple shear 3D tests.""" 2 3import contextlib as cl 4import functools as ft 5from time import process_time 6 7import numpy as np 8import pytest 9 10from pydrex import core as _core 11from pydrex import diagnostics as _diagnostics 12from pydrex import geometry as _geo 13from pydrex import io as _io 14from pydrex import logger as _log 15from pydrex import minerals as _minerals 16from pydrex import stats as _stats 17from pydrex import utils as _utils 18from pydrex import velocity as _velocity 19 20Pool, HAS_RAY = _utils.import_proc_pool() 21if HAS_RAY: 22 import ray # noqa: F401 23 24 from pydrex import distributed as _dstr # noqa: F401 25 26# Subdirectory of `outdir` used to store outputs from these tests. 27SUBDIR = "3d_simple_shear" 28 29 30class TestFraters2021: 31 """Tests inspired by the benchmarks presented in [Fraters & Billen, 2021]. 32 33 [Fraters & Billen, 2021]: https://doi.org/10.1029/2021GC009846 34 35 """ 36 37 class_id = "Fraters2021" 38 39 @classmethod 40 def run( 41 cls, 42 params, 43 timestamps, 44 get_velocity_gradient_initial, 45 get_velocity_gradient_final, 46 switch_time, 47 msg, 48 seed=None, 49 ): 50 """Run simulation with stationary particles in the given velocity gradient. 51 52 The optional RNG `seed` is used for the initial pseudorandom orientations. 53 A prefix `msg` will be printed before each timestep log message if given. 54 Other keyword args are passed to `pydrex.Mineral.update_orientations`. 55 56 Returns a tuple containing one olivine (A-type) and one enstatite mineral. 57 If `params["enstatite_fraction"]` is zero, then the second tuple element will be 58 `None` instead. 59 60 """ 61 62 def get_position(t): 63 return np.full(3, np.nan) 64 65 olivine = _minerals.Mineral( 66 phase=_core.MineralPhase.olivine, 67 fabric=_core.MineralFabric.olivine_A, 68 regime=_core.DeformationRegime.matrix_dislocation, 69 n_grains=params["number_of_grains"], 70 seed=seed, 71 ) 72 if params["enstatite_fraction"] > 0: 73 enstatite = _minerals.Mineral( 74 phase=_core.MineralPhase.enstatite, 75 fabric=_core.MineralFabric.enstatite_AB, 76 regime=_core.DeformationRegime.matrix_dislocation, 77 n_grains=params["number_of_grains"], 78 seed=seed, 79 ) 80 else: 81 enstatite = None 82 deformation_gradient = np.eye(3) # Undeformed initial state. 83 84 for t, time in enumerate(timestamps[:-1], start=1): 85 _log.info( 86 "%s; # %d; step %d/%d (t = %s)", msg, seed, t, len(timestamps) - 1, time 87 ) 88 if time > switch_time: 89 get_velocity_gradient = get_velocity_gradient_final 90 else: 91 get_velocity_gradient = get_velocity_gradient_initial 92 93 if params["enstatite_fraction"] > 0: 94 enstatite.update_orientations( 95 params, 96 deformation_gradient, 97 get_velocity_gradient, 98 pathline=(time, timestamps[t], get_position), 99 ) 100 deformation_gradient = olivine.update_orientations( 101 params, 102 deformation_gradient, 103 get_velocity_gradient, 104 pathline=(time, timestamps[t], get_position), 105 ) 106 _log.debug( 107 "› velocity gradient = %s", 108 get_velocity_gradient(np.nan, np.full(3, np.nan)).flatten(), 109 ) 110 return olivine, enstatite 111 112 @pytest.mark.slow 113 @pytest.mark.parametrize("switch_time_Ma", [0, 1, 2.5, np.inf]) 114 def test_direction_change( 115 self, outdir, seeds, params_Fraters2021, switch_time_Ma, ncpus, ray_session 116 ): 117 """Test a-axis alignment in simple shear with instantaneous geometry change. 118 119 The simulation runs for 5 Ma with a strain rate of 1.58e-14/s, resulting in an 120 accumulated strain invariant of 2.5. 121 122 The initial shear has nonzero du/dz and the final shear has nonzero dv/dx where 123 u is the velocity along X and v the velocity along Y. 124 125 """ 126 # Strain rate in units of strain per year, avoids tiny numbers ∼1e-14. 127 strain_rate = 5e-7 128 # With 500 steps, each step is ~3e11 seconds which is about as small as 129 # geodynamic modelling could feasibly go in most cases. 130 n_timestamps = 500 131 # Solve until D₀t=2.5 ('shear' γ=5). 132 timestamps = np.linspace(0, 5e6, n_timestamps) 133 _seeds = seeds[:500] # 500 runs as per Fraters & Billen, 2021, Figure 5. 134 n_seeds = len(_seeds) 135 136 _, get_velocity_gradient_initial = _velocity.simple_shear_2d( 137 "X", "Z", strain_rate 138 ) 139 _, get_velocity_gradient_final = _velocity.simple_shear_2d( 140 "Y", "X", strain_rate 141 ) 142 143 # Output arrays for mean [100] angles. 144 olA_from_proj_XZ = np.empty((n_seeds, n_timestamps)) 145 olA_from_proj_YX = np.empty_like(olA_from_proj_XZ) 146 ens_from_proj_XZ = np.empty_like(olA_from_proj_XZ) 147 ens_from_proj_YX = np.empty_like(olA_from_proj_XZ) 148 # Output arrays for M-index (CPO strength). 149 olA_strength = np.empty_like(olA_from_proj_XZ) 150 ens_strength = np.empty_like(olA_from_proj_XZ) 151 152 _id = _io.stringify(switch_time_Ma * 1e6) 153 optional_logging = cl.nullcontext() 154 if outdir is not None: 155 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_direction_change_{_id}" 156 optional_logging = _io.logfile_enable(f"{out_basepath}.log") 157 158 with optional_logging: 159 clock_start = process_time() 160 _run = ft.partial( 161 self.run, 162 params_Fraters2021, 163 timestamps, 164 get_velocity_gradient_initial, 165 get_velocity_gradient_final, 166 switch_time_Ma * 1e6, 167 _id, 168 ) 169 with Pool(processes=ncpus) as pool: 170 for s, out in enumerate(pool.imap_unordered(_run, _seeds)): 171 olivine, enstatite = out 172 _log.info("%s; # %d; postprocessing olivine...", _id, _seeds[s]) 173 olA_resampled, _ = _stats.resample_orientations( 174 olivine.orientations, olivine.fractions, seed=_seeds[s] 175 ) 176 olA_mean_vectors = np.array( 177 [ 178 _diagnostics.bingham_average(dcm, axis="a") 179 for dcm in olA_resampled 180 ] 181 ) 182 olA_from_proj_XZ[s, :] = np.array( 183 [ 184 _diagnostics.smallest_angle(v, v - v * [0, 1, 0]) 185 for v in olA_mean_vectors 186 ] 187 ) 188 olA_from_proj_YX[s, :] = np.array( 189 [ 190 _diagnostics.smallest_angle(v, v - v * [0, 0, 1]) 191 for v in olA_mean_vectors 192 ] 193 ) 194 olA_downsampled, _ = _stats.resample_orientations( 195 olivine.orientations, 196 olivine.fractions, 197 seed=_seeds[s], 198 n_samples=1000, 199 ) 200 olA_strength[s, :] = _diagnostics.misorientation_indices( 201 olA_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool 202 ) 203 204 del olivine, olA_resampled, olA_mean_vectors 205 206 _log.info("%s; # %d; postprocessing enstatite...", _id, _seeds[s]) 207 ens_resampled, _ = _stats.resample_orientations( 208 enstatite.orientations, enstatite.fractions, seed=_seeds[s] 209 ) 210 ens_mean_vectors = np.array( 211 [ 212 _diagnostics.bingham_average(dcm, axis="a") 213 for dcm in ens_resampled 214 ] 215 ) 216 ens_from_proj_XZ[s, :] = np.array( 217 [ 218 _diagnostics.smallest_angle(v, v - v * [0, 1, 0]) 219 for v in ens_mean_vectors 220 ] 221 ) 222 ens_from_proj_YX[s, :] = np.array( 223 [ 224 _diagnostics.smallest_angle(v, v - v * [0, 0, 1]) 225 for v in ens_mean_vectors 226 ] 227 ) 228 ens_downsampled, _ = _stats.resample_orientations( 229 enstatite.orientations, 230 enstatite.fractions, 231 seed=_seeds[s], 232 n_samples=1000, 233 ) 234 ens_strength[s, :] = _diagnostics.misorientation_indices( 235 ens_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool 236 ) 237 del enstatite, ens_resampled, ens_mean_vectors 238 239 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start)) 240 _log.info("calculating ensemble averages...") 241 olA_from_proj_XZ_mean = olA_from_proj_XZ.mean(axis=0) 242 olA_from_proj_XZ_err = olA_from_proj_XZ.std(axis=0) 243 244 olA_from_proj_YX_mean = olA_from_proj_YX.mean(axis=0) 245 olA_from_proj_YX_err = olA_from_proj_YX.std(axis=0) 246 247 ens_from_proj_XZ_mean = ens_from_proj_XZ.mean(axis=0) 248 ens_from_proj_XZ_err = ens_from_proj_XZ.std(axis=0) 249 250 ens_from_proj_YX_mean = ens_from_proj_YX.mean(axis=0) 251 ens_from_proj_YX_err = ens_from_proj_YX.std(axis=0) 252 253 if outdir is not None: 254 np.savez( 255 f"{out_basepath}.npz", 256 olA_from_proj_XZ_mean=olA_from_proj_XZ_mean, 257 olA_from_proj_XZ_err=olA_from_proj_XZ_err, 258 olA_from_proj_YX_mean=olA_from_proj_YX_mean, 259 olA_from_proj_YX_err=olA_from_proj_YX_err, 260 ens_from_proj_XZ_mean=ens_from_proj_XZ_mean, 261 ens_from_proj_XZ_err=ens_from_proj_XZ_err, 262 ens_from_proj_YX_mean=ens_from_proj_YX_mean, 263 ens_from_proj_YX_err=ens_from_proj_YX_err, 264 ) 265 np.savez( 266 f"{out_basepath}_strength.npz", 267 olA_mean=olA_strength.mean(axis=0), 268 ens_mean=ens_strength.mean(axis=0), 269 olA_err=olA_strength.std(axis=0), 270 ens_err=ens_strength.std(axis=0), 271 )
SUBDIR =
'3d_simple_shear'
class
TestFraters2021:
31class TestFraters2021: 32 """Tests inspired by the benchmarks presented in [Fraters & Billen, 2021]. 33 34 [Fraters & Billen, 2021]: https://doi.org/10.1029/2021GC009846 35 36 """ 37 38 class_id = "Fraters2021" 39 40 @classmethod 41 def run( 42 cls, 43 params, 44 timestamps, 45 get_velocity_gradient_initial, 46 get_velocity_gradient_final, 47 switch_time, 48 msg, 49 seed=None, 50 ): 51 """Run simulation with stationary particles in the given velocity gradient. 52 53 The optional RNG `seed` is used for the initial pseudorandom orientations. 54 A prefix `msg` will be printed before each timestep log message if given. 55 Other keyword args are passed to `pydrex.Mineral.update_orientations`. 56 57 Returns a tuple containing one olivine (A-type) and one enstatite mineral. 58 If `params["enstatite_fraction"]` is zero, then the second tuple element will be 59 `None` instead. 60 61 """ 62 63 def get_position(t): 64 return np.full(3, np.nan) 65 66 olivine = _minerals.Mineral( 67 phase=_core.MineralPhase.olivine, 68 fabric=_core.MineralFabric.olivine_A, 69 regime=_core.DeformationRegime.matrix_dislocation, 70 n_grains=params["number_of_grains"], 71 seed=seed, 72 ) 73 if params["enstatite_fraction"] > 0: 74 enstatite = _minerals.Mineral( 75 phase=_core.MineralPhase.enstatite, 76 fabric=_core.MineralFabric.enstatite_AB, 77 regime=_core.DeformationRegime.matrix_dislocation, 78 n_grains=params["number_of_grains"], 79 seed=seed, 80 ) 81 else: 82 enstatite = None 83 deformation_gradient = np.eye(3) # Undeformed initial state. 84 85 for t, time in enumerate(timestamps[:-1], start=1): 86 _log.info( 87 "%s; # %d; step %d/%d (t = %s)", msg, seed, t, len(timestamps) - 1, time 88 ) 89 if time > switch_time: 90 get_velocity_gradient = get_velocity_gradient_final 91 else: 92 get_velocity_gradient = get_velocity_gradient_initial 93 94 if params["enstatite_fraction"] > 0: 95 enstatite.update_orientations( 96 params, 97 deformation_gradient, 98 get_velocity_gradient, 99 pathline=(time, timestamps[t], get_position), 100 ) 101 deformation_gradient = olivine.update_orientations( 102 params, 103 deformation_gradient, 104 get_velocity_gradient, 105 pathline=(time, timestamps[t], get_position), 106 ) 107 _log.debug( 108 "› velocity gradient = %s", 109 get_velocity_gradient(np.nan, np.full(3, np.nan)).flatten(), 110 ) 111 return olivine, enstatite 112 113 @pytest.mark.slow 114 @pytest.mark.parametrize("switch_time_Ma", [0, 1, 2.5, np.inf]) 115 def test_direction_change( 116 self, outdir, seeds, params_Fraters2021, switch_time_Ma, ncpus, ray_session 117 ): 118 """Test a-axis alignment in simple shear with instantaneous geometry change. 119 120 The simulation runs for 5 Ma with a strain rate of 1.58e-14/s, resulting in an 121 accumulated strain invariant of 2.5. 122 123 The initial shear has nonzero du/dz and the final shear has nonzero dv/dx where 124 u is the velocity along X and v the velocity along Y. 125 126 """ 127 # Strain rate in units of strain per year, avoids tiny numbers ∼1e-14. 128 strain_rate = 5e-7 129 # With 500 steps, each step is ~3e11 seconds which is about as small as 130 # geodynamic modelling could feasibly go in most cases. 131 n_timestamps = 500 132 # Solve until D₀t=2.5 ('shear' γ=5). 133 timestamps = np.linspace(0, 5e6, n_timestamps) 134 _seeds = seeds[:500] # 500 runs as per Fraters & Billen, 2021, Figure 5. 135 n_seeds = len(_seeds) 136 137 _, get_velocity_gradient_initial = _velocity.simple_shear_2d( 138 "X", "Z", strain_rate 139 ) 140 _, get_velocity_gradient_final = _velocity.simple_shear_2d( 141 "Y", "X", strain_rate 142 ) 143 144 # Output arrays for mean [100] angles. 145 olA_from_proj_XZ = np.empty((n_seeds, n_timestamps)) 146 olA_from_proj_YX = np.empty_like(olA_from_proj_XZ) 147 ens_from_proj_XZ = np.empty_like(olA_from_proj_XZ) 148 ens_from_proj_YX = np.empty_like(olA_from_proj_XZ) 149 # Output arrays for M-index (CPO strength). 150 olA_strength = np.empty_like(olA_from_proj_XZ) 151 ens_strength = np.empty_like(olA_from_proj_XZ) 152 153 _id = _io.stringify(switch_time_Ma * 1e6) 154 optional_logging = cl.nullcontext() 155 if outdir is not None: 156 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_direction_change_{_id}" 157 optional_logging = _io.logfile_enable(f"{out_basepath}.log") 158 159 with optional_logging: 160 clock_start = process_time() 161 _run = ft.partial( 162 self.run, 163 params_Fraters2021, 164 timestamps, 165 get_velocity_gradient_initial, 166 get_velocity_gradient_final, 167 switch_time_Ma * 1e6, 168 _id, 169 ) 170 with Pool(processes=ncpus) as pool: 171 for s, out in enumerate(pool.imap_unordered(_run, _seeds)): 172 olivine, enstatite = out 173 _log.info("%s; # %d; postprocessing olivine...", _id, _seeds[s]) 174 olA_resampled, _ = _stats.resample_orientations( 175 olivine.orientations, olivine.fractions, seed=_seeds[s] 176 ) 177 olA_mean_vectors = np.array( 178 [ 179 _diagnostics.bingham_average(dcm, axis="a") 180 for dcm in olA_resampled 181 ] 182 ) 183 olA_from_proj_XZ[s, :] = np.array( 184 [ 185 _diagnostics.smallest_angle(v, v - v * [0, 1, 0]) 186 for v in olA_mean_vectors 187 ] 188 ) 189 olA_from_proj_YX[s, :] = np.array( 190 [ 191 _diagnostics.smallest_angle(v, v - v * [0, 0, 1]) 192 for v in olA_mean_vectors 193 ] 194 ) 195 olA_downsampled, _ = _stats.resample_orientations( 196 olivine.orientations, 197 olivine.fractions, 198 seed=_seeds[s], 199 n_samples=1000, 200 ) 201 olA_strength[s, :] = _diagnostics.misorientation_indices( 202 olA_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool 203 ) 204 205 del olivine, olA_resampled, olA_mean_vectors 206 207 _log.info("%s; # %d; postprocessing enstatite...", _id, _seeds[s]) 208 ens_resampled, _ = _stats.resample_orientations( 209 enstatite.orientations, enstatite.fractions, seed=_seeds[s] 210 ) 211 ens_mean_vectors = np.array( 212 [ 213 _diagnostics.bingham_average(dcm, axis="a") 214 for dcm in ens_resampled 215 ] 216 ) 217 ens_from_proj_XZ[s, :] = np.array( 218 [ 219 _diagnostics.smallest_angle(v, v - v * [0, 1, 0]) 220 for v in ens_mean_vectors 221 ] 222 ) 223 ens_from_proj_YX[s, :] = np.array( 224 [ 225 _diagnostics.smallest_angle(v, v - v * [0, 0, 1]) 226 for v in ens_mean_vectors 227 ] 228 ) 229 ens_downsampled, _ = _stats.resample_orientations( 230 enstatite.orientations, 231 enstatite.fractions, 232 seed=_seeds[s], 233 n_samples=1000, 234 ) 235 ens_strength[s, :] = _diagnostics.misorientation_indices( 236 ens_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool 237 ) 238 del enstatite, ens_resampled, ens_mean_vectors 239 240 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start)) 241 _log.info("calculating ensemble averages...") 242 olA_from_proj_XZ_mean = olA_from_proj_XZ.mean(axis=0) 243 olA_from_proj_XZ_err = olA_from_proj_XZ.std(axis=0) 244 245 olA_from_proj_YX_mean = olA_from_proj_YX.mean(axis=0) 246 olA_from_proj_YX_err = olA_from_proj_YX.std(axis=0) 247 248 ens_from_proj_XZ_mean = ens_from_proj_XZ.mean(axis=0) 249 ens_from_proj_XZ_err = ens_from_proj_XZ.std(axis=0) 250 251 ens_from_proj_YX_mean = ens_from_proj_YX.mean(axis=0) 252 ens_from_proj_YX_err = ens_from_proj_YX.std(axis=0) 253 254 if outdir is not None: 255 np.savez( 256 f"{out_basepath}.npz", 257 olA_from_proj_XZ_mean=olA_from_proj_XZ_mean, 258 olA_from_proj_XZ_err=olA_from_proj_XZ_err, 259 olA_from_proj_YX_mean=olA_from_proj_YX_mean, 260 olA_from_proj_YX_err=olA_from_proj_YX_err, 261 ens_from_proj_XZ_mean=ens_from_proj_XZ_mean, 262 ens_from_proj_XZ_err=ens_from_proj_XZ_err, 263 ens_from_proj_YX_mean=ens_from_proj_YX_mean, 264 ens_from_proj_YX_err=ens_from_proj_YX_err, 265 ) 266 np.savez( 267 f"{out_basepath}_strength.npz", 268 olA_mean=olA_strength.mean(axis=0), 269 ens_mean=ens_strength.mean(axis=0), 270 olA_err=olA_strength.std(axis=0), 271 ens_err=ens_strength.std(axis=0), 272 )
Tests inspired by the benchmarks presented in [Fraters & Billen, 2021].
@classmethod
def
run( cls, params, timestamps, get_velocity_gradient_initial, get_velocity_gradient_final, switch_time, msg, seed=None):
40 @classmethod 41 def run( 42 cls, 43 params, 44 timestamps, 45 get_velocity_gradient_initial, 46 get_velocity_gradient_final, 47 switch_time, 48 msg, 49 seed=None, 50 ): 51 """Run simulation with stationary particles in the given velocity gradient. 52 53 The optional RNG `seed` is used for the initial pseudorandom orientations. 54 A prefix `msg` will be printed before each timestep log message if given. 55 Other keyword args are passed to `pydrex.Mineral.update_orientations`. 56 57 Returns a tuple containing one olivine (A-type) and one enstatite mineral. 58 If `params["enstatite_fraction"]` is zero, then the second tuple element will be 59 `None` instead. 60 61 """ 62 63 def get_position(t): 64 return np.full(3, np.nan) 65 66 olivine = _minerals.Mineral( 67 phase=_core.MineralPhase.olivine, 68 fabric=_core.MineralFabric.olivine_A, 69 regime=_core.DeformationRegime.matrix_dislocation, 70 n_grains=params["number_of_grains"], 71 seed=seed, 72 ) 73 if params["enstatite_fraction"] > 0: 74 enstatite = _minerals.Mineral( 75 phase=_core.MineralPhase.enstatite, 76 fabric=_core.MineralFabric.enstatite_AB, 77 regime=_core.DeformationRegime.matrix_dislocation, 78 n_grains=params["number_of_grains"], 79 seed=seed, 80 ) 81 else: 82 enstatite = None 83 deformation_gradient = np.eye(3) # Undeformed initial state. 84 85 for t, time in enumerate(timestamps[:-1], start=1): 86 _log.info( 87 "%s; # %d; step %d/%d (t = %s)", msg, seed, t, len(timestamps) - 1, time 88 ) 89 if time > switch_time: 90 get_velocity_gradient = get_velocity_gradient_final 91 else: 92 get_velocity_gradient = get_velocity_gradient_initial 93 94 if params["enstatite_fraction"] > 0: 95 enstatite.update_orientations( 96 params, 97 deformation_gradient, 98 get_velocity_gradient, 99 pathline=(time, timestamps[t], get_position), 100 ) 101 deformation_gradient = olivine.update_orientations( 102 params, 103 deformation_gradient, 104 get_velocity_gradient, 105 pathline=(time, timestamps[t], get_position), 106 ) 107 _log.debug( 108 "› velocity gradient = %s", 109 get_velocity_gradient(np.nan, np.full(3, np.nan)).flatten(), 110 ) 111 return olivine, enstatite
Run simulation with stationary particles in the given velocity gradient.
The optional RNG seed is used for the initial pseudorandom orientations.
A prefix msg will be printed before each timestep log message if given.
Other keyword args are passed to pydrex.Mineral.update_orientations.
Returns a tuple containing one olivine (A-type) and one enstatite mineral.
If params["enstatite_fraction"] is zero, then the second tuple element will be
None instead.
@pytest.mark.slow
@pytest.mark.parametrize('switch_time_Ma', [0, 1, 2.5, np.inf])
def
test_direction_change( self, outdir, seeds, params_Fraters2021, switch_time_Ma, ncpus, ray_session):
113 @pytest.mark.slow 114 @pytest.mark.parametrize("switch_time_Ma", [0, 1, 2.5, np.inf]) 115 def test_direction_change( 116 self, outdir, seeds, params_Fraters2021, switch_time_Ma, ncpus, ray_session 117 ): 118 """Test a-axis alignment in simple shear with instantaneous geometry change. 119 120 The simulation runs for 5 Ma with a strain rate of 1.58e-14/s, resulting in an 121 accumulated strain invariant of 2.5. 122 123 The initial shear has nonzero du/dz and the final shear has nonzero dv/dx where 124 u is the velocity along X and v the velocity along Y. 125 126 """ 127 # Strain rate in units of strain per year, avoids tiny numbers ∼1e-14. 128 strain_rate = 5e-7 129 # With 500 steps, each step is ~3e11 seconds which is about as small as 130 # geodynamic modelling could feasibly go in most cases. 131 n_timestamps = 500 132 # Solve until D₀t=2.5 ('shear' γ=5). 133 timestamps = np.linspace(0, 5e6, n_timestamps) 134 _seeds = seeds[:500] # 500 runs as per Fraters & Billen, 2021, Figure 5. 135 n_seeds = len(_seeds) 136 137 _, get_velocity_gradient_initial = _velocity.simple_shear_2d( 138 "X", "Z", strain_rate 139 ) 140 _, get_velocity_gradient_final = _velocity.simple_shear_2d( 141 "Y", "X", strain_rate 142 ) 143 144 # Output arrays for mean [100] angles. 145 olA_from_proj_XZ = np.empty((n_seeds, n_timestamps)) 146 olA_from_proj_YX = np.empty_like(olA_from_proj_XZ) 147 ens_from_proj_XZ = np.empty_like(olA_from_proj_XZ) 148 ens_from_proj_YX = np.empty_like(olA_from_proj_XZ) 149 # Output arrays for M-index (CPO strength). 150 olA_strength = np.empty_like(olA_from_proj_XZ) 151 ens_strength = np.empty_like(olA_from_proj_XZ) 152 153 _id = _io.stringify(switch_time_Ma * 1e6) 154 optional_logging = cl.nullcontext() 155 if outdir is not None: 156 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_direction_change_{_id}" 157 optional_logging = _io.logfile_enable(f"{out_basepath}.log") 158 159 with optional_logging: 160 clock_start = process_time() 161 _run = ft.partial( 162 self.run, 163 params_Fraters2021, 164 timestamps, 165 get_velocity_gradient_initial, 166 get_velocity_gradient_final, 167 switch_time_Ma * 1e6, 168 _id, 169 ) 170 with Pool(processes=ncpus) as pool: 171 for s, out in enumerate(pool.imap_unordered(_run, _seeds)): 172 olivine, enstatite = out 173 _log.info("%s; # %d; postprocessing olivine...", _id, _seeds[s]) 174 olA_resampled, _ = _stats.resample_orientations( 175 olivine.orientations, olivine.fractions, seed=_seeds[s] 176 ) 177 olA_mean_vectors = np.array( 178 [ 179 _diagnostics.bingham_average(dcm, axis="a") 180 for dcm in olA_resampled 181 ] 182 ) 183 olA_from_proj_XZ[s, :] = np.array( 184 [ 185 _diagnostics.smallest_angle(v, v - v * [0, 1, 0]) 186 for v in olA_mean_vectors 187 ] 188 ) 189 olA_from_proj_YX[s, :] = np.array( 190 [ 191 _diagnostics.smallest_angle(v, v - v * [0, 0, 1]) 192 for v in olA_mean_vectors 193 ] 194 ) 195 olA_downsampled, _ = _stats.resample_orientations( 196 olivine.orientations, 197 olivine.fractions, 198 seed=_seeds[s], 199 n_samples=1000, 200 ) 201 olA_strength[s, :] = _diagnostics.misorientation_indices( 202 olA_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool 203 ) 204 205 del olivine, olA_resampled, olA_mean_vectors 206 207 _log.info("%s; # %d; postprocessing enstatite...", _id, _seeds[s]) 208 ens_resampled, _ = _stats.resample_orientations( 209 enstatite.orientations, enstatite.fractions, seed=_seeds[s] 210 ) 211 ens_mean_vectors = np.array( 212 [ 213 _diagnostics.bingham_average(dcm, axis="a") 214 for dcm in ens_resampled 215 ] 216 ) 217 ens_from_proj_XZ[s, :] = np.array( 218 [ 219 _diagnostics.smallest_angle(v, v - v * [0, 1, 0]) 220 for v in ens_mean_vectors 221 ] 222 ) 223 ens_from_proj_YX[s, :] = np.array( 224 [ 225 _diagnostics.smallest_angle(v, v - v * [0, 0, 1]) 226 for v in ens_mean_vectors 227 ] 228 ) 229 ens_downsampled, _ = _stats.resample_orientations( 230 enstatite.orientations, 231 enstatite.fractions, 232 seed=_seeds[s], 233 n_samples=1000, 234 ) 235 ens_strength[s, :] = _diagnostics.misorientation_indices( 236 ens_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool 237 ) 238 del enstatite, ens_resampled, ens_mean_vectors 239 240 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start)) 241 _log.info("calculating ensemble averages...") 242 olA_from_proj_XZ_mean = olA_from_proj_XZ.mean(axis=0) 243 olA_from_proj_XZ_err = olA_from_proj_XZ.std(axis=0) 244 245 olA_from_proj_YX_mean = olA_from_proj_YX.mean(axis=0) 246 olA_from_proj_YX_err = olA_from_proj_YX.std(axis=0) 247 248 ens_from_proj_XZ_mean = ens_from_proj_XZ.mean(axis=0) 249 ens_from_proj_XZ_err = ens_from_proj_XZ.std(axis=0) 250 251 ens_from_proj_YX_mean = ens_from_proj_YX.mean(axis=0) 252 ens_from_proj_YX_err = ens_from_proj_YX.std(axis=0) 253 254 if outdir is not None: 255 np.savez( 256 f"{out_basepath}.npz", 257 olA_from_proj_XZ_mean=olA_from_proj_XZ_mean, 258 olA_from_proj_XZ_err=olA_from_proj_XZ_err, 259 olA_from_proj_YX_mean=olA_from_proj_YX_mean, 260 olA_from_proj_YX_err=olA_from_proj_YX_err, 261 ens_from_proj_XZ_mean=ens_from_proj_XZ_mean, 262 ens_from_proj_XZ_err=ens_from_proj_XZ_err, 263 ens_from_proj_YX_mean=ens_from_proj_YX_mean, 264 ens_from_proj_YX_err=ens_from_proj_YX_err, 265 ) 266 np.savez( 267 f"{out_basepath}_strength.npz", 268 olA_mean=olA_strength.mean(axis=0), 269 ens_mean=ens_strength.mean(axis=0), 270 olA_err=olA_strength.std(axis=0), 271 ens_err=ens_strength.std(axis=0), 272 )
Test a-axis alignment in simple shear with instantaneous geometry change.
The simulation runs for 5 Ma with a strain rate of 1.58e-14/s, resulting in an accumulated strain invariant of 2.5.
The initial shear has nonzero du/dz and the final shear has nonzero dv/dx where u is the velocity along X and v the velocity along Y.