Verification confirms that the equations are solved correctly (solving the equations right), while validation assesses whether the equations adequately represent the physical system (solving the right equations) [172]. This
section presents verification studies using analytical solutions and the \(\mu \)MAG Standard Problem 4.
5.5.1 The \(\mu \)MAG Standard Problem 4
Standard Problem 4, defined by the Micromagnetic Modeling Activity Group (\(\mu \)MAG) [173], is a widely used benchmark for validating micromagnetic simulation codes. It considers a thin Permalloy film with
dimensions \(500 \times 125 \times 3\) \(\si {\nano \meter \cubed }\) and material parameters: Exchange stiffness \(A = 1.3 \times 10^{-11}\) J/m, Saturation
magnetization \(M_s = 8.0 \times 10^5\) A/m, Gilbert damping \(\alpha = 0.02\), Anisotropy \(K = 0\) (no crystalline anisotropy).
Figure 5.4: Setup for the \(\mu \)MAG Standard Problem 4. (a) Permalloy film geometry (\(500 \times 125 \times 3\) \(\si {\nano \meter \cubed }\)) discretized by an unstructured tetrahedral
mesh with 45391 elements. (b) Initial S-state magnetization configuration showing the characteristic S pattern in the \(xy\)-plane, colored by the \(x\)-component of magnetization.
The initial magnetization is an equilibrium S-state, prepared according to the \(\mu \)MAG specifications by applying and then slowly reducing a saturating field along the \((1,1,1)\) diagonal. To achieve this, a
saturating external field of \(\SI {1.0}{\tesla }\) is applied along the \((1,1,1)\) direction and held constant for \(\SI {0.5}{\nano \second }\) to align the magnetization. This field is then linearly ramped down to zero
over a period of \(\SI {3}{\nano \second }\), after which the system is allowed to relax in zero field for an additional \(\SI {2.5}{\nano \second }\). This procedure eliminates metastable states, yielding the characteristic
"S"-shaped domain pattern required for the benchmark. Following the relaxation, a field of sufficient magnitude to reverse the magnetization is applied. The simulation geometry, discretization, and equilibrium S-state
magnetization configuration are shown in Figure 5.4. The field is applied instantaneously at \(t=0\) to the equilibrium S-state:
The time evolution of the spatially averaged magnetization components, \(\langle m_x \rangle \), \(\langle m_y \rangle \), and \(\langle m_z \rangle \), is then examined as the system approaches equilibrium in the
new field. Reference solutions are available from multiple independent sources [173].
5.5.2 Solver Comparison on Standard Problem 4
Figure 5.5: Magnetization dynamics for Standard Problem 4 (Field 1) computed with different time integrators. All solvers agree with the \(\mu \)MAG reference solutions [174, 175].
Figure 5.5 compares the magnetization dynamics obtained with different solvers for Field 1 of Standard Problem 4. The results are validated against
reference solutions from d’Aquino et al. [174] and Albuquerque et al. [175], both submitted to the \(\mu \)MAG repository. All solvers produce trajectories in excellent
agreement with these reference solutions, validating the correctness of the implementation. The S-state relaxed to an average magnetization of \(\mathbf {m} = (m_x, m_y, m_z) = (0.96849366, 0.12573718,
-0.00088421495)\). The slight deviation from FDM reference values (\(m_x \approx 0.966\), \(m_y \approx 0.127\), \(m_z \approx 0.00000078\)) is attributed to differences in spatial discretization: the FEM provides a
more accurate representation of the magnetization configuration near edges and corners, whereas FDM solvers employ a staircase approximation that affects the local demagnetizing field calculation in these regions.
Figure 5.6: Deviation of computed magnetization trajectories from the \(\mu \)MAG reference solutions: (left) d’Aquino et al. [174], (right) Albuquerque et al.
[175]. All solvers exhibit comparable deviation from both references, with differences attributable to the finite element discretization rather than time integration errors.
Table 5.3 provides quantitative comparison metrics. The maximum deviation from the reference solutions is computed as:
where \(\langle \mathbf {m} \rangle _\text {ref}\) denotes the reference solutions of d’Aquino et al. [174] and Albuquerque et al. [175], with the corresponding
deviations shown in Figure 5.6.
Table 5.3: Quantitative comparison of solvers on \(\mu \)MAG Standard Problem 4. Errors are computed against both reference solutions.
.
\(e_\infty \)
Solver
Configuration
d’Aquino
Albuquerque
Steps
CPU [\(\si {\second }\)]
Speedup
TPS
\(\Delta t = 10^{-13}\,\si {\second }\)
0.10899498
0.12411681
20000
96752
\(1.00\times \)
TPS
\(\Delta t = 10^{-12}\,\si {\second }\)
0.31432466
0.25772080
2000
11962
\(8.09\times \)
RK4
\(\Delta t = 10^{-13}\,\si {\second }\)
0.09044942
0.12825486
20000
78813
\(1.23\times \)
BDF
.
rtol=\(10^{-2}\),
atol=\(10^{-4}\)
0.11185953
0.16341087
1099
40014
\(\mathbf {2.42\times }\)
IMEX
.
rtol=\(10^{-2}\),
atol=\(10^{-4}\)
0.08920000
0.12800000
3818
34510
\(\mathbf {2.80\times }\)
The computational efficiency results reveal that both adaptive implicit methods significantly outperform the fixed-step explicit schemes. The BDF method achieves a 2.4\(\times \) speedup over the TPS baseline despite requiring
approximately 25 linear solver iterations per Newton step. By taking time steps 10–50\(\times \) larger than the explicit CFL limit (1–5 \(\si {\pico \second }\) versus 0.1 \(\si {\pico \second }\)), the
total step count reduces from 20000 to just 1099, more than compensating for the increased per-step cost.
The IMEX method achieves the best overall performance, completing the simulation in 34510 \(\si {\second }\), which results in a 2.8\(\times \) speedup over TPS and 2.3\(\times \) over RK4. Although the time
step remains accuracy-limited to 0.3–0.5 \(\si {\pico \second }\) due to explicit precession (cf. Section 5.4.2), the SPD structure of the implicit
Jacobian enables highly efficient linear solves via AMG. This per-step efficiency advantage compensates for the larger step count (3818 versus 1099 for BDF), demonstrating that IMEX methods can be competitive when the implicit
partition is carefully chosen to yield favorable algebraic structure.
To investigate whether the TPS could achieve comparable efficiency by simply increasing the time step, we tested the scheme with \(\Delta t = 1\) \(\si {\pico \second }\), the same order of magnitude used by the
adaptive BDF solver. While the simulation remains stable (the TPS is unconditionally A-stable for \(\theta = 1\)), the solution exhibits significant phase drift, with \(e_\infty = 0.31\) against d’Aquino and \(e_\infty =
0.26\) against Albuquerque, more than double the error of any other solver. This demonstrates the critical distinction between stability and accuracy: an A-stable method guarantees bounded solutions
regardless of time step, but does not guarantee accurate solutions. The first-order truncation error of the TPS update (5.33) accumulates
as \(\mathcal {O}(N \Delta t) = \mathcal {O}(T)\), where \(N = T/\Delta t\) is the number of steps. Increasing \(\Delta t\) by a factor of 10 increases the per-step error by the same factor, resulting in the observed phase
lag visible in Figure 5.5.
This finding underscores the value of adaptive time-stepping methods. The BDF solver achieves large time steps (\(1\)–\(5\) \(\si {\pico \second }\)) while maintaining accuracy because the embedded error estimator
automatically reduces the step size when truncation error would exceed the specified tolerance. In contrast, a fixed-step method such as TPS requires the user to balance efficiency and accuracy manually, a balance that depends on
the specific problem and is difficult to determine a priori.
These results demonstrate that both fully implicit (BDF) and IMEX methods achieve significant speedup over explicit schemes for precession-dominated problems. The IMEX splitting following Di
Fratta et al. [176] and Mauser [166] successfully addresses exchange stiffness while maintaining SPD structure in the implicit partition, enabling efficient AMG solves that offset the
accuracy-limited time step.
5.5.3 Constraint Preservation and Solver Accuracy
A key aspect of numerical time integration for the LLG equation is the preservation of the unit-sphere constraint \(|\mathbf {m}| = 1\). Figure 5.7 (a) shows
the constraint deviation \(\| |\mathbf {m}| - 1 \|_\infty \) (maximum over all mesh nodes) before renormalization.
For explicit single-step methods, theoretical results establish that a scheme of order \(p\) preserves the magnetization magnitude within \(\mathcal {O}(\Delta t^{p+1})\) error per step [177, 178]. This favorable
behavior exploits the cross-product structure of the LLG equation, which guarantees \(\partial _t\mathbf {m} \perp \mathbf {m}\). Our RK4 results confirm this: with \(\Delta t = 10^{-13}\,\si {\second }\), the
fourth-order method yields per-step drift of \(\mathcal {O}(\Delta t^5) \sim 10^{-65}\), accumulating to \(\sim 10^{-9}\) over the simulation. The TPS method, being first-order in its tangent-plane update, produces
\(\mathcal {O}(\Delta t^2) \sim 10^{-26}\) drift per step, accumulating to \(\sim 10^{-5}\).
The implicit multistep methods (BDF and IMEX) exhibit fundamentally different drift characteristics. The \(\mathcal {O}(\Delta t^{p+1})\) bound does not apply because the multistep formulation combines solution values
from previous steps that each carry minor constraint violations, and the Newton iteration required for the implicit solve explores intermediate states that may lie far from the unit sphere. Unlike the RK4 method, the adaptive BDF
and IMEX schemes do not preserve the constraint order-by-order. Consequently, they exhibit larger oscillating deviations that correlate with the step size (Figure 5.7 (c)). Crucially, however, this norm deviation acts primarily as a scaling artifact rather than a physical error. In the unconstrained formulation, the effective field scales
linearly with magnetization magnitude (\(|\mathbf {H}| \propto |\mathbf {m}|\)). A deviation in norm (e.g., \(|\mathbf {m}| > 1\)) scales the update vector but preserves its direction. Since the projection step
\(\mathbf {m} \leftarrow \mathbf {m} / |\mathbf {m}|\) maps the solution back to the unit sphere radially, this length error is discarded without corrupting the angle error (physical precession). This
geometric property explains why the BDF and IMEX solvers maintain high trajectory accuracy (Figure 5.7 (b)) even when the intermediate norm deviation is
significant.
Figure 5.7: (a) Deviation from the unit-sphere constraint \(|\mathbf {m}|=1\) before renormalization. (b) Difference between solver trajectories and the TPS reference solution. (c) Time step evolution for
adaptive solvers (BDF, IMEX) compared to the fixed time step used by TPS and RK4.
All solvers apply post-step renormalization to restore \(|\mathbf {m}|=1\), ensuring the physical constraint is satisfied regardless of the intermediate drift. The constraint deviation generally decreases as the simulation progresses,
reflecting the reduction in magnetization velocity \(|\partial _t \mathbf {m}|\) as the system relaxes toward equilibrium.
The decreasing constraint deviation for fixed-step methods reflects the diminishing magnetization velocity \(|\partial _t \mathbf {m}| \to 0\) as the system relaxes, thereby reducing the per-step truncation error. In contrast,
the adaptive solvers maintain approximately constant deviation by increasing the time step to consume their full error budget, as evidenced by the growing \(\Delta t\) in Figure 5.7 (c).
Figure 5.7 (b) quantifies the deviation of each solver from the TPS reference solution. All solvers remain within acceptable error bounds, with differences
arising primarily from the time-discretization approach rather than from fundamental inconsistencies in the underlying physics. The deviation from the TPS solution confirms that all implemented solvers converge to the same
physical trajectory within the specified tolerances. Interestingly, the RK4, BDF, and IMEX solvers exhibit nearly identical deviations from TPS (\(e_\infty \approx 5 \times 10^{-2}\)), while differing from each other by only
\(\mathcal {O}(10^{-6})\). This reflects the different underlying formulations: TPS solves for the velocity in the tangent plane with exact constraint enforcement, whereas the other methods integrate the explicit LLG equation in
\(\mathbb {R}^3\) and then renormalize. Both approaches produce physically correct results that match the \(\mu \)MAG reference solutions within spatial discretization error.
5.5.4 Adaptive Time Stepping
A key advantage of the SUNDIALS-based solvers is their automatic time-step adaptation. Figure 5.7 (c) shows the time step history during the
Standard Problem 4 simulation, with the fixed time step of the TPS and RK4 solvers shown for reference.
The BDF solver achieves time steps of 1–5 \(\si {\pico \second }\) during the post-switching relaxation phase, representing a 10–50\(\times \) improvement over the explicit CFL limit. During the rapid switching
transient (\(t \approx 0.1\)–\(0.3\) \(\si {\nano \second }\)), the solver automatically reduces the time step to capture the fast dynamics, then smoothly increases it as the system relaxes. This smooth adaptation
profile reflects the fully implicit treatment: by solving for the magnetization at the future time level, the method effectively filters the high-frequency precession, allowing the error controller to respond to the slowly varying envelope.
The IMEX solver exhibits markedly different behavior, with oscillatory time steps that correlate with the magnetization dynamics. As \(|\mathbf {H}_\text {eff}|\) varies during precession, the local truncation error of the
explicit term fluctuates, causing the adaptive controller to adjust the step size continuously. This results in \(\Delta t \approx 0.15\)–\(0.3\,\si {\pico \second }\) with significant variation, and a step-rejection rate of
approximately 16%. The solver effectively chases the oscillating error estimate rather than settling into a stable large time step.
Note that the tolerances employed here (reltol=\(10^{-2}\), abstol=\(10^{-4}\)) are stricter than those commonly used in micromagnetic simulations. Relaxing to reltol=\(0\) and abstol=\(10^{-3}\) [158] permits larger
time steps at the expense of increased constraint violation, which may be acceptable when only qualitative switching behavior is of interest.
