A comparative study of finite element methods for a class of harmonic map heat flow problems

TL;DR

This paper systematically compares three finite element methods for solving harmonic map heat flow problems from a disk to a sphere, focusing on convergence, efficiency, and stability within a unified framework.

Contribution

It provides a comprehensive comparison and analysis of three finite element discretization methods, including a new stability result for one method.

Findings

All methods show expected convergence rates for smooth solutions

Computational efficiency varies among the methods

A discrete inf-sup stability result is established for one method

Abstract

In this paper, we review and systematically compare three finite element discretization methods for a harmonic map heat flow problem from the unit disk in $\mathbb{R}^2$ to the unit sphere in $\mathbb{R}^3$ in an unified framework. Numerical tests validate the convergence rates in a regime of smooth solutions and are used to compare the methods in terms of computational efficiency. For one of the methods a discrete inf-sup stability result is derived.

Tables

Table 1. Table 1: Error for ( 12 ) with BDF1 and V h 2 V_{h}^{2} ; h = 2 − 14 h=2^{-14} fixed.

$h=2^{-14}$	${\parallel u_h^J-u_h^{\text{ref},J}\parallel }_{L^2\left(I\right)}$	EOC	${\parallel u_h^J-u_h^{\text{ref},J}\parallel }_{H^1\left(I\right)}$	EOC
$\tau =5\cdot {10}^{-2}$	$4.2056e-02$	$-$	$1.3876e-01$	$-$
$\tau =2.5\cdot {10}^{-2}$	$2.3299e-02$	$0.85$	$7.6855e-02$	$0.85$
$\tau =1.25\cdot {10}^{-2}$	$1.2331e-02$	$0.92$	$4.0678e-02$	$0.92$
$\tau =6.25\cdot {10}^{-3}$	$6.3542e-03$	$0.96$	$2.0964e-02$	$0.96$
$\tau =3.125\cdot {10}^{-3}$	$3.2269e-03$	$0.98$	$1.0647e-02$	$0.98$

Table 2. Table 2: Error for ( 12 ) with BDF2 and V h 2 V_{h}^{2} ; h = 2 − 14 h=2^{-14} fixed.

$h=2^{-14}$	${\parallel u_h^J-u_h^{\text{ref},J}\parallel }_{L^2\left(I\right)}$	EOC	${\parallel u_h^J-u_h^{\text{ref},J}\parallel }_{H^1\left(I\right)}$	EOC
$\tau =5\cdot {10}^{-2}$	$2.5176e-02$	$-$	$8.3363e-02$	$-$
$\tau =2.5\cdot {10}^{-2}$	$5.6626e-03$	$2.15$	$1.9086e-02$	$2.13$
$\tau =1.25\cdot {10}^{-2}$	$1.0739e-03$	$2.40$	$3.6513e-03$	$2.39$
$\tau =6.25\cdot {10}^{-3}$	$2.3314e-04$	$2.20$	$7.8868e-04$	$2.21$
$\tau =3.125\cdot {10}^{-3}$	$5.5458e-05$	$2.07$	$1.8704e-04$	$2.08$

Table 3. Table 3: Error for ( 12 ) with V h 1 V_{h}^{1} and BDF2; τ = 10 − 6 \tau=10^{-6} fixed.

$\tau ={10}^{-6}$	${\parallel u_h^J-u_h^{\text{ref},J}\parallel }_{L^2\left(I\right)}$	EOC	${\parallel u_h^J-u_h^{\text{ref},J}\parallel }_{H^1\left(I\right)}$	EOC
$h=2^{-3}$	$1.0893e-03$	$-$	$2.4522e-02$	$-$
$h=2^{-4}$	$2.7888e-04$	$1.97$	$1.2328e-02$	$0.99$
$h=2^{-5}$	$6.9802e-05$	$2.00$	$6.1812e-03$	$1.00$
$h=2^{-6}$	$1.7100e-05$	$2.03$	$3.0995e-03$	$1.00$
$h=2^{-7}$	$3.9710e-06$	$2.11$	$1.5422e-03$	$1.01$

Table 4. Table 4: Error for ( 12 ) with V h 2 V_{h}^{2} and BDF2; τ = 10 − 6 \tau=10^{-6} fixed.

$\tau ={10}^{-6}$	${\parallel u_h^J-u_h^{\text{ref},J}\parallel }_{L^2\left(I\right)}$	EOC	${\parallel u_h^J-u_h^{\text{ref},J}\parallel }_{H^1\left(I\right)}$	EOC
$h=2^{-3}$	$3.5116e-05$	$-$	$2.0124e-03$	$-$
$h=2^{-4}$	$4.1838e-06$	$3.07$	$5.0027e-04$	$2.01$
$h=2^{-5}$	$5.1173e-07$	$3.03$	$1.2469e-04$	$2.00$
$h=2^{-6}$	$6.3310e-08$	$3.01$	$3.1125e-05$	$2.00$
$h=2^{-7}$	$7.8741e-09$	$3.01$	$7.7756e-06$	$2.00$

Table 5. Table 5: Error for PPFEM ( 18 )-( 19 ) with BDF1 and 𝐕 h 2 \mathbf{{V}}_{h}^{2} ; h = 2 − 6 h=2^{-6} fixed.

$h=2^{-6}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1}$	EOC
$\tau =5\cdot {10}^{-2}$	$3.4910e-02$	$-$	$1.4780e-01$	$-$
$\tau =2.5\cdot {10}^{-2}$	$2.3253e-02$	$0.59$	$9.6176e-02$	$0.62$
$\tau =1.25\cdot {10}^{-2}$	$1.3985e-02$	$0.73$	$5.7423e-02$	$0.74$
$\tau =6.25\cdot {10}^{-3}$	$7.8473e-03$	$0.83$	$3.2146e-02$	$0.84$
$\tau =3.125\cdot {10}^{-3}$	$4.2028e-03$	$0.90$	$1.7204e-02$	$0.90$
$\tau =1.5625\cdot {10}^{-3}$	$2.1849e-03$	$0.94$	$8.9428e-03$	$0.94$
$\tau =7.8125\cdot {10}^{-4}$	$1.1158e-03$	$0.97$	$4.5677e-03$	$0.97$

Table 6. Table 6: Error for PPFEM ( 18 )-( 19 ) with BDF2 and 𝐕 h 2 \mathbf{{V}}_{h}^{2} ; h = 2 − 6 h=2^{-6} fixed.

$h=2^{-8}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1}$	EOC
$\tau =5\cdot {10}^{-2}$	$2.5378e-02$	$-$	$1.0438e-01$	$-$
$\tau =2.5\cdot {10}^{-2}$	$7.2493e-03$	$1.81$	$2.9959e-02$	$1.80$
$\tau =1.25\cdot {10}^{-2}$	$1.6907e-03$	$2.10$	$6.9553e-03$	$2.11$
$\tau =6.25\cdot {10}^{-3}$	$4.1608e-04$	$2.02$	$1.7098e-03$	$2.02$
$\tau =3.125\cdot {10}^{-3}$	$1.0539e-04$	$1.98$	$4.3477e-04$	$1.98$

Table 7. Table 7: Error for PPFEM ( 18 )-( 19 ) with 𝐕 h 1 \mathbf{{V}}_{h}^{1} and BDF2; τ = 10 − 6 \tau=10^{-6} fixed.

$\tau ={10}^{-6}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1}$	EOC
$h=2^{-2}$	$3.2225e-02$	$-$	$1.8987e-01$	$-$
$h=2^{-3}$	$8.1866e-03$	$1.98$	$8.5117e-02$	$1.16$
$h=2^{-4}$	$1.8507e-03$	$2.15$	$3.4665e-02$	$1.30$
$h=2^{-5}$	$4.2263e-04$	$2.13$	$1.2632e-02$	$1.45$
$h=2^{-6}$	$1.0192e-04$	$2.05$	$4.5444e-03$	$1.47$

