Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert   equation

Panchi Li; Changjian Xie; Rui Du; Jingrun Chen; Xiao-Ping; Wang

arXiv:1907.11853·math.NA·January 8, 2020

Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert equation

Panchi Li, Changjian Xie, Rui Du, Jingrun Chen, Xiao-Ping, Wang

PDF

TL;DR

This paper introduces two improved Gauss-Seidel projection methods for the Landau-Lifshitz-Gilbert equation that are more efficient and stable, solving fewer heat equations per step while maintaining accuracy and stability.

Contribution

The paper proposes two novel Gauss-Seidel projection methods with unconditional stability that reduce computational effort compared to the original method.

Findings

01

Methods solve fewer heat equations per time step.

02

Achieve first-order accuracy in time and second-order in space.

03

Demonstrate unconditional stability and efficiency in simulations.

Abstract

In this paper, we present two improved Gauss-Seidel projection methods with unconditional stability. The first method updates the gyromagnetic term and the damping term simultaneously and follows by a projection step. The second method introduces two sets of approximate solutions, where we update the gyromagnetic term and the damping term simultaneously for one set of approximate solutions and apply the projection step to the other set of approximate solutions in an alternating manner. Compared to the original Gauss-Seidel projection method which has to solve heat equations $7$ times at each time step, the improved methods solve heat equations $5$ times and $3$ times, respectively. First-order accuracy in time and second-order accuracy in space are verified by examples in both 1D and 3D. In addition, unconditional stability with respect to both the grid size and the damping parameter is…

Tables5

Table 1. Table 1: The number of linear systems of equations to be solved and the number of FFT executions to be used at each step for the original GSPM [ 18 ] , Scheme A, and Scheme B. The savings represent the ratio between costs of two improved schemes over that of the original GSPM.

GSPM Scheme	Number of linear systems	Saving	Execution of FFT	Saving
Original	$7$	$0$	$4$	$0$
Scheme A	$5$	$2 / 7$	$3$	$1 / 4$
Scheme B	$3$	$4 / 7$	$3$	$1 / 4$

Table 2. Table 2 : Recorded CPU time in 1D with respect to the approximation error when only Δ t Δ 𝑡 \Delta t is varied and Δ x = 1 / 100 Δ 𝑥 1 100 \Delta x=1/100 .

	T/1250	T/2500	T/5000	T/10000	Reference
GSPM	7.7882e-01	1.5445e+00	3.1041e+00	6.2196e+00	-
Scheme A	4.8340e-01	9.9000e-01	2.0527e+00	4.4917e+00	-
Scheme B	3.3010e-01	6.3969e-01	1.2281e+00	2.5510e+00	-
ratio-A	0.38	0.36	0.34	0.28	0.29(2/7)
ratio-B	0.58	0.59	0.60	0.59	0.57(4/7)

Table 3. Table 3 : Recorded CPU time in 1D with respect to the approximation error when only Δ x Δ 𝑥 \Delta x is varied and Δ t = 1.0 e − 8 Δ 𝑡 1.0 𝑒 8 \Delta t=1.0e-8 .

	1/100	1/120	1/140	1/160	Reference
GSPM	3.3752e+03	5.2340e+03	9.0334e+03	1.0495e+04	-
Scheme A	2.4391e+03	3.7175e+03	6.5149e+03	8.0429e+03	-
Scheme B	1.4740e+03	2.2448e+03	3.9152e+03	4.8873e+03	-
ratio-A	0.28	0.29	0.28	0.23	0.29(2/7)
ratio-B	0.56	0.57	0.57	0.53	0.57(4/7)

Table 4. Table 4 : Recorded CPU time in 3D with respect to the approximation error when only Δ t Δ 𝑡 \Delta t is varied and the spatial mesh is 128 × 64 × 10 128 64 10 128\times 64\times 10 .

	T/10	T/20	T/40	T/80	Reference
GSPM	3.5188e+01	6.8711e+01	1.4146e+02	2.9769e+02	-
Scheme A	2.3015e+01	4.3920e+01	8.6831e+01	1.7359e+02	-
Scheme B	1.3984e+01	2.6313e+01	5.1928e+01	1.0415e+02	-
ratio-A	0.35	0.36	0.39	0.42	0.29(2/7)
ratio-B	0.60	0.62	0.63	0.65	0.57(4/7)

Table 5. Table 5 : Recorded CPU time in 3D with respect to the approximation error when only the spatial gridsize is varied with Δ x = Δ y = Δ z Δ 𝑥 Δ 𝑦 Δ 𝑧 \Delta x=\Delta y=\Delta z and Δ t = 1.0 e − 09 Δ 𝑡 1.0 𝑒 09 \Delta t=1.0e-09 .

	1/6	1/8	1/10	1/12	Reference
