An Adaptive Step Implicit Midpoint Rule for the Time Integration of Newton’s Linearisations of Non-Linear Problems with Applications in Micromagnetics

The implicit mid-point rule is a Runge–Kutta numerical integrator for the solution of initial value problems, which possesses important properties that are relevant in micromagnetic simulations based on the Landau–Lifshitz–Gilbert equation, because it conserves the magnetization length and accurately reproduces the energy balance (i.e. preserves the geometric properties of the solution). We present an adaptive step size version of the integrator by introducing a suitable local truncation error estimator in the context of a predictor-corrector scheme. We demonstrate on a number of relevant examples that the selected step sizes in the adaptive algorithm are comparable to the widely used adaptive second-order integrators, such as the backward differentiation formula (BDF2) and the trapezoidal rule. The proposed algorithm is suitable for a wider class of non-linear problems, which are linearised by Newton’s method and require the preservation of geometric properties in the numerical solution.

with explicit and implicit methods. The design and development of these algorithms has been the topic of much previous research (see, for example [3,Sect. 3.16], [4][5][6][7]). Many publicly available codes for the solution of IVPs (for example CVODE [8], or the various solvers in Matlab [9]) involve some adaptivity strategy. Adaptive IVP solvers allow the selection of the right step sizes at each stage of the solution process to satisfy the required accuracy, thus removing the need for a "trial-and-error" approach. This also improves the computational efficiency and robustness of the overall algorithm.
In some problems the governing equations dictate that the solution vectors should have constant magnitude, and that the system's energy should either be conserved or accurately reproduced. In such cases IVP solvers should preserve these properties at the discrete level, while computational efficiency is achieved by using the adaptive implicit methods. As an example, we consider micromagnetics [10]. Dynamical micromagnetic models are based on the Landau-Lifshitz-Gilbert (LLG) equation, a time-dependent differential equation modelling the magnetization evolution in a ferromagnetic body. The classical LLG equation models the magnetization at zero temperature, which implies that the magnetization vector length is constant at any point of the domain. The second important property of the LLG equation is the energy conservation in the non-damped case, or equivalently, under a constant applied field the energy decreases at a rate proportional to the damping factor. Solutions of the LLG equation in practical cases often exhibit multiple spatio-temporal scales favouring the use of adaptive implicit time integration schemes [11]. An integration scheme that preserves (or, in the damped case, accurately reproduces) the length/energy properties of the solution is referred to as a geometric integrator. Most commonly used implicit schemes (such as the second-order backward differentiation formula BDF2 or the trapezoidal rule (TR)) do not possess this property. By contrast, the implicit mid-point rule (IMR), combined with an appropriate spatial discretisation in the PDE cases, is a geometric integrator [12,13]. The proposed adaptive version of IMR retains this favourable property. Some alternative numerical techniques aimed at preserving the magnetization length include peiodic renornalisation [14], length constraints imposed through Lagrange multipliers [15], penalty formulations [16], and the self-correcting LLG equation used by Nmag [16]. However, all these formulations change the energy of the system.
The implicit mid-point rule (IMR) is a well-known second order implicit Runge-Kutta method [17, p. 204], [3, p. 262], that has been applied in fixed-step form to micromagnetic problems [12]. In this paper we present an adaptive version of the IMR based on the new local truncation error estimator in the context of a predictor-corrector scheme. We demonstrate the efficiency of the adaptive version of IMR by comparing its accuracy and computational cost with the fixed step version of the method [12], and with the adaptive versions of the BDF2 [3, p. 649] and the TR [3, p. 647] when applied to both ODE and PDE problems that arise in computational micromagetics.
The paper is organised as follows. In Sect. 2 we introduce the LLG equation and describe a method for its discretisation that is used in this paper. Section 3 introduces the IMR and details some of its properties. In Sect. 4 we present the new local truncation error (LTE) estimator for the IMR and introduce the adaptive integrator. Finally, in Sect. 5 we evaluate the effectiveness of the adaptive IMR when applied to ODE and PDE problems associated with the LLG equation. We also compare the adaptive IMR to adaptive BDF2 and the standard version of TR in terms of the computational cost, accuracy and the preservation of geometric properties.

