ELAS3D-Xtal: An OpenMP-accelerated crystal elasticity solver with automated experiment-driven microstructure generation

This paper introduces ELAS3D-Xtal, a high-performance Fortran/OpenMP upgrade of the NIST ELAS3D voxel-based finite element solver for computing 3D elastic fields in polycrystals with defects. The code supports crystal anisotropy by precomputing rotated stiffness tensors from user-specified orientations and solves the equilibrium problem with a matrix-free, OpenMP-parallel preconditioned conjugate-gradient (PCG) method using a point-block Jacobi preconditioner. On a single shared-memory multicore PC, OpenMP threading accelerates the baseline CG solver by ~10X, while the block-preconditioned CG solver achieves 53-61X speedup relative to the serial CG baseline for meshes from 100^3 to 500^3 voxels (scaling to domains up to 800^3 voxels). Accuracy is validated against the analytical Eshelby inclusion solution. ELAS3D-Xtal also integrates microstructure construction, including statistically calibrated polycrystal generation via spatial filtering and parallel voxel-to-grain assignment, direct pore insertion from XCT centroid/radius data, and texture assignment. Full-field phase, orientation, and stress outputs are written in HDF5 to enable scalable post-processing and defect-mechanics workflows. Applications are demonstrated for (i) anisotropy-controlled defect-scale stress fields and (ii) LPBF SS316L microstructures with gas, lack-of-fusion, and keyhole pore morphologies.

💡 Research Summary

ELAS3D‑Xtal is a high‑performance, OpenMP‑accelerated extension of the NIST ELAS3D voxel‑based finite‑element code, designed to compute three‑dimensional elastic fields in polycrystalline materials that contain defects such as pores or lack‑of‑fusion voids. The new solver adds full crystal anisotropy by pre‑computing rotated stiffness tensors for each grain orientation supplied by the user. These tensors are generated from Rodrigues vectors, converted to 3 × 3 rotation matrices, and then assembled into 6 × 6 transformation matrices that rotate the local (crystal) stiffness into the global (sample) frame.

The equilibrium displacement field is obtained by minimizing the total potential energy, which leads to a large sparse linear system Au + b = 0. Because assembling the global stiffness matrix A would be prohibitive for problems with >10⁸ degrees of freedom, ELAS3D‑Xtal employs a matrix‑free preconditioned conjugate‑gradient (PCG) algorithm. The matrix‑vector product is evaluated with a 27‑point stencil, and the implementation uses cache‑blocked loops to improve data locality on modern multi‑core CPUs.

A key innovation is the point‑block Jacobi preconditioner. Instead of the scalar Jacobi scaling used in the legacy code, each node’s three displacement components are treated as a dense 3 × 3 block. The diagonal block is assembled from the contributions of the eight surrounding elements, inverted analytically, and stored as M⁻¹. This block‑diagonal preconditioner dramatically improves convergence for highly anisotropic crystals where the off‑diagonal stiffness terms couple the displacement components.

OpenMP parallelism is applied throughout the preconditioner construction and the PCG iterations. By processing nodes in contiguous chunks (block_size) that fit 50‑70 % of the per‑core L1/L2 cache, the code achieves excellent scaling on shared‑memory machines. Benchmarks on a single workstation show that the baseline CG solver is accelerated roughly 10× by OpenMP threading, while the block‑preconditioned CG achieves a 53‑61× speed‑up over the serial CG baseline for meshes ranging from 100³ to 500³ voxels (10⁶–1.25×10⁸ cells). The solver also scales to domains of 800³ voxels (≈5×10⁸ cells).

Accuracy is validated against the analytical Eshelby inclusion solution, with errors below 1 % for stress and displacement fields.

ELAS3D‑Xtal integrates an automated microstructure generation pipeline. Grain structures are created using spatial filtering to match prescribed grain‑size distributions and aspect‑ratio statistics, followed by a parallel ellipsoidal Voronoi tessellation that assigns voxels to grains. Crystallographic texture can be imposed either randomly or from user‑provided orientation distribution functions. Pore geometries are directly inserted from X‑ray computed tomography (XCT) centroid and radius data; pores are modeled by assigning zero stiffness, effectively creating voids. All microstructure data, as well as full‑field phase, orientation, and stress outputs, are written in HDF5 format, enabling scalable post‑processing and easy interfacing with visualization or data‑analysis tools.

Two application examples illustrate the code’s capabilities. First, a parametric study shows how varying crystal anisotropy controls the stress field around a spherical inclusion, highlighting the importance of correctly representing grain orientations in defect mechanics. Second, realistic laser powder‑bed‑fusion (LPBF) SS316L microstructures are simulated: statistically calibrated grain structures with measured texture are combined with XCT‑derived gas pores, lack‑of‑fusion defects, and keyhole pores. The results demonstrate pronounced stress concentrations that depend on both defect morphology and surrounding grain anisotropy, providing insight into fatigue‑life prediction for additively manufactured components.

In summary, ELAS3D‑Xtal delivers a complete, shared‑memory parallel framework for high‑resolution elastic analysis of polycrystals with defects. By coupling advanced preconditioning, cache‑aware OpenMP parallelism, and an integrated microstructure generator, it enables researchers and engineers to perform large‑scale, physics‑rich simulations on a single workstation, bridging the gap between experimental microstructure characterization and predictive mechanical modeling. Future work will extend the solver to nonlinear elasticity, crystal plasticity, multiphysics coupling (thermal, electrical), and GPU acceleration.