A Gauss-Seidel projection method with the minimal number of updates for stray field in micromagnetic simulations
AA Gauss-Seidel projection method with the minimal number ofupdates for stray field in micromagnetic simulations
Panchi Li a , Zetao Ma a , Rui Du a,b, ∗ , Jingrun Chen a,b, ∗ a School of Mathematical Sciences, Soochow University, Suzhou, 215006, China. b Mathematical Center for Interdisciplinary Research, Soochow University, Suzhou, 215006, China.
A B S T R A C TMagnetization dynamics in magnetic materials is often modeled by the Landau-Lifshitzequation, which is solved numerically in general. In micromagnetic simulations, the compu-tational cost relies heavily on the time-marching scheme and the evaluation of stray field.Explicit marching schemes are efficient but suffer from severe stability constraints, whilenonlinear systems of equations have to be solved in implicit schemes though they are uncon-ditionally stable. A better compromise between stability and efficiency is the semi-implicitscheme, such as the Gauss-Seidel projection method (GSPM) and the second-order backwarddifferentiation formula scheme (BDF2). At each marching step, GSPM solves several linearsystems of equations with constant coefficients and updates the stray field several times,while BDF2 updates the stray field only once but solves a larger linear system of equationswith variable coefficients and a nonsymmetric structure. In this work, we propose a newmethod, dubbed as GSPM-BDF2, by combing the advantages of both GSPM and BDF2.Like GSPM, this method is first-order accurate in time and second-order accurate in space,and is unconditionally stable with respect to the damping parameter. However, GSPM-BDF2updates the stray field only once per time step, leading to an efficiency improvement of about60% than the state-of-the-art GSPM for micromagnetic simulations. For Standard Problem
Keywords:
Micromagnetic simulation, Landau-Lifshitz equation, Gauss-Seidel projection method,Backward differentiation formula, Stray field
1. Introduction
Due to their intrinsic magnetic properties, ferromagnets have been ideal materials formagnetic recording devices in the past several decades [1]. The basic quantity of interest ina ferromagnet is the magnetization, whose dynamics is modeled by the Landau-Lifshitz (LL)equation [2, 3] phenomenologically. The model and its generalizations in the presence ofexternal controls, such as spin current and temperature gradient, have been successfullyused to interpret many interesting experimental observations. Numerically, micromagneticsimulation becomes increasingly important to study the dynamics of magnetization, in ad-dition to experiment and theory. In the LL equation, magnetization dynamics is driven by ∗ Corresponding authors e-mail:
[email protected] (Panchi Li), (Zetao Ma), [email protected] (RuiDu), [email protected] (Jingrun Chen) a r X i v : . [ m a t h . NA ] J a n he gyromagnetic term and the damping term, which are both nonlinear with respect tomagnetization. Moreover, the length of magnetization does not change during its dynamicevolution. These pose challenges in designing efficient and simple numerical methods. Be-sides, the calculation of the stray field in micromagnetic simulations is time-consuming sinceit involves a problem defined over the entire space instead of the ferromagnetic body [4].Explicit schemes such as Runge-Kutta method [5] are favored in the early days, and thefifth-order Cash-Karp Runge-Kutta method is still used in OOMMF [6, 7]. These schemesare simple and efficient in the sense that no linear/nonlinear systems of equations need tobe solved at each marching step, but the temporal stepsize is rather small due to the strongstability restriction of explicit schemes. Implicit methods, such as [8, 9, 10, 11], are provedto be unconditionally stable, and can preserve the length of magnetization automatically.However, a nonlinear system of equations with variable coefficients and a nonsymmetricstructure has to be solved at each step.A nice compromise between efficiency and stability is the semi-implicit scheme. Onenotable example is the Gauss-Seidel projection method (GSPM) [12, 13]. At each step,GSPM only solves linear systems of equations with constant coefficients seven times andupdates the stray field four times. GSPM is tested to be unconditionally stable, and is first-order accurate in time and second-order accurate in space. Recently, it has been found thatthe number of linear systems to be solved and the number of updates for the stray field ateach step can be reduced without sacrificing the stability and accuracy [14, 15]. This leadsto a reduction in computational cost of about 30% in micromagnetic simulations. Othersemi-implicit schemes, such as [16, 17, 18], are second-order accurate in time. In [17], thesemi-implicit schemes are constructed using the backward differentiation formula (BDF) andone-sided interpolation, and the second-order BDF scheme (BDF2) is proved to converge withsecond-order accuracy in both space and time [18]. Note that a projection step is needed insemi-implicit schemes to preserve the length of magnetization at each step.To compare GSPM and BDF2 in details, let us consider a discrete problem with n degreesof freedom (dofs) and the dofs of magnetization is thus 3 n . At each step, GSPM (referredto Scheme A in [14] throughout the paper) solves linear systems of equations with constantcoefficients and n dofs five times and updates the stray field with 3 n dofs three times. BDF2solves linear systems of equations once with variable coefficients, a nonsymmetric structure,and 3 n dofs and updates the stray field once with 3 n dofs. Both schemes are second-orderaccurate in space. GSPM is first-order accurate in time while BDF2 is second-order accuratein time. A natural question arises here: can one design a more efficient method whichcombines the strengths of both GSPM and BDF2 ?In the work, we provide an affirmative answer to the aforementioned question. The newmethod is dubbed as GSPM-BDF2, which solves linear systems of equations with constant2oefficients and n dofs five times and updates the stray field with 3 n dofs only once. Themethod is tested to be unconditionally stable with respect to the damping parameter, first-order accurate in time, and second-order accurate in space. In micromagnetic simulations,GSPM-BDF2 reduces the number of evaluations of the stray field from 3 to 1, yielding about60% reduction of computational time on top of GSPM [14]. For Standard Problem
2. Landau-Lifshitz equation
The dynamics of magnetization M = ( M , M , M ) T in ferromagnetic materials is mod-eled by the phenomenological Landau-Lifshitz equation [2, 3] M t = − γ M × H − γαM s M × ( M × H ) , (1)where γ is the gyromagnetic ratio, α is the dimensionless damping parameter and | M | = M s inthe point-wise sense with M s the saturation magnetization. The local effective field H = − δFδ M is obtained from the Landau-Lifshitz energy functional F [ M ] = 12 (cid:90) Ω (cid:26) AM s |∇ M | + Φ (cid:18) M M s (cid:19) − µ H e · M (cid:27) d x + µ (cid:90) R |∇ U | d x (2)with A the exchange constant and µ the permeability of vacuum. H e is the external magneticfield and Ω is the volume occupied by the material. For a uniaxial material with x -axis theeasy direction, the anisotropy energy can be described by Φ (cid:16) M M s (cid:17) = K u M s ( M + M ) with K u the anisotropy constant. The last term in (2) is the self-interacting energy induced by themagnetization distribution inside the material with U ( x ) = (cid:90) Ω ∇ N ( x − y ) · M ( y )d y , (3)where N ( x − y ) = − π | x − y | is the Newtonian potential. Fast Fourier Transform (FFT) isemployed for the evaluation of the stray field H s = −∇ U in regular-shaped materials [21, 22].For convenience, we nondimensionalize (1) by rescaling variables t → ( µ γM s ) − t and x → Lx with L the diameter of Ω. With m = M /M s , h e = H e /M s , and h s = H s /M s , thedimensionless LL equation reads as m t = − m × ( (cid:15) ∆ m + f ( m )) − α m × ( m × ( (cid:15) ∆ m + f ( m ))) , (4)3here f ( m ) = − Q ( m e + m e ) + h e + h s . (5)Here dimensionless parameters (cid:15) = A/ ( µ M s L ), Q = K u / ( µ M s ), e = (0 , , T , and e = (0 , , T . Neumann boundary condition is used ∂ m ∂ ν (cid:12)(cid:12)(cid:12) ∂ Ω = 0 , (6)where ν is the outward unit normal vector on ∂ Ω.An equivalent form of (1) is the so-called Landau-Lifshitz-Gilbert (LLG) equation M t = − γ M × H + αM s M × M t . (7)In the presence of spin transfer torque (STT), the generalized LLG equation reads as [23] M t = − γ M × H + αM s M × M t − bM s M × ( M × ( j · ∇ ) M ) − bξM s M × ( j · ∇ ) M , (8)where j is the spin polarization current density with magnitude J , b = P µ B (cid:0) eM s (1 + ξ ) (cid:1) − with e the electron charge, µ B the Bohr magneton, and P the spin current polarization. Afterrescaling t → (1 + α )( µ γM s ) − t and x → Lx , the LLG equation can be rewritten into (4)with the dimensionless local field f ( m ) of the following form f ( m ) = − Q ( m e + m e ) + h e + h s + bγM s L m × ( j · ∇ ) m + bξγM s L ( j · ∇ ) m . In particular, when an in-plane current is applied along the x -direction, the local field reducesto f ( m ) = − Q ( m e + m e ) + h e + h s + bJγM s L m × m x + bJξγM s L m x . (9)
3. The proposed method
In this section, we shall describe the proposed method in details. The GSPM [14] isintroduced firstly for completeness and comparison. For the LL equation (4) with field (5)and Neumann boundary condition (6), the fractional step is applied m ∗ − m n ∆ t = (cid:15) ∆ h m ∗ + f ( m n ) , m n +1 − m n ∆ t = − m n × m ∗ − m n ∆ t − α m n × (cid:18) m n × m ∗ − m n ∆ t (cid:19) , (10)where ∆ h represents a discrete approximation to the Laplacian operator. In 3D, we use thesecond-order centered difference∆ h m i,j,k = m i +1 ,j,k − m i,j,k + m i − ,j,k ∆ x + m i,j +1 ,k − m i,j,k + m i,j − ,k ∆ y + m i,j,k +1 − m i,j,k + m i,j,k − ∆ z , (11)4here m i,j,k = m (( i − )∆ x, ( j − )∆ y, ( k − )∆ z ) with i = 0 , , · · · , M, M + 1, j =0 , , · · · , N, N + 1, and k = 0 , , · · · , K, K + 1 being the indices of grid points in x -, y -and z -directions, respectively. For the Neumann boundary condition (6), a second-orderapproximation yields m ,j,k = m ,j,k , m M,j,k = m M +1 ,j,k , j = 1 , · · · , N, k = 1 , · · · , K, m i, ,k = m i, ,k , m i,N,k = m i,N +1 ,k , i = 1 , · · · , M, k = 1 , · · · , K, m i,j, = m i,j, , m i,j,K = m i,j,K +1 , i = 1 , · · · , M, j = 1 , · · · , N. Denote L = ( I − (cid:15) ∆ t ∆ h ) − with I the identity operator. (10) can be rewritten as m ∗ i = L ( m ni + ∆ tf i ( m n )) , i = 1 , , , (12) m n +1 = m n − m n × m ∗ − α m n × ( m n × m ∗ ) . (13)Note that (13) is updated in a Gauss-Seidel manner to avoid the stability issue [12]. It isnice that all linear systems in (12) have constant coefficients with the symmetric and positivedefinite (spd) property, thus the optimal computational cost of solving linear systems hereis O ( n ) where n is the dofs. Due to the homogeneous Neumann boundary condition, thecomputational complexity is O ( n log( n )) if discrete cosine transform is used to solve thelinear systems here [22]. Meanwhile, FFT has also been used to calculate the stray field (3)to reduce the computational complexity to O ( n log( n )) [21]. Thus solving a linear system andupdating the stray field have the same computational complexity in the asymptotic sense.However, the dofs of magnetization is actually 3 n . As a consequence, updating the stray fieldis roughly three times costly compared to solving a linear system in (12). In [14], GSPM solves the LL equation in two steps: • Implicit Gauss-Seidel update g ni = L ( m ni + ∆ tf i ( m n )) , i = 1 , , ,g ∗ i = L ( m ∗ i + ∆ tf i ( m ∗ )) , i = 1 , , (14) m ∗ m ∗ m ∗ = m n − ( m n g n − m n g n ) − α ( m n g n + m n g n + m n g n ) m n + αg n m n − ( m n g ∗ − m ∗ g n ) − α ( m ∗ g ∗ + m n g n + m n g n ) m n + αg n m n − ( m ∗ g ∗ − m ∗ g ∗ ) − α ( m ∗ g ∗ + m ∗ g ∗ + m n g n ) m n + αg n . (15) • Projection onto S m n +11 m n +12 m n +13 = 1 | m ∗ | m ∗ m ∗ m ∗ . (16) Remark 1.
In [16], the semi-implicit BDF1 scheme for the LL equation reads as m n +1 − m n ∆ t = − m n × ( (cid:15) ∆ m n +1 + f ( m n )) − α m n × ( m n × ( (cid:15) ∆ m n +1 + f ( m n ))) . (17)5 projection step is applied after (17) at each step. In fact, GSPM for the LL equation (14) - (16) can be obtained by applying the splitting strategy (12) - (13) to (17) and then updating (13) in the Gauss-Seidel manner. This motivates the current work of applying the splittingstrategy and updating in the Gauss-Seidel manner for the semi-implicit BDF2 scheme. Remark 2.
From (14) , we know that at each time step GSPM solves linear systems ofequations with n dofs five times and updates the stray field with n dofs three times.3.2. GSPM-BDF2 scheme for LL equation In [17, 18], the semi-implicit BDF2 scheme for the LL equation reads as ˆm n +2 − m n +1 + m n ∆ t = − ˜m n +2 × ( (cid:15) ∆ ˆm n +2 + f ( ˜m n +2 )) , − α ˜m n +2 × ( ˜m n +2 × ( (cid:15) ∆ ˆm n +2 + f ( ˜m n +2 ))) (18) m n +2 = 1 | ˆm n +2 | ˆm n +2 (19)with ˜m n +2 = 2 m n +1 − m n . This scheme has second-order accuracy in both space andtime. At each step, it solves one linear system of equations with variable coefficients, anonsymmetric structure, and 3 n dofs. Note that only one update is needed for the stray fieldper step in (18).Similarly, we solve the semi-implicit BDF2 scheme (18) using the splitting strategy andupdating in the Gauss-Seidel manner. The Gauss-Seidel update is applied for the followingequation after splitting ˆm n +2 − m n +1 + m n ∆ t = − ˜m n +2 × L ˜m n +2 − α ˜m n +2 × ( ˜m n +2 × L ˜m n +2 ) . (20)Therefore, in details, GSPM-BDF2 works as follows g n +2 i = L ( ˜ m n +2 i + ∆ tf i ( ˜m n +2 )) , i = 1 , , , m ∗ = 2 m n +11 − m n − ( ˜ m n +22 g n +23 − ˜ m n +23 g n +22 ) − α ( ˜ m n +21 g n +21 + ˜ m n +22 g n +22 + ˜ m n +23 g n +23 ) ˜ m n +21 + αg n +21 , ˜ m ∗ = 2 m ∗ − m n +11 , g ∗ = L ( ˜ m ∗ + ∆ tf ( ˜m n +2 )) , m ∗ = 2 m n +12 − m n − ( ˜ m n +23 g ∗ − ˜ m ∗ g n +23 ) − α ( ˜ m ∗ g ∗ + ˜ m n +22 g n +22 + ˜ m n +23 g n +23 ) ˜ m n +22 + αg n +22 , ˜ m ∗ = 2 m ∗ − m n +12 , g ∗ = L ( ˜ m ∗ + ∆ tf ( ˜m n +2 )) , m ∗ = 2 m n +13 − m n − ( ˜ m ∗ g ∗ − ˜ m ∗ g ∗ ) − α ( ˜ m ∗ g ∗ + ˜ m ∗ g ∗ + ˜ m n +23 g n +23 ) ˜ m n +23 + αg n +23 , m n +21 m n +22 m n +23 = 1 | m ∗ | m ∗ m ∗ m ∗ . In (18), one linear system with variable coefficients, a nonsymmetric structure, and 3 n dofsis solved at each step. This is avoided in GSPM-BDF2 thanks to the splitting strategy. Fivelinear systems with constant coefficients, spd structure, and n dofs are solved. Meanwhile,the Gauss-Seidel update is tested to be unconditionally stable. GSPM-BDF2 inherits the6dvantage of BDF2 that only one update of the stray field is needed with 3 n dofs. Wesummarize the main computational costs of GSPM and GSPM-BDF2 in Table 1 by countingthe number of linear systems of equations and the number of updates for the stray field perstep. Simulations suggest that one update of stray field with 3 n dofs is computationallycomparable to solving three linear systems of equations with n dofs, thus GSPM-BDF2 savesabout 60% computational time over GSPM; see Section 4 for details. Table 1. Main computational costs of GSPM and GSPM-BDF2 in micromagnetic simulations. n ) 3 (3 n )GSPM-BDF2 5 ( n ) 1 (3 n )GSPM-BDF2 is expected to be first-order accurate in time and second-order accuratein space due to the splitting strategy and the Gauss-Seidel update, both shall be verifiednumerically later. Compared to the semi-implicit BDF2 scheme, GSPM-BDF2 only solveslinear systems of equations with constant coefficients and a spd structure. Many efficientlinear solvers can be implemented directly, while the linear system in BDF2 has to be solvedby GMRES which is an iterative solver and is not so efficient as solvers for spd matrices. Remark 3.
There are two improved GSPMs proposed in [14]: Scheme A and Scheme B. Thecurrent work proposes a more efficient method by combining Scheme A and BDF2. SchemeB solves three linear systems of equations with n dofs and updates the stray field with n dofs three times. One may wonder how the combination of Scheme B and BDF2 works.Unfortunately, we find that the stability of the resulting scheme depends on the dampingparameter α for micromagnetic simulations. The underlying reason is that Scheme B updatesthe magnetization in one step and does a projection step in a subsequent step and a delay ofthe update for the stray field leads to the instability with respect to the damping parameter.
4. Numerical simulations
In 1D, we consider the following LL equation m t = − m × m xx − α m × ( m × m xx ) . (21)Choose the exact solution m e = (cos(¯ x ) sin( t ) , sin(¯ x ) sin( t ) , cos( t ))with ¯ x = x (1 − x ) . A forcing term will be added on the right-hand side of (21) with ˆf = m et + m e × m exx + α m e × ( m e × m exx ). The convergence rate with respect to thetemporal step size and the spatial step size is recorded in Table 2.7 able 2. Convergence rates in time and space for the 1D example. The final time T = 1 . e − and the damping parameter α = 0 . . ∆ x = 1 . e − is used for the temporal accuracy and ∆ t = 1 . e − is used for the spatial accuracy, respectively. Temporalaccuracy nt 1000 2000 4000 8000 orderGSPM 6.01e-08 3.01e-08 1.52e-08 7.72e-09 0.99GSPM-BDF2 6.01e-08 3.02e-08 1.52e-08 7.73e-09 0.99Spatialaccuracy nx 10 20 40 80 orderGSPM 1.24e-06 4.24e-07 1.27e-07 4.03e-08 1.66GSPM-BDF2 1.24e-06 4.24e-07 1.27e-07 4.03e-08 1.66In 3D, we choose the exact solution m e = (cos(¯ x ¯ y ¯ z ) sin( t ) , sin(¯ x ¯ y ¯ z ) sin( t ) , cos( t )) , where ¯ x = x (1 − x ) , ¯ y = y (1 − y ) and ¯ z = z (1 − z ) . Uniform discretization overΩ = [0 , is used for the spatial accuracy. Ω = [0 , × [0 , × [0 , .
2] with 128 × × O (∆ t + (∆ x ) ). Moreover, the scheme is tested tobe unconditionally stable by varying the temporal step size. Table 3. Convergence rates in time and space for the 3D example. The final time T = 1 . e − and the damping parameter α = 0 . . ∆ x = 1 . e − is used for the temporal accuracy and ∆ t = 1 . e − is used for the spatial accuracy, respectively. Temporalaccuracy nt 10 20 40 80 orderGSPM 1.00e-06 5.00e-07 2.50e-07 1.25e-07 1.00GSPM-BDF2 1.00e-06 5.00e-07 2.50e-07 1.25e-07 1.00Spatialaccuracy nx=ny=nz 6 8 10 12 orderGSPM 2.91e-14 1.72e-14 1.13e-14 7.92e-15 1.88GSPM-BDF2 2.91e-14 1.72e-14 1.13e-14 7.92e-15 1.88
The original GSPM is found to be conditionally stable with respect to the dampingparameter [13] if the stray field is updated only once per step [12]. Here we examine theperformance of GSPM-BDF2 for small damping parameters. Consider the magnetizationdynamics of a ferromagnet over Ω = 1 µ m × µ m × . µ m and the final time 1 . × × f ( m n ) is used per step. Results are shown in Fig. 1 when α = 0 .
01. The whole simulationtakes 493 .
75 seconds. If we use the GSPM (14)-(16) with three updates of the stray field,we obtain the simulation results in Fig. 2. The whole simulation takes 1225 .
92 seconds.Comparing these two results, we find that updating the stray field only once leads to thenumerical instability with respect to the damping parameter in micromagnetic simulations,although the computational cost is reduced by 60%.Simulation results of GSPM-BDF2 are plotted in Fig. 3 when α = 0 .
01. It is clear thatGSPM-BDF2 is able to produce correct magnetization dynamics though the stray field is8 ( µ m) y ( µ m ) -3-2-10123 (a) Angle profile x ( µ m) y ( µ m ) (b) Magnetization profile Fig. 1. Simulation results of the LL equation using GSPM with only one update of stray fieldat each step. The magnetization on the centered slice of the material in the xy plane is used.Left: a color plot of the angle between the in-plane magnetization and the x axis; Right : anarrow plot of the in-plane magnetization. x ( µ m) y ( µ m ) -1-0.500.51 (a) Angle profile x ( µ m) y ( µ m ) (b) Magnetization profile Fig. 2. Simulation results of the LL equation using GSPM with three updates of stray fieldat each step. The magnetization on the centered slice of the material in the xy plane is used.Left: a color plot of the angle between the in-plane magnetization and the x axis; Right : anarrow plot of the in-plane magnetization. .
39 seconds and saves 58%computational time over GSPM (14)-(16). This study shows that the proposed method isstable with respect to the damping parameter even the number of updates for the stray fieldis minimized. x ( µ m) y ( µ m ) -0.500.511.5 (a) Angle profile x ( µ m) y ( µ m ) (b) Magnetization profile Fig. 3. Simulation results of the LL equation using GSPM-BDF2. The magnetization on thecentered slice of the material in the xy plane is used. Left: a color plot of the angle betweenthe in-plane magnetization and the x axis; Right : an arrow plot of the in-plane magnetization. The setup of this benchmark problem is as follows: film size 500 nm ×
125 nm × x ( µ m) y ( µ m ) × -4 -3-2-1 Fig. 4. Magnetization profile of the initial s-state. This is generated for a random initial stateunder a magnetic field
100 mT along the [1 , , direction with a successive reduction to . Two different magnetic fields are applied: • A field 25 mT, directed 170 degrees counterclockwise from the positive x axis; • A field 36 mT, directed 190 degrees counterclockwise from the positive x axis.10wo different mesh strategies are used: a coarse mesh with cell size 5 nm × × . × . × < m x > = 0 under the external Time (s) × -9 < M > / M s -1-0.500.51
25 mT ; Right: comparison of the averaged < m y > on two different meshes. field of 25 mT is visualized in Fig. 6. The color map is given by the z-component of themagnetization. x ( µ m) y ( µ m ) -0.5-0.4-0.3-0.2-0.1 Fig. 6. Magnetization profile when the averaged < m x > = 0 under the external field of
25 mT .The color map is given by the z -component of the magnetization. In the presence of the 36 mT field, dynamics of the spatially averaged magnetization isplotted in Fig. 7 and the magnetization profile when the averaged < m x > = 0 is visualizedin Fig. 8, respectively.When the field is applied at 170 degrees, the magnetization at the center of the rectanglerotates in the same direction as that at the two ends during the reversal process. Whenthe field is applied at 190 degrees, the magnetization at the center of the rectangle rotatesin the opposite direction as that at the two ends during the reversal process. These arein good agreements with the reports listed in [19]. To check the efficiency, we record thecomputational costs of the proposed method and OOMMF in Table 4. The current work11 ime (s) × -9 < M > / M s -1-0.500.51
36 mT ; Right: comparison of the averaged < m y > on two different meshes. x ( µ m) y ( µ m ) -0.200.2 Fig. 8. Magnetization profile when the averaged < m x > = 0 under the external field of
36 mT .The color map is given by the z -component of the magnetization. Table 4. Computational costs of the proposed method and OOMMF for Standard Problem
Standard Problem
Standard Problem ×
100 nm ×
10 nm; cell size 2nm × × x, y, z ) ∈ Ω, the initial magnetization ischosen as m = g / | g | , where g ( x, y, z ) = [ − y, x, R ] (22)and R = 10 nm. The equilibrium configuration after relaxation is chosen as the initial state(vortex pattern) for all simulations in what follows.There are four sets of parameters used in Standard Problem bJ = 72 . / s and ξ = 0; time (s) × -9 < M x > , < M y > ( U n it: A / m ) × -3-2-1012
35 m / s .Left: dynamics of the spatially averaged magnetization with the result of D. G. Porter as thereference; Right: the final state colored by the x -component of magnetization. bJ = 72 . / s and ξ = 0 . bJ = 71 . / s and ξ = 0 . bJ = 57 . / s and ξ = 0 .
5. 13 ime (s) × -9 < M x > , < M y > ( U n it: A / m ) × -3-2-1012
17 m / s .Left: dynamics of the spatially averaged magnetization with the result of D. G. Porter as thereference; Right: the final state colored by the x -component of magnetization. time (s) × -9 < M x > , < M y > ( U n it: A / m ) × -3-2-1012
64 m / s .Left: dynamics of the spatially averaged magnetization with the result of D. G. Porter as thereference; Right: the final state colored by the x -component of magnetization. ime (s) × -9 < M x > , < M y > ( U n it: A / m ) × -3-2-101234
88 m / s .Left: dynamics of the spatially averaged magnetization with the result of D. G. Porter as thereference; Right: the final state colored by the x -component of magnetization. In the current work, the first component of the spatially averaged magnetization < M x > at10ns is − . × A / m, while the value is − . × A / m produced by OOMMF for case1), 2), and 3). We attribute this discrepancy to different calculations of the stray field. Thishas been observed in [7] that the spatially averaged magnetization simulated by OOMMFand M S are different due to different evaluations of the stray field. In our implementation,the stray field is evaluated in 3D without any simplification. For case 4), a vortex state withdifferent location is observed in our simulation. Both < M y > and < M x > are differentcompared to the result of D. G. Porter. However, if the result of G. Finocchio et al. isused for comparison with the parameters ξ = 0 . bJ = 72 .
45 m / s, our result has a smalldifference in < M x > but no difference in < M y > ; as shown in Fig. 13. This state is found tobe stable and may be caused by different treatments of STT terms in the LLG equation [19].For each case, GSPM-BDF2 takes about 100 seconds while softwares such as NMAG [24],OOMMF [20], and the work in [7] take more than half an hour on the same personal computer.Below we list the computational costs of GSPM-BDF2 and OOMMF in Table 5. The proposedmethod saves 96% computational time over OOMMF. Table 5. Computational costs of the proposed method and OOMMF for Standard Problem
Standard Problem ime (s) × -9 < M x > , < M y > ( U n it: A / m ) × -3-2-101234
45 m / s .Left: dynamics of the spatially averaged magnetization with the result of G. Finocchio et al.as the reference; Right: the final state colored by the x -component of magnetization. The firstcomponent of the spatially averaged magnetization < M x > at
10 ns is − . × A / m .
5. Conclusions
In this work, we propose a new method (GSPM-BDF2 in short) by combing the advantagesof the Gauss-Seidel projection method and the semi-implicit BDF2 scheme. The proposedmethod solves linear systems of equations with constant coefficients and the spd structurefive times and updates the stray field only once. The method is tested to be unconditionallystable with respect to the damping parameter, first-order accurate in time, and second-orderaccurate in space. In micromagnetic simulations, the proposed method reduces the number ofevaluations of the stray field from 3 to 1, yielding about 60% reduction of computational timeon top of GSPM [14]. For Standard Problem
Acknowledgments
P. Li is grateful to Kelong Cheng for helpful discussions during the 18th CSIAM annualmeeting and acknowledges the financial support from the Postgraduate Research & PracticeInnovation Program of Jiangsu Province via grant KYCX20 2711. The work of R. Du was16upported in part by NSFC via grant 11501399. The work of J. Chen was supported byNSFC via grant 11971021.
References [1] I. ˇZuti´c, J. Fabian, S. Das Sarma, Spintronics: Fundamentals and applications, Reviews of ModernPhysics 76 (2004) 323–410.[2] L. D. Landau, E. M. Lifshitz, On the theory of the dispersion of magetic permeability in ferromagneticbodies, “Zeitschrift fur Physik (Sowjetunion),” Physics 8 (1935) 153–169.[3] T. L. Gilbert, A lagrangian formulation of gyromagnetic equation of the magnetization field, PhysicalReview D 100 (1955) 1243–1255.[4] C. Abert, L. Exl, G. Selke, A. Drews, T. Schrefl, Numerical methods for the stray field calculation: Acomparison of recently developed algorithms, Journal of Magnetism and Magnetic Materials 326 (2013)176–185.[5] A. Romeo, G. Finocchio, M. Carpentieri, L. Torres, G. Consolo, B. Azzerboni, A numerical solutionof the magnetization reversal modeling in a permalloy thin film using fifth order Runge-Kutta methodwith adaptive step size control, Physica B: Condensed Matter 403 (2008) 1163–1194.[6] J. R. Cash, A. H. Karp, A variable order Runge-Kutta method for initial value problems with rapidlyvarying right-hand sides, ACM Transaction on Mathematical Software 16 (1999).[7] N. Massoud, K. Benjamin, B. Stellan, F. Matteo, F. Hans, V. Antoine, A. Rolf, B. Markus, M. Ulrich,P. Daniela, M. D. P. F., M. Guido, Proposal for a standard problem for micromagnetic simulationsincluding spin-transfer torque, Journal of Applied Physics 105 (2009) 113914.[8] H. Yamada, N. Hayashi, Implicit solution of the Landau-Lifshitz-Gilbert equation by the Crank-Nicolsonmethod, Journal of the Magnetics Society of Japan 28 (2004) 924–931.[9] B. S¨oren, P. Andreas, Convergence of an implicit finite element method for the Landau-Lifshitz-Gilbertequation, SIAM Journal on Numerical Analysis 44 (2006) 1405–1419.[10] D. Jeong, J. Kim, A Crank-Nicolson scheme for the Landau-Lifshitz equation without damping, Journalof Computational and Applied Mathematics 234 (2010) 613–623.[11] F. Atsushi, I. Tetsuya, T. Masayoshi, Finite difference scheme for the Landau-Lifshitz equation, JapanJournal of Industrial and Applied Mathematics 29 (2012) 83–110.[12] X.-P. Wang, C. J. Garc´ıa-Cervera, W. E, A Gauss-Seidel projection method for micromagnetics simula-tions, Journal of Computational Physics 171 (2001) 357–372.[13] C. J. Garc´ıa-Cervera, W. E, Improved Gauss-Seidel projection method for micromagnetics simulations,IEEE Transactions on Magnetics 39 (2003) 1766–1770.[14] P. Li, C. Xie, R. Du, J. Chen, X.-P. Wang, Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert equation, Journal of Computational Physics 401 (2020) 109046.[15] P. Li, J. Chen, R. Du, X.-P. Wang, Numerical methods for antiferrimagnets, IEEE Transactions onMagnetics 56 (2020) 7200509.[16] C. Ivan, Error estimates for a semi-implicit numerical scheme solving the Landau-Lifshitz equation withan exchange field, IMA Journal of Numerical Analysis (2005) 611–634.[17] C. Xie, C. J. Garc´ıa-Cervera, C. Wang, Z. Zhou, J. Chen, Second-order semi-implicit projection methodsfor micromagnetics simulations, Journal of Computational Physics (2020) 109104.[18] J. Chen, C. Wang, C. Xie, Convergence analysis of a second-order semi-implicit projection method forLandau-Lifshitz equation, arXiv 1902.0974 (2019).[19] Micromagnetic Modeling Activity Group, National institute of standards and technology, , 2000.[20] M. J. Donahue, D. G. Porter, Oommf user’s guide (2019).[21] C. J. Garc´ıa-Cervera, Numerical micromagnetics: a review, Bolet´ın de la Sociedad Espa˜nola deMatem´atica Aplicada. SeMA 39 (2007) 103–135.[22] L. Yang, Current induced domain wall motion: Analysis and simulation, Ph.D Thesis (2008).[23] S. Zhang, Z. Li, Roles of nonequilibrium conduction electrons on the magnetization dynamics of ferro-magnets, Physical Review Letter 93 (2004) 127204.[24] H. Fangohr, T. Fischbacher, M. Franchin, G. Bordignon, J. Generowicz, A. Knittel, M. Walter, M. Albert,Nmag user manual (0.2.1) (2012)., 2000.[20] M. J. Donahue, D. G. Porter, Oommf user’s guide (2019).[21] C. J. Garc´ıa-Cervera, Numerical micromagnetics: a review, Bolet´ın de la Sociedad Espa˜nola deMatem´atica Aplicada. SeMA 39 (2007) 103–135.[22] L. Yang, Current induced domain wall motion: Analysis and simulation, Ph.D Thesis (2008).[23] S. Zhang, Z. Li, Roles of nonequilibrium conduction electrons on the magnetization dynamics of ferro-magnets, Physical Review Letter 93 (2004) 127204.[24] H. Fangohr, T. Fischbacher, M. Franchin, G. Bordignon, J. Generowicz, A. Knittel, M. Walter, M. Albert,Nmag user manual (0.2.1) (2012).