Numerical Analysis and Innovative Simulation
Techniques for Designing Advanced MRAM

5.4 Adaptive Multistep Methods via SUNDIALS

The SUNDIALS (SUite of Nonlinear and DIfferential/ALgebraic equation Solvers) library  [161] provides implementations of adaptive time-stepping methods. Two solvers from this suite are integrated into the present framework: CVODE for implicit multistep methods and ARKStep for IMEX Runge–Kutta methods.

Both solvers operate on the ODE form \(d\mathbf {m}/dt = f(\mathbf {m})\). The FEM discretization developed in Section 4.2.2 transforms the LLG Equation (3.12) into the semi-discrete system (discrete in space, continuous in time):

\begin{equation} M \frac {d\mathbf {m}}{dt} = F(\mathbf {m}), \label {eq::TI_semidiscrete} \end{equation}

where \(M\) is the mass matrix with entries \(M_{ij} = \int _\Omega \phi _i \phi _j \, d\mathbf {x}\), and \(F(\mathbf {m})\) is the weak-form assembly of the LLG RHS:

\begin{equation} F_i(\mathbf {m}) = -\gamma _L \int _\Omega \bigl [ \mathbf {m} \times \mathbf {H}_\text {eff} + \alpha \, \mathbf {m} \times (\mathbf {m} \times \mathbf {H}_\text {eff}) \bigr ] \cdot \boldsymbol {\phi }_i \, d\mathbf {x}, \label {eq::TI_F_assembly} \end{equation}

with \(\mathbf {H}_\text {eff}\) as defined in (3.19). To obtain the explicit ODE form required by SUNDIALS, we employ row-sum mass lumping, replacing \(M\) with a diagonal matrix \(M^L\) having entries \(M_{ii}^L = \sum _j M_{ij}\). This approximation preserves the unit-length constraint at each node and has been proven to converge to a weak solution of the LLG equation  [165, 166]. This yields the nodewise system:

\begin{equation} \frac {d\mathbf {m}}{dt} = f(\mathbf {m}) = (M^L)^{-1} F(\mathbf {m}), \label {eq::TI_mass_lumped_ODE} \end{equation}

enabling efficient explicit evaluation and nodewise constraint enforcement \(|\mathbf {m}_i| = 1\) at each mesh node  [90].

5.4.1 CVODE: Variable-Order BDF

CVODE implements variable-order, variable-step BDF (cf. Section 5.2.3) with automatic error control. At each step, the solver performs a predictor-corrector sequence:

• (i) Prediction: An initial guess \(\mathbf {m}^{(0)}_{n+1}\) is constructed via polynomial extrapolation from the history array.

• (ii) Correction: The predicted solution is refined by solving the nonlinear residual equation \(G(\mathbf {m}_{n+1}) = 0\) using Newton iteration. The residual is defined as:

  \(\seteqnumber{0}{5.}{37}\)

  \begin{equation} G(\mathbf {m}) = \sum _{j=0}^{q} \alpha _j \mathbf {m}_{n+1-j} - \Delta t\, \beta _0\, f(\mathbf {m}), \label {eq::TI_residual_G} \end{equation}

  with coefficients \(\alpha _j\), \(\beta _0\) as defined in (5.16).

• (iii) Adaptation: The local truncation error is estimated to dynamically adjust both the step size \(\Delta t\) and the method order \(q\) (typically \(1 \le q \le 5\)). To ensure unconditional stability for the purely imaginary eigenvalues characteristic of the exchange interaction, the maximum BDF order is restricted to \(q=2\), ensuring A-stability (similar to the strategy in  [158]). This order restriction is essential for the oscillatory dynamics of the LLG equation. The precessional term generates eigenvalues \(\lambda \approx \pm i\omega \) lying close to the imaginary axis, where \(\omega = \gamma _0 |\mathbf {H}_\text {eff}|\) is the Larmor precession frequency. While BDF1 and BDF2 are A-stable (their stability regions include the entire left half-plane), BDF3 and higher orders exclude wedge-shaped regions near the imaginary axis  [156]. When these precessional eigenvalues fall outside the stability region, the adaptive solver is forced to reduce the time step to maintain stability rather than accuracy, thereby negating the efficiency gains of implicit integration. Restricting to order \(q \leq 2\) guarantees unconditional stability for all eigenvalue configurations encountered in micromagnetic simulations.

