Skip to content

Commit

Permalink
fix: GBS now carries over previous orientation
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
adigitoleo committed Jul 20, 2023
1 parent bd8b1d7 commit 058f881
Show file tree
Hide file tree
Showing 3 changed files with 27 additions and 31 deletions.
3 changes: 1 addition & 2 deletions src/pydrex/minerals.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
23 changes: 0 additions & 23 deletions src/pydrex/visualisation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -214,7 +192,6 @@ def _lag_2d_corner_flow(θ):
)


@check_marker_seq
def corner_flow_2d(
x_paths,
z_paths,
Expand Down
32 changes: 26 additions & 6 deletions tests/test_simple_shear_2d.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
):
Expand All @@ -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.
Expand All @@ -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(
Expand Down Expand Up @@ -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,
)

0 comments on commit 058f881

Please sign in to comment.