GSPM	2.1066e+01	9.2615e+01	1.9879e+02	3.7820e+02	-
Scheme A	1.5278e+01	6.5953e+01	1.4215e+02	2.6725e+02	-
Scheme B	8.9698e+00	3.8684e+01	8.4291e+01	1.5977e+02	-
ratio-A	0.27	0.29	0.28	0.29	0.29(2/7)
ratio-B	0.57	0.58	0.58	0.58	0.57(4/7)

Equations83

\frac{\partial M}{\partial t} = - γ M \times H - \frac{γ α}{M _{s}} M \times (M \times H)

\frac{\partial M}{\partial t} = - γ M \times H - \frac{γ α}{M _{s}} M \times (M \times H)

F[\mathbf{M}]=\frac{1}{2}\int_{\Omega}\left\{\frac{A}{M_{s}^{2}}|\nabla\mathbf{M}|^{2}+\Phi\left(\frac{\mathbf{M}}{M_{s}}\right)-2\mu_{0}\mathbf{H}_{e}\cdot\mathbf{M}\right\}\mathrm{d}\mbox{\boldmath$x$}+\frac{\mu_{0}}{2}\int_{\mathbb{R}^{3}}|\nabla U|^{2}\mathrm{d}\mbox{\boldmath$x$},

F[\mathbf{M}]=\frac{1}{2}\int_{\Omega}\left\{\frac{A}{M_{s}^{2}}|\nabla\mathbf{M}|^{2}+\Phi\left(\frac{\mathbf{M}}{M_{s}}\right)-2\mu_{0}\mathbf{H}_{e}\cdot\mathbf{M}\right\}\mathrm{d}\mbox{\boldmath$x$}+\frac{\mu_{0}}{2}\int_{\mathbb{R}^{3}}|\nabla U|^{2}\mathrm{d}\mbox{\boldmath$x$},

U(\mbox{\boldmath$x$})=\int_{\Omega}\nabla N(\mbox{\boldmath$x$}-\mbox{\boldmath$y$})\cdot\mathbf{M}(\mbox{\boldmath$y$})\mathrm{d}\mbox{\boldmath$y$},

U(\mbox{\boldmath$x$})=\int_{\Omega}\nabla N(\mbox{\boldmath$x$}-\mbox{\boldmath$y$})\cdot\mathbf{M}(\mbox{\boldmath$y$})\mathrm{d}\mbox{\boldmath$y$},

\frac{\partial m}{\partial t} = - m \times h - α m \times (m \times h),

\frac{\partial m}{\partial t} = - m \times h - α m \times (m \times h),

h = - Q (m_{2} e_{2} + m_{3} e_{3}) + ϵ Δ m + h_{e} + h_{s}

h = - Q (m_{2} e_{2} + m_{3} e_{3}) + ϵ Δ m + h_{e} + h_{s}

\frac{\partial m}{\partial ν}_{∣ \partial Ω} = 0,

\frac{\partial m}{\partial ν}_{∣ \partial Ω} = 0,

Δ_{h} m_{i, j, k}

Δ_{h} m_{i, j, k}

+ \frac{m _{i, j + 1, k} - 2 m _{i, j, k} + m _{i, j - 1, k}}{Δ y ^{2}}

+ \frac{m _{i, j, k + 1} - 2 m _{i, j, k} + m _{i, j, k - 1}}{Δ z ^{2}},

m_{0, j, k}

m_{0, j, k}

m_{i, 0, k}

m_{i, j, 0}

m_{t} = - m \times Δ m - α m \times (m \times Δ m),

m_{t} = - m \times Δ m - α m \times (m \times Δ m),

m_{t} = - m \times Δ m - α m (m \cdot Δ m) + α Δ m .

m_{t} = - m \times Δ m - α m (m \cdot Δ m) + α Δ m .

m_{t}

m_{t}

m_{t}

\frac{m ^{*} - m ^{n}}{Δ t}

\frac{m ^{*} - m ^{n}}{Δ t}

m^{n + 1}

m^{n + 1} = m^{n} - m^{n} \times (I - Δ t Δ_{h})^{- 1} m^{n},

m^{n + 1} = m^{n} - m^{n} \times (I - Δ t Δ_{h})^{- 1} m^{n},

g_{i}^{n} = (I - Δ t Δ_{h})^{- 1} m_{i}^{n}, i = 1, 2, 3.

g_{i}^{n} = (I - Δ t Δ_{h})^{- 1} m_{i}^{n}, i = 1, 2, 3.

m_{1}^{n + 1} m_{2}^{n + 1} m_{3}^{n + 1} = m_{1}^{n} + (g_{2}^{n} m_{3}^{n} - g_{3}^{n} m_{2}^{n}) m_{2}^{n} + (g_{3}^{n} m_{1}^{n + 1} - g_{1}^{n + 1} m_{3}^{n}) m_{3}^{n} + (g_{1}^{n + 1} m_{2}^{n + 1} - g_{2}^{n + 1} m_{1}^{n + 1}) .