The core computational challenge lies in the Newton correction. Let \(\mathbf {m}^{(k)}\) denote the \(k\)-th iteration. The update \(\delta \mathbf {m} = \mathbf {m}^{(k+1)} - \mathbf {m}^{(k)}\) is obtained by solving the linear system:

\begin{equation} (I - \gamma J) \delta \mathbf {m} = -G(\mathbf {m}^{(k)}), \label {eq::TI_CVODE_Newton} \end{equation}

where \(\gamma = \Delta t \beta _0\) and \(J = \partial f / \partial \mathbf {m}\) is the system Jacobian. For the discretized LLG equation, \(J\) is a dense, non-symmetric matrix of dimension \(3N \times 3N\). Direct factorization of \(I - \gamma J\) is computationally prohibitive (\(\mathcal {O}(N^3)\)). Instead, CVODE employs preconditioned Krylov subspace methods (GMRES), which approximate the solution using only matrix-vector products.

The system Jacobian \(J = \partial f / \partial \mathbf {m}\) inherits contributions from each term in the effective field (3.19):

\begin{equation} J = J_\text {exch} + J_\text {aniso} + J_\text {demag} + J_\text {STT} + \ldots \label {eq::TI_J_decomposition} \end{equation}

The exchange contribution \(J_\text {exch}\), arising from the Laplacian in \(\mathbf {H}_\text {exch}\) (3.27), dominates the spectral radius with eigenvalues scaling as \(\mathcal {O}(h^{-2})\). The demagnetization Jacobian \(J_\text {demag}\) is dense due to the nonlocal nature of \(\mathbf {H}_\text {demag}\), while \(J_\text {aniso}\) is sparse and bounded. This structure motivates a preconditioner that captures the stiff exchange modes while neglecting the expensive demagnetization contribution.

Preconditioning Strategy: The spectral distribution of the iteration matrix governs the convergence rate of GMRES. Without preconditioning, the stiffness of the exchange interaction leads to eigenvalues scaling as \(\mathcal {O}(h^{-2})\) distributed along both the negative real axis (damping) and the imaginary axis (precession). This results in a significant condition number, causing the linear solver to stagnate.

Suess et al.  [158] proposed a physics-based preconditioner within a finite-element framework employing the box scheme  [167] for spatial discretization on tetrahedral meshes. One starts from the exchange contribution to \(f\):

\begin{equation} f_{\text {exch}} = -\gamma _L \mathbf {m} \times \nabla ^2 \mathbf {m} - \gamma _L \alpha \, \mathbf {m} \times (\mathbf {m} \times \nabla ^2 \mathbf {m}), \label {eq::TI_f_exch} \end{equation}

where \(\gamma _L\) denotes the gyromagnetic ratio. Linearizing around the current magnetization \(\mathbf {m}\) and using \(|\mathbf {m}|=1\) yields:

\begin{equation} \frac {\partial f_{\text {exch}}}{\partial \mathbf {m}} \approx -\gamma _L [\mathbf {m}]_\times \nabla ^2 + \gamma _L \alpha \, \nabla ^2, \label {eq::TI_Jexch_linearization} \end{equation}

where \([\mathbf {m}]_\times \) denotes the skew-symmetric matrix representing the cross product with \(\mathbf {m}\):

\begin{equation} [\mathbf {m}]_\times = \begin{pmatrix} 0 & -m_z & m_y \\ m_z & 0 & -m_x \\ -m_y & m_x & 0 \end {pmatrix}, \quad \text {such that } [\mathbf {m}]_\times \mathbf {v} = \mathbf {m} \times \mathbf {v}. \label {eq::TI_cross_matrix} \end{equation}

The original Suess preconditioner takes the form:

\begin{equation} P_{\text {Suess}} = M + \gamma \alpha C_{\text {exch}} K_{\text {diff}} + \gamma C_{\text {exch}} K_{\text {skew}}, \label {eq::TI_Suess_Prec} \end{equation}

where \(M\) is the mass matrix, \(C_{\text {exch}} = 2\gamma _L A_{\text {ex}}/(\mu _0 M_s)\) is the exchange coefficient, and:

