The extrapolated explicit midpoint scheme for variable order and step size controlled integration of the Landau-Lifschitz-Gilbert equation

TL;DR

This paper introduces an efficient, adaptive higher-order integration scheme for the Landau-Lifschitz-Gilbert equation, utilizing extrapolation of the explicit midpoint rule and a piecewise stray field approximation to improve computational efficiency.

Contribution

It presents a novel extrapolated explicit midpoint scheme with adaptive order and step size control for LLG equation integration, incorporating a new stray field approximation method.

Findings

The scheme achieves higher efficiency compared to traditional methods.

Numerical experiments confirm the method's accuracy and computational savings.

Performance improves when stray field computation dominates the cost.

Abstract

A practical and efficient scheme for the higher order integration of the Landau-Lifschitz-Gilbert (LLG) equation is presented. The method is based on extrapolation of the two-step explicit midpoint rule and incorporates adaptive time step and order selection. We make use of a piecewise time-linear stray field approximation to reduce the necessary work per time step. The approximation to the interpolated operator is embedded into the extrapolation process to keep in step with the hierarchic order structure of the scheme. We verify the approach by means of numerical experiments on a standardized NIST problem and compare with a higher order embedded Runge-Kutta formula. The efficiency of the presented approach increases when the stray field computation takes a larger portion of the costs for the effective field evaluation.

Equations

$$\frac{\partial \mathbf{M}}{\partial t}$$

$$\mathbf{e}(t;h):=\mathbf{\eta }(t;h)-\mathbf{M}(t)={\mathbf{c}}_1(t)h^{\beta p}+{\mathbf{c}}_2(t)h^{\beta (p+1)}+\mathcal{O}(h^{\beta (p+2)}),$$

$$\mathbf{\eta }(t;h,h/2):=\mathbf{\eta }(t;h/2)+\frac{\mathbf{\eta }(t;h/2)-\mathbf{\eta }(t;h)}{2^{\beta p}-1}.$$

$$\mathbf{\eta }(t;h,h/2)=\mathbf{M}(t)+\mathcal{O}(h^{\beta (p+1)}).$$

$$\displaystyle e_{tot}(\boldsymbol{m})=\frac{1}{|\Omega|}\Big{(}\frac{A}{\mu_{0}M_{s}^{2}}\int_{\Omega}|\nabla\boldsymbol{m}|^{2}\,\text{d}x-\frac{1}{2}\int_{\Omega}\boldsymbol{m}\cdot\boldsymbol{h}_{s}(\boldsymbol{m})\,\text{d}x-\frac{K_{1}}{\mu_{0}M_{s}^{2}}\int_{\Omega}(\boldsymbol{a}\cdot\boldsymbol{m})^{2}\,\text{d}x-\int_{\Omega}\boldsymbol{m}\cdot\boldsymbol{h}_{ext}\,\,\text{d}x\Big{)},$$

$$\displaystyle\frac{\partial\boldsymbol{m}}{\partial\tau}=-\frac{1}{1+\alpha^{2}}\,\boldsymbol{m}\times\boldsymbol{h}(\boldsymbol{m})-\frac{\alpha}{1+\alpha^{2}}\,\boldsymbol{m}\times\big{(}\boldsymbol{m}\times\boldsymbol{h}(\boldsymbol{m})\big{)},$$

$$\mathbf{h}(\mathbf{m}):=-|\Omega |\frac{\delta e_{tot}}{\delta m}=\frac{2A}{{\mu }_0{M}_s^2}\Delta \mathbf{m}+{\mathbf{h}}_s(\mathbf{m})+\frac{2K_1}{{\mu }_0{M}_s^2}(\mathbf{a}\cdot \mathbf{m})\mathbf{a}+{\mathbf{h}}_{ext}.$$

$$\mathbf{m}(x,0)={\mathbf{m}}^{(0)}(x)$$

$$\frac{\partial \mathbf{m}}{\partial \mathbf{\nu }}=0,x\in \partial \Omega ,$$

$$\displaystyle\begin{aligned} \frac{\partial\boldsymbol{m}_{j}}{\partial\tau}&\,=-\frac{1}{1+\alpha^{2}}\,\boldsymbol{m}_{j}\times\boldsymbol{h}_{j}(\boldsymbol{m})-\frac{\alpha}{1+\alpha^{2}}\,\boldsymbol{m}_{j}\times\big{(}\boldsymbol{m}_{j}\times\boldsymbol{h}_{j}(\boldsymbol{m})\big{)},\quad j=1,\ldots,N,\\ \boldsymbol{m}(0)&\,=\boldsymbol{m}^{(0)},\end{aligned}$$

$$\displaystyle\begin{aligned} \boldsymbol{m}^{\prime}(t)&\,=F\big{(}\boldsymbol{m}(t)\big{)},\\ \boldsymbol{m}(t_{0})&\,=\boldsymbol{m}^{(0)}.\end{aligned}$$