m_{1}^{n + 1} m_{2}^{n + 1} m_{3}^{n + 1} = m_{1}^{n} + (g_{2}^{n} m_{3}^{n} - g_{3}^{n} m_{2}^{n}) m_{2}^{n} + (g_{3}^{n} m_{1}^{n + 1} - g_{1}^{n + 1} m_{3}^{n}) m_{3}^{n} + (g_{1}^{n + 1} m_{2}^{n + 1} - g_{2}^{n + 1} m_{1}^{n + 1}) .

h = ϵ Δ m + \hat{f},

h = ϵ Δ m + \hat{f},

g_{i}^{n}

g_{i}^{n}

g_{i}^{*}

m_{1}^{*} m_{2}^{*} m_{3}^{*} = m_{1}^{n} + (g_{2}^{n} m_{3}^{n} - g_{3}^{n} m_{2}^{n}) m_{2}^{n} + (g_{3}^{n} m_{1}^{*} - g_{1}^{*} m_{3}^{n}) m_{3}^{n} + (g_{1}^{*} m_{2}^{*} - g_{2}^{*} m_{1}^{*}) .

m_{1}^{*} m_{2}^{*} m_{3}^{*} = m_{1}^{n} + (g_{2}^{n} m_{3}^{n} - g_{3}^{n} m_{2}^{n}) m_{2}^{n} + (g_{3}^{n} m_{1}^{*} - g_{1}^{*} m_{3}^{n}) m_{3}^{n} + (g_{1}^{*} m_{2}^{*} - g_{2}^{*} m_{1}^{*}) .

\hat{f}^{*} = - Q (m_{2}^{*} e_{2} + m_{3}^{*} e_{3}) + h_{e} + h_{s}^{n},

\hat{f}^{*} = - Q (m_{2}^{*} e_{2} + m_{3}^{*} e_{3}) + h_{e} + h_{s}^{n},

m_{1}^{**} m_{2}^{**} m_{3}^{**} = m_{1}^{*} + α Δ t (ϵ Δ_{h} m_{1}^{**} + \hat{f}_{1}^{*}) m_{2}^{*} + α Δ t (ϵ Δ_{h} m_{2}^{**} + \hat{f}_{2}^{*}) m_{3}^{*} + α Δ t (ϵ Δ_{h} m_{3}^{**} + \hat{f}_{3}^{*}) .

m_{1}^{**} m_{2}^{**} m_{3}^{**} = m_{1}^{*} + α Δ t (ϵ Δ_{h} m_{1}^{**} + \hat{f}_{1}^{*}) m_{2}^{*} + α Δ t (ϵ Δ_{h} m_{2}^{**} + \hat{f}_{2}^{*}) m_{3}^{*} + α Δ t (ϵ Δ_{h} m_{3}^{**} + \hat{f}_{3}^{*}) .

m_{1}^{n + 1} m_{2}^{n + 1} m_{3}^{n + 1} = \frac{1}{∣ m ^{**} ∣} m_{1}^{**} m_{2}^{**} m_{3}^{**} .

m_{1}^{n + 1} m_{2}^{n + 1} m_{3}^{n + 1} = \frac{1}{∣ m ^{**} ∣} m_{1}^{**} m_{2}^{**} m_{3}^{**} .

g_{i}^{*} = (I - Δ t ϵ Δ_{h})^{- 1} (m_{i}^{*} + Δ t \hat{f}_{i}^{*}), i = 1, 2,

g_{i}^{*} = (I - Δ t ϵ Δ_{h})^{- 1} (m_{i}^{*} + Δ t \hat{f}_{i}^{*}), i = 1, 2,

\hat{f}^{*} = - Q (m_{2}^{*} e_{2} + m_{3}^{*} e_{3}) + h_{e} + h_{s}^{*},

\hat{f}^{*} = - Q (m_{2}^{*} e_{2} + m_{3}^{*} e_{3}) + h_{e} + h_{s}^{*},

g_{i}^{n}

g_{i}^{n}

g_{i}^{*}

m_{1}^{*} m_{2}^{*} m_{3}^{*} = m_{1}^{n} - (m_{2}^{n} g_{3}^{n} - m_{3}^{n} g_{2}^{n}) - α (m_{1}^{n} g_{1}^{n} + m_{2}^{n} g_{2}^{n} + m_{3}^{n} g_{3}^{n}) m_{1}^{n} + α g_{1}^{n} m_{2}^{n} - (m_{3}^{n} g_{1}^{*} - m_{1}^{*} g_{3}^{n}) - α (m_{1}^{*} g_{1}^{*} + m_{2}^{n} g_{2}^{n} + m_{3}^{n} g_{3}^{n}) m_{2}^{n} + α g_{2}^{n} m_{3}^{n} - (m_{1}^{*} g_{2}^{*} - m_{2}^{*} g_{1}^{*}) - α (m_{1}^{*} g_{1}^{*} + m_{2}^{*} g_{2}^{*} + m_{3}^{n} g_{3}^{n}) m_{3}^{n} + α g_{3}^{n} .

