API Reference

Mineral

Store polycrystal texture data for a single mineral phase.

Fields

• phase: MineralPhase
• fabric: MineralFabric
• regime: DeformationRegime
• n_grains: number of grains
• fractions: list of volume fraction snapshots (Vector{Vector{Float64}})
• orientations: list of orientation snapshots (Vector{Array{Float64,3}})
• seed: RNG seed for initial random orientations

Supported mineral phases.

Supported mineral fabrics / CRSS distributions.

Deformation mechanism regimes.

Default parameters for DRex simulations.

Return default parameters as a Dict{Symbol,Any} for easy mutation.

Stiffness tensors (Voigt 6×6), units of GPa.

update_orientations!(mineral, params, deformation_gradient,
                     get_velocity_gradient, pathline; get_regime=nothing)

Update crystalline orientations and grain volume distribution. Returns the updated deformation gradient.

• params — Dict with parameter keys matching default_params()
• deformation_gradient — 3×3 matrix
• get_velocity_gradient(t, x) — callable returning 3×3 velocity gradient
• pathline — (tstart, tend, getposition) where getposition(t) returns 3D position

update_all!(minerals, params, deformation_gradient,
            get_velocity_gradient, pathline; get_regime=nothing)

Update orientations and volume distributions for all mineral phases. Returns the updated deformation gradient tensor.

run_pathlines_batch!(minerals_per_tracer, params, pathlines_data; backend=CPU())

Integrate CPO for n_tracers pathlines using Tsit5 (adaptive Runge-Kutta).

Each outer timestep is integrated with update_all! (Tsit5), which internally uses the same GPU kernel as update_orientations!. Velocity gradients are linearly interpolated between consecutive pathline samples within each interval.

Arguments

• minerals_per_tracer — Vector{Vector{Mineral{T}}} of length n_tracers; each inner vector holds one Mineral per phase in the same order for every tracer.
• params — parameter dict from default_params(), identical for all tracers.
• pathlines_data — Vector of length n_tracers; each element is a named tuple or tuple (timestamps, positions, velocity_gradients):
  • timestamps — Vector{Float64} of length n_steps
  • positions — Vector{Vector} of 3-D positions at each timestamp
  • velocity_gradients — Vector{Matrix} of 3×3 VG matrices at each timestamp
• backend — KernelAbstractions backend (default CPU()).
• snapshot_stride — save an orientation/fraction snapshot every this many steps (default 1 = every step). Use e.g. 10 to reduce snapshot count 10×, keeping only ~51 output frames for 500 inner steps. The returned strains vector is subsampled to match.

Returns Vector{Vector{T}} of per-tracer accumulated strain values (length 1 + cld(n_steps-1, snapshot_stride)).

derivatives!(orientations_diff, fractions_diff,
             regime, phase, fabric, n_grains,
             orientations, fractions,
             strain_rate, velocity_gradient, deformation_gradient_spin,
             stress_exponent, deformation_exponent,
             nucleation_efficiency, gbm_mobility, volume_fraction)

Compute derivatives of orientation and volume distribution (in-place, minimally allocating).

• orientations_diff — pre-allocated n_grains×3×3 output array for rotation rates
• fractions_diff — pre-allocated n_grains output vector for volume fraction changes
• orientations — n_grains×3×3 input orientation matrices
• fractions — n_grains input volume fractions

The inner grain-level computations are fully allocation-free using StaticArrays.

get_crss(phase, fabric) -> SVector{4,Float64}

Get normalised Critical Resolved Shear Stress for the given mineral phase and fabric. Returns a static vector (allocation-free).

voigt_averages(minerals, phase_assemblage, phase_fractions; elastic_tensors=StiffnessTensors())

Calculate elastic tensors as Voigt averages of a collection of minerals. Returns a 6×6 Voigt matrix.

voigt_decompose(matrix)

Decompose elastic tensor (6×6 Voigt matrix) into dilatational and deviatoric stiffness tensors. Returns (dij, vij) where dij = Cijkk and vij = Cikjk.

voigt_to_elastic_tensor(matrix)

Create 4th-order elastic tensor from a 6×6 Voigt matrix.

elastic_tensor_to_voigt(tensor)

Create a 6×6 Voigt matrix from a 4th-order elastic tensor.

voigt_matrix_to_vector(matrix)

Create the 21-component Voigt vector from the 6×6 Voigt matrix.

voigt_vector_to_matrix(vector)

Create the 6×6 Voigt matrix from the 21-component Voigt vector.

rotate_tensor(tensor, rotation)

Rotate a 4th-order tensor using a 3×3 rotation matrix.

polar_decompose(M; left=true)

Compute polar decomposition. If left=true, returns (R, V) with M = VR. If left=false, returns (R, U) with M = RU.

upper_tri_to_symmetric(A)

Create symmetric matrix using the upper triangle of the input matrix.

misorientation_index(orientations, system; bins=nothing)

Calculate M-index for polycrystal orientations.

bingham_average(orientations; axis="a")