$$\displaystyle\begin{aligned} \boldsymbol{m}^{(1)}&\,=\boldsymbol{m}^{(0)}+h\,F(\boldsymbol{m}^{(0)})\\ \boldsymbol{m}^{(\nu+1)}&\,=\boldsymbol{m}^{(\nu-1)}+2h\,F(\boldsymbol{m}^{(\nu)})\quad\nu=1,2,\ldots,n\\ \boldsymbol{m}_{h}(t)&\,=\frac{1}{4}(\boldsymbol{m}^{(n-1)}+2\boldsymbol{m}^{(n)}+\boldsymbol{m}^{(n+1)})=\frac{1}{2}\big{(}\boldsymbol{m}^{(n)}+\boldsymbol{m}^{(n-1)}+h\,F(\boldsymbol{m}^{(n)})\big{)}.\end{aligned}$$

$$\begin{array}{rl}\mathbf{m}(t_0+\nu h)-{\mathbf{m}}_h(t_0+\nu h)& =\sum _{j=1}^{\ell }{\mathbf{a}}_j(t_0+\nu h)h^{2j}+h^{2j+2}\mathbf{c}(t_0+\nu h;h),\nu \text{even}\\ \mathbf{m}(t_0+\nu h)-{\mathbf{m}}_h(t_0+\nu h)& =\sum _{j=1}^{\ell }{\mathbf{b}}_j(t_0+\nu h)h^{2j}+h^{2j+2}\tilde{\mathbf{c}}(t_0+\nu h;h),\nu \text{odd},\end{array}$$

$${\mathbf{\mu }}_{\ell ,m}={\mathbf{\mu }}_{\ell ,m-1}+\frac{{\mathbf{\mu }}_{\ell ,m-1}-{\mathbf{\mu }}_{\ell -1,m-1}}{(n_{\ell }/n_{\ell -m+1}{)}^2-1}m=1,\dots ,\ell .$$

$${\epsilon }_{\ell ,m-1}:=\Vert {\mathbf{\mu }}_{\ell ,m-1}-{\mathbf{\mu }}_{\ell ,m}\Vert =\mathcal{O}(H^{2m-1}).$$

$$\displaystyle\begin{aligned} \boldsymbol{m}^{\prime}(t)&\,=F\big{(}\boldsymbol{m}(t);\mathcal{D}\big{)},\\ \boldsymbol{m}(0)&\,=\boldsymbol{m}^{(0)},\end{aligned}$$

$${\mathcal{D}}_H\mathbf{m}(\nu h_1)$$

$$\displaystyle\begin{aligned} \boldsymbol{m}^{\prime}(t)&\,=F\big{(}\boldsymbol{m}(t);\mathcal{D}_{H}\big{)},\\ \boldsymbol{m}(0)&\,=\boldsymbol{m}^{(0)},\end{aligned}$$

$${\mathcal{D}}_H{\mathbf{m}}^{(\nu )}$$

$${\mathbf{m}}_1^{(n_1/2)}$$

$$\mathcal{D}{\mathbf{m}}_1^{(n_1/2)}$$

$$\displaystyle H_{\ell}=0.94\cdot H\cdot\Big{(}0.65\,\frac{tol}{err_{\ell}}\Big{)}^{1/(2\ell-1)}.$$

$$W_{\ell }=f_{sf}(2\ell +1)+(1-f_{sf})(1+\sum _{\nu =1}^{\ell }{n}_{\nu })=f_{sf}(2\ell +1)+(1-f_{sf})(2^{\ell +1}-1),$$

$$⟨\alpha ⟩=-h\frac{{\sum }_j\Delta {\epsilon }_j^{(\nu )}}{{\sum }_j|\Delta {\mathbf{m}}_j^{(\nu )}{|}^2},$$

$$\alpha =-\frac{{\int }_{\Omega }\text{d}\epsilon /\text{d}\tau }{{\int }_{\Omega }(\text{d}\mathbf{m}/\text{d}\tau {)}^2},$$

Full text

The extrapolated explicit midpoint scheme for variable order and step size controlled integration of the Landau-Lifschitz-Gilbert equation

Lukas Exl lukas.exl@univie.ac.at Fak. Mathematik, Univ. Wien, 1090 Vienna, Austria.

Inst. of Solid State Physics, TU Wien, 1040 Vienna, Austria.

Norbert J. Mauser

Wolfgang Pauli Institute c/o Fak. Mathematik, Univ. Wien, 1090 Vienna, Austria.

Thomas Schrefl

Center of Integrated Sensor Systems, Danube Univ. Krems, 2700 Wr. Neustadt, Austria

Dieter Suess