Table 8. Table 8: Error for PPFEM ( 18 )-( 19 ) with 𝐕 h 2 \mathbf{{V}}_{h}^{2} and BDF2; τ = 10 − 6 \tau=10^{-6} fixed.

$\tau ={10}^{-6}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2\left(I\right)}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1\left(I\right)}$	EOC
$h=2^{-2}$	$8.2075e-02$	$-$	$3.4101e-01$	$-$
$h=2^{-3}$	$3.6998e-02$	$1.15$	$1.5203e-01$	$1.17$
$h=2^{-4}$	$2.5116e-03$	$3.88$	$1.0399e-02$	$3.87$
$h=2^{-5}$	$1.6563e-04$	$3.92$	$6.9987e-04$	$3.89$
$h=2^{-6}$	$1.1221e-05$	$3.88$	$4.8919e-05$	$3.84$

Table 9. Table 9: Error for TFEM (( 27 )-( 28 ) and ( 26 )) with BDF1 and 𝐕 h 2 \mathbf{{V}}_{h}^{2} ; h = 2 − 6 h=2^{-6} fixed.

$h=2^{-6}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1}$	EOC
$\tau =5\cdot {10}^{-2}$	$1.0814e-01$	$-$	$4.3607e-01$	$-$
$\tau =2.5\cdot {10}^{-2}$	$6.1103e-02$	$0.82$	$2.4512e-01$	$0.83$
$\tau =1.25\cdot {10}^{-2}$	$3.2653e-02$	$0.90$	$1.3072e-01$	$0.91$
$\tau =6.25\cdot {10}^{-3}$	$1.6902e-02$	$0.95$	$6.7616e-02$	$0.95$
$\tau =3.125\cdot {10}^{-3}$	$8.6006e-03$	$0.97$	$3.4404e-02$	$0.97$
$\tau =1.5625\cdot {10}^{-3}$	$4.3381e-03$	$0.99$	$1.7355e-02$	$0.99$
$\tau =7.8125\cdot {10}^{-4}$	$2.1783e-03$	$0.99$	$8.7161e-03$	$0.99$

Table 10. Table 10: Error for TFEM (( 27 )-( 28 ) and ( 26 )) with BDF2 and 𝐕 h 2 \mathbf{{V}}_{h}^{2} ; h = 2 − 6 h=2^{-6} fixed.

$h=2^{-6}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1}$	EOC
$\tau =5\cdot {10}^{-2}$	$7.4457e-02$	$-$	$2.9762e-01$	$-$
$\tau =2.5\cdot {10}^{-2}$	$2.0902e-02$	$1.83$	$8.4048e-02$	$1.82$
$\tau =1.25\cdot {10}^{-2}$	$5.2568e-03$	$1.99$	$2.2060e-02$	$1.93$
$\tau =6.25\cdot {10}^{-3}$	$1.3969e-03$	$1.91$	$6.3743e-03$	$1.80$
$\tau =3.125\cdot {10}^{-3}$	$3.6834e-04$	$1.92$	$1.8829e-03$	$1.76$

Table 11. Table 11: Error for TFEM (( 27 )-( 28 ) and ( 26 )) with 𝐕 h 1 \mathbf{{V}}_{h}^{1} and BDF2; τ = 10 − 6 \tau=10^{-6} fixed.

$\tau ={10}^{-6}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1}$	EOC
$h=2^{-2}$	$3.7887e-02$	$-$	$2.0071e-01$	$-$
$h=2^{-3}$	$9.2785e-03$	$2.03$	$8.7882e-02$	$1.19$
$h=2^{-4}$	$2.0418e-03$	$2.18$	$3.5116e-02$	$1.32$
$h=2^{-5}$	$4.6612e-04$	$2.13$	$1.2721e-02$	$1.46$
$h=2^{-6}$	$1.1216e-04$	$2.06$	$4.5600e-03$	$1.48$

Table 12. Table 12: Error for TFEM (( 27 )-( 28 ) and ( 26 )) with 𝐕 h 2 \mathbf{{V}}_{h}^{2} and BDF2; τ = 10 − 6 \tau=10^{-6} fixed.

$\tau ={10}^{-6}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1}$	EOC
$h=2^{-2}$	$5.5044e-03$	$-$	$4.4757e-02$	$-$
$h=2^{-3}$	$6.0176e-04$	$3.19$	$7.0017e-03$	$2.68$
$h=2^{-4}$	$7.0011e-05$	$3.10$	$1.1024e-03$	$2.67$
$h=2^{-5}$	$8.3761e-06$	$3.06$	$1.8427e-04$	$2.58$
$h=2^{-6}$	$1.0403e-06$	$3.01$	$3.4234e-05$	$2.43$

Table 13. Table 13: Error for Algorithm 1 with ϵ = 10 − 10 \epsilon=10^{-10} ; τ = 10 − 6 \tau=10^{-6} fixed.

$\tau ={10}^{-6}$	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{L^2}$	EOC	${\parallel {𝐮}_h^J-{𝐮}_h^{\text{ref},J}\parallel }_{H^1}$	EOC
$h=2^{-2}$	$6.3361e-02$	$-$	$3.0327e-01$	$-$
$h=2^{-3}$	$1.5403e-02$	$2.04$	$1.0703e-01$	$1.50$
$h=2^{-4}$	$3.6554e-03$	$2.08$	$3.9342e-02$	$1.44$
$h=2^{-5}$	$9.2286e-04$	$1.99$	$1.3897e-02$	$1.50$
$h=2^{-6}$	$2.3430e-04$	$1.98$	$5.0386e-03$	$1.46$

Equations

$${\partial }_t\mathbf{u}=△\mathbf{u}+|\nabla \mathbf{u}{|}^2\mathbf{u},\mathbf{u}(0,\cdot )={\mathbf{u}}_0,\mathbf{u}(t,\cdot {)}_{|\partial \Omega }=({\mathbf{u}}_0{)}_{|\partial \Omega },t\in (0,T].$$

$$D$$

$$S^2$$

$$\begin{array}{rl}{\partial }_t\mathbf{u}& =△\mathbf{u}+|\nabla \mathbf{u}{|}^2\mathbf{u}\text{ on }D,t\in (0,T],\\ \mathbf{u}(0,x)& ={\mathbf{u}}_0(x)\text{ for }x\in D,\end{array}$$

$$\mathbf{u}(t,x)$$

$$E(\mathbf{u}):=\frac{1}{2}{\int }_D{\left|\nabla \mathbf{u}\right|}^{2}dx$$

$${\mathbf{u}}_0(r,\psi )=\left(\begin{array}{c}\cos \psi \sin u_0(r)\\ \sin \psi \sin u_0(r)\\ \cos u(r)\end{array}\right)$$

$$\mathbf{u}(t,r,\psi )=\left(\begin{array}{c}\cos \psi \sin u(t,r)\\ \sin \psi \sin u(t,r)\\ \cos u(t,r)\end{array}\right).$$

$$\begin{array}{rl}{\partial }_{t}u& ={\partial }_{rr}u+\frac{1}{r}{\partial }_{r}u-\frac{\sin (2u)}{2r^2}\text{ for }r\in (0,1),t\in (0,T],\\ u(0,r)& =u_0(r)\text{ for }r\in (0,1),\end{array}$$

$$u(t,0)=u_0(0),u(t,1)=u_0(1)\text{ for }t\in [0,T].$$

$$E(u)=\pi {\int }_0^{1}r({\partial }_{r}u{)}^2+\frac{{\sin }^2(u)}{r}dr.$$