m_{1}^{*} m_{2}^{*} m_{3}^{*} = m_{1}^{n} - (m_{2}^{n} g_{3}^{n} - m_{3}^{n} g_{2}^{n}) - α (m_{1}^{n} g_{1}^{n} + m_{2}^{n} g_{2}^{n} + m_{3}^{n} g_{3}^{n}) m_{1}^{n} + α g_{1}^{n} m_{2}^{n} - (m_{3}^{n} g_{1}^{*} - m_{1}^{*} g_{3}^{n}) - α (m_{1}^{*} g_{1}^{*} + m_{2}^{n} g_{2}^{n} + m_{3}^{n} g_{3}^{n}) m_{2}^{n} + α g_{2}^{n} m_{3}^{n} - (m_{1}^{*} g_{2}^{*} - m_{2}^{*} g_{1}^{*}) - α (m_{1}^{*} g_{1}^{*} + m_{2}^{*} g_{2}^{*} + m_{3}^{n} g_{3}^{n}) m_{3}^{n} + α g_{3}^{n} .

m_{1}^{n + 1} m_{2}^{n + 1} m_{3}^{n + 1} = \frac{1}{∣ m ^{*} ∣} m_{1}^{*} m_{2}^{*} m_{3}^{*} .

m_{1}^{n + 1} m_{2}^{n + 1} m_{3}^{n + 1} = \frac{1}{∣ m ^{*} ∣} m_{1}^{*} m_{2}^{*} m_{3}^{*} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert equation

Panchi Li

Changjian Xie

Rui Du

Jingrun Chen

Xiao-Ping Wang

School of Mathematical Sciences, Soochow University, Suzhou, 215006, China.

Mathematical Center for Interdisciplinary Research, Soochow University, Suzhou, 215006, China.

Department of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China

Abstract

Micromagnetic simulation is an important tool to study various dynamic behaviors of magnetic order in ferromagnetic materials. The underlying model is the Landau-Lifshitz-Gilbert equation, where the magnetization dynamics is driven by the gyromagnetic torque term and the Gilbert damping term. Numerically, considerable progress has been made in the past decades. One of the most popular methods is the Gauss-Seidel projection method developed by Xiao-Ping Wang, Carlos García-Cervera, and Weinan E in 2001. It first solves a set of heat equations with constant coefficients and updates the gyromagnetic term in the Gauss-Seidel manner, and then solves another set of heat equations with constant coefficients for the damping term. Afterwards, a projection step is applied to preserve the length constraint in the pointwise sense. This method has been verified to be unconditionally stable numerically and successfully applied to study magnetization dynamics under various controls.

In this paper, we present two improved Gauss-Seidel projection methods with unconditional stability. The first method updates the gyromagnetic term and the damping term simultaneously and follows by a projection step. The second method introduces two sets of approximate solutions, where we update the gyromagnetic term and the damping term simultaneously for one set of approximate solutions and apply the projection step to the other set of approximate solutions in an alternating manner. Compared to the original Gauss-Seidel projection method which has to solve heat equations $7$ times at each time step, the improved methods solve heat equations $5$ times and $3$ times, respectively. First-order accuracy in time and second-order accuracy in space are verified by examples in both 1D and 3D. In addition, unconditional stability with respect to both the grid size and the damping parameter is confirmed numerically. Application of both methods to a realistic material is also presented with hysteresis loops and magnetization profiles. Compared with the original method, the recorded running times suggest that savings of both methods are about $2/7$ and $4/7$ for the same accuracy requirement, respectively.

keywords:

\KWDLandau-Lifshitz-Gilbert equation, Gauss-Seidel projection method, unconditional stability, micromagnetic simulation

MSC:

[2000] 35Q99 , 65Z05 , 65M06

1 Introduction

In ferromagnetic materials, the intrinsic magnetic order, known as magnetization $\mbox{\boldmath$ M $}=(M_{1},M_{2},M_{3})^{T}$ , is modeled by the following Landau-Lifshitz-Gilbert (LLG) equation [1, 2, 3]

[TABLE]

with $\gamma$ the gyromagnetic ratio and $|\mathbf{M}|=M_{s}$ the saturation magnetization. On the right-hand side of (1), the first term is the gyromagnetic term and the second term is the Gilbert damping term with $\alpha$ the dimensionless damping coefficient [2]. Note that the gyromagnetic term is a conservative term, whereas the damping term is a dissipative term. The local field $\boldsymbol{\mathcal{H}}=-\frac{\delta F}{\delta\mathbf{M}}$ is computed from the Landau-Lifshitz energy functional

[TABLE]