CD-Laboratory for Advanced Magnetic Sensing and Materials, TU Wien, 1040 Vienna, Austria

Abstract. A practical and efficient scheme for the higher order integration of the Landau-Lifschitz-Gilbert (LLG) equation is presented. The method is based on extrapolation of the two-step explicit midpoint rule and incorporates adaptive time step and order selection. We make use of a piecewise time-linear stray field approximation to reduce the necessary work per time step. The approximation to the interpolated operator is embedded into the extrapolation process to keep in step with the hierarchic order structure of the scheme. We verify the approach by means of numerical experiments on a standardized NIST problem and compare with a higher order embedded Runge-Kutta formula. The efficiency of the presented approach increases when the stray field computation takes a larger portion of the costs for the effective field evaluation.

Keywords: Landau-Lifschitz-Gilbert equation, extrapolation method, explicit midpoint scheme, variable

order method, micromagnetics

1 Introduction

Micromagnetics is a continuum theory of ferromagnetic materials located between classical Maxwell’s theory of electromagnetism and quantum theory [1]. A ferromagnetic system is described by the total magnetic energy of its magnetic distribution, which is modeled as a continuous vector field within the magnetic material. Typical length scales, that can be resolved by micromagnetic models, are in the range of a few nanometers to micrometers, which is too large for atomistic spin dynamics. On the other hand, these length scales are large enough for computer simulations of magnetic data storage systems like hard discs [2, 3] or random access memory [4] and high performance permanent magnets [5, 6].

The fundamental equation for dynamic processes of the magnetization in a magnetic body $\Omega\subset\mathbb{R}^{3}$, a vector field $\boldsymbol{M}(x,t)=M_{s}\boldsymbol{m}(x,t),\,|\boldsymbol{m}(x,t)|=1$ depending on the position $x\in\Omega$ and time $t\in\mathbb{R}$, is the Landau-Lifschitz-Gilbert equation. It is given in explicit form as [7]

[TABLE]

where $\gamma_{0}$ is the gyromagnetic ratio, $\alpha$ the damping constant and $\boldsymbol{H}$ the effective field, which is the sum of nonlocal and local fields such as the stray field and the exchange field, respectively.

Typically, equation (1) is numerically treated by a spatial semi-discrete approach [8, 9, 10]. We mention here, that in recent years also lower order finite element methods for the LLG equation were developed along with convergence analysis of weak solutions [11, 12, 13]. The computational main difficulty for numerics of the LLG equation arises from the expensive right hand side evaluation, mostly due to the nonlocal part in the effective field, namely the stray field. Several numerical methods were developed for the stray field calculation [14]. They either rely on a scalar potential or a field-based approach and scale at best linearly with the number of discrete magnetic spins or computational units. Nevertheless, the amount of computational costs for this calculation is typically $80-90\%$ of that of the (total) effective field. Hence, it is desirable to develop numerical schemes for micromagnetics that try to avoid excessive field evaluations, while also maintaining accuracy and efficiency. In the large $\alpha$ case, equation (1) degenerates to a steepest descent method for minimizing the total energy, owing to the Lyapunov structure of the LLG equation [8]. In this case, steepest descent methods [15] and conjugate gradient variants [16] were recently developed, which already require fairly optimal amounts of field evaluations. In moderately damped cases high accuracy and large time steps can be achieved by either higher order non-stiff integrators, as Runge-Kutta methods [10] or implicit schemes such as the midpoint method [8] or backward differentiation formulas (BDF) [9]. Other methods, as semi-analytic and geometric integration and projected Gauss-Seidel, can be found in the review [17] and references therein. Implicit schemes require special treatment of the (non)linear systems of equations, which have to be solved each time step. These systems are basically dense owing to the nonlocal stray field. Efficiency will also strongly rely on the successful application of preconditioners, which might also have to be recomputed during integration [9]. Typical time step lengths reached by implicit second order methods are in the range of picoseconds, while those of explicit higher order schemes lie in the range of several femtoseconds. Hence, implicit schemes will need fewer derivative evaluations for establishing the time step equations, but shift the computational task to the numerical treatment of the (non)linear systems. These systems should be solved accurately and efficient and might require additional field evaluations as well. On the other hand, higher order explicit schemes require larger amounts of field evaluations per step, for instance, the classical $4$th order Runge-Kutta scheme with $5$th order local error estimate requires $6$ evaluations per time step and an $8$th order Dormand-Prince formula with $7$th order local error estimate requires already $13$ evaluations per time step [18]. These explicit methods also incorporate adaptive step size selection, which provides them with additional efficiency and robustness. Equation (1) is non-stiff for largely homogeneous materials and simple geometries [19], but might only get stiffer if grain structures are also modeled [9]. For instance, OOMMF [10], likely the most widely used micromagnetic simulation package, uses explicit (non-stiff) embedded Runge-Kutta formulas of different selectable order for the integration routines of the spatially semi-discretized equation (1). We will construct an explicit higher order scheme for (1) that is especially cheap in terms of stray field evaluations, while maintaining higher order properties for iterates and local error estimates. This is achieved by exploiting extrapolation for the Gragg method [20], also known as explicit midpoint scheme. The meta-principle of (Richardson) extrapolation applies to computed quantities, which depend on a parameter like a mesh or step size. Consider, for instance, a spatially semi-discretized version of (1) and a prescribed initial magnetization. Now, consider the error $\boldsymbol{e}(t;h)$ of a numerical approximation $\boldsymbol{\eta}(t;h)$ of the magnetization $\boldsymbol{M}(t)$ at some time $t$ obtained from an iteration scheme (some ODE solver) that uses a step size $h$. If the error possess an asymptotic expansion in $h$