The Landau-Lifshitz-Gilbert Equation
Micromagnetics is a semi-continuum theory of ferromagnetic materials that covers the lengthscales between the atomistic spin dynamics and the classical Maxwell's theory. Its aim is to describe the behaviour of a magnetic body, expressed by its magnetization vector m, under various external conditions and internal material properties. This can be achieved either by calculating the minimum energy of a magnetic body [10], or by modelling the magnetization dynamics. In the latter case the time evolution of the magnetization is modelled by a differential equation. Given a space of k-continuous functions C k (·), the evolution of the magnetization vector denoted by m ∈ [C 2 ( )] 3 × C 1 (T ) in a magnetic body ⊂ R d , d = 2, 3 with a smooth boundary ∂ over a time interval T follows the Landau-Lifshitz (LL) equation, which, in non-dimensional form, reads and the initial condition is the effective magnetic field, and α is an empirical damping parameter. In (2)n denotes the outward normal vector on ∂ . The first term on the right-hand side of (1) determines the precession of m, while the second term is responsible for damping. A mathematically equivalent form of (1), referred to as the Landau-Lifshitz-Gilbert (LLG) equation is frequently used, which in non-dimensional form reads The effective field h in (1) and (4) where h ap is the applied (external) field, h ex = ∇ 2 m is the is a phenomenological exchange field that represents the effect of quantum mechanical exchange, which acts to align the magnetic moments of neighbouring atoms in ferromagnetic materials, h ms = ∇φ is the magnetostatic (demagnetisation or stray) field expressed in terms of the magnetostatic potential φ, and h mc = k 1 (m ·ê)ê is the magnetocrystaline anisotropy field associated with the easy magnetization axisê of a ferromagnetic material [18, p. 306]. The magnetostatic potential φ(x), x ∈ R 3 is related to the magnetization vector through the equation where * is the exterior of a ferromagnetic body . Boundary conditions associated with (6) are of the form on ∂ , where φ int and φ ext are the values of the potential just inside and outside of , respectively, and φ(x) → 0 for |x| → ∞.
The fact that (6) is posed over an infinite domain, together with the BC (9) at infinity cause difficulties with the discretisation of this part of the problem.

Finite Element Discretisation
The weak formulation of (4) in a residual form is to find where H 1 ( ) is the standard Sobolev space. The discrete problem is obtained by restricting this space to a finite element space S h associated with the subdivision of into a set of quadrilateral elements and adopting piecewise bi/trilinear polynomial approximation m h for the solution m over the elements, i.e. S h = span{ψ k } N k=1 , N = N i + N b (N i is the number of nodes in the interior of and N b the number of nodes on ∂ ). Then we have with m k = [m k,x m k,y m k,z ] t representing the values of m at node k. This process results in a semi-discrete problem, which is fully discretised approximating the time derivatives in (10). We use the IMR method in this context (see Sect. 3). The linearisation of the non-linear discrete problem at each time step is done using the Newton method, leading to the solution of a sequence of linear algebraic systems of the form J δ m = −r ,m +1 =m +δ m , = 0, 1, . . .. The Jacobian matrix with elements J mk = ∂r m ∂m k can be written in block form whereM is a matrix similar to a mass matrix [19, p. 55], while the blocks K * are sums of mass and Laplacian-like matrices (for detailed expressions see [13, pp. 97-102]). The Jacobian matrices are assembled from the elemental contributions with their coefficients calculated two ways: either using standard Gaussian quadratures [19, p. 30], or nodal (Newton-Cotes) quadratures [20]. The former approach is standard in FE computations and gives optimal order of accuracy per amount of computational work. However, it does not lead to the conservation of the magnetization length. 1 The nodal quadrature (also known in micromagnetic literature as the "reduced integration" [21]) for an integral over an element e is defined by where ϕ (e) is the local basis function associated with the local node and V ( e ) is the set of nodes of an element e . The nodal quadratures use element's nodes as the knots, thus allowing the point-wise magnetization length conservation when the weak formulation of the problem is used as the basis for spatial discretisation [20]. In order for this procedure to be effective, we need to set a tight tolerance N in the stopping criterion of Newton's method ( e.g. N < 10 −10 ). Notice that J in (12) is a block skew-symmetric matrix. Linear systems with the coefficient matrix J are solved either directly, or by preconditioned Krylov solvers. In our experiments reported in Sects. 5.2 and 5.3 GMRES preconditioned by ILU(1) method is used [22].

