From aac2397fbacf39a86696d69d669c37dd4118463d Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Tue, 11 Jul 2023 13:41:59 +1000 Subject: [PATCH 01/30] remove: Remove outdated copy of voigt averaging method --- src/pydrex/stats.py | 42 ------------------------------------------ 1 file changed, 42 deletions(-) diff --git a/src/pydrex/stats.py b/src/pydrex/stats.py index dd9edb55..466bc720 100644 --- a/src/pydrex/stats.py +++ b/src/pydrex/stats.py @@ -2,52 +2,10 @@ import numpy as np from pydrex import geometry as _geo -from pydrex import minerals as _minerals -from pydrex import tensors as _tensors _RNG = np.random.default_rng(seed=8845) -def average_stiffness(minerals, config): - """Calculate average elastic tensor from a list of `minerals`. - - The `config` dictionary must contain volume fractions of all occurring mineral phases, - indexed by keys of the format `"_fraction"`. - - """ - n_grains = minerals[0].n_grains - assert np.all( - [m.n_grains == n_grains for m in minerals[1:]] - ), "cannot average minerals with varying grain counts" - elastic_tensors = [] - for phase in _minerals.MineralPhase: - if phase == _minerals.MineralPhase.olivine: - elastic_tensors.append( - _tensors.Voigt_to_elastic_tensor(_minerals.OLIVINE_STIFFNESS) - ) - elif phase == _minerals.MineralPhase.enstatite: - elastic_tensors.append( - _tensors.Voigt_to_elastic_tensor(_minerals.ENSTATITE_STIFFNESS) - ) - - average_tensor = np.zeros((3, 3, 3, 3)) - for n in n_grains: - for mineral in minerals: - if mineral.phase == _minerals.MineralPhase.olivine: - average_tensor += ( - _tensors.rotated_tensor(mineral.orientations[n, ...].transpose()) - * mineral.fractions[n] - * config["olivine_fraction"] - ) - elif mineral.phase == _minerals.MineralPhase.enstatite: - average_tensor += ( - _tensors.rotated_tensor(minerals.orientations[n, ...].transpose()) - * mineral.fractions[n] - * config["enstatite_fraction"] - ) - return _tensors.elastic_tensor_to_Voigt(average_tensor) - - def resample_orientations(orientations, fractions, n_samples=None, rng=_RNG): """Generate new samples from `orientations` weighed by the volume distribution. From e93bc48f7ef3f03b680ae91baaa14728c3391628 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Tue, 11 Jul 2023 15:31:40 +1000 Subject: [PATCH 02/30] test: Clean up simple shear 2d tests and compare to Kaminski 2001 --- src/pydrex/minerals.py | 4 + src/pydrex/visualisation.py | 6 +- tests/test_simple_shear_2d.py | 395 +++++++++++++--------------------- 3 files changed, 159 insertions(+), 246 deletions(-) diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index e098ff28..bd3539d7 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -13,6 +13,7 @@ from pydrex import io as _io from pydrex import logger as _log from pydrex import tensors as _tensors +from pydrex import stats as _stats OLIVINE_STIFFNESS = np.array( [ @@ -202,6 +203,8 @@ class Mineral: - `n_grains` (int) — number of grains in the aggregate - `fractions` (list of arrays) — grain volume fractions for each texture snapshot - `orientations` (list of arrays) — grain orientation matrices for each texture snapshot + - `rng` (`numpy.random.Generator`) — random number generator used for the isotropic + initial condition when `fractions_init` or `orientations_init` are not provided """ @@ -214,6 +217,7 @@ class Mineral: orientations_init: np.ndarray = None fractions: list = field(default_factory=list) orientations: list = field(default_factory=list) + rng: np.random.Generator = _stats._RNG def __str__(self): # String output, used for str(self) and f"{self}", etc. diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index b51bce1a..1c5bdf4e 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -167,21 +167,21 @@ def simple_shear_stationary_2d( _io.data("thirdparty") / "Kaminski2001_GBMshear.scsv" ) lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M0) / 200, + np.asarray(data_Kaminski2001.equivalent_strain_M0) / 200 * timestop, data_Kaminski2001.angle_M0, alpha=0.33, label=r"$M^{\ast}=0$", ) colors.append(lines[0].get_color()) lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M50) / 200, + np.asarray(data_Kaminski2001.equivalent_strain_M50) / 200 * timestop, data_Kaminski2001.angle_M50, alpha=0.33, label=r"$M^{\ast}=50$", ) colors.append(lines[0].get_color()) lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M200) / 200, + np.asarray(data_Kaminski2001.equivalent_strain_M200) / 200 * timestop, data_Kaminski2001.angle_M200, alpha=0.33, label=r"$M^{\ast}=200$", diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 7aee7ef6..5582c12d 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -1,32 +1,27 @@ -"""> PyDRex: 2D simple shear tests. - -NOTE: In scipy, rotations are represented as a matrix -that transforms [1, 0, 0] to the new a-axis vector. -We need the inverse of that rotation, -which represents the change of coordinates -from the grain-local to the global (Eulerian) frame. - -""" +"""> PyDRex: 2D simple shear tests.""" import contextlib as cl import numpy as np -from scipy.spatial.transform import Rotation -from pydrex import core as _core from pydrex import diagnostics as _diagnostics from pydrex import logger as _log from pydrex import minerals as _minerals from pydrex import stats as _stats from pydrex import visualisation as _vis +# Subdirectory of `outdir` used to store outputs from these tests. +SUBDIR = "2d_simple_shear" + + +class TestOlivineA: + """Tests for A-type olivine polycrystals in 2D simple shear.""" -class TestSinglePolycrystalOlivineA: - """Tests for a single A-type olivine polycrystal in 2D simple shear.""" + class_id = "olivineA" def get_position(self, t): return np.zeros(3) - def test_shearYZ_initQ1( + def test_dvdx_GBM( self, params_Kaminski2001_fig5_solid, # GBM = 0 params_Kaminski2001_fig5_shortdash, # GBM = 50 @@ -34,163 +29,122 @@ def test_shearYZ_initQ1( rng, outdir, ): - """Test clockwise a-axis rotation around X. + r"""Test a-axis alignment to shear in the Y direction. - Initial condition: randomised a-axis orientation within the first - quadrant (between +Y and +Z axes) and steady flow with velocity gradient + Velocity gradient: + $$\bm{L} = \begin{bmatrix} 0 & 0 & 0 \cr 2 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$ - 0 0 0 - L = 0 0 2 - 0 0 0 - - Orientations set up for slip on (010)[100]. - Tests the effect of grain boundary migration, - similar to Fig. 5 in [Kaminski 2001]. - - [Kaminski 2001]: https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9 + Tests the effect of grain boundary migration, similar to Fig. 5 in + [Kaminski 2001](https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9). """ - strain_rate_scale = 2e-5 - timescale = 1 / (strain_rate_scale / 2) - n_grains = 1000 + strain_rate = 1 def get_velocity_gradient(x): - Δv = np.zeros((3, 3)) - Δv[1, 2] = strain_rate_scale - return Δv + grad_v = np.zeros((3, 3)) + grad_v[1, 0] = 2 * strain_rate + return grad_v - # Initial orientations with a-axis in first quadrant of the YZ plane, - # and c-axis along the +X direction (important!) - # This means that the starting average is the same, - # which is convenient for comparing the texture evolution. - orientations_init = ( - Rotation.from_euler( - "zxz", - [[x * np.pi / 2, np.pi / 2, np.pi / 2] for x in rng.random(n_grains)], - ) - .inv() - .as_matrix() - ) - # Uncomment to check a-axis vectors, should be near [0, a, a]. - # assert False, f"{orientations_init[0:10, 0, :]}" - # Uncomment to check c-axis vectors, should be near [1, 0, 0]. - # assert False, f"{orientations_init[0:10, 2, :]}" + shear_direction = [0, 1, 0] + timestamps = np.linspace(0, 2e3, 500) - # One mineral to test each value of grain boundary mobility. - minerals = [ - _minerals.Mineral( - _core.MineralPhase.olivine, - _core.MineralFabric.olivine_A, - _core.DeformationRegime.dislocation, - n_grains=n_grains, - fractions_init=np.full(n_grains, 1 / n_grains), - orientations_init=orientations_init, - ) - for i in range(3) - ] - - # Optional plotting setup. + # Optional logging and plotting setup. + optional_logging = cl.nullcontext() if outdir is not None: + out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_GBM" + optional_logging = _log.logfile_enable(f"{out_basepath}.log") + θ_fse = [45] labels = [] angles = [] indices = [] - for mineral, params in zip( - minerals, - ( - params_Kaminski2001_fig5_solid, # GBM = 0 - params_Kaminski2001_fig5_shortdash, # GBM = 50 - params_Kaminski2001_fig5_longdash, # GBM = 200 - ), - ): - time = 0 - timestep = 0.025 - timestop = 1 - deformation_gradient = np.eye(3) # Undeformed initial state. - while time < timestop * timescale: - deformation_gradient = mineral.update_orientations( - params, - deformation_gradient, - get_velocity_gradient, - pathline=(time, time + timestep * timescale, self.get_position), - ) - time += timestep * timescale - - n_timesteps = len(mineral.orientations) - misorient_angles = np.zeros(n_timesteps) - misorient_indices = np.zeros(n_timesteps) - # Loop over first dimension (time steps) of orientations. - for idx, matrices in enumerate(mineral.orientations): - orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx] - ) - direction_mean = _diagnostics.bingham_average( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - ) - misorient_angles[idx] = _diagnostics.smallest_angle( - direction_mean, [0, 1, 0] - ) - misorient_indices[idx] = _diagnostics.misorientation_index( - orientations_resampled + with optional_logging: + for p, params in enumerate( + ( + params_Kaminski2001_fig5_solid, # GBM = 0 + params_Kaminski2001_fig5_shortdash, # GBM = 50 + params_Kaminski2001_fig5_longdash, # GBM = 200 + ), + ): + mineral = _minerals.Mineral( + n_grains=params["number_of_grains"], rng=rng ) + deformation_gradient = np.eye(3) # Undeformed initial state. + for t, time in enumerate(timestamps[:-1], start=1): + _log.info( + "calculating CPO at %s (t=%s) with velocity gradient: %s", + self.get_position(None), + time, + get_velocity_gradient(None).flatten(), + ) + deformation_gradient = mineral.update_orientations( + params, + deformation_gradient, + get_velocity_gradient, + pathline=(time, time + timestamps[t], self.get_position), + ) + fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) + _log.info("strain = %s", fse_λ) + if p == 0: + θ_fse.append( + _diagnostics.smallest_angle(fse_v, shear_direction) + ) - # Check for uncorrupted record of initial condition. - assert np.isclose(misorient_angles[0], 45.0, atol=5.0) - assert misorient_indices[0] < 0.71 - # Check for mostly smoothly decreasing misalignment. - angles_diff = np.diff(misorient_angles) - assert np.max(angles_diff) < 3.6 - assert np.min(angles_diff) > -7.5 - assert np.sum(angles_diff) < -24.0 - - # Check alignment and texture strength (half way and final value). - halfway = round(n_timesteps / 2) - match params["gbm_mobility"]: - case 0: - np.testing.assert_allclose( - misorient_angles, - misorient_angles[0] - * np.exp( - np.linspace(0, timestop, n_timesteps) - * (np.cos(2 * np.deg2rad(misorient_angles[0])) - 1) - ), - atol=5.2, - rtol=0.1, + n_timesteps = len(mineral.orientations) + misorient_angles = np.zeros(n_timesteps) + misorient_indices = np.zeros(n_timesteps) + # Loop over first dimension (time steps) of orientations. + for idx, matrices in enumerate(mineral.orientations): + orientations_resampled, _ = _stats.resample_orientations( + matrices, mineral.fractions[idx] + ) + direction_mean = _diagnostics.bingham_average( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], ) - assert np.isclose(misorient_angles[halfway], 25, atol=2.0) - assert np.isclose(misorient_angles[-1], 17.0, atol=1.0) - assert np.isclose(misorient_indices[halfway], 0.925, atol=0.075) - assert np.isclose(misorient_indices[-1], 0.975, atol=0.05) - case 50: - assert np.isclose(misorient_angles[halfway], 15, atol=2.0) - assert np.isclose(misorient_angles[-1], 10, atol=1.0) - assert np.isclose(misorient_indices[halfway], 0.925, atol=0.075) - assert np.isclose(misorient_indices[-1], 0.97, atol=0.03) - case 200: - assert np.isclose(misorient_angles[halfway], 9, atol=1.5) - assert np.isclose(misorient_angles[-1], 7, atol=1.0) - assert np.isclose(misorient_indices[halfway], 0.975, atol=0.05) - assert np.isclose(misorient_indices[-1], 0.99, atol=0.03) + misorient_angles[idx] = _diagnostics.smallest_angle( + direction_mean, shear_direction + ) + # misorient_indices[idx] = _diagnostics.misorientation_index( + # orientations_resampled + # ) - # Optionally store plotting metadata. - if outdir is not None: - labels.append(f"$M^∗$ = {params['gbm_mobility']}") - angles.append(misorient_angles) - indices.append(misorient_indices) + # Optionally store plotting metadata. + if outdir is not None: + labels.append(f"$M^∗$ = {params['gbm_mobility']}") + # angles.append(misorient_angles) + angles.append( + [ + # _diagnostics.smallest_angle( + # _diagnostics.anisotropy(x)[1][2, :], shear_direction + # ) + np.abs( + np.rad2deg( + np.arcsin(_diagnostics.anisotropy(x)[1][2, 2]) + ) + ) + for x in _minerals.voigt_averages([mineral], params) + ] + ) + indices.append(misorient_indices) + mineral.save( + f"{out_basepath}.npz", + postfix=f"M{params['gbm_mobility']}", + ) # Optionally plot figure. if outdir is not None: _vis.simple_shear_stationary_2d( angles, indices, - timestop=timestop, - savefile=f"{outdir}/simple_shearYZ_stationary_olivineA_initQ1.png", + timestop=timestamps[-1], + savefile=f"{out_basepath}.png", markers=("o", "v", "s"), labels=labels, + θ_fse=θ_fse, ) - def test_shearXZ_initQ1( + def test_dudz_GBS( self, params_Kaminski2004_fig4_triangles, # GBS = 0.4 params_Kaminski2004_fig4_squares, # GBS = 0.2 @@ -198,89 +152,63 @@ def test_shearXZ_initQ1( rng, outdir, ): - """Test clockwise a-axis rotation around Y. + r"""Test a-axis alignment to shear in X direction. - Initial condition: randomised a-axis orientation within the first - quadrant (between +X and +Z axes) and steady flow with velocity gradient - - 0 0 2 - L = 0 0 0 - 0 0 0 - - Orientations set up for slip on (010)[100]. + Velocity gradient: + $$\bm{L} = \begin{bmatrix} 0 & 0 & 2 \cr 0 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$ """ - strain_rate_scale = 2e-5 - timescale = 1 / (strain_rate_scale / 2) - n_grains = 1000 + strain_rate = 1 def get_velocity_gradient(x): - Δv = np.zeros((3, 3)) - Δv[0, 2] = strain_rate_scale - return Δv + grad_v = np.zeros((3, 3)) + grad_v[0, 2] = 2 * strain_rate + return grad_v - # Initial orientations with a-axis in first quadrant of the XZ plane, - # and c-axis along the -Y direction (important!) - # This means that the starting average is the same, - # which is convenient for comparing the texture evolution. - orientations_init = ( - Rotation.from_euler( - "zxz", - [[x * np.pi / 2, np.pi / 2, 0] for x in rng.random(n_grains)], - ) - .inv() - .as_matrix() - ) - # Uncomment to check a-axis vectors, should be near [a, 0, a]. - # assert False, f"{orientations_init[0:10, 0, :]}" - # Uncomment to check c-axis vectors, should be near [0, -1, 0]. - # assert False, f"{orientations_init[0:10, 2, :]}" - - # One mineral to test each grain boundary sliding threshold. - minerals = [ - _minerals.Mineral( - _core.MineralPhase.olivine, - _core.MineralFabric.olivine_A, - _core.DeformationRegime.dislocation, - n_grains=n_grains, - fractions_init=np.full(n_grains, 1 / n_grains), - orientations_init=orientations_init, - ) - for i in range(3) - ] + shear_direction = [1, 0, 0] + timestamps = np.linspace(0, 2e3, 500) # Optional plotting and logging setup. optional_logging = cl.nullcontext() if outdir is not None: - optional_logging = _log.logfile_enable( - f"{outdir}/simple_shearXZ_initQ1.log" - ) + out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dudz_GBS" + optional_logging = _log.logfile_enable(f"{out_basepath}.log") + θ_fse = [45] labels = [] angles = [] indices = [] with optional_logging: - for mineral, params in zip( - minerals, + for p, params in enumerate( ( params_Kaminski2004_fig4_triangles, # GBS = 0.4 params_Kaminski2004_fig4_squares, # GBS = 0.2 params_Kaminski2004_fig4_circles, # GBS = 0 ), ): - time = 0 - timestep = 0.025 - timestop = 1 + mineral = _minerals.Mineral( + n_grains=params["number_of_grains"], rng=rng + ) deformation_gradient = np.eye(3) # Undeformed initial state. - while time < timestop * timescale: + for t, time in enumerate(timestamps[:-1], start=1): + _log.info( + "calculating CPO at %s (t=%s) with velocity gradient: %s", + self.get_position(None), + time, + get_velocity_gradient(None).flatten(), + ) deformation_gradient = mineral.update_orientations( params, deformation_gradient, get_velocity_gradient, - pathline=(time, time + timestep * timescale, self.get_position), - # atol=1e-99, + pathline=(time, time + timestamps[t], self.get_position), ) - time += timestep * timescale + fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) + _log.info("strain = %s", fse_λ) + if p == 0: + θ_fse.append( + _diagnostics.smallest_angle(fse_v, shear_direction) + ) n_timesteps = len(mineral.orientations) misorient_angles = np.zeros(n_timesteps) @@ -297,60 +225,41 @@ def get_velocity_gradient(x): misorient_angles[idx] = _diagnostics.smallest_angle( direction_mean, [1, 0, 0] ) - misorient_indices[idx] = _diagnostics.misorientation_index( - orientations_resampled - ) - - # Check for uncorrupted record of initial condition. - assert np.isclose(misorient_angles[0], 45.0, atol=5.0) - assert misorient_indices[0] < 0.71 - # Check for mostly smoothly decreasing misalignment. - angles_diff = np.diff(misorient_angles) - assert np.max(angles_diff) < 3.6 - assert np.min(angles_diff) > -7.5 - assert np.sum(angles_diff) < -24.0 - # Check that recrystallization is causing faster alignment. - np.testing.assert_array_less( - misorient_angles - 3.8, # Tolerance for GBM onset latency. - misorient_angles[0] - * np.exp( - np.linspace(0, timestop, n_timesteps) - * (np.cos(2 * np.deg2rad(misorient_angles[0])) - 1) - ), - ) - - # Check alignment and texture strength (half way and final value). - halfway = round(n_timesteps / 2) - match params["gbs_threshold"]: - case 0: - assert np.isclose(misorient_angles[halfway], 11, atol=1.5) - assert np.isclose(misorient_angles[-1], 8, atol=1.25) - assert np.isclose(misorient_indices[halfway], 0.975, atol=0.075) - assert np.isclose(misorient_indices[-1], 0.99, atol=0.03) - case 0.2: - assert np.isclose(misorient_angles[halfway], 13, atol=2.075) - assert np.isclose(misorient_angles[-1], 11, atol=1.5) - assert np.isclose(misorient_indices[halfway], 0.755, atol=0.08) - assert np.isclose(misorient_indices[-1], 0.755, atol=0.075) - case 0.4: - assert np.isclose(misorient_angles[halfway], 19, atol=2.0) - assert np.isclose(misorient_angles[-1], 16, atol=2.5) - assert misorient_indices[halfway] < 0.75 - assert misorient_indices[-1] < 0.71 + # misorient_indices[idx] = _diagnostics.misorientation_index( + # orientations_resampled + # ) # Optionally store plotting metadata. if outdir is not None: labels.append(f"$f_{{gbs}}$ = {params['gbs_threshold']}") - angles.append(misorient_angles) + # angles.append(misorient_angles) + angles.append( + [ + # _diagnostics.smallest_angle( + # _diagnostics.anisotropy(x)[1][2, :], shear_direction + # ) + np.abs( + np.rad2deg( + np.arcsin(_diagnostics.anisotropy(x)[1][2, 2]) + ) + ) + for x in _minerals.voigt_averages([mineral], params) + ] + ) indices.append(misorient_indices) + mineral.save( + f"{out_basepath}.npz", + postfix=f"X{params['gbs_threshold']}", + ) # Optionally plot figure. if outdir is not None: _vis.simple_shear_stationary_2d( angles, indices, - timestop=timestop, - savefile=f"{outdir}/simple_shearXZ_stationary_olivineA_initQ1.png", + timestop=timestamps[-1], + savefile=f"{out_basepath}.png", markers=("o", "v", "s"), labels=labels, + θ_fse=θ_fse, ) From 0fd741648235e6fccd314284b5ff5438dd128022 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Tue, 11 Jul 2023 16:30:56 +1000 Subject: [PATCH 03/30] test: Fix bad conditional --- tests/test_simple_shear_2d.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 5582c12d..825151fb 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -85,7 +85,7 @@ def get_velocity_gradient(x): ) fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) _log.info("strain = %s", fse_λ) - if p == 0: + if p == 0 and outdir is not None: θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) @@ -205,7 +205,7 @@ def get_velocity_gradient(x): ) fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) _log.info("strain = %s", fse_λ) - if p == 0: + if p == 0 and outdir is not None: θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) From 71497a20f371528b5e63a0b2d00ea07cd0834e0e Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Wed, 12 Jul 2023 03:00:54 +1000 Subject: [PATCH 04/30] fix: Fix unused rng attribute --- src/pydrex/minerals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index bd3539d7..636aff28 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -253,7 +253,7 @@ def __post_init__(self): self.fractions_init = np.full(self.n_grains, 1.0 / self.n_grains) if self.orientations_init is None: self.orientations_init = Rotation.random( - self.n_grains, random_state=1 + self.n_grains, random_state=self.rng ).as_matrix() # Copy the initial values to the storage lists. From f6e66c5450f4cbe70696403a843cab534e9ad96e Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Wed, 12 Jul 2023 03:02:42 +1000 Subject: [PATCH 05/30] test: Use shorter time for CPO integration --- tests/test_simple_shear_2d.py | 68 +++++++++++++++++------------------ 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 825151fb..6ea9a2fd 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -46,7 +46,7 @@ def get_velocity_gradient(x): return grad_v shear_direction = [0, 1, 0] - timestamps = np.linspace(0, 2e3, 500) + timestamps = np.linspace(0, 1, 500) # Optional logging and plotting setup. optional_logging = cl.nullcontext() @@ -91,23 +91,23 @@ def get_velocity_gradient(x): ) n_timesteps = len(mineral.orientations) - misorient_angles = np.zeros(n_timesteps) misorient_indices = np.zeros(n_timesteps) - # Loop over first dimension (time steps) of orientations. - for idx, matrices in enumerate(mineral.orientations): - orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx] - ) - direction_mean = _diagnostics.bingham_average( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - ) - misorient_angles[idx] = _diagnostics.smallest_angle( - direction_mean, shear_direction - ) - # misorient_indices[idx] = _diagnostics.misorientation_index( - # orientations_resampled - # ) + # misorient_angles = np.zeros(n_timesteps) + # # Loop over first dimension (time steps) of orientations. + # for idx, matrices in enumerate(mineral.orientations): + # orientations_resampled, _ = _stats.resample_orientations( + # matrices, mineral.fractions[idx] + # ) + # direction_mean = _diagnostics.bingham_average( + # orientations_resampled, + # axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + # ) + # misorient_angles[idx] = _diagnostics.smallest_angle( + # direction_mean, shear_direction + # ) + # misorient_indices[idx] = _diagnostics.misorientation_index( + # orientations_resampled + # ) # Optionally store plotting metadata. if outdir is not None: @@ -166,7 +166,7 @@ def get_velocity_gradient(x): return grad_v shear_direction = [1, 0, 0] - timestamps = np.linspace(0, 2e3, 500) + timestamps = np.linspace(0, 1, 500) # Optional plotting and logging setup. optional_logging = cl.nullcontext() @@ -211,23 +211,23 @@ def get_velocity_gradient(x): ) n_timesteps = len(mineral.orientations) - misorient_angles = np.zeros(n_timesteps) misorient_indices = np.zeros(n_timesteps) - # Loop over first dimension (time steps) of orientations. - for idx, matrices in enumerate(mineral.orientations): - orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx] - ) - direction_mean = _diagnostics.bingham_average( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - ) - misorient_angles[idx] = _diagnostics.smallest_angle( - direction_mean, [1, 0, 0] - ) - # misorient_indices[idx] = _diagnostics.misorientation_index( - # orientations_resampled - # ) + # misorient_angles = np.zeros(n_timesteps) + # # Loop over first dimension (time steps) of orientations. + # for idx, matrices in enumerate(mineral.orientations): + # orientations_resampled, _ = _stats.resample_orientations( + # matrices, mineral.fractions[idx] + # ) + # direction_mean = _diagnostics.bingham_average( + # orientations_resampled, + # axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + # ) + # misorient_angles[idx] = _diagnostics.smallest_angle( + # direction_mean, [1, 0, 0] + # ) + # misorient_indices[idx] = _diagnostics.misorientation_index( + # orientations_resampled + # ) # Optionally store plotting metadata. if outdir is not None: From 3bc31ff9612fc27527d592796157d7e69da2d0a1 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Wed, 12 Jul 2023 03:44:41 +1000 Subject: [PATCH 06/30] fix: Add missing datafile --- .../thirdparty/Kaminski2001_GBMshear.scsv | 125 ++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv diff --git a/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv b/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv new file mode 100644 index 00000000..769d4df9 --- /dev/null +++ b/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv @@ -0,0 +1,125 @@ +--- +# Digitized data from Figure 5 in Kaminski & Ribe 2001, https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9 +# This data records the mean angle of the a-axes of an olivine aggregate from the horizontal shear direction. +# The flow field is a static simple shear and the aggregate starts with all grain orientations restricted to the first quadrant +# in the 2D deformation plane (so the mean angle always starts at 45°). + +schema: + delimiter: ',' + missing: '-' + fields: + - name: equivalent_strain_M0 + type: float + fill: NaN + unit: "200 x max_strain x time" + - name: angle_M0 + type: float + fill: NaN + unit: degree + - name: equivalent_strain_M50 + type: float + fill: NaN + unit: "200 x max_strain x time" + - name: angle_M50 + type: float + fill: NaN + unit: degree + - name: equivalent_strain_M200 + type: float + fill: NaN + unit: "200 x max_strain x time" + - name: angle_M200 + type: float + fill: NaN + unit: degree +--- +equivalent_strain_M0, angle_M0,equivalent_strain_M50, angle_M50,equivalent_strain_M200, angle_M200 + 0.3076923076923084, 44.93959731543627, 1.098613804496157, 44.77430481378372, 1.5111685699920994, 44.72085812156877 + 1.9237233354880399,44.718921206880175, 3.367665014723837, 44.37236387335087, 3.7802197802197792, 44.29219383502844 + 4.19277454571572,44.382124404232606, 5.5748330101271275, 44.00933456156514, 5.911752735282144, 43.91091652503244 + 6.4618257559434,44.026404780609305, 7.390073978309269, 43.76427140580321, 7.390073978309269, 43.64428903552475 + 8.730876966171083, 43.71627845769182, 10.381096028154849, 43.30748138181449, 10.277957336780863, 43.0062852881645 + 10.999928176398763, 43.3868987470251, 13.243194713782946, 42.80741200283247, 13.475256769374413, 42.10306095185686 + 13.330862601450834,43.052174367136885, 15.385908209437625, 42.32208498611974, 15.744307979602093, 41.497695356360985 + 15.538030596854123, 42.7215155332013, 16.98197227608992, 42.024527036765505, 17.291388350211875, 41.01667512642642 + 17.78645406880701,42.402471503142664, 19.972994325935495, 41.529599759366846, 24.40795805501687, 38.67374665953426 + 20.07613301730948,42.073610733697606, 23.17029375852905, 40.951884646476046, 26.573870573870565, 37.89593367548361 + 22.34518422753716, 41.76929181271859, 26.47073188249658, 40.37596926913943, 27.811534870358393, 37.43520137361432 + 24.61423543776484, 41.4584283988153, 29.564892623716148, 39.7928549495861, 34.30927242691948, 35.035008593634714 + 26.88328664799252,41.134475999063454, 32.65905336493571, 39.21693957224949, 36.47518494577317, 34.24606997161278 + 29.1523378582202, 40.80725135284947, 35.85635279752926, 38.62302683937109, 37.60971055088701, 33.81773290971867 + 31.421389068447887, 40.4898434460219, 39.05365223012281, 38.00031833762587, 43.79803203332614, 31.350095494324965 + 33.69044027867557, 40.16916329273219, 42.147812971342375, 37.417204018072546, 45.55138978668389, 30.569210234429306 + 35.95949148890325, 39.85502763236676, 48.542411836529475, 36.157389130148694, 46.685915391797735, 30.050286482974958 + 38.22854269913093, 39.54089197200133, 51.636572577749035, 35.545479041728534, 52.15226603461896, 27.35218293133803 + 40.497593909358606, 39.22021181871162, 54.7307333189686, 34.89037530000813, 53.69934640522875, 26.551943372427214 + 42.766645119586286, 38.91262065127047, 57.721755368814186, 34.2766654760338, 54.83387201034259, 25.970885893507234 + 45.035696329813966,38.611573976753604, 60.712777418659755, 33.63775935430098, 59.99080657904186, 23.145301073449442 + 47.304747540041646,38.294166069926035, 63.80693815987932, 32.91786513263021, 61.537886949651636, 22.319993769319733 + 49.573798750269326, 37.98984714894703, 66.7979602097249, 32.24596385907082, 62.5692738633915, 21.759504696733202 + 51.842849960497006, 37.68552822796802, 69.78898225957047, 31.53566822702232, 67.9030176173033, 18.90563831796692 + 54.11190117072469, 37.38120930698901, 72.78000430941606, 30.860167482354576, 68.49974861739567, 18.304012432713485 + 56.40387209014658, 37.08125338129285, 75.77102635926164, 30.16067026363114, 69.99525964231844, 17.81551278229403 + 58.65000359118005, 36.78893269734169, 78.65890971773322, 29.46237286861049, 75.77102635926164, 14.781672848110055 + 61.00156575450691, 36.47414258768383, 81.75307045895279, 28.728080762506302, 77.31810672987142, 14.002644458230606 + 63.18810601163541, 36.19665608769437, 84.64095381742439, 28.022584425268942, 78.45263233498525, 13.396733488324372 + 65.45715722186308, 35.90869839902606, 87.52883717589599, 27.302690203598168, 83.81584428643251, 10.6422810735031 + 67.72620843209077, 35.61092397097133, 90.41672053436757, 26.58279598192739, 85.36292465704229, 9.827257972540117 + 69.99525964231844,35.316421789378744, 93.30460389283917, 25.891697529123448, 86.49745026215612, 9.24534347668957 + 72.26431085254613, 35.03828084009685, 96.29562594268475, 25.11601150527319, 91.86066221360338, 6.348626251395267 + 74.5333620627738, 34.75359539789068, 99.28664799253032, 24.387118605831528, 93.40774258421315, 5.530389336942683 + 76.80241327300149, 34.46890995568451, 102.07139265962793, 23.6402283508481, 94.439129497953, 5.00396668734593 + 79.07146448322918,34.177680020554064, 105.06241470947352, 22.868141798106194, 100.52431228901816, 2.361627669167966 + 81.34051569345685,33.909355810658596, 107.9502980679451, 22.097854980918466, 102.79336349924583, 1.6744559121185958 + 83.60956690368454, 33.64103160076312, 110.79692594986709, 21.33476710594745, 104.13416648710763, 1.3505035123667497 + 85.87861811391221, 33.36289065148123, 113.62292609351431, 20.56268055320554, 113.7260647848883, 0.08850712680145989 + 88.1476693241399, 33.07820520927506, 116.40767076061192, 19.794193471571987, 115.99511599511597, 0.06560140156648941 + 90.41672053436757,32.816425492303864, 119.08927673633553, 19.04550348103438, 117.74847374847371, 0.02469832078973866 + 92.68577174459526, 32.55137352887054, 121.97716009480712, 18.232023010546406, 120.7394957983193,-0.08928493097479873 + 94.99607843137254, 32.28075874645157, 124.76190476190473, 17.44913804447944, 123.00854700854697,-0.08928493097479873 + 97.22387416505062, 32.01799735554174, 127.54664942900234, 16.641056780653997, 125.27759821877467,-0.08928493097479873 + 99.49292537527829, 31.75948988503269, 130.12511671335199, 15.901365467887274, 127.54664942900234,-0.10564616328550613 + 101.76197658550598,31.520615893296476, 133.01300007182357, 15.037492401882346, 129.81570063923002, -0.1612743531418772 + 104.03102779573365,31.232658204628166, 135.79774473892118, 14.207814311406779, 132.08475184945772, -0.1449131208311769 + 106.30007900596134,30.987239719967675, 138.4793507146448, 13.41593066756893, 133.32241614594554,-0.08928493097479873 + 108.56913021618902,30.725460002996485, 141.16095669036844, 12.559256543780705, 144.255117431588, -0.1612743531418772 + 110.8381814264167,30.473497025411714, 143.73942397471805, 11.728978541453742, 146.52416864181566,-0.14818536729332266 + 113.0247216835452,30.236259156906574, 146.4210299504417, 10.89810062727539, 148.79321985204336,-0.15472986021759993 + 115.37628384687206, 29.96302657731789, 149.1026359261653, 10.097218305666658, 151.06227106227104,-0.08274043805052145 + 117.64533505709974, 29.74051381789238, 151.88738059326292, 9.283737835178684, 153.3313222724987,-0.08928493097479873 + 119.91438626732742,29.495095333231887, 154.67212526036053, 8.48345542542134, 155.60037348272638,-0.08928493097479873 + 122.18343747755509,29.266038080882094, 157.45686992745814, 7.692771605286275, 157.86942469295408,-0.08928493097479873 + 124.45248868778279,29.017347349759465, 160.24161459455573, 6.920685052544364, 160.13847590318176,-0.08928493097479873 + 126.72153989801046,28.807923576182514, 163.12949795302734, 6.119802730935625, 162.38460740421522,-0.08928493097479873 + 128.99059110823814,28.549416105673462, 165.91424262012495, 5.381911153723081, 164.6765783236371,-0.08928493097479873 + 131.2596423184658,28.326903346247953, 168.9052646699705, 4.608024865426998, 166.9456295338648,-0.08928493097479873 + 133.5286935286935, 28.09784609389816, 171.69000933706812, 3.9277248259481183, 169.21468074409248,-0.08928493097479873 + 135.79774473892118,27.878605580934785, 174.7841700782877, 3.222228488710762, 171.48373195432015,-0.08928493097479873 + 138.06679594914885,27.646276082122853, 177.87833081950726, 2.567124746990352, 173.75278316454782,-0.08928493097479873 + 140.33584715937656, 27.4204910762352, 180.9724915607268, 1.9696125430036133, 176.02183437477552,-0.09910167036122175 + 142.60489836960423, 27.21106730265825, 184.16979099332036, 1.3685008679085158, 178.2908855850032,-0.08928493097479873 + 144.8739495798319, 26.98201005030846, 187.26395173453994, 0.8735735905098565, 180.55993679523087,-0.08928493097479873 + 147.14300079005957,26.779130769655787, 191.3158288956608, 0.2655200782772269, 182.82898800545857,-0.08928493097479873 + 149.41205200028728,26.569706996078832, 194.20002872944045, 0.06969170964415383, 185.09803921568624,-0.08928493097479873 + 151.57796451914095,26.381225599859572, 197.40592305298182,-0.02796060838803527, 187.36709042591391,-0.10891840974763767 + 154.1564318034906, 26.12795372368995, 199.1936603701309, -0.0652884569191059, 189.6361416361416,-0.11546290267192205 + 156.42548301371826, 25.94143567534797, -, -, 191.9051928463693,-0.08928493097479873 + 158.69453422394596, 25.72546740884674, -, -, 194.17424405659696,-0.05329021989126659 + 160.96358543417364,25.525860374656208, -, -, 196.44329526682463,-0.07292369866409842 + 163.2326366444013,25.336070079852096, -, -, -, - + 165.501687854629,25.136463045661564, -, -, -, - + 167.77073906485668, 24.95321724378173, -, -, -, - + 170.03979027508436, 24.77651593482618, -, -, -, - + 172.30884148531203,24.576908900635644, -, -, -, - + 174.57789269553973,24.396935345217948, -, -, -, - + 176.8469439057674,24.203872803951693, -, -, -, - + 179.0128564246211,24.030770966104498, -, -, -, - + 181.38504632622275, 23.86028692542701, -, -, -, - + 183.65409753645045,23.686857862933593, -, -, -, - + 185.92314874667812, 23.50361206105376, -, -, -, - + 188.1921999569058,23.349816477333185, -, -, -, - + 190.46125116713347, 23.17311516837763, -, -, -, - + 192.73030237736117,23.002958352346358, -, -, -, - + 194.99935358758884,22.836073782777223, -, -, -, - + 197.26840479781652, 22.66591696674595, -, -, -, - + 198.8154851684263,22.551388340571055, -, -, -, - From eb9a65cf591c09c7ebdcb310945210eac6395f42 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Wed, 12 Jul 2023 06:14:41 +1000 Subject: [PATCH 07/30] test: Use even shorter time for CPO integration --- tests/test_simple_shear_2d.py | 64 +++++++++++++++++------------------ 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 6ea9a2fd..576cfc82 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -46,7 +46,7 @@ def get_velocity_gradient(x): return grad_v shear_direction = [0, 1, 0] - timestamps = np.linspace(0, 1, 500) + timestamps = np.linspace(0, 0.005, 500) # Optional logging and plotting setup. optional_logging = cl.nullcontext() @@ -92,40 +92,40 @@ def get_velocity_gradient(x): n_timesteps = len(mineral.orientations) misorient_indices = np.zeros(n_timesteps) - # misorient_angles = np.zeros(n_timesteps) - # # Loop over first dimension (time steps) of orientations. - # for idx, matrices in enumerate(mineral.orientations): - # orientations_resampled, _ = _stats.resample_orientations( - # matrices, mineral.fractions[idx] - # ) - # direction_mean = _diagnostics.bingham_average( - # orientations_resampled, - # axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - # ) - # misorient_angles[idx] = _diagnostics.smallest_angle( - # direction_mean, shear_direction - # ) - # misorient_indices[idx] = _diagnostics.misorientation_index( - # orientations_resampled - # ) + misorient_angles = np.zeros(n_timesteps) + # Loop over first dimension (time steps) of orientations. + for idx, matrices in enumerate(mineral.orientations): + orientations_resampled, _ = _stats.resample_orientations( + matrices, mineral.fractions[idx] + ) + direction_mean = _diagnostics.bingham_average( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + ) + misorient_angles[idx] = _diagnostics.smallest_angle( + direction_mean, shear_direction + ) + misorient_indices[idx] = _diagnostics.misorientation_index( + orientations_resampled + ) # Optionally store plotting metadata. if outdir is not None: labels.append(f"$M^∗$ = {params['gbm_mobility']}") - # angles.append(misorient_angles) - angles.append( - [ - # _diagnostics.smallest_angle( - # _diagnostics.anisotropy(x)[1][2, :], shear_direction - # ) - np.abs( - np.rad2deg( - np.arcsin(_diagnostics.anisotropy(x)[1][2, 2]) - ) - ) - for x in _minerals.voigt_averages([mineral], params) - ] - ) + angles.append(misorient_angles) + # angles.append( + # [ + # # _diagnostics.smallest_angle( + # # _diagnostics.anisotropy(x)[1][2, :], shear_direction + # # ) + # np.abs( + # np.rad2deg( + # np.arcsin(_diagnostics.anisotropy(x)[1][2, 2]) + # ) + # ) + # for x in _minerals.voigt_averages([mineral], params) + # ] + # ) indices.append(misorient_indices) mineral.save( f"{out_basepath}.npz", @@ -166,7 +166,7 @@ def get_velocity_gradient(x): return grad_v shear_direction = [1, 0, 0] - timestamps = np.linspace(0, 1, 500) + timestamps = np.linspace(0, 0.005, 500) # Optional plotting and logging setup. optional_logging = cl.nullcontext() From 0f1f170de1ab0009f7c0bc4da22c9c981a75f06e Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Wed, 12 Jul 2023 10:55:02 +1000 Subject: [PATCH 08/30] dev: Add CLI tool to plot polefigures from npz file --- tools/pfplot.py | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 tools/pfplot.py diff --git a/tools/pfplot.py b/tools/pfplot.py new file mode 100644 index 00000000..90e4edd3 --- /dev/null +++ b/tools/pfplot.py @@ -0,0 +1,10 @@ +import sys +from pydrex import visualisation as _vis + +_vis.polefigures( + sys.argv[1], + postfix=sys.argv[2], + i_range=range(*[int(s) for s in sys.argv[3].split(":")]), + density=bool(int(sys.argv[4])), + ref_axes=sys.argv[5], +) From 7091984a2dcd686cbe212c1bf01dcbe4f2852ed7 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Wed, 12 Jul 2023 16:20:24 +1000 Subject: [PATCH 09/30] test: Use Fraters & Billen, 2021 parameters for CPO integration --- tests/test_simple_shear_2d.py | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 576cfc82..a1d203f8 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -38,7 +38,7 @@ def test_dvdx_GBM( [Kaminski 2001](https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9). """ - strain_rate = 1 + strain_rate = 1.5844e-14 # Strain rate from Fraters & Billen, 2021. def get_velocity_gradient(x): grad_v = np.zeros((3, 3)) @@ -46,7 +46,9 @@ def get_velocity_gradient(x): return grad_v shear_direction = [0, 1, 0] - timestamps = np.linspace(0, 0.005, 500) + # Solve texture until 20 Ma. + timestamps = np.linspace(0, 20 * 1e6 * 31556952, 100) + n_timesteps = len(timestamps) # Optional logging and plotting setup. optional_logging = cl.nullcontext() @@ -72,8 +74,9 @@ def get_velocity_gradient(x): deformation_gradient = np.eye(3) # Undeformed initial state. for t, time in enumerate(timestamps[:-1], start=1): _log.info( - "calculating CPO at %s (t=%s) with velocity gradient: %s", - self.get_position(None), + "step %s/%s (t=%s) with velocity gradient: %s", + t, + n_timesteps, time, get_velocity_gradient(None).flatten(), ) @@ -84,13 +87,12 @@ def get_velocity_gradient(x): pathline=(time, time + timestamps[t], self.get_position), ) fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) - _log.info("strain = %s", fse_λ) + _log.info("strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) if p == 0 and outdir is not None: θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) - n_timesteps = len(mineral.orientations) misorient_indices = np.zeros(n_timesteps) misorient_angles = np.zeros(n_timesteps) # Loop over first dimension (time steps) of orientations. @@ -105,9 +107,9 @@ def get_velocity_gradient(x): misorient_angles[idx] = _diagnostics.smallest_angle( direction_mean, shear_direction ) - misorient_indices[idx] = _diagnostics.misorientation_index( - orientations_resampled - ) + # misorient_indices[idx] = _diagnostics.misorientation_index( + # orientations_resampled + # ) # Optionally store plotting metadata. if outdir is not None: @@ -158,7 +160,7 @@ def test_dudz_GBS( $$\bm{L} = \begin{bmatrix} 0 & 0 & 2 \cr 0 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$ """ - strain_rate = 1 + strain_rate = 1.5844e-14 # Strain rate from Fraters & Billen, 2021. def get_velocity_gradient(x): grad_v = np.zeros((3, 3)) @@ -166,7 +168,9 @@ def get_velocity_gradient(x): return grad_v shear_direction = [1, 0, 0] - timestamps = np.linspace(0, 0.005, 500) + # Solve texture until 20 Ma. + timestamps = np.linspace(0, 20 * 1e6 * 31556952, 100) + n_timesteps = len(timestamps) # Optional plotting and logging setup. optional_logging = cl.nullcontext() @@ -192,8 +196,9 @@ def get_velocity_gradient(x): deformation_gradient = np.eye(3) # Undeformed initial state. for t, time in enumerate(timestamps[:-1], start=1): _log.info( - "calculating CPO at %s (t=%s) with velocity gradient: %s", - self.get_position(None), + "step %s/%s (t=%s) with velocity gradient: %s", + t, + n_timesteps, time, get_velocity_gradient(None).flatten(), ) @@ -204,13 +209,12 @@ def get_velocity_gradient(x): pathline=(time, time + timestamps[t], self.get_position), ) fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) - _log.info("strain = %s", fse_λ) + _log.info("strain √λ-1=%s", fse_λ) if p == 0 and outdir is not None: θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) - n_timesteps = len(mineral.orientations) misorient_indices = np.zeros(n_timesteps) # misorient_angles = np.zeros(n_timesteps) # # Loop over first dimension (time steps) of orientations. From 42bf59ab8e5154dc687769bf9b1703118dc1ea57 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 17 Jul 2023 17:03:40 +1000 Subject: [PATCH 10/30] test: Fix 2d simple shear tests! Much reproducibility! Such convergence! Very wow! --- src/pydrex/utils.py | 7 +++ src/pydrex/visualisation.py | 13 +++- tests/test_simple_shear_2d.py | 114 ++++++++++++++++------------------ 3 files changed, 73 insertions(+), 61 deletions(-) create mode 100644 src/pydrex/utils.py diff --git a/src/pydrex/utils.py b/src/pydrex/utils.py new file mode 100644 index 00000000..b3332b8b --- /dev/null +++ b/src/pydrex/utils.py @@ -0,0 +1,7 @@ +"""> PyDRex: Miscellaneous utility methods.""" +import numpy as np + + +def angle_fse_simpleshear(strain): + """Get angle of FSE long axis anticlockwise from the X axis in simple shear.""" + return np.rad2deg(np.arctan(np.sqrt(strain**2 + 1) + strain)) diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index 1c5bdf4e..98d6d226 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -212,7 +212,18 @@ def simple_shear_stationary_2d( ) if θ_fse is not None: - ax_mean.plot(timestamps, θ_fse, "r--", label="FSE") + ax_mean.plot(timestamps, θ_fse, "r--", alpha=0.33, label="calc. FSE") + ax_mean.plot( + timestamps, + [ + 90 - np.rad2deg( + np.arctan(np.sqrt(5e-6**2 * t**2 + 1) + 5e-6 * t) + ) + for t in timestamps + ], + "k", + label="true FSE", + ) if labels is not None: ax_mean.legend() diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index a1d203f8..595bf9db 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -8,6 +8,7 @@ from pydrex import minerals as _minerals from pydrex import stats as _stats from pydrex import visualisation as _vis +from pydrex import utils as _utils # Subdirectory of `outdir` used to store outputs from these tests. SUBDIR = "2d_simple_shear" @@ -38,17 +39,23 @@ def test_dvdx_GBM( [Kaminski 2001](https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9). """ - strain_rate = 1.5844e-14 # Strain rate from Fraters & Billen, 2021. + strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. + timestamps = np.linspace(0, 2e5, 501) # Solve until D₀t=1 ('shear strain' γ=2). + n_timesteps = len(timestamps) - 1 def get_velocity_gradient(x): + # It is independent of time or position in this test. grad_v = np.zeros((3, 3)) grad_v[1, 0] = 2 * strain_rate return grad_v - shear_direction = [0, 1, 0] - # Solve texture until 20 Ma. - timestamps = np.linspace(0, 20 * 1e6 * 31556952, 100) - n_timesteps = len(timestamps) + shear_direction = [0, 1, 0] # Used to calculate the angular diagnostics. + + # Theoretical FSE axis for simple shear. + # We want the angle from the Y axis (shear direction), so subtract from 90. + θ_fse_eq = [ + 90 - _utils.angle_fse_simpleshear(t * strain_rate) for t in timestamps + ] # Optional logging and plotting setup. optional_logging = cl.nullcontext() @@ -84,11 +91,12 @@ def get_velocity_gradient(x): params, deformation_gradient, get_velocity_gradient, - pathline=(time, time + timestamps[t], self.get_position), + pathline=(time, timestamps[t], self.get_position), ) fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) - _log.info("strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) + _log.info("› strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) if p == 0 and outdir is not None: + # θ_fse.append(90 - np.rad2deg(np.arctan2(fse_v[1], fse_v[0]))) θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) @@ -107,33 +115,23 @@ def get_velocity_gradient(x): misorient_angles[idx] = _diagnostics.smallest_angle( direction_mean, shear_direction ) - # misorient_indices[idx] = _diagnostics.misorientation_index( - # orientations_resampled - # ) + misorient_indices[idx] = _diagnostics.misorientation_index( + orientations_resampled + ) # Optionally store plotting metadata. if outdir is not None: labels.append(f"$M^∗$ = {params['gbm_mobility']}") angles.append(misorient_angles) - # angles.append( - # [ - # # _diagnostics.smallest_angle( - # # _diagnostics.anisotropy(x)[1][2, :], shear_direction - # # ) - # np.abs( - # np.rad2deg( - # np.arcsin(_diagnostics.anisotropy(x)[1][2, 2]) - # ) - # ) - # for x in _minerals.voigt_averages([mineral], params) - # ] - # ) indices.append(misorient_indices) mineral.save( f"{out_basepath}.npz", postfix=f"M{params['gbm_mobility']}", ) + # Check that FSE is correct. + np.testing.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-6, atol=0) + # Optionally plot figure. if outdir is not None: _vis.simple_shear_stationary_2d( @@ -160,17 +158,23 @@ def test_dudz_GBS( $$\bm{L} = \begin{bmatrix} 0 & 0 & 2 \cr 0 & 0 & 0 \cr 0 & 0 & 0 \end{bmatrix}$$ """ - strain_rate = 1.5844e-14 # Strain rate from Fraters & Billen, 2021. + strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. + timestamps = np.linspace(0, 2e5, 501) # Solve until D₀t=1 ('shear strain' γ=2). + n_timesteps = len(timestamps) - 1 def get_velocity_gradient(x): + # It is independent of time or position in this test. grad_v = np.zeros((3, 3)) grad_v[0, 2] = 2 * strain_rate return grad_v - shear_direction = [1, 0, 0] - # Solve texture until 20 Ma. - timestamps = np.linspace(0, 20 * 1e6 * 31556952, 100) - n_timesteps = len(timestamps) + shear_direction = [1, 0, 0] # Used to calculate the angular diagnostics. + + # Theoretical FSE axis for simple shear. + # We want the angle from the Y axis (shear direction), so subtract from 90. + θ_fse_eq = [ + 90 - _utils.angle_fse_simpleshear(t * strain_rate) for t in timestamps + ] # Optional plotting and logging setup. optional_logging = cl.nullcontext() @@ -206,56 +210,46 @@ def get_velocity_gradient(x): params, deformation_gradient, get_velocity_gradient, - pathline=(time, time + timestamps[t], self.get_position), + pathline=(time, timestamps[t], self.get_position), ) fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) - _log.info("strain √λ-1=%s", fse_λ) + _log.info("› strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) if p == 0 and outdir is not None: θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) misorient_indices = np.zeros(n_timesteps) - # misorient_angles = np.zeros(n_timesteps) - # # Loop over first dimension (time steps) of orientations. - # for idx, matrices in enumerate(mineral.orientations): - # orientations_resampled, _ = _stats.resample_orientations( - # matrices, mineral.fractions[idx] - # ) - # direction_mean = _diagnostics.bingham_average( - # orientations_resampled, - # axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - # ) - # misorient_angles[idx] = _diagnostics.smallest_angle( - # direction_mean, [1, 0, 0] - # ) - # misorient_indices[idx] = _diagnostics.misorientation_index( - # orientations_resampled - # ) + misorient_angles = np.zeros(n_timesteps) + # Loop over first dimension (time steps) of orientations. + for idx, matrices in enumerate(mineral.orientations): + orientations_resampled, _ = _stats.resample_orientations( + matrices, mineral.fractions[idx] + ) + direction_mean = _diagnostics.bingham_average( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + ) + misorient_angles[idx] = _diagnostics.smallest_angle( + direction_mean, [1, 0, 0] + ) + misorient_indices[idx] = _diagnostics.misorientation_index( + orientations_resampled + ) # Optionally store plotting metadata. if outdir is not None: labels.append(f"$f_{{gbs}}$ = {params['gbs_threshold']}") - # angles.append(misorient_angles) - angles.append( - [ - # _diagnostics.smallest_angle( - # _diagnostics.anisotropy(x)[1][2, :], shear_direction - # ) - np.abs( - np.rad2deg( - np.arcsin(_diagnostics.anisotropy(x)[1][2, 2]) - ) - ) - for x in _minerals.voigt_averages([mineral], params) - ] - ) + angles.append(misorient_angles) indices.append(misorient_indices) mineral.save( f"{out_basepath}.npz", postfix=f"X{params['gbs_threshold']}", ) + # Check that FSE is correct. + np.testing.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-6, atol=0) + # Optionally plot figure. if outdir is not None: _vis.simple_shear_stationary_2d( From 43b5fea842ab7b185822cf470870b98a58e969ae Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Tue, 18 Jul 2023 12:39:58 +1000 Subject: [PATCH 11/30] test: Add asserts for GBM test angles --- src/pydrex/utils.py | 6 ++ src/pydrex/visualisation.py | 72 +++++-------------- tests/test_simple_shear_2d.py | 126 ++++++++++++++++++++++++---------- 3 files changed, 114 insertions(+), 90 deletions(-) diff --git a/src/pydrex/utils.py b/src/pydrex/utils.py index b3332b8b..043d83af 100644 --- a/src/pydrex/utils.py +++ b/src/pydrex/utils.py @@ -2,6 +2,12 @@ import numpy as np +def skip_nans(a): + """Skip NaN values in array.""" + a = np.asarray(a) + return a[~np.isnan(a)] + + def angle_fse_simpleshear(strain): """Get angle of FSE long axis anticlockwise from the X axis in simple shear.""" return np.rad2deg(np.arctan(np.sqrt(strain**2 + 1) + strain)) diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index 98d6d226..9538c8cc 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -143,90 +143,52 @@ def _get_marker_and_label(data, seq_index, markers, labels=None): @check_marker_seq def simple_shear_stationary_2d( + strains, + target_angles, angles, indices, - timestop, savefile="pydrex_simple_shear_stationary_2d.png", markers=("."), - labels=None, θ_fse=None, + labels=None, ): """Plot diagnostics for stationary A-type olivine 2D simple shear box tests.""" fig = plt.figure(figsize=(5, 8), dpi=300) grid = fig.add_gridspec(2, 1, hspace=0.05) ax_mean = fig.add_subplot(grid[0]) ax_mean.set_ylabel("Mean angle ∈ [0, 90]°") - ax_mean.axhline(0, color=plt.rcParams["axes.edgecolor"]) ax_mean.tick_params(labelbottom=False) + ax_mean.set_xlim((strains[0], strains[-1])) + ax_mean.set_ylim((0, 60)) ax_strength = fig.add_subplot(grid[1], sharex=ax_mean) + ax_strength.set_xlim((strains[0], strains[-1])) + ax_strength.set_ylim((0, 1)) ax_strength.set_ylabel("Texture strength (M-index)") - ax_strength.set_xlabel(r"Strain ($\dot{ε}_0 t$)") + ax_strength.set_xlabel(r"Strain ($\dot{D}_0 t = γ/2$)") - colors = [] - data_Kaminski2001 = _io.read_scsv( - _io.data("thirdparty") / "Kaminski2001_GBMshear.scsv" - ) - lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M0) / 200 * timestop, - data_Kaminski2001.angle_M0, - alpha=0.33, - label=r"$M^{\ast}=0$", - ) - colors.append(lines[0].get_color()) - lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M50) / 200 * timestop, - data_Kaminski2001.angle_M50, - alpha=0.33, - label=r"$M^{\ast}=50$", - ) - colors.append(lines[0].get_color()) - lines = ax_mean.plot( - np.asarray(data_Kaminski2001.equivalent_strain_M200) / 200 * timestop, - data_Kaminski2001.angle_M200, - alpha=0.33, - label=r"$M^{\ast}=200$", - ) - colors.append(lines[0].get_color()) - - for i, (misorient_angles, misorient_indices) in enumerate(zip(angles, indices)): - timestamps = np.linspace(0, timestop, len(misorient_angles)) + for i, (θ_target, θ, m_indices) in enumerate(zip(target_angles, angles, indices)): marker, label = _get_marker_and_label(angles, i, markers, labels) - ax_mean.plot( - timestamps, - misorient_angles, - marker, - markersize=5, - alpha=0.33, - label=label, - color=colors[i], - ) + lines = ax_mean.plot(strains, θ_target, alpha=0.66, label=label) + color = lines[0].get_color() + ax_mean.plot(strains, θ, marker, markersize=5, alpha=0.33, color=color) ax_strength.plot( - timestamps, - misorient_indices, + strains, + m_indices, marker, markersize=5, alpha=0.33, + color=color, label=label, - color=colors[i], ) if θ_fse is not None: - ax_mean.plot(timestamps, θ_fse, "r--", alpha=0.33, label="calc. FSE") ax_mean.plot( - timestamps, - [ - 90 - np.rad2deg( - np.arctan(np.sqrt(5e-6**2 * t**2 + 1) + 5e-6 * t) - ) - for t in timestamps - ], - "k", - label="true FSE", + strains, θ_fse, linestyle=(0, (5, 5)), alpha=0.66, label="FSE" ) - if labels is not None: ax_mean.legend() + ax_strength.legend() fig.savefig(_io.resolve_path(savefile)) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 595bf9db..bd6ea84e 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -2,13 +2,16 @@ import contextlib as cl import numpy as np +from numpy import testing as nt +from scipy.interpolate import PchipInterpolator from pydrex import diagnostics as _diagnostics +from pydrex import io as _io from pydrex import logger as _log from pydrex import minerals as _minerals from pydrex import stats as _stats -from pydrex import visualisation as _vis from pydrex import utils as _utils +from pydrex import visualisation as _vis # Subdirectory of `outdir` used to store outputs from these tests. SUBDIR = "2d_simple_shear" @@ -40,8 +43,10 @@ def test_dvdx_GBM( """ strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. - timestamps = np.linspace(0, 2e5, 501) # Solve until D₀t=1 ('shear strain' γ=2). - n_timesteps = len(timestamps) - 1 + timestamps = np.linspace(0, 2e5, 201) # Solve until D₀t=1 ('shear strain' γ=2). + n_timesteps = len(timestamps) + i_first_cpo = 50 # First index where Bingham averages are sufficiently stable. + i_strain_50p = [0, 50, 100, 150, 200] # Indices for += 50% strain. def get_velocity_gradient(x): # It is independent of time or position in this test. @@ -51,21 +56,15 @@ def get_velocity_gradient(x): shear_direction = [0, 1, 0] # Used to calculate the angular diagnostics. - # Theoretical FSE axis for simple shear. - # We want the angle from the Y axis (shear direction), so subtract from 90. - θ_fse_eq = [ - 90 - _utils.angle_fse_simpleshear(t * strain_rate) for t in timestamps - ] - - # Optional logging and plotting setup. + # Output setup with optional logging and data series labels. + θ_fse = [45] + angles = [] + indices = [] optional_logging = cl.nullcontext() if outdir is not None: out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_GBM" optional_logging = _log.logfile_enable(f"{out_basepath}.log") - θ_fse = [45] labels = [] - angles = [] - indices = [] with optional_logging: for p, params in enumerate( @@ -83,7 +82,7 @@ def get_velocity_gradient(x): _log.info( "step %s/%s (t=%s) with velocity gradient: %s", t, - n_timesteps, + n_timesteps - 1, time, get_velocity_gradient(None).flatten(), ) @@ -96,7 +95,6 @@ def get_velocity_gradient(x): fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) _log.info("› strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) if p == 0 and outdir is not None: - # θ_fse.append(90 - np.rad2deg(np.arctan2(fse_v[1], fse_v[0]))) θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) @@ -129,21 +127,78 @@ def get_velocity_gradient(x): postfix=f"M{params['gbm_mobility']}", ) - # Check that FSE is correct. - np.testing.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-6, atol=0) + # Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`. + _log.info("interpolating target misorientation angles...") + strains = timestamps * strain_rate + data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2001_GBMshear.scsv") + cs_M0 = PchipInterpolator( + _utils.skip_nans(data.equivalent_strain_M0) / 200, + _utils.skip_nans(data.angle_M0), + ) + cs_M50 = PchipInterpolator( + _utils.skip_nans(data.equivalent_strain_M50) / 200, + _utils.skip_nans(data.angle_M50), + ) + cs_M200 = PchipInterpolator( + _utils.skip_nans(data.equivalent_strain_M200) / 200, + _utils.skip_nans(data.angle_M200), + ) + target_angles = [cs_M0(strains), cs_M50(strains), cs_M200(strains)] # Optionally plot figure. if outdir is not None: _vis.simple_shear_stationary_2d( + strains, + target_angles, angles, indices, - timestop=timestamps[-1], - savefile=f"{out_basepath}.png", + savefile=f"{out_basepath}.pdf", markers=("o", "v", "s"), - labels=labels, θ_fse=θ_fse, + labels=labels, ) + # Check that FSE is correct. + # First, get theoretical FSE axis for simple shear. + # We want the angle from the Y axis (shear direction), so subtract from 90. + θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] + nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) + + # Check Bingham angles, ignoring the first portion. + # Average orientations of near-isotropic distributions are unstable. + nt.assert_allclose( + θ_fse[i_first_cpo:], angles[0][i_first_cpo:], rtol=0.1, atol=0 + ) + nt.assert_allclose( + target_angles[0][i_first_cpo:], angles[0][i_first_cpo:], rtol=0.1, atol=0 + ) + nt.assert_allclose( + target_angles[1][i_first_cpo:], angles[1][i_first_cpo:], rtol=0, atol=5.7 + ) + nt.assert_allclose( + target_angles[2][i_first_cpo:], angles[2][i_first_cpo:], rtol=0, atol=5.5 + ) + + # Check texture strength (M-index) at strains of 0%, 50%, 100%, 150% & 200%. + nt.assert_allclose( + [0.29, 0.28, 0.26, 0.24, 0.24], + indices[1].take(i_strain_50p), + rtol=0, + atol=0.02, + ) + nt.assert_allclose( + [0.29, 0.27, 0.26, 0.42, 0.6], + indices[1].take(i_strain_50p), + rtol=0, + atol=0.02, + ) + nt.assert_allclose( + [0.29, 0.3, 0.52, 0.83, 0.89], + indices[2].take(i_strain_50p), + rtol=0, + atol=0.02, + ) + def test_dudz_GBS( self, params_Kaminski2004_fig4_triangles, # GBS = 0.4 @@ -159,8 +214,10 @@ def test_dudz_GBS( """ strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. - timestamps = np.linspace(0, 2e5, 501) # Solve until D₀t=1 ('shear strain' γ=2). - n_timesteps = len(timestamps) - 1 + timestamps = np.linspace(0, 2e5, 201) # Solve until D₀t=1 ('shear strain' γ=2). + n_timesteps = len(timestamps) + i_first_cpo = 50 # First index where Bingham averages are sufficiently stable. + i_strain_50p = [0, 50, 100, 150, 200] # Indices for += 50% strain. def get_velocity_gradient(x): # It is independent of time or position in this test. @@ -170,21 +227,15 @@ def get_velocity_gradient(x): shear_direction = [1, 0, 0] # Used to calculate the angular diagnostics. - # Theoretical FSE axis for simple shear. - # We want the angle from the Y axis (shear direction), so subtract from 90. - θ_fse_eq = [ - 90 - _utils.angle_fse_simpleshear(t * strain_rate) for t in timestamps - ] - - # Optional plotting and logging setup. + # Output setup with optional logging and data series labels. + θ_fse = [45] + angles = [] + indices = [] optional_logging = cl.nullcontext() if outdir is not None: out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dudz_GBS" optional_logging = _log.logfile_enable(f"{out_basepath}.log") - θ_fse = [45] labels = [] - angles = [] - indices = [] with optional_logging: for p, params in enumerate( @@ -202,7 +253,7 @@ def get_velocity_gradient(x): _log.info( "step %s/%s (t=%s) with velocity gradient: %s", t, - n_timesteps, + n_timesteps - 1, time, get_velocity_gradient(None).flatten(), ) @@ -248,14 +299,19 @@ def get_velocity_gradient(x): ) # Check that FSE is correct. - np.testing.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-6, atol=0) + # First get theoretical FSE axis for simple shear. + # We want the angle from the Y axis (shear direction), so subtract from 90. + θ_fse_eq = [ + 90 - _utils.angle_fse_simpleshear(t * strain_rate) for t in timestamps + ] + nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) # Optionally plot figure. if outdir is not None: _vis.simple_shear_stationary_2d( + timestamps[-1] * strain_rate, angles, indices, - timestop=timestamps[-1], savefile=f"{out_basepath}.png", markers=("o", "v", "s"), labels=labels, From c3f9899173794789e999164c61b0fa6b6407762a Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Thu, 20 Jul 2023 23:58:17 +1000 Subject: [PATCH 12/30] feat: Allow using strain data to set polefigure suptitles --- src/pydrex/cli.py | 15 ++++++++++++++- src/pydrex/visualisation.py | 37 +++++++++++++++++++++++++------------ 2 files changed, 39 insertions(+), 13 deletions(-) diff --git a/src/pydrex/cli.py b/src/pydrex/cli.py index 8b8d3055..8f39612f 100644 --- a/src/pydrex/cli.py +++ b/src/pydrex/cli.py @@ -6,6 +6,7 @@ import numpy as np +from pydrex import io as _io from pydrex import exceptions as _err from pydrex import logger as _log from pydrex import minerals as _minerals @@ -31,7 +32,9 @@ def __call__(self): if args.range is None: i_range = None else: - i_range = range(*(int(s) for s in args.range.split(":"))) + start, stop_ex, step = (int(s) for s in args.range.split(":")) + # Make command line start:stop:step stop-inclusive, it's more intuitive. + i_range = range(start, stop_ex + step, step) density_kwargs = {"kernel": args.kernel} if args.smoothing is not None: @@ -58,6 +61,7 @@ def __call__(self): i_range=i_range, density=args.density, savefile=args.out, + strains=_io.read_scsv(args.scsv).strain, **density_kwargs, ) except (argparse.ArgumentError, ValueError, _err.Error) as e: @@ -73,6 +77,15 @@ def _get_args(self) -> argparse.Namespace: help="range of strain indices to be plotted, in the format start:stop:step", default=None, ) + parser.add_argument( + "-f", + "--scsv", + help=( + "path to SCSV file with a column named 'strain'" + + " that lists shear strain percentages for each strain index" + ), + default=None, + ) parser.add_argument( "-p", "--postfix", diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index 9538c8cc..31f7c31a 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -23,7 +23,13 @@ def polefigures( - orientations, ref_axes, i_range, density=False, savefile="polefigures.png", **kwargs + orientations, + ref_axes, + i_range, + density=False, + savefile="polefigures.png", + strains=None, + **kwargs, ): """Plot pole figures of a series of (Nx3x3) orientation matrix stacks. @@ -42,18 +48,25 @@ def polefigures( grid = fig.add_gridspec( 4, n_orientations, height_ratios=((1, 3, 3, 3)), hspace=0, wspace=0.2 ) - fig_time = fig.add_subfigure(grid[0, :]) + fig_strain = fig.add_subfigure(grid[0, :]) first_row = 1 - fig_time.suptitle( - f"N ⋅ (max strain) / {i_range.stop}", x=0.5, y=0.85, fontsize="small" - ) - ax_time = fig_time.add_subplot(111) - ax_time.set_frame_on(False) - ax_time.grid(False) - ax_time.yaxis.set_visible(False) - ax_time.xaxis.set_tick_params(labelsize="x-small", length=0) - ax_time.set_xticks(list(i_range)) - ax_time.set_xlim((-i_range.step / 2, i_range.stop - i_range.step / 2)) + ax_strain = fig_strain.add_subplot(111) + ax_strain.set_xlim((i_range.start, i_range.stop)) + + if strains is None: + fig_strain.suptitle( + f"N ⋅ (max strain) / {i_range.stop}", x=0.5, y=0.85, fontsize="small" + ) + ax_strain.set_xticks(list(i_range)) + else: + fig_strain.suptitle("strain (%)", x=0.5, y=0.85, fontsize="small") + ax_strain.set_xticks(strains[i_range.start : i_range.stop : i_range.step]) + + ax_strain.set_frame_on(False) + ax_strain.grid(False) + ax_strain.yaxis.set_visible(False) + ax_strain.xaxis.set_tick_params(labelsize="x-small", length=0) + ax_strain.set_xlim((-i_range.step / 2, i_range.stop - i_range.step / 2)) fig100 = fig.add_subfigure( grid[first_row, :], edgecolor=plt.rcParams["grid.color"], linewidth=1 From 5e5c6e739ac1e8509ed5f3914ade267da5c71a86 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 00:00:32 +1000 Subject: [PATCH 13/30] test: Plot symmetry diagnostic instead of strength and tweak asserts --- src/pydrex/diagnostics.py | 2 +- src/pydrex/visualisation.py | 28 ++++++++-------- tests/test_simple_shear_2d.py | 63 +++++++++++++++++++++++------------ 3 files changed, 57 insertions(+), 36 deletions(-) diff --git a/src/pydrex/diagnostics.py b/src/pydrex/diagnostics.py index 7331010a..458dcad5 100644 --- a/src/pydrex/diagnostics.py +++ b/src/pydrex/diagnostics.py @@ -240,7 +240,7 @@ def coaxial_index(orientations, axis1="b", axis2="a"): The BA index of [Mainprice et al. 2015](https://doi.org/10.1144/SP409.8) is derived from the eigenvalue `symmetry` diagnostics and measures point vs girdle - symmetry in the aggregate. $BA \in [0, 1]$ where $BA = 0$ corresponds to a perfect + symmetry in the aggregate. $BA ∈ [0, 1]$ where $BA = 0$ corresponds to a perfect axial girdle texture and $BA = 1$ represents a point symmetry texture assuming that the random component $R$ is negligible. May be less susceptible to random fluctuations compared to the raw eigenvalue diagnostics. diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index 31f7c31a..5173660e 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -159,7 +159,7 @@ def simple_shear_stationary_2d( strains, target_angles, angles, - indices, + point100_symmetry, savefile="pydrex_simple_shear_stationary_2d.png", markers=("."), θ_fse=None, @@ -173,21 +173,23 @@ def simple_shear_stationary_2d( ax_mean.tick_params(labelbottom=False) ax_mean.set_xlim((strains[0], strains[-1])) ax_mean.set_ylim((0, 60)) - ax_strength = fig.add_subplot(grid[1], sharex=ax_mean) - ax_strength.set_xlim((strains[0], strains[-1])) - ax_strength.set_ylim((0, 1)) - ax_strength.set_ylabel("Texture strength (M-index)") - ax_strength.set_xlabel(r"Strain ($\dot{D}_0 t = γ/2$)") - - for i, (θ_target, θ, m_indices) in enumerate(zip(target_angles, angles, indices)): + ax_symmetry = fig.add_subplot(grid[1], sharex=ax_mean) + ax_symmetry.set_xlim((strains[0], strains[-1])) + ax_symmetry.set_ylim((0, 1)) + ax_symmetry.set_ylabel(r"Texture symmetry ($P_{[100]}$)") + ax_symmetry.set_xlabel(r"Strain ($\dot{D}_0 t = γ/2$)") + + for i, (θ_target, θ, point100) in enumerate( + zip(target_angles, angles, point100_symmetry) + ): marker, label = _get_marker_and_label(angles, i, markers, labels) lines = ax_mean.plot(strains, θ_target, alpha=0.66, label=label) color = lines[0].get_color() ax_mean.plot(strains, θ, marker, markersize=5, alpha=0.33, color=color) - ax_strength.plot( + ax_symmetry.plot( strains, - m_indices, + point100, marker, markersize=5, alpha=0.33, @@ -196,12 +198,10 @@ def simple_shear_stationary_2d( ) if θ_fse is not None: - ax_mean.plot( - strains, θ_fse, linestyle=(0, (5, 5)), alpha=0.66, label="FSE" - ) + ax_mean.plot(strains, θ_fse, linestyle=(0, (5, 5)), alpha=0.66, label="FSE") if labels is not None: ax_mean.legend() - ax_strength.legend() + ax_symmetry.legend() fig.savefig(_io.resolve_path(savefile)) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index bd6ea84e..7c9e9e04 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -59,7 +59,7 @@ def get_velocity_gradient(x): # Output setup with optional logging and data series labels. θ_fse = [45] angles = [] - indices = [] + point100_symmetry = [] optional_logging = cl.nullcontext() if outdir is not None: out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_GBM" @@ -99,8 +99,11 @@ def get_velocity_gradient(x): _diagnostics.smallest_angle(fse_v, shear_direction) ) - misorient_indices = np.zeros(n_timesteps) - misorient_angles = np.zeros(n_timesteps) + assert len(mineral.orientations) == n_timesteps + assert len(mineral.fractions) == n_timesteps + assert len(θ_fse) == n_timesteps + texture_symmetry = np.zeros(n_timesteps) + mean_angles = np.zeros(n_timesteps) # Loop over first dimension (time steps) of orientations. for idx, matrices in enumerate(mineral.orientations): orientations_resampled, _ = _stats.resample_orientations( @@ -110,18 +113,19 @@ def get_velocity_gradient(x): orientations_resampled, axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], ) - misorient_angles[idx] = _diagnostics.smallest_angle( + mean_angles[idx] = _diagnostics.smallest_angle( direction_mean, shear_direction ) - misorient_indices[idx] = _diagnostics.misorientation_index( - orientations_resampled - ) + texture_symmetry[idx] = _diagnostics.symmetry( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + )[0] # Optionally store plotting metadata. if outdir is not None: labels.append(f"$M^∗$ = {params['gbm_mobility']}") - angles.append(misorient_angles) - indices.append(misorient_indices) + angles.append(mean_angles) + point100_symmetry.append(texture_symmetry) mineral.save( f"{out_basepath}.npz", postfix=f"M{params['gbm_mobility']}", @@ -147,12 +151,29 @@ def get_velocity_gradient(x): # Optionally plot figure. if outdir is not None: + schema = { + "delimiter": ",", + "missing": "-", + "fields": [ + { + "name": "strain", + "type": "integer", + "unit": "percent", + "fill": 999999, + } + ], + } + _io.save_scsv( + f"{out_basepath}_strains.scsv", + schema, + [[int(γ * 200) for γ in strains]], + ) _vis.simple_shear_stationary_2d( strains, target_angles, angles, - indices, - savefile=f"{out_basepath}.pdf", + point100_symmetry, + savefile=f"{out_basepath}.png", markers=("o", "v", "s"), θ_fse=θ_fse, labels=labels, @@ -179,24 +200,24 @@ def get_velocity_gradient(x): target_angles[2][i_first_cpo:], angles[2][i_first_cpo:], rtol=0, atol=5.5 ) - # Check texture strength (M-index) at strains of 0%, 50%, 100%, 150% & 200%. + # Check point symmetry of [100] at strains of 0%, 50%, 100%, 150% & 200%. nt.assert_allclose( - [0.29, 0.28, 0.26, 0.24, 0.24], - indices[1].take(i_strain_50p), + [0.015, 0.11, 0.19, 0.27, 0.34], + point100_symmetry[0].take(i_strain_50p), rtol=0, - atol=0.02, + atol=0.015, ) nt.assert_allclose( - [0.29, 0.27, 0.26, 0.42, 0.6], - indices[1].take(i_strain_50p), + [0.015, 0.15, 0.33, 0.57, 0.72], + point100_symmetry[1].take(i_strain_50p), rtol=0, - atol=0.02, + atol=0.015, ) nt.assert_allclose( - [0.29, 0.3, 0.52, 0.83, 0.89], - indices[2].take(i_strain_50p), + [0.015, 0.22, 0.64, 0.86, 0.91], + point100_symmetry[2].take(i_strain_50p), rtol=0, - atol=0.02, + atol=0.015, ) def test_dudz_GBS( From dc3e5831eabe93f3161e2863ac68842514aa5555 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 02:38:52 +1000 Subject: [PATCH 14/30] test: Fix GBS test and header in GBM datafile --- .../thirdparty/Kaminski2001_GBMshear.scsv | 3 +- tests/test_simple_shear_2d.py | 85 +++++++++++++------ 2 files changed, 59 insertions(+), 29 deletions(-) diff --git a/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv b/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv index 769d4df9..ff5180d2 100644 --- a/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv +++ b/src/pydrex/data/thirdparty/Kaminski2001_GBMshear.scsv @@ -1,8 +1,7 @@ --- # Digitized data from Figure 5 in Kaminski & Ribe 2001, https://doi.org/10.1016%2Fs0012-821x%2801%2900356-9 # This data records the mean angle of the a-axes of an olivine aggregate from the horizontal shear direction. -# The flow field is a static simple shear and the aggregate starts with all grain orientations restricted to the first quadrant -# in the 2D deformation plane (so the mean angle always starts at 45°). +# The flow field is a static simple shear and the aggregate starts with random (i.e. isotropically distributed) grain orientations. schema: delimiter: ',' diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 7c9e9e04..17efa6a8 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -94,14 +94,11 @@ def get_velocity_gradient(x): ) fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) _log.info("› strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) - if p == 0 and outdir is not None: + if p == 0: θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) - assert len(mineral.orientations) == n_timesteps - assert len(mineral.fractions) == n_timesteps - assert len(θ_fse) == n_timesteps texture_symmetry = np.zeros(n_timesteps) mean_angles = np.zeros(n_timesteps) # Loop over first dimension (time steps) of orientations. @@ -132,7 +129,7 @@ def get_velocity_gradient(x): ) # Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`. - _log.info("interpolating target misorientation angles...") + _log.info("interpolating target CPO angles...") strains = timestamps * strain_rate data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2001_GBMshear.scsv") cs_M0 = PchipInterpolator( @@ -235,7 +232,7 @@ def test_dudz_GBS( """ strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. - timestamps = np.linspace(0, 2e5, 201) # Solve until D₀t=1 ('shear strain' γ=2). + timestamps = np.linspace(0, 5e5, 201) # Solve until D₀t=2.5 ('shear' γ=5). n_timesteps = len(timestamps) i_first_cpo = 50 # First index where Bingham averages are sufficiently stable. i_strain_50p = [0, 50, 100, 150, 200] # Indices for += 50% strain. @@ -251,7 +248,7 @@ def get_velocity_gradient(x): # Output setup with optional logging and data series labels. θ_fse = [45] angles = [] - indices = [] + point100_symmetry = [] optional_logging = cl.nullcontext() if outdir is not None: out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dudz_GBS" @@ -286,13 +283,13 @@ def get_velocity_gradient(x): ) fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) _log.info("› strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) - if p == 0 and outdir is not None: + if p == 0: θ_fse.append( _diagnostics.smallest_angle(fse_v, shear_direction) ) - misorient_indices = np.zeros(n_timesteps) - misorient_angles = np.zeros(n_timesteps) + texture_symmetry = np.zeros(n_timesteps) + mean_angles = np.zeros(n_timesteps) # Loop over first dimension (time steps) of orientations. for idx, matrices in enumerate(mineral.orientations): orientations_resampled, _ = _stats.resample_orientations( @@ -302,39 +299,73 @@ def get_velocity_gradient(x): orientations_resampled, axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], ) - misorient_angles[idx] = _diagnostics.smallest_angle( - direction_mean, [1, 0, 0] - ) - misorient_indices[idx] = _diagnostics.misorientation_index( - orientations_resampled + mean_angles[idx] = _diagnostics.smallest_angle( + direction_mean, shear_direction ) + texture_symmetry[idx] = _diagnostics.symmetry( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + )[0] # Optionally store plotting metadata. if outdir is not None: labels.append(f"$f_{{gbs}}$ = {params['gbs_threshold']}") - angles.append(misorient_angles) - indices.append(misorient_indices) + angles.append(mean_angles) + point100_symmetry.append(texture_symmetry) mineral.save( f"{out_basepath}.npz", postfix=f"X{params['gbs_threshold']}", ) - # Check that FSE is correct. - # First get theoretical FSE axis for simple shear. - # We want the angle from the Y axis (shear direction), so subtract from 90. - θ_fse_eq = [ - 90 - _utils.angle_fse_simpleshear(t * strain_rate) for t in timestamps - ] - nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) + # Interpolate Kaminski & Ribe, 2001 data to get target angles at `strains`. + _log.info("interpolating target CPO angles...") + strains = timestamps * strain_rate + data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2004_GBSshear.scsv") + cs_X0 = PchipInterpolator( + _utils.skip_nans(data.dimensionless_time_X0), + 45 + _utils.skip_nans(data.angle_X0), + ) + cs_X0d2 = PchipInterpolator( + _utils.skip_nans(data.dimensionless_time_X0d2), + 45 + _utils.skip_nans(data.angle_X0d2), + ) + cs_X0d4 = PchipInterpolator( + _utils.skip_nans(data.dimensionless_time_X0d4), + 45 + _utils.skip_nans(data.angle_X0d4), + ) + target_angles = [cs_X0(strains), cs_X0d2(strains), cs_X0d4(strains)] # Optionally plot figure. if outdir is not None: + schema = { + "delimiter": ",", + "missing": "-", + "fields": [ + { + "name": "strain", + "type": "integer", + "unit": "percent", + "fill": 999999, + } + ], + } + _io.save_scsv( + f"{out_basepath}_strains.scsv", + schema, + [[int(γ * 200) for γ in strains]], + ) _vis.simple_shear_stationary_2d( - timestamps[-1] * strain_rate, + strains, + target_angles, angles, - indices, + point100_symmetry, savefile=f"{out_basepath}.png", markers=("o", "v", "s"), labels=labels, - θ_fse=θ_fse, ) + + # Check that FSE is correct. + # First, get theoretical FSE axis for simple shear. + # We want the angle from the Y axis (shear direction), so subtract from 90. + θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] + nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) From 31669abaf73544ce02e39214f8feb3f09ba3f434 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 03:45:13 +1000 Subject: [PATCH 15/30] fix: GBS now carries over previous orientation Fixes the grain boundary sliding routine to match what is defined in the original fortran, where dotacs and rot are set to zero for grains that are sliding, i.e. there is zero rotation _rate_, but not resetting the whole orientation matrix to the initial one. Now simply carries over the previous orientation. Also adds symmetry asserts to the GBS test. --- src/pydrex/minerals.py | 3 +-- src/pydrex/visualisation.py | 23 ----------------------- tests/test_simple_shear_2d.py | 32 ++++++++++++++++++++++++++------ 3 files changed, 27 insertions(+), 31 deletions(-) diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index 636aff28..61321252 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -362,8 +362,7 @@ def apply_gbs(orientations, fractions, config): len(np.nonzero(mask)) / len(fractions), ) # No rotation: carry over previous orientations. - # TODO: Should we really be resetting to initial orientations here? - orientations[mask, :, :] = self.orientations[0][mask, :, :] + orientations[mask, :, :] = self.orientations[-1][mask, :, :] fractions[mask] = config["gbs_threshold"] / self.n_grains fractions /= fractions.sum() _log.debug( diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index 5173660e..edee6db8 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -125,27 +125,6 @@ def polefigures( fig.savefig(_io.resolve_path(savefile)) -def check_marker_seq(func): - """Raises a `ValueError` if number of markers and data series don't match. - - The decorated function is expected to take the data as the first positional - argument, and a keyword argument called `markers`. - - """ - - @ft.wraps(func) - def wrapper(data, *args, **kwargs): - markers = kwargs["markers"] - if len(data) % len(markers) != 0: - raise ValueError( - "Number of data series must be divisible by number of markers." - + f" You've supplied {len(data)} data series and {len(markers)} markers." - ) - func(data, *args, **kwargs) - - return wrapper - - def _get_marker_and_label(data, seq_index, markers, labels=None): marker = markers[int(seq_index / (len(data) / len(markers)))] label = None @@ -154,7 +133,6 @@ def _get_marker_and_label(data, seq_index, markers, labels=None): return marker, label -@check_marker_seq def simple_shear_stationary_2d( strains, target_angles, @@ -214,7 +192,6 @@ def _lag_2d_corner_flow(θ): ) -@check_marker_seq def corner_flow_2d( x_paths, z_paths, diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 17efa6a8..0aab6065 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -219,9 +219,9 @@ def get_velocity_gradient(x): def test_dudz_GBS( self, - params_Kaminski2004_fig4_triangles, # GBS = 0.4 - params_Kaminski2004_fig4_squares, # GBS = 0.2 params_Kaminski2004_fig4_circles, # GBS = 0 + params_Kaminski2004_fig4_squares, # GBS = 0.2 + params_Kaminski2004_fig4_triangles, # GBS = 0.4 rng, outdir, ): @@ -232,10 +232,10 @@ def test_dudz_GBS( """ strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. - timestamps = np.linspace(0, 5e5, 201) # Solve until D₀t=2.5 ('shear' γ=5). + timestamps = np.linspace(0, 5e5, 251) # Solve until D₀t=2.5 ('shear' γ=5). n_timesteps = len(timestamps) i_first_cpo = 50 # First index where Bingham averages are sufficiently stable. - i_strain_50p = [0, 50, 100, 150, 200] # Indices for += 50% strain. + i_strain_100p = [0, 50, 100, 150, 200] # Indices for += 100% strain. def get_velocity_gradient(x): # It is independent of time or position in this test. @@ -258,9 +258,9 @@ def get_velocity_gradient(x): with optional_logging: for p, params in enumerate( ( - params_Kaminski2004_fig4_triangles, # GBS = 0.4 - params_Kaminski2004_fig4_squares, # GBS = 0.2 params_Kaminski2004_fig4_circles, # GBS = 0 + params_Kaminski2004_fig4_squares, # GBS = 0.2 + params_Kaminski2004_fig4_triangles, # GBS = 0.4 ), ): mineral = _minerals.Mineral( @@ -369,3 +369,23 @@ def get_velocity_gradient(x): # We want the angle from the Y axis (shear direction), so subtract from 90. θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) + + # Check point symmetry of [100] at strains of 0%, 100%, 200%, 300% & 400%. + nt.assert_allclose( + [0.015, 0.52, 0.86, 0.93, 0.94], + point100_symmetry[0].take(i_strain_100p), + rtol=0, + atol=0.015, + ) + nt.assert_allclose( + [0.015, 0.4, 0.72, 0.77, 0.79], + point100_symmetry[1].take(i_strain_100p), + rtol=0, + atol=0.015, + ) + nt.assert_allclose( + [0.015, 0.39, 0.58, 0.61, 0.62], + point100_symmetry[2].take(i_strain_100p), + rtol=0, + atol=0.015, + ) From 89b53905d93bc2d227b40593cc23310a572e1818 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 05:47:48 +1000 Subject: [PATCH 16/30] test: Fix GBS test asserts and add commented versions for SCCS --- src/pydrex/visualisation.py | 2 -- tests/test_simple_shear_2d.py | 25 ++++++++++++++++++++++--- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index edee6db8..8867e116 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -1,6 +1,4 @@ """> PyDRex: Visualisation functions for test outputs and examples.""" -import functools as ft - import numpy as np from matplotlib import projections as mproj from matplotlib import pyplot as plt diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 0aab6065..89bdf1d1 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -118,6 +118,16 @@ def get_velocity_gradient(x): axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], )[0] + # Uncomment to use SCCS axis (hexagonal symmetry) for the angle instead. + # mean_angles = np.array( + # [ + # _diagnostics.smallest_angle( + # _diagnostics.anisotropy(v)[1][2, :], shear_direction + # ) + # for v in _minerals.voigt_averages([mineral], params) + # ] + # ) + # Optionally store plotting metadata. if outdir is not None: labels.append(f"$M^∗$ = {params['gbm_mobility']}") @@ -234,7 +244,6 @@ def test_dudz_GBS( strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. timestamps = np.linspace(0, 5e5, 251) # Solve until D₀t=2.5 ('shear' γ=5). n_timesteps = len(timestamps) - i_first_cpo = 50 # First index where Bingham averages are sufficiently stable. i_strain_100p = [0, 50, 100, 150, 200] # Indices for += 100% strain. def get_velocity_gradient(x): @@ -307,6 +316,16 @@ def get_velocity_gradient(x): axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], )[0] + # Uncomment to use SCCS axis (hexagonal symmetry) for the angle instead. + # mean_angles = np.array( + # [ + # _diagnostics.smallest_angle( + # _diagnostics.anisotropy(v)[1][2, :], shear_direction + # ) + # for v in _minerals.voigt_averages([mineral], params) + # ] + # ) + # Optionally store plotting metadata. if outdir is not None: labels.append(f"$f_{{gbs}}$ = {params['gbs_threshold']}") @@ -378,13 +397,13 @@ def get_velocity_gradient(x): atol=0.015, ) nt.assert_allclose( - [0.015, 0.4, 0.72, 0.77, 0.79], + [0.015, 0.42, 0.71, 0.77, 0.79], point100_symmetry[1].take(i_strain_100p), rtol=0, atol=0.015, ) nt.assert_allclose( - [0.015, 0.39, 0.58, 0.61, 0.62], + [0.015, 0.36, 0.57, 0.6, 0.62], point100_symmetry[2].take(i_strain_100p), rtol=0, atol=0.015, From a6f9afe70c323b888d26a97c2ff62dc65f9aa376 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 20 Jul 2023 20:00:56 +0000 Subject: [PATCH 17/30] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/pydrex/cli.py | 2 +- src/pydrex/minerals.py | 2 +- tools/pfplot.py | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/pydrex/cli.py b/src/pydrex/cli.py index 8f39612f..911d709c 100644 --- a/src/pydrex/cli.py +++ b/src/pydrex/cli.py @@ -6,8 +6,8 @@ import numpy as np -from pydrex import io as _io from pydrex import exceptions as _err +from pydrex import io as _io from pydrex import logger as _log from pydrex import minerals as _minerals from pydrex import stats as _stats diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index 61321252..1e1627f2 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -12,8 +12,8 @@ from pydrex import exceptions as _err from pydrex import io as _io from pydrex import logger as _log -from pydrex import tensors as _tensors from pydrex import stats as _stats +from pydrex import tensors as _tensors OLIVINE_STIFFNESS = np.array( [ diff --git a/tools/pfplot.py b/tools/pfplot.py index 90e4edd3..06563aa0 100644 --- a/tools/pfplot.py +++ b/tools/pfplot.py @@ -1,4 +1,5 @@ import sys + from pydrex import visualisation as _vis _vis.polefigures( From 71a98b6b874429846467009136e7e24fde7a8490 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 06:57:55 +1000 Subject: [PATCH 18/30] fix: Add missing data file for GBS --- .../thirdparty/Kaminski2004_GBSshear.scsv | 61 +++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 src/pydrex/data/thirdparty/Kaminski2004_GBSshear.scsv diff --git a/src/pydrex/data/thirdparty/Kaminski2004_GBSshear.scsv b/src/pydrex/data/thirdparty/Kaminski2004_GBSshear.scsv new file mode 100644 index 00000000..57721214 --- /dev/null +++ b/src/pydrex/data/thirdparty/Kaminski2004_GBSshear.scsv @@ -0,0 +1,61 @@ +--- +# Digitized data from Figure 4 in Kaminski & Ribe 2004, https://doi.org/10.1111%2Fj.1365-246x.2004.02308.x +# This data records the mean angle of the a-axes of an olivine aggregate in the shear plane. +# The flow field is a static simple shear and the aggregate starts with random (i.e. isotropically distributed) grain orientations. + +schema: + delimiter: ',' + missing: '-' + fields: + - name: dimensionless_time_X0 + type: float + fill: NaN + unit: "max_strain x time" + - name: angle_X0 + type: float + fill: NaN + unit: degree + - name: dimensionless_time_X0d2 + type: float + fill: NaN + unit: "max_strain x time" + - name: angle_X0d2 + type: float + fill: NaN + unit: degree + - name: dimensionless_time_X0d4 + type: float + fill: NaN + unit: "max_strain x time" + - name: angle_X0d4 + type: float + fill: NaN + unit: degree +--- +dimensionless_time_X0, angle_X0,dimensionless_time_X0d2, angle_X0d2,dimensionless_time_X0d4, angle_X0d4 + 0,0.13513513513513775, 0,0.06756756756757198, -0.004098360655737765,3.552713678800501e-15 + 0.09836065573770492, 0.2702702702702737, 0.09836065573770492, 0.3378378378378404, 0.09631147540983603, 0.27027027027027195 + 0.19672131147540978,-0.2702702702702666, 0.19672131147540978,-0.2702702702702684, 0.194672131147541, -0.4054054054054035 + 0.3032786885245901, -2.297297297297293, 0.3012295081967214,-2.2972972972972965, 0.3012295081967213, -2.3648648648648614 + 0.40163934426229503, -9.189189189189186, 0.3995901639344262, -9.054054054054053, 0.39754098360655715, -8.108108108108098 + 0.49999999999999994,-20.945945945945944, 0.49999999999999994,-19.256756756756758, 0.4959016393442623, -16.486486486486484 + 0.598360655737705, -35, 0.598360655737705,-31.486486486486484, 0.5963114754098366, -26.081081081081045 + 0.6967213114754098, -42.70270270270271, 0.6987704918032793, -38.04054054054055, 0.6946721311475408, -31.351351351351347 + 0.7950819672131146, -45.00000000000001, 0.7971311475409832, -39.93243243243243, 0.7930327868852456, -33.04054054054057 + 0.8975409836065573, -45.4054054054054, 0.8954918032786883, -40.74324324324325, 0.8934426229508197, -33.44594594594596 + 1, -45.4054054054054, 1.0020491803278693, -41.14864864864866, 0.9979508196721318, -34.05405405405407 + 1.0983606557377048, -45.4054054054054, 1.100409836065574, -41.2837837837838, 1.0963114754098358, -34.662162162162176 + 1.1987704918032789, -45.20270270270271, 1.1987704918032789, -41.35135135135137, 1.196721311475411, -35.40540540540541 + 1.2991803278688523, -45.13513513513514, 1.2991803278688525, -41.35135135135137, 1.2950819672131149, -35.810810810810814 + 1.3975409836065575, -45.00000000000001, 1.3975409836065575, -41.55405405405406, 1.395491803278689, -36.418918918918926 + 1.4959016393442621, -45.00000000000001, 1.4959016393442617, -41.75675675675676, 1.4938524590163944, -36.6891891891892 + 1.6024590163934427, -45.00000000000001, 1.6024590163934427, -42.1621621621622, 1.5983606557377048, -36.89189189189188 + 1.7008196721311473, -45.00000000000001, 1.7008196721311477,-42.162162162162126, 1.6967213114754094, -36.891891891891895 + 1.7991803278688525, -45.00000000000001, 1.7991803278688525,-42.567567567567586, 1.797131147540984, -37.2972972972973 + 1.8975409836065575, -45.00000000000001, 1.897540983606557, -42.77027027027028, 1.8954918032786876, -37.50000000000001 + 1.9979508196721312, -45.2027027027027, 1.9979508196721303, -42.97297297297298, 1.9959016393442626, -37.905405405405396 + 2.102459016393443, -45.13513513513514, 2.096311475409836, -43.37837837837841, 2.0942622950819674, -38.31081081081081 + 2.1967213114754096, -45.20270270270271, 2.1946721311475406,-43.378378378378365, 2.192622950819672, -38.51351351351353 + 2.3012295081967222, -45.2027027027027, 2.3012295081967222, -43.37837837837841, 2.2971311475409837, -38.918918918918926 + 2.3995901639344264, -45.2027027027027, 2.3995901639344273, -43.58108108108105, 2.397540983606557, -39.12162162162162 + 2.500000000000001, -45.2027027027027, 2.500000000000001,-43.783783783783804, 2.495901639344262, -39.52702702702706 From 062a27aa197143cb2d4cd89170f7bdd0451fcd43 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 08:03:15 +1000 Subject: [PATCH 19/30] fix: Fix incorrect strain ticks in polefigures --- src/pydrex/visualisation.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index 8867e116..8bafee0b 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -49,22 +49,29 @@ def polefigures( fig_strain = fig.add_subfigure(grid[0, :]) first_row = 1 ax_strain = fig_strain.add_subplot(111) - ax_strain.set_xlim((i_range.start, i_range.stop)) if strains is None: fig_strain.suptitle( f"N ⋅ (max strain) / {i_range.stop}", x=0.5, y=0.85, fontsize="small" ) + ax_strain.set_xlim( + (i_range.start - i_range.step / 2, i_range.stop - i_range.step / 2) + ) ax_strain.set_xticks(list(i_range)) else: fig_strain.suptitle("strain (%)", x=0.5, y=0.85, fontsize="small") ax_strain.set_xticks(strains[i_range.start : i_range.stop : i_range.step]) + ax_strain.set_xlim( + ( + strains[i_range.start] - strains[i_range.step] / 2, + strains[i_range.stop - i_range.step] + strains[i_range.step] / 2, + ) + ) ax_strain.set_frame_on(False) ax_strain.grid(False) ax_strain.yaxis.set_visible(False) ax_strain.xaxis.set_tick_params(labelsize="x-small", length=0) - ax_strain.set_xlim((-i_range.step / 2, i_range.stop - i_range.step / 2)) fig100 = fig.add_subfigure( grid[first_row, :], edgecolor=plt.rcParams["grid.color"], linewidth=1 From 2c08eb1990f189d9a49178a84ce11f1702e781d5 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 08:46:48 +1000 Subject: [PATCH 20/30] test: Fix incorrecto outdir conditional --- tests/test_simple_shear_2d.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 89bdf1d1..dfa91e50 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -128,11 +128,11 @@ def get_velocity_gradient(x): # ] # ) - # Optionally store plotting metadata. + # Store data series and optional plotting metadata. + angles.append(mean_angles) + point100_symmetry.append(texture_symmetry) if outdir is not None: labels.append(f"$M^∗$ = {params['gbm_mobility']}") - angles.append(mean_angles) - point100_symmetry.append(texture_symmetry) mineral.save( f"{out_basepath}.npz", postfix=f"M{params['gbm_mobility']}", @@ -326,11 +326,11 @@ def get_velocity_gradient(x): # ] # ) - # Optionally store plotting metadata. + # Store data series and optional plotting metadata. + angles.append(mean_angles) + point100_symmetry.append(texture_symmetry) if outdir is not None: labels.append(f"$f_{{gbs}}$ = {params['gbs_threshold']}") - angles.append(mean_angles) - point100_symmetry.append(texture_symmetry) mineral.save( f"{out_basepath}.npz", postfix=f"X{params['gbs_threshold']}", From 284a17c0a1780cde64d7c7cc7ca972441946e19d Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 09:13:07 +1000 Subject: [PATCH 21/30] test: Fix tolerance value --- tests/test_simple_shear_2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index dfa91e50..e52e7ce9 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -195,7 +195,7 @@ def get_velocity_gradient(x): # Check Bingham angles, ignoring the first portion. # Average orientations of near-isotropic distributions are unstable. nt.assert_allclose( - θ_fse[i_first_cpo:], angles[0][i_first_cpo:], rtol=0.1, atol=0 + θ_fse[i_first_cpo:], angles[0][i_first_cpo:], rtol=0.11, atol=0 ) nt.assert_allclose( target_angles[0][i_first_cpo:], angles[0][i_first_cpo:], rtol=0.1, atol=0 From 90a845a65794de180e4fd8fdb557b174a3f0b2cc Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 10:38:18 +1000 Subject: [PATCH 22/30] test: Fix minor tolerance issue. --- tests/test_simple_shear_2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index e52e7ce9..fb553349 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -198,7 +198,7 @@ def get_velocity_gradient(x): θ_fse[i_first_cpo:], angles[0][i_first_cpo:], rtol=0.11, atol=0 ) nt.assert_allclose( - target_angles[0][i_first_cpo:], angles[0][i_first_cpo:], rtol=0.1, atol=0 + target_angles[0][i_first_cpo:], angles[0][i_first_cpo:], rtol=0.11, atol=0 ) nt.assert_allclose( target_angles[1][i_first_cpo:], angles[1][i_first_cpo:], rtol=0, atol=5.7 From a3156be9c5311a860bf39008230b0221c5e7b1a4 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 11:44:20 +1000 Subject: [PATCH 23/30] dev: Don't show pytest "short summary" duplicate errors in CI --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8cdd9d41..33a56bd4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,4 +28,4 @@ jobs: sudo apt-get install -y python3-pip python3 -m pip install --upgrade pip python3 -m pip install "$PWD"[test] - python3 -m pytest + python3 -m pytest -rN From 7cc021161ab307d883b1845946d2496a144035dc Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Fri, 21 Jul 2023 11:45:17 +1000 Subject: [PATCH 24/30] refactor: Rename skip_nans -> remove_nans --- src/pydrex/utils.py | 4 ++-- tests/test_simple_shear_2d.py | 24 ++++++++++++------------ 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/pydrex/utils.py b/src/pydrex/utils.py index 043d83af..52b23a6e 100644 --- a/src/pydrex/utils.py +++ b/src/pydrex/utils.py @@ -2,8 +2,8 @@ import numpy as np -def skip_nans(a): - """Skip NaN values in array.""" +def remove_nans(a): + """Remove NaN values from array.""" a = np.asarray(a) return a[~np.isnan(a)] diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index fb553349..e5030371 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -143,16 +143,16 @@ def get_velocity_gradient(x): strains = timestamps * strain_rate data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2001_GBMshear.scsv") cs_M0 = PchipInterpolator( - _utils.skip_nans(data.equivalent_strain_M0) / 200, - _utils.skip_nans(data.angle_M0), + _utils.remove_nans(data.equivalent_strain_M0) / 200, + _utils.remove_nans(data.angle_M0), ) cs_M50 = PchipInterpolator( - _utils.skip_nans(data.equivalent_strain_M50) / 200, - _utils.skip_nans(data.angle_M50), + _utils.remove_nans(data.equivalent_strain_M50) / 200, + _utils.remove_nans(data.angle_M50), ) cs_M200 = PchipInterpolator( - _utils.skip_nans(data.equivalent_strain_M200) / 200, - _utils.skip_nans(data.angle_M200), + _utils.remove_nans(data.equivalent_strain_M200) / 200, + _utils.remove_nans(data.angle_M200), ) target_angles = [cs_M0(strains), cs_M50(strains), cs_M200(strains)] @@ -341,16 +341,16 @@ def get_velocity_gradient(x): strains = timestamps * strain_rate data = _io.read_scsv(_io.data("thirdparty") / "Kaminski2004_GBSshear.scsv") cs_X0 = PchipInterpolator( - _utils.skip_nans(data.dimensionless_time_X0), - 45 + _utils.skip_nans(data.angle_X0), + _utils.remove_nans(data.dimensionless_time_X0), + 45 + _utils.remove_nans(data.angle_X0), ) cs_X0d2 = PchipInterpolator( - _utils.skip_nans(data.dimensionless_time_X0d2), - 45 + _utils.skip_nans(data.angle_X0d2), + _utils.remove_nans(data.dimensionless_time_X0d2), + 45 + _utils.remove_nans(data.angle_X0d2), ) cs_X0d4 = PchipInterpolator( - _utils.skip_nans(data.dimensionless_time_X0d4), - 45 + _utils.skip_nans(data.angle_X0d4), + _utils.remove_nans(data.dimensionless_time_X0d4), + 45 + _utils.remove_nans(data.angle_X0d4), ) target_angles = [cs_X0(strains), cs_X0d2(strains), cs_X0d4(strains)] From a1ffdf638c418ab28b021a56fc320ffa31d1387e Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Sat, 22 Jul 2023 14:35:27 +1000 Subject: [PATCH 25/30] remove: [skip ci] Remove old pfplot script, use pydrex-polefigures --- tools/pfplot.py | 11 ----------- 1 file changed, 11 deletions(-) delete mode 100644 tools/pfplot.py diff --git a/tools/pfplot.py b/tools/pfplot.py deleted file mode 100644 index 06563aa0..00000000 --- a/tools/pfplot.py +++ /dev/null @@ -1,11 +0,0 @@ -import sys - -from pydrex import visualisation as _vis - -_vis.polefigures( - sys.argv[1], - postfix=sys.argv[2], - i_range=range(*[int(s) for s in sys.argv[3].split(":")]), - density=bool(int(sys.argv[4])), - ref_axes=sys.argv[5], -) From 839ba3192d0e1d23545c2b98cae33ef58c855a34 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 24 Jul 2023 11:40:40 +1000 Subject: [PATCH 26/30] fix: Fix incorrect RNG seeding and add --runslow for long tests --- src/pydrex/data/rng/seeds.scsv | 1010 ++++++++++++++++++++++++++++++++ src/pydrex/minerals.py | 6 +- src/pydrex/stats.py | 7 +- src/pydrex/visualisation.py | 5 + tests/conftest.py | 46 +- tests/test_corner_flow_2d.py | 6 +- tests/test_geometry.py | 3 +- tests/test_simple_shear_2d.py | 442 +++++++------- 8 files changed, 1302 insertions(+), 223 deletions(-) create mode 100644 src/pydrex/data/rng/seeds.scsv diff --git a/src/pydrex/data/rng/seeds.scsv b/src/pydrex/data/rng/seeds.scsv new file mode 100644 index 00000000..8c030997 --- /dev/null +++ b/src/pydrex/data/rng/seeds.scsv @@ -0,0 +1,1010 @@ +--- +schema: + delimiter: ',' + missing: '-' + fields: + - name: seeds + type: integer + fill: 0 +--- +seeds +101 +121 +1218 +3052 +4232 +4557 +5675 +8236 +8711 +10113 +10325 +12198 +12467 +12762 +13178 +13515 +14756 +15181 +16535 +16608 +16913 +19182 +20342 +20442 +20519 +20692 +20893 +20909 +21232 +21400 +21590 +24094 +26169 +26974 +28635 +31173 +31719 +33355 +34472 +34851 +36199 +43410 +43879 +44099 +44409 +44886 +45928 +46708 +46908 +47611 +48561 +49754 +50018 +50036 +50296 +50569 +50669 +51406 +51623 +52535 +53478 +54061 +54882 +55276 +56346 +56394 +57961 +59776 +60397 +60411 +60955 +63519 +63566 +64142 +64206 +67640 +70280 +70618 +71666 +71838 +72127 +74830 +75257 +75269 +75509 +76437 +76509 +77002 +77468 +79002 +79073 +80782 +80962 +83466 +84510 +85192 +86036 +86446 +87099 +87126 +90646 +90949 +91002 +91827 +92409 +92420 +94664 +95791 +97951 +98873 +100710 +101899 +104164 +104963 +105813 +106577 +107646 +107784 +108242 +109640 +110784 +111304 +111478 +112182 +113363 +113535 +113663 +114704 +115470 +116042 +116335 +117661 +118109 +120051 +120765 +121729 +124533 +126017 +126657 +129519 +129878 +130272 +130950 +131319 +132147 +132333 +133142 +134203 +134561 +134679 +136585 +136866 +137027 +138114 +138160 +139623 +141670 +142614 +142630 +142677 +143678 +144881 +144972 +145400 +145478 +145637 +146239 +147742 +148686 +149713 +150432 +152078 +152272 +152827 +153859 +154469 +154522 +154661 +157064 +159107 +161425 +161575 +162130 +162194 +163472 +163863 +165201 +165590 +166203 +167019 +168486 +170180 +171159 +172981 +173861 +174256 +175802 +176530 +176625 +176661 +176762 +177884 +177902 +178679 +179030 +179426 +179772 +180748 +182886 +183850 +184169 +185118 +185422 +185920 +186416 +187608 +188228 +191526 +191953 +193226 +194848 +195211 +199081 +199740 +200064 +200969 +201228 +201897 +202641 +203371 +205585 +206823 +207006 +207651 +207694 +209018 +212606 +213511 +213662 +214205 +215936 +216028 +217889 +218486 +218898 +219503 +220767 +221923 +224888 +227310 +227658 +227774 +229028 +229228 +229427 +232214 +233523 +234135 +234590 +235239 +236275 +236990 +237639 +240855 +241246 +241308 +241519 +242014 +244336 +245199 +247561 +247885 +250933 +251730 +254091 +254312 +254421 +254663 +255256 +255659 +255767 +255813 +255850 +257203 +257382 +259837 +260294 +260583 +261521 +261725 +261746 +262517 +262591 +262643 +263132 +263585 +263664 +264112 +266552 +267007 +267096 +267664 +269002 +271054 +271783 +271879 +272072 +272103 +273326 +275228 +275445 +276563 +278368 +278855 +278877 +280582 +280807 +283301 +284463 +284933 +287994 +289416 +291901 +292819 +295332 +295940 +296427 +297562 +298877 +300906 +301643 +302043 +302942 +304463 +305273 +306387 +310707 +311713 +313520 +314883 +315141 +317174 +320123 +320152 +320651 +322318 +323683 +324677 +325364 +326698 +327262 +327333 +327578 +329023 +329785 +332494 +332962 +334043 +335832 +336802 +336919 +338138 +338987 +343049 +345528 +346444 +346585 +347563 +349501 +352775 +355260 +356715 +357839 +358841 +358865 +358948 +359893 +360175 +360396 +361291 +362095 +362290 +363809 +364865 +365604 +370152 +371380 +372149 +373716 +373928 +374838 +374909 +378247 +379464 +379554 +381361 +382202 +382449 +382599 +382866 +383357 +383533 +383853 +385205 +385469 +386029 +388625 +389027 +390437 +391075 +391332 +391520 +393397 +393906 +395389 +398151 +398261 +398460 +399329 +399601 +402816 +403551 +403990 +404522 +405987 +412442 +412630 +412656 +413271 +413467 +415170 +415344 +415773 +415936 +418239 +418388 +418701 +418778 +419633 +420532 +421043 +421263 +423522 +423915 +425898 +427400 +427946 +428117 +429464 +429593 +433933 +433983 +434484 +435212 +435348 +435604 +437148 +437283 +438641 +439831 +441009 +441581 +441979 +442138 +442804 +444584 +447019 +447986 +449235 +449525 +450820 +450822 +452296 +452971 +454420 +454661 +455744 +455876 +456950 +458972 +459948 +462615 +463811 +464004 +465753 +466349 +467139 +467531 +467584 +467710 +468484 +470558 +470682 +470934 +471145 +471826 +471955 +472864 +473529 +475202 +476465 +476753 +477252 +477384 +478311 +479818 +481452 +482138 +483000 +483973 +484116 +484217 +484872 +485865 +487222 +491172 +491999 +493057 +493696 +493824 +494459 +498480 +499571 +500049 +500346 +500445 +502323 +502950 +504119 +504713 +504941 +506016 +507663 +511945 +512633 +513989 +516960 +517294 +517801 +517834 +522321 +523270 +523637 +526204 +526431 +526893 +526906 +526973 +528611 +530070 +530178 +530407 +530639 +530947 +531128 +531275 +531590 +533010 +533197 +533276 +534678 +537643 +538678 +539939 +540030 +542130 +542633 +544457 +546530 +548646 +548717 +548805 +552227 +553607 +557117 +558482 +558904 +559903 +560888 +561368 +563245 +563294 +563520 +564397 +564819 +565582 +566792 +567812 +569815 +572072 +572878 +576465 +576587 +577557 +579867 +580260 +580682 +581478 +582940 +583940 +584145 +586959 +588493 +588532 +589990 +590776 +593155 +593430 +593632 +594885 +595898 +598620 +599420 +601715 +601872 +602178 +602591 +602611 +602888 +605984 +606798 +607131 +607931 +608085 +610722 +613000 +613276 +613315 +614909 +615460 +615779 +617764 +617963 +619214 +620203 +622414 +623978 +624052 +624286 +624650 +625527 +626657 +626716 +626825 +630264 +630760 +631079 +631264 +631575 +631733 +633212 +634957 +635967 +636516 +641406 +641606 +642197 +643329 +643397 +643632 +644194 +644308 +644785 +645815 +647888 +648320 +649396 +650704 +652522 +654155 +654313 +654704 +655820 +658800 +658832 +659286 +660245 +660659 +661129 +664010 +664975 +667304 +667599 +668089 +668210 +668477 +668915 +669689 +670024 +673294 +673510 +676008 +678502 +679143 +679305 +679762 +680229 +682241 +682719 +684165 +685298 +686964 +690741 +691328 +691402 +692105 +693619 +696799 +700232 +700421 +700540 +701239 +701297 +701498 +702498 +702840 +702962 +702971 +703363 +703402 +704854 +706082 +706428 +707107 +707957 +711738 +712288 +712397 +712757 +713260 +714453 +714743 +715834 +716660 +716975 +718168 +718628 +719645 +720116 +721325 +721367 +722486 +724389 +724722 +725961 +726835 +732226 +733141 +736511 +737244 +741208 +741645 +742280 +743755 +745060 +745957 +746212 +746979 +747513 +749077 +749419 +749989 +751348 +755830 +756270 +757168 +757597 +759115 +760318 +761602 +762269 +763085 +763215 +763560 +764092 +765464 +766537 +766679 +767822 +768656 +769114 +771565 +771988 +775845 +775887 +776480 +777256 +777739 +779012 +779863 +779897 +779912 +780362 +780999 +781265 +781550 +781796 +782246 +783556 +787035 +787654 +788261 +788442 +788591 +789326 +791016 +791183 +792031 +793039 +794748 +795684 +795955 +796421 +797516 +799256 +799777 +801081 +803588 +807751 +807914 +808120 +809399 +814791 +814922 +815033 +815346 +817360 +818647 +818951 +820052 +820449 +826627 +827412 +830216 +830482 +830580 +831032 +832097 +832245 +839261 +839381 +839789 +841781 +843469 +844115 +844932 +845283 +845317 +846331 +847041 +847373 +847699 +848528 +851434 +851465 +851466 +851774 +852268 +852789 +853654 +855316 +856162 +856547 +858407 +859504 +860219 +860654 +862215 +862307 +864936 +865908 +866458 +866753 +867911 +870193 +871661 +873418 +877441 +878009 +878255 +880540 +881472 +881988 +882393 +883930 +884255 +885117 +885118 +887174 +887953 +888765 +889304 +889352 +890564 +893397 +894275 +894633 +896104 +897298 +898800 +899154 +900886 +901107 +901859 +901993 +902375 +907005 +907141 +909813 +911359 +912036 +912699 +912766 +913740 +914224 +915912 +916088 +916528 +919320 +922633 +924463 +926970 +927283 +927983 +929035 +932962 +933185 +933715 +933810 +934681 +934725 +934727 +940356 +941238 +942410 +942480 +943924 +945225 +945769 +946655 +947763 +948215 +948388 +949448 +951693 +952402 +952635 +952872 +953753 +959636 +960626 +963786 +964389 +965439 +965750 +966015 +967207 +967844 +968522 +968848 +968990 +969958 +969985 +971247 +971737 +972461 +972976 +974358 +974404 +974628 +977057 +977819 +977995 +978427 +978718 +980127 +980330 +981505 +982017 +982332 +983276 +983368 +984405 +984453 +985561 +985918 +986581 +986698 +986809 +988366 +989661 +989866 +990988 +992476 +992937 +993197 +993759 +994253 +994818 +995242 diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index 1e1627f2..a77136ee 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -203,7 +203,7 @@ class Mineral: - `n_grains` (int) — number of grains in the aggregate - `fractions` (list of arrays) — grain volume fractions for each texture snapshot - `orientations` (list of arrays) — grain orientation matrices for each texture snapshot - - `rng` (`numpy.random.Generator`) — random number generator used for the isotropic + - `seed` (`int`) — seed used by the random number generator to set up the isotropic initial condition when `fractions_init` or `orientations_init` are not provided """ @@ -217,7 +217,7 @@ class Mineral: orientations_init: np.ndarray = None fractions: list = field(default_factory=list) orientations: list = field(default_factory=list) - rng: np.random.Generator = _stats._RNG + seed: int = None def __str__(self): # String output, used for str(self) and f"{self}", etc. @@ -253,7 +253,7 @@ def __post_init__(self): self.fractions_init = np.full(self.n_grains, 1.0 / self.n_grains) if self.orientations_init is None: self.orientations_init = Rotation.random( - self.n_grains, random_state=self.rng + self.n_grains, random_state=self.seed ).as_matrix() # Copy the initial values to the storage lists. diff --git a/src/pydrex/stats.py b/src/pydrex/stats.py index 466bc720..6ac42545 100644 --- a/src/pydrex/stats.py +++ b/src/pydrex/stats.py @@ -3,17 +3,16 @@ from pydrex import geometry as _geo -_RNG = np.random.default_rng(seed=8845) - -def resample_orientations(orientations, fractions, n_samples=None, rng=_RNG): +def resample_orientations(orientations, fractions, n_samples=None, seed=None): """Generate new samples from `orientations` weighed by the volume distribution. If the optional number of samples `n_samples` is not specified, it will be set to the number of original "grains" (length of `fractions`). - The argument `rng` can be used to specify a custom random number generator. + The argument `seed` can be used to seed the random number generator. """ + rng = np.random.default_rng(seed=seed) if n_samples is None: n_samples = len(fractions) sort_ascending = np.argsort(fractions) diff --git a/src/pydrex/visualisation.py b/src/pydrex/visualisation.py index 8bafee0b..b67e97fd 100644 --- a/src/pydrex/visualisation.py +++ b/src/pydrex/visualisation.py @@ -143,6 +143,7 @@ def simple_shear_stationary_2d( target_angles, angles, point100_symmetry, + angles_err=None, savefile="pydrex_simple_shear_stationary_2d.png", markers=("."), θ_fse=None, @@ -170,6 +171,10 @@ def simple_shear_stationary_2d( lines = ax_mean.plot(strains, θ_target, alpha=0.66, label=label) color = lines[0].get_color() ax_mean.plot(strains, θ, marker, markersize=5, alpha=0.33, color=color) + if angles_err is not None: + ax_mean.fill_between( + strains, θ - angles_err[i], θ + angles_err[i], alpha=0.22, color=color + ) ax_symmetry.plot( strains, point100, diff --git a/tests/conftest.py b/tests/conftest.py index bd232377..e17842df 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,15 +2,31 @@ import matplotlib import pytest from _pytest.logging import LoggingPlugin, _LiveLoggingStreamHandler -from numpy import random as rn from pydrex import logger as _log from pydrex import mock as _mock +from pydrex import io as _io matplotlib.use("Agg") # Stop matplotlib from looking for $DISPLAY in env. _log.quiet_aliens() # Stop imported modules from spamming the logs. +# Set up custom pytest CLI arguments. +def pytest_addoption(parser): + parser.addoption( + "--outdir", + metavar="DIR", + default=None, + help="output directory in which to store PyDRex figures/logs", + ) + parser.addoption( + "--runslow", + action="store_true", + default=False, + help="run slow tests (HPC cluster recommended, large memory requirement)", + ) + + # The default pytest logging plugin always creates its own handlers... class PytestConsoleLogger(LoggingPlugin): """Pytest plugin that allows linking up a custom console logger.""" @@ -34,10 +50,12 @@ def pytest_runtest_teardown(self, item): yield from self._runtest_for(item, "teardown") -# Hook up our logging plugin last, -# it relies on terminalreporter and capturemanager. @pytest.hookimpl(trylast=True) def pytest_configure(config): + config.addinivalue_line("markers", "slow: mark test as slow to run") + + # Hook up our logging plugin last, + # it relies on terminalreporter and capturemanager. if config.option.verbose > 0: terminal_reporter = config.pluginmanager.get_plugin("terminalreporter") capture_manager = config.pluginmanager.get_plugin("capturemanager") @@ -50,13 +68,13 @@ def pytest_configure(config): ) -def pytest_addoption(parser): - parser.addoption( - "--outdir", - metavar="DIR", - default=None, - help="output directory in which to store PyDRex figures/logs", - ) +def pytest_collection_modifyitems(config, items): + if config.getoption("--runslow"): + return # Don't skip slow tests. + skip_slow = pytest.mark.skip(reason="need --runslow option to run") + for item in items: + if "slow" in item.keywords: + item.add_marker(skip_slow) @pytest.fixture(scope="session") @@ -66,7 +84,7 @@ def outdir(request): @pytest.fixture(scope="function") def console_handler(request): - if request.config.option.verbose > 0: + if request.config.option.verbose > 0: # Show console logs if -v/--verbose given. return request.config.pluginmanager.get_plugin( "pytest-console-logger" ).log_cli_handler @@ -124,6 +142,6 @@ def ref_axes(request): @pytest.fixture -def rng(): - """A seeded RNG for tests to have (more) reproducible results.""" - return rn.default_rng(seed=8816) +def seeds(): + """1000 unique seeds for ensemble runs that need an RNG seed.""" + return _io.read_scsv(_io.data("rng") / "seeds.scsv").seeds diff --git a/tests/test_corner_flow_2d.py b/tests/test_corner_flow_2d.py index 70d1b5ac..7f94d5c2 100644 --- a/tests/test_corner_flow_2d.py +++ b/tests/test_corner_flow_2d.py @@ -81,7 +81,7 @@ class TestOlivineA: def test_corner_prescribed_init_isotropic( self, params_Kaminski2001_fig5_shortdash, - rng, + seed, outdir, ): """Test CPO evolution in prescribed 2D corner flow. @@ -96,7 +96,7 @@ def test_corner_prescribed_init_isotropic( domain_height = 2.0e5 # Represents the depth of olivine-spinel transition. domain_width = 1.0e6 n_grains = 2000 - orientations_init = Rotation.random(n_grains, random_state=rng).as_matrix() + orientations_init = Rotation.random(n_grains, random_state=seed).as_matrix() n_timesteps = 20 # Number of places along the pathline to compute CPO. # Z-values at the end of each pathline. z_ends = list(map(lambda x: x * domain_height, (-0.1, -0.3, -0.54, -0.78))) @@ -191,7 +191,7 @@ def _get_velocity_gradient(point): # Loop over first dimension (time steps) of orientations. for idx, matrices in enumerate(mineral.orientations): orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx], rng=rng + matrices, mineral.fractions[idx], seed=seed ) direction_mean = _diagnostics.bingham_average( orientations_resampled, diff --git a/tests/test_geometry.py b/tests/test_geometry.py index e6c24cf6..24189396 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -26,7 +26,7 @@ def test_poles_example(hkl, ref_axes): np.testing.assert_allclose(ref_data["zvals"], zvals, atol=1e-16, rtol=0) -def test_lambert_equal_area(rng): +def test_lambert_equal_area(seed): """Test Lambert equal area projection.""" x, y = np.mgrid[-1:1:11j, -1:1:11j] x_flat = [j for i in x for j in i] @@ -35,6 +35,7 @@ def test_lambert_equal_area(rng): x_disk, y_disk = _geo.shirley_concentric_squaredisk(x_flat, y_flat) # Project onto the unit sphere by adding z = ± (1 - r). # Then project back onto the disk using Lambert equal-area, should be the same. + rng = np.random.default_rng(seed=seed) sign = rng.integers(low=0, high=2, size=len(x_disk)) x_laea, y_laea = _geo.lambert_equal_area( x_disk, y_disk, (-1) ** sign * (1 - (x_disk**2 + y_disk**2)) diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index e5030371..46ede930 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -1,8 +1,9 @@ """> PyDRex: 2D simple shear tests.""" import contextlib as cl +import pytest import numpy as np -from numpy import testing as nt +# from numpy import testing as nt from scipy.interpolate import PchipInterpolator from pydrex import diagnostics as _diagnostics @@ -25,12 +26,13 @@ class TestOlivineA: def get_position(self, t): return np.zeros(3) - def test_dvdx_GBM( + @pytest.mark.slow + def test_dvdx_GBM_ensemble( self, params_Kaminski2001_fig5_solid, # GBM = 0 params_Kaminski2001_fig5_shortdash, # GBM = 50 params_Kaminski2001_fig5_longdash, # GBM = 200 - rng, + seeds, outdir, ): r"""Test a-axis alignment to shear in the Y direction. @@ -44,9 +46,9 @@ def test_dvdx_GBM( """ strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. timestamps = np.linspace(0, 2e5, 201) # Solve until D₀t=1 ('shear strain' γ=2). - n_timesteps = len(timestamps) - i_first_cpo = 50 # First index where Bingham averages are sufficiently stable. - i_strain_50p = [0, 50, 100, 150, 200] # Indices for += 50% strain. + n_timestamps = len(timestamps) # Number of steps + 1. + # i_first_cpo = 50 # First index where Bingham averages are sufficiently stable. + # i_strain_50p = [0, 50, 100, 150, 200] # Indices for += 50% strain. def get_velocity_gradient(x): # It is independent of time or position in this test. @@ -57,12 +59,12 @@ def get_velocity_gradient(x): shear_direction = [0, 1, 0] # Used to calculate the angular diagnostics. # Output setup with optional logging and data series labels. - θ_fse = [45] - angles = [] - point100_symmetry = [] + θ_fse = np.empty((len(seeds), n_timestamps)) + angles = np.empty((3, len(seeds), n_timestamps)) + point100_symmetry = np.empty_like(angles) optional_logging = cl.nullcontext() if outdir is not None: - out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_GBM" + out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dvdx_GBM_ensemble" optional_logging = _log.logfile_enable(f"{out_basepath}.log") labels = [] @@ -74,63 +76,73 @@ def get_velocity_gradient(x): params_Kaminski2001_fig5_longdash, # GBM = 200 ), ): - mineral = _minerals.Mineral( - n_grains=params["number_of_grains"], rng=rng - ) - deformation_gradient = np.eye(3) # Undeformed initial state. - for t, time in enumerate(timestamps[:-1], start=1): - _log.info( - "step %s/%s (t=%s) with velocity gradient: %s", - t, - n_timesteps - 1, - time, - get_velocity_gradient(None).flatten(), + for s, seed in enumerate(seeds): + mineral = _minerals.Mineral( + n_grains=params["number_of_grains"], seed=seed ) - deformation_gradient = mineral.update_orientations( - params, - deformation_gradient, - get_velocity_gradient, - pathline=(time, timestamps[t], self.get_position), - ) - fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) - _log.info("› strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) - if p == 0: - θ_fse.append( - _diagnostics.smallest_angle(fse_v, shear_direction) + deformation_gradient = np.eye(3) # Undeformed initial state. + θ_fse[s, 0] = 45 + for t, time in enumerate(timestamps[:-1], start=1): + _log.info( + "$M^∗=%s$; #%s/%s; step %s/%s (t=%s)", + params["gbm_mobility"], + s + 1, + len(seeds), + t, + n_timestamps - 1, + time, ) - - texture_symmetry = np.zeros(n_timesteps) - mean_angles = np.zeros(n_timesteps) - # Loop over first dimension (time steps) of orientations. - for idx, matrices in enumerate(mineral.orientations): - orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx] - ) - direction_mean = _diagnostics.bingham_average( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - ) - mean_angles[idx] = _diagnostics.smallest_angle( - direction_mean, shear_direction + deformation_gradient = mineral.update_orientations( + params, + deformation_gradient, + get_velocity_gradient, + pathline=(time, timestamps[t], self.get_position), + ) + _, fse_v = _diagnostics.finite_strain(deformation_gradient) + _log.info( + "› velocity gradient = %s\n› strain D₀t=%s", + get_velocity_gradient(None).flatten(), + strain_rate * time, + ) + if p == 0: + θ_fse[s, t] = _diagnostics.smallest_angle( + fse_v, shear_direction + ) + + texture_symmetry = np.zeros(n_timestamps) + # mean_angles = np.zeros(n_timestamps) + # Loop over first dimension (time steps) of orientations. + for idx, matrices in enumerate(mineral.orientations): + orientations_resampled, _ = _stats.resample_orientations( + matrices, mineral.fractions[idx], seed=seed + ) + # direction_mean = _diagnostics.bingham_average( + # orientations_resampled, + # axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + # ) + # mean_angles[idx] = _diagnostics.smallest_angle( + # direction_mean, shear_direction + # ) + texture_symmetry[idx] = _diagnostics.symmetry( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + )[0] + + # Use SCCS axis (hexagonal symmetry) for the angle instead (opt). + mean_angles = np.array( + [ + _diagnostics.smallest_angle( + _diagnostics.anisotropy(v)[1][2, :], shear_direction + ) + for v in _minerals.voigt_averages([mineral], params) + ] ) - texture_symmetry[idx] = _diagnostics.symmetry( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - )[0] - - # Uncomment to use SCCS axis (hexagonal symmetry) for the angle instead. - # mean_angles = np.array( - # [ - # _diagnostics.smallest_angle( - # _diagnostics.anisotropy(v)[1][2, :], shear_direction - # ) - # for v in _minerals.voigt_averages([mineral], params) - # ] - # ) - - # Store data series and optional plotting metadata. - angles.append(mean_angles) - point100_symmetry.append(texture_symmetry) + + # Store data series and optional plotting metadata. + angles[p, s, :] = mean_angles + point100_symmetry[p, s, :] = texture_symmetry + + # Update labels and store the last mineral of the ensemble for polefigs. if outdir is not None: labels.append(f"$M^∗$ = {params['gbm_mobility']}") mineral.save( @@ -156,7 +168,11 @@ def get_velocity_gradient(x): ) target_angles = [cs_M0(strains), cs_M50(strains), cs_M200(strains)] - # Optionally plot figure. + # Take ensemble means and optionally plot figure. + result_angles = angles.mean(axis=1) + result_angles_err = angles.std(axis=1) + result_point100_symmetry = point100_symmetry.mean(axis=1) + result_θ_fse = θ_fse.mean(axis=0) if outdir is not None: schema = { "delimiter": ",", @@ -178,61 +194,75 @@ def get_velocity_gradient(x): _vis.simple_shear_stationary_2d( strains, target_angles, - angles, - point100_symmetry, + result_angles, + result_point100_symmetry, + angles_err=result_angles_err, savefile=f"{out_basepath}.png", markers=("o", "v", "s"), - θ_fse=θ_fse, + θ_fse=result_θ_fse, labels=labels, ) - # Check that FSE is correct. - # First, get theoretical FSE axis for simple shear. - # We want the angle from the Y axis (shear direction), so subtract from 90. - θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] - nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) - - # Check Bingham angles, ignoring the first portion. - # Average orientations of near-isotropic distributions are unstable. - nt.assert_allclose( - θ_fse[i_first_cpo:], angles[0][i_first_cpo:], rtol=0.11, atol=0 - ) - nt.assert_allclose( - target_angles[0][i_first_cpo:], angles[0][i_first_cpo:], rtol=0.11, atol=0 - ) - nt.assert_allclose( - target_angles[1][i_first_cpo:], angles[1][i_first_cpo:], rtol=0, atol=5.7 - ) - nt.assert_allclose( - target_angles[2][i_first_cpo:], angles[2][i_first_cpo:], rtol=0, atol=5.5 - ) - - # Check point symmetry of [100] at strains of 0%, 50%, 100%, 150% & 200%. - nt.assert_allclose( - [0.015, 0.11, 0.19, 0.27, 0.34], - point100_symmetry[0].take(i_strain_50p), - rtol=0, - atol=0.015, - ) - nt.assert_allclose( - [0.015, 0.15, 0.33, 0.57, 0.72], - point100_symmetry[1].take(i_strain_50p), - rtol=0, - atol=0.015, - ) - nt.assert_allclose( - [0.015, 0.22, 0.64, 0.86, 0.91], - point100_symmetry[2].take(i_strain_50p), - rtol=0, - atol=0.015, - ) - - def test_dudz_GBS( + # # Check that FSE is correct. + # # First, get theoretical FSE axis for simple shear. + # # We want the angle from the Y axis (shear direction), so subtract from 90. + # θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] + # nt.assert_allclose(result_θ_fse, θ_fse_eq, rtol=1e-7, atol=0) + + # # Check Bingham angles, ignoring the first portion. + # # Average orientations of near-isotropic distributions are unstable. + # nt.assert_allclose( + # result_θ_fse[i_first_cpo:], + # result_angles[0][i_first_cpo:], + # rtol=0.11, + # atol=0, + # ) + # nt.assert_allclose( + # target_angles[0][i_first_cpo:], + # result_angles[0][i_first_cpo:], + # rtol=0.11, + # atol=0, + # ) + # nt.assert_allclose( + # target_angles[1][i_first_cpo:], + # result_angles[1][i_first_cpo:], + # rtol=0, + # atol=5.7, + # ) + # nt.assert_allclose( + # target_angles[2][i_first_cpo:], + # result_angles[2][i_first_cpo:], + # rtol=0, + # atol=5.5, + # ) + + # # Check point symmetry of [100] at strains of 0%, 50%, 100%, 150% & 200%. + # nt.assert_allclose( + # [0.015, 0.11, 0.19, 0.27, 0.34], + # result_point100_symmetry[0].take(i_strain_50p), + # rtol=0, + # atol=0.015, + # ) + # nt.assert_allclose( + # [0.015, 0.15, 0.33, 0.57, 0.72], + # result_point100_symmetry[1].take(i_strain_50p), + # rtol=0, + # atol=0.015, + # ) + # nt.assert_allclose( + # [0.015, 0.22, 0.64, 0.86, 0.91], + # result_point100_symmetry[2].take(i_strain_50p), + # rtol=0, + # atol=0.015, + # ) + + @pytest.mark.slow + def test_dudz_GBS_ensemble( self, params_Kaminski2004_fig4_circles, # GBS = 0 params_Kaminski2004_fig4_squares, # GBS = 0.2 params_Kaminski2004_fig4_triangles, # GBS = 0.4 - rng, + seeds, outdir, ): r"""Test a-axis alignment to shear in X direction. @@ -243,7 +273,7 @@ def test_dudz_GBS( """ strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. timestamps = np.linspace(0, 5e5, 251) # Solve until D₀t=2.5 ('shear' γ=5). - n_timesteps = len(timestamps) + n_timestamps = len(timestamps) i_strain_100p = [0, 50, 100, 150, 200] # Indices for += 100% strain. def get_velocity_gradient(x): @@ -255,12 +285,12 @@ def get_velocity_gradient(x): shear_direction = [1, 0, 0] # Used to calculate the angular diagnostics. # Output setup with optional logging and data series labels. - θ_fse = [45] - angles = [] - point100_symmetry = [] + θ_fse = np.empty((len(seeds), n_timestamps)) + angles = np.empty((3, len(seeds), n_timestamps)) + point100_symmetry = np.empty_like(angles) optional_logging = cl.nullcontext() if outdir is not None: - out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dudz_GBS" + out_basepath = f"{outdir}/{SUBDIR}/{self.class_id}_dudz_GBS_ensemble" optional_logging = _log.logfile_enable(f"{out_basepath}.log") labels = [] @@ -272,63 +302,73 @@ def get_velocity_gradient(x): params_Kaminski2004_fig4_triangles, # GBS = 0.4 ), ): - mineral = _minerals.Mineral( - n_grains=params["number_of_grains"], rng=rng - ) - deformation_gradient = np.eye(3) # Undeformed initial state. - for t, time in enumerate(timestamps[:-1], start=1): - _log.info( - "step %s/%s (t=%s) with velocity gradient: %s", - t, - n_timesteps - 1, - time, - get_velocity_gradient(None).flatten(), + for s, seed in enumerate(seeds): + mineral = _minerals.Mineral( + n_grains=params["number_of_grains"], seed=seed ) - deformation_gradient = mineral.update_orientations( - params, - deformation_gradient, - get_velocity_gradient, - pathline=(time, timestamps[t], self.get_position), - ) - fse_λ, fse_v = _diagnostics.finite_strain(deformation_gradient) - _log.info("› strain √λ-1=%s (D₀t=%s)", fse_λ, strain_rate * time) - if p == 0: - θ_fse.append( - _diagnostics.smallest_angle(fse_v, shear_direction) + deformation_gradient = np.eye(3) # Undeformed initial state. + θ_fse[s, 0] = 45 + for t, time in enumerate(timestamps[:-1], start=1): + _log.info( + "X = %s; # %s/%s; step %s/%s (t = %s)", + params["gbs_threshold"], + s + 1, + len(seeds), + t, + n_timestamps - 1, + time, ) - - texture_symmetry = np.zeros(n_timesteps) - mean_angles = np.zeros(n_timesteps) - # Loop over first dimension (time steps) of orientations. - for idx, matrices in enumerate(mineral.orientations): - orientations_resampled, _ = _stats.resample_orientations( - matrices, mineral.fractions[idx] - ) - direction_mean = _diagnostics.bingham_average( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - ) - mean_angles[idx] = _diagnostics.smallest_angle( - direction_mean, shear_direction + deformation_gradient = mineral.update_orientations( + params, + deformation_gradient, + get_velocity_gradient, + pathline=(time, timestamps[t], self.get_position), + ) + _, fse_v = _diagnostics.finite_strain(deformation_gradient) + _log.info( + "› velocity gradient = %s", + get_velocity_gradient(None).flatten(), + ) + _log.info("› strain D₀t = %s", strain_rate * time) + if p == 0: + θ_fse[s, t] = _diagnostics.smallest_angle( + fse_v, shear_direction + ) + + texture_symmetry = np.zeros(n_timestamps) + # mean_angles = np.zeros(n_timestamps) + # Loop over first dimension (time steps) of orientations. + for idx, matrices in enumerate(mineral.orientations): + orientations_resampled, _ = _stats.resample_orientations( + matrices, mineral.fractions[idx], seed=seed + ) + # direction_mean = _diagnostics.bingham_average( + # orientations_resampled, + # axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + # ) + # mean_angles[idx] = _diagnostics.smallest_angle( + # direction_mean, shear_direction + # ) + texture_symmetry[idx] = _diagnostics.symmetry( + orientations_resampled, + axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], + )[0] + + # Use SCCS axis (hexagonal symmetry) for the angle instead (opt). + mean_angles = np.array( + [ + _diagnostics.smallest_angle( + _diagnostics.anisotropy(v)[1][2, :], shear_direction + ) + for v in _minerals.voigt_averages([mineral], params) + ] ) - texture_symmetry[idx] = _diagnostics.symmetry( - orientations_resampled, - axis=_minerals.OLIVINE_PRIMARY_AXIS[mineral.fabric], - )[0] - - # Uncomment to use SCCS axis (hexagonal symmetry) for the angle instead. - # mean_angles = np.array( - # [ - # _diagnostics.smallest_angle( - # _diagnostics.anisotropy(v)[1][2, :], shear_direction - # ) - # for v in _minerals.voigt_averages([mineral], params) - # ] - # ) - - # Store data series and optional plotting metadata. - angles.append(mean_angles) - point100_symmetry.append(texture_symmetry) + + # Store data series and optional plotting metadata. + angles[p, s, :] = mean_angles + point100_symmetry[p, s, :] = texture_symmetry + + # Update labels and store the last mineral of the ensemble for polefigs. if outdir is not None: labels.append(f"$f_{{gbs}}$ = {params['gbs_threshold']}") mineral.save( @@ -354,7 +394,11 @@ def get_velocity_gradient(x): ) target_angles = [cs_X0(strains), cs_X0d2(strains), cs_X0d4(strains)] - # Optionally plot figure. + # Take ensemble means and optionally plot figure. + result_angles = angles.mean(axis=1) + result_angles_err = angles.std(axis=1) + result_point100_symmetry = point100_symmetry.mean(axis=1) + result_θ_fse = θ_fse.mean(axis=0) if outdir is not None: schema = { "delimiter": ",", @@ -376,35 +420,37 @@ def get_velocity_gradient(x): _vis.simple_shear_stationary_2d( strains, target_angles, - angles, - point100_symmetry, + result_angles, + result_point100_symmetry, + angles_err=result_angles_err, savefile=f"{out_basepath}.png", markers=("o", "v", "s"), + θ_fse=result_θ_fse, labels=labels, ) - # Check that FSE is correct. - # First, get theoretical FSE axis for simple shear. - # We want the angle from the Y axis (shear direction), so subtract from 90. - θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] - nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) - - # Check point symmetry of [100] at strains of 0%, 100%, 200%, 300% & 400%. - nt.assert_allclose( - [0.015, 0.52, 0.86, 0.93, 0.94], - point100_symmetry[0].take(i_strain_100p), - rtol=0, - atol=0.015, - ) - nt.assert_allclose( - [0.015, 0.42, 0.71, 0.77, 0.79], - point100_symmetry[1].take(i_strain_100p), - rtol=0, - atol=0.015, - ) - nt.assert_allclose( - [0.015, 0.36, 0.57, 0.6, 0.62], - point100_symmetry[2].take(i_strain_100p), - rtol=0, - atol=0.015, - ) + # # Check that FSE is correct. + # # First, get theoretical FSE axis for simple shear. + # # We want the angle from the Y axis (shear direction), so subtract from 90. + # θ_fse_eq = [90 - _utils.angle_fse_simpleshear(strain) for strain in strains] + # nt.assert_allclose(θ_fse, θ_fse_eq, rtol=1e-7, atol=0) + + # # Check point symmetry of [100] at strains of 0%, 100%, 200%, 300% & 400%. + # nt.assert_allclose( + # [0.015, 0.52, 0.86, 0.93, 0.94], + # point100_symmetry[0].take(i_strain_100p), + # rtol=0, + # atol=0.015, + # ) + # nt.assert_allclose( + # [0.015, 0.42, 0.71, 0.77, 0.79], + # point100_symmetry[1].take(i_strain_100p), + # rtol=0, + # atol=0.015, + # ) + # nt.assert_allclose( + # [0.015, 0.36, 0.57, 0.6, 0.62], + # point100_symmetry[2].take(i_strain_100p), + # rtol=0, + # atol=0.015, + # ) From c578c4cb37240aa2774d36e024657f2910228a08 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 24 Jul 2023 01:44:14 +0000 Subject: [PATCH 27/30] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/conftest.py | 2 +- tests/test_simple_shear_2d.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index e17842df..0acab459 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -3,9 +3,9 @@ import pytest from _pytest.logging import LoggingPlugin, _LiveLoggingStreamHandler +from pydrex import io as _io from pydrex import logger as _log from pydrex import mock as _mock -from pydrex import io as _io matplotlib.use("Agg") # Stop matplotlib from looking for $DISPLAY in env. _log.quiet_aliens() # Stop imported modules from spamming the logs. diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index 46ede930..b902bd71 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -1,8 +1,9 @@ """> PyDRex: 2D simple shear tests.""" import contextlib as cl -import pytest import numpy as np +import pytest + # from numpy import testing as nt from scipy.interpolate import PchipInterpolator From 95b67e494f91589a15b56abb46752c2b7b999366 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 24 Jul 2023 11:49:53 +1000 Subject: [PATCH 28/30] tests: Add missing seed fixture --- tests/conftest.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/conftest.py b/tests/conftest.py index 0acab459..6679cbe7 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -145,3 +145,8 @@ def ref_axes(request): def seeds(): """1000 unique seeds for ensemble runs that need an RNG seed.""" return _io.read_scsv(_io.data("rng") / "seeds.scsv").seeds + +@pytest.fixture +def seed(): + """Default seed for test RNG.""" + return 8816 From 43cf3dc46b6106fe15bcf7ba8faf958a7222cae6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 24 Jul 2023 01:51:30 +0000 Subject: [PATCH 29/30] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/conftest.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/conftest.py b/tests/conftest.py index 6679cbe7..7519e103 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -146,6 +146,7 @@ def seeds(): """1000 unique seeds for ensemble runs that need an RNG seed.""" return _io.read_scsv(_io.data("rng") / "seeds.scsv").seeds + @pytest.fixture def seed(): """Default seed for test RNG.""" From a7876c845a9c20a7c8acdd099ce0ccb0cdcd73f8 Mon Sep 17 00:00:00 2001 From: adigitoleo Date: Mon, 24 Jul 2023 11:53:02 +1000 Subject: [PATCH 30/30] fix: Fix pre-commit issues --- src/pydrex/minerals.py | 1 - tests/test_simple_shear_2d.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/pydrex/minerals.py b/src/pydrex/minerals.py index a77136ee..8276a8bd 100644 --- a/src/pydrex/minerals.py +++ b/src/pydrex/minerals.py @@ -12,7 +12,6 @@ from pydrex import exceptions as _err from pydrex import io as _io from pydrex import logger as _log -from pydrex import stats as _stats from pydrex import tensors as _tensors OLIVINE_STIFFNESS = np.array( diff --git a/tests/test_simple_shear_2d.py b/tests/test_simple_shear_2d.py index b902bd71..d82ecb89 100644 --- a/tests/test_simple_shear_2d.py +++ b/tests/test_simple_shear_2d.py @@ -275,7 +275,7 @@ def test_dudz_GBS_ensemble( strain_rate = 5e-6 # Strain rate from Fraters & Billen, 2021, fig. 3. timestamps = np.linspace(0, 5e5, 251) # Solve until D₀t=2.5 ('shear' γ=5). n_timestamps = len(timestamps) - i_strain_100p = [0, 50, 100, 150, 200] # Indices for += 100% strain. + # i_strain_100p = [0, 50, 100, 150, 200] # Indices for += 100% strain. def get_velocity_gradient(x): # It is independent of time or position in this test.