$$\left({\partial }_{t}u,v\right)+\left({\partial }_{r}u,{\partial }_{r}v\right)-\left(\frac{1}{r}{\partial }_{r}u,v\right)+\left(\frac{\sin (2u)}{2r^2},v\right)=0,t\in (0,T],$$

$$\frac{\sin (2u)}{2r^2}=\left(\frac{\sin (2u)}{2u{r}^2}\right)u.$$

$$V_h^p$$

$$V_{h,u_0}^p$$

$${\partial }_{t}u(t_j,\cdot )$$

$$\begin{array}{rl}{\hat{u}}^{j,(1)}& :=u^{j-1}\text{(BDF1)},\\ {\hat{u}}^{j,(2)}& :=2u^{j-1}-u^{j-2}\text{(BDF2)}.\end{array}$$

$$\frac{{\delta }_0^{(k)}}{\tau }\left(u_h^{j+1},v_h\right)$$

$$+\left(\frac{\sin (2{\hat{u}}_h^{j+1,(k)})}{2r^2{\hat{u}}_h^{j+1,(k)}}{u}_h^{j+1},v_h\right)=-\frac{1}{\tau }\left(\sum _{i=1}^k{\delta }_i^{(k)}{u}_h^{j-i},v_h\right).$$

$$u_0(r)=\frac{1}{2}\pi r^2.$$

$$V_h^p$$

$${\mathbf{V}}_h^p$$

$${\mathbf{V}}_{h,{\mathbf{u}}_0}^p$$

$$F(u)(x,y)=\left(\begin{array}{c}x/r\sin u(r)\\ y/r\sin u(r)\\ \cos u(r)\end{array}\right),r=\sqrt{x^2+y^2},(x,y)\in D.$$

$${\mathbf{u}}_h^{\text{ref}}={\mathcal{I}}_h\left(F\left(u_h^{\text{ref}}\right)\right)$$

$$\left({\partial }_t\mathbf{u},\mathbf{v}\right)+\left(\nabla \mathbf{u},\nabla \mathbf{v}\right)-\left({\left|\nabla \mathbf{u}\right|}^2\mathbf{u},\mathbf{v}\right)=0,t\in (0,T],$$

$$|\mathbf{u}|=1\text{ a.e. in }[0,T]\times D,$$

$$\frac{{\delta }_0^{(k)}}{\tau }\left({\tilde{\mathbf{u}}}_h^{j+1},{\mathbf{v}}_h\right)+\left(\nabla {\tilde{\mathbf{u}}}_h^{j+1},\nabla {\mathbf{v}}_h\right)-\left({\left|\nabla {\hat{\mathbf{u}}}_h^{j+1,(k)}\right|}^2{\tilde{\mathbf{u}}}_h^{j+1},{\mathbf{v}}_h\right)$$

$$=-\frac{1}{\tau }\sum _{i=1}^k\left({\delta }_i^{(k)}{\mathbf{u}}_h^{j+1-i},{\mathbf{v}}_h\right),$$

$${\mathbf{u}}_h^{j+1}(z)=\frac{{\tilde{\mathbf{u}}}_h^{j+1}(z)}{\left|{\tilde{\mathbf{u}}}_h^{j+1}(z)\right|}\text{for all}z\in {\mathcal{N}}_h,$$

$${\mathbf{u}}_0(r,\psi )=\left(\begin{array}{c}\cos \psi \sin u_0(r)\\ \sin \psi \sin u_0(r)\\ \cos u_0(r)\end{array}\right),$$

$$\frac{1}{2}{\partial }_t\left({\left|\mathbf{u}\right|}^2\right)={\partial }_t\mathbf{u}\cdot \mathbf{u}$$

$$\mathbf{T}(\mathbf{u}):=\{\mathbf{v}\in L^2(D;{\mathbb{R}}^3)]:\mathbf{u}\cdot \mathbf{v}=0\text{ a.e.}\}$$

$$\left({\partial }_t\mathbf{u}(t,\cdot ),\mathbf{v}\right)+\left(\nabla \mathbf{u}(t,\cdot ),\nabla \mathbf{v}\right)=0,t\in (0,T],$$

Full text

1 Introduction

Let $\Omega\subset\mathbb{R}^{2}$ be a Lipschitz domain and $S^{2}$ the unit sphere in $\mathbb{R}^{3}$. In this paper, we study numerical methods for the discretization of the following harmonic map heat flow (HMHF) problem. Given an initial condition $\mathbf{u}_{0}:\Omega\to S^{2}$, determine $\mathbf{u}(t,\cdot):\Omega\to S^{2}$ such that

[TABLE]

This problem is obtained as the $L^{2}$ gradient flow of the Dirichlet energy $\int_{\Omega}|\nabla\mathbf{v}|^{2}\,dx$ on vector fields $\mathbf{v}:\Omega\to\mathbb{R}^{3}$ that satisfy a pointwise unit length constraint. Unit lenght minimizers of this Dirichlet energy are called harmonic maps. The HMHF problem is also closely related to the Landau-Lifshitz-Gilbert (LLG) equation. The HMHF equation can be considered as the limit of the LLG equation where the precessional term vanishes and only damping is left [26]. Harmonic maps, HMHF and LLG equations have numerous practical applications, for example, in the modeling of ferromagnetic materials or of liquid crystals, cf. e.g. [23, 25, 30, 21].

There is an extensive mathematical literature in which topics related to well-posedness, weak formulations, regularity, blow-up phenomena and convergence of solutions of the HMHF problem to harmonic maps are studied, cf. e.g [31, 32, 16, 19, 22, 17, 14, 34].

In this paper, we focus on finite element discretization methods for the HMHF problem. Early work on the development and analysis of numerical methods for HMHF or LLG problems is found in [10, 9, 30, 11, 4, 5, 15]. In recent years, there has been a renewed interest in the numerical analysis of methods for this problem class [3, 8, 7, 1, 2, 28, 29].

The main difficulties in solving HMHF (or LLG) numerically are the non-linear term and the fact that the solution has to satisfy the unit length constraint. In the literature, three different finite element based approaches for treating the unit length constraint and the nonlinearity have been proposed. We briefly outline the basic ideas of these approaches and refer to Section 3 for more details.

In the first approach [30, Chapter 4], [10, 9, 6, 18], the unit length constraint is enforced weakly by an additional step of renormalizing the numerical solution at every mesh point. Prohl [30, Chapter 4] shows that the renormalization step can be interpreted as a penalizing term in the equation. This renormalization leads to a nonlinear system that has to be solved in each time step, cf.  [9]. In [30, Section 4.3.2] and [6, 18], the nonlinearity is treated by a simple extrapolation in time in a semi-implicit Euler scheme which results in a linear system in each time step.

In the second approach [4, 5, 3, 1, 2], the unit length constraint is reformulated as an orthogonality condition between the solution and its time derivative. This leads to a solution space tangential to the sphere and the finite element space used in the discretization mimics this. In this approach one uses semi-implict BDF time stepping schemes in which a straightforward extrapolation is used to treat the nonlinearity.

In the third approach [11], the double cross product reformulation of the HMHF (or LLG) equation is used, cf. (32) below. This reformulation has the nice property that the unit length constraint is implicit in the formulation, meaning that its solution preserves the unit length of the initial condition. In the finite element spatial discretization, the norm preserving structure is mirrored. The reformulation, however, introduces a nonlinearity in the equation that is much more severe than the nonlinearity in the original equation. This much stronger nonlinearity is treated by fixed point iteration in each time step. This approach has also been applied to the $p$-harmonic map heat flow [12].

