Advanced Features

This page catalogues the cross-cutting numerical and algorithmic techniques implemented in OpenPhase — techniques that are used by more than one module and are worth knowing about when you build a custom simulation. Each section names the feature, summarises what it does, cites the module(s) where the implementation lives, and points to the detailed module page.

Active parameter tracking

A straightforward multiphase-field evaluation scales as $\mathcal{O}(n^3)$ in the number of phase fields $n$ — one pairwise contribution per pair per point. Since the double-obstacle potential produces $\sin$-shaped profiles with finitely supported interfaces, every grid point actually touches only a small subset of phase fields. OpenPhase exploits this by maintaining an indicator function ${\chi }_{\alpha }$ that masks out points where ${\varphi }_{\alpha }=0\wedge {\mathrm{\nabla }}^2{\varphi }_{\alpha }=0$:

$$\begin{array}{}\text{(1)}& {\chi }_{\alpha }=\{\begin{array}{ll}1& \text{if }{\varphi }_{\alpha }\ne 0\text{or}{\mathrm{\nabla }}^2{\varphi }_{\alpha }\ne 0,\\ 0& \text{otherwise.}\end{array}\end{array}$$

The phase-field update loop skips the masked-out combinations, so the effective cost is proportional to the number of active phase fields per point (typically two or three at interfaces, zero in the bulk). The full derivation is on the Phase Field page; the code lives in the PhaseField and DoubleObstacle modules.

Interface regularization

Strong bulk driving forces distort the steady phase-field profile, which in turn distorts interface kinetics. OpenPhase addresses this with a regularisation field $\kappa$ computed by the InterfaceRegularization module and fed into the optional overload DoubleObstacle::CalculatePhaseFieldIncrements(Phi, IP, dG, Kappa). The regulariser smooths transient profile deviations while letting the interface move at the velocity dictated by its driving force.

Anti-trapping current

Diffuse-interface diffusion models exhibit a spurious "solute trapping" at moderate interface velocities. OpenPhase ships the Karma / Eiken–Steinbach anti-trapping current as a switchable flux inside the binary-alloy diffusion solver:

$$\begin{array}{}\text{(2)}& \frac{\partial c}{\partial t}=\mathrm{\nabla }\cdot (D(\varphi )\mathrm{\nabla }c)-\mathrm{\nabla }\cdot {\mathbf{j}}_{\text{AT}}.\end{array}$$

Set $EnableAntiTrapping: Yes in the @EquilibriumPartitionDiffusionBinary block; the flux ${\mathbf{j}}_{\text{AT}}$ is computed internally from ${\varphi }_{\alpha }$ and the per-phase diffusivities.

Multiphase-junction partitioning

At junctions where three or more phases meet, standard binary- partitioning relations no longer close. OpenPhase's Thermodynamics page derives the generalised partitioning scheme used internally: an equal-undercooling condition across every pair simultaneously, with a per-pair partitioning coefficient

$$\begin{array}{}\text{(3)}& K_{\alpha \beta }=m_{\alpha \beta }/m_{\beta \alpha }.\end{array}$$

The scheme is used under the hood by the binary-diffusion solver and by the Grand Potential Solver.

Spectral solvers

Two mechanical / field solvers use FFT-based spectral methods:

• ElasticitySolverSpectral solves mechanical equilibrium with a homogeneous / inhomogeneous split that sidesteps the Fourier-periodic-only limitation. Use the accompanying MechanicalLoads for time-scheduled loading.
• ElectricSolverSpectral solves the Poisson equation for the electric potential in frequency space with cost $\mathcal{O}(N\log N)$ where $N=N_x{N}_y{N}_z$.

Both depend on FFTW3 (MKL FFTW on Windows).

Lattice Boltzmann for fluid flow

The FluidDynamics module provides a D3Q27 Lattice Boltzmann solver with:

• Single-phase incompressible flow (BGK-LBM).
• Multi-phase flow via pseudo-potential methods (Benzi, Kupershtokh).
• Phase-field coupling for diffuse-interface solid-fluid boundaries.
• Solid-fluid interaction (bounce-back, drag, rigid-body motion).
• Solid-solid short-range repulsion (Standard / Wang).
• Thermally compressible flow (low-Mach-number extension).
• Reactive flow coupling to species transport and heat sources.

Cantera-backed reactive chemistry

If the library was built with CANTERA_PATH set at CMake time, the ReactiveFlows module loads a Cantera mechanism file ($ReactionMechanism, typically a .yaml) and couples the computed species and energy source terms into SpeciesTransport and EnergyTransport.

Moving reference frame

Directional-solidification-style setups use MovingFrame to shift the entire simulation domain (every Storage3D field) along the chosen axis when a selected trigger phase crosses a chosen grid position. The shift keeps the interesting region inside the computational box without growing the box itself.

Faceted interface anisotropy

Beyond the standard cubic and hexagonal anisotropy models, OpenPhase implements a faceted family for materials with strong faceting tendencies:

$$\begin{array}{}\text{(4)}& \sigma (\theta )={\sigma }_F\sqrt{{\sin }^2\theta +{ϵ}_{EF}^2{\cos }^2\theta }\end{array}$$

where $\theta$ is the inclination angle to the nearest facet family, ${\sigma }_F$ the facet energy, and ${ϵ}_{EF}$ the facet anisotropy parameter. The facet-selection logic and the enum values FACETED, FACETEDFULL, HEXBOETTGER, HEXSUN, HEXYANG, CUBIC, CUBICFULL, ISO live on InterfaceProperties.

Temperature control

Temperature supports an arbitrary list of control modes ($ControlMode_<n>) that ramp or relax the temperature between holding plateaus:

• LINEAR — constant-rate cooling / heating up to a target.
• NEWTON — Newtonian relaxation towards a target.
• NONE — no evolution.

Each record has its own $ActivationTime_<n> so chained experiments (isothermal hold, quench, re-heat) are expressed in the input file.

Latent-heat coupling

Phase-change energy is released / absorbed via an antisymmetric latent-heat matrix on Temperature:

$$\begin{array}{}\text{(5)}& L_{\alpha \beta }=-L_{\beta \alpha }.\end{array}$$

$LatentHeatMode selects OFF, LOCAL, or GLOBAL coupling; per-phase $VolumetricHeatCapacity_<n> and the matrix $LatentHeat_<α>_<β> complete the specification. LOCAL deposits the latent heat at the grid point where the phase transition occurs; GLOBAL distributes it over the domain.

Driving-force limiting

The shared DrivingForce accumulator caps each per-pair contribution to a physically stable maximum. $Limit_<a>_<b> (alias $CutOff_<a>_<b>) sets the cap as a fraction of the driving force that would destabilise the interface profile; the default of 0.95 is a reasonable starting point. With $Average: Yes the accumulated driving force is smoothed across the diffuse interface before consumption; with $bUnify: Yes multi- physics contributions to the same pair are unified into a single value.

Robust rotation extraction

Crystal-plasticity and mechanical coupling in phase-field simulations require a clean separation of the rotational and strain parts of a deformation. OpenPhase implements an iterative quaternion-based extraction that is numerically stable in the range of deformation gradients that occur in practice:

/*
 * A Robust Method to Extract the Rotational Part of Deformations
 * Müller et al., SCA / MIG 2016:
 * https://animation.rwth-aachen.de/media/papers/2016-MIG-StableRotation.pdf
 *
 * The input quaternion can be:
 *   1. The previous iterate of a time-dependent solve.
 *   2. Quat(Mat) / |Quat(Mat)| on first use.
 *   3. (1, 0, 0, 0) as a last resort.
 *
 * Convergence: maxIter of 20 suffices for every case we have tested.
 */
void Tools::ExtractRotation(dMatrix3x3 &Mat, Quaternion &Quat,
                            const size_t maxIter)
{
    for (size_t iter = 0; iter < maxIter; iter++)
    {
        dMatrix3x3 Rot = Quat.getRotationMatrix();

        std::vector<dVector3> Rot_col = Col(Rot);
        std::vector<dVector3> Mat_col = Col(Mat);

        double b = 1.0 / fabs(Rot_col[0]*Mat_col[0] +
                              Rot_col[1]*Mat_col[1] +
                              Rot_col[2]*Mat_col[2]) + 1.0e-9;

        dVector3 omega = (Rot_col[0].cross(Mat_col[0]) +
                          Rot_col[1].cross(Mat_col[1]) +
                          Rot_col[2].cross(Mat_col[2])) * b;

        double Angle = omega.length();
        if (Angle < 1.0e-9) break;

        dVector3 Axis = omega * (1.0 / Angle);
        Quaternion tempQuat;
        tempQuat.set(Axis, Angle);
        Quat = tempQuat * Quat;
        Quat.normalize();
    }
}

The algorithm:

1. Start with an initial quaternion.
2. Convert it to a rotation matrix.
3. Compute the axis-angle correction between the current rotation and the target matrix via cross products of column vectors.
4. Apply the correction to the quaternion; normalise.
5. Stop when the angular correction falls below ${10}^{-9}$ or after maxIter iterations.

The same Tools::ExtractRotation is used by crystal-plasticity slip- system rotations and by the crystallographic misorientation writers on Orientations.

Strain measures at large deformation

The default St. Venant–Kirchhoff hyperelastic model is built on the Green-Lagrange strain, which has clean physical meaning but is unstable in compression. At large compressive loads, the Mechanics module offers two alternatives that share the same library code path — select them in the input:

• Hencky (logarithmic) strain — ${\mathbf{E}}_{\mathrm{Hencky}}=\frac{1}{2}\ln ({\mathbf{F}}^{\mathsf{T}}\mathbf{F})$; correct physical behaviour but expensive.
• Bažant (2) — ${\mathbf{E}}_{\mathrm{Bazant}}=\frac{1}{4}({\mathbf{F}}^{\mathsf{T}}\mathbf{F}-({\mathbf{F}}^{\mathsf{T}}\mathbf{F}{)}^{-1})$; very close to Hencky at large deformations but without the logarithm. Practical default for compression-heavy simulations.

Hybrid parallelism

Every field-carrying module allocates via GridParameters and uses the halo-exchange machinery on BoundaryConditions. OpenMP is on by default; MPI is enabled with SETTINGS=mpi-parallel or -DENABLE_MPI=ON. 3D cartesian decomposition is selected by $MPI3D: Yes in the @GridParameters block, with $Ncx, $Ncy, $Ncz declaring the number of processes along each axis.

Spectral solvers (mechanics and electrics) switch to FFTW3-MPI automatically in MPI builds.

Copy-protection gating

The academic edition is built with -DNOSENTINEL and unlocked; the commercial edition uses the Sentinel licence library (external/Sentinel/). The ENABLE_SENTINEL CMake option gates this path; static linking is only supported when dynamic OpenMP + Sentinel is disabled (ENABLE_DYNAMIC_LINKING=ON).

See also

• Phase Field — the multiphase-field equations and the derivation of active-parameter tracking.
• Thermodynamics — driving-force derivation and partitioning relations.
• Mechanics — continuum-mechanics derivations for the hyperelastic solver and the strain-measure discussion.
• Interface Properties — the anisotropy model enumeration and the Arrhenius mobility coupling.