Hybrid Finite-Boundary Element Method for the Exchange Field
Spatial discretisation of (6) with BCS (7)-(9) is problematic, as the domain * is infinite. The application of FEM in this context would require so-called "infinite elements" [23]. An alternative is to truncate the domain * or to use asymptotic BCs [24]. The asymptotic accuracy of such approaches is usually low. An alternative approach is to apply the hybrid finite/boundary element method, in which the external domain is replaced by a dipole layer placed on ∂ , which simulates the correct BC (9) [25]. Standard FEM is used to calculate the potentials in .
We use the linearity of the problem (6)-(9) to introduce two new potentials u(x) and v(x), such that φ = u + v, and u is defined to be non-zero only in . Therefore, u(x) satisfies Thus, u is calculated in the domain using the FEM and the double layer charge distribution on ∂ simulates the potential v in * and is solved by BEM. After setting the conditions for v [13, p. 114-116], we obtain the following problem: subject to the BCs where G[ · ] is defined as In (20) G(x, y) denotes the Green's function of the Laplacian, and γ (x) is the fractional (space) angle, where for any smooth part of ∂ we have γ = 1 2 . After discretisation, Eq. (18) becomes a matrix-vector productφ = Gū, where the elements of the BEM matrix G are obtained as In (21) ψ j is the basis function, d = 2, 3 is the spatial dimension, and δ ji is the Kronecker delta.
The accuracy of the FEM/BEM method is better than that of the asymptotic BCs method [26]. The main drawback of the hybrid method is that it produces a dense coefficient matrix of size N b , leading to both the computational and storage cost of O(N 2 b ) (in the discrete version of (18) the BEM matrix needs to be multiplied by a vector to impose the BCs). This cost can be prohibitive when modelling micromagnetic problems posed over thin domains, where N b = O(N ). This problem was circumvented by the approximation of the BEM matrix by low rank hierarchical matrices, as implemented in the library HLib [27,28], reducing both the storage and computational cost from

The Implicit Mid-Point Rule
Consider a system of N coupled IVPs y = f(t, y(t)), for t ∈ T = [0, t max ], where y(t) ∈ [C 1 (T )] N is a vector function, y 0 is the initial condition, and f is a locally Lipshitz-continuous function. We denote an approximation to y(t n ) by y n and write n = t n+1 − t n for the size of n-th time step, and introduce three commonly used second-order implicit methods for numerical integration of (22). The IMR is defined by [17, p. 204] With the notation t n+1/2 = t n+1 +t n 2 = t n + n 2 , y n+1/2 = y n+1 +y n 2 for the mid-point approximations of t and y, (23) becomes y n+1 = y n + n f(t n+1/2 , y n+1/2 ).