where $A$ is the exchange constant, $\frac{A}{M_{s}^{2}}|\nabla\mathbf{M}|$ is the exchange interaction energy; $\Phi\left(\frac{\mathbf{M}}{M_{s}}\right)$ is the anisotropy energy, and for simplicity the material is assumed to be uniaxial with $\Phi\left(\frac{\mathbf{M}}{M_{s}}\right)=\frac{K_{u}}{M_{s}^{2}}(M_{2}^{2}+M_{3}^{2})$ with $K_{u}$ the anisotropy constant; $-2\mu_{0}\mathbf{H}_{e}\cdot\mathbf{M}$ is the Zeeman energy due to the external field with $\mu_{0}$ the permeability of vacuum. $\Omega$ is the volume occupied by the material. The last term in (2) is the energy resulting from the field induced by the magnetization distribution inside the material. This stray field $\boldsymbol{\mathbf{H}}_{s}=-\nabla U$ where $U(\mbox{\boldmath$ x $})$ satisfies

[TABLE]

where $N(\mbox{\boldmath$ x $}-\mbox{\boldmath$ y $})=-\frac{1}{4\pi}\frac{1}{|\mbox{\boldmath$ x $}-\mbox{\boldmath$ y $}|}$ is the Newtonian potential.

For convenience, we rescale the original LLG equation (1) by changes of variables $t\rightarrow(\mu_{0}\gamma M_{s})^{-1}t$ and $x\rightarrow Lx$ with $L$ the diameter of $\Omega$ . Define $\mathbf{m}=\mbox{\boldmath$ M $}/M_{s}$ and $\mathbf{h}=M_{s}\mathcal{H}$ . The dimensionless LLG equation reads as

[TABLE]

where

[TABLE]

with dimensionless parameters $Q=K_{u}/(\mu_{0}M_{s}^{2})$ and $\epsilon=A/(\mu_{0}M_{s}^{2}L^{2})$ . Here $\mathbf{e}_{2}=(0,1,0)$ , $\mathbf{e}_{3}=(0,0,1)$ . Neumann boundary condition is used

[TABLE]

where $\boldsymbol{\nu}$ is the outward unit normal vector on $\partial\Omega$ .

The LLG equation is a weakly nonlinear equation. In the absence of Gilbert damping, $\alpha=0$ , equation (4) is a degenerate equation of parabolic type and is related to the sympletic flow of harmonic maps [4]. In the large damping limit, $\alpha\to\infty$ , equation (4) is related to the heat flow for harmonic maps [5]. It is easy to check that $\lvert\mathbf{m}\rvert=1$ in the pointwise sense in the evolution. All these properties possesses interesting challenges for designing numerical methods to solve the LLG equation. Meanwhile, micromagnetic simulation is an important tool to study magnetization dynamics of magnetic materials [3, 6]. Over the past decades, there has been increasing progress on numerical methods for the LLG equation; see [7, 8, 9] for reviews and references therein. Finite difference method and finite element method have been used for the spatial discretization.

For the temporal discretization, there are explicit schemes such as Runge-Kutta methods [10, 11]. Their stepsizes are subject to strong stability constraint. Another issue is that the length of magnetization cannot be preserved and thus a projection step is needed. Implicit schemes [12, 13, 14] are unconditionally stable and usually can preserve the length of magnetization automatically. The difficulty of implicit schemes is how to solve a nonlinear system of equations at each step. Therefore, semi-implicit methods [15, 16, 17, 18, 19] provide a compromise between stability and the difficult for solving the equation at each step. A projection step is also needed to preserve the length of magnetization.

Among the semi-implicit schemes, the most popular one is the Gauss-Seidel projection method (GSPM) proposed by Wang, García-Cervera, and E [15, 18]. GSPM first solves a set of heat equations with constant coefficients and updates the gyromagnetic term in the Gauss-Seidel manner, and then solves another set of heat equations with constant coefficients for the damping term. Afterwards, a projection step is applied to preserve the length of magnetization. GSPM is first-order accurate in time and has been verified to be unconditionally stable numerically.

In this paper, we present two improved Gauss-Seidel projection methods with unconditional stability. The first method updates the gyromagnetic term and the damping term simultaneously and follows by a projection step. The second method introduces two sets of approximate solutions, where we update the gyromagnetic term and the damping term simultaneously for one set of approximate solutions and apply the projection step to the other set of approximate solutions in an alternating manner. Compared to the original Gauss-Seidel projection method, which solves heat equations $7$ times at each time step, the improved methods solve heat equations $5$ times and $3$ times, respectively. First-order accuracy in time and second-order accuracy in space are verified by examples in both 1D and 3D. In addition, unconditional stability with respect to both the grid size and the damping parameter is confirmed numerically. Application of both methods to a realistic material is also presented with hysteresis loops and magnetization profiles. Compared with the original method, the recorded running times suggest that savings of both methods are about $2/7$ and $4/7$ for the same accuracy requirement, respectively.

The rest of the paper is organized as follows. For completeness and comparison, we first introduce GSPM in Section 2. Two improved GSPMs are presented in Section 3. Detailed numerical tests are given in Section 4, including accuracy check and efficiency check in both 1D and 3D, unconditional stability with respect to both the grid size and the damping parameter, hysteresis loops, and magnetization profiles. Conclusions are drawn in Section 5.