Compute Bingham average (antipodally symmetric mean) of orientation matrices. orientations is N×3×3. Returns the mean direction as a 3-vector.

finite_strain(deformation_gradient)

Extract largest principal strain value and maximum extension axis from the 3×3 deformation gradient. Returns (strain, axis).

symmetry_pgr(orientations; axis="a")

Compute Point, Girdle, Random symmetry diagnostics. Returns (P, G, R).

elasticity_components(voigt_matrices)

Calculate elasticity decompositions for N×6×6 array of Voigt matrices. Returns a Dict with keys: bulkmodulus, shearmodulus, percentanisotropy, percenthexagonal, percenttetragonal, percentorthorhombic, percentmonoclinic, percenttriclinic, hexagonal_axis.

coaxial_index(orientations; axis1="b", axis2="a")

Calculate coaxial "BA" index. Returns value in [0, 1].

smallest_angle(vector, axis; plane=nothing)

Get smallest angle (degrees) between a unit vector and a bidirectional axis.

simple_shear_2d(direction, deformation_plane, strain_rate)

Return (velocity, velocity_gradient) callables for 2D simple shear.

• direction — "X", "Y", or "Z": velocity vector direction
• deformation_plane — "X", "Y", or "Z": direction of velocity gradient
• strain_rate — 1/2 × magnitude of the largest eigenvalue of the velocity gradient

cell_2d(horizontal, vertical, velocity_edge; edge_length=2.0)

Return (velocity, velocity_gradient) callables for a steady-state 2D Stokes convection cell.

The cell is centered at (0,0) with velocity: u = U cos(πx/d) sin(πz/d) ĥ - U sin(πx/d) cos(πz/d) v̂

• horizontal — "X", "Y", or "Z"
• vertical — "X", "Y", or "Z"
• velocity_edge — velocity magnitude at cell edge center
• edge_length — cell edge length (default 2.0)

corner_2d(horizontal, vertical, plate_speed)

Return (velocity, velocity_gradient) callables for steady-state 2D corner flow.

Velocity field: u = (2U/π)(arctan(x/(-z)) + xz/(x²+z²)) x̂ + (2U/π)(z²/(x²+z²)) ẑ

• horizontal — "X", "Y", or "Z"
• vertical — "X", "Y", or "Z"
• plate_speed — upper boundary ("plate") speed

get_pathline(final_location, get_velocity, get_velocity_gradient,
             min_coords, max_coords; max_strain=10.0, regular_steps=nothing)

Determine the pathline for a particle in a steady-state flow by integrating backwards in time from final_location.

Returns (timestamps, position_fn) where position_fn(t) gives the 3D position at time t (negative times, since integration is backwards).

• final_location — 3D coordinates of the final (exit) location
• get_velocity(t, x) — velocity callable
• get_velocity_gradient(t, x) — velocity gradient callable
• min_coords — lower-bound coordinates of the domain
• max_coords — upper-bound coordinates of the domain
• max_strain — target strain at the final location (default 10)
• regular_steps — if set, return regularly-spaced timestamps (default: solver timestamps)

poles(orientations, ref_axes="xz", hkl=[1,0,0])

Extract 3D vectors of crystallographic directions from orientation matrices. orientations is an N×3×3 array.

lambert_equal_area(xvals, yvals, zvals)

Project axial data from a 3D sphere onto a 2D disk using Lambert equal-area projection.

to_cartesian(ϕ, θ; r=1.0)

Convert spherical to Cartesian coordinates. ϕ = longitude (azimuth), θ = colatitude (inclination).

to_spherical(x, y, z)

Convert Cartesian to spherical coordinates. Returns (r, ϕ, θ).

resample_orientations(orientations, fractions; n_samples=nothing, seed=nothing)

Return new samples from orientations weighted by the volume distribution.

Accepts either:

• orientations::Vector{Array{Float64,3}}, fractions::Vector{Vector{Float64}} (list of N×3×3 arrays, list of N-vectors — matches Python API)
• orientations::AbstractArray{Float64,4}, fractions::AbstractArray{Float64,2} (P×M×3×3 and P×M arrays)

misorientation_angles(q1_array, q2_array)

Calculate minimum misorientation angles for collections of rotation quaternions. q1array has shape (N, A, 4), q2array has shape (N, B, 4). Returns N-element vector of minimum misorientation angles in degrees.

misorientation_hist(orientations, system; bins=nothing)

Calculate misorientation histogram for polycrystal orientations. Returns (counts, bin_edges) as from a normalised histogram.

Uses only proper rotation symmetry operations (not reflections) to match the theoretical random distribution of Skemer et al. (2005). Since proper rotations form a group, the minimum over cross-combinations reduces to iterating only nsym per pair (fix one side, vary the other).

The inner loop is allocation-free: symmetry-applied quaternions are precomputed per grain (n × nsym tuples) and the histogram is built on-the-fly.

misorientations_random(low, high, system)

Expected count of random misorientation angles in (low, high) degrees.

strain_increment(dt, velocity_gradient)

Calculate strain increment for a given time increment and velocity gradient. Returns tensorial strain increment ε = γ/2.