Local Truncation Errors and the Order of Convergence
The local truncation error (LTE) of an integration scheme is the error due to a single integration step [1, p. 40]: whereŷ n+1 is an approximation to y(t n+1 ) obtained by the integration method with all the history values y n , y n−1 , . . . being exact (i.e. y n−i = y(t n−i )). The LTE expression for the IMR is more complex than those of TR or BDF2 due to the approximation y ≈ y n +y n+1 and Substituting (28) and (29) into gives The LTE expression (31) has two terms: terms similar to (I) feature in LTEs of many other second-order integrators [cf. (32) and (33) and

Properties of the IMR
The IMR has two important properties that makes it a method of choice for the time integration of semi-discretised PDEs: • It is A-stable [17, p. 43,251], i.e. it can be applied to stiff problems.
• It requires only one evaluation of the function f per step (a one-step method, unlike many other Runge-Kutta methods).
These properties of the IMR method are shared with TR and BDF2. However, the IMR does not introduce spurious damping of oscillatory solutions (TR does possess this property, although it is lost in the modified versions of it [4,30], introduced to avoid the locking effect in the step size selection caused by the accumulation of the round-off errors). All backward differentiation formulae are damping [3, p. 265]. When applied to the time integration of the LLG equation, the IMR has geometric properties, related to the conservation of the magnetization length and reproduces a correct energy behaviour of the system. It is also worth noticing that when applied to the integration of linear problems, IMR and standard TR will produce exactly the same numerical solutions, i.e. the two methods are equivalent (by virtue of

Adaptive Implicit Mid-Point Rule
Implicit time integration of ODEs or semi-discretised PDEs is a computationally intensive task. Numerically effective integrators employ adaptivity, which enables an integrator to take steps of the size governed, within the prescribed tolerance, by the physics of the problem, rather than artificial (method-specific) constraints. Adaptive step selection algorithms rely on computable estimates of the LTE. In this section we discuss the difficulties in applying the commonly used LTE estimation schemes to IMR and introduce the new approach.

Local Truncation Error Estimation
Typical LTE estimators used in implicit IVP integrators (such as TR or BDF2) are based on the solution estimate y E n+1 ∼ y(t n+1 ) computed by an explicit method with the same asymptotic order of accuracy, using the history values computed by an implicit method. This is referred to as the predictor step. Algebraic rearrangements of the analytical expressions for the LTEs of an explicit and an implicit method lead to a computable estimate of the implicit LTE with a desired order of accuracy. This approach is known as Milne's device [3, p. 707-716]. It is not straightforward to apply this approach to estimating the LTE of the IMR, due to the presence of the term II in (31). The difficulty of constructing a Milne device in this case lies in finding a suitable predictor which LTE involves the term II from (31). In the case of linear or linearised problems (such as the Simo-Armero linearisation of the Navier-Stokes equations [31]), the IMR is numerically equivalent to TR, which implies that the use of the standard AB2 predictor will produce the same behaviour (in terms of the number and the size of the selected steps) for both AB2-TR and AB2-IMR integrators. This is no longer the case when f is non-linear and we focus on this scenario in the remainder of this paper.
Adaptivity in explicit Runge-Kutta methods is achieved by finding a pair of methods of different orders and to obtain a LTE estimate by comparing the two solutions [17, p. 165]. The effectiveness of this approach relies on the existence of pairs of methods that share the most of their function evaluation points -the so-called embedded methods (the examples include the Dormand-Prince pair of order 4/5 deployed in Matlab's function ode45 [9]). To deploy this approach in the context of the IMR would require a third order implicit method that involve a function evaluation in addition to f(t n+1/2 , y n+1/2 ). However, the function evaluation f(t n+1/2 , y n+1/2 ) at the "mid-point" of the interval [t n , t n+1 ] leads to the cancellation of the second-order terms in (28)- (29), and an additional function evaluation at a point in [t n , t n+1 ] would break this symmetry. To restore it, we would require two additional function evaluations.
In order to bypass these problems, we consider explicit third order methods as predictors in the non-linear case. Suitable candidates are AB3 [1, p. 127], RK3 [1, p. 73], and the eBDF3 method [17, p. 378]. The first two of these methods are computationally more expensive and/or require the storage of more history values than the AB2 method. This is, however, not the case with the eBDF3 method. In order to be a competitive predictor, it should lead either to the overall reduction in the number of time steps in non-linear cases for a set LTE tolerance, or a better accuracy for a fixed number of steps performed compared to the AB2 predictor. In particular, the AB2 method requires the storage of the function values (the solution derivatives) at two history points. In the non-linear case this ammounts to either storing the Jacobians at these points or assembling them on the fly. We notice that the function evaluations at multiple points also introduce additional implementation complications in ODE/PDE cases where the time derivative is given implicitly (such as in the LLG equation). The situation is even worse for AB3 and RK3 methods that require the computation or storage of the function values at three history points. Neither of these predictors is self-starting, which is a drawback for a self-starting IMR method. The eBDF3 method is not self-starting either, but requires the storage of three history solution values, rather than its derivatives, and only one function evaluation (at n ). This makes it a competitive alternative to the AB2 method and the other two third order methods. The computational cost of applying the eBDF3 predictor consists of one sparse matrix-vector multiplication (where we assume that a dense matrix block obtained in the BEM discretisation of the magnetostatic field is sparsified prior to its use) and, in PDE cases, the solution of one trivial linear system with a diagonal mass matrix (assuming that reduced quadratures are used). The eBDF3 predictor is therefore judiciously chosen among the other alternatives (AB2, AB3, RK3) due to its low computational and/or storage overhead.
The eBDF3 method is unstable for fixed step sizes [17, p. 378], however this is not an issue when used as a predictor, as we only perform a single step of it using IMR-computed history values. Thus, the eBDF3 method can be viewed as an extrapolation of y at t n+1 using the history values.
The LTE estimate of a generic explicit third-order method at the time step n = t n+1 − t n is given by: and the LTE of the IMR can be written as Expressing y(t n+1 ) from (34) and substituting to (35) gives i.e. we have an estimate:

The Variable Step eBDF3 Method
Both explicit and implicit BDF methods are derived from a divided difference representation of an interpolation polynomial p(t) through y n+1 and a number of history values. Setting Differentiating (38) yields Equation (39) can be simplified by noticing that the product on the right-hand side evaluated at t = t n (for i = 1) is zero. This gives For k = 1 and k = 2 (40) reduces to the well-known forward Euler and the leapfrog (explicit mid-point rule) methods, respectively [3, p. 715]. For larger values of k the complexity of (40) increases significantly. The eBDF3 method is obtained for k = 3: The rearrangement of (41) for y n+1 was done using the symbolic library SymPy [32,33] resulting in y n+1 = bf(t n , y n ) + c n y n + c n−1 y n−1 + c n−2 y n−2 , where b = n+1 n ( n + n−1 ) For constant step size ( n = n−1 = n+1 = ) (43) reduces to the scheme given in [17, p. 364 The computational cost of (43) can be high for linear scalar IVPs. However, if non-linear IVP systems are solved in the context of the method of lines, the cost of computing (43) is small compared to the (repeated) solution of a linear system required by an implicit corrector.
The LTE estimate for IMR (37) is used to predict the next step size via the standard heuristics that involves the user-prescribed tolerance ε T [3, p. 268]: where · is a suitable norm of the LTE (e.g. a mass matrix norm, see [3, p. 708]). In the example from Sect. 5.1 we do not introduce any further restrictions in the step selection while in the cases from Sects. 5.2 and 5.3 we make the criterion (44) more robust by adding the following two heuristical rules, which are in line with the recommendations from [19, p. 432]: i) if n+1 / n is smaller than 0.7, we adopt n = n /2 and repeat the calculation for the step n (the step cancellation), and ii) we limit the step size growth factor n+1 / n by 4 to improve numerical stability. We notice that some public IVP solvers, such as CVODE [8, Sect. 2.1] keep the step size constant as long as feasible. This allows the use of the inexact Newton's method, which does not require the assembly of the Jacobians at every iteration, thus reducing the computational cost. However, if the associated linear systems are solved by a preconditioned Krylov method, the gains in convergence speed obtained using exact Newton's method can outweigh the savings obtained by not assembling the Jacobian matrices [34, p. 128]. The adaptive highorder explicit method from [7] also selects fairly constant time steps over long time intervals.

Numerical Examples
In this section we present numerical experiments that demonstrate the effectiveness, accuracy, and geometric properties of the proposed adaptive IMR scheme. The first case study is an ODE system modelling a magnetization reversal in a small spherical (both isotropic and anisotropic) ferromagnetic particle, while the two further cases are the PDE problems which involve spatial variation of the magnetization vector and are discretised using finite element method and the method of lines is applied to the semi-discretised system of IVPs [35]. The implementation of the ODE example in Sect. 5.1 is done in Matlab, while the PDE examples are implemented in oomph-lib, an object-oriented multi-physics finite element library [36].

The Magnetic Reversal of a Spherical Particle
This is a widely used model-problem from micromagnetics which is represented by an ODE system. We consider the magnetization of a small spherical particle made of either an isotropic or anisotropic ferromagnetic material under a spatially uniform applied field h ap . The evolution of the magnetization vector m ∈ C 1 (T ) follows the LL Eq. (1) with h = h ap + k 1 (m · e)e is the effective field. In case of an isotropic material (k 1 = 0) h has only the applied field component [18, p. 306]. For a sufficiently small sphere the magnetization is spatially uniform, i.e. m = m(t) and the problem (1) is a system of three ODEs for the unknown Cartesian components of m. The problem has an analytical solution [37], which describes the process of the magnetization switching and can be expressed in spherical polar coordinates (θ, φ) as θ = cos −1 (m z /1), φ = tan −1 (m y /m x ), where θ is the angle between m and h ap . The energy of the system is in this case given by E = −m · h ap . In the ideal non-damped case (α = 0) the energy remains conserved throughout the magnetization reversal (i.e. the process of the realignment of m from the initial configuration m 0 to the external field h ap ) -the property that needs to be replicated by the time integrator. For anisotropic materials (k 1 > 0) we study the behaviour of the time integrators for the realistic values of the phenomenological anisotropy parameter k 1 ∈ [0, 4]. a) Isotropic case (k 1 = 0). We consider the damping value α = 0.01. In this case the problem (1) is only mildly non-linear. Thus, we would expect that the AB2-IMR method performs well, and its step size selection to be comparable to that obtained with a thirdorder predictor. We choose h ap = (0, 0, −1.1) t , which is a sufficient magnetic field strength  Table 1 The number of the time steps N t and the minimum of magnetization length |m| min for the spherical particle reversal problem with α = 0.01 and k 1 = 0 integrated over the time interval T = [0, 1000] using GL-BDF2, AB2-TR, AB2-IMR and eBDF3-IMR schemes to reverse the initial magnetization m 0 =m 0 / m 0 , wherem 0 = (0.01, 0, 1) t . In all our experiments it took 1-2 iterations for the exact Newton method to converge to a prescribed tolerance ε N = 10 −14 . The system (1) is integrated over the interval T = [0, 1000], which is sufficient to reach a steady state after the reversal. The time evolution of the Cartesian components of the magnetization vector m = (m x , m y , m z ) t is presented in Fig. 1. We first examine the efficiency of the step size selection procedure and the magnetization length conservation properties of different integration schemes as a function of the LTE tolerance T . The results are summarised in Table 1. From these figures we see that the adaptive step IMR method conserves |m| to the level of Newton's tolerance ε N , while this is not the case for the BDF2 and TR methods. TR is more accurate in this context than BDF2, which can also be seen from Fig. 2, which depicts the time evolution of the magnetization length for the three time integrators. Comparing the number of time steps required to complete the integration, it follows that AB2-TR and AB2-IMR have the same behaviour (which is expected in the case of linear and mildly non-linear problems). The combination eBDF3-IMR leads to approximately 40% higher number of time steps than the AB2 predictor for a fixed value of T . The time step evolution for various integrators is presented in Fig. 3. Next we test the accuracy of the proposed time integrators. To this end, we record the time t 1 for which m z (t 1 ) = 0 (cf. Fig. 1), i.e. the time t 1 for which m·k = 0. The reference solution for this problem is calculated by both adaptive TR and IMR with a tight tolerance ( T = 10 −10 ) as t s = 481.72. In Table 2 we present the accuracy results when all the tested methods take the same number of steps to perform the simulation over the interval T = [0, 490]. This is achieved by adjusting the LTE tolerance T for each method separately. The rationale behind this experiment is to evaluate the quality of the integration points distribution for each of the predictors considered.
From these results we see that the combination eBDF3-IMR is as accurate as AB2-TR. This is due to the character of the problem, which is only weakly non-linear, implying a near Table 2 The switching times t 1 for which m z (t 1 ) = 0 with fixed number of time steps N t for the spherical particle reversal problem with α = 0.01 and k 1 = 0  0). This case involves material anisotropy, expressed by an additional field component (the anisotropic field). We test the efficiency and accuracy of the proposed time integrators when the non-linearity in the problem, controlled by the parameter k 1 , increases. In particular, there is no equivalence between the TR and the IMR method in these cases, and we expect that the third order predictor for IMR leads to a better accuracy when the number of time steps is fixed, when compared to AB2-TR and AB2-IMR methods. As in the previous example, we choose h ap = (0, 0, −1.1) t , m 0 =m 0 / m 0 , wherem 0 = (0.01, 0, 1) t , and select the direction of the material anisotropy asê = e/ e where e = (1, −0.3, 0) t . The remaining problem and algorithm parameters are the same as in the isotropic case. The size of the time interval required to complete the reversal varies for different values of k 1 and is selected to be sufficiently long for the system to reach a steady state. For any non-zero value of k 1 the steady-state position of the magnetization vector after the reversal is no longer aligned with h ap . The time evolution of the magnetization vector components with k 1 = 4 is presented in Fig. 4.
In Table 3 we present the numbers of time steps N required to complete the integration for the LTE tolerance T = 10 −5 and the bounds for the magnetization length |m| for different values of k 1 and different integrators. From the table we see that the increase in k 1 leads to the poorer magnetization length preservation (although the used level of toletance T gives fairly good values for all integrators, looser values of T will cause problems for BDF2 and TR, with notable loss of accuracy in |m|, see Fig. 5 for the case T = 10 −4 ). As k 1 is increased, the combination AB2-IMR requires fewer steps than AB2-TR.
The time step size history for three different combinations of time integrators is presented in Fig. 6. Comparing to the isotropic case, where a smooth variation in step sizes is observed in all time integrators, we see that the step size selection exhibits a considerably more dynamic behaviour, which is the consequence of a more complex dynamics in the problem.
Next, we report the accuracy of the first zero crossing in the z component of the magnetization for k 1 = 4 (i.e. the smallest time t 1 for which m z (t 1 ) = 0). We do this with respect to two different criteria: fixed number of time steps and fixed execution time for all tested integrators. The first criterion reflects the quality of a predictor, i.e. how the distribution of a prescribed number of the integration points that it produces affects the solution accuracy. The second criterion tests whether the computational overhead associated with adaptive step  Table 3 The number of the time steps N t and the magnetization length interval [|m| min , |m| max ] for (1) for the anisotropic spherical particle reversal problem with α = 0.01 as a function of k 1 with T = 10 −5 . The integration is done over the interval T = [0, t k 1 ]   Fig. 6 The time step size evolution for GL-BDF2, AB2-TR and eBDF3-IMR methods for the anisotropic spherical particle reversal problem with α = 0.01, k 1 = 4 and T = 10 −5 selection leads to a more accurate method than computationally cheaper fixed step counterpart within a set amount of wall clock time.
For the first criterion we consider three different levels of precision, determined by a fixed number of steps N 1 required to cover the interval [0, 150] for all methods. These levels are N 1 = 4142, 8967, and 19,336, obtained by setting T in AB2-TR method to 10 −5 , 10 −6 , and 10 −7 , respectively, and adjusting it in other integrators to match this number of steps. The reference solution t ref 1 = 145.038 is obtained by the eBDF3-IMR and AB2-TR methods with T = 10 −10 . The results are presented in Table 4. From the results we observe that Table 4 The computed minimum values of t 1 for which m z (t 1 ) = 0 for the anisotropic spherical particle reversal problem with α = 0.01, k 1 = 4 as a function of the number of time steps N 1  Table 5 The computed minimum values of t 1 for which m z (t 1 ) = 0 for the anisotropic spherical particle reversal problem with α = 0. the eBDF3-IMR combination is the most accurate, although for very tight tolerances T the differences to AB2-TR and AB2-IMR methods become less pronounced. For the second criterion (fixed execution time) we compare only the adaptive and the fixed step versions of the IMR integrator. We consider again three different levels of precision for adaptive IMR and adjust the number of time steps for the fixed step variant so that both methods have the same execution time to perform the integration over the time interval [0, 500]. In Table 5 we report the values of t 1 at which the first zero crossing m z (t 1 ) occurs. From these figures it can be concluded that fixed step IMR completes more steps in a set amount of time (due to an increased computational cost per step in the adaptive method), but still has lower accuracy, justifying the additional costs of the adaptive version. We remark that the number of steps reported in Table 5 suggest that the computational cost of an adaptive step is only 5−10% larger than that of a fixed step. This is expected in the case of non-linear problems, with the cost of computing a predictor consists of a matrix-vector multiplication, while the cost of the corrector is dominated by multiple linear system solves that are costly.
Finally, in Table 6 we report the energy E = −m · (h ap + h mc ) at the time t = 600 of the simulation with k 1 = 4 as the function of the LTE tolerance T . This time is sufficient for the magnetization vector to reach the steady-state after the reversal. We observe from the results that for the larger values of T there is a notable error in the computed energy (more than 5% for the BDF2 method with T = 10 −4 ). By contrast, the computed energy remains constant with the IMR method, irrespectively of the LTE tolerance.