[TABLE]

we could recompute the approximation with reduced step size, e.g., halved $h/2$, and establish a new extrapolated approximation according to

[TABLE]

For the new approximation the lowest error term is canceled, that is

[TABLE]

This is especially efficient if $\beta>1$, which is true, with $\beta=2$, for symmetric methods [21, 22]. Natural candidates are the midpoint scheme or the trapezoidal rule, which are both implicit and second order in time. Due to the implicit nature, the error expansion of such methods only holds within the numerical accuracy of the solutions of the (non)linear systems. On the other hand, the Gragg method is a symmetric explicit two-step scheme, which is therefore ideal for establishing an exact extrapolation approach for the LLG equation (1). This is done in a triangular Aitken-Neville scheme for polynomial extrapolation, which offers a natural way for adaptive step size and order selection via computationally available local error estimates and the hierarchic order structure. The well-known Gragg-Bulirsch-Stoer (GBS) algorithm [23] for general non-stiff initial value problems is based on the Gragg method and rational function extrapolation. However, it turned out that polynomial extrapolation is almost always more effective [24]. While extrapolation methods for initial value problems are designed for highly accurate nuemrical solutions, the drawback is the increased amount of derivative evaluations because of successive step doubling. In this paper we construct higher order schemes for (1) via polynomial extrapolation of the Gragg method and save expensive stray field evaluations, while simultaneously maintaining the order properties for the iterates and the local error estimates. This is achieved by treating a version of equation (1) with time-linear stray field, where the computational realization of the linear interpolation is incorporated in the extrapolation procedure. We combine the resulting hierarchic structure of higher order schemes in an interplaying step size and order adaptive procedure.

In the following two sections we will clarify the problem setting and give details to the extrapolated Gragg method. Section 4 explains the approach for taming the complexity of the extrapolation scheme. A further section is dedicated to the adaptive step size and order selection. Finally, we validate the method in terms of accuracy and efficiency on variations of the NIST $\mu$MAG Standard problem $\#4$ [25] and also compare it to a higher order Dormand-Prince formula.

2 Problem setting

Let $\Omega\subset\mathbb{R}^{3}$ denote a magnet and $\boldsymbol{m}:\,\Omega\rightarrow\mathbb{R}^{3}$ the reduced (dimensionless) magnetization. The magnetic Gibbs free energy (in dimensionless form) is given by [7]

[TABLE]

that is the sum of exchange-, demagnetizing-, (uniaxial/first order) anisotropy- and external energy, respectively. Here $\mu_{0}$ is the vacuum permeability, $M_{s}$ the saturation magnetization, $A$ the exchange constant, $K_{1}$ the first magnetocrystalline anisotropy constant and $\boldsymbol{a}$ the unit vector parallel to the easy axis. Further, $\boldsymbol{h}_{s}$ is the (dimensionless) stray field, which is defined by the magnetostatic Maxwell equation $-\nabla\cdot\boldsymbol{h}_{s}=\nabla\cdot\boldsymbol{m}\,\,\text{in}\,\,\mathbb{R}^{3}$ and $\boldsymbol{h}_{ext}:=\boldsymbol{H}_{ext}/M_{s}$ is the (dimensionless) external field.

The Landau-Lifschitz-Gilbert (LLG) equation [7, 1, 26] describes the time evolution of the magnetization and is given in a dimensionless and explicit form as

[TABLE]

where $\alpha>0$ is the (dimensionless) damping constant and $\boldsymbol{m}=\boldsymbol{m}(x,t):\,\Omega\times[0,T]\rightarrow\mathbb{R}^{3}$ the time-dependent magnetization. The parameter $\tau$ in equation (6) is dimensionless owing to the relation $\tau=M_{s}\gamma_{0}\,t$ to the physical time $t$, where $\gamma_{0}$ is the gyromagnetic ratio. The effective field $\boldsymbol{h}$ is defined via the functional derivative of the energy

[TABLE]