2 Gauss-Seidel projection method for Landau-Lifshitz-Gilbert equation

Before the introduction of the GSPM [15, 18], we first use the finite difference method for spatial discretization. Figure 1 shows a schematic picture of spatial grids in 1D. Let $i=0,1,\cdots,M,M+1$ , $j=0,1,\cdots,N,N+1$ , and $k=0,1,\cdots,K,K+1$ be the indices of grid points in 3D.

Second-order centered difference for $\Delta\mathbf{m}$ reads as

[TABLE]

where $\mathbf{m}_{i,j,k}=\mathbf{m}((i-\frac{1}{2})\Delta x,(j-\frac{1}{2})\Delta y,(k-\frac{1}{2})\Delta z)$ . For the Neumann boundary condition, a second-order approximation yields

[TABLE]

To illustrate the main ideas, we first consider the following simplified equation

[TABLE]

which can be rewritten as

[TABLE]

We split (8) into two equations

[TABLE]

However, (9) is still nonlinear. Therefore, we consider a fractional step scheme to solve (9)

[TABLE]

or

[TABLE]

where $I$ is the identity matrix. This scheme is subject to strong stability constraint, and thus the implicit Gauss-Seidel scheme is introduced to overcome this issue. Let

[TABLE]

We then have

[TABLE]

This scheme solve (9) with unconditional stability. (10) is linear heat equation which can be solved easily. However, the splitting scheme (9) - (10) cannot preserve $\lvert\mbox{\boldmath$ m $}\rvert=1$ , and thus a projection step needs to be added.

For the full LLG equation (4), the GSPM works as follows. Define

[TABLE]

where $\mathbf{\hat{f}}=-Q(m_{2}\mathbf{e}_{2}+m_{3}\mathbf{e}_{3})+\mathbf{h}_{e}+\mathbf{h}_{s}$ .

The original GSPM [15] solves the equation (4) in three steps:

Implicit Gauss-Seidel

[TABLE]

2.

Heat flow without constraints

[TABLE]

3.

Projection onto $S^{2}$

[TABLE]

Here the numerical stability of the original GSPM [15] was founded to be independent of gridsizes but depend on the damping parameter $\alpha$ . This issue was solved in [18] by replacing (14) and (16) with

[TABLE]

and

[TABLE]

respectively. Update of the stray field is done using fast Fourier transform [15]. It is easy to see that the GSPM solves $7$ linear systems of equations with constant coefficients and updates the stray field using FFT $6$ times at each step.

3 Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert equation

Based on the description of the original GSPM in Section 2, we introduce two improved GSPMs for LLG equation. The first improvement updates both the gyromagnetic term and the damping term simultaneously, termed as Scheme A. The second improvement introduces two sets of approximate solution with one set for implicit Gauss-Seidel step and the other set for projection in an alternating manner, termed as Scheme B. Details are given in below.

3.1 Scheme A

The main improvement of Scheme A over the original GSPM is the combination of (13) - (17), or (9) - (10).

Implicit-Gauss-Seidel

[TABLE]

2.

Projection onto $S^{2}$

[TABLE]

It is easy to see that Scheme A solves $5$ linear systems of equations with constant coefficients and uses FFT $5$ times at each step.

3.2 Scheme B

The main improvement of Scheme B over Scheme A is the introduction of two sets of approximate solutions, one for (1) - (20) and the other for (21) and the update of these two sets of solutions in an alternating manner.

Given the initialized $\mathbf{g}^{0}$

[TABLE]

Scheme B works as follows

Implicit Gauss-Seidel

[TABLE]

2.

Projection onto $S^{2}$

[TABLE]

Here one set of approximate solution $\{\mbox{\boldmath$ m $}^{*}\}$ is updated in the implicit Gauss-Seidel step and the other set of approximate solution $\{\mbox{\boldmath$ m $}^{n+1}\}$ is updated in the projection step. Note that (23) is defined only for $\{\mbox{\boldmath$ m $}^{*}\}$ which can be used in two successive temporal steps, and thus only $3$ linear systems of equations with constant coefficients are solved at each step and $3$ FFT executions are used for the stray field. The length of magnetization can be preserved in the time evolution.

The computational cost of GSPM and its improvements comes from solving the linear systems of equations with constant coefficients. To summarize, we list the number of linear systems of equations to be solved and the number of FFT executions to be used at each step for the original GSPM [18], Scheme A, and Scheme B in table 1. The savings represent the ratio between costs of two improved schemes over that of the original GSPM.

4 Numerical Experiments

In this section, we compare the original GSPM [15, 18], Scheme A, and Scheme B via a series of examples in both 1D and 3D, including accuracy check and efficiency check, unconditional stability with respect to both the grid size and the damping parameter, hysteresis loops, and magnetization profiles. For convenience, we define

[TABLE]

for $i=$ A and B, which quantifies the improved efficiency of Scheme A and Scheme B over the original GSPM [15, 18].

