- 33class TestPreliminaries:
- 34 """Preliminary tests to check that various auxiliary routines are working."""
- 35
- 36 def test_strain_increment(self):
- 37 """Test for accumulating strain via strain increment calculations."""
- 38 _, get_velocity_gradient = _velocity.simple_shear_2d("X", "Z", 1)
- 39 timestamps = np.linspace(0, 1, 10) # Solve until D₀t=1 (tensorial strain).
- 40 strains_inc = np.zeros_like(timestamps)
- 41 L = get_velocity_gradient(np.nan, Ŋ([0e0, 0e0, 0e0]))
- 42 for i, ε in enumerate(strains_inc[1:]):
- 43 strains_inc[i + 1] = strains_inc[i] + _utils.strain_increment(
- 44 timestamps[1] - timestamps[0],
- 45 L,
- 46 )
- 47 # For constant timesteps, check strains == positive_timestamps * strain_rate.
- 48 nt.assert_allclose(strains_inc, timestamps, atol=6e-16, rtol=0)
- 49
- 50 # Same thing, but for strain rate similar to experiments.
- 51 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", 1e-5)
- 52 timestamps = np.linspace(0, 1e6, 10) # Solve until D₀t=10 (tensorial strain).
- 53 strains_inc = np.zeros_like(timestamps)
- 54 L = get_velocity_gradient(np.nan, Ŋ([0e0, 0e0, 0e0]))
- 55 for i, ε in enumerate(strains_inc[1:]):
- 56 strains_inc[i + 1] = strains_inc[i] + _utils.strain_increment(
- 57 timestamps[1] - timestamps[0],
- 58 L,
- 59 )
- 60 nt.assert_allclose(strains_inc, timestamps * 1e-5, atol=5e-15, rtol=0)
- 61
- 62 # Again, but this time the particle will move (using get_pathline).
- 63 # We use a 400km x 400km box and a strain rate of 1e-15 s⁻¹.
- 64 get_velocity, get_velocity_gradient = _velocity.simple_shear_2d("X", "Z", 1e-15)
- 65 timestamps, get_position = _paths.get_pathline(
- 66 Ŋ([1e5, 0e0, 1e5]),
- 67 get_velocity,
- 68 get_velocity_gradient,
- 69 Ŋ([-2e5, 0e0, -2e5]),
- 70 Ŋ([2e5, 0e0, 2e5]),
- 71 2,
- 72 regular_steps=10,
- 73 )
- 74 positions = [get_position(t) for t in timestamps]
- 75 velocity_gradients = [get_velocity_gradient(np.nan, Ŋ(x)) for x in positions]
- 76
- 77 # Check that polycrystal is experiencing steady velocity gradient.
- 78 nt.assert_array_equal(
- 79 velocity_gradients, np.full_like(velocity_gradients, velocity_gradients[0])
- 80 )
- 81 # Check that positions are changing as expected.
- 82 xdiff = np.diff(Ŋ([x[0] for x in positions]))
- 83 zdiff = np.diff(Ŋ([x[2] for x in positions]))
- 84 assert xdiff[0] > 0
- 85 assert zdiff[0] == 0
- 86 nt.assert_allclose(xdiff, np.full_like(xdiff, xdiff[0]), rtol=0, atol=1e-10)
- 87 nt.assert_allclose(zdiff, np.full_like(zdiff, zdiff[0]), rtol=0, atol=1e-10)
- 88 strains_inc = np.zeros_like(timestamps)
- 89 for t, time in enumerate(timestamps[:-1], start=1):
- 90 strains_inc[t] = strains_inc[t - 1] + (
- 91 _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
- 92 )
- 93 # fig, ax, _, _ = _vis.pathline_box2d(
- 94 # None,
- 95 # get_velocity,
- 96 # "xz",
- 97 # strains_inc,
- 98 # positions,
- 99 # ".",
-100 # Ŋ([-2e5, -2e5]),
-101 # Ŋ([2e5, 2e5]),
-102 # [20, 20],
-103 # )
-104 # fig.savefig("/tmp/fig.png")
-105 nt.assert_allclose(
-106 strains_inc,
-107 (timestamps - timestamps[0]) * 1e-15,
-108 atol=5e-15,
-109 rtol=0,
-110 )
+ 32class TestPreliminaries:
+ 33 """Preliminary tests to check that various auxiliary routines are working."""
+ 34
+ 35 def test_strain_increment(self):
+ 36 """Test for accumulating strain via strain increment calculations."""
+ 37 _, get_velocity_gradient = _velocity.simple_shear_2d("X", "Z", 1)
+ 38 timestamps = np.linspace(0, 1, 10) # Solve until D₀t=1 (tensorial strain).
+ 39 strains_inc = np.zeros_like(timestamps)
+ 40 L = get_velocity_gradient(np.nan, Ŋ([0e0, 0e0, 0e0]))
+ 41 for i, ε in enumerate(strains_inc[1:]):
+ 42 strains_inc[i + 1] = strains_inc[i] + _utils.strain_increment(
+ 43 timestamps[1] - timestamps[0],
+ 44 L,
+ 45 )
+ 46 # For constant timesteps, check strains == positive_timestamps * strain_rate.
+ 47 nt.assert_allclose(strains_inc, timestamps, atol=6e-16, rtol=0)
+ 48
+ 49 # Same thing, but for strain rate similar to experiments.
+ 50 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", 1e-5)
+ 51 timestamps = np.linspace(0, 1e6, 10) # Solve until D₀t=10 (tensorial strain).
+ 52 strains_inc = np.zeros_like(timestamps)
+ 53 L = get_velocity_gradient(np.nan, Ŋ([0e0, 0e0, 0e0]))
+ 54 for i, ε in enumerate(strains_inc[1:]):
+ 55 strains_inc[i + 1] = strains_inc[i] + _utils.strain_increment(
+ 56 timestamps[1] - timestamps[0],
+ 57 L,
+ 58 )
+ 59 nt.assert_allclose(strains_inc, timestamps * 1e-5, atol=5e-15, rtol=0)
+ 60
+ 61 # Again, but this time the particle will move (using get_pathline).
+ 62 # We use a 400km x 400km box and a strain rate of 1e-15 s⁻¹.
+ 63 get_velocity, get_velocity_gradient = _velocity.simple_shear_2d("X", "Z", 1e-15)
+ 64 timestamps, get_position = _paths.get_pathline(
+ 65 Ŋ([1e5, 0e0, 1e5]),
+ 66 get_velocity,
+ 67 get_velocity_gradient,
+ 68 Ŋ([-2e5, 0e0, -2e5]),
+ 69 Ŋ([2e5, 0e0, 2e5]),
+ 70 2,
+ 71 regular_steps=10,
+ 72 )
+ 73 positions = [get_position(t) for t in timestamps]
+ 74 velocity_gradients = [get_velocity_gradient(np.nan, Ŋ(x)) for x in positions]
+ 75
+ 76 # Check that polycrystal is experiencing steady velocity gradient.
+ 77 nt.assert_array_equal(
+ 78 velocity_gradients, np.full_like(velocity_gradients, velocity_gradients[0])
+ 79 )
+ 80 # Check that positions are changing as expected.
+ 81 xdiff = np.diff(Ŋ([x[0] for x in positions]))
+ 82 zdiff = np.diff(Ŋ([x[2] for x in positions]))
+ 83 assert xdiff[0] > 0
+ 84 assert zdiff[0] == 0
+ 85 nt.assert_allclose(xdiff, np.full_like(xdiff, xdiff[0]), rtol=0, atol=1e-10)
+ 86 nt.assert_allclose(zdiff, np.full_like(zdiff, zdiff[0]), rtol=0, atol=1e-10)
+ 87 strains_inc = np.zeros_like(timestamps)
+ 88 for t, time in enumerate(timestamps[:-1], start=1):
+ 89 strains_inc[t] = strains_inc[t - 1] + (
+ 90 _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
+ 91 )
+ 92 # fig, ax, _, _ = _vis.pathline_box2d(
+ 93 # None,
+ 94 # get_velocity,
+ 95 # "xz",
+ 96 # strains_inc,
+ 97 # positions,
+ 98 # ".",
+ 99 # Ŋ([-2e5, -2e5]),
+100 # Ŋ([2e5, 2e5]),
+101 # [20, 20],
+102 # )
+103 # fig.savefig("/tmp/fig.png")
+104 nt.assert_allclose(
+105 strains_inc,
+106 (timestamps - timestamps[0]) * 1e-15,
+107 atol=5e-15,
+108 rtol=0,
+109 )
36 def test_strain_increment(self):
- 37 """Test for accumulating strain via strain increment calculations."""
- 38 _, get_velocity_gradient = _velocity.simple_shear_2d("X", "Z", 1)
- 39 timestamps = np.linspace(0, 1, 10) # Solve until D₀t=1 (tensorial strain).
- 40 strains_inc = np.zeros_like(timestamps)
- 41 L = get_velocity_gradient(np.nan, Ŋ([0e0, 0e0, 0e0]))
- 42 for i, ε in enumerate(strains_inc[1:]):
- 43 strains_inc[i + 1] = strains_inc[i] + _utils.strain_increment(
- 44 timestamps[1] - timestamps[0],
- 45 L,
- 46 )
- 47 # For constant timesteps, check strains == positive_timestamps * strain_rate.
- 48 nt.assert_allclose(strains_inc, timestamps, atol=6e-16, rtol=0)
- 49
- 50 # Same thing, but for strain rate similar to experiments.
- 51 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", 1e-5)
- 52 timestamps = np.linspace(0, 1e6, 10) # Solve until D₀t=10 (tensorial strain).
- 53 strains_inc = np.zeros_like(timestamps)
- 54 L = get_velocity_gradient(np.nan, Ŋ([0e0, 0e0, 0e0]))
- 55 for i, ε in enumerate(strains_inc[1:]):
- 56 strains_inc[i + 1] = strains_inc[i] + _utils.strain_increment(
- 57 timestamps[1] - timestamps[0],
- 58 L,
- 59 )
- 60 nt.assert_allclose(strains_inc, timestamps * 1e-5, atol=5e-15, rtol=0)
- 61
- 62 # Again, but this time the particle will move (using get_pathline).
- 63 # We use a 400km x 400km box and a strain rate of 1e-15 s⁻¹.
- 64 get_velocity, get_velocity_gradient = _velocity.simple_shear_2d("X", "Z", 1e-15)
- 65 timestamps, get_position = _paths.get_pathline(
- 66 Ŋ([1e5, 0e0, 1e5]),
- 67 get_velocity,
- 68 get_velocity_gradient,
- 69 Ŋ([-2e5, 0e0, -2e5]),
- 70 Ŋ([2e5, 0e0, 2e5]),
- 71 2,
- 72 regular_steps=10,
- 73 )
- 74 positions = [get_position(t) for t in timestamps]
- 75 velocity_gradients = [get_velocity_gradient(np.nan, Ŋ(x)) for x in positions]
- 76
- 77 # Check that polycrystal is experiencing steady velocity gradient.
- 78 nt.assert_array_equal(
- 79 velocity_gradients, np.full_like(velocity_gradients, velocity_gradients[0])
- 80 )
- 81 # Check that positions are changing as expected.
- 82 xdiff = np.diff(Ŋ([x[0] for x in positions]))
- 83 zdiff = np.diff(Ŋ([x[2] for x in positions]))
- 84 assert xdiff[0] > 0
- 85 assert zdiff[0] == 0
- 86 nt.assert_allclose(xdiff, np.full_like(xdiff, xdiff[0]), rtol=0, atol=1e-10)
- 87 nt.assert_allclose(zdiff, np.full_like(zdiff, zdiff[0]), rtol=0, atol=1e-10)
- 88 strains_inc = np.zeros_like(timestamps)
- 89 for t, time in enumerate(timestamps[:-1], start=1):
- 90 strains_inc[t] = strains_inc[t - 1] + (
- 91 _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
- 92 )
- 93 # fig, ax, _, _ = _vis.pathline_box2d(
- 94 # None,
- 95 # get_velocity,
- 96 # "xz",
- 97 # strains_inc,
- 98 # positions,
- 99 # ".",
-100 # Ŋ([-2e5, -2e5]),
-101 # Ŋ([2e5, 2e5]),
-102 # [20, 20],
-103 # )
-104 # fig.savefig("/tmp/fig.png")
-105 nt.assert_allclose(
-106 strains_inc,
-107 (timestamps - timestamps[0]) * 1e-15,
-108 atol=5e-15,
-109 rtol=0,
-110 )
+ 35 def test_strain_increment(self):
+ 36 """Test for accumulating strain via strain increment calculations."""
+ 37 _, get_velocity_gradient = _velocity.simple_shear_2d("X", "Z", 1)
+ 38 timestamps = np.linspace(0, 1, 10) # Solve until D₀t=1 (tensorial strain).
+ 39 strains_inc = np.zeros_like(timestamps)
+ 40 L = get_velocity_gradient(np.nan, Ŋ([0e0, 0e0, 0e0]))
+ 41 for i, ε in enumerate(strains_inc[1:]):
+ 42 strains_inc[i + 1] = strains_inc[i] + _utils.strain_increment(
+ 43 timestamps[1] - timestamps[0],
+ 44 L,
+ 45 )
+ 46 # For constant timesteps, check strains == positive_timestamps * strain_rate.
+ 47 nt.assert_allclose(strains_inc, timestamps, atol=6e-16, rtol=0)
+ 48
+ 49 # Same thing, but for strain rate similar to experiments.
+ 50 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", 1e-5)
+ 51 timestamps = np.linspace(0, 1e6, 10) # Solve until D₀t=10 (tensorial strain).
+ 52 strains_inc = np.zeros_like(timestamps)
+ 53 L = get_velocity_gradient(np.nan, Ŋ([0e0, 0e0, 0e0]))
+ 54 for i, ε in enumerate(strains_inc[1:]):
+ 55 strains_inc[i + 1] = strains_inc[i] + _utils.strain_increment(
+ 56 timestamps[1] - timestamps[0],
+ 57 L,
+ 58 )
+ 59 nt.assert_allclose(strains_inc, timestamps * 1e-5, atol=5e-15, rtol=0)
+ 60
+ 61 # Again, but this time the particle will move (using get_pathline).
+ 62 # We use a 400km x 400km box and a strain rate of 1e-15 s⁻¹.
+ 63 get_velocity, get_velocity_gradient = _velocity.simple_shear_2d("X", "Z", 1e-15)
+ 64 timestamps, get_position = _paths.get_pathline(
+ 65 Ŋ([1e5, 0e0, 1e5]),
+ 66 get_velocity,
+ 67 get_velocity_gradient,
+ 68 Ŋ([-2e5, 0e0, -2e5]),
+ 69 Ŋ([2e5, 0e0, 2e5]),
+ 70 2,
+ 71 regular_steps=10,
+ 72 )
+ 73 positions = [get_position(t) for t in timestamps]
+ 74 velocity_gradients = [get_velocity_gradient(np.nan, Ŋ(x)) for x in positions]
+ 75
+ 76 # Check that polycrystal is experiencing steady velocity gradient.
+ 77 nt.assert_array_equal(
+ 78 velocity_gradients, np.full_like(velocity_gradients, velocity_gradients[0])
+ 79 )
+ 80 # Check that positions are changing as expected.
+ 81 xdiff = np.diff(Ŋ([x[0] for x in positions]))
+ 82 zdiff = np.diff(Ŋ([x[2] for x in positions]))
+ 83 assert xdiff[0] > 0
+ 84 assert zdiff[0] == 0
+ 85 nt.assert_allclose(xdiff, np.full_like(xdiff, xdiff[0]), rtol=0, atol=1e-10)
+ 86 nt.assert_allclose(zdiff, np.full_like(zdiff, zdiff[0]), rtol=0, atol=1e-10)
+ 87 strains_inc = np.zeros_like(timestamps)
+ 88 for t, time in enumerate(timestamps[:-1], start=1):
+ 89 strains_inc[t] = strains_inc[t - 1] + (
+ 90 _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
+ 91 )
+ 92 # fig, ax, _, _ = _vis.pathline_box2d(
+ 93 # None,
+ 94 # get_velocity,
+ 95 # "xz",
+ 96 # strains_inc,
+ 97 # positions,
+ 98 # ".",
+ 99 # Ŋ([-2e5, -2e5]),
+100 # Ŋ([2e5, 2e5]),
+101 # [20, 20],
+102 # )
+103 # fig.savefig("/tmp/fig.png")
+104 nt.assert_allclose(
+105 strains_inc,
+106 (timestamps - timestamps[0]) * 1e-15,
+107 atol=5e-15,
+108 rtol=0,
+109 )
113class TestOlivineA:
-114 """Tests for stationary A-type olivine polycrystals in 2D simple shear."""
-115
-116 class_id = "olivineA"
-117
-118 @classmethod
-119 def run(
-120 cls,
-121 params,
-122 timestamps,
-123 strain_rate,
-124 get_velocity_gradient,
-125 shear_direction,
-126 seed=None,
-127 return_fse=None,
-128 get_position=lambda t: np.zeros(3), # Stationary particles by default.
-129 ):
-130 """Reusable logic for 2D olivine (A-type) simple shear tests.
-131
-132 Returns a tuple with the mineral and the FSE angle (or `None` if `return_fse` is
-133 `None`).
-134
-135 """
-136 mineral = _minerals.Mineral(
-137 phase=_core.MineralPhase.olivine,
-138 fabric=_core.MineralFabric.olivine_A,
-139 regime=_core.DeformationRegime.matrix_dislocation,
-140 n_grains=params["number_of_grains"],
-141 seed=seed,
-142 )
-143 deformation_gradient = np.eye(3) # Undeformed initial state.
-144 θ_fse = np.empty_like(timestamps)
-145 θ_fse[0] = 45
-146
-147 for t, time in enumerate(timestamps[:-1], start=1):
-148 # Set up logging message depending on dynamic parameter and seeds.
-149 msg_start = (
-150 f"N = {params['number_of_grains']}; "
-151 + f"λ∗ = {params['nucleation_efficiency']}; "
-152 + f"X = {params['gbs_threshold']}; "
-153 + f"M∗ = {params['gbm_mobility']}; "
-154 )
-155 if seed is not None:
-156 msg_start += f"# {seed}; "
-157
-158 _log.info(msg_start + "step %s/%s (t = %s)", t, len(timestamps) - 1, time)
-159
-160 deformation_gradient = mineral.update_orientations(
-161 params,
-162 deformation_gradient,
-163 get_velocity_gradient,
-164 pathline=(time, timestamps[t], get_position),
-165 )
-166 _log.debug(
-167 "› velocity gradient = %s",
-168 get_velocity_gradient(np.nan, np.full(3, np.nan)).flatten(),
-169 )
-170 _log.debug("› strain D₀t = %.2f", strain_rate * timestamps[t])
-171 _log.debug(
-172 "› grain fractions: median = %s, max = %s, min = %s",
-173 np.median(mineral.fractions[-1]),
-174 np.max(mineral.fractions[-1]),
-175 np.min(mineral.fractions[-1]),
-176 )
-177 if return_fse:
-178 _, fse_v = _diagnostics.finite_strain(deformation_gradient)
-179 θ_fse[t] = _diagnostics.smallest_angle(fse_v, shear_direction)
-180 else:
-181 θ_fse = None
-182
-183 return mineral, θ_fse
-184
-185 @classmethod
-186 def interp_GBM_Kaminski2001(cls, strains):
-187 """Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`."""
-188 _log.info("interpolating target CPO angles...")
-189 data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2001_GBMshear.scsv")
-190 cs_M0 = PchipInterpolator(
-191 _utils.remove_nans(data.equivalent_strain_M0) / 200,
-192 _utils.remove_nans(data.angle_M0),
-193 )
-194 cs_M50 = PchipInterpolator(
-195 _utils.remove_nans(data.equivalent_strain_M50) / 200,
-196 _utils.remove_nans(data.angle_M50),
-197 )
-198 cs_M200 = PchipInterpolator(
-199 _utils.remove_nans(data.equivalent_strain_M200) / 200,
-200 _utils.remove_nans(data.angle_M200),
-201 )
-202 return [cs_M0(strains), cs_M50(strains), cs_M200(strains)]
-203
-204 @classmethod
-205 def interp_GBM_FortranDRex(cls, strains):
-206 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
-207 _log.info("interpolating target CPO angles...")
-208 data = _io.read_scsv(_io.data("drexF90") / "olA_D1E4_dt50_X0_L5.scsv")
-209 data_strains = np.linspace(0, 1, 200)
-210 cs_M0 = PchipInterpolator(data_strains, _utils.remove_nans(data.M0_angle))
-211 cs_M50 = PchipInterpolator(data_strains, _utils.remove_nans(data.M50_angle))
-212 cs_M200 = PchipInterpolator(data_strains, _utils.remove_nans(data.M200_angle))
-213 return [cs_M0(strains), cs_M50(strains), cs_M200(strains)]
-214
-215 @classmethod
-216 def interp_GBS_FortranDRex(cls, strains):
-217 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
-218 _log.info("interpolating target CPO angles...")
-219 data = _io.read_scsv(_io.data("thirdparty") / "a_axis_GBS_fortran.scsv")
-220 data_strains = np.linspace(0, 1, 200)
-221 cs_X0 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X0))
-222 cs_X20 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X20))
-223 cs_X40 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X40))
-224 return [cs_X0(strains), cs_X20(strains), cs_X40(strains)]
-225
-226 @classmethod
-227 def interp_GBS_long_FortranDRex(cls, strains):
-228 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
-229 _log.info("interpolating target CPO angles...")
-230 data = _io.read_scsv(_io.data("thirdparty") / "a_axis_GBS_long_fortran.scsv")
-231 data_strains = np.linspace(0, 2.5, 500)
-232 cs_X0 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X0))
-233 cs_X20 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X20))
-234 cs_X40 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X40))
-235 return [cs_X0(strains), cs_X20(strains), cs_X40(strains)]
-236
-237 @classmethod
-238 def interp_GBS_Kaminski2004(cls, strains):
-239 """Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`."""
-240 _log.info("interpolating target CPO angles...")
-241 data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2004_GBSshear.scsv")
-242 cs_X0 = PchipInterpolator(
-243 _utils.remove_nans(data.dimensionless_time_X0),
-244 45 + _utils.remove_nans(data.angle_X0),
-245 )
-246 cs_X0d2 = PchipInterpolator(
-247 _utils.remove_nans(data.dimensionless_time_X0d2),
-248 45 + _utils.remove_nans(data.angle_X0d2),
-249 )
-250 cs_X0d4 = PchipInterpolator(
-251 _utils.remove_nans(data.dimensionless_time_X0d4),
-252 45 + _utils.remove_nans(data.angle_X0d4),
-253 )
-254 return [cs_X0(strains), cs_X0d2(strains), cs_X0d4(strains)]
-255
-256 @pytest.mark.skipif(_utils.in_ci("win32"), reason="Unable to allocate memory")
-257 def test_zero_recrystallisation(self, seed):
-258 """Check that M*=0 is a reliable switch to turn off recrystallisation."""
-259 params = _core.DefaultParams().as_dict()
-260 params["gbm_mobility"] = 0
-261 strain_rate = 1
-262 timestamps = np.linspace(0, 1, 25) # Solve until D₀t=1 (tensorial strain).
-263 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-264 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-265 mineral, _ = self.run(
-266 params,
-267 timestamps,
-268 strain_rate,
-269 get_velocity_gradient,
-270 shear_direction,
-271 seed=seed,
-272 )
-273 for fractions in mineral.fractions[1:]:
-274 nt.assert_allclose(fractions, mineral.fractions[0], atol=1e-15, rtol=0)
-275
-276 @pytest.mark.skipif(_utils.in_ci("win32"), reason="Unable to allocate memory")
-277 @pytest.mark.parametrize("gbm_mobility", [50, 100, 150])
-278 def test_grainsize_median(self, seed, gbm_mobility):
-279 """Check that M*={50,100,150}, λ*=5 causes decreasing grain size median."""
-280 params = _core.DefaultParams().as_dict()
-281 params["gbm_mobility"] = gbm_mobility
-282 params["nucleation_efficiency"] = 5
-283 strain_rate = 1
-284 timestamps = np.linspace(0, 1, 25) # Solve until D₀t=1 (tensorial strain).
-285 n_timestamps = len(timestamps)
-286 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-287 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-288 mineral, _ = self.run(
-289 params,
-290 timestamps,
-291 strain_rate,
-292 get_velocity_gradient,
-293 shear_direction,
-294 seed=seed,
-295 )
-296 medians = np.empty(n_timestamps)
-297 for i, fractions in enumerate(mineral.fractions):
-298 medians[i] = np.median(fractions)
-299
-300 # The first diff is positive (~1e-6) before nucleation sets in.
-301 nt.assert_array_less(np.diff(medians)[1:], np.full(n_timestamps - 2, 0))
-302
-303 @pytest.mark.slow
-304 @pytest.mark.parametrize("gbs_threshold", [0, 0.2, 0.4])
-305 @pytest.mark.parametrize("nucleation_efficiency", [3, 5, 10])
-306 def test_dvdx_ensemble(
-307 self, outdir, seeds_nearX45, ncpus, gbs_threshold, nucleation_efficiency
-308 ):
-309 r"""Test a-axis alignment to shear in Y direction (init. SCCS near 45° from X).
-310
-311 Velocity gradient:
-312 $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$
-313
-314 """
-315 strain_rate = 1
-316 timestamps = np.linspace(0, 1, 201) # Solve until D₀t=1 ('shear' γ=2).
-317 n_timestamps = len(timestamps)
-318 # Use `seeds` instead of `seeds_nearX45` if you have even more RAM and CPU time.
-319 _seeds = seeds_nearX45
-320 n_seeds = len(_seeds)
-321
-322 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-323 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-324
-325 gbm_mobilities = [0, 50, 125, 200]
-326 markers = ("x", "*", "d", "s")
-327
-328 _id = f"X{_io.stringify(gbs_threshold)}_L{_io.stringify(nucleation_efficiency)}"
-329 # Output setup with optional logging and data series labels.
-330 θ_fse = np.empty_like(timestamps)
-331 angles = np.empty((len(gbm_mobilities), n_seeds, n_timestamps))
-332 optional_logging = cl.nullcontext()
-333 if outdir is not None:
-334 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_ensemble_{_id}"
-335 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
-336 labels = []
-337
-338 with optional_logging:
-339 clock_start = process_time()
-340 for m, gbm_mobility in enumerate(gbm_mobilities):
-341 if m == 0:
-342 return_fse = True
-343 else:
-344 return_fse = False
-345
-346 params = {
-347 "phase_assemblage": (_core.MineralPhase.olivine,),
-348 "phase_fractions": (1.0,),
-349 "enstatite_fraction": 0.0,
-350 "stress_exponent": 1.5,
-351 "deformation_exponent": 3.5,
-352 "gbm_mobility": gbm_mobility,
-353 "gbs_threshold": gbs_threshold,
-354 "nucleation_efficiency": nucleation_efficiency,
-355 "number_of_grains": 5000,
-356 "initial_olivine_fabric": "A",
-357 }
-358
-359 _run = ft.partial(
-360 self.run,
-361 params,
-362 timestamps,
-363 strain_rate,
-364 get_velocity_gradient,
-365 shear_direction,
-366 return_fse=return_fse,
-367 )
-368 with Pool(processes=ncpus) as pool:
-369 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
-370 mineral, fse_angles = out
-371 angles[m, s, :] = [
-372 _diagnostics.smallest_angle(v, shear_direction)
-373 for v in _diagnostics.elasticity_components(
-374 _minerals.voigt_averages([mineral], params)
-375 )["hexagonal_axis"]
-376 ]
-377 # Save the whole mineral for the first seed only.
-378 if outdir is not None and s == 0:
-379 postfix = (
-380 f"M{_io.stringify(gbm_mobility)}"
-381 + f"_X{_io.stringify(gbs_threshold)}"
-382 + f"_L{_io.stringify(nucleation_efficiency)}"
-383 )
-384 mineral.save(f"{out_basepath}.npz", postfix=postfix)
-385 if return_fse:
-386 θ_fse += fse_angles
-387
-388 if return_fse:
-389 θ_fse /= n_seeds
-390
-391 if outdir is not None:
-392 labels.append(f"$M^∗$ = {gbm_mobility}")
-393
-394 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-395
-396 # Take ensemble means and optionally plot figure.
-397 strains = timestamps * strain_rate
-398 _log.info("postprocessing results for %s", _id)
-399 result_angles = angles.mean(axis=1)
-400 result_angles_err = angles.std(axis=1)
-401
-402 if outdir is not None:
-403 schema = {
-404 "delimiter": ",",
-405 "missing": "-",
-406 "fields": [
-407 {
-408 "name": "strain",
-409 "type": "integer",
-410 "unit": "percent",
-411 "fill": 999999,
-412 }
-413 ],
-414 }
-415 np.savez(
-416 f"{out_basepath}.npz",
-417 angles=result_angles,
-418 angles_err=result_angles_err,
-419 )
-420 _io.save_scsv(
-421 f"{out_basepath}_strains.scsv",
-422 schema,
-423 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
-424 )
-425 fig, ax, colors = _vis.alignment(
-426 None,
-427 strains,
-428 result_angles,
-429 markers,
-430 labels,
-431 err=result_angles_err,
-432 θ_max=60,
-433 θ_fse=θ_fse,
-434 )
-435 fig.savefig(_io.resolve_path(f"{out_basepath}.pdf"))
-436
-437 @pytest.mark.slow
-438 def test_dvdx_GBM(self, outdir, seeds_nearX45, ncpus):
-439 r"""Test a-axis alignment to shear in Y direction (init. SCCS near 45° from X).
-440
-441 Velocity gradient:
-442 $$
-443 \bm{L} = 10^{-4} ×
-444 \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}
-445 $$
-446
-447 Results are compared to the Fortran 90 output.
-448
-449 """
-450 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-451 strain_rate = 1e-4
-452 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-453 timestamps = np.linspace(0, 1e4, 51) # Solve until D₀t=1 ('shear' γ=2).
-454 i_strain_40p = 10 # Index of 40% strain, lower strains are not relevant here.
-455 i_strain_100p = 25 # Index of 100% strain, when M*=0 matches FSE.
-456 params = _core.DefaultParams().as_dict()
-457 params["gbs_threshold"] = 0 # No GBS, to match the Fortran parameters.
-458 gbm_mobilities = (0, 10, 50, 125, 200) # Must be in ascending order.
-459 markers = ("x", ".", "*", "d", "s")
-460 # Use `seeds` instead of `seeds_nearX45` if you have even more RAM and CPU time.
-461 _seeds = seeds_nearX45
-462 n_seeds = len(_seeds)
-463 angles = np.empty((len(gbm_mobilities), n_seeds, len(timestamps)))
-464 θ_fse = np.zeros_like(timestamps)
-465 strains = timestamps * strain_rate
-466 M0_drexF90, M50_drexF90, M200_drexF90 = self.interp_GBM_FortranDRex(strains)
-467
-468 optional_logging = cl.nullcontext()
-469 if outdir is not None:
-470 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_mobility"
-471 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
-472 labels = []
-473
-474 with optional_logging:
-475 clock_start = process_time()
-476 for m, gbm_mobility in enumerate(gbm_mobilities):
-477 if m == 0:
-478 return_fse = True
-479 else:
-480 return_fse = False
-481 params["gbm_mobility"] = gbm_mobility
-482
-483 _run = ft.partial(
-484 self.run,
-485 params,
-486 timestamps,
-487 strain_rate,
-488 get_velocity_gradient,
-489 shear_direction,
-490 return_fse=True,
-491 )
-492 with Pool(processes=ncpus) as pool:
-493 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
-494 mineral, fse_angles = out
-495 angles[m, s, :] = [
-496 _diagnostics.smallest_angle(v, shear_direction)
-497 for v in _diagnostics.elasticity_components(
-498 _minerals.voigt_averages([mineral], params)
-499 )["hexagonal_axis"]
-500 ]
-501 # Save the whole mineral for the first seed only.
-502 if outdir is not None and s == 0:
-503 mineral.save(
-504 f"{out_basepath}.npz",
-505 postfix=f"M{_io.stringify(gbm_mobility)}",
-506 )
-507 if return_fse:
-508 θ_fse += fse_angles
-509
-510 if return_fse:
-511 θ_fse /= n_seeds
-512
-513 if outdir is not None:
-514 labels.append(f"$M^∗$ = {params['gbm_mobility']}")
-515
-516 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-517
-518 # Take ensemble means and optionally plot figure.
-519 _log.info("postprocessing results...")
-520 result_angles = angles.mean(axis=1)
-521 result_angles_err = angles.std(axis=1)
-522
-523 if outdir is not None:
-524 schema = {
-525 "delimiter": ",",
-526 "missing": "-",
-527 "fields": [
-528 {
-529 "name": "strain",
-530 "type": "integer",
-531 "unit": "percent",
-532 "fill": 999999,
-533 }
-534 ],
-535 }
-536 _io.save_scsv(
-537 f"{out_basepath}_strains.scsv",
-538 schema,
-539 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
-540 )
-541 np.savez(
-542 f"{out_basepath}_angles.npz",
-543 angles=result_angles,
-544 err=result_angles_err,
-545 )
-546 fig, ax, colors = _vis.alignment(
-547 None,
-548 strains,
-549 result_angles,
-550 markers,
-551 labels,
-552 err=result_angles_err,
-553 θ_max=60,
-554 θ_fse=θ_fse,
-555 )
-556 ax.plot(strains, M0_drexF90, c=colors[0])
-557 ax.plot(strains, M50_drexF90, c=colors[2])
-558 ax.plot(strains, M200_drexF90, c=colors[4])
-559 _vis.show_Skemer2016_ShearStrainAngles(
-560 ax,
-561 ["Z&K 1200 C", "Z&K 1300 C"],
-562 ["v", "^"],
-563 ["k", "k"],
-564 ["none", None],
-565 [
-566 "Zhang & Karato, 1995\n(1473 K)",
-567 "Zhang & Karato, 1995\n(1573 K)",
-568 ],
-569 _core.MineralFabric.olivine_A,
-570 )
-571 # There is a lot of stuff on this legend, so put it outside the axes.
-572 # These values might need to be tweaked depending on the font size, etc.
-573 _legend = _utils.redraw_legend(ax, fig=fig, bbox_to_anchor=(1.66, 0.99))
-574 fig.savefig(
-575 _io.resolve_path(f"{out_basepath}.pdf"),
-576 bbox_extra_artists=(_legend,),
-577 bbox_inches="tight",
-578 )
-579
-580 # Check that GBM speeds up the alignment between 40% and 100% strain.
-581 _log.info("checking grain orientations...")
-582 for i, θ in enumerate(result_angles[:-1], start=1):
-583 nt.assert_array_less(
-584 result_angles[i][i_strain_40p:i_strain_100p],
-585 θ[i_strain_40p:i_strain_100p],
-586 )
-587
-588 # Check that M*=0 matches FSE (±1°) past 100% strain.
-589 nt.assert_allclose(
-590 result_angles[0][i_strain_100p:],
-591 θ_fse[i_strain_100p:],
-592 atol=1,
-593 rtol=0,
-594 )
-595
-596 # Check that results match Fortran output past 40% strain.
-597 nt.assert_allclose(
-598 result_angles[0][i_strain_40p:],
-599 M0_drexF90[i_strain_40p:],
-600 atol=0,
-601 rtol=0.1, # At 40% strain the match is worse than at higher strain.
-602 )
-603 nt.assert_allclose(
-604 result_angles[2][i_strain_40p:],
-605 M50_drexF90[i_strain_40p:],
-606 atol=1,
-607 rtol=0,
-608 )
-609 nt.assert_allclose(
-610 result_angles[4][i_strain_40p:],
-611 M200_drexF90[i_strain_40p:],
-612 atol=1.5,
-613 rtol=0,
-614 )
-615
-616 @pytest.mark.slow
-617 def test_GBM_calibration(self, outdir, seeds, ncpus):
-618 r"""Compare results for various values of $$M^∗$$ to A-type olivine data.
-619
-620 Velocity gradient:
-621 $$
-622 \bm{L} = 10^{-4} ×
-623 \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}
-624 $$
-625
-626 Unlike `test_dvdx_GBM`,
-627 grain boudary sliding is enabled here (see `_core.DefaultParams`).
-628 Data are provided by [Skemer & Hansen, 2016](http://dx.doi.org/10.1016/j.tecto.2015.12.003).
-629
-630 """
-631 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-632 strain_rate = 1
-633 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-634 timestamps = np.linspace(0, 3.2, 65) # Solve until D₀t=3.2 ('shear' γ=6.4).
-635 params = _core.DefaultParams().as_dict()
-636 params["number_of_grains"] = 5000
-637 gbm_mobilities = (0, 10, 50, 125) # Must be in ascending order.
-638 markers = ("x", "*", "1", ".")
-639 # Uses 100 seeds by default; use all 1000 if you have more RAM and CPU time.
-640 _seeds = seeds[:100]
-641 n_seeds = len(_seeds)
-642 angles = np.empty((len(gbm_mobilities), n_seeds, len(timestamps)))
-643 θ_fse = np.zeros_like(timestamps)
-644 strains = timestamps * strain_rate
-645
-646 optional_logging = cl.nullcontext()
-647 if outdir is not None:
-648 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_calibration"
-649 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
-650 labels = []
-651
-652 with optional_logging:
-653 clock_start = process_time()
-654 for m, gbm_mobility in enumerate(gbm_mobilities):
-655 return_fse = True if m == 0 else False
-656 params["gbm_mobility"] = gbm_mobility
-657 _run = ft.partial(
-658 self.run,
-659 params,
-660 timestamps,
-661 strain_rate,
-662 get_velocity_gradient,
-663 shear_direction,
-664 return_fse=return_fse,
-665 )
-666 with Pool(processes=ncpus) as pool:
-667 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
-668 mineral, fse_angles = out
-669 angles[m, s, :] = [
-670 _diagnostics.smallest_angle(v, shear_direction)
-671 for v in _diagnostics.elasticity_components(
-672 _minerals.voigt_averages([mineral], params)
-673 )["hexagonal_axis"]
-674 ]
-675 # Save the whole mineral for the first seed only.
-676 if outdir is not None and s == 0:
-677 mineral.save(
-678 f"{out_basepath}.npz",
-679 postfix=f"M{_io.stringify(gbm_mobility)}",
-680 )
-681 if return_fse:
-682 θ_fse += fse_angles
-683
-684 if return_fse:
-685 θ_fse /= n_seeds
-686
-687 if outdir is not None:
-688 labels.append(f"$M^∗$ = {params['gbm_mobility']}")
-689
-690 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-691
-692 # Take ensemble means and optionally plot figure.
-693 _log.info("postprocessing results...")
-694 result_angles = angles.mean(axis=1)
-695 result_angles_err = angles.std(axis=1)
-696
-697 if outdir is not None:
-698 schema = {
-699 "delimiter": ",",
-700 "missing": "-",
-701 "fields": [
-702 {
-703 "name": "strain",
-704 "type": "integer",
-705 "unit": "percent",
-706 "fill": 999999,
-707 }
-708 ],
-709 }
-710 _io.save_scsv(
-711 f"{out_basepath}_strains.scsv",
-712 schema,
-713 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
-714 )
-715 np.savez(
-716 _io.resolve_path(f"{out_basepath}_ensemble_means.npz"),
-717 angles=result_angles,
-718 err=result_angles_err,
-719 )
-720 fig = _vis.figure(
-721 figsize=(_vis.DEFAULT_FIG_WIDTH * 3, _vis.DEFAULT_FIG_HEIGHT)
-722 )
-723 fig, ax, colors = _vis.alignment(
-724 fig.add_subplot(),
-725 strains,
-726 result_angles,
-727 markers,
-728 labels,
-729 err=result_angles_err,
-730 θ_max=80,
-731 θ_fse=θ_fse,
-732 )
-733 (
+ 112class TestOlivineA:
+113 """Tests for stationary A-type olivine polycrystals in 2D simple shear."""
+114
+115 class_id = "olivineA"
+116
+117 @classmethod
+118 def run(
+119 cls,
+120 params,
+121 timestamps,
+122 strain_rate,
+123 get_velocity_gradient,
+124 shear_direction,
+125 seed=None,
+126 return_fse=None,
+127 get_position=lambda t: np.zeros(3), # Stationary particles by default.
+128 ):
+129 """Reusable logic for 2D olivine (A-type) simple shear tests.
+130
+131 Returns a tuple with the mineral and the FSE angle (or `None` if `return_fse` is
+132 `None`).
+133
+134 """
+135 mineral = _minerals.Mineral(
+136 phase=_core.MineralPhase.olivine,
+137 fabric=_core.MineralFabric.olivine_A,
+138 regime=_core.DeformationRegime.matrix_dislocation,
+139 n_grains=params["number_of_grains"],
+140 seed=seed,
+141 )
+142 deformation_gradient = np.eye(3) # Undeformed initial state.
+143 θ_fse = np.empty_like(timestamps)
+144 θ_fse[0] = 45
+145
+146 for t, time in enumerate(timestamps[:-1], start=1):
+147 # Set up logging message depending on dynamic parameter and seeds.
+148 msg_start = (
+149 f"N = {params['number_of_grains']}; "
+150 + f"λ∗ = {params['nucleation_efficiency']}; "
+151 + f"X = {params['gbs_threshold']}; "
+152 + f"M∗ = {params['gbm_mobility']}; "
+153 )
+154 if seed is not None:
+155 msg_start += f"# {seed}; "
+156
+157 _log.info(msg_start + "step %s/%s (t = %s)", t, len(timestamps) - 1, time)
+158
+159 deformation_gradient = mineral.update_orientations(
+160 params,
+161 deformation_gradient,
+162 get_velocity_gradient,
+163 pathline=(time, timestamps[t], get_position),
+164 )
+165 _log.debug(
+166 "› velocity gradient = %s",
+167 get_velocity_gradient(np.nan, np.full(3, np.nan)).flatten(),
+168 )
+169 _log.debug("› strain D₀t = %.2f", strain_rate * timestamps[t])
+170 _log.debug(
+171 "› grain fractions: median = %s, max = %s, min = %s",
+172 np.median(mineral.fractions[-1]),
+173 np.max(mineral.fractions[-1]),
+174 np.min(mineral.fractions[-1]),
+175 )
+176 if return_fse:
+177 _, fse_v = _diagnostics.finite_strain(deformation_gradient)
+178 θ_fse[t] = _diagnostics.smallest_angle(fse_v, shear_direction)
+179 else:
+180 θ_fse = None
+181
+182 return mineral, θ_fse
+183
+184 @classmethod
+185 def interp_GBM_Kaminski2001(cls, strains):
+186 """Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`."""
+187 _log.info("interpolating target CPO angles...")
+188 data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2001_GBMshear.scsv")
+189 cs_M0 = PchipInterpolator(
+190 _utils.remove_nans(data.equivalent_strain_M0) / 200,
+191 _utils.remove_nans(data.angle_M0),
+192 )
+193 cs_M50 = PchipInterpolator(
+194 _utils.remove_nans(data.equivalent_strain_M50) / 200,
+195 _utils.remove_nans(data.angle_M50),
+196 )
+197 cs_M200 = PchipInterpolator(
+198 _utils.remove_nans(data.equivalent_strain_M200) / 200,
+199 _utils.remove_nans(data.angle_M200),
+200 )
+201 return [cs_M0(strains), cs_M50(strains), cs_M200(strains)]
+202
+203 @classmethod
+204 def interp_GBM_FortranDRex(cls, strains):
+205 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
+206 _log.info("interpolating target CPO angles...")
+207 data = _io.read_scsv(_io.data("drexF90") / "olA_D1E4_dt50_X0_L5.scsv")
+208 data_strains = np.linspace(0, 1, 200)
+209 cs_M0 = PchipInterpolator(data_strains, _utils.remove_nans(data.M0_angle))
+210 cs_M50 = PchipInterpolator(data_strains, _utils.remove_nans(data.M50_angle))
+211 cs_M200 = PchipInterpolator(data_strains, _utils.remove_nans(data.M200_angle))
+212 return [cs_M0(strains), cs_M50(strains), cs_M200(strains)]
+213
+214 @classmethod
+215 def interp_GBS_FortranDRex(cls, strains):
+216 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
+217 _log.info("interpolating target CPO angles...")
+218 data = _io.read_scsv(_io.data("thirdparty") / "a_axis_GBS_fortran.scsv")
+219 data_strains = np.linspace(0, 1, 200)
+220 cs_X0 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X0))
+221 cs_X20 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X20))
+222 cs_X40 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X40))
+223 return [cs_X0(strains), cs_X20(strains), cs_X40(strains)]
+224
+225 @classmethod
+226 def interp_GBS_long_FortranDRex(cls, strains):
+227 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
+228 _log.info("interpolating target CPO angles...")
+229 data = _io.read_scsv(_io.data("thirdparty") / "a_axis_GBS_long_fortran.scsv")
+230 data_strains = np.linspace(0, 2.5, 500)
+231 cs_X0 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X0))
+232 cs_X20 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X20))
+233 cs_X40 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X40))
+234 return [cs_X0(strains), cs_X20(strains), cs_X40(strains)]
+235
+236 @classmethod
+237 def interp_GBS_Kaminski2004(cls, strains):
+238 """Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`."""
+239 _log.info("interpolating target CPO angles...")
+240 data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2004_GBSshear.scsv")
+241 cs_X0 = PchipInterpolator(
+242 _utils.remove_nans(data.dimensionless_time_X0),
+243 45 + _utils.remove_nans(data.angle_X0),
+244 )
+245 cs_X0d2 = PchipInterpolator(
+246 _utils.remove_nans(data.dimensionless_time_X0d2),
+247 45 + _utils.remove_nans(data.angle_X0d2),
+248 )
+249 cs_X0d4 = PchipInterpolator(
+250 _utils.remove_nans(data.dimensionless_time_X0d4),
+251 45 + _utils.remove_nans(data.angle_X0d4),
+252 )
+253 return [cs_X0(strains), cs_X0d2(strains), cs_X0d4(strains)]
+254
+255 @pytest.mark.skipif(_utils.in_ci("win32"), reason="Unable to allocate memory")
+256 def test_zero_recrystallisation(self, seed):
+257 """Check that M*=0 is a reliable switch to turn off recrystallisation."""
+258 params = _core.DefaultParams().as_dict()
+259 params["gbm_mobility"] = 0
+260 strain_rate = 1
+261 timestamps = np.linspace(0, 1, 25) # Solve until D₀t=1 (tensorial strain).
+262 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+263 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+264 mineral, _ = self.run(
+265 params,
+266 timestamps,
+267 strain_rate,
+268 get_velocity_gradient,
+269 shear_direction,
+270 seed=seed,
+271 )
+272 for fractions in mineral.fractions[1:]:
+273 nt.assert_allclose(fractions, mineral.fractions[0], atol=1e-15, rtol=0)
+274
+275 @pytest.mark.skipif(_utils.in_ci("win32"), reason="Unable to allocate memory")
+276 @pytest.mark.parametrize("gbm_mobility", [50, 100, 150])
+277 def test_grainsize_median(self, seed, gbm_mobility):
+278 """Check that M*={50,100,150}, λ*=5 causes decreasing grain size median."""
+279 params = _core.DefaultParams().as_dict()
+280 params["gbm_mobility"] = gbm_mobility
+281 params["nucleation_efficiency"] = 5
+282 strain_rate = 1
+283 timestamps = np.linspace(0, 1, 25) # Solve until D₀t=1 (tensorial strain).
+284 n_timestamps = len(timestamps)
+285 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+286 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+287 mineral, _ = self.run(
+288 params,
+289 timestamps,
+290 strain_rate,
+291 get_velocity_gradient,
+292 shear_direction,
+293 seed=seed,
+294 )
+295 medians = np.empty(n_timestamps)
+296 for i, fractions in enumerate(mineral.fractions):
+297 medians[i] = np.median(fractions)
+298
+299 # The first diff is positive (~1e-6) before nucleation sets in.
+300 nt.assert_array_less(np.diff(medians)[1:], np.full(n_timestamps - 2, 0))
+301
+302 @pytest.mark.slow
+303 @pytest.mark.parametrize("gbs_threshold", [0, 0.2, 0.4])
+304 @pytest.mark.parametrize("nucleation_efficiency", [3, 5, 10])
+305 def test_dvdx_ensemble(
+306 self, outdir, seeds_nearX45, ncpus, gbs_threshold, nucleation_efficiency
+307 ):
+308 r"""Test a-axis alignment to shear in Y direction (init. SCCS near 45° from X).
+309
+310 Velocity gradient:
+311 $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$
+312
+313 """
+314 strain_rate = 1
+315 timestamps = np.linspace(0, 1, 201) # Solve until D₀t=1 ('shear' γ=2).
+316 n_timestamps = len(timestamps)
+317 # Use `seeds` instead of `seeds_nearX45` if you have even more RAM and CPU time.
+318 _seeds = seeds_nearX45
+319 n_seeds = len(_seeds)
+320
+321 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+322 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+323
+324 gbm_mobilities = [0, 50, 125, 200]
+325 markers = ("x", "*", "d", "s")
+326
+327 _id = f"X{_io.stringify(gbs_threshold)}_L{_io.stringify(nucleation_efficiency)}"
+328 # Output setup with optional logging and data series labels.
+329 θ_fse = np.empty_like(timestamps)
+330 angles = np.empty((len(gbm_mobilities), n_seeds, n_timestamps))
+331 optional_logging = cl.nullcontext()
+332 if outdir is not None:
+333 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_ensemble_{_id}"
+334 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
+335 labels = []
+336
+337 with optional_logging:
+338 clock_start = process_time()
+339 for m, gbm_mobility in enumerate(gbm_mobilities):
+340 if m == 0:
+341 return_fse = True
+342 else:
+343 return_fse = False
+344
+345 params = {
+346 "phase_assemblage": (_core.MineralPhase.olivine,),
+347 "phase_fractions": (1.0,),
+348 "enstatite_fraction": 0.0,
+349 "stress_exponent": 1.5,
+350 "deformation_exponent": 3.5,
+351 "gbm_mobility": gbm_mobility,
+352 "gbs_threshold": gbs_threshold,
+353 "nucleation_efficiency": nucleation_efficiency,
+354 "number_of_grains": 5000,
+355 "initial_olivine_fabric": "A",
+356 }
+357
+358 _run = ft.partial(
+359 self.run,
+360 params,
+361 timestamps,
+362 strain_rate,
+363 get_velocity_gradient,
+364 shear_direction,
+365 return_fse=return_fse,
+366 )
+367 with Pool(processes=ncpus) as pool:
+368 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
+369 mineral, fse_angles = out
+370 angles[m, s, :] = [
+371 _diagnostics.smallest_angle(v, shear_direction)
+372 for v in _diagnostics.elasticity_components(
+373 _minerals.voigt_averages([mineral], params)
+374 )["hexagonal_axis"]
+375 ]
+376 # Save the whole mineral for the first seed only.
+377 if outdir is not None and s == 0:
+378 postfix = (
+379 f"M{_io.stringify(gbm_mobility)}"
+380 + f"_X{_io.stringify(gbs_threshold)}"
+381 + f"_L{_io.stringify(nucleation_efficiency)}"
+382 )
+383 mineral.save(f"{out_basepath}.npz", postfix=postfix)
+384 if return_fse:
+385 θ_fse += fse_angles
+386
+387 if return_fse:
+388 θ_fse /= n_seeds
+389
+390 if outdir is not None:
+391 labels.append(f"$M^∗$ = {gbm_mobility}")
+392
+393 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
+394
+395 # Take ensemble means and optionally plot figure.
+396 strains = timestamps * strain_rate
+397 _log.info("postprocessing results for %s", _id)
+398 result_angles = angles.mean(axis=1)
+399 result_angles_err = angles.std(axis=1)
+400
+401 if outdir is not None:
+402 schema = {
+403 "delimiter": ",",
+404 "missing": "-",
+405 "fields": [
+406 {
+407 "name": "strain",
+408 "type": "integer",
+409 "unit": "percent",
+410 "fill": 999999,
+411 }
+412 ],
+413 }
+414 np.savez(
+415 f"{out_basepath}.npz",
+416 angles=result_angles,
+417 angles_err=result_angles_err,
+418 )
+419 _io.save_scsv(
+420 f"{out_basepath}_strains.scsv",
+421 schema,
+422 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
+423 )
+424 fig, ax, colors = _vis.alignment(
+425 None,
+426 strains,
+427 result_angles,
+428 markers,
+429 labels,
+430 err=result_angles_err,
+431 θ_max=60,
+432 θ_fse=θ_fse,
+433 )
+434 fig.savefig(_io.resolve_path(f"{out_basepath}.pdf"))
+435
+436 @pytest.mark.slow
+437 def test_dvdx_GBM(self, outdir, seeds_nearX45, ncpus):
+438 r"""Test a-axis alignment to shear in Y direction (init. SCCS near 45° from X).
+439
+440 Velocity gradient:
+441 $$
+442 \bm{L} = 10^{-4} ×
+443 \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}
+444 $$
+445
+446 Results are compared to the Fortran 90 output.
+447
+448 """
+449 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+450 strain_rate = 1e-4
+451 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+452 timestamps = np.linspace(0, 1e4, 51) # Solve until D₀t=1 ('shear' γ=2).
+453 i_strain_40p = 10 # Index of 40% strain, lower strains are not relevant here.
+454 i_strain_100p = 25 # Index of 100% strain, when M*=0 matches FSE.
+455 params = _core.DefaultParams().as_dict()
+456 params["gbs_threshold"] = 0 # No GBS, to match the Fortran parameters.
+457 gbm_mobilities = (0, 10, 50, 125, 200) # Must be in ascending order.
+458 markers = ("x", ".", "*", "d", "s")
+459 # Use `seeds` instead of `seeds_nearX45` if you have even more RAM and CPU time.
+460 _seeds = seeds_nearX45
+461 n_seeds = len(_seeds)
+462 angles = np.empty((len(gbm_mobilities), n_seeds, len(timestamps)))
+463 θ_fse = np.zeros_like(timestamps)
+464 strains = timestamps * strain_rate
+465 M0_drexF90, M50_drexF90, M200_drexF90 = self.interp_GBM_FortranDRex(strains)
+466
+467 optional_logging = cl.nullcontext()
+468 if outdir is not None:
+469 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_mobility"
+470 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
+471 labels = []
+472
+473 with optional_logging:
+474 clock_start = process_time()
+475 for m, gbm_mobility in enumerate(gbm_mobilities):
+476 if m == 0:
+477 return_fse = True
+478 else:
+479 return_fse = False
+480 params["gbm_mobility"] = gbm_mobility
+481
+482 _run = ft.partial(
+483 self.run,
+484 params,
+485 timestamps,
+486 strain_rate,
+487 get_velocity_gradient,
+488 shear_direction,
+489 return_fse=True,
+490 )
+491 with Pool(processes=ncpus) as pool:
+492 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
+493 mineral, fse_angles = out
+494 angles[m, s, :] = [
+495 _diagnostics.smallest_angle(v, shear_direction)
+496 for v in _diagnostics.elasticity_components(
+497 _minerals.voigt_averages([mineral], params)
+498 )["hexagonal_axis"]
+499 ]
+500 # Save the whole mineral for the first seed only.
+501 if outdir is not None and s == 0:
+502 mineral.save(
+503 f"{out_basepath}.npz",
+504 postfix=f"M{_io.stringify(gbm_mobility)}",
+505 )
+506 if return_fse:
+507 θ_fse += fse_angles
+508
+509 if return_fse:
+510 θ_fse /= n_seeds
+511
+512 if outdir is not None:
+513 labels.append(f"$M^∗$ = {params['gbm_mobility']}")
+514
+515 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
+516
+517 # Take ensemble means and optionally plot figure.
+518 _log.info("postprocessing results...")
+519 result_angles = angles.mean(axis=1)
+520 result_angles_err = angles.std(axis=1)
+521
+522 if outdir is not None:
+523 schema = {
+524 "delimiter": ",",
+525 "missing": "-",
+526 "fields": [
+527 {
+528 "name": "strain",
+529 "type": "integer",
+530 "unit": "percent",
+531 "fill": 999999,
+532 }
+533 ],
+534 }
+535 _io.save_scsv(
+536 f"{out_basepath}_strains.scsv",
+537 schema,
+538 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
+539 )
+540 np.savez(
+541 f"{out_basepath}_angles.npz",
+542 angles=result_angles,
+543 err=result_angles_err,
+544 )
+545 fig, ax, colors = _vis.alignment(
+546 None,
+547 strains,
+548 result_angles,
+549 markers,
+550 labels,
+551 err=result_angles_err,
+552 θ_max=60,
+553 θ_fse=θ_fse,
+554 )
+555 ax.plot(strains, M0_drexF90, c=colors[0])
+556 ax.plot(strains, M50_drexF90, c=colors[2])
+557 ax.plot(strains, M200_drexF90, c=colors[4])
+558 _vis.show_Skemer2016_ShearStrainAngles(
+559 ax,
+560 ["Z&K 1200 C", "Z&K 1300 C"],
+561 ["v", "^"],
+562 ["k", "k"],
+563 ["none", None],
+564 [
+565 "Zhang & Karato, 1995\n(1473 K)",
+566 "Zhang & Karato, 1995\n(1573 K)",
+567 ],
+568 _core.MineralFabric.olivine_A,
+569 )
+570 # There is a lot of stuff on this legend, so put it outside the axes.
+571 # These values might need to be tweaked depending on the font size, etc.
+572 _legend = _utils.redraw_legend(ax, fig=fig, bbox_to_anchor=(1.66, 0.99))
+573 fig.savefig(
+574 _io.resolve_path(f"{out_basepath}.pdf"),
+575 bbox_extra_artists=(_legend,),
+576 bbox_inches="tight",
+577 )
+578
+579 # Check that GBM speeds up the alignment between 40% and 100% strain.
+580 _log.info("checking grain orientations...")
+581 for i, θ in enumerate(result_angles[:-1], start=1):
+582 nt.assert_array_less(
+583 result_angles[i][i_strain_40p:i_strain_100p],
+584 θ[i_strain_40p:i_strain_100p],
+585 )
+586
+587 # Check that M*=0 matches FSE (±1°) past 100% strain.
+588 nt.assert_allclose(
+589 result_angles[0][i_strain_100p:],
+590 θ_fse[i_strain_100p:],
+591 atol=1,
+592 rtol=0,
+593 )
+594
+595 # Check that results match Fortran output past 40% strain.
+596 nt.assert_allclose(
+597 result_angles[0][i_strain_40p:],
+598 M0_drexF90[i_strain_40p:],
+599 atol=0,
+600 rtol=0.1, # At 40% strain the match is worse than at higher strain.
+601 )
+602 nt.assert_allclose(
+603 result_angles[2][i_strain_40p:],
+604 M50_drexF90[i_strain_40p:],
+605 atol=1,
+606 rtol=0,
+607 )
+608 nt.assert_allclose(
+609 result_angles[4][i_strain_40p:],
+610 M200_drexF90[i_strain_40p:],
+611 atol=1.5,
+612 rtol=0,
+613 )
+614
+615 @pytest.mark.slow
+616 def test_GBM_calibration(self, outdir, seeds, ncpus):
+617 r"""Compare results for various values of $$M^∗$$ to A-type olivine data.
+618
+619 Velocity gradient:
+620 $$
+621 \bm{L} = 10^{-4} ×
+622 \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}
+623 $$
+624
+625 Unlike `test_dvdx_GBM`,
+626 grain boudary sliding is enabled here (see `_core.DefaultParams`).
+627 Data are provided by [Skemer & Hansen, 2016](http://dx.doi.org/10.1016/j.tecto.2015.12.003).
+628
+629 """
+630 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+631 strain_rate = 1
+632 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+633 timestamps = np.linspace(0, 3.2, 65) # Solve until D₀t=3.2 ('shear' γ=6.4).
+634 params = _core.DefaultParams().as_dict()
+635 params["number_of_grains"] = 5000
+636 gbm_mobilities = (0, 10, 50, 125) # Must be in ascending order.
+637 markers = ("x", "*", "1", ".")
+638 # Uses 100 seeds by default; use all 1000 if you have more RAM and CPU time.
+639 _seeds = seeds[:100]
+640 n_seeds = len(_seeds)
+641 angles = np.empty((len(gbm_mobilities), n_seeds, len(timestamps)))
+642 θ_fse = np.zeros_like(timestamps)
+643 strains = timestamps * strain_rate
+644
+645 optional_logging = cl.nullcontext()
+646 if outdir is not None:
+647 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_calibration"
+648 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
+649 labels = []
+650
+651 with optional_logging:
+652 clock_start = process_time()
+653 for m, gbm_mobility in enumerate(gbm_mobilities):
+654 return_fse = True if m == 0 else False
+655 params["gbm_mobility"] = gbm_mobility
+656 _run = ft.partial(
+657 self.run,
+658 params,
+659 timestamps,
+660 strain_rate,
+661 get_velocity_gradient,
+662 shear_direction,
+663 return_fse=return_fse,
+664 )
+665 with Pool(processes=ncpus) as pool:
+666 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
+667 mineral, fse_angles = out
+668 angles[m, s, :] = [
+669 _diagnostics.smallest_angle(v, shear_direction)
+670 for v in _diagnostics.elasticity_components(
+671 _minerals.voigt_averages([mineral], params)
+672 )["hexagonal_axis"]
+673 ]
+674 # Save the whole mineral for the first seed only.
+675 if outdir is not None and s == 0:
+676 mineral.save(
+677 f"{out_basepath}.npz",
+678 postfix=f"M{_io.stringify(gbm_mobility)}",
+679 )
+680 if return_fse:
+681 θ_fse += fse_angles
+682
+683 if return_fse:
+684 θ_fse /= n_seeds
+685
+686 if outdir is not None:
+687 labels.append(f"$M^∗$ = {params['gbm_mobility']}")
+688
+689 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
+690
+691 # Take ensemble means and optionally plot figure.
+692 _log.info("postprocessing results...")
+693 result_angles = angles.mean(axis=1)
+694 result_angles_err = angles.std(axis=1)
+695
+696 if outdir is not None:
+697 schema = {
+698 "delimiter": ",",
+699 "missing": "-",
+700 "fields": [
+701 {
+702 "name": "strain",
+703 "type": "integer",
+704 "unit": "percent",
+705 "fill": 999999,
+706 }
+707 ],
+708 }
+709 _io.save_scsv(
+710 f"{out_basepath}_strains.scsv",
+711 schema,
+712 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
+713 )
+714 np.savez(
+715 _io.resolve_path(f"{out_basepath}_ensemble_means.npz"),
+716 angles=result_angles,
+717 err=result_angles_err,
+718 )
+719 fig = _vis.figure(
+720 figsize=(_vis.DEFAULT_FIG_WIDTH * 3, _vis.DEFAULT_FIG_HEIGHT)
+721 )
+722 fig, ax, colors = _vis.alignment(
+723 fig.add_subplot(),
+724 strains,
+725 result_angles,
+726 markers,
+727 labels,
+728 err=result_angles_err,
+729 θ_max=80,
+730 θ_fse=θ_fse,
+731 )
+732 (
+733 _,
734 _,
735 _,
-736 _,
-737 data_Skemer2016,
-738 indices,
-739 ) = _vis.show_Skemer2016_ShearStrainAngles(
-740 ax,
-741 [
-742 "Z&K 1200 C",
-743 "Z&K 1300 C",
-744 "Skemer 2011",
-745 "Hansen 2014",
-746 "Warren 2008",
-747 "Webber 2010",
-748 "H&W 2015",
-749 ],
-750 ["v", "^", "o", "s", "v", "o", "s"],
-751 ["k", "k", "k", "k", "k", "k", "k"],
-752 ["none", "none", "none", "none", None, None, None],
-753 [
-754 "Zhang & Karato, 1995 (1473 K)",
-755 "Zhang & Karato, 1995 (1573 K)",
-756 "Skemer et al., 2011 (1500 K)",
-757 "Hansen et al., 2014 (1473 K)",
-758 "Warren et al., 2008",
-759 "Webber et al., 2010",
-760 "Hansen & Warren, 2015",
-761 ],
-762 fabric=_core.MineralFabric.olivine_A,
-763 )
-764 _legend = _utils.redraw_legend(ax, loc="upper right", ncols=3)
-765 fig.savefig(
-766 _io.resolve_path(f"{out_basepath}.pdf"),
-767 bbox_extra_artists=(_legend,),
-768 bbox_inches="tight",
-769 )
-770 r2vals = []
-771 for angles in result_angles:
-772 _angles = PchipInterpolator(strains, angles)
-773 r2 = np.sum(
-774 [
-775 (a - b) ** 2
-776 for a, b in zip(
-777 _angles(
-778 np.take(data_Skemer2016.shear_strain, indices) / 200
-779 ),
-780 np.take(data_Skemer2016.angle, indices),
-781 )
-782 ]
-783 )
-784 r2vals.append(r2)
-785 _log.info(
-786 "Sums of squared residuals (r-values) for each M∗: %s", r2vals
-787 )
-788
-789 @pytest.mark.big
-790 def test_dudz_pathline(self, outdir, seed):
-791 """Test alignment of olivine a-axis for a polycrystal advected on a pathline."""
-792 test_id = "dudz_pathline"
-793 optional_logging = cl.nullcontext()
-794 if outdir is not None:
-795 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}"
-796 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
-797
-798 with optional_logging:
-799 shear_direction = Ŋ([1, 0, 0], dtype=np.float64)
-800 strain_rate = 1e-15 # Moderate, realistic shear in the upper mantle.
-801 get_velocity, get_velocity_gradient = _velocity.simple_shear_2d(
-802 "X", "Z", strain_rate
-803 )
-804 n_timesteps = 10
-805 timestamps, get_position = _paths.get_pathline(
-806 Ŋ([1e5, 0e0, 1e5]),
-807 get_velocity,
-808 get_velocity_gradient,
-809 Ŋ([-2e5, 0e0, -2e5]),
-810 Ŋ([2e5, 0e0, 2e5]),
-811 2,
-812 regular_steps=n_timesteps,
-813 )
-814 positions = [get_position(t) for t in timestamps]
-815 velocity_gradients = [
-816 get_velocity_gradient(np.nan, Ŋ(x)) for x in positions
-817 ]
-818
-819 params = _core.DefaultParams().as_dict()
-820 params["gbm_mobility"] = 10
-821 params["number_of_grains"] = 5000
-822 olA = _minerals.Mineral(n_grains=params["number_of_grains"], seed=seed)
-823 deformation_gradient = np.eye(3)
-824 strains = np.zeros_like(timestamps)
-825 for t, time in enumerate(timestamps[:-1], start=1):
-826 strains[t] = strains[t - 1] + (
-827 _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
-828 )
-829 _log.info("step %d/%d (ε = %.2f)", t, len(timestamps) - 1, strains[t])
-830 deformation_gradient = olA.update_orientations(
-831 params,
-832 deformation_gradient,
-833 get_velocity_gradient,
-834 pathline=(time, timestamps[t], get_position),
-835 )
-836
-837 if outdir is not None:
-838 olA.save(f"{out_basepath}.npz")
-839
-840 orient_resampled, fractions_resampled = _stats.resample_orientations(
-841 olA.orientations, olA.fractions, seed=seed
-842 )
-843 # About 36GB, 26 min needed with float64. GitHub macos runner has 14GB.
-844 misorient_indices = _diagnostics.misorientation_indices(
-845 orient_resampled,
-846 _geo.LatticeSystem.orthorhombic,
-847 ncpus=3,
-848 )
-849 cpo_vectors = np.zeros((n_timesteps + 1, 3))
-850 cpo_angles = np.zeros(n_timesteps + 1)
-851 for i, matrices in enumerate(orient_resampled):
-852 cpo_vectors[i] = _diagnostics.bingham_average(
-853 matrices,
-854 axis=_minerals.OLIVINE_PRIMARY_AXIS[olA.fabric],
-855 )
-856 cpo_angles[i] = _diagnostics.smallest_angle(
-857 cpo_vectors[i], shear_direction
-858 )
-859
-860 # Check for mostly decreasing CPO angles (exclude initial condition).
-861 _log.debug("cpo angles: %s", cpo_angles)
-862 nt.assert_array_less(np.diff(cpo_angles[1:]), np.ones(n_timesteps - 1))
-863 # Check for increasing CPO strength (M-index).
-864 _log.debug("cpo strengths: %s", misorient_indices)
-865 nt.assert_array_less(
-866 np.full(n_timesteps, -0.01), np.diff(misorient_indices)
-867 )
-868 # Check that last angle is <5° (M*=125) or <10° (M*=10).
-869 assert cpo_angles[-1] < 10
-870
-871 if outdir is not None:
-872 fig, ax, _, _ = _vis.steady_box2d(
-873 None,
-874 (get_velocity, [20, 20]),
-875 (positions, Ŋ([-2e5, -2e5]), Ŋ([2e5, 2e5])),
-876 "xz",
-877 (cpo_vectors, misorient_indices),
-878 strains,
-879 )
-880 fig.savefig(f"{out_basepath}.pdf")
+736 data_Skemer2016,
+737 indices,
+738 ) = _vis.show_Skemer2016_ShearStrainAngles(
+739 ax,
+740 [
+741 "Z&K 1200 C",
+742 "Z&K 1300 C",
+743 "Skemer 2011",
+744 "Hansen 2014",
+745 "Warren 2008",
+746 "Webber 2010",
+747 "H&W 2015",
+748 ],
+749 ["v", "^", "o", "s", "v", "o", "s"],
+750 ["k", "k", "k", "k", "k", "k", "k"],
+751 ["none", "none", "none", "none", None, None, None],
+752 [
+753 "Zhang & Karato, 1995 (1473 K)",
+754 "Zhang & Karato, 1995 (1573 K)",
+755 "Skemer et al., 2011 (1500 K)",
+756 "Hansen et al., 2014 (1473 K)",
+757 "Warren et al., 2008",
+758 "Webber et al., 2010",
+759 "Hansen & Warren, 2015",
+760 ],
+761 fabric=_core.MineralFabric.olivine_A,
+762 )
+763 _legend = _utils.redraw_legend(ax, loc="upper right", ncols=3)
+764 fig.savefig(
+765 _io.resolve_path(f"{out_basepath}.pdf"),
+766 bbox_extra_artists=(_legend,),
+767 bbox_inches="tight",
+768 )
+769 r2vals = []
+770 for angles in result_angles:
+771 _angles = PchipInterpolator(strains, angles)
+772 r2 = np.sum(
+773 [
+774 (a - b) ** 2
+775 for a, b in zip(
+776 _angles(
+777 np.take(data_Skemer2016.shear_strain, indices) / 200
+778 ),
+779 np.take(data_Skemer2016.angle, indices),
+780 )
+781 ]
+782 )
+783 r2vals.append(r2)
+784 _log.info(
+785 "Sums of squared residuals (r-values) for each M∗: %s", r2vals
+786 )
+787
+788 @pytest.mark.big
+789 def test_dudz_pathline(self, outdir, seed):
+790 """Test alignment of olivine a-axis for a polycrystal advected on a pathline."""
+791 test_id = "dudz_pathline"
+792 optional_logging = cl.nullcontext()
+793 if outdir is not None:
+794 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}"
+795 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
+796
+797 with optional_logging:
+798 shear_direction = Ŋ([1, 0, 0], dtype=np.float64)
+799 strain_rate = 1e-15 # Moderate, realistic shear in the upper mantle.
+800 get_velocity, get_velocity_gradient = _velocity.simple_shear_2d(
+801 "X", "Z", strain_rate
+802 )
+803 n_timesteps = 10
+804 timestamps, get_position = _paths.get_pathline(
+805 Ŋ([1e5, 0e0, 1e5]),
+806 get_velocity,
+807 get_velocity_gradient,
+808 Ŋ([-2e5, 0e0, -2e5]),
+809 Ŋ([2e5, 0e0, 2e5]),
+810 2,
+811 regular_steps=n_timesteps,
+812 )
+813 positions = [get_position(t) for t in timestamps]
+814 velocity_gradients = [
+815 get_velocity_gradient(np.nan, Ŋ(x)) for x in positions
+816 ]
+817
+818 params = _core.DefaultParams().as_dict()
+819 params["gbm_mobility"] = 10
+820 params["number_of_grains"] = 5000
+821 olA = _minerals.Mineral(n_grains=params["number_of_grains"], seed=seed)
+822 deformation_gradient = np.eye(3)
+823 strains = np.zeros_like(timestamps)
+824 for t, time in enumerate(timestamps[:-1], start=1):
+825 strains[t] = strains[t - 1] + (
+826 _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
+827 )
+828 _log.info("step %d/%d (ε = %.2f)", t, len(timestamps) - 1, strains[t])
+829 deformation_gradient = olA.update_orientations(
+830 params,
+831 deformation_gradient,
+832 get_velocity_gradient,
+833 pathline=(time, timestamps[t], get_position),
+834 )
+835
+836 if outdir is not None:
+837 olA.save(f"{out_basepath}.npz")
+838
+839 orient_resampled, fractions_resampled = _stats.resample_orientations(
+840 olA.orientations, olA.fractions, seed=seed
+841 )
+842 # About 36GB, 26 min needed with float64. GitHub macos runner has 14GB.
+843 misorient_indices = _diagnostics.misorientation_indices(
+844 orient_resampled,
+845 _geo.LatticeSystem.orthorhombic,
+846 ncpus=3,
+847 )
+848 cpo_vectors = np.zeros((n_timesteps + 1, 3))
+849 cpo_angles = np.zeros(n_timesteps + 1)
+850 for i, matrices in enumerate(orient_resampled):
+851 cpo_vectors[i] = _diagnostics.bingham_average(
+852 matrices,
+853 axis=_minerals.OLIVINE_PRIMARY_AXIS[olA.fabric],
+854 )
+855 cpo_angles[i] = _diagnostics.smallest_angle(
+856 cpo_vectors[i], shear_direction
+857 )
+858
+859 # Check for mostly decreasing CPO angles (exclude initial condition).
+860 _log.debug("cpo angles: %s", cpo_angles)
+861 nt.assert_array_less(np.diff(cpo_angles[1:]), np.ones(n_timesteps - 1))
+862 # Check for increasing CPO strength (M-index).
+863 _log.debug("cpo strengths: %s", misorient_indices)
+864 nt.assert_array_less(
+865 np.full(n_timesteps, -0.01), np.diff(misorient_indices)
+866 )
+867 # Check that last angle is <5° (M*=125) or <10° (M*=10).
+868 assert cpo_angles[-1] < 10
+869
+870 if outdir is not None:
+871 fig, ax, _, _ = _vis.steady_box2d(
+872 None,
+873 (get_velocity, [20, 20]),
+874 (positions, Ŋ([-2e5, -2e5]), Ŋ([2e5, 2e5])),
+875 "xz",
+876 (cpo_vectors, misorient_indices),
+877 strains,
+878 )
+879 fig.savefig(f"{out_basepath}.pdf")
118 @classmethod
-119 def run(
-120 cls,
-121 params,
-122 timestamps,
-123 strain_rate,
-124 get_velocity_gradient,
-125 shear_direction,
-126 seed=None,
-127 return_fse=None,
-128 get_position=lambda t: np.zeros(3), # Stationary particles by default.
-129 ):
-130 """Reusable logic for 2D olivine (A-type) simple shear tests.
-131
-132 Returns a tuple with the mineral and the FSE angle (or `None` if `return_fse` is
-133 `None`).
-134
-135 """
-136 mineral = _minerals.Mineral(
-137 phase=_core.MineralPhase.olivine,
-138 fabric=_core.MineralFabric.olivine_A,
-139 regime=_core.DeformationRegime.matrix_dislocation,
-140 n_grains=params["number_of_grains"],
-141 seed=seed,
-142 )
-143 deformation_gradient = np.eye(3) # Undeformed initial state.
-144 θ_fse = np.empty_like(timestamps)
-145 θ_fse[0] = 45
-146
-147 for t, time in enumerate(timestamps[:-1], start=1):
-148 # Set up logging message depending on dynamic parameter and seeds.
-149 msg_start = (
-150 f"N = {params['number_of_grains']}; "
-151 + f"λ∗ = {params['nucleation_efficiency']}; "
-152 + f"X = {params['gbs_threshold']}; "
-153 + f"M∗ = {params['gbm_mobility']}; "
-154 )
-155 if seed is not None:
-156 msg_start += f"# {seed}; "
-157
-158 _log.info(msg_start + "step %s/%s (t = %s)", t, len(timestamps) - 1, time)
-159
-160 deformation_gradient = mineral.update_orientations(
-161 params,
-162 deformation_gradient,
-163 get_velocity_gradient,
-164 pathline=(time, timestamps[t], get_position),
-165 )
-166 _log.debug(
-167 "› velocity gradient = %s",
-168 get_velocity_gradient(np.nan, np.full(3, np.nan)).flatten(),
-169 )
-170 _log.debug("› strain D₀t = %.2f", strain_rate * timestamps[t])
-171 _log.debug(
-172 "› grain fractions: median = %s, max = %s, min = %s",
-173 np.median(mineral.fractions[-1]),
-174 np.max(mineral.fractions[-1]),
-175 np.min(mineral.fractions[-1]),
-176 )
-177 if return_fse:
-178 _, fse_v = _diagnostics.finite_strain(deformation_gradient)
-179 θ_fse[t] = _diagnostics.smallest_angle(fse_v, shear_direction)
-180 else:
-181 θ_fse = None
-182
-183 return mineral, θ_fse
+ 117 @classmethod
+118 def run(
+119 cls,
+120 params,
+121 timestamps,
+122 strain_rate,
+123 get_velocity_gradient,
+124 shear_direction,
+125 seed=None,
+126 return_fse=None,
+127 get_position=lambda t: np.zeros(3), # Stationary particles by default.
+128 ):
+129 """Reusable logic for 2D olivine (A-type) simple shear tests.
+130
+131 Returns a tuple with the mineral and the FSE angle (or `None` if `return_fse` is
+132 `None`).
+133
+134 """
+135 mineral = _minerals.Mineral(
+136 phase=_core.MineralPhase.olivine,
+137 fabric=_core.MineralFabric.olivine_A,
+138 regime=_core.DeformationRegime.matrix_dislocation,
+139 n_grains=params["number_of_grains"],
+140 seed=seed,
+141 )
+142 deformation_gradient = np.eye(3) # Undeformed initial state.
+143 θ_fse = np.empty_like(timestamps)
+144 θ_fse[0] = 45
+145
+146 for t, time in enumerate(timestamps[:-1], start=1):
+147 # Set up logging message depending on dynamic parameter and seeds.
+148 msg_start = (
+149 f"N = {params['number_of_grains']}; "
+150 + f"λ∗ = {params['nucleation_efficiency']}; "
+151 + f"X = {params['gbs_threshold']}; "
+152 + f"M∗ = {params['gbm_mobility']}; "
+153 )
+154 if seed is not None:
+155 msg_start += f"# {seed}; "
+156
+157 _log.info(msg_start + "step %s/%s (t = %s)", t, len(timestamps) - 1, time)
+158
+159 deformation_gradient = mineral.update_orientations(
+160 params,
+161 deformation_gradient,
+162 get_velocity_gradient,
+163 pathline=(time, timestamps[t], get_position),
+164 )
+165 _log.debug(
+166 "› velocity gradient = %s",
+167 get_velocity_gradient(np.nan, np.full(3, np.nan)).flatten(),
+168 )
+169 _log.debug("› strain D₀t = %.2f", strain_rate * timestamps[t])
+170 _log.debug(
+171 "› grain fractions: median = %s, max = %s, min = %s",
+172 np.median(mineral.fractions[-1]),
+173 np.max(mineral.fractions[-1]),
+174 np.min(mineral.fractions[-1]),
+175 )
+176 if return_fse:
+177 _, fse_v = _diagnostics.finite_strain(deformation_gradient)
+178 θ_fse[t] = _diagnostics.smallest_angle(fse_v, shear_direction)
+179 else:
+180 θ_fse = None
+181
+182 return mineral, θ_fse
185 @classmethod
-186 def interp_GBM_Kaminski2001(cls, strains):
-187 """Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`."""
-188 _log.info("interpolating target CPO angles...")
-189 data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2001_GBMshear.scsv")
-190 cs_M0 = PchipInterpolator(
-191 _utils.remove_nans(data.equivalent_strain_M0) / 200,
-192 _utils.remove_nans(data.angle_M0),
-193 )
-194 cs_M50 = PchipInterpolator(
-195 _utils.remove_nans(data.equivalent_strain_M50) / 200,
-196 _utils.remove_nans(data.angle_M50),
-197 )
-198 cs_M200 = PchipInterpolator(
-199 _utils.remove_nans(data.equivalent_strain_M200) / 200,
-200 _utils.remove_nans(data.angle_M200),
-201 )
-202 return [cs_M0(strains), cs_M50(strains), cs_M200(strains)]
+ 184 @classmethod
+185 def interp_GBM_Kaminski2001(cls, strains):
+186 """Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`."""
+187 _log.info("interpolating target CPO angles...")
+188 data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2001_GBMshear.scsv")
+189 cs_M0 = PchipInterpolator(
+190 _utils.remove_nans(data.equivalent_strain_M0) / 200,
+191 _utils.remove_nans(data.angle_M0),
+192 )
+193 cs_M50 = PchipInterpolator(
+194 _utils.remove_nans(data.equivalent_strain_M50) / 200,
+195 _utils.remove_nans(data.angle_M50),
+196 )
+197 cs_M200 = PchipInterpolator(
+198 _utils.remove_nans(data.equivalent_strain_M200) / 200,
+199 _utils.remove_nans(data.angle_M200),
+200 )
+201 return [cs_M0(strains), cs_M50(strains), cs_M200(strains)]
204 @classmethod
-205 def interp_GBM_FortranDRex(cls, strains):
-206 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
-207 _log.info("interpolating target CPO angles...")
-208 data = _io.read_scsv(_io.data("drexF90") / "olA_D1E4_dt50_X0_L5.scsv")
-209 data_strains = np.linspace(0, 1, 200)
-210 cs_M0 = PchipInterpolator(data_strains, _utils.remove_nans(data.M0_angle))
-211 cs_M50 = PchipInterpolator(data_strains, _utils.remove_nans(data.M50_angle))
-212 cs_M200 = PchipInterpolator(data_strains, _utils.remove_nans(data.M200_angle))
-213 return [cs_M0(strains), cs_M50(strains), cs_M200(strains)]
+ 203 @classmethod
+204 def interp_GBM_FortranDRex(cls, strains):
+205 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
+206 _log.info("interpolating target CPO angles...")
+207 data = _io.read_scsv(_io.data("drexF90") / "olA_D1E4_dt50_X0_L5.scsv")
+208 data_strains = np.linspace(0, 1, 200)
+209 cs_M0 = PchipInterpolator(data_strains, _utils.remove_nans(data.M0_angle))
+210 cs_M50 = PchipInterpolator(data_strains, _utils.remove_nans(data.M50_angle))
+211 cs_M200 = PchipInterpolator(data_strains, _utils.remove_nans(data.M200_angle))
+212 return [cs_M0(strains), cs_M50(strains), cs_M200(strains)]
215 @classmethod
-216 def interp_GBS_FortranDRex(cls, strains):
-217 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
-218 _log.info("interpolating target CPO angles...")
-219 data = _io.read_scsv(_io.data("thirdparty") / "a_axis_GBS_fortran.scsv")
-220 data_strains = np.linspace(0, 1, 200)
-221 cs_X0 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X0))
-222 cs_X20 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X20))
-223 cs_X40 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X40))
-224 return [cs_X0(strains), cs_X20(strains), cs_X40(strains)]
+ 214 @classmethod
+215 def interp_GBS_FortranDRex(cls, strains):
+216 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
+217 _log.info("interpolating target CPO angles...")
+218 data = _io.read_scsv(_io.data("thirdparty") / "a_axis_GBS_fortran.scsv")
+219 data_strains = np.linspace(0, 1, 200)
+220 cs_X0 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X0))
+221 cs_X20 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X20))
+222 cs_X40 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X40))
+223 return [cs_X0(strains), cs_X20(strains), cs_X40(strains)]
226 @classmethod
-227 def interp_GBS_long_FortranDRex(cls, strains):
-228 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
-229 _log.info("interpolating target CPO angles...")
-230 data = _io.read_scsv(_io.data("thirdparty") / "a_axis_GBS_long_fortran.scsv")
-231 data_strains = np.linspace(0, 2.5, 500)
-232 cs_X0 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X0))
-233 cs_X20 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X20))
-234 cs_X40 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X40))
-235 return [cs_X0(strains), cs_X20(strains), cs_X40(strains)]
+ 225 @classmethod
+226 def interp_GBS_long_FortranDRex(cls, strains):
+227 """Interpolate angles produced using 'tools/drex_forward_simpleshear.f90'."""
+228 _log.info("interpolating target CPO angles...")
+229 data = _io.read_scsv(_io.data("thirdparty") / "a_axis_GBS_long_fortran.scsv")
+230 data_strains = np.linspace(0, 2.5, 500)
+231 cs_X0 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X0))
+232 cs_X20 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X20))
+233 cs_X40 = PchipInterpolator(data_strains, _utils.remove_nans(data.a_mean_X40))
+234 return [cs_X0(strains), cs_X20(strains), cs_X40(strains)]
237 @classmethod
-238 def interp_GBS_Kaminski2004(cls, strains):
-239 """Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`."""
-240 _log.info("interpolating target CPO angles...")
-241 data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2004_GBSshear.scsv")
-242 cs_X0 = PchipInterpolator(
-243 _utils.remove_nans(data.dimensionless_time_X0),
-244 45 + _utils.remove_nans(data.angle_X0),
-245 )
-246 cs_X0d2 = PchipInterpolator(
-247 _utils.remove_nans(data.dimensionless_time_X0d2),
-248 45 + _utils.remove_nans(data.angle_X0d2),
-249 )
-250 cs_X0d4 = PchipInterpolator(
-251 _utils.remove_nans(data.dimensionless_time_X0d4),
-252 45 + _utils.remove_nans(data.angle_X0d4),
-253 )
-254 return [cs_X0(strains), cs_X0d2(strains), cs_X0d4(strains)]
+ 236 @classmethod
+237 def interp_GBS_Kaminski2004(cls, strains):
+238 """Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`."""
+239 _log.info("interpolating target CPO angles...")
+240 data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2004_GBSshear.scsv")
+241 cs_X0 = PchipInterpolator(
+242 _utils.remove_nans(data.dimensionless_time_X0),
+243 45 + _utils.remove_nans(data.angle_X0),
+244 )
+245 cs_X0d2 = PchipInterpolator(
+246 _utils.remove_nans(data.dimensionless_time_X0d2),
+247 45 + _utils.remove_nans(data.angle_X0d2),
+248 )
+249 cs_X0d4 = PchipInterpolator(
+250 _utils.remove_nans(data.dimensionless_time_X0d4),
+251 45 + _utils.remove_nans(data.angle_X0d4),
+252 )
+253 return [cs_X0(strains), cs_X0d2(strains), cs_X0d4(strains)]
256 @pytest.mark.skipif(_utils.in_ci("win32"), reason="Unable to allocate memory")
-257 def test_zero_recrystallisation(self, seed):
-258 """Check that M*=0 is a reliable switch to turn off recrystallisation."""
-259 params = _core.DefaultParams().as_dict()
-260 params["gbm_mobility"] = 0
-261 strain_rate = 1
-262 timestamps = np.linspace(0, 1, 25) # Solve until D₀t=1 (tensorial strain).
-263 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-264 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-265 mineral, _ = self.run(
-266 params,
-267 timestamps,
-268 strain_rate,
-269 get_velocity_gradient,
-270 shear_direction,
-271 seed=seed,
-272 )
-273 for fractions in mineral.fractions[1:]:
-274 nt.assert_allclose(fractions, mineral.fractions[0], atol=1e-15, rtol=0)
+ 255 @pytest.mark.skipif(_utils.in_ci("win32"), reason="Unable to allocate memory")
+256 def test_zero_recrystallisation(self, seed):
+257 """Check that M*=0 is a reliable switch to turn off recrystallisation."""
+258 params = _core.DefaultParams().as_dict()
+259 params["gbm_mobility"] = 0
+260 strain_rate = 1
+261 timestamps = np.linspace(0, 1, 25) # Solve until D₀t=1 (tensorial strain).
+262 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+263 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+264 mineral, _ = self.run(
+265 params,
+266 timestamps,
+267 strain_rate,
+268 get_velocity_gradient,
+269 shear_direction,
+270 seed=seed,
+271 )
+272 for fractions in mineral.fractions[1:]:
+273 nt.assert_allclose(fractions, mineral.fractions[0], atol=1e-15, rtol=0)
276 @pytest.mark.skipif(_utils.in_ci("win32"), reason="Unable to allocate memory")
-277 @pytest.mark.parametrize("gbm_mobility", [50, 100, 150])
-278 def test_grainsize_median(self, seed, gbm_mobility):
-279 """Check that M*={50,100,150}, λ*=5 causes decreasing grain size median."""
-280 params = _core.DefaultParams().as_dict()
-281 params["gbm_mobility"] = gbm_mobility
-282 params["nucleation_efficiency"] = 5
-283 strain_rate = 1
-284 timestamps = np.linspace(0, 1, 25) # Solve until D₀t=1 (tensorial strain).
-285 n_timestamps = len(timestamps)
-286 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-287 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-288 mineral, _ = self.run(
-289 params,
-290 timestamps,
-291 strain_rate,
-292 get_velocity_gradient,
-293 shear_direction,
-294 seed=seed,
-295 )
-296 medians = np.empty(n_timestamps)
-297 for i, fractions in enumerate(mineral.fractions):
-298 medians[i] = np.median(fractions)
-299
-300 # The first diff is positive (~1e-6) before nucleation sets in.
-301 nt.assert_array_less(np.diff(medians)[1:], np.full(n_timestamps - 2, 0))
+ 275 @pytest.mark.skipif(_utils.in_ci("win32"), reason="Unable to allocate memory")
+276 @pytest.mark.parametrize("gbm_mobility", [50, 100, 150])
+277 def test_grainsize_median(self, seed, gbm_mobility):
+278 """Check that M*={50,100,150}, λ*=5 causes decreasing grain size median."""
+279 params = _core.DefaultParams().as_dict()
+280 params["gbm_mobility"] = gbm_mobility
+281 params["nucleation_efficiency"] = 5
+282 strain_rate = 1
+283 timestamps = np.linspace(0, 1, 25) # Solve until D₀t=1 (tensorial strain).
+284 n_timestamps = len(timestamps)
+285 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+286 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+287 mineral, _ = self.run(
+288 params,
+289 timestamps,
+290 strain_rate,
+291 get_velocity_gradient,
+292 shear_direction,
+293 seed=seed,
+294 )
+295 medians = np.empty(n_timestamps)
+296 for i, fractions in enumerate(mineral.fractions):
+297 medians[i] = np.median(fractions)
+298
+299 # The first diff is positive (~1e-6) before nucleation sets in.
+300 nt.assert_array_less(np.diff(medians)[1:], np.full(n_timestamps - 2, 0))
303 @pytest.mark.slow
-304 @pytest.mark.parametrize("gbs_threshold", [0, 0.2, 0.4])
-305 @pytest.mark.parametrize("nucleation_efficiency", [3, 5, 10])
-306 def test_dvdx_ensemble(
-307 self, outdir, seeds_nearX45, ncpus, gbs_threshold, nucleation_efficiency
-308 ):
-309 r"""Test a-axis alignment to shear in Y direction (init. SCCS near 45° from X).
-310
-311 Velocity gradient:
-312 $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$
-313
-314 """
-315 strain_rate = 1
-316 timestamps = np.linspace(0, 1, 201) # Solve until D₀t=1 ('shear' γ=2).
-317 n_timestamps = len(timestamps)
-318 # Use `seeds` instead of `seeds_nearX45` if you have even more RAM and CPU time.
-319 _seeds = seeds_nearX45
-320 n_seeds = len(_seeds)
-321
-322 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-323 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-324
-325 gbm_mobilities = [0, 50, 125, 200]
-326 markers = ("x", "*", "d", "s")
-327
-328 _id = f"X{_io.stringify(gbs_threshold)}_L{_io.stringify(nucleation_efficiency)}"
-329 # Output setup with optional logging and data series labels.
-330 θ_fse = np.empty_like(timestamps)
-331 angles = np.empty((len(gbm_mobilities), n_seeds, n_timestamps))
-332 optional_logging = cl.nullcontext()
-333 if outdir is not None:
-334 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_ensemble_{_id}"
-335 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
-336 labels = []
-337
-338 with optional_logging:
-339 clock_start = process_time()
-340 for m, gbm_mobility in enumerate(gbm_mobilities):
-341 if m == 0:
-342 return_fse = True
-343 else:
-344 return_fse = False
-345
-346 params = {
-347 "phase_assemblage": (_core.MineralPhase.olivine,),
-348 "phase_fractions": (1.0,),
-349 "enstatite_fraction": 0.0,
-350 "stress_exponent": 1.5,
-351 "deformation_exponent": 3.5,
-352 "gbm_mobility": gbm_mobility,
-353 "gbs_threshold": gbs_threshold,
-354 "nucleation_efficiency": nucleation_efficiency,
-355 "number_of_grains": 5000,
-356 "initial_olivine_fabric": "A",
-357 }
-358
-359 _run = ft.partial(
-360 self.run,
-361 params,
-362 timestamps,
-363 strain_rate,
-364 get_velocity_gradient,
-365 shear_direction,
-366 return_fse=return_fse,
-367 )
-368 with Pool(processes=ncpus) as pool:
-369 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
-370 mineral, fse_angles = out
-371 angles[m, s, :] = [
-372 _diagnostics.smallest_angle(v, shear_direction)
-373 for v in _diagnostics.elasticity_components(
-374 _minerals.voigt_averages([mineral], params)
-375 )["hexagonal_axis"]
-376 ]
-377 # Save the whole mineral for the first seed only.
-378 if outdir is not None and s == 0:
-379 postfix = (
-380 f"M{_io.stringify(gbm_mobility)}"
-381 + f"_X{_io.stringify(gbs_threshold)}"
-382 + f"_L{_io.stringify(nucleation_efficiency)}"
-383 )
-384 mineral.save(f"{out_basepath}.npz", postfix=postfix)
-385 if return_fse:
-386 θ_fse += fse_angles
-387
-388 if return_fse:
-389 θ_fse /= n_seeds
-390
-391 if outdir is not None:
-392 labels.append(f"$M^∗$ = {gbm_mobility}")
-393
-394 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-395
-396 # Take ensemble means and optionally plot figure.
-397 strains = timestamps * strain_rate
-398 _log.info("postprocessing results for %s", _id)
-399 result_angles = angles.mean(axis=1)
-400 result_angles_err = angles.std(axis=1)
-401
-402 if outdir is not None:
-403 schema = {
-404 "delimiter": ",",
-405 "missing": "-",
-406 "fields": [
-407 {
-408 "name": "strain",
-409 "type": "integer",
-410 "unit": "percent",
-411 "fill": 999999,
-412 }
-413 ],
-414 }
-415 np.savez(
-416 f"{out_basepath}.npz",
-417 angles=result_angles,
-418 angles_err=result_angles_err,
-419 )
-420 _io.save_scsv(
-421 f"{out_basepath}_strains.scsv",
-422 schema,
-423 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
-424 )
-425 fig, ax, colors = _vis.alignment(
-426 None,
-427 strains,
-428 result_angles,
-429 markers,
-430 labels,
-431 err=result_angles_err,
-432 θ_max=60,
-433 θ_fse=θ_fse,
-434 )
-435 fig.savefig(_io.resolve_path(f"{out_basepath}.pdf"))
+ 302 @pytest.mark.slow
+303 @pytest.mark.parametrize("gbs_threshold", [0, 0.2, 0.4])
+304 @pytest.mark.parametrize("nucleation_efficiency", [3, 5, 10])
+305 def test_dvdx_ensemble(
+306 self, outdir, seeds_nearX45, ncpus, gbs_threshold, nucleation_efficiency
+307 ):
+308 r"""Test a-axis alignment to shear in Y direction (init. SCCS near 45° from X).
+309
+310 Velocity gradient:
+311 $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$
+312
+313 """
+314 strain_rate = 1
+315 timestamps = np.linspace(0, 1, 201) # Solve until D₀t=1 ('shear' γ=2).
+316 n_timestamps = len(timestamps)
+317 # Use `seeds` instead of `seeds_nearX45` if you have even more RAM and CPU time.
+318 _seeds = seeds_nearX45
+319 n_seeds = len(_seeds)
+320
+321 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+322 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+323
+324 gbm_mobilities = [0, 50, 125, 200]
+325 markers = ("x", "*", "d", "s")
+326
+327 _id = f"X{_io.stringify(gbs_threshold)}_L{_io.stringify(nucleation_efficiency)}"
+328 # Output setup with optional logging and data series labels.
+329 θ_fse = np.empty_like(timestamps)
+330 angles = np.empty((len(gbm_mobilities), n_seeds, n_timestamps))
+331 optional_logging = cl.nullcontext()
+332 if outdir is not None:
+333 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_ensemble_{_id}"
+334 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
+335 labels = []
+336
+337 with optional_logging:
+338 clock_start = process_time()
+339 for m, gbm_mobility in enumerate(gbm_mobilities):
+340 if m == 0:
+341 return_fse = True
+342 else:
+343 return_fse = False
+344
+345 params = {
+346 "phase_assemblage": (_core.MineralPhase.olivine,),
+347 "phase_fractions": (1.0,),
+348 "enstatite_fraction": 0.0,
+349 "stress_exponent": 1.5,
+350 "deformation_exponent": 3.5,
+351 "gbm_mobility": gbm_mobility,
+352 "gbs_threshold": gbs_threshold,
+353 "nucleation_efficiency": nucleation_efficiency,
+354 "number_of_grains": 5000,
+355 "initial_olivine_fabric": "A",
+356 }
+357
+358 _run = ft.partial(
+359 self.run,
+360 params,
+361 timestamps,
+362 strain_rate,
+363 get_velocity_gradient,
+364 shear_direction,
+365 return_fse=return_fse,
+366 )
+367 with Pool(processes=ncpus) as pool:
+368 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
+369 mineral, fse_angles = out
+370 angles[m, s, :] = [
+371 _diagnostics.smallest_angle(v, shear_direction)
+372 for v in _diagnostics.elasticity_components(
+373 _minerals.voigt_averages([mineral], params)
+374 )["hexagonal_axis"]
+375 ]
+376 # Save the whole mineral for the first seed only.
+377 if outdir is not None and s == 0:
+378 postfix = (
+379 f"M{_io.stringify(gbm_mobility)}"
+380 + f"_X{_io.stringify(gbs_threshold)}"
+381 + f"_L{_io.stringify(nucleation_efficiency)}"
+382 )
+383 mineral.save(f"{out_basepath}.npz", postfix=postfix)
+384 if return_fse:
+385 θ_fse += fse_angles
+386
+387 if return_fse:
+388 θ_fse /= n_seeds
+389
+390 if outdir is not None:
+391 labels.append(f"$M^∗$ = {gbm_mobility}")
+392
+393 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
+394
+395 # Take ensemble means and optionally plot figure.
+396 strains = timestamps * strain_rate
+397 _log.info("postprocessing results for %s", _id)
+398 result_angles = angles.mean(axis=1)
+399 result_angles_err = angles.std(axis=1)
+400
+401 if outdir is not None:
+402 schema = {
+403 "delimiter": ",",
+404 "missing": "-",
+405 "fields": [
+406 {
+407 "name": "strain",
+408 "type": "integer",
+409 "unit": "percent",
+410 "fill": 999999,
+411 }
+412 ],
+413 }
+414 np.savez(
+415 f"{out_basepath}.npz",
+416 angles=result_angles,
+417 angles_err=result_angles_err,
+418 )
+419 _io.save_scsv(
+420 f"{out_basepath}_strains.scsv",
+421 schema,
+422 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
+423 )
+424 fig, ax, colors = _vis.alignment(
+425 None,
+426 strains,
+427 result_angles,
+428 markers,
+429 labels,
+430 err=result_angles_err,
+431 θ_max=60,
+432 θ_fse=θ_fse,
+433 )
+434 fig.savefig(_io.resolve_path(f"{out_basepath}.pdf"))
437 @pytest.mark.slow
-438 def test_dvdx_GBM(self, outdir, seeds_nearX45, ncpus):
-439 r"""Test a-axis alignment to shear in Y direction (init. SCCS near 45° from X).
-440
-441 Velocity gradient:
-442 $$
-443 \bm{L} = 10^{-4} ×
-444 \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}
-445 $$
-446
-447 Results are compared to the Fortran 90 output.
-448
-449 """
-450 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-451 strain_rate = 1e-4
-452 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-453 timestamps = np.linspace(0, 1e4, 51) # Solve until D₀t=1 ('shear' γ=2).
-454 i_strain_40p = 10 # Index of 40% strain, lower strains are not relevant here.
-455 i_strain_100p = 25 # Index of 100% strain, when M*=0 matches FSE.
-456 params = _core.DefaultParams().as_dict()
-457 params["gbs_threshold"] = 0 # No GBS, to match the Fortran parameters.
-458 gbm_mobilities = (0, 10, 50, 125, 200) # Must be in ascending order.
-459 markers = ("x", ".", "*", "d", "s")
-460 # Use `seeds` instead of `seeds_nearX45` if you have even more RAM and CPU time.
-461 _seeds = seeds_nearX45
-462 n_seeds = len(_seeds)
-463 angles = np.empty((len(gbm_mobilities), n_seeds, len(timestamps)))
-464 θ_fse = np.zeros_like(timestamps)
-465 strains = timestamps * strain_rate
-466 M0_drexF90, M50_drexF90, M200_drexF90 = self.interp_GBM_FortranDRex(strains)
-467
-468 optional_logging = cl.nullcontext()
-469 if outdir is not None:
-470 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_mobility"
-471 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
-472 labels = []
-473
-474 with optional_logging:
-475 clock_start = process_time()
-476 for m, gbm_mobility in enumerate(gbm_mobilities):
-477 if m == 0:
-478 return_fse = True
-479 else:
-480 return_fse = False
-481 params["gbm_mobility"] = gbm_mobility
-482
-483 _run = ft.partial(
-484 self.run,
-485 params,
-486 timestamps,
-487 strain_rate,
-488 get_velocity_gradient,
-489 shear_direction,
-490 return_fse=True,
-491 )
-492 with Pool(processes=ncpus) as pool:
-493 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
-494 mineral, fse_angles = out
-495 angles[m, s, :] = [
-496 _diagnostics.smallest_angle(v, shear_direction)
-497 for v in _diagnostics.elasticity_components(
-498 _minerals.voigt_averages([mineral], params)
-499 )["hexagonal_axis"]
-500 ]
-501 # Save the whole mineral for the first seed only.
-502 if outdir is not None and s == 0:
-503 mineral.save(
-504 f"{out_basepath}.npz",
-505 postfix=f"M{_io.stringify(gbm_mobility)}",
-506 )
-507 if return_fse:
-508 θ_fse += fse_angles
-509
-510 if return_fse:
-511 θ_fse /= n_seeds
-512
-513 if outdir is not None:
-514 labels.append(f"$M^∗$ = {params['gbm_mobility']}")
-515
-516 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-517
-518 # Take ensemble means and optionally plot figure.
-519 _log.info("postprocessing results...")
-520 result_angles = angles.mean(axis=1)
-521 result_angles_err = angles.std(axis=1)
-522
-523 if outdir is not None:
-524 schema = {
-525 "delimiter": ",",
-526 "missing": "-",
-527 "fields": [
-528 {
-529 "name": "strain",
-530 "type": "integer",
-531 "unit": "percent",
-532 "fill": 999999,
-533 }
-534 ],
-535 }
-536 _io.save_scsv(
-537 f"{out_basepath}_strains.scsv",
-538 schema,
-539 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
-540 )
-541 np.savez(
-542 f"{out_basepath}_angles.npz",
-543 angles=result_angles,
-544 err=result_angles_err,
-545 )
-546 fig, ax, colors = _vis.alignment(
-547 None,
-548 strains,
-549 result_angles,
-550 markers,
-551 labels,
-552 err=result_angles_err,
-553 θ_max=60,
-554 θ_fse=θ_fse,
-555 )
-556 ax.plot(strains, M0_drexF90, c=colors[0])
-557 ax.plot(strains, M50_drexF90, c=colors[2])
-558 ax.plot(strains, M200_drexF90, c=colors[4])
-559 _vis.show_Skemer2016_ShearStrainAngles(
-560 ax,
-561 ["Z&K 1200 C", "Z&K 1300 C"],
-562 ["v", "^"],
-563 ["k", "k"],
-564 ["none", None],
-565 [
-566 "Zhang & Karato, 1995\n(1473 K)",
-567 "Zhang & Karato, 1995\n(1573 K)",
-568 ],
-569 _core.MineralFabric.olivine_A,
-570 )
-571 # There is a lot of stuff on this legend, so put it outside the axes.
-572 # These values might need to be tweaked depending on the font size, etc.
-573 _legend = _utils.redraw_legend(ax, fig=fig, bbox_to_anchor=(1.66, 0.99))
-574 fig.savefig(
-575 _io.resolve_path(f"{out_basepath}.pdf"),
-576 bbox_extra_artists=(_legend,),
-577 bbox_inches="tight",
-578 )
-579
-580 # Check that GBM speeds up the alignment between 40% and 100% strain.
-581 _log.info("checking grain orientations...")
-582 for i, θ in enumerate(result_angles[:-1], start=1):
-583 nt.assert_array_less(
-584 result_angles[i][i_strain_40p:i_strain_100p],
-585 θ[i_strain_40p:i_strain_100p],
-586 )
-587
-588 # Check that M*=0 matches FSE (±1°) past 100% strain.
-589 nt.assert_allclose(
-590 result_angles[0][i_strain_100p:],
-591 θ_fse[i_strain_100p:],
-592 atol=1,
-593 rtol=0,
-594 )
-595
-596 # Check that results match Fortran output past 40% strain.
-597 nt.assert_allclose(
-598 result_angles[0][i_strain_40p:],
-599 M0_drexF90[i_strain_40p:],
-600 atol=0,
-601 rtol=0.1, # At 40% strain the match is worse than at higher strain.
-602 )
-603 nt.assert_allclose(
-604 result_angles[2][i_strain_40p:],
-605 M50_drexF90[i_strain_40p:],
-606 atol=1,
-607 rtol=0,
-608 )
-609 nt.assert_allclose(
-610 result_angles[4][i_strain_40p:],
-611 M200_drexF90[i_strain_40p:],
-612 atol=1.5,
-613 rtol=0,
-614 )
+ 436 @pytest.mark.slow
+437 def test_dvdx_GBM(self, outdir, seeds_nearX45, ncpus):
+438 r"""Test a-axis alignment to shear in Y direction (init. SCCS near 45° from X).
+439
+440 Velocity gradient:
+441 $$
+442 \bm{L} = 10^{-4} ×
+443 \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}
+444 $$
+445
+446 Results are compared to the Fortran 90 output.
+447
+448 """
+449 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+450 strain_rate = 1e-4
+451 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+452 timestamps = np.linspace(0, 1e4, 51) # Solve until D₀t=1 ('shear' γ=2).
+453 i_strain_40p = 10 # Index of 40% strain, lower strains are not relevant here.
+454 i_strain_100p = 25 # Index of 100% strain, when M*=0 matches FSE.
+455 params = _core.DefaultParams().as_dict()
+456 params["gbs_threshold"] = 0 # No GBS, to match the Fortran parameters.
+457 gbm_mobilities = (0, 10, 50, 125, 200) # Must be in ascending order.
+458 markers = ("x", ".", "*", "d", "s")
+459 # Use `seeds` instead of `seeds_nearX45` if you have even more RAM and CPU time.
+460 _seeds = seeds_nearX45
+461 n_seeds = len(_seeds)
+462 angles = np.empty((len(gbm_mobilities), n_seeds, len(timestamps)))
+463 θ_fse = np.zeros_like(timestamps)
+464 strains = timestamps * strain_rate
+465 M0_drexF90, M50_drexF90, M200_drexF90 = self.interp_GBM_FortranDRex(strains)
+466
+467 optional_logging = cl.nullcontext()
+468 if outdir is not None:
+469 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_mobility"
+470 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
+471 labels = []
+472
+473 with optional_logging:
+474 clock_start = process_time()
+475 for m, gbm_mobility in enumerate(gbm_mobilities):
+476 if m == 0:
+477 return_fse = True
+478 else:
+479 return_fse = False
+480 params["gbm_mobility"] = gbm_mobility
+481
+482 _run = ft.partial(
+483 self.run,
+484 params,
+485 timestamps,
+486 strain_rate,
+487 get_velocity_gradient,
+488 shear_direction,
+489 return_fse=True,
+490 )
+491 with Pool(processes=ncpus) as pool:
+492 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
+493 mineral, fse_angles = out
+494 angles[m, s, :] = [
+495 _diagnostics.smallest_angle(v, shear_direction)
+496 for v in _diagnostics.elasticity_components(
+497 _minerals.voigt_averages([mineral], params)
+498 )["hexagonal_axis"]
+499 ]
+500 # Save the whole mineral for the first seed only.
+501 if outdir is not None and s == 0:
+502 mineral.save(
+503 f"{out_basepath}.npz",
+504 postfix=f"M{_io.stringify(gbm_mobility)}",
+505 )
+506 if return_fse:
+507 θ_fse += fse_angles
+508
+509 if return_fse:
+510 θ_fse /= n_seeds
+511
+512 if outdir is not None:
+513 labels.append(f"$M^∗$ = {params['gbm_mobility']}")
+514
+515 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
+516
+517 # Take ensemble means and optionally plot figure.
+518 _log.info("postprocessing results...")
+519 result_angles = angles.mean(axis=1)
+520 result_angles_err = angles.std(axis=1)
+521
+522 if outdir is not None:
+523 schema = {
+524 "delimiter": ",",
+525 "missing": "-",
+526 "fields": [
+527 {
+528 "name": "strain",
+529 "type": "integer",
+530 "unit": "percent",
+531 "fill": 999999,
+532 }
+533 ],
+534 }
+535 _io.save_scsv(
+536 f"{out_basepath}_strains.scsv",
+537 schema,
+538 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
+539 )
+540 np.savez(
+541 f"{out_basepath}_angles.npz",
+542 angles=result_angles,
+543 err=result_angles_err,
+544 )
+545 fig, ax, colors = _vis.alignment(
+546 None,
+547 strains,
+548 result_angles,
+549 markers,
+550 labels,
+551 err=result_angles_err,
+552 θ_max=60,
+553 θ_fse=θ_fse,
+554 )
+555 ax.plot(strains, M0_drexF90, c=colors[0])
+556 ax.plot(strains, M50_drexF90, c=colors[2])
+557 ax.plot(strains, M200_drexF90, c=colors[4])
+558 _vis.show_Skemer2016_ShearStrainAngles(
+559 ax,
+560 ["Z&K 1200 C", "Z&K 1300 C"],
+561 ["v", "^"],
+562 ["k", "k"],
+563 ["none", None],
+564 [
+565 "Zhang & Karato, 1995\n(1473 K)",
+566 "Zhang & Karato, 1995\n(1573 K)",
+567 ],
+568 _core.MineralFabric.olivine_A,
+569 )
+570 # There is a lot of stuff on this legend, so put it outside the axes.
+571 # These values might need to be tweaked depending on the font size, etc.
+572 _legend = _utils.redraw_legend(ax, fig=fig, bbox_to_anchor=(1.66, 0.99))
+573 fig.savefig(
+574 _io.resolve_path(f"{out_basepath}.pdf"),
+575 bbox_extra_artists=(_legend,),
+576 bbox_inches="tight",
+577 )
+578
+579 # Check that GBM speeds up the alignment between 40% and 100% strain.
+580 _log.info("checking grain orientations...")
+581 for i, θ in enumerate(result_angles[:-1], start=1):
+582 nt.assert_array_less(
+583 result_angles[i][i_strain_40p:i_strain_100p],
+584 θ[i_strain_40p:i_strain_100p],
+585 )
+586
+587 # Check that M*=0 matches FSE (±1°) past 100% strain.
+588 nt.assert_allclose(
+589 result_angles[0][i_strain_100p:],
+590 θ_fse[i_strain_100p:],
+591 atol=1,
+592 rtol=0,
+593 )
+594
+595 # Check that results match Fortran output past 40% strain.
+596 nt.assert_allclose(
+597 result_angles[0][i_strain_40p:],
+598 M0_drexF90[i_strain_40p:],
+599 atol=0,
+600 rtol=0.1, # At 40% strain the match is worse than at higher strain.
+601 )
+602 nt.assert_allclose(
+603 result_angles[2][i_strain_40p:],
+604 M50_drexF90[i_strain_40p:],
+605 atol=1,
+606 rtol=0,
+607 )
+608 nt.assert_allclose(
+609 result_angles[4][i_strain_40p:],
+610 M200_drexF90[i_strain_40p:],
+611 atol=1.5,
+612 rtol=0,
+613 )
616 @pytest.mark.slow
-617 def test_GBM_calibration(self, outdir, seeds, ncpus):
-618 r"""Compare results for various values of $$M^∗$$ to A-type olivine data.
-619
-620 Velocity gradient:
-621 $$
-622 \bm{L} = 10^{-4} ×
-623 \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}
-624 $$
-625
-626 Unlike `test_dvdx_GBM`,
-627 grain boudary sliding is enabled here (see `_core.DefaultParams`).
-628 Data are provided by [Skemer & Hansen, 2016](http://dx.doi.org/10.1016/j.tecto.2015.12.003).
-629
-630 """
-631 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
-632 strain_rate = 1
-633 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
-634 timestamps = np.linspace(0, 3.2, 65) # Solve until D₀t=3.2 ('shear' γ=6.4).
-635 params = _core.DefaultParams().as_dict()
-636 params["number_of_grains"] = 5000
-637 gbm_mobilities = (0, 10, 50, 125) # Must be in ascending order.
-638 markers = ("x", "*", "1", ".")
-639 # Uses 100 seeds by default; use all 1000 if you have more RAM and CPU time.
-640 _seeds = seeds[:100]
-641 n_seeds = len(_seeds)
-642 angles = np.empty((len(gbm_mobilities), n_seeds, len(timestamps)))
-643 θ_fse = np.zeros_like(timestamps)
-644 strains = timestamps * strain_rate
-645
-646 optional_logging = cl.nullcontext()
-647 if outdir is not None:
-648 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_calibration"
-649 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
-650 labels = []
-651
-652 with optional_logging:
-653 clock_start = process_time()
-654 for m, gbm_mobility in enumerate(gbm_mobilities):
-655 return_fse = True if m == 0 else False
-656 params["gbm_mobility"] = gbm_mobility
-657 _run = ft.partial(
-658 self.run,
-659 params,
-660 timestamps,
-661 strain_rate,
-662 get_velocity_gradient,
-663 shear_direction,
-664 return_fse=return_fse,
-665 )
-666 with Pool(processes=ncpus) as pool:
-667 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
-668 mineral, fse_angles = out
-669 angles[m, s, :] = [
-670 _diagnostics.smallest_angle(v, shear_direction)
-671 for v in _diagnostics.elasticity_components(
-672 _minerals.voigt_averages([mineral], params)
-673 )["hexagonal_axis"]
-674 ]
-675 # Save the whole mineral for the first seed only.
-676 if outdir is not None and s == 0:
-677 mineral.save(
-678 f"{out_basepath}.npz",
-679 postfix=f"M{_io.stringify(gbm_mobility)}",
-680 )
-681 if return_fse:
-682 θ_fse += fse_angles
-683
-684 if return_fse:
-685 θ_fse /= n_seeds
-686
-687 if outdir is not None:
-688 labels.append(f"$M^∗$ = {params['gbm_mobility']}")
-689
-690 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-691
-692 # Take ensemble means and optionally plot figure.
-693 _log.info("postprocessing results...")
-694 result_angles = angles.mean(axis=1)
-695 result_angles_err = angles.std(axis=1)
-696
-697 if outdir is not None:
-698 schema = {
-699 "delimiter": ",",
-700 "missing": "-",
-701 "fields": [
-702 {
-703 "name": "strain",
-704 "type": "integer",
-705 "unit": "percent",
-706 "fill": 999999,
-707 }
-708 ],
-709 }
-710 _io.save_scsv(
-711 f"{out_basepath}_strains.scsv",
-712 schema,
-713 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
-714 )
-715 np.savez(
-716 _io.resolve_path(f"{out_basepath}_ensemble_means.npz"),
-717 angles=result_angles,
-718 err=result_angles_err,
-719 )
-720 fig = _vis.figure(
-721 figsize=(_vis.DEFAULT_FIG_WIDTH * 3, _vis.DEFAULT_FIG_HEIGHT)
-722 )
-723 fig, ax, colors = _vis.alignment(
-724 fig.add_subplot(),
-725 strains,
-726 result_angles,
-727 markers,
-728 labels,
-729 err=result_angles_err,
-730 θ_max=80,
-731 θ_fse=θ_fse,
-732 )
-733 (
+ 615 @pytest.mark.slow
+616 def test_GBM_calibration(self, outdir, seeds, ncpus):
+617 r"""Compare results for various values of $$M^∗$$ to A-type olivine data.
+618
+619 Velocity gradient:
+620 $$
+621 \bm{L} = 10^{-4} ×
+622 \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}
+623 $$
+624
+625 Unlike `test_dvdx_GBM`,
+626 grain boudary sliding is enabled here (see `_core.DefaultParams`).
+627 Data are provided by [Skemer & Hansen, 2016](http://dx.doi.org/10.1016/j.tecto.2015.12.003).
+628
+629 """
+630 shear_direction = Ŋ([0, 1, 0], dtype=np.float64)
+631 strain_rate = 1
+632 _, get_velocity_gradient = _velocity.simple_shear_2d("Y", "X", strain_rate)
+633 timestamps = np.linspace(0, 3.2, 65) # Solve until D₀t=3.2 ('shear' γ=6.4).
+634 params = _core.DefaultParams().as_dict()
+635 params["number_of_grains"] = 5000
+636 gbm_mobilities = (0, 10, 50, 125) # Must be in ascending order.
+637 markers = ("x", "*", "1", ".")
+638 # Uses 100 seeds by default; use all 1000 if you have more RAM and CPU time.
+639 _seeds = seeds[:100]
+640 n_seeds = len(_seeds)
+641 angles = np.empty((len(gbm_mobilities), n_seeds, len(timestamps)))
+642 θ_fse = np.zeros_like(timestamps)
+643 strains = timestamps * strain_rate
+644
+645 optional_logging = cl.nullcontext()
+646 if outdir is not None:
+647 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_calibration"
+648 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
+649 labels = []
+650
+651 with optional_logging:
+652 clock_start = process_time()
+653 for m, gbm_mobility in enumerate(gbm_mobilities):
+654 return_fse = True if m == 0 else False
+655 params["gbm_mobility"] = gbm_mobility
+656 _run = ft.partial(
+657 self.run,
+658 params,
+659 timestamps,
+660 strain_rate,
+661 get_velocity_gradient,
+662 shear_direction,
+663 return_fse=return_fse,
+664 )
+665 with Pool(processes=ncpus) as pool:
+666 for s, out in enumerate(pool.imap_unordered(_run, _seeds)):
+667 mineral, fse_angles = out
+668 angles[m, s, :] = [
+669 _diagnostics.smallest_angle(v, shear_direction)
+670 for v in _diagnostics.elasticity_components(
+671 _minerals.voigt_averages([mineral], params)
+672 )["hexagonal_axis"]
+673 ]
+674 # Save the whole mineral for the first seed only.
+675 if outdir is not None and s == 0:
+676 mineral.save(
+677 f"{out_basepath}.npz",
+678 postfix=f"M{_io.stringify(gbm_mobility)}",
+679 )
+680 if return_fse:
+681 θ_fse += fse_angles
+682
+683 if return_fse:
+684 θ_fse /= n_seeds
+685
+686 if outdir is not None:
+687 labels.append(f"$M^∗$ = {params['gbm_mobility']}")
+688
+689 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
+690
+691 # Take ensemble means and optionally plot figure.
+692 _log.info("postprocessing results...")
+693 result_angles = angles.mean(axis=1)
+694 result_angles_err = angles.std(axis=1)
+695
+696 if outdir is not None:
+697 schema = {
+698 "delimiter": ",",
+699 "missing": "-",
+700 "fields": [
+701 {
+702 "name": "strain",
+703 "type": "integer",
+704 "unit": "percent",
+705 "fill": 999999,
+706 }
+707 ],
+708 }
+709 _io.save_scsv(
+710 f"{out_basepath}_strains.scsv",
+711 schema,
+712 [[int(D * 200) for D in strains]], # Shear strain % is 200 * D₀.
+713 )
+714 np.savez(
+715 _io.resolve_path(f"{out_basepath}_ensemble_means.npz"),
+716 angles=result_angles,
+717 err=result_angles_err,
+718 )
+719 fig = _vis.figure(
+720 figsize=(_vis.DEFAULT_FIG_WIDTH * 3, _vis.DEFAULT_FIG_HEIGHT)
+721 )
+722 fig, ax, colors = _vis.alignment(
+723 fig.add_subplot(),
+724 strains,
+725 result_angles,
+726 markers,
+727 labels,
+728 err=result_angles_err,
+729 θ_max=80,
+730 θ_fse=θ_fse,
+731 )
+732 (
+733 _,
734 _,
735 _,
-736 _,
-737 data_Skemer2016,
-738 indices,
-739 ) = _vis.show_Skemer2016_ShearStrainAngles(
-740 ax,
-741 [
-742 "Z&K 1200 C",
-743 "Z&K 1300 C",
-744 "Skemer 2011",
-745 "Hansen 2014",
-746 "Warren 2008",
-747 "Webber 2010",
-748 "H&W 2015",
-749 ],
-750 ["v", "^", "o", "s", "v", "o", "s"],
-751 ["k", "k", "k", "k", "k", "k", "k"],
-752 ["none", "none", "none", "none", None, None, None],
-753 [
-754 "Zhang & Karato, 1995 (1473 K)",
-755 "Zhang & Karato, 1995 (1573 K)",
-756 "Skemer et al., 2011 (1500 K)",
-757 "Hansen et al., 2014 (1473 K)",
-758 "Warren et al., 2008",
-759 "Webber et al., 2010",
-760 "Hansen & Warren, 2015",
-761 ],
-762 fabric=_core.MineralFabric.olivine_A,
-763 )
-764 _legend = _utils.redraw_legend(ax, loc="upper right", ncols=3)
-765 fig.savefig(
-766 _io.resolve_path(f"{out_basepath}.pdf"),
-767 bbox_extra_artists=(_legend,),
-768 bbox_inches="tight",
-769 )
-770 r2vals = []
-771 for angles in result_angles:
-772 _angles = PchipInterpolator(strains, angles)
-773 r2 = np.sum(
-774 [
-775 (a - b) ** 2
-776 for a, b in zip(
-777 _angles(
-778 np.take(data_Skemer2016.shear_strain, indices) / 200
-779 ),
-780 np.take(data_Skemer2016.angle, indices),
-781 )
-782 ]
-783 )
-784 r2vals.append(r2)
-785 _log.info(
-786 "Sums of squared residuals (r-values) for each M∗: %s", r2vals
-787 )
+736 data_Skemer2016,
+737 indices,
+738 ) = _vis.show_Skemer2016_ShearStrainAngles(
+739 ax,
+740 [
+741 "Z&K 1200 C",
+742 "Z&K 1300 C",
+743 "Skemer 2011",
+744 "Hansen 2014",
+745 "Warren 2008",
+746 "Webber 2010",
+747 "H&W 2015",
+748 ],
+749 ["v", "^", "o", "s", "v", "o", "s"],
+750 ["k", "k", "k", "k", "k", "k", "k"],
+751 ["none", "none", "none", "none", None, None, None],
+752 [
+753 "Zhang & Karato, 1995 (1473 K)",
+754 "Zhang & Karato, 1995 (1573 K)",
+755 "Skemer et al., 2011 (1500 K)",
+756 "Hansen et al., 2014 (1473 K)",
+757 "Warren et al., 2008",
+758 "Webber et al., 2010",
+759 "Hansen & Warren, 2015",
+760 ],
+761 fabric=_core.MineralFabric.olivine_A,
+762 )
+763 _legend = _utils.redraw_legend(ax, loc="upper right", ncols=3)
+764 fig.savefig(
+765 _io.resolve_path(f"{out_basepath}.pdf"),
+766 bbox_extra_artists=(_legend,),
+767 bbox_inches="tight",
+768 )
+769 r2vals = []
+770 for angles in result_angles:
+771 _angles = PchipInterpolator(strains, angles)
+772 r2 = np.sum(
+773 [
+774 (a - b) ** 2
+775 for a, b in zip(
+776 _angles(
+777 np.take(data_Skemer2016.shear_strain, indices) / 200
+778 ),
+779 np.take(data_Skemer2016.angle, indices),
+780 )
+781 ]
+782 )
+783 r2vals.append(r2)
+784 _log.info(
+785 "Sums of squared residuals (r-values) for each M∗: %s", r2vals
+786 )
789 @pytest.mark.big
-790 def test_dudz_pathline(self, outdir, seed):
-791 """Test alignment of olivine a-axis for a polycrystal advected on a pathline."""
-792 test_id = "dudz_pathline"
-793 optional_logging = cl.nullcontext()
-794 if outdir is not None:
-795 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}"
-796 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
-797
-798 with optional_logging:
-799 shear_direction = Ŋ([1, 0, 0], dtype=np.float64)
-800 strain_rate = 1e-15 # Moderate, realistic shear in the upper mantle.
-801 get_velocity, get_velocity_gradient = _velocity.simple_shear_2d(
-802 "X", "Z", strain_rate
-803 )
-804 n_timesteps = 10
-805 timestamps, get_position = _paths.get_pathline(
-806 Ŋ([1e5, 0e0, 1e5]),
-807 get_velocity,
-808 get_velocity_gradient,
-809 Ŋ([-2e5, 0e0, -2e5]),
-810 Ŋ([2e5, 0e0, 2e5]),
-811 2,
-812 regular_steps=n_timesteps,
-813 )
-814 positions = [get_position(t) for t in timestamps]
-815 velocity_gradients = [
-816 get_velocity_gradient(np.nan, Ŋ(x)) for x in positions
-817 ]
-818
-819 params = _core.DefaultParams().as_dict()
-820 params["gbm_mobility"] = 10
-821 params["number_of_grains"] = 5000
-822 olA = _minerals.Mineral(n_grains=params["number_of_grains"], seed=seed)
-823 deformation_gradient = np.eye(3)
-824 strains = np.zeros_like(timestamps)
-825 for t, time in enumerate(timestamps[:-1], start=1):
-826 strains[t] = strains[t - 1] + (
-827 _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
-828 )
-829 _log.info("step %d/%d (ε = %.2f)", t, len(timestamps) - 1, strains[t])
-830 deformation_gradient = olA.update_orientations(
-831 params,
-832 deformation_gradient,
-833 get_velocity_gradient,
-834 pathline=(time, timestamps[t], get_position),
-835 )
-836
-837 if outdir is not None:
-838 olA.save(f"{out_basepath}.npz")
-839
-840 orient_resampled, fractions_resampled = _stats.resample_orientations(
-841 olA.orientations, olA.fractions, seed=seed
-842 )
-843 # About 36GB, 26 min needed with float64. GitHub macos runner has 14GB.
-844 misorient_indices = _diagnostics.misorientation_indices(
-845 orient_resampled,
-846 _geo.LatticeSystem.orthorhombic,
-847 ncpus=3,
-848 )
-849 cpo_vectors = np.zeros((n_timesteps + 1, 3))
-850 cpo_angles = np.zeros(n_timesteps + 1)
-851 for i, matrices in enumerate(orient_resampled):
-852 cpo_vectors[i] = _diagnostics.bingham_average(
-853 matrices,
-854 axis=_minerals.OLIVINE_PRIMARY_AXIS[olA.fabric],
-855 )
-856 cpo_angles[i] = _diagnostics.smallest_angle(
-857 cpo_vectors[i], shear_direction
-858 )
-859
-860 # Check for mostly decreasing CPO angles (exclude initial condition).
-861 _log.debug("cpo angles: %s", cpo_angles)
-862 nt.assert_array_less(np.diff(cpo_angles[1:]), np.ones(n_timesteps - 1))
-863 # Check for increasing CPO strength (M-index).
-864 _log.debug("cpo strengths: %s", misorient_indices)
-865 nt.assert_array_less(
-866 np.full(n_timesteps, -0.01), np.diff(misorient_indices)
-867 )
-868 # Check that last angle is <5° (M*=125) or <10° (M*=10).
-869 assert cpo_angles[-1] < 10
-870
-871 if outdir is not None:
-872 fig, ax, _, _ = _vis.steady_box2d(
-873 None,
-874 (get_velocity, [20, 20]),
-875 (positions, Ŋ([-2e5, -2e5]), Ŋ([2e5, 2e5])),
-876 "xz",
-877 (cpo_vectors, misorient_indices),
-878 strains,
-879 )
-880 fig.savefig(f"{out_basepath}.pdf")
+ 788 @pytest.mark.big
+789 def test_dudz_pathline(self, outdir, seed):
+790 """Test alignment of olivine a-axis for a polycrystal advected on a pathline."""
+791 test_id = "dudz_pathline"
+792 optional_logging = cl.nullcontext()
+793 if outdir is not None:
+794 out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_{test_id}"
+795 optional_logging = _io.logfile_enable(f"{out_basepath}.log")
+796
+797 with optional_logging:
+798 shear_direction = Ŋ([1, 0, 0], dtype=np.float64)
+799 strain_rate = 1e-15 # Moderate, realistic shear in the upper mantle.
+800 get_velocity, get_velocity_gradient = _velocity.simple_shear_2d(
+801 "X", "Z", strain_rate
+802 )
+803 n_timesteps = 10
+804 timestamps, get_position = _paths.get_pathline(
+805 Ŋ([1e5, 0e0, 1e5]),
+806 get_velocity,
+807 get_velocity_gradient,
+808 Ŋ([-2e5, 0e0, -2e5]),
+809 Ŋ([2e5, 0e0, 2e5]),
+810 2,
+811 regular_steps=n_timesteps,
+812 )
+813 positions = [get_position(t) for t in timestamps]
+814 velocity_gradients = [
+815 get_velocity_gradient(np.nan, Ŋ(x)) for x in positions
+816 ]
+817
+818 params = _core.DefaultParams().as_dict()
+819 params["gbm_mobility"] = 10
+820 params["number_of_grains"] = 5000
+821 olA = _minerals.Mineral(n_grains=params["number_of_grains"], seed=seed)
+822 deformation_gradient = np.eye(3)
+823 strains = np.zeros_like(timestamps)
+824 for t, time in enumerate(timestamps[:-1], start=1):
+825 strains[t] = strains[t - 1] + (
+826 _utils.strain_increment(timestamps[t] - time, velocity_gradients[t])
+827 )
+828 _log.info("step %d/%d (ε = %.2f)", t, len(timestamps) - 1, strains[t])
+829 deformation_gradient = olA.update_orientations(
+830 params,
+831 deformation_gradient,
+832 get_velocity_gradient,
+833 pathline=(time, timestamps[t], get_position),
+834 )
+835
+836 if outdir is not None:
+837 olA.save(f"{out_basepath}.npz")
+838
+839 orient_resampled, fractions_resampled = _stats.resample_orientations(
+840 olA.orientations, olA.fractions, seed=seed
+841 )
+842 # About 36GB, 26 min needed with float64. GitHub macos runner has 14GB.
+843 misorient_indices = _diagnostics.misorientation_indices(
+844 orient_resampled,
+845 _geo.LatticeSystem.orthorhombic,
+846 ncpus=3,
+847 )
+848 cpo_vectors = np.zeros((n_timesteps + 1, 3))
+849 cpo_angles = np.zeros(n_timesteps + 1)
+850 for i, matrices in enumerate(orient_resampled):
+851 cpo_vectors[i] = _diagnostics.bingham_average(
+852 matrices,
+853 axis=_minerals.OLIVINE_PRIMARY_AXIS[olA.fabric],
+854 )
+855 cpo_angles[i] = _diagnostics.smallest_angle(
+856 cpo_vectors[i], shear_direction
+857 )
+858
+859 # Check for mostly decreasing CPO angles (exclude initial condition).
+860 _log.debug("cpo angles: %s", cpo_angles)
+861 nt.assert_array_less(np.diff(cpo_angles[1:]), np.ones(n_timesteps - 1))
+862 # Check for increasing CPO strength (M-index).
+863 _log.debug("cpo strengths: %s", misorient_indices)
+864 nt.assert_array_less(
+865 np.full(n_timesteps, -0.01), np.diff(misorient_indices)
+866 )
+867 # Check that last angle is <5° (M*=125) or <10° (M*=10).
+868 assert cpo_angles[-1] < 10
+869
+870 if outdir is not None:
+871 fig, ax, _, _ = _vis.steady_box2d(
+872 None,
+873 (get_velocity, [20, 20]),
+874 (positions, Ŋ([-2e5, -2e5]), Ŋ([2e5, 2e5])),
+875 "xz",
+876 (cpo_vectors, misorient_indices),
+877 strains,
+878 )
+879 fig.savefig(f"{out_basepath}.pdf")
19
20Pool, HAS_RAY = _utils.import_proc_pool()
21if HAS_RAY:
- 22 import ray
+ 22 import ray # noqa: F401
23
- 24 from pydrex import distributed as _dstr
+ 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"
@@ -287,77 +287,83 @@
192 ]
193 )
194 olA_downsampled, _ = _stats.resample_orientations(
-195 olivine.orientations, olivine.fractions, seed=_seeds[s], n_samples=1000
-196 )
-197 olA_strength[s, :] = _diagnostics.misorientation_indices(
-198 olA_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool
+195 olivine.orientations,
+196 olivine.fractions,
+197 seed=_seeds[s],
+198 n_samples=1000,
199 )
-200
-201 del olivine, olA_resampled, olA_mean_vectors
-202
-203 _log.info("%s; # %d; postprocessing enstatite...", _id, _seeds[s])
-204 ens_resampled, _ = _stats.resample_orientations(
-205 enstatite.orientations, enstatite.fractions, seed=_seeds[s]
-206 )
-207 ens_mean_vectors = np.array(
-208 [
-209 _diagnostics.bingham_average(dcm, axis="a")
-210 for dcm in ens_resampled
-211 ]
-212 )
-213 ens_from_proj_XZ[s, :] = np.array(
-214 [
-215 _diagnostics.smallest_angle(v, v - v * [0, 1, 0])
-216 for v in ens_mean_vectors
-217 ]
-218 )
-219 ens_from_proj_YX[s, :] = np.array(
-220 [
-221 _diagnostics.smallest_angle(v, v - v * [0, 0, 1])
-222 for v in ens_mean_vectors
-223 ]
-224 )
-225 ens_downsampled, _ = _stats.resample_orientations(
-226 enstatite.orientations, enstatite.fractions, seed=_seeds[s], n_samples=1000
+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_strength[s, :] = _diagnostics.misorientation_indices(
-229 ens_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool
-230 )
-231 del enstatite, ens_resampled, ens_mean_vectors
-232
-233 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-234 _log.info("calculating ensemble averages...")
-235 olA_from_proj_XZ_mean = olA_from_proj_XZ.mean(axis=0)
-236 olA_from_proj_XZ_err = olA_from_proj_XZ.std(axis=0)
-237
-238 olA_from_proj_YX_mean = olA_from_proj_YX.mean(axis=0)
-239 olA_from_proj_YX_err = olA_from_proj_YX.std(axis=0)
-240
-241 ens_from_proj_XZ_mean = ens_from_proj_XZ.mean(axis=0)
-242 ens_from_proj_XZ_err = ens_from_proj_XZ.std(axis=0)
+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 ens_from_proj_YX_mean = ens_from_proj_YX.mean(axis=0)
-245 ens_from_proj_YX_err = ens_from_proj_YX.std(axis=0)
+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 if outdir is not None:
-248 np.savez(
-249 f"{out_basepath}.npz",
-250 olA_from_proj_XZ_mean=olA_from_proj_XZ_mean,
-251 olA_from_proj_XZ_err=olA_from_proj_XZ_err,
-252 olA_from_proj_YX_mean=olA_from_proj_YX_mean,
-253 olA_from_proj_YX_err=olA_from_proj_YX_err,
-254 ens_from_proj_XZ_mean=ens_from_proj_XZ_mean,
-255 ens_from_proj_XZ_err=ens_from_proj_XZ_err,
-256 ens_from_proj_YX_mean=ens_from_proj_YX_mean,
-257 ens_from_proj_YX_err=ens_from_proj_YX_err,
-258 )
-259 np.savez(
-260 f"{out_basepath}_strength.npz",
-261 olA_mean=olA_strength.mean(axis=0),
-262 ens_mean=ens_strength.mean(axis=0),
-263 olA_err=olA_strength.std(axis=0),
-264 ens_err=ens_strength.std(axis=0),
-265 )
+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 )
193 ]
194 )
195 olA_downsampled, _ = _stats.resample_orientations(
-196 olivine.orientations, olivine.fractions, seed=_seeds[s], n_samples=1000
-197 )
-198 olA_strength[s, :] = _diagnostics.misorientation_indices(
-199 olA_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool
+196 olivine.orientations,
+197 olivine.fractions,
+198 seed=_seeds[s],
+199 n_samples=1000,
200 )
-201
-202 del olivine, olA_resampled, olA_mean_vectors
-203
-204 _log.info("%s; # %d; postprocessing enstatite...", _id, _seeds[s])
-205 ens_resampled, _ = _stats.resample_orientations(
-206 enstatite.orientations, enstatite.fractions, seed=_seeds[s]
-207 )
-208 ens_mean_vectors = np.array(
-209 [
-210 _diagnostics.bingham_average(dcm, axis="a")
-211 for dcm in ens_resampled
-212 ]
-213 )
-214 ens_from_proj_XZ[s, :] = np.array(
-215 [
-216 _diagnostics.smallest_angle(v, v - v * [0, 1, 0])
-217 for v in ens_mean_vectors
-218 ]
-219 )
-220 ens_from_proj_YX[s, :] = np.array(
-221 [
-222 _diagnostics.smallest_angle(v, v - v * [0, 0, 1])
-223 for v in ens_mean_vectors
-224 ]
-225 )
-226 ens_downsampled, _ = _stats.resample_orientations(
-227 enstatite.orientations, enstatite.fractions, seed=_seeds[s], n_samples=1000
+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_strength[s, :] = _diagnostics.misorientation_indices(
-230 ens_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool
-231 )
-232 del enstatite, ens_resampled, ens_mean_vectors
-233
-234 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-235 _log.info("calculating ensemble averages...")
-236 olA_from_proj_XZ_mean = olA_from_proj_XZ.mean(axis=0)
-237 olA_from_proj_XZ_err = olA_from_proj_XZ.std(axis=0)
-238
-239 olA_from_proj_YX_mean = olA_from_proj_YX.mean(axis=0)
-240 olA_from_proj_YX_err = olA_from_proj_YX.std(axis=0)
-241
-242 ens_from_proj_XZ_mean = ens_from_proj_XZ.mean(axis=0)
-243 ens_from_proj_XZ_err = ens_from_proj_XZ.std(axis=0)
+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 ens_from_proj_YX_mean = ens_from_proj_YX.mean(axis=0)
-246 ens_from_proj_YX_err = ens_from_proj_YX.std(axis=0)
+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 if outdir is not None:
-249 np.savez(
-250 f"{out_basepath}.npz",
-251 olA_from_proj_XZ_mean=olA_from_proj_XZ_mean,
-252 olA_from_proj_XZ_err=olA_from_proj_XZ_err,
-253 olA_from_proj_YX_mean=olA_from_proj_YX_mean,
-254 olA_from_proj_YX_err=olA_from_proj_YX_err,
-255 ens_from_proj_XZ_mean=ens_from_proj_XZ_mean,
-256 ens_from_proj_XZ_err=ens_from_proj_XZ_err,
-257 ens_from_proj_YX_mean=ens_from_proj_YX_mean,
-258 ens_from_proj_YX_err=ens_from_proj_YX_err,
-259 )
-260 np.savez(
-261 f"{out_basepath}_strength.npz",
-262 olA_mean=olA_strength.mean(axis=0),
-263 ens_mean=ens_strength.mean(axis=0),
-264 olA_err=olA_strength.std(axis=0),
-265 ens_err=ens_strength.std(axis=0),
-266 )
+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 )
@@ -836,77 +848,83 @@
193 ]
194 )
195 olA_downsampled, _ = _stats.resample_orientations(
-196 olivine.orientations, olivine.fractions, seed=_seeds[s], n_samples=1000
-197 )
-198 olA_strength[s, :] = _diagnostics.misorientation_indices(
-199 olA_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool
+196 olivine.orientations,
+197 olivine.fractions,
+198 seed=_seeds[s],
+199 n_samples=1000,
200 )
-201
-202 del olivine, olA_resampled, olA_mean_vectors
-203
-204 _log.info("%s; # %d; postprocessing enstatite...", _id, _seeds[s])
-205 ens_resampled, _ = _stats.resample_orientations(
-206 enstatite.orientations, enstatite.fractions, seed=_seeds[s]
-207 )
-208 ens_mean_vectors = np.array(
-209 [
-210 _diagnostics.bingham_average(dcm, axis="a")
-211 for dcm in ens_resampled
-212 ]
-213 )
-214 ens_from_proj_XZ[s, :] = np.array(
-215 [
-216 _diagnostics.smallest_angle(v, v - v * [0, 1, 0])
-217 for v in ens_mean_vectors
-218 ]
-219 )
-220 ens_from_proj_YX[s, :] = np.array(
-221 [
-222 _diagnostics.smallest_angle(v, v - v * [0, 0, 1])
-223 for v in ens_mean_vectors
-224 ]
-225 )
-226 ens_downsampled, _ = _stats.resample_orientations(
-227 enstatite.orientations, enstatite.fractions, seed=_seeds[s], n_samples=1000
+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_strength[s, :] = _diagnostics.misorientation_indices(
-230 ens_downsampled, _geo.LatticeSystem.orthorhombic, pool=pool
-231 )
-232 del enstatite, ens_resampled, ens_mean_vectors
-233
-234 _log.info("elapsed CPU time: %s", np.abs(process_time() - clock_start))
-235 _log.info("calculating ensemble averages...")
-236 olA_from_proj_XZ_mean = olA_from_proj_XZ.mean(axis=0)
-237 olA_from_proj_XZ_err = olA_from_proj_XZ.std(axis=0)
-238
-239 olA_from_proj_YX_mean = olA_from_proj_YX.mean(axis=0)
-240 olA_from_proj_YX_err = olA_from_proj_YX.std(axis=0)
-241
-242 ens_from_proj_XZ_mean = ens_from_proj_XZ.mean(axis=0)
-243 ens_from_proj_XZ_err = ens_from_proj_XZ.std(axis=0)
+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 ens_from_proj_YX_mean = ens_from_proj_YX.mean(axis=0)
-246 ens_from_proj_YX_err = ens_from_proj_YX.std(axis=0)
+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 if outdir is not None:
-249 np.savez(
-250 f"{out_basepath}.npz",
-251 olA_from_proj_XZ_mean=olA_from_proj_XZ_mean,
-252 olA_from_proj_XZ_err=olA_from_proj_XZ_err,
-253 olA_from_proj_YX_mean=olA_from_proj_YX_mean,
-254 olA_from_proj_YX_err=olA_from_proj_YX_err,
-255 ens_from_proj_XZ_mean=ens_from_proj_XZ_mean,
-256 ens_from_proj_XZ_err=ens_from_proj_XZ_err,
-257 ens_from_proj_YX_mean=ens_from_proj_YX_mean,
-258 ens_from_proj_YX_err=ens_from_proj_YX_err,
-259 )
-260 np.savez(
-261 f"{out_basepath}_strength.npz",
-262 olA_mean=olA_strength.mean(axis=0),
-263 ens_mean=ens_strength.mean(axis=0),
-264 olA_err=olA_strength.std(axis=0),
-265 ens_err=ens_strength.std(axis=0),
-266 )
+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 )