Equation (6) is supplemented with the initial condition

[TABLE]

and the boundary condition

[TABLE]

where $\boldsymbol{\nu}$ is the outward unit normal on the boundary $\partial{\Omega}$. The LLG equation preserves the magnitude of the initial magnetization as can be seen by scalar multiplication with $\boldsymbol{m}$.

We treat the spatially semi-discretized LLG equation in the spirit of ordinary differential equations as in [8], where $\boldsymbol{m}$ is represented by a discrete mesh vector in each of the $N$ computational units, e.g., cubical/rectangular computational cells. Also, a linear finite element approach on tetrahedral meshes can lead to a similar system of ordinary differential equations [9]. The discrete mesh vectors are collected in a long vector of size $3N$, as well as the discrete effective field components evaluated at the discrete magnetization. The exchange field is discretized by symmetric second order finite differences where the boundary condition (9) is taken into account. The stray field is nonlocal and computed by the algorithm described in [27], which is based on a scalar potential and accelerated by FFT as introduced in [28]. The other components in (7) permit evident discrete local representations. In this sense, we treat equation (6) as a system of ordinary differential equations for the $3N$ components of the long discrete magnetization vector

[TABLE]

where $\boldsymbol{m}=(\boldsymbol{m}_{1}^{T},\boldsymbol{m}_{2}^{T},\ldots,\boldsymbol{m}_{N}^{T})^{T}\in\mathbb{R}^{3N}$.

Numerical algorithms for the LLG equation, which do not preserve the unit norm constraint of the magnetization inherently, consider renormalization of the discrete magnetic spins after each iteration or if some accuracy tolerance is violated. We present our method without renormalization, hence the deviation from the unit norm constraint may also serve as a measure of accuracy. However, there is no limitation to it in the forthcoming method, so that renormalization could be incorporated.

3 The extrapolated Gragg method

Our scheme is based on the explicit two-step midpoint rule (Gragg method) [20, 24] for the initial value problem (10) given in the abbreviated form

[TABLE]

Let the desired approximation to (11) at $t=t_{0}+H,\,H>0$ be denoted with $\boldsymbol{m}_{h}(t)$ where $h:=H/n$ and $n$ an even number. Gragg’s midpoint rule reads

[TABLE]

The method is consistent of order $2$ and the equivalent one-step scheme is symmetric [21]. It therefore possesses an asymptotic error expansion in even powers of $h$, provided the function $F(\boldsymbol{m})$ is sufficiently smooth, that is

[TABLE]

where the expansions are different for even and odd indices $\nu$. There holds $\boldsymbol{a}_{j}(t_{0})=0$ in the even case, but $\boldsymbol{b}_{j}(t_{0})\neq 0$ in the odd case. The existence of an error expansion in powers of $h^{2}$ is crucial for the efficiency of the Richardson extrapolation based method in the forthcoming. Due to Gragg [20] the (first order) explicit Euler starting step is enough for guaranteeing the $h^{2}$-expansion. The averaging in the last line of (12) is a smoothing step, which should originally reduce ’weak stability’ by eliminating the lowest corresponding error term. It is actually not needed for that purpose, if the Gragg method is applied together with extrapolation, which cancels theses error terms anyway. Omitting the smoothing would save one $F$-evaluation and also maintains the asymptotic expansions, where in this case we would simply have $\boldsymbol{m}_{h}(t):=\boldsymbol{m}^{(n)}$. In the algorithm for the LLG equation we will use the smoothing step at no additional cost, since the $F$-evaluation at the interval end is provided for different reason. The well known Gragg-Bulirsch-Stoer (GBS) algorithm [23] for general non-stiff initial value problems (11) is based on Gragg’s midpoint scheme and rational extrapolation. However, it turned out that polynomial extrapolation is more efficient [24]. The extrapolated Gragg method is built up of approximations $\boldsymbol{\mu}_{\ell,1}:=\boldsymbol{m}_{h_{\ell}}(t),\,\,\ell=1,\,2,\ldots$, where $h_{\ell}=H/n_{\ell}$ and $n_{\ell}$ an increasing sequence of even numbers, e.g. the Romberg (power-two) sequence $\{2,4,8,16,32,64,\ldots$}, see Fig. 1.

Owing to the $h^{2}$-expansion (13) the Aitken-Neville algorithm on level $\ell$ leads to

[TABLE]

The $\boldsymbol{\mu}_{\ell,m}$ represent explicit Runge-Kutta (ERK) methods of order $2m$, that is $\boldsymbol{\mu}_{\ell,m}=\boldsymbol{m}(t)+\mathcal{O}(H^{2m})$. Thus, also stability behavior is that of ERK methods. Here, the formula $\boldsymbol{\mu}_{\ell,m-1}$ is the most accurate approximation that can be associated with a computationally available error estimate of order $2m-1$ [22]