Although for all three approaches one finds papers in which results of numerical experiments with methods based on these are presented, we are not aware of any literature in which a systematic comparison of these three very different techniques is presented.

The main contribution of this paper is twofold. Firstly, we perform a systematic numerical study and compare convergence and efficiency properties of the above-mentioned methods. Secondly, the paper has review character in the sense that we present these methods in an unified framework. For the comparative study we restrict to the smooth regime, i.e. we do not consider HMHF with finite time blow up. In the smooth regime, for a given (standard) finite element space and consistency order of the time stepping scheme it is known what the optimal rates of convergence (w.r.t. mesh size and time step size) of the different methods are, which then allows a fair comparison of these methods. In the blow up case, singularities develop and a fair comparison of computational efficiencies becomes much more difficult. For the HMHF, exact solutions are not known, except for very special cases. Hence, for the method comparison we have to determine reliable and highly accurate reference solutions that can be used as sufficiently accurate approximations of the exact solution. We therefore restrict to a class of initial data in (1), for which global in time smoothness is guaranteed and, due to rotational symmetry, the problem can be transformed to a spatially one-dimensional problem. Using the latter property we can determine very accurate reference solutions by numerically solving this lower dimensional problem.

As a further contribution of this paper we present a discrete inf-sup stability result for the second approach, cf. Lemma 2 below. This result is relevant for an efficient implementation of the method.

The rest of the paper is organized as follows. In Section 2, we describe the harmonic map heat flow problem and its formulation for a radially symmetric case. Section 3 describes the finite element discretizations for this radially symmetric case and for the general harmonic map heat flow problem in the respective subsections. In each subsection, we also present numerical results to validate the convergence rates. The source code used in these experiments can be found in [27] and is based on the software package Netgen/NGSolve. At the end of Section 3 we draw conclusions and discuss computational efficiency aspects.

2 The harmonic map heat flow problem

Let

[TABLE]

be the open unit disk in $\mathbb{R}^{2}$ and the unit sphere in $\mathbb{R}^{3}$, respectively. Then, given $\mathbf{u}_{0}\in C^{2}(\bar{D},S^{2})$ and some $T>0$, we consider the following initial-boundary value problem

[TABLE]

which is the harmonic map heat flow problem. We prescribe the Dirichlet boundary condition

[TABLE]

The PDE (1) is the $L^{2}$-gradient flow of the energy

[TABLE]

for vector fields satisfying a pointwise unit length condition. Stationary solutions of (1) are called harmonic maps and are critical points of (3).

In this paper we study three different finite element discretization methods for (1) in a regime where we have a global in time smooth solution. For computing discretization errors we need a sufficiently accurate approximation of the exact solution of (1). This motivates why we restrict to a class of initial data in (1), for which global in time smoothness is guaranteed and, due to rotational symmetry, the problem can be transformed to a spatially one-dimensional problem. Due to the latter property we can determine very accurate reference solutions by numerically solving this lower dimensional problem, cf. Section 3.1 below. This class of smooth rotationally symmetric solutions is obtained from the following fundamental result.

Theorem 1** (Chang-Ding [13]).**

Using polar coordinates $(r,\psi)$ on the disk $D$, let the initial condition have the form

[TABLE]

where $u_{0}\in C^{2}([0,1])$ with $\|u_{0}\|_{L^{\infty}}\leq\pi$ is given, then there exists a unique global smooth solution to (1) of the form

[TABLE]

The spherically or radially symmetric form of the solution (5) reduces the harmonic map heat flow (1) to the one-dimensional problem

[TABLE]

with $u_{0}(0)=0$ where we prescribe the Dirichlet boundary conditions

[TABLE]

We will call (6) the radially symmetric harmonic map heat flow (RSHMHF). The energy of solutions of (6) is then given by [33, Equation (1.4)]

[TABLE]

3 Finite element discretization methods

In this section, we describe discretization methods for the HMHF problem. We apply the method of lines approach in which for spatial discretization we use finite element methods and for discretization in time a low order BDF (BDF1 or BDF2) method is used. In Section 3.1, we consider the RSHMHF (6). We propose a basic method for which optimal order convergence rates in $L^{2}$- and $H^{1}$-norms are demonstrated numerically. We use this method to construct highly accurate numerical reference solutions for the HMHF problem (1). These reference solutions are used to compute discretization errors for the three methods that are treated in Section 3.2.

In Section 3.2 we consider three very different finite element discretizations for HMHF (1) based on the schemes presented in [30, Chapter 4], [3] and [11]. These methods differ in how the unit length constraint of HMHF is treated. We will study convergence rates and efficiency aspect of these methods. The source code for the implementation of all the finite element methods can be found in [27].

3.1 Method for the radially symmetric case

We propose a very elementary method of lines discretization for the RSHMHF (6) using linear or quadratic finite elements in space and BDF1 or BDF2 for time discretization. Optimal order convergence of this method is demonstrated below. In [20], a finite difference based scheme with a moving mesh ansatz is developed for RSHMHF in the singular regime with blow-up of solutions in finite time. We are not aware of other literature in which discretization schemes for RSHMHF are systematically studied. We note that due to the strong nonlinearity of RSHMHF an error analysis of discretization methods for this problem is not straightforward, cf. the recent work [28, 29].

Let an initial condition $u_{0}\in C^{2}(\bar{I})$, $I:=(0,1)$, with $u_{0}(0)=0$ and $\|u_{0}\|_{L^{\infty}}\leq\pi$ be given. We introduce the solution space $H^{1}_{u_{0}}(I):=\{\,v\in H^{1}(I)~|~v(0)=0,~v(1)=u_{0}(1)\,\}$.

A variational formulation of (6)-(7) is as follows: find u\in C^{1}\big{(}[0,T];H^{1}_{u_{0}}(I)\big{)} with $u(0,\cdot)=u_{0}$ such that for all $v\in H^{1}_{0}(I)$

[TABLE]

where $(\cdot,\cdot)$ denotes the $L^{2}$-scalar product on $I$. Below in the time discretization, we use a linearization of the nonlinear term in (9) based on

[TABLE]

For simplicity, we use an equidistant partitioning of $\bar{I}$ into $N$ subintervals of equal length $h=1/N$, denoted by $\mathcal{T}_{h}$. We introduce the finite element spaces of continuous and piecewise polynomial functions of degree less than or equal to $p$:

[TABLE]

In the experiments below we restrict to $p=1$ and $p=2$.

Let $\tau$ be a constant time step and $J\in\mathbb{N}$ such that $J\tau=T$. For BDF1 (implicit Euler) and BDF2, we use the standard compact representation

[TABLE]

where the superscript $k\in\{1,2\}$ in $\delta_{i}^{(k)}$ corresponds to BDF$k$. The extrapolation used in the linearization is given by

[TABLE]

Let $\mathcal{I}_{h}$ be the nodal interpolation operator from $C^{0}(\bar{I},\mathbb{R})$ to $V_{h}^{p}$. We fix $k\in\{1,2\}$ and a timestep $\tau$. The discretization is as follows: Given $u_{h}^{0}=\mathcal{I}_{h}(u_{0})$ and $u_{h}^{1}\in V_{h,u_{0}}^{p}$ (only for $k=2$), for $j=k,\ldots,J-1$, find $u_{h}^{j+1}\in V_{h,u_{0}}^{p}$ such that for all $v_{h}\in V^{p}_{h,0}$

[TABLE]

For $k=2$ the initial value $u_{h}^{1}$ is determined by applying one step of the BDF1 scheme. An error analysis of the scheme (12) for $k=1$ is given in [29]. In [28], an error analysis of a finite difference scheme for (6) is presented.

Numerical results

We consider the RSHMHF for $t\in[0,10^{-1}]$ and with initial condition