• (i) \(K_{\text {diff}}\): Symmetric positive-definite matrix arising from the damping term, mathematically equivalent to a vector diffusion operator.

• (ii) \(K_{\text {skew}}\): Skew-symmetric matrix arising from the precessional term \(-\mathbf {m} \times \nabla ^2 \mathbf {m}\), capturing the rotation of the magnetization.

In the box scheme, the effective field at each node \(i\) is computed as \(H^k_{i,\text {eff}} \approx -m_i^{-1}\,\partial E_t / \partial u^k_i\), where \(m_i = \int _{V_i} J_s(\mathbf {x})\,\mathrm {d}V\) is the magnetic moment associated with the control volume \(V_i\) surrounding node \(i\). This yields a naturally diagonal weighting with entries \(m_i\), enabling the second derivatives of the total energy to be computed analytically on an element-by-element basis, from which the system Jacobian can be constructed explicitly  [158]. The stray field contribution is omitted from the preconditioner to preserve sparsity. The resulting preconditioner system is solved using Biconjugate Gradient Stabilized (BiCGSTAB) method  [168] with relaxed incomplete LU (RILU) factorization  [169].

Adaptation from Box Scheme to Galerkin Finite Elements: In contrast to the box scheme employed by Suess et al., the present work uses a standard Galerkin finite-element discretization with row-sum mass lumping (cf. Section 5.5.4). This change in spatial discretization requires adapting the preconditioner construction. The semi-discrete system (5.35) involves the consistent mass matrix \(M\), which is sparse but non-diagonal. Following standard practice  [90, 165], we employ row-sum mass lumping to obtain a diagonal approximation \(M^L\) with entries \(M^L_{ii} = \sum _j M_{ij}\). The preconditioner is then constructed in the dimensionless form appropriate for the lumped system:

\begin{equation} P = I + \gamma (M^L)^{-1} \bigl ( \alpha C_{\text {exch}} K_{\text {diff}} + C_{\text {exch}} K_{\text {skew}} \bigr ). \label {eq::TI_Prec_lumped} \end{equation}

In both the box scheme and the lumped Galerkin formulation, the diagonal entries (\(m_i\) or \(M^L_{ii}\), respectively) vary across nodes on non-uniform meshes. However, in the box scheme, the \(m_i^{-1}\) weighting is intrinsic to the effective field computation, and the Jacobian is assembled directly in this weighted form. In the Galerkin approach, the weighting arises as an a posteriori left-multiplication \((M^L)^{-1} K_{\text {skew}}\) that destroys the skew-symmetric structure of the assembled stiffness matrix. Specifically, letting \(D = (M^L)^{-1}\) denote the non-uniform diagonal scaling, the product \(D K_{\text {skew}}\) is no longer skew-symmetric:

\begin{equation} (D K_{\text {skew}})^T = K_{\text {skew}}^T D = -K_{\text {skew}} D \neq -D K_{\text {skew}}, \label {eq::TI_skew_destruction} \end{equation}

unless \(D\) is a scalar multiple of the identity. The resulting matrix is neither symmetric nor skew-symmetric, with eigenvalues scattered in the complex plane. This causes BiCGSTAB to stagnate or diverge, particularly at larger time steps, where the skew contribution dominates  [170]. For this reason, the full Suess preconditioner with the skew-symmetric term is not used in this work.

Practical Implementation: Instead, we adopt a simplified preconditioner that omits the problematic skew-symmetric term entirely. The original Suess preconditioner  [158] was developed using the box scheme for spatial discretization, in which the \(m_i^{-1}\) weighting is incorporated during Jacobian assembly rather than applied as a post-hoc scaling. In the present Galerkin formulation, the analogous scaling \((M^L)^{-1}\) is applied after assembly of the global stiffness matrices, which disrupts the skew-symmetric structure of the precessional contribution. Standard Galerkin finite-element discretization with mass lumping produces a diagonal matrix \(M^L\) whose entries vary significantly on unstructured meshes. Following the same approach as Suess et al.  [158], who validated the original preconditioner against the \(\mu \)MAG standard problems, the simplified preconditioner adopted in this work reads:

\begin{equation} P = I + \gamma \alpha C_{\text {exch}} (M^L)^{-1} K_{\text {diff}}. \label {eq::TI_Prec_damping_only} \end{equation}