4.1 Accuracy Test

Example 4.1 (1D case)

In 1D, we choose the exact solution over the unit interval $\Omega=[0,1]$

[TABLE]

which satisfies

[TABLE]

with $\bar{x}=x^{2}(1-x)^{2}$ , and $\mathbf{f}=\mathbf{m}_{et}+\mathbf{m}_{e}\times\mathbf{m}_{exx}+\alpha\mathbf{m}_{e}\times(\mathbf{m}_{e}\times\mathbf{m}_{exx})$ . Parameters are $\alpha=0.00001$ and $T=5.0e-2$ .

We first show the error $\|\mathbf{m}_{e}-\mathbf{m}_{h}\|_{\infty}$ with $\mathbf{m}_{h}$ being the numerical solution with respect to the temporal stepsize $\Delta t$ and the spatial stepsize $\Delta x$ . As shown in fig. 2(a) and fig. 2(c), suggested by the least squares fitting, both first-order accuracy in time and second-order accuracy in space are observed. Meanwhile, we record the CPU time as a function of accuracy (error) by varying the temporal stepsize and the spatial stepsize in fig. 2(b) and fig. 2(d), table 2 and table 3, respectively. In addition, from table 2 and table 3, the saving of Scheme A over GSPM is about $2/7$ , which equals $1-5/7$ , and the saving of Scheme B over GSPM is about $4/7$ , respectively. This observation is in good agreement with the number of linear systems being solved at each step for these three methods, as shown in table 1.

Example 4.2 (3D case)

In 3D, we choose the exact solution over $\Omega=[0,2]\times[0,1]\times[0,0.2]$

[TABLE]

which satisfies

[TABLE]

with $\bar{x}=x^{2}(1-x)^{2}$ , $\bar{y}=y^{2}(1-y)^{2}$ , $\bar{z}=z^{2}(1-z)^{2}$ and $\mathbf{f}=\mathbf{m}_{et}+\mathbf{m}_{e}\times\Delta\mathbf{m}_{e}+\alpha\mathbf{m}_{e}\times(\mathbf{m}_{e}\times\Delta\mathbf{m}_{e})$ . Parameters are $T=1.0e-05$ and $\alpha=0.01$ .

Like in the 1D case, we first show the error $\|\mathbf{m}_{e}-\mathbf{m}_{h}\|_{\infty}$ with $\mathbf{m}_{h}$ being the numerical solution with respect to the temporal stepsize $\Delta t$ and the spatial stepsize $\Delta x$ . As shown in fig. 3(a) and fig. 3(c), suggested by the least squares fitting, both first-order accuracy in time and second-order accuracy in space are observed. Meanwhile, we record the CPU time as a function of accuracy (error) by varying the temporal stepsize and the spatial stepsize in fig. 3(b) and fig. 3(d), table 4 and table 5, respectively. In addition, from table 4 and table 5, the saving of Scheme A over GSPM is about $2/7$ , and the saving of Scheme B over GSPM is about $4/7$ , respectively. This observation is in good agreement with the number of linear systems being solved at each step for these three methods, as shown in table 1.

It worths mentioning that all these three methods are tested to be unconditionally stable with respect to the spatial gridsize and the temporal stepsize.

4.2 Micromagnetic Simulations

To compare the performance of Scheme A and Scheme B with GSPM, we have carried out micromagnetic simulations of the full LLG equation with realistic material parameters. In all our following simulations, we consider a thin film ferromagnet of size $\Omega=1\;\mu\textrm{m}\times 1\;\mu\textrm{m}\times 0.02\;\mu\textrm{m}$ with the spatial gridsize $4\;\textrm{nm}\times 4\;\textrm{nm}\times 4\;\textrm{nm}$ and the temporal stepsize ${\Delta t}=1$ picosecond. The demagnetization field (stray field) is calculated via FFT [15, 18].

4.2.1 Comparison of hysteresis loops

The hysteresis loop is calculated in the following way. First, a positive external field $H_{0}=\mu_{0}H$ is applied and the system is allowed to reach a stable state. Afterwards, the external field is reduced by a certain amount and the system is relaxed to a stable state again. The process continues until the external field attains a negative field of strength $H_{0}$ . Then the external field starts to increase and the system relaxes until the initial applied external field $H_{0}$ is approached. In the hysteresis loop, we can monitor the magnetization dynamics and plot the average magnetization at the stable state as a function of the strength of the external field. The stopping criterion for a steady state is that the relative change of the total energy is less than $10^{-7}$ . The applied field is parallel to the $x$ axis. The initial state we take is the uniform state and the damping parameter $\alpha=0.1$ .

In fig. 4, we compare the average magnetization in the hysteresis loop simulated by GSPM, Scheme A and Scheme B. Profiles of the average magnetization of these three methods are in quantitative agreements with approximately the same switch field $9\;(\pm 0.4)\;\mathrm{mT}$ .

4.2.2 Comparison of magnetization profiles