[TABLE]

In practice, however, the formula $\boldsymbol{\mu}_{\ell,\ell}$ is taken as the numerical approximation of order $2\ell$ for a prescribed (or determined) level $\ell$ and the error estimate (15) (with $m=\ell$) is used for step size control in the notion of local extrapolation. Note that, due to the different expansions (13) for odd and even indices, the corresponding orders of the error and error estimate for the extrapolated formula would be reduced by one in the case of odd step numbers $n_{\ell}$.

One practically relevant feature of the extrapolated Gragg method is the possibility of adapting the level $\ell$ (and hence the order) in accordance with the step size during computation. We will briefly describe this procedure for the order and step size adapted time integration of the LLG equation in section 5.

The extrapolation via the Aitken-Neville scheme requires, each level $\ell$, the renewed evaluation of the explicit midpoint rule with increased number of steps $n_{\ell}$. The number of right hand side evaluations up to level $\ell$ (with smoothing) is $1+\sum_{\nu=1}^{\ell}n_{\nu}$, where $F(\boldsymbol{m}^{(0)})$ is only computed for $\ell=1$. This complexity is exponentially increasing in the case of the Romberg sequence, which turned out to be most effective in our tests (including step size and order control) compared to different choices, like for instance the harmonic sequences $\{2,4,6,8,10,12,\ldots\}$. To tame this complexity, we will treat the expensive nonlocal stray field differently from the rest of the effective field components.

4 Taming the complexity of the extrapolation

The computational effort for derivative evaluation in numerical schemes for the LLG equation is dominated by the nonlocal part of the effective field. Typically, the computational effort of the stray field computation amounts to $80-90\%$ [14] of the total effective field, depending on the spatial discretization scheme, the numerical scheme for the stray field computation and other local field components (especially concerning the Laplacian for the exchange field). This makes plain extrapolation schemes for the LLG equation inefficient. On the other hand, extrapolation delivers naturally local error estimates and hence the possibility to incorporate adaptive time step selection, which is necessary to make the integration scheme practical and more robust. In addition, the extrapolated Gragg method offers a hierarchy of accurate higher order schemes and the opportunity to adapt the order as well. However, the linearity of the stray field operator with respect to the magnetization offers a way out to tame the complexity. We perform the full explicit midpoint scheme (12) for level $\ell=1$ ($n_{1}=2,\,h_{1}=H/2$), which leads to second order approximations $\boldsymbol{m}^{(1)}=\boldsymbol{m}(H/2)+\mathcal{O}(h_{1}^{2})$ and $\boldsymbol{m}^{(2)}=\boldsymbol{m}(H)+\mathcal{O}(h_{1}^{2})$ according to (13). Let us denote the (linear) discrete stray field operator as $\mathcal{D}$, that is $\boldsymbol{h}_{s}(\boldsymbol{m}):=\mathcal{D}\boldsymbol{m}$. We address the discretized LLG equation in the form ($t_{0}=0$)

[TABLE]

where we emphasize the dependence of the right hand side on the stray field operator $\mathcal{D}$. We now define a piecewise linear stray field on $[0,H/2]\cup[H/2,H]$ and denote the corresponding operator with $\mathcal{D}_{H}$. Hence, there shall hold the interpolation condition at exact (unknown) solution values

[TABLE]

while for intermediate times $t\in[0,H]$ the approximation is second order. Our discretized LLG equation (11) takes now the approximate form

[TABLE]

where $F$ depends on the linearly interpolated stray field instead of the operator $\mathcal{D}$. A computationally realizable approximation of the time-linearized operator $\mathcal{D}_{H}$ is first obtained from the computation at level $1$, where we evaluate the stray field by using the operator $\mathcal{D}$ at the iterates $\boldsymbol{m}^{(0)},\boldsymbol{m}^{(1)}$ and $\boldsymbol{m}^{(2)}$. There holds

[TABLE]

Note that, due to the $h_{1}^{2}$-error expansions of the iterates (13) and the linearity of the operators $\mathcal{D}$ and $\mathcal{D}_{H}$, the above error also involves only even powers of $h_{1}$, that is

[TABLE]