This matrix is symmetric positive-definite, enabling the use of the PCG method with AMG preconditioning via BoomerAMG from the HYPRE library  [171]. The AMG solver is configured for 3D vector problems using the systems interface with nodal coarsening. A loose relative tolerance of \(10^{-1}\) is employed, as the preconditioner needs only to provide an approximate solution as the outer GMRES iteration in CVODE corrects any residual error.

While this preconditioner does not capture the precessional dynamics, it effectively clusters the stiff exchange eigenvalues, reducing the GMRES iteration count from \(\mathcal {O}(100)\) (unpreconditioned) to \(\mathcal {O}(20)\) per Newton step. The trade-off, approximately 30–40% more linear iterations compared to the full Suess preconditioner, is acceptable given the guaranteed convergence and robustness on meshes of arbitrary quality. The correctness of this approach is validated against the \(\mu \)MAG Standard Problem 4 in Section 5.5.2, where all solvers using this preconditioner reproduce the published reference solutions. Table 5.1 summarizes the preconditioner variants and their applicability.

Alternative: Incomplete LU Factorization. An alternative approach employs incomplete LU (ILU) factorization, which directly factorizes the non-symmetric preconditioner matrix without requiring structural modifications. Unlike AMG, which relies on smoothing properties that degrade for skew-symmetric operators, ILU factorization handles indefinite systems naturally by computing approximate \(L\) and \(U\) factors. The Euclid parallel ILU implementation in HYPRE  [171] provides a scalable option for moderate problem sizes. With ILU preconditioning, the full preconditioner including \(K_{\text {skew}}\) and \(J_{\text {STT}}\) can be employed:

\begin{equation} P_{\text {ILU}} = I + \gamma \alpha C_{\text {exch}} (M^L)^{-1} K_{\text {diff}} + \gamma C_{\text {exch}} (M^L)^{-1} K_{\text {skew}} + \gamma (M^L)^{-1} J_{\text {STT}}. \label {eq::TI_Prec_ILU} \end{equation}

This approach recovers the full Jacobian approximation at the cost of increased fill-in compared to AMG. For problems where STT dominates, the improved Jacobian approximation can reduce the number of outer GMRES iterations sufficiently to offset the higher preconditioner cost.

Extension for Spin-Transfer Torque: The standard preconditioner assumes the system is driven towards equilibrium by damping. However, current-induced switching introduces a non-conservative STT that acts as negative damping, pumping energy into the system and destabilizing the initial state.

In the switching regime, the STT terms dominate the dynamics. If the preconditioner ignores these terms, the linear solver will operate with an incorrect model of the system stiffness, leading to slow convergence. The total derivative of the STT with respect to magnetization is:

\begin{equation} \frac {d\boldsymbol {\tau }_{\text {STT}}}{d\mathbf {m}} = \frac {\partial \boldsymbol {\tau }_{\text {STT}}}{\partial \mathbf {m}}\bigg |_{\mathbf {S}} + \frac {\partial \boldsymbol {\tau }_{\text {STT}}}{\partial \mathbf {S}} \cdot \frac {\partial \mathbf {S}}{\partial \mathbf {m}}. \label {eq::TI_Full_Jacobian} \end{equation}

The second term represents the feedback of magnetization changes on the spin-accumulation through the drift-diffusion system. Computing this term exactly would require differentiating through the solution of the spin transport PDE, effectively solving an adjoint problem for each magnetization degree of freedom, rendering it computationally prohibitive.

We therefore adopt the frozen-\(\mathbf {S}\) approximation:

\begin{equation} J_{\text {STT}} \approx \frac {\partial \boldsymbol {\tau }_{\text {STT}}}{\partial \mathbf {m}}\bigg |_{\mathbf {S}=\text {const}}. \label {eq::TI_Frozen_S} \end{equation}

Explicitly, for the damping-like torque \(\boldsymbol {\tau }_{\text {DL}} = c_{\text {DL}}\, \mathbf {m} \times (\mathbf {m} \times \mathbf {S})\) and field-like torque \(\boldsymbol {\tau }_{\text {FL}} = c_{\text {FL}}\, \mathbf {m} \times \mathbf {S}\):

