diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md new file mode 100644 index 00000000..8ded361b --- /dev/null +++ b/.github/copilot-instructions.md @@ -0,0 +1,31 @@ +# File formatting + +Use LF (`\n`) line endings, not CRLF (`\r\n`), in every text file you create, +including on Windows. + +# Testing + +When running pytest directly, always load the project's pytest plugin: + +```console +python -m pytest -p matflow.pytest_plugin matflow/tests/unit/test_cli.py -q +``` + +Use the configured Python interpreter and select the tests relevant to the change. +The plugin registers options required by the test fixtures, including `--slurm` +and `--repeat`. Without it, collection fails. The `matflow test` command loads +the plugin automatically. + +When running integration tests, always pass `--configure-python-env` to configure +the matflow Python environment using the active Python environment: + +```console +matflow test --integration --configure-python-env +``` + +When invoking pytest directly, include both `-p matflow.pytest_plugin` and +`--configure-python-env`. + +Before running integration tests, read +`.github/instructions/local-testing.instructions.md` if it exists for +machine-specific test setup. This optional file is local and must not be committed. diff --git a/.github/instructions/python-typechecking.instructions.md b/.github/instructions/python-typechecking.instructions.md new file mode 100644 index 00000000..0bca84ea --- /dev/null +++ b/.github/instructions/python-typechecking.instructions.md @@ -0,0 +1,12 @@ +--- +description: "Use when adding or modifying Python code, including tests. Requires mypy validation before completing changes." +applyTo: "**/*.py" +--- + +# Python type checking + +- Run mypy on new or modified Python files, including tests, before completing a code change. +- Use the configured project environment and existing mypy configuration. For focused validation, run `poetry run mypy ` from the repository root. +- Fix type errors introduced by the change with appropriate annotations and narrowing, not blanket ignores or unnecessary casts. +- For changes affecting shared types or APIs, also run the CI command: `poetry run mypy hpcflow`. +- If validation is blocked or pre-existing errors remain, report the exact limitation rather than claiming the type check passed. diff --git a/docs/source/user/how_to/task_schemas.rst b/docs/source/user/how_to/task_schemas.rst index 6a850fde..6bf42ba1 100644 --- a/docs/source/user/how_to/task_schemas.rst +++ b/docs/source/user/how_to/task_schemas.rst @@ -1,6 +1,170 @@ Task schemas ------------ +Generating multiple load cases +============================== + +The built-in ``generate_load_cases`` task generates a list of independent uniaxial +``LoadCase`` objects with identical loading parameters: + +.. code-block:: yaml + + name: multiple_load_cases + tasks: + - schema: generate_load_cases_uniaxial + inputs: + num_load_cases: 3 + total_time: 100 + num_increments: 200 + direction: x + target_def_grad_rate: 1.0e-3 + +``num_load_cases`` must be a positive integer. The ``load_case`` schema output has +``multiplicity: "<>"``, allowing downstream tasks to source +the individual load cases. With a multiplicity of one, the output remains a +single-member list, in accordance with the framework's multiplicity semantics. + +Sampling EBSD textures with Python MTEX +====================================== + +The ``sample_texture_from_CTF_file_mtex_python`` and +``sample_texture_from_CRC_file_mtex_python`` schemas estimate an orientation density +from one EBSD phase and sample ``num_orientations`` orientations from it. They require +the Python edition of MTEX (``mtex>=1.0``) in ``python_env``; MATLAB is not needed. +The existing ``*_mtex`` schemas still use MATLAB. + +.. code-block:: yaml + + tasks: + - schema: sample_texture_from_CTF_file_mtex_python + inputs: + CTF_file_path: C:\data\texture.ctf + specimen_symmetry: '1' + EBSD_phase: Al + num_orientations: 1000 + num_orientation_sets: 3 + kernel_halfwidth_deg: 10.0 + +For CRC data, use ``sample_texture_from_CRC_file_mtex_python`` and +``CRC_file_path`` instead. The matching ``.cpr`` file must be beside the ``.crc`` +file, with the same stem. + +``num_orientation_sets`` defaults to 1 and must be a positive integer. The density +is estimated once, then independently sampled for each set using the run's seeded +random generator. Each set contains ``num_orientations`` orientations. +The ``orientations`` output has multiplicity ``<>``: +its value is a list of orientation sets, including when the count is 1. +Downstream tasks sourcing ``orientations`` receive individual sets. + +All inputs and the ``orientations`` output are passed directly to/from the script, +without JSON/HDF5 handoff files. There is no ``execution_mode`` input. +Sampling uses ``MATFLOW_RUN_RANDOM_SEED`` for reproducibility within Python MTEX; +samples are not expected to match MATLAB bit for bit. ``kernel_halfwidth_deg`` +sets the density-estimation kernel halfwidth in degrees and defaults to ``10.0``, +preserving the previous smoothing choice. It must be finite, greater than zero, +and at most 180 degrees. The halfwidth is passed explicitly to MTEX, so it no +longer warns about selecting a default for correlated EBSD pixels. + +``EBSD_reference_frame_transformation`` defaults to ``euler_to_spatial`` (also +accepted as ``convertEuler2SpatialReferenceFrame``), using the reader's Oxford +Euler correction. ``spatial_to_euler`` (or +``convertSpatial2EulerReferenceFrame``) leaves the Euler angles unchanged, as does +a null transformation: spatial coordinates are not part of this task's output. +The distinction matters because passing ``None`` as the reader's correction asks +MTEX to apply its default Oxford correction, while an explicit identity rotation +leaves the file's Euler angles unchanged. + +``EBSD_rotation`` optionally applies a rigid map rotation specified by three +``euler_angles_deg``. By default, both map positions and orientations are rotated. +``keep_XY: true`` keeps map positions fixed and rotates orientations; +``keep_euler: true`` keeps orientations fixed and rotates map positions. The flags +cannot both be true. Since this task samples only orientations, ``keep_euler: true`` +does not change the sampled texture. + +The returned orientations use scalar-first inverse MTEX quaternions and the +crystal frame's unit-cell alignment. A frame whose Cartesian axes cannot all be +described by MatFlow lattice directions is rejected explicitly. + +MATLAB MTEX compatibility +~~~~~~~~~~~~~~~~~~~~~~~~ + +The ``sample_texture_from_CTF_file_mtex`` and +``sample_texture_from_CRC_file_mtex`` schemas accept the same sampling inputs, +including ``kernel_halfwidth_deg`` with the same default and valid range, +reference-frame aliases, and rotation flags. ``spatial_to_euler`` (or +``convertSpatial2EulerReferenceFrame``) leaves Euler angles unchanged in both +implementations; ``euler_to_spatial`` applies the Oxford correction. They retain +JSON input/HDF5 output and the ``script``, ``compile``, and ``precompiled`` +execution modes. Their output continues to use scalar-first inverse MTEX +quaternions and crystal alignment. Both implementations use the run's random +seed, but do not produce identical samples. + +The MATLAB schemas sample one orientation set per task. ``num_orientation_sets`` +is currently supported only by the Python MTEX schemas. + +Both implementations read ``MATFLOW_RUN_NUM_THREADS`` from the run's +``num_threads`` resource (default 1). Python sets MTEX's ``threads`` preference; +MATLAB sets its computational thread limit using ``maxNumCompThreads``. MATLAB +MTEX's separate NFFT/MEX libraries may have their own threading limits. + +Use ``execution_mode: script`` or recompile with ``execution_mode: compile`` to +use the updated MATLAB behavior. Existing precompiled executables must be rebuilt +from the updated scripts; changing the schema does not update those binaries. + +Plotting orientations with Python MTEX +===================================== + +The ``visualise_orientations_pole_figure_mtex_python`` and +``visualise_orientations_inverse_pole_figure_mtex_python`` schemas are Python +alternatives to the corresponding MATLAB ``*_mtex`` schemas. They require Python +MTEX (``mtex>=1.0``) in ``python_env`` and pass inputs directly to the script, +without JSON/HDF5 handoff files or an ``execution_mode`` input. + +For example, append these tasks after an orientation-sampling task: + +.. code-block:: yaml + + - schema: visualise_orientations_pole_figure_mtex_python + inputs: + crystal_symmetry: cubic + pole_figure_directions: [[0, 0, 1], [1, 0, 1], [1, 1, 1]] + use_contours: false + IPF_reference_direction: z + + - schema: visualise_orientations_inverse_pole_figure_mtex_python + inputs: + crystal_symmetry: cubic + IPF_reference_directions: [x, y, z] + use_contours: true + +Each task consumes a single ``orientations`` value. When a sampling task produces +multiple orientation sets, output multiplicity creates one plotting element per +set automatically. + +Both quaternion component orders are supported. MatFlow's inverse quaternions +are converted back to MTEX's crystal-to-specimen convention, with the input +unit-cell alignment preserved. Euler-angle inputs are not supported and raise +an explicit error. As in the MATLAB tasks, the crystal frame is constructed from +``crystal_symmetry`` and the alignment, using MTEX's default lattice parameters. + +``use_contours: false`` draws scatter plots; pole-figure markers use IPF colours +for ``IPF_reference_direction`` (x, y or z, default z). Inverse-pole-figure markers +are uncoloured, as in the MATLAB task, with one panel per +``IPF_reference_directions`` entry (default ``[z]``). Both scatter tasks also save +an IPF colour key. ``use_contours: true`` instead draws filled density contours +and colourbars, without a colour-key image. For pole figures, +``colourbar_limits`` optionally supplies two finite, increasing values, and +``use_one_colourbar: true`` uses a shared colour scale and a single bar. These +colourbar options apply only to contour plots. + +The files are saved as ``pole_figure.png`` or ``inverse_pole_figure.png`` and, +for scatter plots, ``IPF_key.png``. The existing file labels +``mtex_pole_figures``, ``mtex_inverse_pole_figures`` and ``mtex_IPF_key`` are +retained, so the images are saved as workflow artifacts just as for MATLAB. +Rendering is headless, uses x pointing east, and closes its figures after saving. +``MATFLOW_RUN_NUM_THREADS`` controls MTEX's thread preference, as for Python +texture sampling. Figure styling and density contours may differ from MATLAB. + Action command strings ====================== diff --git a/matflow/data/data_manifests/data.json b/matflow/data/data_manifests/data.json index daf6cc1f..2246bcc7 100644 --- a/matflow/data/data_manifests/data.json +++ b/matflow/data/data_manifests/data.json @@ -46,5 +46,8 @@ }, "Zr_cooling_data": { "zip_file": "Zr_cooling_data.zip" - } + }, + "posterior.txt": { + "zip_file": "posterior.zip" + } } diff --git a/matflow/data/scripts/generate_load_case.py b/matflow/data/scripts/generate_load_case.py new file mode 100644 index 00000000..789c6ff4 --- /dev/null +++ b/matflow/data/scripts/generate_load_case.py @@ -0,0 +1,20 @@ +from __future__ import annotations + +from matflow.param_classes.load import LoadCase + + +def generate_load_case( + total_time: float | int, + num_increments: int, + direction: str, + target_def_grad_rate: float, +) -> dict[str, LoadCase]: + """Generate a single-step uniaxial load case.""" + return { + "load_case": LoadCase.uniaxial( + total_time=total_time, + num_increments=num_increments, + direction=direction, + target_def_grad_rate=target_def_grad_rate, + ) + } diff --git a/matflow/data/scripts/generate_load_cases.py b/matflow/data/scripts/generate_load_cases.py new file mode 100644 index 00000000..edb088bd --- /dev/null +++ b/matflow/data/scripts/generate_load_cases.py @@ -0,0 +1,31 @@ +from __future__ import annotations + +from matflow.data.scripts.generate_load_case import generate_load_case +from matflow.param_classes.load import LoadCase + + +def generate_load_cases( + num_load_cases: int, + total_time: float | int, + num_increments: int, + direction: str, + target_def_grad_rate: float, +) -> dict[str, list[LoadCase]]: + """Generate independent uniaxial load cases with the same loading parameters.""" + if ( + isinstance(num_load_cases, bool) + or not isinstance(num_load_cases, int) + or num_load_cases < 1 + ): + raise ValueError("num_load_cases must be a positive integer.") + return { + "load_case": [ + generate_load_case( + total_time=total_time + idx, + num_increments=num_increments, + direction=direction, + target_def_grad_rate=target_def_grad_rate, + )["load_case"] + for idx in range(num_load_cases) + ] + } diff --git a/matflow/data/scripts/generate_load_cases_random_exponential.py b/matflow/data/scripts/generate_load_cases_random_exponential.py new file mode 100644 index 00000000..85996240 --- /dev/null +++ b/matflow/data/scripts/generate_load_cases_random_exponential.py @@ -0,0 +1,172 @@ +from __future__ import annotations +import os +from typing import Literal +import warnings + +import numpy as np +from numpy.typing import ArrayLike + +from matflow.param_classes.load import LoadCase, LoadStep + + +def generate_load_cases_random_exponential( + num_load_cases: int, + total_time: float | int, + num_steps: int, + num_increments: int = 1, + conserve_volume: bool = True, + seed: int | None | Literal["resource_value"] = None, + max_tries: int | None = None, + k_range: ArrayLike = (0, 0.1), + delta_f_range: ArrayLike = (-0.1, 0.1), +) -> dict[str, list[LoadCase] | float]: + """Sample component-wise monotonic exponential deformation-gradient paths. + + Each candidate follows ``F(t) = I + deltaF * (1 - exp(-k*t))``, with + independent components sampled uniformly from ``k_range`` for ``k`` and + ``delta_f_range`` for ``deltaF``. If requested, normalise each gradient to unit + determinant before checking monotonicity of every component. + + The path includes the identity at time zero and ``num_steps`` subsequent + samples. Each subsequent sample defines the target of a load step of duration + ``total_time / num_steps``, with ``num_increments`` increments per step. + + Parameters + ---------- + num_load_cases : int + Number of accepted load cases to generate. Must be a positive integer. + total_time : float | int + Total simulation time for each load case. Must be positive and finite. + num_steps : int + Number of equally spaced load steps per load case. Must be a positive integer. + Each step targets the deformation gradient at its end time. + num_increments : int, optional + Number of simulation increments within each load step. Must be a positive integer. + By default, 1. + conserve_volume : bool, optional + Normalise each sampled deformation gradient to unit determinant before checking + component-wise monotonicity. By default, True. + seed : int | None | Literal["resource_value"], optional + Seed for reproducible random sampling. By default, None (unseeded sampling). + Use "resource_value" to read MATFLOW_RUN_RANDOM_SEED. If the environment + variable is unset, warn and use None. + max_tries : int | None, optional + Maximum total number of candidate paths sampled across all load cases, including + accepted and rejected paths. Must be a positive integer if specified. By default, + None, which uses ``100 * num_load_cases``. Raises RuntimeError if too few load + cases are accepted within this limit. + k_range : array-like, optional + Finite (min, max) sampling bounds for every component of ``k``. Must be + non-negative and ordered. By default, (0, 0.1). + delta_f_range : array-like, optional + Finite, ordered (min, max) sampling bounds for every component of ``deltaF``. By + default, (-0.1, 0.1). Equal bounds fix the sampled value for either range. + Candidates with non-positive or non-finite determinants are rejected. + + Returns + ------- + dict + ``load_case`` contains the accepted load cases. ``monotonicity`` is the fraction + of attempted paths found to be monotonic (accepted cases divided by total + attempts), between 0 and 1. + """ + + if seed == "resource_value": + resource_seed = os.environ.get("MATFLOW_RUN_RANDOM_SEED") + if resource_seed is None: + seed = None + warnings.warn( + "MATFLOW_RUN_RANDOM_SEED is not set; using seed=None for unseeded " + "sampling.", + UserWarning, + stacklevel=2, + ) + else: + seed = int(resource_seed) + rng = np.random.default_rng(seed) + + for name, value in ( + ("num_load_cases", num_load_cases), + ("num_steps", num_steps), + ("num_increments", num_increments), + ): + if isinstance(value, bool) or not isinstance(value, int) or value < 1: + raise ValueError(f"{name} must be a positive integer.") + + if not np.isfinite(total_time) or total_time <= 0: + raise ValueError("total_time must be positive and finite.") + + if max_tries is None: + max_tries = 100 * num_load_cases + + if isinstance(max_tries, bool) or not isinstance(max_tries, int) or max_tries < 1: + raise ValueError("max_tries must be a positive integer.") + + ranges = {} + for name, bounds in (("k_range", k_range), ("delta_f_range", delta_f_range)): + err_msg = f"{name} must contain finite, ordered (min, max) bounds." + try: + bounds = np.asarray(bounds, dtype=float) + except (TypeError, ValueError) as exc: + raise ValueError(err_msg) from exc + if ( + bounds.shape != (2,) + or not np.all(np.isfinite(bounds)) + or bounds[0] > bounds[1] + ): + raise ValueError(err_msg) + ranges[name] = bounds + if ranges["k_range"][0] < 0: + raise ValueError("k_range bounds must be non-negative.") + + times = np.linspace(0, total_time, num_steps + 1) + load_cases: list[LoadCase] = [] + for num_tries in range(1, max_tries + 1): + + k = rng.uniform(*ranges["k_range"], size=(3, 3)) + delta_f = rng.uniform(*ranges["delta_f_range"], size=(3, 3)) + gradients = np.eye(3) + delta_f * ( + 1 - np.exp(-k * times[:, np.newaxis, np.newaxis]) + ) + determinants = np.linalg.det(gradients) + + if not np.all(np.isfinite(determinants) & (determinants > 0)): + continue + + if conserve_volume: + gradients /= np.cbrt(determinants)[:, np.newaxis, np.newaxis] + + # check each component always increasing or always decreasing: + differences = np.diff(gradients, axis=0) + if not np.all( + np.all(differences <= 0, axis=0) | np.all(differences >= 0, axis=0) + ): + continue + + load_cases.append( + LoadCase( + steps=[ + LoadStep( + total_time=total_time / num_steps, + num_increments=num_increments, + target_def_grad=np.ma.masked_array(gradient, mask=False), + stress=np.ma.masked_array(np.zeros((3, 3)), mask=True), + ) + for gradient in gradients[1:] + ] + ) + ) + if len(load_cases) == num_load_cases: + monotonicity = len(load_cases) / num_tries + print( + f"Generated {len(load_cases)} load cases; monotonicity rate was: " + f"{monotonicity}." + ) + return { + "load_case": load_cases, + "monotonicity": monotonicity, + } + + raise RuntimeError( + "Could not sample enough monotonic load cases within " f"{max_tries} attempts." + ) diff --git a/matflow/data/scripts/mtex/plot_orientations.py b/matflow/data/scripts/mtex/plot_orientations.py new file mode 100644 index 00000000..054b2e93 --- /dev/null +++ b/matflow/data/scripts/mtex/plot_orientations.py @@ -0,0 +1,170 @@ +from __future__ import annotations + +from contextlib import contextmanager +import os +from typing import Any, Iterator + +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt +import mtex # type: ignore[import-untyped] +import numpy as np + +from matflow.param_classes.orientations import ( + OrientationRepresentationType, + Orientations, + QuatOrder, +) + + +def _mtex_orientations(orientations: Orientations, crystal_symmetry: str) -> Any: + if orientations.representation.type is not OrientationRepresentationType.QUATERNION: + raise ValueError("MTEX plotting requires quaternion orientations.") + data = np.array(orientations.data, dtype=float, copy=True) + if ( + data.ndim != 2 + or data.shape[1] != 4 + or data.shape[0] == 0 + or not np.all(np.isfinite(data)) + or np.any(np.linalg.norm(data, axis=1) == 0) + ): + raise ValueError( + "orientations.data must contain finite, nonzero (N, 4) quaternions." + ) + if orientations.representation.quat_order is QuatOrder.VECTOR_SCALAR: + data = np.roll(data, 1, axis=1) + elif orientations.representation.quat_order is not QuatOrder.SCALAR_VECTOR: + raise ValueError("orientations must specify a quaternion component order.") + # MatFlow stores the inverse of MTEX's crystal-to-specimen rotation. + data[:, 1:] *= -1 + + alignment = [] + for axis in ("x", "y", "z"): + direction = getattr(orientations.unit_cell_alignment, axis) + if direction is None: + raise ValueError( + "orientations must specify all three unit-cell alignment axes." + ) + name = direction.name.lower().replace("_star", "*") + alignment.append(f"{axis.upper()}||{name}") + cs = mtex.crystalFrame(crystal_symmetry, alignment=alignment) + return mtex.orientation(data, cs, mtex.specimenFrame("1")) + + +def _reference_direction(name: str) -> Any: + if name not in ("x", "y", "z"): + raise ValueError(f"Unknown IPF reference direction: {name!r}; use x, y or z.") + return mtex.vector3d(*np.eye(3)[("x", "y", "z").index(name)]) + + +@contextmanager +def _figure() -> Iterator[Any]: + figure = mtex.newMtexFigure() + try: + yield figure + finally: + plt.close(figure.fig) + + +def _annotate_axes(cs: Any) -> None: + labels = ["x", "y", "z"] + extra = [] + for name in ("a", "b", "c"): + direction = getattr(cs, f"{name}Axis") + unit = direction.normalize().data + matches = np.flatnonzero(np.all(np.isclose(np.eye(3), unit, atol=1e-10), axis=1)) + if matches.size: + labels[int(matches[0])] += f"/{name}" + else: + extra.append((direction, name)) + mtex.annotate(mtex.vector3d(np.eye(3)), label=labels, backgroundColor="white") + for direction, name in extra: + mtex.annotate( + mtex.vector3d(direction.data), label=[name], backgroundColor="white" + ) + + +def plot_orientations( + orientations: Orientations, + crystal_symmetry: str, + use_contours: bool, + pole_figure_directions: list[list[float]] | None = None, + IPF_reference_direction: str = "z", + IPF_reference_directions: list[str] | None = None, + colourbar_limits: list[float] | None = None, + use_one_colourbar: bool = False, +) -> None: + """Write pole figures or inverse pole figures and, for scatter plots, an IPF key. + + Supplying ``pole_figure_directions`` selects pole figures; otherwise plot inverse + pole figures for ``IPF_reference_directions`` (z by default). Inputs are passed + directly by the corresponding schemas; no orientation handoff files are used. + """ + ori = _mtex_orientations(orientations, crystal_symmetry) + reference = _reference_direction(IPF_reference_direction) + if pole_figure_directions is not None: + indices = np.asarray(pole_figure_directions, dtype=float) + if ( + indices.ndim != 2 + or indices.shape[1] != 3 + or indices.shape[0] == 0 + or not np.all(np.isfinite(indices)) + or np.any(np.linalg.norm(indices, axis=1) == 0) + ): + raise ValueError( + "pole_figure_directions must contain finite, nonzero triplets." + ) + directions = [mtex.Miller(*h, ori.CS) for h in indices] + plot = mtex.plotPF + filename = "pole_figure.png" + else: + names = ["z"] if IPF_reference_directions is None else IPF_reference_directions + if not names: + raise ValueError("IPF_reference_directions must not be empty.") + directions = [_reference_direction(name) for name in names] + plot = mtex.plotIPF + filename = "inverse_pole_figure.png" + + if colourbar_limits is not None: + limits = np.asarray(colourbar_limits, dtype=float) + if ( + limits.shape != (2,) + or not np.all(np.isfinite(limits)) + or limits[0] >= limits[1] + ): + raise ValueError( + "colourbar_limits must contain two finite, increasing values." + ) + if num_threads := os.environ.get("MATFLOW_RUN_NUM_THREADS"): + mtex.setMTEXpref("threads", int(num_threads)) + + key = mtex.ipfColorKey(ori.CS, ipfDirection=reference) + with _figure() as figure: + options = { + "antipodal": True, + "how2plot": mtex.plottingConvention(), + } + if use_contours: + color_range = ( + colourbar_limits + if colourbar_limits is not None + else "equal" + if use_one_colourbar + else "individual" + ) + plot(ori, directions, contourf=True, colorRange=color_range, **options) + figure.colorbar(multiple=not use_one_colourbar) + elif pole_figure_directions is not None: + plot(ori, key.orientation2color(ori), directions, **options) + else: + plot(ori, directions, **options) + if pole_figure_directions is not None: + _annotate_axes(ori.CS) + mtex.saveFigure(filename, fig=figure) + + if not use_contours: + with _figure() as figure: + key.plot(how2plot=mtex.plottingConvention()) + mtex.saveFigure("IPF_key.png", fig=figure) diff --git a/matflow/data/scripts/mtex/sample_texture_CRC.m b/matflow/data/scripts/mtex/sample_texture_CRC.m index bb649dcd..b0b870c6 100644 --- a/matflow/data/scripts/mtex/sample_texture_CRC.m +++ b/matflow/data/scripts/mtex/sample_texture_CRC.m @@ -8,34 +8,57 @@ function sample_texture_CRC(inputs_JSON_path, outputs_HDF5_path) referenceFrameTransformation = all_args.EBSD_reference_frame_transformation; specimenSym = all_args.specimen_symmetry; phase = all_args.EBSD_phase; - rotation = all_args.EBSD_rotation; + ebsdRotation = all_args.EBSD_rotation; numOrientations = all_args.num_orientations; + kernelHalfwidthDeg = all_args.kernel_halfwidth_deg; - EBSD_orientations = get_EBSD_orientations_from_CRC_file(CRCFilePath, referenceFrameTransformation, specimenSym, phase, rotation); - ODF = calcDensity(EBSD_orientations); + validateattributes(numOrientations, {'numeric'}, {'real', 'scalar', 'finite', 'integer', 'positive'}, mfilename, 'num_orientations'); + validateattributes(kernelHalfwidthDeg, {'numeric'}, {'real', 'scalar', 'finite', '>', 0, '<=', 180}, mfilename, 'kernel_halfwidth_deg'); + + numThreads = getenv('MATFLOW_RUN_NUM_THREADS'); + if ~isempty(numThreads) + numThreads = str2double(numThreads); + validateattributes(numThreads, {'numeric'}, {'real', 'scalar', 'finite', 'integer', 'positive'}, mfilename, 'MATFLOW_RUN_NUM_THREADS'); + maxNumCompThreads(numThreads); + end + + EBSD_orientations = get_EBSD_orientations_from_CRC_file(CRCFilePath, referenceFrameTransformation, specimenSym, phase, ebsdRotation); + ODF = calcDensity(EBSD_orientations, 'halfwidth', kernelHalfwidthDeg * degree); orientations = discreteSample(ODF, numOrientations); export_orientations_HDF5(orientations, outputs_HDF5_path); end -function EBSD_orientations = get_EBSD_orientations_from_CRC_file(CRCFilePath, referenceFrameTransformation, specimenSym, phase, rotation) +function EBSD_orientations = get_EBSD_orientations_from_CRC_file(CRCFilePath, referenceFrameTransformation, specimenSym, phase, ebsdRotation) if isempty(referenceFrameTransformation) referenceFrameTransformation = {}; else - referenceFrameTransformation = {referenceFrameTransformation}; + switch referenceFrameTransformation + case {'euler_to_spatial', 'convertEuler2SpatialReferenceFrame'} + referenceFrameTransformation = {'convertEuler2SpatialReferenceFrame'}; + case {'spatial_to_euler', 'convertSpatial2EulerReferenceFrame'} + referenceFrameTransformation = {}; + otherwise + error('matflow:InvalidReferenceFrame', 'Unknown EBSD reference frame transformation: %s', referenceFrameTransformation); + end end ebsd = loadEBSD_crc(CRCFilePath, referenceFrameTransformation{:}); - if ~isempty(rotation) + if ~isempty(ebsdRotation) % Note that using rotate appears to remove the non-indexed phase. - rotationEulersRad = num2cell(rotation.euler_angles_deg * degree); - rot = rotation('euler', rotationEulersRad{:}); - - if isfield(rotation, 'keep_XY') + validateattributes(ebsdRotation.euler_angles_deg, {'numeric'}, {'real', 'vector', 'numel', 3, 'finite'}, mfilename, 'EBSD_rotation.euler_angles_deg'); + rotationEulersRad = num2cell(ebsdRotation.euler_angles_deg * degree); + rot = rotation.byEuler(rotationEulersRad{:}); + + keepXY = isfield(ebsdRotation, 'keep_XY') && ebsdRotation.keep_XY; + keepEuler = isfield(ebsdRotation, 'keep_euler') && ebsdRotation.keep_euler; + if keepXY && keepEuler + error('matflow:InvalidRotation', 'EBSD_rotation.keep_XY and keep_euler cannot both be true.'); + elseif keepXY ebsd = rotate(ebsd, rot, 'keepXY'); - elseif isfield(rotation, 'keep_euler') + elseif keepEuler ebsd = rotate(ebsd, rot, 'keepEuler'); else ebsd = rotate(ebsd, rot); @@ -44,6 +67,10 @@ function sample_texture_CRC(inputs_JSON_path, outputs_HDF5_path) end EBSD_orientations = ebsd(phase).orientations; + oriData = [EBSD_orientations.a, EBSD_orientations.b, EBSD_orientations.c, EBSD_orientations.d]; + if isempty(oriData) || any(~isfinite(oriData(:))) + error('matflow:InvalidPhaseOrientations', 'Phase %s has no finite indexed orientations.', phase); + end EBSD_orientations.SS = specimenSymmetry(specimenSym); end diff --git a/matflow/data/scripts/mtex/sample_texture_CTF.m b/matflow/data/scripts/mtex/sample_texture_CTF.m index 88cb3320..334e107c 100644 --- a/matflow/data/scripts/mtex/sample_texture_CTF.m +++ b/matflow/data/scripts/mtex/sample_texture_CTF.m @@ -8,34 +8,57 @@ function sample_texture_CTF(inputs_JSON_path, outputs_HDF5_path) referenceFrameTransformation = all_args.EBSD_reference_frame_transformation; specimenSym = all_args.specimen_symmetry; phase = all_args.EBSD_phase; - rotation = all_args.EBSD_rotation; + ebsdRotation = all_args.EBSD_rotation; numOrientations = all_args.num_orientations; + kernelHalfwidthDeg = all_args.kernel_halfwidth_deg; - EBSD_orientations = get_EBSD_orientations_from_CTF_file(CTFFilePath, referenceFrameTransformation, specimenSym, phase, rotation); - ODF = calcDensity(EBSD_orientations); + validateattributes(numOrientations, {'numeric'}, {'real', 'scalar', 'finite', 'integer', 'positive'}, mfilename, 'num_orientations'); + validateattributes(kernelHalfwidthDeg, {'numeric'}, {'real', 'scalar', 'finite', '>', 0, '<=', 180}, mfilename, 'kernel_halfwidth_deg'); + + numThreads = getenv('MATFLOW_RUN_NUM_THREADS'); + if ~isempty(numThreads) + numThreads = str2double(numThreads); + validateattributes(numThreads, {'numeric'}, {'real', 'scalar', 'finite', 'integer', 'positive'}, mfilename, 'MATFLOW_RUN_NUM_THREADS'); + maxNumCompThreads(numThreads); + end + + EBSD_orientations = get_EBSD_orientations_from_CTF_file(CTFFilePath, referenceFrameTransformation, specimenSym, phase, ebsdRotation); + ODF = calcDensity(EBSD_orientations, 'halfwidth', kernelHalfwidthDeg * degree); orientations = discreteSample(ODF, numOrientations); export_orientations_HDF5(orientations, outputs_HDF5_path); end -function EBSD_orientations = get_EBSD_orientations_from_CTF_file(CTFFilePath, referenceFrameTransformation, specimenSym, phase, rotation) +function EBSD_orientations = get_EBSD_orientations_from_CTF_file(CTFFilePath, referenceFrameTransformation, specimenSym, phase, ebsdRotation) if isempty(referenceFrameTransformation) referenceFrameTransformation = {}; else - referenceFrameTransformation = {referenceFrameTransformation}; + switch referenceFrameTransformation + case {'euler_to_spatial', 'convertEuler2SpatialReferenceFrame'} + referenceFrameTransformation = {'convertEuler2SpatialReferenceFrame'}; + case {'spatial_to_euler', 'convertSpatial2EulerReferenceFrame'} + referenceFrameTransformation = {}; + otherwise + error('matflow:InvalidReferenceFrame', 'Unknown EBSD reference frame transformation: %s', referenceFrameTransformation); + end end ebsd = loadEBSD_ctf(CTFFilePath, referenceFrameTransformation{:}); - if ~isempty(rotation) + if ~isempty(ebsdRotation) % Note that using rotate appears to remove the non-indexed phase. - rotationEulersRad = num2cell(rotation.euler_angles_deg * degree); - rot = rotation('euler', rotationEulersRad{:}); - - if isfield(rotation, 'keep_XY') + validateattributes(ebsdRotation.euler_angles_deg, {'numeric'}, {'real', 'vector', 'numel', 3, 'finite'}, mfilename, 'EBSD_rotation.euler_angles_deg'); + rotationEulersRad = num2cell(ebsdRotation.euler_angles_deg * degree); + rot = rotation.byEuler(rotationEulersRad{:}); + + keepXY = isfield(ebsdRotation, 'keep_XY') && ebsdRotation.keep_XY; + keepEuler = isfield(ebsdRotation, 'keep_euler') && ebsdRotation.keep_euler; + if keepXY && keepEuler + error('matflow:InvalidRotation', 'EBSD_rotation.keep_XY and keep_euler cannot both be true.'); + elseif keepXY ebsd = rotate(ebsd, rot, 'keepXY'); - elseif isfield(rotation, 'keep_euler') + elseif keepEuler ebsd = rotate(ebsd, rot, 'keepEuler'); else ebsd = rotate(ebsd, rot); @@ -44,6 +67,10 @@ function sample_texture_CTF(inputs_JSON_path, outputs_HDF5_path) end EBSD_orientations = ebsd(phase).orientations; + oriData = [EBSD_orientations.a, EBSD_orientations.b, EBSD_orientations.c, EBSD_orientations.d]; + if isempty(oriData) || any(~isfinite(oriData(:))) + error('matflow:InvalidPhaseOrientations', 'Phase %s has no finite indexed orientations.', phase); + end EBSD_orientations.SS = specimenSymmetry(specimenSym); end diff --git a/matflow/data/scripts/mtex/sample_texture_EBSD.py b/matflow/data/scripts/mtex/sample_texture_EBSD.py new file mode 100644 index 00000000..bfdda92c --- /dev/null +++ b/matflow/data/scripts/mtex/sample_texture_EBSD.py @@ -0,0 +1,142 @@ +from __future__ import annotations + +import os +from pathlib import Path +from typing import Literal + +import mtex # type: ignore[import-untyped] +import numpy as np +from numpy.typing import ArrayLike +from typing_extensions import NotRequired, TypedDict + + +class EBSDRotation(TypedDict): + """Rigid map rotation with optional preservation of positions or orientations. + + ``keep_XY`` preserves map positions while rotating orientations. + ``keep_euler`` preserves orientations while rotating map positions. + With neither flag set, both are rotated; the flags cannot both be true. + """ + + euler_angles_deg: ArrayLike + keep_XY: NotRequired[bool] + keep_euler: NotRequired[bool] + + +def sample_texture_EBSD( + specimen_symmetry: str, + EBSD_phase: str, + num_orientations: int, + EBSD_reference_frame_transformation: str | None = "euler_to_spatial", + EBSD_rotation: EBSDRotation | None = None, + CTF_file_path: str | Path | None = None, + CRC_file_path: str | Path | None = None, + kernel_halfwidth_deg: float = 10.0, + num_orientation_sets: int = 1, + sampling_method: Literal["discreteSample", "optimalSample"] = "discreteSample", +) -> dict: + """Estimate one EBSD density and return a list of independently sampled sets. + + The list contains ``num_orientation_sets`` members, even when this is 1. + Sampling advances a single generator seeded from the run's random seed. + ``optimalSample`` optimizes each seeded random sample to better match the density. + """ + if (CTF_file_path is None) == (CRC_file_path is None): + raise ValueError("Specify exactly one of CTF_file_path and CRC_file_path.") + if sampling_method not in ("discreteSample", "optimalSample"): + raise ValueError("sampling_method must be 'discreteSample' or 'optimalSample'.") + if ( + not isinstance(num_orientations, (int, np.integer)) + or isinstance(num_orientations, bool) + or num_orientations <= 0 + ): + raise ValueError("num_orientations must be a positive integer.") + if ( + not isinstance(num_orientation_sets, (int, np.integer)) + or isinstance(num_orientation_sets, bool) + or num_orientation_sets <= 0 + ): + raise ValueError("num_orientation_sets must be a positive integer.") + if ( + isinstance(kernel_halfwidth_deg, bool) + or not np.isfinite(kernel_halfwidth_deg) + or not 0 < kernel_halfwidth_deg <= 180 + ): + raise ValueError("kernel_halfwidth_deg must be finite and in (0, 180] degrees.") + + transform = EBSD_reference_frame_transformation + if transform in ("euler_to_spatial", "convertEuler2SpatialReferenceFrame"): + correction = None + elif transform in ( + None, + "", + "spatial_to_euler", + "convertSpatial2EulerReferenceFrame", + ): + # None would request MTEX's default Oxford correction; identity preserves Euler angles. + correction = mtex.rotation.id + else: + raise ValueError(f"Unknown EBSD reference frame transformation: {transform!r}") + + file_path = CTF_file_path if CTF_file_path is not None else CRC_file_path + ebsd = mtex.loadEBSD( + file_path, + interface="ctf" if CTF_file_path is not None else "crc", + grid=False, + eulerCorrection=correction, + ) + if EBSD_rotation is not None: + angles = np.asarray(EBSD_rotation["euler_angles_deg"], dtype=float) + if angles.shape != (3,) or not np.all(np.isfinite(angles)): + raise ValueError( + "EBSD_rotation.euler_angles_deg must contain three finite angles." + ) + rot = mtex.rotation.byEuler(*(angles * mtex.degree)) + ebsd = ebsd.rotate( + rot, + keepXY=EBSD_rotation.get("keep_XY", False), + keepEuler=EBSD_rotation.get("keep_euler", False), + ) + + source = ebsd[EBSD_phase].orientations + if source.size == 0 or not np.all(np.isfinite(source.data)): + raise ValueError(f"Phase {EBSD_phase!r} has no finite indexed orientations.") + source.SS = mtex.specimenFrame(specimen_symmetry) + + if num_threads := os.environ.get("MATFLOW_RUN_NUM_THREADS"): + # MTEX otherwise counts logical CPUs on Windows, oversubscribing FINUFFT. + mtex.setMTEXpref("threads", int(num_threads)) + + rng = np.random.default_rng(int(os.environ["MATFLOW_RUN_RANDOM_SEED"])) + density = mtex.calcDensity(source, halfwidth=kernel_halfwidth_deg * mtex.degree) + alignment = source.CS.alignment() + if not alignment: + unit_cell_alignment = {"x": "a", "y": "b", "z": "c"} + else: + unit_cell_alignment = { + axis.lower(): direction.replace("*", "_star") + for axis, direction in (item.split("||") for item in alignment) + } + if set(unit_cell_alignment) != {"x", "y", "z"}: + raise ValueError( + f"The crystal frame alignment cannot be represented by MatFlow: {alignment}" + ) + + orientation_sets = [] + for _ in range(num_orientation_sets): + # always use discreteSample first, to ensure the RNG is used (otherwise we'd get + # the same orientations for each set): + sampled = mtex.discreteSample(density, num_orientations, rng=rng) + if sampling_method == "optimalSample": + sampled = mtex.optimalSample(density, sampled) + # MTEX maps crystal to specimen; MatFlow stores the inverse quaternion. + data = np.array(sampled.data, dtype=float, copy=True).reshape(num_orientations, 4) + data[:, 1:] *= -1 + orientation_sets.append( + { + "data": data, + "representation": {"type": "quaternion", "quat_order": "scalar_vector"}, + "unit_cell_alignment": unit_cell_alignment, + } + ) + return {"orientations": orientation_sets} diff --git a/matflow/data/template_components/parameters.yaml b/matflow/data/template_components/parameters.yaml index 42e0d1e5..59850ffc 100644 --- a/matflow/data/template_components/parameters.yaml +++ b/matflow/data/template_components/parameters.yaml @@ -1,3 +1,9 @@ +- type: monotonicity + _validation: + - path: [] + doc: Fraction of attempted load-case paths found to be component-wise monotonic. + condition: { value.is_instance: [float] } + - type: VE_size _validation: - path: [] @@ -188,6 +194,24 @@ and: - value.is_instance: [dict] - key.in: ["euler_angles_deg", "keep_XY", "keep_euler"] +- type: kernel_halfwidth_deg + _validation: + - path: [] + doc: Orientation density kernel halfwidth in degrees; must be finite and in (0, 180]. + condition: { value.is_instance: [int, float] } +- type: num_orientation_sets + _validation: + - path: [] + doc: Positive integer number of independently sampled orientation sets. + condition: { value.is_instance: [int] } +- type: sampling_method + _validation: + - path: [] + doc: MTEX density sampling method; optimalSample optimizes seeded random starting orientations. + condition: + and: + - value.is_instance: [str] + - value.in: [discreteSample, optimalSample] - type: num_orientations _validation: - path: [] diff --git a/matflow/data/template_components/task_schemas.yaml b/matflow/data/template_components/task_schemas.yaml index 1a836d5a..1aa765eb 100644 --- a/matflow/data/template_components/task_schemas.yaml +++ b/matflow/data/template_components/task_schemas.yaml @@ -204,7 +204,10 @@ environment: damask_parse_env - objective: sample_texture - doc: Extract the crystal orientations from a CTF file. + doc: > + Sample an orientation density estimated from a CTF file using MATLAB MTEX. + The kernel halfwidth defaults to 10 degrees; the run's num_threads resource + sets the MATLAB computation thread limit. method: from_CTF_file implementation: mtex inputs: @@ -215,6 +218,8 @@ - parameter: specimen_symmetry - parameter: EBSD_phase - parameter: num_orientations + - parameter: kernel_halfwidth_deg + default_value: 10.0 - parameter: EBSD_reference_frame_transformation default_value: euler_to_spatial - parameter: EBSD_rotation @@ -279,7 +284,11 @@ value.equal_to: precompiled - objective: sample_texture - doc: Extract the crystal orientations from a channel reconstruction file. + doc: > + Sample an orientation density estimated from a CRC file using MATLAB MTEX. + A CPR file with the same stem must be beside the CRC file. + The kernel halfwidth defaults to 10 degrees; the run's num_threads resource + sets the MATLAB computation thread limit. method: from_CRC_file implementation: mtex inputs: @@ -290,6 +299,8 @@ - parameter: specimen_symmetry - parameter: EBSD_phase - parameter: num_orientations + - parameter: kernel_halfwidth_deg + default_value: 10.0 - parameter: EBSD_reference_frame_transformation default_value: euler_to_spatial - parameter: EBSD_rotation @@ -353,6 +364,77 @@ condition: value.equal_to: precompiled +- objective: sample_texture + doc: > + Sample a density estimated from a CTF file using Python MTEX (mtex>=1.0). + Choose discreteSample (default) or optimalSample with sampling_method. + Inputs and outputs are passed directly to and from the Python script. + method: from_CTF_file + implementation: mtex_python + inputs: + - parameter: CTF_file_path + - parameter: specimen_symmetry + - parameter: EBSD_phase + - parameter: num_orientations + - parameter: num_orientation_sets + default_value: 1 + - parameter: sampling_method + default_value: discreteSample + - parameter: kernel_halfwidth_deg + default_value: 10.0 + - parameter: EBSD_reference_frame_transformation + default_value: euler_to_spatial + - parameter: EBSD_rotation + default_value: null + outputs: + - parameter: orientations + multiplicity: <> + actions: + - script: <> + script_data_in: direct + script_data_out: direct + script_exe: python_script + environments: + - scope: + type: any + environment: python_env + +- objective: sample_texture + doc: > + Sample a density estimated from a CRC file using Python MTEX (mtex>=1.0). + Choose discreteSample (default) or optimalSample with sampling_method. + A CPR file with the same stem must be beside the CRC file. + Inputs and outputs are passed directly to and from the Python script. + method: from_CRC_file + implementation: mtex_python + inputs: + - parameter: CRC_file_path + - parameter: specimen_symmetry + - parameter: EBSD_phase + - parameter: num_orientations + - parameter: num_orientation_sets + default_value: 1 + - parameter: sampling_method + default_value: discreteSample + - parameter: kernel_halfwidth_deg + default_value: 10.0 + - parameter: EBSD_reference_frame_transformation + default_value: euler_to_spatial + - parameter: EBSD_rotation + default_value: null + outputs: + - parameter: orientations + multiplicity: <> + actions: + - script: <> + script_data_in: direct + script_data_out: direct + script_exe: python_script + environments: + - scope: + type: any + environment: python_env + - objective: sample_texture doc: Extract the crystal orientations from an ODF material file. method: from_ODF_mat_file @@ -647,6 +729,53 @@ condition: value.equal_to: precompiled +- objective: visualise_orientations + doc: Plot pole figures with Python MTEX using direct orientation inputs. + method: pole_figure + implementation: mtex_python + inputs: + - parameter: orientations + - parameter: crystal_symmetry + - parameter: pole_figure_directions + - parameter: use_contours + - parameter: IPF_reference_direction + default_value: z + - parameter: colourbar_limits + default_value: null + - parameter: use_one_colourbar + default_value: false + save_files: [mtex_pole_figures, mtex_IPF_key] + actions: + - script: <> + script_data_in: direct + script_exe: python_script + environments: + - scope: + type: any + environment: python_env + output_files: [mtex_pole_figures, mtex_IPF_key] + +- objective: visualise_orientations + doc: Plot inverse pole figures with Python MTEX using direct orientation inputs. + method: inverse_pole_figure + implementation: mtex_python + inputs: + - parameter: orientations + - parameter: crystal_symmetry + - parameter: IPF_reference_directions + default_value: ["z"] + - parameter: use_contours + save_files: [mtex_inverse_pole_figures, mtex_IPF_key] + actions: + - script: <> + script_data_in: direct + script_exe: python_script + environments: + - scope: + type: any + environment: python_env + output_files: [mtex_inverse_pole_figures, mtex_IPF_key] + - objective: visualise_multi_orientations # Note: this should be made redundant by some sort of built-in hpcflow functionality to # source from a group (with some rules about how the parameters are merged?). @@ -734,6 +863,93 @@ inputs: - parameter: load_case +- objective: generate_load_case + doc: Generate a uniaxial load case. + web_doc: false + method: uniaxial + inputs: + - parameter: total_time + default_value: 100 + - parameter: num_increments + default_value: 200 + - parameter: direction + default_value: x + - parameter: target_def_grad_rate + default_value: 1.0e-3 + outputs: + - parameter: load_case + actions: + - script: <> + script_data_in: direct + script_data_out: direct + script_exe: python_script + environments: + - scope: + type: any + environment: damask_parse_env + +- objective: generate_load_cases + doc: Generate a list of independent, identical uniaxial load cases. + web_doc: false + method: uniaxial + inputs: + - parameter: num_load_cases + - parameter: total_time + default_value: 100 + - parameter: num_increments + default_value: 200 + - parameter: direction + default_value: x + - parameter: target_def_grad_rate + default_value: 1.0e-3 + outputs: + - parameter: load_case + multiplicity: <> + actions: + - script: <> + script_data_in: direct + script_data_out: direct + script_exe: python_script + environments: + - scope: + type: any + environment: damask_parse_env + +- objective: generate_load_cases + doc: Generate monotonic random exponential deformation-gradient load cases. + method: random_exponential + inputs: + - parameter: num_load_cases + - parameter: total_time + default_value: 100 + - parameter: num_steps + default_value: 100 + - parameter: num_increments + default_value: 1 + - parameter: conserve_volume + default_value: true + - parameter: seed + default_value: resource_value + - parameter: max_tries + default_value: null + - parameter: k_range + default_value: [0, 0.1] + - parameter: delta_f_range + default_value: [-0.1, 0.1] + outputs: + - parameter: load_case + multiplicity: <> + - parameter: monotonicity + actions: + - script: <> + script_data_in: direct + script_data_out: direct + script_exe: python_script + environments: + - scope: + type: any + environment: damask_parse_env + - objective: define_orientations inputs: - parameter: orientations diff --git a/matflow/data/workflows/generate_ROM_training_data.yaml b/matflow/data/workflows/generate_ROM_training_data.yaml new file mode 100644 index 00000000..1cdade3e --- /dev/null +++ b/matflow/data/workflows/generate_ROM_training_data.yaml @@ -0,0 +1,194 @@ +doc: + - > + Generate crystal plasticity training data for training a reduced order model. For + each load case, multiple microstructures (grain morphologies using Dream3D via a + specified grain size distribution, and orientations sampled from EBSD data via a + CTF file using MTEX) are simulated. + + Key outputs: + - volume-averaged Cauchy stress per time increment + - deformation gradient per time increment + + The following variables are available for specification on the command line: + - num_load_cases: number of load cases to consider + - num_microstructures_per_load_case: number of microstructures to simulate for each load case. + + Example submission: + `matflow go demo-workflow generate_ROM_training_data --var num_load_cases 100 --var num_microstructures_per_load_case 5` + + Recommended resources: as specified, use a resources `any` scope for + the DAMASK task, such that all actions (pre-/post-processing and the simulation + itself) run in a single jobscript. Otherwise, a large number of simulation files + may exist on disk concurrently. + + This workflow will typically submit three jobscripts (on scheduled systems): + 1. Non-array job to sample orientations and then generate load cases + 2. Array-job (of maximum size `max_array_size`) to run the Dream3D microstructure + generation + 3. Array-job (of maximum size `max_array_size`) to run the DAMASK task + + Once the workflow has finished, zip the workflow or consolidate the array outputs + for best read performance: + `matflow zip /path/to/workflow` + `matflow workflow /path/to/workflow consolidate-array-outputs` + + In a Python, load array outputs like this: + + ```py + wk = mf.Workflow("/path/to/workflow.zip") + arr_out = wk.get_array_outputs("VE_response.phase_data.vol_avg_def_grad.data") + arr_out.data # numpy array of shape (num_simulations, num_load_increments, 3, 3) + ``` + + From the CLI, export array outputs to Numpy .npy files like this: + + ```bash + matflow workflow /path/to/workflow.zip export-array-output VE_response.phase_data.vol_avg_def_grad.data vol_avg_def_grad.npy + ``` + + TODO: update CTF/CRC file! + TODO: update material_parameters/phase label etc + TODO: validation: check we get similar pole figures with Python and Matlab CTF orientation sampling + +name: gen_ROM_data_Nl_<>_Nm_<> + +inputs: + grid_size: [8, 8, 8] + num_cores_DAMASK: 1 + time_texture: 5 # minutes to sample *all* orientations (may take ~1hr for 25_000) + time_load: 10 # minutes to generate *all* load cases + time_DAMASK: 10 # minutes per single DAMASK simulation + time_Dream3D: 5 # minutes per single Dream3D invocation + material_parameters: <> + max_array_size: 1000 # maximum concurrency + slurm_partition: ukaea-icl + + # specify these in the submission command (e.g. `--var num_load_cases 100`): + num_load_cases: <> + num_microstructures_per_load_case: <> + + # calculated inputs; do not change: + num_simulations: <> * <>)>> + num_stacked_elements: <> // <> + 1)>> + time_DAMASK_all: <> * <>)>> + time_Dream3D_all: <> * <>)>> + +resources: + any: + max_array_size: <> + scheduler_args: + directives: + --partition: <> + +tasks: + - schema: sample_texture_from_CTF_file_mtex_python + resources: + any: + scheduler_args: + directives: + --time: <> + <>)>> + inputs: + CTF_file_path: <> + EBSD_reference_frame_transformation: convertEuler2SpatialReferenceFrame + specimen_symmetry: "1" + EBSD_phase: Al + EBSD_rotation: null + sampling_method: optimalSample + num_orientations: 1000 + num_orientation_sets: <> + + - schema: generate_load_cases_random_exponential + inputs: + total_time: 100 + num_steps: 100 + num_increments: 1 # per load step + conserve_volume: true + resources: + any: + scheduler_args: + directives: + --time: <> + <>)>> + + - schema: generate_volume_element_from_statistics + resources: + any: + scheduler_args: + directives: + --time: <> + save_files: false # disable saving Dream3D XDMF and HDF5 files + inputs: + grid_size: <> + size: [1, 1, 1] + phase_statistics: + - type: primary + name: Al + crystal_structure: cubic + volume_fraction: 1.00 # `volume_fraction` for all phases must sum to 1! + size_distribution: + ESD_mean: 0.5 # 0.12 # Specify either `ESD_mean` or `ESD_log_mean` + ESD_log_stddev: 0.05 + num_bins: 10 # Specify either `num_bins` or `bin_step_size + + - schema: simulate_VE_loading_damask + resources: + any: + num_cores: <> + scheduler_args: + directives: + --time: <> + repeats: + number: <> + nesting_order: 2 + nesting_order: + inputs.load_case: 1 + inputs.volume_element: 2.5 + multi_path_sequences: + - paths: + - inputs.damask_phases.Al.mechanical.plastic.a_sl + - inputs.damask_phases.Al.mechanical.plastic.xi_0_sl.0 + - inputs.damask_phases.Al.mechanical.plastic.xi_inf_sl.0 + - inputs.damask_phases.Al.mechanical.plastic.h_0_sl-sl + nesting_order: 2.5 + values::from_load_txt: + file_path: <> + num_samples: <> + inputs: + homogenization: + SX: + mechanical: { type: "pass" } + N_constituents: 1 + damask_phases: + Al: + lattice: cF + mechanical: + output: [F, P, F_e, F_p, L_p, O] + elastic: + type: Hooke + C_11: 106750000000 + C_12: 60410000000 + C_44: 28340000000 + plastic: + type: phenopowerlaw + N_sl: [12] + a_sl: 2.25 + atol_xi: 1 + dot_gamma_0_sl: 0.001 + h_0_sl-sl: 75.0e+6 + h_sl-sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4] + n_sl: 20 + output: [xi_sl] + xi_0_sl: [31.0e+6] + xi_inf_sl: [63.0e+6] + damask_post_processing: + - name: add_stress_Cauchy + args: { P: P, F: F } + VE_response_data: + phase_data: + - field_name: sigma + phase_name: Al + out_name: vol_avg_stress + transforms: [{ mean_along_axes: 1 }] + - field_name: F + phase_name: Al + out_name: vol_avg_def_grad + transforms: [{ mean_along_axes: 1 }] diff --git a/matflow/param_classes/load.py b/matflow/param_classes/load.py index 9e9d090f..56ceeb4c 100644 --- a/matflow/param_classes/load.py +++ b/matflow/param_classes/load.py @@ -1344,6 +1344,46 @@ def random_3D(cls, **kwargs) -> Self: """ return cls(steps=[LoadStep.random_3D(**kwargs)]) + @classmethod + def random_exponential( + cls, + total_time: int | float, + num_steps: int, + num_increments: int = 1, + conserve_volume: bool = True, + seed: int | None = None, + max_tries: int | None = None, + k_range: ArrayLike = (0, 0.1), + delta_f_range: ArrayLike = (-0.1, 0.1), + ) -> Self: + """A monotonic random exponential deformation-gradient load case. + + Sample ``F(t) = I + deltaF * (1 - exp(-k*t))``, with independent + components sampled uniformly from ``k_range`` (default (0, 0.1)) and + ``delta_f_range`` (default (-0.1, 0.1)). Bounds must be finite and ordered; + ``k_range`` must be non-negative. Equal bounds fix the sampled value. + Optionally conserve volume, then reject paths with non-monotonic + components. Generate ``num_steps`` load steps over ``total_time``, each + with ``num_increments`` increments. Use ``seed`` for reproducible sampling. + ``max_tries`` limits candidate paths and defaults to 100 for this single case. + """ + from matflow.data.scripts.generate_load_cases_random_exponential import ( + generate_load_cases_random_exponential, + ) + + case = generate_load_cases_random_exponential( + num_load_cases=1, + total_time=total_time, + num_steps=num_steps, + num_increments=num_increments, + conserve_volume=conserve_volume, + seed=seed, + max_tries=max_tries, + k_range=k_range, + delta_f_range=delta_f_range, + )["load_case"][0] + return cls(steps=case.steps) + @classmethod def random_polynomial(cls, **kwargs) -> Self: """A multi-step load case in which each deformation gradient component follows diff --git a/matflow/tests/data/define_load.yaml b/matflow/tests/data/define_load.yaml index 532e1bd5..799a2b58 100644 --- a/matflow/tests/data/define_load.yaml +++ b/matflow/tests/data/define_load.yaml @@ -32,3 +32,17 @@ tasks: - [-1.0e-7, 1.0e-7] # c_3 num_increments: 20 seed: 12345 + + # A monotonic exponential deformation-gradient path, normalised to conserve + # volume, with multiple increments per load step and a fixed seed: + - schema: define_load_case + inputs: + load_case::random_exponential: + total_time: 100 + num_steps: 100 + num_increments: 4 + conserve_volume: true + seed: 12345 + max_tries: 200 + k_range: [0.01, 0.08] + delta_f_range: [-0.05, 0.05] diff --git a/matflow/tests/param_classes/test_load.py b/matflow/tests/param_classes/test_load.py index a91b7521..5c510188 100644 --- a/matflow/tests/param_classes/test_load.py +++ b/matflow/tests/param_classes/test_load.py @@ -1,15 +1,226 @@ import sys from pathlib import Path +from unittest.mock import Mock import pytest -import requests import numpy as np import matflow as mf from matflow.param_classes import LoadCase +from matflow.data.scripts.generate_load_case import generate_load_case +from matflow.data.scripts.generate_load_cases import generate_load_cases +from matflow.data.scripts.generate_load_cases_random_exponential import ( + generate_load_cases_random_exponential, +) from matflow.tests.utils import make_test_data_YAML_workflow +def test_generate_load_case(): + inputs = { + "total_time": 100, + "num_increments": 200, + "direction": "x", + "target_def_grad_rate": 1e-3, + } + result = generate_load_case(**inputs) + + assert isinstance(result["load_case"], LoadCase) + assert result["load_case"] == LoadCase.uniaxial(**inputs) + + +@pytest.mark.parametrize("num_load_cases", [1, 3]) +@pytest.mark.parametrize("direction", ["x", "y", "z"]) +def test_generate_load_cases(num_load_cases: int, direction: str) -> None: + result = generate_load_cases( + num_load_cases=num_load_cases, + total_time=10, + num_increments=20, + direction=direction, + target_def_grad_rate=2e-3, + ) + assert set(result) == {"load_case"} + load_cases = result["load_case"] + assert len(load_cases) == num_load_cases + expected = lambda idx: LoadCase.uniaxial( + total_time=10 + idx, + num_increments=20, + direction=direction, + target_def_grad_rate=2e-3, + ) + assert all( + isinstance(case, LoadCase) and case == expected(idx) + for idx, case in enumerate(load_cases) + ) + assert len({id(case) for case in load_cases}) == num_load_cases + assert len({id(case.steps[0]) for case in load_cases}) == num_load_cases + + +@pytest.mark.parametrize("num_load_cases", [0, -1, True]) +def test_generate_load_cases_invalid_count(num_load_cases: int) -> None: + with pytest.raises(ValueError, match="num_load_cases must be a positive integer"): + generate_load_cases( + num_load_cases=num_load_cases, + total_time=100, + num_increments=200, + direction="x", + target_def_grad_rate=1e-3, + ) + + +def test_generate_load_cases_schema() -> None: + template = mf.WorkflowTemplate.from_YAML_string( + "name: generate_load_cases_test\n" + "tasks:\n" + " - schema: generate_load_cases_uniaxial\n" + " inputs:\n" + " num_load_cases: 3\n" + ) + output = template.tasks[0].all_schema_outputs[0] + assert output.parameter.typ == "load_case" + assert output.multiplicity == "<>" + assert output.resolve_multiplicity(3) == 3 + + +@pytest.mark.parametrize("conserve_volume", [True, False]) +def test_generate_load_cases_random_exponential(conserve_volume): + args = { + "num_load_cases": 3, + "total_time": 100, + "num_steps": 100, + "num_increments": 4, + "conserve_volume": conserve_volume, + "seed": 123, + } + result = generate_load_cases_random_exponential(**args) + assert set(result) == {"load_case", "monotonicity"} + assert 0 < result["monotonicity"] <= 1 + if not conserve_volume: + assert result["monotonicity"] == 1 + cases = result["load_case"] + assert len(cases) == args["num_load_cases"] + + for case in cases: + assert isinstance(case, LoadCase) + assert len(case.steps) == args["num_steps"] + gradients = np.array([np.eye(3)] + [s.target_def_grad.data for s in case.steps]) + differences = np.diff(gradients, axis=0) + assert np.all(np.all(differences <= 0, axis=0) | np.all(differences >= 0, axis=0)) + if conserve_volume: + np.testing.assert_allclose(np.linalg.det(gradients), 1) + else: + assert not np.allclose(np.linalg.det(gradients), 1) + assert sum(step.total_time for step in case.steps) == args["total_time"] + for step in case.steps: + assert step.num_increments == args["num_increments"] + assert step.total_time == args["total_time"] / args["num_steps"] + assert not np.any(step.target_def_grad.mask) + assert np.all(step.stress.mask) + + assert cases == generate_load_cases_random_exponential(**args)["load_case"] + assert ( + cases + != generate_load_cases_random_exponential(**{**args, "seed": 124})["load_case"] + ) + + +@pytest.mark.parametrize( + "kwargs, parameter", + [({"num_load_cases": value}, "num_load_cases") for value in [0, -1, True, 1.5]] + + [({"num_steps": value}, "num_steps") for value in [0, -1, True, 1.5]] + + [({"num_increments": value}, "num_increments") for value in [0, -1, True, 1.5]] + + [({"max_tries": value}, "max_tries") for value in [0, -1, True, 1.5]] + + [({"total_time": value}, "total_time") for value in [0, -1, np.inf, np.nan]], +) +def test_generate_load_cases_random_exponential_invalid_inputs(kwargs, parameter): + args = {"num_load_cases": 1, "total_time": 100, "num_steps": 100} + with pytest.raises(ValueError, match=parameter): + generate_load_cases_random_exponential(**{**args, **kwargs}) + + +def test_generate_load_cases_random_exponential_monotonicity(monkeypatch): + rng = Mock(wraps=np.random.default_rng(123)) + monkeypatch.setattr(np.random, "default_rng", lambda seed: rng) + result = generate_load_cases_random_exponential( + num_load_cases=3, total_time=100, num_steps=100, seed=123 + ) + # Each attempted path samples k and deltaF, including rejected paths. + attempts = rng.uniform.call_count // 2 + assert attempts > 3 + assert result["monotonicity"] == 3 / attempts + + +@pytest.mark.parametrize("parameter", ["k_range", "delta_f_range"]) +@pytest.mark.parametrize("bounds", [[], [0], [[0, 1]], [1, 0], [0, np.inf], [np.nan, 1]]) +def test_generate_load_cases_random_exponential_invalid_range(parameter, bounds): + with pytest.raises(ValueError, match=parameter): + generate_load_cases_random_exponential( + num_load_cases=1, total_time=100, num_steps=10, **{parameter: bounds} + ) + + +def test_generate_load_cases_random_exponential_sampling_ranges(): + args = {"num_load_cases": 1, "total_time": 100, "num_steps": 10, "seed": 123} + with pytest.raises(ValueError, match="k_range"): + generate_load_cases_random_exponential(**args, k_range=(-0.1, 0.1)) + # Zero decay or zero amplitude must produce an undeformed path. + for ranges in ( + {"k_range": [0, 0], "delta_f_range": [0.2, 0.3]}, + {"k_range": [0.02, 0.03], "delta_f_range": [0, 0]}, + ): + result = generate_load_cases_random_exponential(**args, **ranges) + assert result["monotonicity"] == 1 + for step in result["load_case"][0].steps: + np.testing.assert_array_equal(step.target_def_grad.data, np.eye(3)) + + +def test_generate_load_cases_random_exponential_invalid_determinants(): + # Large negative amplitudes make these candidates physically invalid. + with pytest.raises(RuntimeError, match="within 2 attempts"): + generate_load_cases_random_exponential( + num_load_cases=1, + total_time=100, + num_steps=10, + k_range=[0.1, 0.1], + delta_f_range=[-2, -2], + max_tries=2, + ) + + +def test_generate_load_cases_random_exponential_max_tries(): + # With volume conservation disabled, every candidate is monotonic. + args = { + "num_load_cases": 3, + "total_time": 100, + "num_steps": 10, + "conserve_volume": False, + "seed": 123, + } + with pytest.raises(RuntimeError, match="within 2 attempts"): + generate_load_cases_random_exponential(**args, max_tries=2) + cases = generate_load_cases_random_exponential(**args, max_tries=3)["load_case"] + assert len(cases) == 3 + assert all(step.num_increments == 1 for case in cases for step in case.steps) + + +def test_generate_load_case_random_exponential_schema(): + template = mf.WorkflowTemplate.from_YAML_string( + "name: generate_load_cases_test\n" + "tasks:\n" + " - schema: generate_load_cases_random_exponential\n" + " inputs:\n" + " num_load_cases: 3\n" + ) + task = template.tasks[0] + output = task.all_schema_outputs[0] + assert output.parameter.typ == "load_case" + assert output.multiplicity == "<>" + assert output.resolve_multiplicity(3) == 3 + + monotonicity_output = task.all_schema_outputs[1] + assert monotonicity_output.parameter.typ == "monotonicity" + assert monotonicity_output.multiplicity is None + + def test_single_multistep_uniaxial(): args = { "num_increments": 200, @@ -37,6 +248,30 @@ def test_load_case_yaml_init(tmp_path: Path, load_case_1: LoadCase): seed=12345, ) + # third task uses the same factory syntax for a random exponential path: + load_case_exp = wk.tasks.define_load_case_3.elements[0].inputs.load_case.value + args = { + "total_time": 100, + "num_steps": 100, + "num_increments": 4, + "conserve_volume": True, + "seed": 12345, + "max_tries": 200, + "k_range": [0.01, 0.08], + "delta_f_range": [-0.05, 0.05], + } + assert load_case_exp == LoadCase.random_exponential(**args) + assert ( + load_case_exp + == generate_load_cases_random_exponential(num_load_cases=1, **args)["load_case"][ + 0 + ] + ) + assert load_case_exp.num_steps == args["num_steps"] + assert all( + step.num_increments == args["num_increments"] for step in load_case_exp.steps + ) + def test_load_case_from_npz_file(): npz_file_path = mf.get_demo_data_file_path("load_cases.npz") diff --git a/matflow/tests/test_plot_orientations_mtex_python.py b/matflow/tests/test_plot_orientations_mtex_python.py new file mode 100644 index 00000000..6d3fc1da --- /dev/null +++ b/matflow/tests/test_plot_orientations_mtex_python.py @@ -0,0 +1,248 @@ +from pathlib import Path +from unittest.mock import patch + +import numpy as np +import pytest + +import matflow as mf +from matflow.param_classes.orientations import ( + LatticeDirection, + OrientationRepresentation, + Orientations, + QuatOrder, + UnitCellAlignment, +) + +mtex = pytest.importorskip("mtex", reason="Python MTEX is an optional task dependency") + +from matflow.data.scripts.mtex.plot_orientations import ( + _mtex_orientations, + plot_orientations, +) +from matflow.tests.test_sample_texture_mtex_python import ebsd_file + +import matplotlib.pyplot as plt +from PIL import Image + + +@pytest.fixture +def orientations() -> Orientations: + data = np.random.default_rng(42).normal(size=(24, 4)) + data /= np.linalg.norm(data, axis=1, keepdims=True) + return Orientations( + data=data, + representation=OrientationRepresentation.quaternion(), + unit_cell_alignment=UnitCellAlignment( + x=LatticeDirection.A, y=LatticeDirection.B, z=LatticeDirection.C + ), + ) + + +@pytest.mark.parametrize("order", list(QuatOrder)) +@pytest.mark.parametrize("symmetry", ["cubic", "hexagonal"]) +def test_orientation_conversion( + orientations: Orientations, order: QuatOrder, symmetry: str +) -> None: + expected = orientations.data.copy() + expected[:, 1:] *= -1 + if order is QuatOrder.VECTOR_SCALAR: + orientations.data = np.roll(orientations.data, -1, axis=1) + orientations.representation = OrientationRepresentation.quaternion(order) + if symmetry == "hexagonal": + orientations.unit_cell_alignment = UnitCellAlignment( + x=LatticeDirection.A_STAR, y=LatticeDirection.B, z=LatticeDirection.C + ) + original = orientations.data.copy() + result = _mtex_orientations(orientations, symmetry) + np.testing.assert_allclose(result.data, expected) + np.testing.assert_array_equal(orientations.data, original) + if symmetry == "hexagonal": + assert result.CS.alignment() == ["X||a*", "Y||b", "Z||c"] + else: + assert result.CS.alignment() == [] + assert result.SS.sym.name == mtex.specimenFrame("1").sym.name + + +@pytest.mark.parametrize("data", [[], [[0, 0, 0, 0]], [[np.nan, 0, 0, 1]], [[1, 0, 0]]]) +def test_invalid_quaternions(orientations: Orientations, data: list) -> None: + orientations.data = np.asarray(data) + with pytest.raises(ValueError, match="nonzero.*quaternions"): + _mtex_orientations(orientations, "cubic") + + +def test_unsupported_representation(orientations: Orientations) -> None: + orientations.representation = OrientationRepresentation.euler() + with pytest.raises(ValueError, match="requires quaternion"): + _mtex_orientations(orientations, "cubic") + + +def test_missing_alignment(orientations: Orientations) -> None: + orientations.unit_cell_alignment.x = None + with pytest.raises(ValueError, match="all three"): + _mtex_orientations(orientations, "cubic") + + +@pytest.mark.parametrize( + "kwargs, message", + [ + ({"pole_figure_directions": []}, "pole_figure_directions"), + ({"pole_figure_directions": [[0, 0, 0]]}, "pole_figure_directions"), + ({"pole_figure_directions": [[1, 2]]}, "pole_figure_directions"), + ({"pole_figure_directions": [[0, np.inf, 1]]}, "pole_figure_directions"), + ({"IPF_reference_direction": "bad"}, "reference direction"), + ({"IPF_reference_directions": []}, "must not be empty"), + ({"IPF_reference_directions": ["bad"]}, "reference direction"), + ({"colourbar_limits": [1.0, 0.0]}, "colourbar_limits"), + ({"colourbar_limits": [0.0, np.nan]}, "colourbar_limits"), + ], +) +def test_invalid_plot_options( + orientations: Orientations, kwargs: dict, message: str +) -> None: + with pytest.raises(ValueError, match=message): + plot_orientations(orientations, "cubic", False, **kwargs) + + +def assert_png(path: Path) -> None: + with Image.open(path) as image: + assert image.format == "PNG" + assert image.width > 100 and image.height > 100 + pixels = np.asarray(image.convert("RGB")) + assert np.std(pixels) > 5 + + +@pytest.mark.parametrize("inverse", [False, True]) +@pytest.mark.parametrize("contours", [False, True]) +@pytest.mark.parametrize("symmetry", ["cubic", "hexagonal"]) +def test_render_figures( + orientations: Orientations, + tmp_path: Path, + monkeypatch: pytest.MonkeyPatch, + inverse: bool, + contours: bool, + symmetry: str, +) -> None: + monkeypatch.chdir(tmp_path) + monkeypatch.setenv("MATFLOW_RUN_NUM_THREADS", "1") + if symmetry == "hexagonal": + orientations.unit_cell_alignment.x = LatticeDirection.A_STAR + before = plt.get_fignums() + with patch.object(mtex, "setMTEXpref", wraps=mtex.setMTEXpref) as threads: + plot_orientations( + orientations, + symmetry, + contours, + pole_figure_directions=None if inverse else [[0, 0, 1], [1, 1, 1]], + IPF_reference_directions=["x", "z"], + ) + threads.assert_called_once_with("threads", 1) + assert plt.get_fignums() == before + assert_png(tmp_path / ("inverse_pole_figure.png" if inverse else "pole_figure.png")) + if contours: + assert not (tmp_path / "IPF_key.png").exists() + else: + assert_png(tmp_path / "IPF_key.png") + + +@pytest.mark.parametrize("one_bar", [False, True]) +@pytest.mark.parametrize("limits", [None, [0.0, 5.0]]) +def test_colourbar_options( + orientations: Orientations, + tmp_path: Path, + monkeypatch: pytest.MonkeyPatch, + one_bar: bool, + limits: list[float] | None, +) -> None: + monkeypatch.chdir(tmp_path) + figures = [] + save = mtex.saveFigure + + def capture(filename, fig): + figures.append(fig) + save(filename, fig=fig) + + with patch.object(mtex, "saveFigure", side_effect=capture): + plot_orientations( + orientations, + "cubic", + True, + pole_figure_directions=[[0, 0, 1], [1, 1, 1]], + colourbar_limits=limits, + use_one_colourbar=one_bar, + ) + assert len(figures) == 1 + assert len(figures[0].colorbars) == (1 if one_bar else 2) + ranges = [ax._mtexSphericalPlot.mappable.get_clim() for ax in figures[0].children] + if limits is not None: + assert all(value == tuple(limits) for value in ranges) + elif one_bar: + assert ranges[0] == ranges[1] + + +def test_figure_closed_after_error( + orientations: Orientations, tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + monkeypatch.chdir(tmp_path) + before = plt.get_fignums() + with patch.object(mtex, "saveFigure", side_effect=OSError("cannot write figure")): + with pytest.raises(OSError, match="cannot write"): + plot_orientations(orientations, "cubic", False) + assert plt.get_fignums() == before + + +@pytest.mark.parametrize("method", ["pole_figure", "inverse_pole_figure"]) +def test_plot_schema(method: str) -> None: + schema = mf.task_schemas.get(f"visualise_orientations_{method}_mtex_python") + inputs = {item.parameter.typ for item in schema.inputs} + assert "execution_mode" not in inputs + action = schema.actions[0] + assert set(action.script_data_in) == {f"inputs.{name}" for name in inputs} + assert all(value["format"] == "direct" for value in action.script_data_in.values()) + assert action.script_exe == "python_script" + + +@pytest.mark.integration +@pytest.mark.parametrize("contours", [False, True]) +def test_sample_and_plot_workflow( + ebsd_file: Path, tmp_path: Path, contours: bool +) -> None: + file_type = ebsd_file.suffix[1:].upper() + template = mf.WorkflowTemplate.from_YAML_string( + "name: sample_and_plot_python\n" + "resources:\n" + " any:\n" + " random_seed: 1234\n" + "tasks:\n" + f" - schema: sample_texture_from_{file_type}_file_mtex_python\n" + " inputs:\n" + f" {file_type}_file_path: '{ebsd_file}'\n" + " specimen_symmetry: '1'\n" + " EBSD_phase: Al\n" + " num_orientations: 24\n" + " num_orientation_sets: 2\n" + " - schema: visualise_orientations_pole_figure_mtex_python\n" + " inputs:\n" + " crystal_symmetry: cubic\n" + " pole_figure_directions: [[0, 0, 1], [1, 1, 1]]\n" + f" use_contours: {str(contours).lower()}\n" + " - schema: visualise_orientations_inverse_pole_figure_mtex_python\n" + " inputs:\n" + " crystal_symmetry: cubic\n" + " IPF_reference_directions: [x, z]\n" + f" use_contours: {str(contours).lower()}\n" + ) + workflow = mf.Workflow.from_template(template, path=tmp_path) + assert workflow.tasks[1].num_elements == workflow.tasks[2].num_elements == 2 + workflow.submit(wait=True, add_to_known=False) + for run in workflow.get_all_EARs(): + assert run.status.name == "success" + reopened = mf.Workflow(workflow.path) + for run in reopened.get_all_EARs()[1:]: + files = run.output_files + if run.action.output_files[0].label == "mtex_pole_figures": + path = files.mtex_pole_figures.value + else: + path = files.mtex_inverse_pole_figures.value + assert_png(Path(path)) + if not contours: + assert_png(Path(files.mtex_IPF_key.value)) diff --git a/matflow/tests/test_sample_texture_mtex.py b/matflow/tests/test_sample_texture_mtex.py new file mode 100644 index 00000000..34275fc4 --- /dev/null +++ b/matflow/tests/test_sample_texture_mtex.py @@ -0,0 +1,49 @@ +from pathlib import Path + +import pytest + +import matflow as mf + + +@pytest.mark.parametrize("file_type", ["CTF", "CRC"]) +@pytest.mark.parametrize("halfwidth", [None, 25.0]) +@pytest.mark.parametrize("execution_mode", ["script", "compile", "precompiled"]) +def test_sample_texture_matlab_schema( + file_type: str, halfwidth: float | None, execution_mode: str +) -> None: + schema_name = f"sample_texture_from_{file_type}_file_mtex" + template = mf.WorkflowTemplate.from_YAML_string( + "name: sample_texture_matlab\n" + "tasks:\n" + f" - schema: {schema_name}\n" + " inputs:\n" + f" execution_mode: {execution_mode}\n" + f" {file_type}_file_path: texture.{file_type.lower()}\n" + " specimen_symmetry: '1'\n" + " EBSD_phase: Al\n" + " num_orientations: 12\n" + + ("" if halfwidth is None else f" kernel_halfwidth_deg: {halfwidth}\n") + ) + assert template.tasks[0].all_schema_outputs[0].parameter.typ == "orientations" + schema = mf.task_schemas.get(schema_name) + assert schema.outputs[0].multiplicity is None + inputs = {item.parameter.typ: item for item in schema.inputs} + defaults = { + name: item.default_value.value + for name, item in inputs.items() + if item.default_value is not None + } + assert defaults["kernel_halfwidth_deg"] == 10.0 + assert "num_orientation_sets" not in inputs + assert defaults["execution_mode"] == "script" + assert defaults["EBSD_reference_frame_transformation"] == "euler_to_spatial" + assert defaults["EBSD_rotation"] is None + + expected_inputs = {f"inputs.{name}": {"format": "json"} for name in inputs} + expected_outputs = {"outputs.orientations": {"format": "hdf5"}} + assert len(schema.actions) == 4 + for action in schema.actions[1:3]: + assert action.script_data_in == expected_inputs + assert action.script_data_out == expected_outputs + assert schema.actions[3].program_data_in == expected_inputs + assert schema.actions[3].program_data_out == expected_outputs diff --git a/matflow/tests/test_sample_texture_mtex_python.py b/matflow/tests/test_sample_texture_mtex_python.py new file mode 100644 index 00000000..e68b2aaf --- /dev/null +++ b/matflow/tests/test_sample_texture_mtex_python.py @@ -0,0 +1,374 @@ +from pathlib import Path +from typing import Literal +from unittest.mock import patch +import warnings + +import numpy as np +import pytest + +import matflow as mf +from matflow.param_classes.orientations import ( + LatticeDirection, + Orientations, + QuatOrder, + UnitCellAlignment, +) + +mtex = pytest.importorskip("mtex", reason="Python MTEX is an optional task dependency") + +from matflow.data.scripts.mtex.sample_texture_EBSD import ( + EBSDRotation, + sample_texture_EBSD, +) + + +@pytest.fixture(params=["ctf", "crc"]) +def ebsd_file(tmp_path: Path, request: pytest.FixtureRequest) -> Path: + angles = np.array([[10, 20, 30], [40, 50, 60], [0, 0, 0], [70, 80, 90]]) + phases = [1, 2, 0, 1] + if request.param == "ctf": + path = tmp_path / "texture.ctf" + header = ( + "Channel Text File\n" + "XCells\t4\nYCells\t1\nXStep\t1\nYStep\t1\n" + "Phases\t2\n" + "4;4;4\t90;90;90\tAl\t11\t225\n" + "3;3;5\t90;90;120\tMg\t9\t194\n" + "Phase\tX\tY\tEuler1\tEuler2\tEuler3\tBC\n" + ) + rows = [ + "\t".join(str(v) for v in [phase, i, 0, *euler, 100]) + for i, (phase, euler) in enumerate(zip(phases, angles)) + ] + path.write_text(header + "\n".join(rows) + "\n", encoding="utf-8", newline="\n") + else: + path = tmp_path / "texture.crc" + header = ( + "[Job]\nXCells=4\nYCells=1\nGridDistX=1\nGridDistY=1\n" + "[Fields]\nCount=3\nField1=3\nField2=4\nField3=5\n" + "[Phases]\nCount=2\n" + "[Phase1]\nStructureName=Al\nLaueGroup=11\n" + "a=4\nb=4\nc=4\nAlpha=90\nBeta=90\nGamma=90\n" + "[Phase2]\nStructureName=Mg\nLaueGroup=9\n" + "a=3\nb=3\nc=5\nAlpha=90\nBeta=90\nGamma=120\n" + ) + path.with_suffix(".cpr").write_text(header, encoding="utf-8", newline="\n") + data = np.zeros(4, dtype=[("phase", "i1"), ("angles", " dict: + path_key = "CTF_file_path" if path.suffix == ".ctf" else "CRC_file_path" + return sample_texture_EBSD( + **{ + path_key: path, + "specimen_symmetry": "cubic", + "EBSD_phase": "Al", + "num_orientations": 12, + **kwargs, + } + ) + + +@pytest.mark.parametrize("phase", ["Al", "Mg"]) +def test_sample_texture_real_readers( + ebsd_file: Path, monkeypatch: pytest.MonkeyPatch, phase: str +) -> None: + monkeypatch.setenv("MATFLOW_RUN_RANDOM_SEED", "1234") + first = sample_file(ebsd_file, EBSD_phase=phase)["orientations"][0] + second = sample_file(ebsd_file, EBSD_phase=phase)["orientations"][0] + ori = Orientations(**first) + assert ori.data.shape == (12, 4) + assert np.all(np.isfinite(ori.data)) + np.testing.assert_allclose(np.linalg.norm(ori.data, axis=1), 1) + np.testing.assert_array_equal(first["data"], second["data"]) + assert ori.representation.quat_order is QuatOrder.SCALAR_VECTOR + alignment = ( + UnitCellAlignment( + x=LatticeDirection.A, y=LatticeDirection.B, z=LatticeDirection.C + ) + if phase == "Al" + else UnitCellAlignment( + x=LatticeDirection.A_STAR, y=LatticeDirection.B, z=LatticeDirection.C + ) + ) + assert ori.unit_cell_alignment == alignment + monkeypatch.setenv("MATFLOW_RUN_RANDOM_SEED", "5678") + third = sample_file(ebsd_file, EBSD_phase=phase)["orientations"][0] + assert not np.allclose(first["data"], third["data"]) + + +@pytest.mark.parametrize( + "transform", + [ + "euler_to_spatial", + "convertEuler2SpatialReferenceFrame", + "spatial_to_euler", + "convertSpatial2EulerReferenceFrame", + None, + ], +) +@pytest.mark.parametrize("keep", [None, "keep_XY", "keep_euler"]) +def test_sample_texture_frames_and_rotation( + ebsd_file: Path, + monkeypatch: pytest.MonkeyPatch, + transform: str | None, + keep: Literal["keep_XY", "keep_euler"] | None, +) -> None: + monkeypatch.setenv("MATFLOW_RUN_RANDOM_SEED", "1234") + rotation: EBSDRotation = {"euler_angles_deg": [15, 25, 35]} + if keep is not None: + rotation[keep] = True + with patch.object( + mtex, "calcDensity", side_effect=lambda source, halfwidth: source + ) as density: + with patch.object( + mtex, "discreteSample", side_effect=lambda source, n, rng: source + ): + result = sample_file( + ebsd_file, + num_orientations=2, + EBSD_reference_frame_transformation=transform, + EBSD_rotation=rotation, + ) + + expected = mtex.rotation.byEuler( + np.array([10, 70]) * mtex.degree, + np.array([20, 80]) * mtex.degree, + np.array([30, 90]) * mtex.degree, + ) + if transform in ("euler_to_spatial", "convertEuler2SpatialReferenceFrame"): + expected = mtex.rotation.byAxisAngle(mtex.vector3d(0, 1, 0), np.pi) * expected + if keep != "keep_euler": + expected = ( + mtex.rotation.byEuler(*(np.array([15, 25, 35]) * mtex.degree)) * expected + ) + source = density.call_args.args[0] + np.testing.assert_allclose(source.data, expected.data, atol=1e-7) + assert source.SS.sym.name == mtex.specimenFrame("cubic").sym.name + np.testing.assert_allclose( + result["orientations"][0]["data"], expected.inv().data, atol=1e-7 + ) + + +@pytest.mark.parametrize( + "kwargs,match", + [ + ({"num_orientations": 0}, "positive integer"), + ({"num_orientations": -1}, "positive integer"), + ({"sampling_method": "bad"}, "sampling_method"), + ({"EBSD_reference_frame_transformation": "bad"}, "reference frame"), + ({"EBSD_rotation": {"euler_angles_deg": [1, 2]}}, "three finite angles"), + ({"EBSD_rotation": {"euler_angles_deg": [1, 2, np.nan]}}, "three finite angles"), + ( + { + "EBSD_rotation": { + "euler_angles_deg": [0, 0, 0], + "keep_XY": True, + "keep_euler": True, + } + }, + "not both", + ), + ({"EBSD_phase": "missing"}, "phase"), + ], +) +def test_sample_texture_invalid_inputs(ebsd_file: Path, kwargs: dict, match: str) -> None: + with pytest.raises(ValueError, match=match): + sample_file(ebsd_file, **kwargs) + + +@pytest.mark.parametrize("halfwidth", [0, -1, 181, np.nan, np.inf, -np.inf, True]) +def test_sample_texture_invalid_halfwidth(ebsd_file: Path, halfwidth: float) -> None: + with pytest.raises(ValueError, match="kernel_halfwidth_deg"): + sample_file(ebsd_file, kernel_halfwidth_deg=halfwidth) + + +@pytest.mark.parametrize("count", [0, -1, True, 1.5, None]) +def test_sample_texture_invalid_set_count(ebsd_file: Path, count: object) -> None: + with pytest.raises( + ValueError, match="num_orientation_sets must be a positive integer" + ): + sample_file(ebsd_file, num_orientation_sets=count) + + +@pytest.mark.parametrize("count", [1, 3]) +def test_sample_texture_sets( + ebsd_file: Path, monkeypatch: pytest.MonkeyPatch, count: int +) -> None: + monkeypatch.setenv("MATFLOW_RUN_RANDOM_SEED", "1234") + with patch.object(mtex, "calcDensity", wraps=mtex.calcDensity) as density: + first = sample_file(ebsd_file, num_orientation_sets=count)["orientations"] + density.assert_called_once() + second = sample_file(ebsd_file, num_orientation_sets=count)["orientations"] + assert len(first) == count + for member, repeated in zip(first, second): + assert Orientations(**member).data.shape == (12, 4) + np.testing.assert_array_equal(member["data"], repeated["data"]) + for member in first[1:]: + assert not np.array_equal(first[0]["data"], member["data"]) + single = sample_file(ebsd_file)["orientations"][0] + np.testing.assert_array_equal(first[0]["data"], single["data"]) + + +def test_sample_texture_optimal_sets( + ebsd_file: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + monkeypatch.setenv("MATFLOW_RUN_RANDOM_SEED", "1234") + monkeypatch.setenv("MATFLOW_RUN_NUM_THREADS", "1") + kwargs = {"sampling_method": "optimalSample", "num_orientation_sets": 2} + first = sample_file(ebsd_file, **kwargs)["orientations"] + second = sample_file(ebsd_file, **kwargs)["orientations"] + assert len(first) == 2 + for member, repeated in zip(first, second): + data = Orientations(**member).data + assert data.shape == (12, 4) + assert np.all(np.isfinite(data)) + np.testing.assert_allclose(np.linalg.norm(data, axis=1), 1) + np.testing.assert_array_equal(member["data"], repeated["data"]) + assert not np.allclose(first[0]["data"], first[1]["data"]) + monkeypatch.setenv("MATFLOW_RUN_RANDOM_SEED", "5678") + third = sample_file(ebsd_file, sampling_method="optimalSample")["orientations"][0] + assert not np.allclose(first[0]["data"], third["data"]) + + +@pytest.mark.parametrize("halfwidth", [None, 5.0, 25.0]) +def test_sample_texture_halfwidth( + ebsd_file: Path, monkeypatch: pytest.MonkeyPatch, halfwidth: float | None +) -> None: + monkeypatch.setenv("MATFLOW_RUN_RANDOM_SEED", "1234") + kwargs = {} if halfwidth is None else {"kernel_halfwidth_deg": halfwidth} + with patch.object(mtex, "calcDensity", wraps=mtex.calcDensity) as density: + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + sample_file(ebsd_file, specimen_symmetry="1", **kwargs) + assert not any("halfwidth" in str(warning.message) for warning in caught) + assert density.call_args.kwargs["halfwidth"] == pytest.approx( + (10.0 if halfwidth is None else halfwidth) * np.pi / 180 + ) + + +def test_sample_texture_requires_one_file() -> None: + with pytest.raises(ValueError, match="exactly one"): + sample_texture_EBSD("cubic", "Al", 12) + with pytest.raises(ValueError, match="exactly one"): + sample_texture_EBSD( + "cubic", "Al", 12, CTF_file_path="a.ctf", CRC_file_path="a.crc" + ) + + +def test_crc_requires_cpr(ebsd_file: Path) -> None: + if ebsd_file.suffix != ".crc": + pytest.skip("CRC-specific error") + ebsd_file.with_suffix(".cpr").unlink() + with pytest.raises(ValueError, match="no .cpr beside"): + sample_file(ebsd_file) + + +@pytest.mark.parametrize("num_threads", [None, "3"]) +def test_sample_texture_threads( + ebsd_file: Path, monkeypatch: pytest.MonkeyPatch, num_threads: str | None +) -> None: + monkeypatch.setenv("MATFLOW_RUN_RANDOM_SEED", "1234") + if num_threads is None: + monkeypatch.delenv("MATFLOW_RUN_NUM_THREADS", raising=False) + else: + monkeypatch.setenv("MATFLOW_RUN_NUM_THREADS", num_threads) + with patch.object(mtex, "setMTEXpref") as set_pref: + sample_file(ebsd_file, specimen_symmetry="1") + if num_threads is None: + set_pref.assert_not_called() + else: + set_pref.assert_called_once_with("threads", 3) + + +@pytest.mark.parametrize("file_type", ["CTF", "CRC"]) +@pytest.mark.parametrize("sampling_method", ["discreteSample", "optimalSample"]) +def test_sample_texture_schema(file_type: str, sampling_method: str) -> None: + template = mf.WorkflowTemplate.from_YAML_string( + "name: sample_texture_python\n" + "tasks:\n" + f" - schema: sample_texture_from_{file_type}_file_mtex_python\n" + " inputs:\n" + f" {file_type}_file_path: texture.{file_type.lower()}\n" + " specimen_symmetry: cubic\n" + " EBSD_phase: Al\n" + " num_orientations: 12\n" + f" sampling_method: {sampling_method}\n" + ) + assert template.tasks[0].all_schema_outputs[0].parameter.typ == "orientations" + schema = mf.task_schemas.get(f"sample_texture_from_{file_type}_file_mtex_python") + assert schema.outputs[0].multiplicity == "<>" + count_input = next( + item for item in schema.inputs if item.parameter.typ == "num_orientation_sets" + ) + assert count_input.default_value.value == 1 + method_input = next( + item for item in schema.inputs if item.parameter.typ == "sampling_method" + ) + assert method_input.default_value.value == "discreteSample" + action = schema.actions[0] + assert set(action.script_data_in) == { + f"inputs.{file_type}_file_path", + "inputs.specimen_symmetry", + "inputs.EBSD_phase", + "inputs.num_orientations", + "inputs.num_orientation_sets", + "inputs.sampling_method", + "inputs.kernel_halfwidth_deg", + "inputs.EBSD_reference_frame_transformation", + "inputs.EBSD_rotation", + } + assert all(value["format"] == "direct" for value in action.script_data_in.values()) + assert action.script_data_out == {"outputs.orientations": {"format": "direct"}} + + +@pytest.mark.integration +@pytest.mark.parametrize("halfwidth", [None, 25.0]) +@pytest.mark.parametrize("count", [1, 3]) +def test_sample_texture_workflow( + ebsd_file: Path, tmp_path: Path, halfwidth: float | None, count: int +) -> None: + file_type = ebsd_file.suffix[1:].upper() + template = mf.WorkflowTemplate.from_YAML_string( + "name: sample_texture_python\n" + "resources:\n" + " any:\n" + " random_seed: 1234\n" + "tasks:\n" + f" - schema: sample_texture_from_{file_type}_file_mtex_python\n" + " inputs:\n" + f" {file_type}_file_path: '{ebsd_file}'\n" + " specimen_symmetry: '1'\n" + " EBSD_phase: Al\n" + " num_orientations: 12\n" + + ("" if count == 1 else f" num_orientation_sets: {count}\n") + + ("" if halfwidth is None else f" kernel_halfwidth_deg: {halfwidth}\n") + ) + template.tasks.append( + mf.Task( + mf.TaskSchema( + objective="consume_orientation_set", + inputs=[mf.SchemaInput("orientations")], + ) + ) + ) + workflow = mf.Workflow.from_template(template, path=tmp_path) + assert workflow.tasks[1].num_elements == count + workflow.submit(wait=True, add_to_known=False) + output = workflow.tasks[0].elements[0].outputs.orientations.value + assert isinstance(output, list) + assert len(output) == count + for data in output: + member = Orientations(**data) + assert member.data.shape == (12, 4) + np.testing.assert_allclose(np.linalg.norm(member.data, axis=1), 1) + reopened = mf.Workflow(workflow.path) + persisted = reopened.tasks[0].elements[0].outputs.orientations.value + for member, data in zip(persisted, output): + assert Orientations(**member) == Orientations(**data) + for idx, element in enumerate(reopened.tasks[1].elements): + assert element.inputs.orientations.value == Orientations(**output[idx]) diff --git a/poetry.lock b/poetry.lock index 261a37ee..34c5a73a 100644 --- a/poetry.lock +++ b/poetry.lock @@ -1546,14 +1546,14 @@ numpy = ">=1.21.2" [[package]] name = "hpcflow" -version = "0.4.0a6" +version = "0.4.0a9" description = "Computational workflow management" optional = false python-versions = "<3.14,>=3.10" groups = ["main"] files = [ - {file = "hpcflow-0.4.0a6-py3-none-any.whl", hash = "sha256:d2f6e085d5dca5a2c31ef86fa57158059e1fbb5cbe4f53e623247d1caa1c27a6"}, - {file = "hpcflow-0.4.0a6.tar.gz", hash = "sha256:3dd6a58c3b2bc0ccd452db43cab8ba4aacf2ddbe745418d230680d8a7841fd54"}, + {file = "hpcflow-0.4.0a9-py3-none-any.whl", hash = "sha256:4605f2d633f00a19bda63c9f1e6d064b245de241d8146dabbe0ffcc93902156c"}, + {file = "hpcflow-0.4.0a9.tar.gz", hash = "sha256:300bbf9d0cf0f526cc2914fc6b53d5bc964a6437d4d1bc61f0534176d102109a"}, ] [package.dependencies] @@ -5706,4 +5706,4 @@ test = ["pytest"] [metadata] lock-version = "2.1" python-versions = ">=3.10,<3.14" -content-hash = "5ffccfb0ab6108623da5043a0dbd3f224401f407ee12b1ef372d21d137774b0f" +content-hash = "91dbb84d07b5574fe244b7a4104e2a80659170638f9cbb3632e05c3c163ad095" diff --git a/pyproject.toml b/pyproject.toml index 33c827b9..a95e2049 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,7 +20,7 @@ pytest = {version = "^9.0.2", optional = true} matplotlib = "^3.7" typing-extensions = "^4.12.2" skops = "^0.13.0" -hpcflow = "0.4.0a6" +hpcflow = "0.4.0a9" [tool.poetry.group.dev.dependencies] pylint = "^2.12.2"