It is tested that GSPM in [15] was unstable with a very small damping parameter $\alpha$ and was resolved in [18]. This section is devoted to the unconditional stability of Scheme A and Scheme B with respect to $\alpha$ . We consider a thin film ferromagnet of size $1\;\mu\textrm{m}\times 1\;\mu\textrm{m}\times 0.02\;\mu\textrm{m}$ with the spatial gridsize $4\;\textrm{nm}\times 4\;\textrm{nm}\times 4\;\textrm{nm}$ and the temporal stepsize is $1$ picosecond. Following [18], we consider the full LLG equation with $\alpha=0.1$ and $\alpha=0.01$ and without the external field. The initial state is $\mbox{\boldmath$ m $}_{0}=(0,1,0)$ if $x\in[0,L_{x}/5]\cup[4L_{x}/5,L_{x}]$ and $\mbox{\boldmath$ m $}_{0}=(1,0,0)$ otherwise. The final time is $10\;\textrm{ns}$ . In Figures 5, 6 and 7, we present a color plot of the angle between the in-plane magnetization and the $x$ axis, and an arrow plot of the in-plane magnetization for the original GSPM [15], Scheme A, and Scheme B, respectively. In these figures, $\alpha=0.1$ is presented in the top row and $\alpha=0.01$ is presented in the bottom row; a color plot of the angle between the in-palne magnetization and the $x$ axis is presented in the left column and an arrow plot of the in-plane magnetization is presented in the right column.

5 Conclusion

In this paper, based on the original Gauss-Seidel projection methods, we present two improved Gauss-Seidel projection methods with the first-order accuracy in time and the second-order accuracy in space. The first method updates the gyromagnetic term and the damping term simultaneously and follows by a projection step, which requires to solve heat equations $5$ times at each time step. The second method introduces two sets of approximate solutions, where we update the gyromagnetic term and the damping term simultaneously for one set of approximate solutions and apply the projection step to the other set of approximate solutions in an alternating manner. Therefore, only $3$ heat equations are needed to be solved at each step. Compared to the original Gauss-Seidel projection method, which solves heat equations $7$ times at each step, savings of these two improved methods are about $2/7$ and $4/7$ , which is verified by both 1D and 3D examples for the same accuracy requirement. In addition, unconditional stability with respect to both the grid size and the damping parameter is confirmed numerically. Application of both methods to a realistic material is also presented with hysteresis loops and magnetization profiles.

Acknowledgments

This work is supported in part by the grants NSFC 21602149 (J. Chen), NSFC 11501399 (R. Du), the Hong Kong Research Grants Council (GRF grants 16302715, 16324416, 16303318 and NSFC-RGC joint research grant N-HKUST620/15) (X.-P. Wang), and the Innovation Program for postgraduates in Jiangsu province via grant KYCX19_1947 (C. Xie).

Bibliography19

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Landau and Lifshitz [1935] L. Landau, E. Lifshitz, On the theory of the dispersion of magetic permeability in ferromagnetic bodies, Phys. Z. Sowjetunion 8 (1935) 153–169.
2Gilbert [1955] T. Gilbert, A lagrangian formulation of gyromagnetic equation of the magnetization field, Phys. Rev. 100 (1955) 1243–1255.
3Jr. [1963] W. F. B. Jr., Micromagnetics, Interscience Tracts on Physics and Astronomy, 1963.
4Sulem et al. [1986] P. Sulem, C. Sulem, C. Bardos, On the continuous limit limit for a system of classical spins, Comm. Math. Phys. 107 (1986) 431–454.
5Struwe [1988] M. Struwe, On the evolution of harmonic maps in higher dimensions, J. Differential Geom. 28 (1988) 485–502.
6Žutić et al. [2004] I. Žutić, J. Fabian, S. Das Sarma, Spintronics: Fundamentals and applications, Rev. Mod. Phys. 76 (2004) 323–410.
7Kruzik and Prohl [2006] M. Kruzik, A. Prohl, Recent developments in the modeling, analysis, and numerics of ferromagnetism, SIAM Rev. 48 (2006) 439–483.
8Cimrák [2008] I. Cimrák, A survey on the numerics and computations for the Landau-Lifshitz equation of micromagnetism, Arch. Comput. Methods Eng. 15 (2008) 277–309.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert equation

Abstract

keywords:

MSC:

1 Introduction

2 Gauss-Seidel projection method for Landau-Lifshitz-Gilbert equation

3 Two improved Gauss-Seidel projection methods for Landau-Lifshitz-Gilbert equation

3.1 Scheme A

3.2 Scheme B

4 Numerical Experiments

4.1 Accuracy Test

Example 4.1** (1D case)**

Example 4.2** (3D case)**

4.2 Micromagnetic Simulations

4.2.1 Comparison of hysteresis loops

4.2.2 Comparison of magnetization profiles

5 Conclusion

Acknowledgments

Example 4.1 (1D case)

Example 4.2 (3D case)