where $\boldsymbol{e}_{j}=\boldsymbol{a}_{j}$ for $\nu=0,2$ with $\boldsymbol{a}_{j}(0)=0$ and $\boldsymbol{e}_{j}=\boldsymbol{b}_{j}$ for $\nu=1$. This means that the interpolation conditions (17) hold approximately ($\nu=1,2$) for the computational realization of $\mathcal{D}_{H}$ with an error involving only powers of $h_{1}^{2}$. This shall make us aware of the possibility of exploiting efficient extrapolation for the values $\mathcal{D}\boldsymbol{m}^{(\nu)},\,\nu=1,2$ to establish more accurate approximations to $\mathcal{D}_{H}\boldsymbol{m}^{(\nu)},\,\nu=1,2$. As step sequence we choose the Romberg sequence, hence the amount of steps is doubled and the step size halved from one to the next extrapolation level, compare with Fig. 1. Each level we perform the Gragg method (12) with smoothing, where we save stray field evaluations by using interpolated values from the current approximate version of $\mathcal{D}_{H}$, i.e. we are solving (18). Renewed evaluations of the stray field are only necessary at $t=H/2$ and $t=H$ followed by the computation of a new line in the Aitken-Neville scheme (14) for both, the current iterate $\boldsymbol{m}_{h_{\ell}}(H)$ and the stray field values $\mathcal{D}\boldsymbol{m}^{(n_{\ell}/2)}$ and $\mathcal{D}\boldsymbol{m}^{(n_{\ell})}$. Note that the midpoint in level $1$ has odd parity, while for all subsequent levels it is even. We therefore use for $\ell=1$ a centered average of the stray field at the midpoint. This is

[TABLE]

where the error term involves only even powers of $h_{1}$. Note that no further evaluations of $\mathcal{D}$ are needed here, since

[TABLE]

where the values on the right hand side are already available.

5 Step size and order control

According to [24, 22] we take for level $\ell$ the order $2\ell-1$ error estimate (15) of $\boldsymbol{\mu}_{\ell,\ell-1}$ for the numerical approximation $\boldsymbol{\mu}_{\ell,\ell}$ in the notion of local extrapolation. We remark, that also the error estimates of the extrapolation of the stray field values are available and can be incorporated in several different ways. One possibility is to simply establish a weighted sum of relative error estimates. As usual, we require the dominant term $err_{\ell}\approx C\,H^{2\ell-1}$ in the error estimate for a given basic step size $H$ to reach a tolerance $tol\approx C\widetilde{H}^{2\ell-1}$ obtained from an adapted step size $\widetilde{H}$. This, together with incorporated safety factors, yields the empirically optimal choice [18] for an adapted step size $H_{\ell}$ at level $\ell$

[TABLE]

Equation (23) is used for determining a next step size within a convergence monitor for the three subsequent levels $\ell,\ell+1$ and $\ell+2$, which determine whether the current approximation is accepted or rejected and the order increased or decreased [18]. The tool for measuring the necessity and efficiency for order and step size adaption is the reduction of work per time step size $W_{\ell}/H_{\ell}$, where the work $W_{\ell}$ measures the effort for computing the numerical approximation $\boldsymbol{\mu}_{\ell,\ell}$. An adapted choice for the step size and order shall reduce the work per time step size. We define $W_{\ell}$ as the weighted sum of stray field evaluations and effective field evaluations up to level $\ell$. The latter one uses already computed stray field evaluations and, hence, can be understood as the amount of evaluations of all other field components except the stray field. As a weighting factor we take $f_{sf}=0.8-0.9$, which shall be a rough estimate of the portion of the costs for the stray field compared to the total field. Hence, the work is defined as

[TABLE]

where $n_{\nu}=2^{\nu}$ for the Romberg sequence. Note that the stray field part of the work increases linearly with the level, while the other part increases exponentially.

6 Numerics

We look at the NIST $\mu$MAG Standard problem $\#4$ [25]. The geometry is a magnetic plate of size $500\times 125\times 3$ nm with material parameters of permalloy: $A=1.3\times 10^{-11}$ J/m, $M_{s}=8.0\times 10^{5}$ A/m, $\alpha=0.02$. The initial state is an equilibrium s-state, obtained after applying and slowly reducing a saturating field along the diagonal direction $[1,1,1]$ to zero. Then an external field of magnitude $25$mA is applied with an angle of $170^{\circ}$ c.c.w. from the positive $x$ axis. We use different spatial discretizations, where the finest is built of $1$nm cubes. Errors are computed in the relative Euclidean norm, where we weighted error estimates from the stray field extrapolation marginally with two percent. However, investigation of the decrease of the error estimates of the stray field values at the mid- and endpoint, $H/2$ and $H$ respectively, showed analogue decay rate and magnitudes as for the error estimates of the magnetization iterates at the endpoint. Simulations are performed by subdividing the time interval into $1$ps subintervals and data were captured at every $1$ps simulated time. Some measures like the time step sizes, the extrapolation level or the numerical damping parameter were recorded within the $1$ps subintervals and archived as averaged values. Fig. 2 shows the time evolution of the averaged magnetization components for $1$ns simulated time and the $1$nm cube discretization. Computations in Fig. 2 were performed with $f_{sf}=0.85$ and a tolerance of $tol=1.0\times 10^{-12}$.