\begin{align} \frac {\partial \boldsymbol {\tau }_{\text {DL}}}{\partial \mathbf {m}}\bigg |_{\mathbf {S}} &= -c_{\text {DL}} \bigl [ \mathbf {m} \otimes \mathbf {S} + (\mathbf {m} \cdot \mathbf {S})\, I \bigr ], \label {eq::TI_J_DL}\\[4pt] \frac {\partial \boldsymbol {\tau }_{\text {FL}}}{\partial \mathbf {m}}\bigg |_{\mathbf {S}} &= c_{\text {FL}}\, [\mathbf {S}]_\times . \label {eq::TI_J_FL} \end{align} Where \([\mathbf {S}]_\times \) denotes the skew-symmetric cross-product matrix of \(\mathbf {S}\). This approximation neglects the feedback term \((\partial \boldsymbol {\tau }/\partial \mathbf {S})(\partial \mathbf {S}/\partial \mathbf {m})\) in the Jacobian. Physically, this relies on the same time-scale separation that motivates the steady-state assumption \(\partial _t \mathbf {S} = 0\) in the drift-diffusion formulation (see Section 6.1): the spin-flip relaxation time \(\tau _\mathrm {sf}\) is of the order of picoseconds, whereas magnetization dynamics evolve on nanosecond time scales  [35]. Consequently, \(\mathbf {S}\) equilibrates quasi-instantaneously with respect to \(\mathbf {m}\) and can be recomputed from \(\mathbf {m}\) at the beginning of each time step and held constant during the Newton iteration that advances \(\mathbf {m}\).

Both the damping-like torque contribution (5.51) and the field-like term (5.52) are incorporated into the preconditioner. While the damping-like term \(J_{\text {DL}}\) is symmetric, the field-like contribution \(J_{\text {FL}} = c_{\text {FL}}[\mathbf {S}]_\times \) is skew-symmetric. This skew-symmetry introduces the same structural problem encountered with \(K_{\text {skew}}\): AMG smoothers, designed for elliptic operators, degrade when applied to matrices with significant skew-symmetric components.

Consequently, when spin-transfer torque is active in the BDF solver, the preconditioner employs ILU factorization rather than AMG. Unlike algebraic multigrid, ILU directly factorizes the non-symmetric matrix without relying on smoothing properties, handling the skew-symmetric \(J_{\text {FL}}\) naturally. The Euclid parallel ILU implementation from HYPRE  [171] is employed with fill level \(k=1\).

The extended preconditioner for BDF with spin-transfer torque reads:

\begin{equation} P_{\text {BDF}} = I + \gamma \alpha C_{\text {exch}} (M^L)^{-1} K_{\text {diff}} - \gamma (M^L)^{-1} J_{\text {STT}}, \label {eq::TI_Prec_with_STT} \end{equation}

where \(J_{\text {STT}} = J_{\text {DL}} + J_{\text {FL}}\) comprises the damping-like (symmetric) and field-like (skew-symmetric) Jacobian contributions defined in Equations (5.51) and (5.52). Note that \(K_{\text {skew}}\) is omitted from the preconditioner despite appearing in the true Jacobian: numerical experiments showed that including \(K_{\text {skew}}\) after mass lumping degraded convergence rather than improving it (cf. the discussion following Equation (5.46)). The outer linear system is solved using GMRES preconditioned by ILU.

The numerical validation of this extended preconditioner for spin-transfer torque switching is presented in Section 6.8, following the derivation of the spin-accumulation model.

5.4.2 ARKStep: IMEX Runge–Kutta

ARKStep implements Additive Runge–Kutta (ARK) methods for the IMEX splitting introduced in Section 5.2.4. The effectiveness of IMEX strategies for micromagnetic simulations, treating exchange implicitly while evaluating lower-order contributions explicitly, has been analyzed rigorously by Mauser et al.  [166]. The formulation applies diagonally implicit tableaux to \(f_I\) and explicit tableaux to \(f_E\) simultaneously:

\begin{align} \mathbf {k}_i^E &= f_E\left ( t_n + c_i \Delta t,\; \mathbf {m}_n + \Delta t \sum _{j=1}^{i-1} a_{ij}^E \mathbf {k}_j^E + \Delta t \sum _{j=1}^{i-1} a_{ij}^I \mathbf {k}_j^I \right ), \label {eq::TI_ARK_E}\\ \mathbf {k}_i^I &= f_I\left ( t_n + c_i \Delta t,\; \mathbf {m}_n + \Delta t \sum _{j=1}^{i-1} a_{ij}^E \mathbf {k}_j^E + \Delta t \sum _{j=1}^{i} a_{ij}^I \mathbf {k}_j^I \right ), \label {eq::TI_ARK_I} \end{align} with the solution advanced via:

\begin{equation} \mathbf {m}_{n+1} = \mathbf {m}_n + \Delta t \sum _{i=1}^{s} \bigl ( b_i^E \mathbf {k}_i^E + b_i^I \mathbf {k}_i^I \bigr ), \label {eq::TI_ARK_update} \end{equation}

where \(b_i^E\), \(b_i^I\) are the quadrature weights from the respective Butcher tableaux.

Implicit Preconditioning: The computational bottleneck in IMEX schemes is the solution of the nonlinear systems arising in the implicit stages (5.55). For the LLG equation, the implicit partition \(f_I\) (5.19) isolates the exchange damping term. Because the exchange field depends linearly on \(\mathbf {m}\), the Jacobian of the implicit function in weak form is:

\begin{equation} J_{\text {IMEX}} = M - \gamma \frac {\partial f_I}{\partial \mathbf {m}} = M + \gamma \alpha C_{\text {exch}} K_{\text {stiff}}, \label {eq::TI_JIMEX} \end{equation}

where \(M\) is the mass matrix and \(K_{\text {stiff}}\) is the stiffness matrix associated with the vector Laplacian of the magnetization. Applying mass lumping (\(M \to M^L\), \(M^L_{ii} = \sum _j M_{ij}\)) and left-multiplying by \((M^L)^{-1}\) yields the dimensionless form suitable for iterative solution:

\begin{equation} (M^L)^{-1} J_{\text {IMEX}} \approx I + \gamma \alpha C_{\text {exch}} (M^L)^{-1} K_{\text {stiff}}. \label {eq::TI_JIMEX_lumped} \end{equation}

Unlike the full Newton Jacobian in CVODE (5.53), this operator is symmetric positive-definite (SPD) when the precessional cross-product and spin-transfer torque are moved to \(f_E\). This structural property enables highly efficient linear solvers such as PCG with AMG preconditioning, which are optimal for elliptic diffusion problems. The extension to spin-transfer torque simulations preserves this advantage: since \(\mathbf {J}_{\text {STT}}\) does not appear in the implicit partition, it is omitted from the IMEX preconditioner, maintaining the SPD structure regardless of spin-torque magnitude.

However, a fundamental limitation arises from the explicit treatment of precession. The partition \(f_E\) contains the full precession term \(-\gamma _L \mathbf {m} \times \mathbf {H}_\text {eff}\), which oscillates at the Larmor frequency \(\omega = \gamma _0 |\mathbf {H}_\text {eff}| \sim 10^{10}\)–\(10^{11}\,\si {\radian \per \second }\). While the IMEX splitting successfully removes the exchange-driven CFL stability restriction, the time step remains constrained by the truncation error of this explicit term. To accurately resolve the rapid oscillations, the local error must remain small, requiring \(\Delta t \lesssim C/\omega \sim 0.1\)–\(1\) \(\si {\pico \second }\), comparable to the explicit stability limit. This accuracy constraint, rather than stability, becomes the limiting factor.

Evidence for this behavior appears in the ARKStep solver statistics: throughout our Standard Problem 4 simulations, num_exp_limited_stepfails \(= 0\) (no stability-induced rejections) while num_acc_limited_steps equals the total step count (all steps accuracy-limited). The solver never encounters the explicit stability boundary, while the truncation error of the explicit precession term determines every time step.

A further consequence is oscillatory time-step behavior. As the magnetization precesses, the effective field magnitude \(|\mathbf {H}_\text {eff}|\) varies periodically, causing the truncation error and hence the allowed time step to fluctuate, see Figure 5.7 (c). The adaptive controller continuously adjusts the step size, chasing the oscillating error estimate. This contrasts sharply with the smooth adaptation observed in fully implicit BDF methods, where the implicit treatment of precession filters the high frequency dynamics, allowing the error controller to respond to the slowly varying solution. Table 5.2 summarizes the time integrators implemented in this work.