[TABLE]

Note that this initial condition satisfies $\|u_{0}(r)\|_{L^{\infty}}<\pi$ for all $r\in[0,1]$ and thus Theorem 1 ensures that there is no blow-up and the exact solution is globally smooth. As a reference solution, we take the numerical solution $u_{h}^{\text{ref}}$ computed using (12) with $h=2^{-14}$, piecewise quadratic finite elements ($p=2$), $\tau=10^{-6}$ and BDF2 ($k=2$). We compute errors in $u_{h}^{J}\approx u(0.1,\cdot)$ in the $L^{2}$ and $H^{1}$ norms, i.e., $\lVert u_{h}^{J}-u_{h}^{\text{ref}}\rVert_{L^{2}(I)}$ and $\lVert u_{h}^{J}-u_{h}^{\text{ref}}\rVert_{H^{1}(I)}$. Results are presented in the Tables 1–4.

The results in Tables 1 and 2 show that the scheme (12) has optimal order of convergence in the time step $\tau$ for BDF1 and BDF2.

In Tables 3 and 4 we observe optimal order convergence rates with respect to $h$, in $L^{2}$- and $H^{1}$-norms, for $p=1$ and $p=2$. Below the scheme (12) is used to determine a highly accurate reference solution for the HMHF problem (1), cf. Remark 1.

3.2 Methods for the harmonic map heat flow

In this section, we study discretization methods for the HMHF problem (1). They are based on the works of Prohl [30, Section 4.3.2], Akrivis *et al. *[3] and Bartels-Prohl [11] and differ in how the unit length constraint of HMHF is treated. In [30, Section 4.3.2], the unit length constraint is weakly enforced by normalizing the numerical solution at every mesh point after each time step. A different way, proposed in [11], is to use the double cross product reformulation of HMHF, which automatically preserves the unit length of the initial condition, without using an additional unit length constraint. The method studied in [3] is based on the approach introduced in [5, 4] and exploits the fact that preserving the unit length constraint can be formulated as the time derivative of the solution being orthogonal to the solution itself. Motivated by this, one constructs an approximate tangent space to the sphere as the solution space for the time derivative to enforce the unit length. This idea of solving for the approximate time derivative in the finite element discretization goes back to [5, 4] where a pointwise orthogonality in the mesh points is used to define a discrete tangent space. In [3], an $L^{2}$-projection of the orthogonality condition into the finite element space is used.

Since the domain $D$ has a curved boundary, for the case of a quadratic finite element space $V_{h}^{2}$ we will use isoparametric elements, which is essential for obtaining optimal order of convergence, cf. [24]. We denote by ${\mathcal{T}}_{h}$ a regular and quasi-uniform triangulation of ${D}$ with maximum mesh parameter $h$. Note that ${\mathcal{T}}_{h}$ is not related to the interval partitioning (also denoted by ${\mathcal{T}}_{h}$) used in Section 3.1 above. The vertices of the boundary triangles are assumed to lie on $\partial D$. The set of vertices in $\mathcal{T}_{h}$ is denoted by $\mathbf{\mathcal{N}}_{{h}}$. We introduce finite element spaces of continuous piecewise polynomial functions of maximal degree $p$:

[TABLE]

In the numerical experiments, we restrict to $p=1$ and $p=2$. For $p=2$ the isoparametric variants of these spaces are used.

Remark 1**.**

We define the mapping $F:C(I)\rightarrow C(D;\mathbb{R}^{3})$ by

[TABLE]

This mapping is used to map a reference solution $u_{h}^{\text{ref}}$, cf. Section 3.1, to a corresponding $\mathbf{u}^{\text{ref}}_{h}$ for HMHF as follows. Given $u_{h}^{\text{ref}}$ at a fixed time $t$ we compute

[TABLE]

where $\mathcal{I}_{h}$ is the nodal interpolation operator from $C^{0}(\bar{D};\mathbb{R}^{3})$ to $\mathbf{V}_{h}^{p}$.

3.2.1 Treatment of unit length constraint by pointwise projection

We recall the method from [30, Section 4.3.2]. A variational formulation of (1) reads as follows: Given $\mathbf{u}_{0}\in H^{1}({D};\mathbb{R}^{3})$ with $\left|\mathbf{u}_{0}\right|=1$ (a.e), find \mathbf{u}\in C^{1}\big{(}[0,T];H^{1}({D};\mathbb{R}^{3})\big{)} with $\mathbf{u}(0,\cdot)=\mathbf{u}_{0}$ and $\mathbf{u}_{|\partial D}=(\mathbf{u}_{0})_{|\partial D}$ such that for all $\mathbf{v}\in H^{1}_{0}(D;\mathbb{R}^{3})$

[TABLE]

where $(\cdot,\cdot)$ denotes the $L^{2}$-scalar product on $D$. Let the time step $\tau$ and $J\in\mathbb{N}$ be such that $J\tau=T$. Based on this variational formulation we obtain the following discrete problem, cf. [30, Section 4.3.2], with given $k\in\{1,2\}$ and time step $\tau$: Given $\mathbf{u}_{h}^{0}=\mathcal{I}_{h}(\mathbf{u}_{0})$ and $\mathbf{u}_{h}^{1}\in\mathbf{V}_{h,\mathbf{u}_{0}}^{p}$ (only for $k=2$), for $j=k-1,\ldots,J-1$, find ${\mathbf{u}}_{h}^{j+1}\in\mathbf{{V}}^{p}_{{h},\mathbf{u}_{0}}$ such that for all $\mathbf{v}_{h}\in\mathbf{V}^{p}_{{h},0}$:

[TABLE]

where $\hat{\mathbf{u}}^{j+1,(k)}_{h}$ is the extrapolation given by the formula (11) and $\delta_{i}^{(k)}$ are the coefficients of the time stepping scheme (10). For $k=2$, the initial function $\mathbf{u}_{h}^{1}$ is determined by applying one step of the BDF1 ($k=1$) variant of the scheme (18)-(19). Note that the norm constraint is enforced by an additional pointwise normalization step. We will call the scheme (18)-(19) the pointwise projection finite element method (PPFEM) for the harmonic map heat flow.

Remark 2**.**

Under the assumption that the solution is sufficiently smooth and a parameter constraint $\tau^{-1/2}=o(h^{-1})$ is satisfied, it is shown in [30, Theorem 4.11] that a variant of this scheme, modified for the Landau-Lifschitz-Gilbert equation, converges with first order in $\tau$ for BDF1 $(k=1)$. The author also proves for linear finite elements first order convergence in $h$ in the $H^{1}$-norm. Recently, in [18] a slight modification of PPFEM is proposed, using a lumped mass finite element method. For this lumped mass PPFEM it is proved that on rectangular meshes and under the assumption $h^{p+1}=\mathcal{O}(\tau)$, one obtains optimal convergence order in $h$ and in $\tau$ in the $L^{2}$-norm, using BDF1 for time discretization and linear or higher order conforming finite element spaces. In the case of triangular meshes, one order in space is lost. The paper [18] also provides numerical results for the convergence behavior using a two dimensional rectangular domain (hence, no spherical symmetry). The lumped mass technique is also used in the method treated in [11], cf. Section 3.2.3 below. In [6], the discretization method described here is analyzed for LLG and they derive, for $p\geq 2$, optimal convergence order $\mathcal{O}(\tau+h^{p+1})$ in the $L^{2}$-norm under the time step condition $\tau=\mathcal{O}(h)$.

Numerical results

We consider the HMHF problem (1) for $t\in[0,10^{-1}]$ and with initial condition

[TABLE]