apply_gbs!(orientations, fractions, gbs_threshold, orientations_prev, n_grains)

Apply grain boundary sliding for small grains. Modifies orientations and fractions in-place.

LaMEMSnapshot

One timestep of LaMEM output: velocity gradient tensor and (optionally) velocity on a structured rectilinear grid.

CPOTracer

A Lagrangian particle that carries crystallographic preferred orientation (CPO) through a LaMEM simulation.

compute_cpo_from_lamem(sim_name, sim_dir;
                       target_positions=nothing, initial_positions=nothing,
                       output_dir="$(sim_name)_tracers",
                       skip_initial_steps=5, start_step=nothing, end_step=nothing,
                       steady_state_step=nothing, steady_state_duration=nothing,
                       steady_state_n_steps=100,
                       drex_params=default_params(), n_grains=1000, seed=42,
                       n_substeps=5, fabric=olivine_A, vel_grad_field=:vel_gr_tensor)

All-in-one function: load LaMEM output, (optionally backtrack positions), evolve CPO forward, and write a Paraview time-series (.pvd + .vtp files).

Returns (tracers, snapshots_used).

Exactly one of target_positions or initial_positions must be provided:

• target_positions — [x,y,z] km where you want CPO at the last used snapshot; positions are backtracked to find the material source.
• initial_positions — [x,y,z] km of material at the first used snapshot; CPO is evolved forward from here directly.

Time-dependent mode (default)

Reads velocity gradients from every snapshot in the usable window.

• skip_initial_steps — drop the first N snapshots (spin-up transients)
• start_step / end_step — optional sub-window within the usable range

Steady-state mode

Set steady_state_step and steady_state_duration to use a single snapshot's velocity field for the entire integration (constant flow assumption).

• steady_state_step — index into the usable snapshot window for the reference field
• steady_state_duration — total integration time (Myr)
• steady_state_n_steps — number of output steps (default 100)

Requires LaMEM and GeophysicalModelGenerator to be loaded.

load_snapshots(sim_name, dir=""; vel_grad_field=:vel_gr_tensor, load_velocity=true)
→ Vector{LaMEMSnapshot}

Read all timesteps of a LaMEM simulation. Requires LaMEM and GeophysicalModelGenerator to be loaded.

create_tracers(positions; n_grains=1000, seed=42,
               phase_assemblage=[olivine,enstatite], fabric=olivine_A)
→ Vector{CPOTracer}

Create CPO tracers with random initial orientations at the given positions. Requires LaMEM and GeophysicalModelGenerator to be loaded.

backtrack_positions(target_positions, snapshots; n_substeps=5)
→ Vector{Vector{Float64}}

Advect positions backward through the snapshot sequence to find source positions at snapshots[1]. Requires LaMEM and GeophysicalModelGenerator to be loaded.

evolve_cpo!(tracers, params, snapshots; advect=true, n_substeps=5)

Evolve CPO on all tracers through the LaMEM snapshot sequence. Requires LaMEM and GeophysicalModelGenerator to be loaded.

run_cpo_at_locations(target_positions, snapshots, drex_params;
                     n_grains=1000, seed=42, n_substeps=5, fabric=olivine_A)
→ Vector{CPOTracer}

High-level function: compute CPO at known target locations at the final snapshot. Backtracks positions to find material source, then runs CPO forward. Requires LaMEM and GeophysicalModelGenerator to be loaded.

trilinear_interpolate(x_grid, y_grid, z_grid, field_3d, point) → Float64

Trilinear interpolation of a scalar 3D field at a point. Boundary-clamped. Requires LaMEM and GeophysicalModelGenerator to be loaded.

interpolate_vel_grad(snap, point) → Matrix{Float64}(3,3)

Interpolate the 3×3 velocity gradient at an arbitrary point from a snapshot. Requires LaMEM and GeophysicalModelGenerator to be loaded.

interpolate_velocity(snap, point) → Vector{Float64}(3)

Interpolate the velocity vector at an arbitrary point from a snapshot. Requires LaMEM and GeophysicalModelGenerator to be loaded.

make_velocity_gradient_func(snap_prev, snap_next) → Function

Return a (t, x) → 3×3 Matrix callable that linearly interpolates the velocity gradient in time and trilinearly in space between two snapshots. Requires LaMEM and GeophysicalModelGenerator to be loaded.

make_velocity_func(snap_prev, snap_next) → Function

Return a (t, x) → Vector callable that linearly interpolates velocity in time and trilinearly in space between two snapshots. Requires LaMEM and GeophysicalModelGenerator to be loaded.

read_scsv(filepath) -> NamedTuple

Read an SCSV file and return a NamedTuple whose field names match the column names. Each field is a Vector of the appropriate type.

save_scsv(filepath, schema, data; comments=nothing)

Write data columns to an SCSV file.

• schema — Dict with keys "delimiter", "missing", "fields"
• data — Vector of column vectors, one per field

write_scsv_header(io, schema; comments=nothing)

Write the YAML header of an SCSV file to io.