The PDE Example with an Analytical Solution
This model problem is defined in the domain ⊂ R 2 , and is modelled by the non-dimensional form of the LLG equation (4). In this example we consider only the exchange coupling, hence h = ∇ 2 m. We report on the behaviour of the adaptive IMR when applied to a problem whose analytical solution has a non-trivial spatial variation. A more detailed description can be found in [38,39]. We consider an infinite two-dimensional domain, which is reduced to a finite counterpart by applying periodic BCs. For chosen constants c ∈ R and k ∈ R 3 , definê Then the components of the magnetization are given by For large intervals of the values of c and k (46) represents a damped periodic precession tending to a steady state m = [0 0 ± 1] t as t → ∞, reducing to a simple wave in space and time for α → 0. The domain = [0, 1] 2 ⊂ R 2 is used in our simulations with periodic BCs. Spatial discretisation is performed by finite elements (uniform quadrilateral elements with the resolution 80 × 80 resulting in 1681 nodes). Our spatial discretisation is similar to that used in [20]. The difference between the two approaches is that instead of a provably convergent, but restrictive linearisation procedure (the convergence is proved under a time step size restriction = O(h 2 )), we use the standard Newton method with a tolerance linked to the desired level of error in the magnetization length and/or energy, leaving to the adaptive eBDF3-IMR procedure to select the appropriate time step sizes. The initial condition is set according to (46) for t = 0 and h ap = 0. The solution parameters are k = [2π, 2π] t , c = 0.1π and h ap = 0.01. For the energy conservation experiments we set α = 0.
The snapshot of the solution at t = 0.1 is given in Fig. 7. In Fig. 8 we depict the error in the nodal magnetization where m j,n are the FE approximations at the nodes j of the magnetization vector at t = t n and V is the set of all nodes in the mesh. We consider four different configurations: BDF2, TR and IMR with Gaussian integration, and IMR with nodal (reduced) integration. When IMR with nodal quadrature is used the error (47) remains small and controlled by the adopted The integral in (48) can be evaluated exactly by both the Gaussian and the nodal quadratures when piecewise linear elements are used, since ∇m is a constant in each element. The results are summarised in Fig. 9. We observe that BDF2 does not conserve energy, but both TR and IMR with Gaussian quadratures perform well for this particular problem. The best results are obtained, however, when IMR and nodal quadratures in FE discretisation are used. Accurate results for the length and energy conservation with Gaussian quadratures are not observed for any of the time integrators in more general cases (see, for example, [13,Sect. 8.2] for the test-problem with a non-uniform external field).