where $(r,\psi)$ are the polar coordinates on the disk and $u_{0}$ as given by (13). The numerical solution $\mathbf{u}_{h}^{j},j=0,\dots,J$, is determined by PPFEM (18)-(19). For computing errors we use the reference solution $u_{h}^{\text{ref}}$ computed using the method explained in Section 3.1 (with $h=2^{-14}$, $V_{h}^{2}$, BDF2 with $\tau=10^{-6}$). The reference solution $\mathbf{u}_{h}^{\text{ref}}$ is constructed from $u_{h}^{\text{ref}}$ by (17). We compute the errors in $\mathbf{u}_{h}^{J}\approx\mathbf{u}(0.1,\cdot)$ in the $L^{2}$ and $H^{1}$ norms, i.e., $\lVert\mathbf{u}_{h}^{J}-\mathbf{u}_{h}^{\text{ref}}\rVert_{L^{2}}$ and $\lVert\mathbf{u}_{h}^{J}-\mathbf{u}_{h}^{\text{ref}}\rVert_{H^{1}}$.

Results are presented in the Tables 5–8. In the Tables  5–6 we consider quadratic finite elements ($V_{h}^{2}$) with a fixed sufficiently small $h$ and vary $\tau$ to study the accuracy of the time discretization scheme.

Table 5 shows that for a small enough time step we reach a rate of convergence with optimal order 1 for BDF1 with respect to the time step $\tau$ in $L^{2}$ and $H^{1}$ norms. The results in Table 6 show the optimal convergence rate of order 2 for BDF2 with respect to $\tau$ in $L^{2}$ and $H^{1}$ norms.

In the Tables  7–8 we consider BDF2 with a fixed sufficiently small $\tau$ and vary $h$ to study the accuracy of the finite element discretization.

Table 7 shows optimal order of convergence in the $L^{2}$ norm with respect to $h$, while convergence rate in the $H^{1}$ norm is slightly above the optimal rate of 1 for the piecewise linear finite element space. Table 8 shows a significantly higher convergence rate than expected for quadratic finite elements. Furthermore, the error reductions in the $H^{1}$-norm and in the $L^{2}$-norm are approximately the same. Repeating the numerical experiment with a different initial condition $u_{0}(r)=0.5\pi\left(\sin(2\pi r)+r\right)$ we observe a very similar behavior in the error reductions. We have no satisfactory explanation for these observations.

3.2.2 Treatment of unit length constraint using a tangent space approach

From

[TABLE]

it follows that if the initial condition $\mathbf{u}(0,\cdot)=\mathbf{u}_{0}$ satisfies $|\mathbf{u}_{0}|=1$ on $D$ then the solution $\mathbf{u}$ satisfies the unit length constraint $|\mathbf{u}(t,\cdot)|=1$ on $D$ for $t\in[0,T]$ iff $\partial_{t}\mathbf{u}\cdot\mathbf{u}=0$ on $[0,T]\times D$. In [4, 5, 3], this elementary observation is used to construct a discretization method for the Landau-Lifshitz-Gilbert equation. In [4, 5], a discrete tangent space is constructed where the orthogonality is satisfied at each mesh point. Here, we apply the method studied in [3] to the harmonic map heat flow problem, where an $L^{2}$ projection of the orthogonality condition onto the finite element space is used, cf. (22).

Let

[TABLE]

be the tangent space at some $\mathbf{u}\in H^{1}({D},\mathbb{R}^{3})$. Note that if $\mathbf{u}$ is the solution of (1), then for $\mathbf{v}\in\mathbf{T}(\mathbf{u})$ we have $|\nabla\mathbf{u}|^{2}\mathbf{u}\cdot\mathbf{v}=0$. Using this and the elementary observation above leads to the following variational formulation. Given $\mathbf{u}_{0}\in H^{1}({D};\mathbb{R}^{3})$ with $\left|\mathbf{u}_{0}\right|=1$ on $D$, find \mathbf{u}\in C^{1}\big{(}[0,T];H^{1}({D},\mathbb{R}^{3})\big{)} such that $\mathbf{u}(0,\cdot)=\mathbf{u}_{0}$ and a time derivative $\partial_{t}\mathbf{u}(t,\cdot)\in\mathbf{T}(\mathbf{u}(t,\cdot))$ such that $(\partial_{t}\mathbf{u})_{|\partial D}=0$ and satisfies for all $\mathbf{v}\in\mathbf{T}(\mathbf{u}(t,\cdot))\cap H^{1}_{0}({D},\mathbb{R}^{3})$

[TABLE]

where $(\cdot,\cdot)$ is again the $L^{2}$ inner product on $D$. In the discretization, we use a discretized tangent space where the orthogonality is reduced to orthogonality in the finite element space as follows:

[TABLE]

We can write the approximation of the time derivative $\dot{\mathbf{u}}_{h}^{j}\approx\partial_{t}\mathbf{u}(t_{j},\cdot)$ by using (10) for a given $k\in\{1,2\}$ as

[TABLE]

Since the tangent space depends itself on the solution of the current time step, we extrapolate the solution from previous time steps to construct the discrete tangent space. The extrapolation is given by, cf. (11),

[TABLE]

Thus, we obtain the following scheme, cf. [3, Equation (2.6)]. Let the time step $\tau$ and $J\in\mathbb{N}$ be such that $J\tau=T$. Given $\tau>0$ and $p,k\in\{1,2\}$, $\mathbf{u}_{h}^{0}=\mathcal{I}_{h}(\mathbf{u}_{0})$ and $\mathbf{u}_{h}^{1}\in\mathbf{V}_{h,\mathbf{u}_{0}}^{p}$ (only for $k=2$), for $j=k-1,\ldots,J-1$, find the approximate time derivative $\dot{\mathbf{u}}_{h}^{j+1}\in\mathbf{T}^{p}_{h,0}\left(\hat{\mathbf{u}}_{h}^{j+1,(k)}\right)$ such that for all $\mathbf{v}_{h}\in\mathbf{T}^{p}_{{h,0}}\left(\hat{\mathbf{u}}_{h}^{j+1,(k)}\right)$

[TABLE]

and the update step is given by

[TABLE]

which follows directly from (23). For $k=2$ the initial function $\mathbf{u}_{h}^{1}$ is determined by applying one step of the BDF1 variant of the scheme. The scheme (25)-(26) is analyzed in [3].

Remark 3**.**

In [3, Theorem 3.1], for a variant of this scheme, modified for the LLG equation with Neumann boundary condition, an error bound as in (31) is derived (under a mild CFL-type condition) for BDF$k$ with $k\in\{1,2\}$. For BDF$k$ with $3\leq k\leq 5$, optimal $H^{1}$ norm error bounds of order $\tau^{k}+h^{p}$ are derived for $p\geq 2$, using a much stronger CFL type condition $\tau=\mathcal{O}(h)$ and with an additional assumption concerning the size of the damping parameter in the LLG equation, cf.  [3, Theorem 3.2]. The paper also contains results of numerical results with exact solutions for the LLG equation on a three-dimensional cube with a thin layer in $z$-direction to confirm the theoretical discretization error bounds.

Concerning implementation of this scheme we note the following. It is inconvenient and in general computationally relatively expensive to construct a basis of the finite element space $\mathbf{T}_{h,0}^{p}(\mathbf{u})$. This can be circumvented using a Lagrange multiplier approach, as proposed in [3]. We introduce a Lagrange multiplier $\lambda_{h}\in V_{h,0}^{p}$ to treat the constraint $\mathbf{u}_{h}\cdot\mathbf{v}_{h}=0$ in (22). This leads to the following saddle point problem: Given a time step $\tau>0$ (with $\tau J=T$) and $p,k\in\{1,2\}$, $\mathbf{u}_{h}^{0}=\mathcal{I}_{h}(\mathbf{u}_{0})$, $\mathbf{u}_{h}^{1}\in\mathbf{V}_{h,\mathbf{u}_{0}}^{p}$ (only for $k=2$), for $j=k-1,\ldots,J-1$, find $\left(\dot{\mathbf{u}}_{h}^{j+1},\lambda_{h}\right)\in\left(\mathbf{{V}}_{h,0}^{p},{V}_{h,0}^{p}\right)$ such that for all $\left(\mathbf{v}_{h},w_{h}\right)\in\left(\mathbf{{V}}_{h,0}^{p},{V}_{h,0}^{p}\right)$