In Tab. 1 we give statistics of these computations for tolerances of $1.0\times 10^{-10}$ and $1.0\times 10^{-12}$ including the spent work w.r.t. $f_{sf}=0.8,0.85$ and $0.9$, number of function evaluations, average extrapolation level and step size, number of rejected steps and absolute maximum error of unit norm constraint (LLG preserves the magnitude of the moments). No renormalization is performed.

We also compared our method with the Dopri8 (Dormand-Prince) method from [18], which is an $8$th order embedded Runge-Kutta method using a $7$th order estimate for the error, results are also included in Tab. 1. Averaged magnetization components coincide with our method with absolute error in the range of $1.0\times 10^{-6}-1.0\times 10^{-4}$. At time $t=0.138$ns, approximately the moment where $\langle m_{x}\rangle$ crosses zero for the first time, the discrete magnetization configurations were captured for the Runge-Kutta and our method, giving a calculated maximum absolute deviation on the entire $1$nm grid of about $1.0\times 10^{-6}$.

We consider one more error measure: the relative error of the numerical damping parameter $|\langle\alpha\rangle-\alpha|/\alpha$ at a time $t_{\nu}$ according to [29]

[TABLE]

where $j$ is the node index and we approximate $\Delta\varepsilon_{j}^{(\nu)}=\varepsilon^{(\nu)}_{j}-\varepsilon^{(\nu-1)}_{j}\approx-\boldsymbol{h}^{(\nu-1/2)}_{j}\cdot\Delta\boldsymbol{m}_{j}^{(\nu)}=-\boldsymbol{h}^{(\nu-1/2)}_{j}\cdot(\boldsymbol{m}^{(\nu)}_{j}-\boldsymbol{m}^{(\nu-1)}_{j})$ and use for the field at the midpoint $\boldsymbol{h}^{(\nu-1/2)}_{j}\approx(\boldsymbol{h}^{(\nu)}_{j}+\boldsymbol{h}^{(\nu-1)}_{j})/2$. Note, however, that this is itself a first order approximation to the analytical expression

[TABLE]

where $\varepsilon$ is the energy density. The errors are plotted in Fig. 3 associated with the computations of Fig. 2, where the average step size was about $195$fs. The computations in Tab. 1 are repeated on a coarser mesh consisting of about $2\times 2\times 1$nm prisms (mesh size $250\times 64\times 3$), see Fig. 4 and Tab. 2.

As a second numerical test we change parameters of the original setting of Standard problem $\#4$. Results in Fig.5 were obtained by changing the damping parameter to $\alpha=0.2$ and shows propagation of averaged magnetization with $1$nm discretization, the averaged time steps (within the $1$ps subintervals), the average extrapolation levels and work per (reduced) time step. One can recognize the interplay between order and step size adaption, while the work per time step remains roughly unchanged.

Now we change the anisotropy constant to $K_{1}=1.0\times 10^{6}$ J/m3 with the easy axis $\boldsymbol{a}=[1,1,1]/\sqrt{3}$, while maintaining all other original parameters. Fig. 6 shows computation results with $1$nm discretization for the propagation of the averaged $y$-component (the others oscillate similarly) and the deviation from the corresponding values obtained from the Runge-Kutta method. In Fig. 7 we give averaged time steps (within the $1$ps subintervals), the average extrapolation levels and work per (reduced) time step, associated with Fig. 6.

From the test examples one can recognize an advantage of the extrapolation method in terms of required work. This aspect gets more significant when the true $f_{sf}$ gets larger, that is, stray field computation increasingly dominates other computational costs. In all tests the average step sizes are clearly larger and stray field computations are fewer. Moreover, it is noticeable that due to the simultaneous order and step size control there are actually no rejected (wasted) steps.

7 Conclusions

We developed a step size and order adaptive solver for the Landau-Lifschitz-Gilbert equation. The method uses extrapolation of the symmetric explicit midpoint scheme, which possess an asymptotic error expansion in even powers of the step size parameter. The necessary number of expensive stray field evaluations is reduced to linear dependence on the order of the method. This is achieved by a piecewise time-linear stray field approximation. We show how to efficiently extrapolate this approximation by utilizing the $h^{2}$-expansion and the linearity of the stray field operator. Numerical experiments indicate that the proposed scheme gets more and more efficient, compared to conventional methods as higher order Runge-Kutta, when the stray field computation takes a larger portion of the costs for the effective field evaluation. This is more likely the case in field-based stray field approaches.

Acknowledgments

Financial support by the Austrian Science Foundation (FWF) under grant No F41 (SFB ’VICOM’), grant No F65 (SFB ’Complexity in PDEs’) and grant No W1245 (DK ’Nonlinear PDEs’) and the Wiener Wissenschafts- und TechnologieFonds (WWTF) project No MA16-066 (’SEQUEX’). The computational results have been achieved using the Vienna Scientific Cluster (VSC).