The NIST MAG Standard Problem #4
The μMAG standard problem #4 [40] is a widely used benchmark for testing various micromagnetic methods and models. It involves modelling the reversal of the magnetization in a thin film of a permalloy material under an applied field. The problem is modelled by a PDE system (4)-(9), where we consider an isotrpic case (k 1 = 0). The damping parameter is α = 0.02 and we consider two cases of applied field: μ 0 h ap,1 = [−24. This choice allows the system to reach the initial condition quickly. The domain is discretised by a structured uniform mesh of hexahedral elements with a single layer of elements in the z-direction. Then we apply the hybrid FEM/BEM method described in Sect. 2.3. The parameters used in the eBDF3-IMR method are 0 = 10 −4 and T = 10 −5 . Newton's method tolerance is N = 10 −11 . We comment briefly on the computational cost of the FEM/BEM method. At each time step the eBDF3 predictor requires one sparse matrix-vector multiplication (recall that we use hierarchical approximation of the BEM matrix) at a cost O (N + N b log N b ) and a trivial solution of a linear system with a mass coefficient matrix with a cost O(N ) (when reduced integration is used, the mass matrix is diagonal by construction [13, p. 104]). The Newton iteration in the corrector takes 2-3 steps to converge [13,Fig. 8.24]. The cost of assembling the Jacobian is O(N ), as the BEM matrix does not need to be reassembled. The Krylov solver is preconditioned by a general algebraic, rather than a problem-specific, preconditioner, leading to an inevitable growth in iteration counts as N is increased. For our simulations GMRES converged in 10-25 iterations. Therefore, we pessimistically estimate the computational cost per time step of the eBDF3-IMR method as 150C mv , where C mv is the cost of a sparse matrix-vector multiplication.
We compare the accuracy of the solutions obtained by the adaptive IMR scheme to the results from [12] obtained by the IMR scheme with fixed step sizes and finite difference method applied to spatial discretisation (magnetostatic effects are handled by the FFT and monolithically coupled to the LLG model). In Fig. 10 we present the mean magnetization and the time step sizes n generated by the proposed adaptive scheme for the field h ap,1 , while Fig. 11 plots the same results for h ap,2 .
The results obtained by adaptive IMR (in red) are superimposed to those in green plotted using the raw data from [12]. The two solutions agree reasonably well, except in the interval [100, 200] for h ap,2 . It is worth noticing that this is the interval where all the solutions submited to μMAG repository [40] exhibit different behaviour. The upshot is that the proposed adaptive algorithm correctly identifies this difficult part of the time interval and responds by cutting the  Fig. 10 The evolution of the Cartesian components of the mean magnetization and the time step size n for the μMAG problem #4 for h ap,1 obtained with adaptive IMR and FEM (red) and the solution plotted using the raw data from [12] (green). One non-dimensional time unit is equal to 5.65 ps (Color figure online)  Fig. 11 The evolution of the Cartesian components of the mean magnetization and the time step size n for the μMAG problem #4 for h ap,2 obtained with adaptive IMR and FEM (red) and the solution plotted using the raw data from [12] (green). One non-dimensional time unit is equal to 5.65 ps (Color figure online) stepsize to a value significantly below the fixed steps adopted in [12]. After the oscillations in m are damped out for t > 200, the step size increases steadily. In Fig. 12 we present the mean value of m y in the critical part of the time interval calculated with two different spatial resolutions (the number of grid points in x direction is equal to the number assigned to the quadrature rule in the figure legend) and two different types of numerical quadratures. When compared to adaptive BDF2 and TR methods, IMR selects significantly smaller steps during the initial relaxation phase (for t < 0). This effect may be attributed to order reduction observed in IMR when applied to very stiff problems. For t > 0 time steps selected by all three methods are commensurate. We notice that recently introduced adaptive highorder explicit midpoint method [7] keeps the step size fairly constant both in the highlighted time interval [100, 200] and during the entire simulation for the case h ap,1 . The time step of 200 f s reported in [7] is between 3 and 20 times smaller than the step sizes taken by the eBDF3-IMR method for t > 0 (cf. Fig. 10).
The comparison of the computational cost per time step between the explicit method from [7] and the coupled FEM/BEM method is not straightforward. Results from Table 1 in [7] indicate approximately between 25 and 40 total (effective and stray) field evaluations per time step. If we assume that the cost of one field computation is proportional to C mv , we can conclude that the overall computational times of both methods are broadly comparable in this particular case. 2

Fig. 12
The mean value of m y for the μMAG problem #4 with h app,2 for adaptive IMR with different spatial resolutions and different numerical quadratures. One non-dimensional time unit is equal to 5.65 ps The maximum pointwise error (47) in the magnetization length observed during the simulation was in both cases 4 · 10 −10 , in line with the adopted Newton's tolerance.

Conclusions
In this paper we have introduced an adaptive IMR method. Using the example domain of micromagnetic problems we have shown that the method conserves the magnetisation length and accurately reproduces the energy balance (i.e. it preserves the geometric properties of the solution even for loose LTE tolerances). The method uses a computationally efficient and reliable method for the LTE estimation in the IMR and is based on the use of a third-order explicit method as the predictor. In our experiments the selected step sizes of the adaptive IMR were commensurate to the other commonly used second-order methods, while the geometric integration properties were preserved. We also remark that geometric integrator properties of IMR are relevant for other problems (for example, if the Navier-Stokes problem is solved, the kinetic energy u t u will be correctly reproduced independently of the step size, where u is the fluid velocity [3, p. 651], provided that an appropriate spatial discretisation is used, see [41]). This makes the proposed adaptive integrator suitable choice for a much wider class of problems than studied in this paper.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.