[TABLE]

Given $\dot{\mathbf{u}}_{h}^{j+1}$ we determine $\mathbf{u}_{h}^{j+1}$ as in (26). We call (27)-(28) together with (26) the tangential finite element method (TFEM) for HMHF.

Concerning well-posedness of (27)-(28) and the relation between this scheme and (25) we note the following. Let $\left(\dot{\mathbf{u}}_{h}^{j+1},\lambda_{h}\right)\in\left(\mathbf{{V}}_{h,0}^{p},{V}_{h,0}^{p}\right)$be a solution of (27)-(28); then (28) implies $\dot{\mathbf{u}}_{h}^{j+1}\in\mathbf{T}^{p}_{h,0}\left(\hat{\mathbf{u}}_{h}^{j+1,(k)}\right)$ and, restricting in (27) to $\mathbf{v}_{h}\in\mathbf{T}^{p}_{h,0}\left(\hat{\mathbf{u}}_{h}^{j+1,(k)}\right)$, then yields that $\dot{\mathbf{u}}_{h}^{j+1}$ solves (25). Hence, if the discrete saddle point problem (27)-(28) has a unique solution, then its $\mathbf{u}$-solution is the same as that of (25).

In [3] the well-posedness of (27)-(28) is not addressed. In the subsection below we show that under reasonable assumptions the saddle point formulation is well-posed.

Well-posedness of the saddle point formulation

The bilinear form $(\mathbf{w}_{h},\mathbf{v}_{h})\to\left(\mathbf{w}_{h},\mathbf{v}_{h}\right)+\frac{\tau}{\delta_{0}^{(k)}}\left(\nabla\mathbf{w}_{h}^{j+1},\nabla\mathbf{v}_{h}\right)$ used in (27) is elliptic on $\mathbf{{V}}_{h,0}^{p}$. Hence, we have well-posedness of (27)-(28) iff the bilinear form $\mathbf{b}(\hat{\mathbf{u}}_{h}^{j+1,(k)};\cdot,\cdot)$ has the discrete inf-sup property. We have the following result:

Lemma 2**.**

Assume that the solution $\mathbf{u}$ of (1) is sufficiently smooth and for the discrete approximation $(\mathbf{u}_{h}^{j})_{0\leq j\leq J}$ of (25)-(26) we have uniform convergence $\mathbf{u}_{h}^{j}\to\mathbf{u}(t_{j},\cdot)$ in the following sense: For arbitrary $\varepsilon\in(0,1]$ there exist $\varepsilon_{1}>0,\varepsilon_{2}>0$ (depending on $\mathbf{u}$) such that

[TABLE]

Then there exists $\beta=\beta(\mathbf{u})>0$, independent of $h$ and $\tau$, such that for $k\in\{0,1\}$ and $\tau$, $h$ sufficiently small we have:

[TABLE]

Proof.

A proof is given in the Appendix A. ∎

The result (30) yields well-posedness of (27)-(28), provided assumption (29) is satisfied and $\tau$ and $h$ are sufficiently small. Moreover, the discrete inf-sup constant $\beta$ in (30) is independent of the discretization parameters $h$ and $\tau$, which implies a control of the condition number of the Schur complement of the symmetric indefinite system matrix corresponding to the saddle point problem (27)-(28).

Remark 4**.**

Concerning the assumption (29) we note the following. In [3, Theorem 3.1], for a variant of the scheme (25)-(26) applied to the LLG equation, an optimal error bound of the form

[TABLE]

is derived under the very mild CFL-type condition $\tau^{k}\leq\bar{c}h^{\frac{1}{2}}$ for a sufficiently small constant $\bar{c}$. Assume that (31) also holds for the scheme applied to the HMHF problem and assume that a slightly stronger mesh-parameter condition is satisfied, namely $\tau^{k}\leq ch^{\frac{1}{2}+\epsilon}$ holds for some constants $c>0$ and $\epsilon>0$. Then we obtain, with $I_{h}$ the nodal interpolation in the finite element space $\mathbf{V}_{h,0}^{p}$:

[TABLE]

which implies that the assumption (29) is satisfied.

Numerical results

A reference solution $\mathbf{u}_{h}^{\text{ref}}$ is determined in the same way as in Section 3.2.1. The discrete approximations $\mathbf{u}_{h}^{j}$ are computed using TFEM (27)-(28) and (26). Results are shown in Tables 9–12.

Similar to PPFEM, Table 9 shows for BDF1 an optimal convergence rate of order 1 in $\tau$, both in the $L^{2}$- and $H^{1}$-norm. In Table 10 we observe the optimal convergence rate of order 2 for BDF2 with respect to $\tau$ in $L^{2}$ and slightly less than 2 in the $H^{1}$-norm.

Table 11 shows an optimal convergence rate of order 2 in the $L^{2}$ norm with respect to $h$, while the convergence rate in the $H^{1}$-norm is slightly above the optimal rate of 1, again similar to PPFEM, cf. Table 7. Table 12 also shows an optimal convergence rate of order 3 in the $L^{2}$ norm with respect to $h$, while the convergence rate in the $H^{1}$-norm is slightly above the optimal rate of 2.

Summarizing, the results for PPFEM and TFEM are very similar, except for BDF2 combined with $\mathbf{V}_{h}^{2}$ in the Tables 8 and 12.

3.2.3 Treatment of unit length constraint based on a constraint preserving formulation

In view of the vector identity

[TABLE]

we can write

[TABLE]

Because of $\mathbf{u}\in\mathbb{S}^{2}$, we have $|\mathbf{u}|=1$ and $\mathbf{u}\cdot\nabla\mathbf{u}=0$ and obtain

[TABLE]

Thus, we can rewrite (1) as

[TABLE]

Note that multiplying both sides of (32) with $\mathbf{u}$ results in $\partial_{t}\mathbf{u}\cdot\mathbf{u}=0$, and thus

[TABLE]

Since $|\mathbf{u}_{0}|=1$ on $D$ holds, this implies that a solution of (32) preserves the unit length property $\left|\mathbf{u}\right|=1$ on $D$. The method presented in [11] is based on the double cross formulation of the harmonic map heat flow (32). A variational formulation reads: Given $\mathbf{u}_{0}\in H^{1}({D};\mathbb{R}^{3})$ with $\left|\mathbf{u}_{0}\right|=1$ (a.e.), find $\mathbf{u}\in C^{1}([0,T];H^{1}({D};\mathbb{R}^{3})$ with $\mathbf{u}(0,\cdot)=\mathbf{u}_{0}$ and $\mathbf{u}_{|\partial D}=(\mathbf{u}_{0})_{|\partial D}$ such that for all $\mathbf{v}\in H^{1}_{0}({D};\mathbb{R}^{3})$

[TABLE]

Let $\{\mathbf{\phi}_{z}:z\in\mathbf{\mathcal{N}}_{h}\}$ be the nodal basis of the linear finite element space $\mathbf{V}_{h}^{1}$. For $\mathbf{u},\mathbf{v}\in C^{0}(\bar{D};\mathbb{R}^{3})$, we define

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the Euclidean scalar product in $\mathbb{R}^{3}$. The corresponding seminorm is denoted by $\lVert\mathbf{u}\rVert_{h}:=(\mathbf{u},\mathbf{u})_{h}^{\frac{1}{2}}$. On $\mathbf{V}_{h}^{1}$ this defines a norm that is equivalent to the $L^{2}$-norm [11, Equation (2.3)]. The discrete Laplacian $\tilde{\mathop{}\!\mathbin{\bigtriangleup}}_{h}:H^{1}(D;\mathbb{R}^{3})\rightarrow\mathbf{V}_{h}^{1}$ is given by

[TABLE]

We define

[TABLE]

Following [11], we use BDF1 for time discretization. For a given time step $\tau$ (with $\tau J=T$) the resulting discretization reads: Given $\mathbf{u}_{h}^{0}=\mathcal{I}_{h}(\mathbf{u}_{0})$, for $j=0,\ldots,J-1$, find $\mathbf{u}_{h}^{j+1}\in\mathbf{V}_{h,\mathbf{u}_{0}}^{1}$ such that for all $\mathbf{v}_{h}\in\mathbf{V}^{1}_{h,0}$

[TABLE]

Testing (34) with $\mathbf{v}_{h}=\bar{\mathbf{u}}_{h}^{j+1/2}(z)\mathbf{\phi}_{z}$ for $z\in{\mathbf{\mathcal{N}}_{h}}$ results in

[TABLE]

which is a discrete analogon of the length preservation property (33). It implies that the constraint is pointwise satisfied at every mesh point $z\in{\mathbf{\mathcal{N}}_{h}}$ provided the initial condition satisfies the norm constraint. Due to the cross product the discrete problem (34) for the unknown $\mathbf{u}_{h}^{j+1}$ is strongly nonlinear. In [11], the following iterative linearization method is proposed. We introduce $\mathbf{w}^{j+1}_{h}:=\bar{\mathbf{u}}_{h}^{j+1/2}$ and rewrite (34) as

[TABLE]

for all $\mathbf{v}_{h}\in\mathbf{V}^{1}_{h,0}$. To this nonlinear problem a fixed point iteration (Algorithm 1) is applied, cf. [11, Algorithm 4.1].

Remark 5**.**

In [11], for the case of a Neumann boundary condition, a result on (unconditional) weak convergence (for $h\to 0$) of the solution of (34) to the (weak) solution of the harmonic map heat flow (1) is derived. No result for the rate of convergence rate is available. Furthermore, in [11, Theorem 4.1], it is proved that Algorithm 1 converges provided the strong CFL condition $\tau=\mathcal{O}(h^{2})$ is satisfied. No rigorous results on rate of convergence of Algorithm 1 are known.

Numerical results

A reference solution $\mathbf{u}_{h}^{\text{ref}}$ is determined in the same way as in Section 3.2.1. The discrete approximations $\mathbf{u}_{h}^{j}$ are computed using Algorithm 1.

Table 13 shows an optimal convergence rate of order 2 in the $L^{2}$ norm with respect to $h$, while the convergence rate in the $H^{1}$ norm is slightly higher than the optimal one. The results are similar to those of PPFEM, cf. Table 7 and of TFEM, cf. Table 11. We use a tolerance parameter $\epsilon=10^{-10}$; in this example the fixed point algorithm 1 then needs on average 3 iterations per time step to satisfy the tolerance criterion. Numerical experiments show that the CFL condition $\tau=\mathcal{O}(h^{2})$ that is needed in the analysis is essential. For example, the algorithm 1 does not convergence for our example is we use a scaling $\tau\sim h$.

Because of this CFL time step restriction $\tau=\mathcal{O}(h^{2})$, it is not possible to study the convergence rate depending on $\tau$ for a fixed $h$ that is small enough such that the spatial error does not dominate.

3.3 Comparison of methods

All three methods show similar convergence behavior with BDF1 und piecewise linear finite element spaces.

Although the nonlinear scheme (34) converges unconditionally, the proposed Algorithm 1 suffers from the severe time step restriction. Higher order BDF schemes do not satisfy the relation (35) and it is not clear how to extend the scheme to higher order BDF methods. Algorithm 1 is computationally the most expensive among the three methods since in each time step one needs to solve a fixed point iteration which typically takes some iterations to converge for a given suitable tolerance. In our comparative study Algorithm 1 turns out to be the least efficient one.

PPFEM ((18)-(19)) is the simplest one to implement. One has to solve a linear system of dimension equal to the dimension of the underlying finite element space in each time step. The renormalization at every mesh point afterwards is computationally negligible in comparison to solving the discrete system. Altogether, in each time step, this is significantly less computationally expensive compared to TFEM and Algorithm 1. PPFEM with BDF2 combined with piecewise quadratic finite elements shows a significantly higher rate of convergence in the mesh size than expected.

TFEM ((27)-(28) and (26)) performs as expected from the theoretical results in terms of convergence behavior in mesh size and time step. One drawback of TFEM is the higher computational time compared to PPFEM since due to the Lagrange multiplier used in the saddle point formulation the dimension of the linear equation, which has to be solved in each time step, is increased.

Acknowledgements The authors acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 442047500 – through the Collaborative Research Center “Sparsity and Singular Structures” (SFB 1481).

Appendix A Proof of Lemma 2

Take $w_{h}\in V_{h,0}^{p}$. The constants below (all denoted by $c$) may depend on $\mathbf{u}$, but are independent of $w_{h}$, $\tau$, $h$. Let $I_{h}:C(\bar{D};\mathbb{R}^{3})\to\mathbf{V}_{h,0}^{p}$ be the componentwise nodal interpolation in the finite element space. We write $\mathbf{u}(t_{j})=\mathbf{u}(t_{j},\cdot)$. For $K\in\mathcal{T}_{h}$ we obtain, using $|w_{h}|_{H^{p+1}(K)}=0$ and an inverse inequality:

[TABLE]

From this we obtain (with $c=c(\mathbf{u})$)

[TABLE]

We restrict to $\tau\in(0,\varepsilon_{1}]$, $h\in(0,\varepsilon_{2}]$ with suitable $\varepsilon_{1}>0$, $\varepsilon_{2}>0$. Define $\delta_{\tau,h}:=\max_{0\leq j\leq J}\|\mathbf{u}_{h}^{j}-\mathbf{u}(t_{j})\|_{L^{\infty}}$. Due to assumption (29) the size of $\delta_{\tau,h}$ is controlled by $\varepsilon_{1}$ and $\varepsilon_{2}$. Using $|\mathbf{u}(t_{n})|=1$ on $\Omega$ for all $0\leq n\leq J$, we obtain $\left||\mathbf{u}_{h}^{n}|-1\right|\leq\delta_{\tau,h}$ and $\left||2\mathbf{u}_{h}^{n}-\mathbf{u}_{h}^{n-1}|-1\right|\leq 3\delta_{\tau,h}+c\tau$. Hence, for the extrapolations defined in (24) and with $\varepsilon_{1},\varepsilon_{2}$ sufficiently small we get

[TABLE]

Since $w_{h}=0$ on the boundary of the triangulation we have $\mathbf{v}_{h}=I_{h}(\mathbf{u}(t_{j})w_{h})\in\mathbf{V}_{h,0}^{p}$. Using (36) and (37) we obtain

[TABLE]

For $\varepsilon_{1}$ and $\varepsilon_{2}$ sufficiently small we have $1-c(\delta_{\tau,h}+\tau+h)\geq\tfrac{1}{2}$. Finally note that $\|I_{h}(\mathbf{u}(t_{j})w_{h})\|_{L^{2}}\leq c\|w_{h}\|_{L^{2}}$ holds.