High Order Hierarchical Divergence-free Constrained Transport $H(div)$ Finite Element Method for Magnetic Induction Equation

TL;DR

This paper introduces a high-order hierarchical divergence-free finite element method for the magnetic induction equation, improving the enforcement of the divergence-free condition of magnetic fields in MHD simulations.

Contribution

It proposes a novel hierarchical basis approach using interior functions of BDM_p elements to enforce divergence-free magnetic fields in finite element solutions.

Findings

Effective enforcement of divergence-free condition in 3D magnetic induction simulations

High-order accuracy achieved with hierarchical basis functions

Numerical results demonstrate improved solution quality

Abstract

In this paper, we will use the interior functions of an hierarchical basis for high order $BDM_p$ elements to enforce the divergence-free condition of a magnetic field $B$ approximated by the H(div) $BDM_p$ basis. The resulting constrained finite element method can be used to solve magnetic induction equation in MHD equations. The proposed procedure is based on the fact that the scalar $(p-1)$-th order polynomial space on each element can be decomposed as an orthogonal sum of the subspace defined by the divergence of the interior functions of the $p$-th order $BDM_p$ basis and the constant function. Therefore, the interior functions can be used to remove element-wise all higher order terms except the constant in the divergence error of the finite element solution of $B$-field. The constant terms from each element can be then easily corrected using a first order H(div) basis globally. Numerical results for a 3-D magnetic induction equation show the effectiveness of the proposed method in enforcing divergence-free condition of the magnetic field.

Tables

Table 1. Table 2.1: Decomposition of the ( ℙ p ​ ( K ) ) 3 superscript subscript ℙ 𝑝 𝐾 3 \left(\mathbb{P}_{p}(K)\right)^{3} B ​ D ​ M p 𝐵 𝐷 subscript 𝑀 𝑝 BDM_{p} tetrahedral finite element space.

Decomposition	Dimension
Edge-based face functions	$12p$
Face bubble functions	$2\left(p-2\right)\left(p-1\right)$
Edge-based interior functions	$6\left(p-1\right)$
Face-based interior functions	$4\left(p-2\right)\left(p-1\right)$
Interior bubble functions	$\left(p-3\right)\left(p-2\right)\left(p-1\right)/2$
Total	$\left(p+1\right)\left(p+2\right)\left(p+3\right)/2=dim{\left(P_p\left(K\right)\right)}^3$

Table 2. Table 4.1: The errors and the order of convergence, by the P 3 subscript 𝑃 3 P_{3} element, for ( 4.1 ) and divergence-free corrections.

	${\Vert B-B_h\Vert }_0$	$h^n$	${\Vert \mathrm{div}\left(B-B_h\right)\Vert }_0={\Vert \mathrm{div}{B}_h\Vert }_0$	$h^n$	$dim{V}_3$
2	0.86157111	0.0	2.52192914	0.0	1920
3	0.31631008	1.4	1.27803825	1.0	15360
4	0.05191410	2.6	0.44154994	1.5	122880
5	0.00635164	3.0	0.07211612	2.6	983040
For $\mathrm{div}{B}_h=0$ corrected solution with both (3.13) and (3.17).
2	0.95693378	0.0	0.00000032	0.0	1920
3	0.32810096	1.5	0.00000012	1.4	15360
4	0.05713843	2.5	0.00000040	0.0	122880
5	0.00698271	3.0	0.00000063	0.0	983040

Table 3. Table 4.2: The errors and the order of convergence of P 3 subscript 𝑃 3 P_{3} solutions, for ( 4.3 ) with global P 3 − H ​ ( d ​ i ​ v ) subscript 𝑃 3 𝐻 𝑑 𝑖 𝑣 P_{3}-H(div) and local interior modes only corrections.

	${\Vert B-B_h\Vert }_0$	$h^n$	${\Vert \mathrm{div}\left(B-B_h\right)\Vert }_0={\Vert \mathrm{div}{B}_h\Vert }_0$	$h^n$	$dim{V}_3$
For a global $\mathrm{div}{B}_h=0$ corrected solution with $P_3-H\left(div\right)$ basis.
2	0.84092712	0.0	0.00000484	0.0	1920
3	0.31432784	1.4	0.00000476	0.0	15360
4	0.05155935	2.6	0.00000487	0.0	122880
5	0.00618254	3.0	0.00000447	0.0	983040
For local $\mathrm{div}{B}_h=0$ corrected solution only with (3.13).
2	0.93913487	0.0	1.57251299	0.0	1920
3	0.32826239	1.5	1.03345712	0.6	15360
4	0.05715300	2.5	0.32677858	1.7	122880
5	0.00665786	3.1	0.02944335	3.5	983040

Equations

$${\lambda }_0:=1-\xi -\eta -\zeta ,{\lambda }_1:=\xi ,{\lambda }_2:=\eta ,{\lambda }_3:=\zeta .$$

$${\tau }^{{\mathbf{e}}_j}:={\tau }^{[j_1,j_2]}={\mathbf{v}}_{j_2}-{\mathbf{v}}_{j_1},j_1<j_2.$$

$${\gamma }_{{\mathbf{e}}_j}:={\lambda }_{j_2}-{\lambda }_{j_1},j_1<j_2.$$

$${\mathbf{e}}_j:=[j_1,j_2],j_1=0,1,2,j_1<j_2\le 3,j=j_1+j_2+\mathrm{sign}(j_1),$$

$${\mathbf{f}}_{j_1}:=[j_2,j_3,j_4],0\le \{j_1,j_2,j_3,j_4\}\le 3,j_2<j_3<j_4.$$

$${\mathbf{n}}^{{\mathbf{f}}_{j_k}}\cdot {\Phi }_{\mathbf{e}[k_1,k_2]}^{{\mathbf{f}}_{j_1},i}=0,j_k\ne j_1,$$

$${\tilde{\Phi }}_{\mathbf{e}[k_1,k_2]}^{{\mathbf{f}}_{j_1},0}={\lambda }_{k_1}\nabla {\lambda }_{k_2}\times \nabla {\lambda }_{k_3},$$

$${\tilde{\Phi }}_{\mathbf{e}[k_1,k_2]}^{{\mathbf{f}}_{j_1},1}={\lambda }_{k_1}{\lambda }_{k_2}\nabla {\lambda }_{k_3}\times \nabla {\lambda }_{k_1},$$

$${\tilde{\Phi }}_{\mathbf{e}[k_1,k_2]}^{{\mathbf{f}}_{j_1},i+1}$$

$$={\ell }_i({\gamma }_{{\mathbf{e}}_k})\left[{\lambda }_{k_1}{\lambda }_{k_2}\nabla {\lambda }_{k_3}\times \nabla {\lambda }_{k_1}\right]+{\ell }_{i-1}({\gamma }_{{\mathbf{e}}_k})\left[{\lambda }_{k_1}\nabla {\lambda }_{k_2}\times \nabla {\lambda }_{k_3}\right],$$

$${\mathbf{n}}^{{\mathbf{f}}_{j_k}}\cdot {\Phi }_{m,n}^{{\mathbf{f}}_{j_1}}=0,j_k\ne j_1.$$

$${\Phi }_{m,n}^{{\mathbf{f}}_{j_1}}={\lambda }_{j_2}{\lambda }_{j_3}{\lambda }_{j_4}{L}_{m,n}\frac{\nabla {\lambda }_{j_3}\times \nabla {\lambda }_{j_4}}{|\nabla {\lambda }_{j_3}\times \nabla {\lambda }_{j_4}|},$$

$$L_{m,n}=(1-{\lambda }_{j_2}{)}^m(1-{\lambda }_{j_2}-{\lambda }_{j_3}{)}^n{P}_m^{(2n+3,2)}\left(\frac{2{\lambda }_{j_3}}{1-{\lambda }_{j_2}}-1\right)P_n^{(0,2)}\left(\frac{2{\lambda }_{j_4}}{1-{\lambda }_{j_2}-{\lambda }_{j_3}}-1\right)$$

$$0\le \{m,n\},m+n\le p-3.$$

$$<{\Phi }_{m_1,n_1}^{{\mathbf{f}}_{j_1}},{\Phi }_{m_2,n_2}^{{\mathbf{f}}_{j_1}}>{|}_{K^3}={\delta }_{m_1{m}_2}{\delta }_{n_1{n}_2},$$

$${\mathbf{n}}^{{\mathbf{f}}_j}\cdot {\Phi }^{\mathbf{t}}=0,j=\{0,1,2,3\}.$$

$${\tau }^{{\mathbf{e}}_j}\cdot {\Phi }_{\mathbf{e}[k_1,k_2]}^{\mathbf{t},i}=0,{\mathbf{e}}_j\ne {\mathbf{e}}_k,$$

$${\Phi }_{\mathbf{e}[k_1,k_2]}^{\mathbf{t},i}={\lambda }_{k_1}{\lambda }_{k_2}\left\{(1-{\lambda }_{k_1}{)}^i{P}_i^{(1,2)}\left(\frac{2{\lambda }_{k_2}}{1-{\lambda }_{k_1}}-1\right)\right\}\frac{{\tau }^{{\mathbf{e}}_k}}{|{\tau }^{{\mathbf{e}}_k}|},$$

$$<{\Phi }_{\mathbf{e}[k_1,k_2]}^{\mathbf{t},m},{\Phi }_{\mathbf{e}[k_1,k_2]}^{\mathbf{t},n}>{|}_{K^3}={\delta }_{mn},\{m,n\}=0,1,\cdots ,p-2.$$

$${\mathbf{n}}^{{\mathbf{f}}_{j_k}}\times {\Phi }_{m,n}^{\mathbf{t},{\mathbf{f}}_{j_1}}=0,j_k\ne j_1.$$

$${\tau }^{{\mathbf{e}}_k}\cdot {\Phi }_{m,n}^{\mathbf{t},{\mathbf{f}}_{j_1}}=0.$$

$${\Phi }_{m,n}^{\mathbf{t},{\mathbf{f}}_{j_1}^1}={\lambda }_{j_2}{\lambda }_{j_3}{\lambda }_{j_4}{L}_{mn}\frac{{\tau }^{[j_2,j_3]}}{\left|{\tau }^{[j_2,j_3]}\right|},{\Phi }_{m,n}^{\mathbf{t},{\mathbf{f}}_{j_1}^2}={\lambda }_{j_2}{\lambda }_{j_3}{\lambda }_{j_4}{L}_{mn}\frac{{\tau }^{[j_2,j_4]}}{\left|{\tau }^{[j_2,j_4]}\right|}$$

$$L_{mn}=(1-{\lambda }_{j_2}{)}^m(1-{\lambda }_{j_2}-{\lambda }_{j_3}{)}^n{P}_m^{(2n+3,2)}\left(\frac{2{\lambda }_{j_3}}{1-{\lambda }_{j_2}}-1\right)P_n^{(0,2)}\left(\frac{2{\lambda }_{j_4}}{1-{\lambda }_{j_2}-{\lambda }_{j_3}}-1\right)$$

$$<{\Phi }_{m_1,n_1}^{\mathbf{t},{\mathbf{f}}_{j_1}^i},{\Phi }_{m_2,n_2}^{\mathbf{t},{\mathbf{f}}_{j_1}^i}>{|}_{K^3}={\delta }_{m_1{m}_2}{\delta }_{n_1{n}_2},$$

$${\Phi }_{\ell ,m,n}^{\mathbf{t},{\vec{e}}_i}={\lambda }_0{\lambda }_1{\lambda }_2{\lambda }_3{L}_{lmn}{\vec{e}}_i,i=1,2,3,$$

$$L_{lmn}$$

$$\cdot P_m^{(2n+5,2)}\left(\frac{2{\lambda }_2}{1-{\lambda }_1}-1\right)P_n^{(2,2)}\left(\frac{2{\lambda }_3}{1-{\lambda }_1-{\lambda }_2}-1\right)$$

$$0\le \{\ell ,m,n\},\ell +m+n\le p-4.$$

$$<{\Phi }_{{\ell }_1,m_1,n_1}^{\mathbf{t},{\vec{e}}_i},{\Phi }_{{\ell }_2,m_2,n_2}^{\mathbf{t},{\vec{e}}_j}>{|}_{K^3}={\delta }_{{\ell }_1{\ell }_2}{\delta }_{m_1{m}_2}{\delta }_{n_1{n}_2},$$

$$0\le \{{\ell }_1,{\ell }_2,m_1,m_2,n_1,n_2\},{\ell }_1+m_1+n_1,{\ell }_2+m_2+n_2\le p-4,\{i,j\}=1,2,3.$$

$${\Sigma }_{\mathrm{int}}=\mathrm{span}\{{\Phi }_i{\}}_{i=1}^{n_i},n_i=\frac{1}{2}(p-1)(p+1)(p+2).$$

Full text

\ams

65M60,76W05

1 Introduction

Numerical modeling of magneto-hydrodynamic fluids has shown that the observance of the zero divergence of the magnetic field plays an important role in reproducing the correct physics in plasmas [3]. Various numerical techniques have been devised to ensure the computed magnetic field to be divergence-free [6]. In the early work of [3] a projection approach was used to correct the magnetic field to have a zero divergence. A more natural way to satisfy this constraint is through a class of the so-called constrained transport (CT) numerical methods based on the ideas in [5]. In most CT algorithms for the MHD, the surface averaged magnetic flux over the surface of a 3-D element is used to represent the magnetic field so normal continuity of the magnetic field can be assured while the volume averaged conserved quantities are used for mass, momentum, and energy variables.

In this paper, we will propose a high order transport finite element method using a recently developed high order hierarchical basis for the $BDM_{p}$ element [4] for the magnetic induction equation in the MHD problems. The divergence condition on the $B$ field is enforced through corrections with interior functions in the basis set such that the global divergence-free condition will be satisfied.

The paper is organized as follows. In section 2, we will present the hierarchical H(div) basis functions in various modes (edge, face, and interior). In section 3, we will first characterize the divergence of the interior basis functions for the hierarchical H(div) basis, then, we will introduce a two-step procedure to remove non-zero divergence in the finite element solution. Numerical test of the proposed procedure will be carried out for a 3-D magnetic induction equation in Section 4. Finally, a conclusion is given in Section 5.

2 Basis functions for the tetrahedral element

In this section we present hierarchical shape functions proposed in [2] for the  H(div)-conforming tetrahedral $BDM_{p}$ element on the canonical reference 3-simplex . The shape functions are grouped into several categories based upon their geometrical entities on the reference 3-simplex [1]. The basis functions in each category are constructed so that they are also orthonormal within each category on the reference element.

Any point in the 3-simplex $K^{3}$ is uniquely located in terms of the local coordinate system $(\xi,\eta,\zeta)$. The vertexes are numbered as $\mathbf{v}_{0}(0,0,0),$ $\mathbf{v}_{1}(1,0,0),$ $\mathbf{v}_{2}(0,1,0),$ $\mathbf{v}_{3}(0,0,1)$. The barycentric coordinates are given as

[TABLE]

The directed tangent on a generic edge $\mathbf{e}_{j}=[j_{1},j_{2}]$ is defined as

[TABLE]

The edge is parameterized as

[TABLE]

A generic edge can be uniquely identified with

[TABLE]

where $\mathrm{sign}(0)=0$. Each face on the 3-simplex can be identified by the associated three vertexes, and is uniquely defined as

[TABLE]

The standard bases in $\mathbb{R}^{n}$ are noted as $\vec{e}_{i}$, $i=1,\cdots,n$, and $n=\{2,3\}$.

2.1 Face functions

The face functions are further grouped into two categories: edge-based face functions and face bubble functions.

Edge-based face functions:

These functions are associated with the three edges of a certain face $\mathbf{f}_{j_{1}}$, and by construction all have non-zero normal components only on the associated face $\mathbf{f}_{j_{1}}$, i.e.,

[TABLE]

where $\mathbf{n}^{\mathbf{f}_{j_{k}}}$ is the unit outward normal vector to face $\mathbf{f}_{j_{k}}$.

Using the idea of recursion from [1], independent edge-based face functions are proposed in [2] as follows.

For $p=1$, for each edge we have one face function for this edge

[TABLE]

and for $p=2$, one additional new basis function can be constructed as

[TABLE]

which can be shown to satisfy the condition (2.6), and and for $p\geq 3,$ the basis functions are given by

[TABLE]

for $\;i=1,\cdots,p-2.$ It can be shown numerically that there are exactly $p$ functions that are independent and whose normal component is non-zero only on the associated edge $\mathbf{e}_{k}$.

Face bubble functions:

The face bubble functions which belong to each specific group are associated with a particular face $\mathbf{f}_{j_{1}}$. They vanish on all edges of the reference 3-simplex $K^{3}$, and the normal components of which vanish on other three faces, i.e.,

[TABLE]

The explicit formula is given as

[TABLE]

[TABLE]

and

[TABLE]

By construction the face bubble functions share an orthonormal property on the reference 3-simplex $K^{3}$:

[TABLE]

for $0\leq\{m_{1},m_{2},n_{1},n_{2}\},m_{1}+n_{1},m_{2}+n_{2}\leq p-3.$

2.2 Interior functions

Interior functions will have zero components on all four faces while some of them may still have non-zero tangential components (Edge-based and face-based interior functions defined below).

The interior functions are classified into three categories: edge-based, face-based and bubble interior functions. By construction the normal component of each interior function vanishes on all faces of the reference 3-simplex $K^{3}$, i.e.,

[TABLE]

Edge-based interior functions:

The tangential component of each edge-based function does not vanish on the associated only edge $\mathbf{e}_{k}:=[k_{1},k_{2}]$ but vanishes all other five edges, i.e.,

[TABLE]

where $\tau^{\mathbf{e}_{j}}$ is the directed tangent along the edge $\mathbf{e}_{j}:=[j_{1},j_{2}]$. The shape functions are given as

[TABLE]

where $i=0,1,\cdots,p-2.$

Again one can prove the orthonormal property of edge-based interior functions:

[TABLE]

Face-based interior functions:

These functions which are associated with a particular face $\mathbf{f}_{j_{1}}$ have non-zero tangential components on their associated face only, and have no contribution to the tangential components on all other three faces, i.e.,

[TABLE]

Further each face-based interior function vanishes on all the edges of the 3-simplex $K^{3}$, i.e.,

[TABLE]

The formulas of these functions are given as

[TABLE]

[TABLE]

where $0\leq\{m,n\},m+n\leq p-3$. The face-based interior functions enjoy the orthonormal property on the reference 3-simplex $K^{3}$:

[TABLE]

for $i=\{1,2\},0\leq\{m_{1},m_{2},n_{1},n_{2}\},m_{1}+n_{1},m_{2}+n_{2}\leq p-3.$

Interior bubble functions:

The interior bubble functions vanish on the entire boundary $\partial K^{3}$ of the reference 3-simplex $K^{3}$. The formulas of these functions are given as

[TABLE]

[TABLE]

where

[TABLE]

Again, one can show the orthonormal property of the interior bubble functions

[TABLE]

where

[TABLE]

In Table 1 we summarize the decomposition of the space $\left(\mathbb{P}_{p}(K)\right)^{3}$ for the $\mathbf{\mathcal{H}}(\mathbf{div})$-conforming tetrahedral $BDM_{p}$ element.

3 Enforcing Divergence-free condition for $B$-field

The magnetic field B in the MHD equations is assumed to be divergence free, however, the time evolution from a fully discretized finite element method for the magnetic equation will render the divergence of $B$ to be non-zero at later time. There are many ways to remove the non-zero divergence in the magnetic field such as the projection method through Helmholtz decomposition. In this paper, we will use the interior functions in the H(div) basis set to correct the non-zero divergence element by element. Due to the vanishing property of the normal components of the interior basis functions, such a local correction will still keep the corrected finite element solution in H(div) globally. The ability of using only the interior functions to reduce the non-divergence error in the magnetic field is based on the following result.

First let us denote the subspace spanned all the interior functions defined in (2.16), (2.20), and (2.22) as

[TABLE]

Lemma 3.1**.**

$\operatorname{div}\Sigma_{\rm int}=P_{p-1}(K)\backslash\{1\}=\{1\}^{\perp}.$**

Proof 3.2**.**

We will prove the result by subspace inclusion argument. First, we will show $\operatorname{div}\Sigma_{\rm int}\subset\{1\}^{\perp}$. Take any function $\Phi\in\Sigma_{\rm int}$, we have

[TABLE]

due to the fact the the normal component of interior function vanishes on the edge of the element K, thus $\operatorname{div}\Sigma_{\rm int}\subset\{1\}^{\perp}.$

On the other hand, we will show $\{1\}^{\perp}\subset\operatorname{div}\Sigma_{\rm int}$ by showing that $(\operatorname{div}\Sigma_{\rm int})^{\perp}\subset\{1\}$, instead. Let $v\in P_{p-1}(K)$ and $v\in(\operatorname{div}\Sigma_{\rm int})^{\perp}$, then we have

[TABLE]

which gives with an integration by parts

[TABLE]

where the vector field $\nabla v\in\left(P_{p-2}(K)\right)^{3}.$

*Now take the tangential vectors of the three edges sharing the common vertex $\mathbf{v}_{0}$,$\tau^{[0,1]}$,

$\tau^{[0,2]}$,$\tau^{[0,3]},$ we can see easily that the following vector functions are also interior functions (with zero normal components on all faces),*

[TABLE]

where the scalar functions $\ g_{i}\ ,i=1,2,3$ are polynomials of degree $(p-2).$ Next, we construct three bi-orthogonal vectors $s_{j}$ with respect to $\tau^{[0,i]},i=1,2,3$ with the following property

[TABLE]

We can express the vector field $\nabla v$ using the basis vector $s_{j}$ as follows

[TABLE]

which will be substituted into (3.4), resulting in

[TABLE]

By setting $\Phi=\lambda_{0}\lambda_{1}f_{1}\tau^{[0,1]}$ in the above identity and using the bi-orthogonality property of (3.6), we have

[TABLE]

which implies that $f_{1}=0$ as both $\lambda_{0},\lambda_{1}$ are positive inside $\overline{K}$. Similar arguement will show that $f_{2}=f_{3}=0$ also hold.

Therefore, we have $\nabla v=0,$ namely, $v=$const., thus $(\operatorname{div}\Sigma_{\rm int})^{\perp}\subset\{1\}.$

Algorithm: we propose a two-step algorithm to remove the non-divergence in the numerical solution for the magnetic field $B$.

• •

Step 1 (local correction) Element-wise removal of high order terms in $\operatorname{div}B$.

Due to Lemma 3.1, we can use the interior function in $\Sigma_{\rm int}$ to remove higher order terms in $\operatorname{div}B$. The remaining component in $\operatorname{div}B$ will be a constant on each element. We proceed to finding a vector function

[TABLE]

such that

[TABLE]

[TABLE]

namely, $div(B+\Phi)\in(\operatorname{div}\Sigma_{\rm int})^{\perp}=\{P_{p-1}(K)\setminus P_{0}\}^{\perp}$, which gives the following linear system for the unique $\Phi\in\Sigma_{\rm int}$

[TABLE]

where $v\in P_{p-1}(K)\setminus P_{0}$. Here (3.13) is a local mixed finite element approximation to the following Poisson equation,

[TABLE]

• •

Step 2 (global correction) remove the constant term in $\operatorname{div}B_{1}$ in the whole domain.

Due to the result of lemma 1, we will have a residual constant term left in the corrected magnetic field $B_{1}$, which can only be removed by a global correction with the first order H(div) basis defined in Section 2. We proceed as follows by finding a second function $\Phi\in H(div,\Omega)$ using the first order $H(div)$ basis functions defined in (2.7),

[TABLE]

such that

[TABLE]

[TABLE]

resulting into the following linear system for $\Phi$, in a similar argument for equations (3.13),

[TABLE]

4 Numerical results

We will solve a magnetic induction field equation on the unit cube $\Omega=[0,1]^{3}$,

[TABLE]

with a periodic $B\cdot\mathbf{n}$ boundary condition is considered

[TABLE]

In (4.1),

[TABLE]

Also for (4.1), the initial condition $B(0,\mathbf{x})$ is given by the exact solution

[TABLE]

We note that due to the initial condition $\operatorname{div}B(0,\mathbf{x})=0$, (4.1) ensures

[TABLE]

We discritize equation (4.1) by a $P_{3}$ $H(\operatorname{div})$ mixed finite element in Section 2, i.e. the space $V_{3}$, on the uniform tetrahedral grids depicted in Fig. 4.1. For time evolution, we use the characteristic method. We compute the solution $B(1/2)$ with 100 time steps with $dt=0.005$. The $L^{2}$ errors are listed in Table 4.1, along with the order of convergence, which is $3$ as expected for the $P_{3}$ mixed finite element. We can see that the numerical solutions are not divergence-free any more as the divergence-free condition for $B_{h}(t,\mathbf{x})$ is not enforced at the time discretization level. Following the correction method proposed in Section 3, we correct the solution $B_{h}(1/2,\mathbf{x})$ by direct solution of both (3.13) and (3.17). The error and the convergence order are listed also in Table 4.1 at the bottom. We can see that the divergence of corrected solution is almost zero, up to the computer accuracy. Moreover, the correction does not change the order of convergence of the finite element solution, though the L2 norm of the error for the corrected solution is slightly larger.

To see the effect of the two correction steps on the solution, we plot in Figure 4.2 the solution $B_{1}(1/2)$ at $t=1/2$ and the errors before and after the correction on a plane at $z=0.485.$ We also compare the new method with the existing, global correction method with the full $P_{3}-H(div)$ basis. To make a comparison to the existing full divergence-zero correction method, we list in Table 4.2 (top part) the errors of the divergence-corrected solution with $P_{3}-H(div)$ basis in (3.17). The correction equations are solved by an iterative Uzawa method, stopped when the $\operatorname{div}$ norm is less than $0.000005$. As expected, we see that the computation complexity is much higher than our new, two-step correction method. We notice that from Table 4.2 that the corrected solution now has a smaller error while achieving an almost zero divergence. In doing the two-step correction (3.13) and (3.17), the work for the global correction (3.17) is of several orders higher than the local correction step. So, we might wish do just the first step correction (3.13). However, the corrected solution is not divergence-free anymore and we list in the bottom half of Table 4.2 the errors and the order of convergence for such one-step corrected solutions.

5 Conclusion

In this paper, we have proposed an efficient correction procedure to ensure the divergence free condition of the magnetic field. The correction is done in two steps: the first step can be done locally on each element which removes the high order terms in the divergence error of the B-field; the second step is a global one which removes the remaining constant term in the divergence error on each element. Numerical results have shown the effectiveness of the proposed method in enforcing divergence-free condition for a magnetic induction equation while maintaining the accuracy of the solution itself.

Acknowledgement

The author (W.C) acknowledges the support of the US Army Office of Research (Grant No. W911NF-14-1-0297) and US National Science Foundation (Grant No. DMS-1315128) and the National Natural Science Foundation of China (No. 91330110) for the work in this paper. The research of the author (J.H.) is supported by NSFC projects 11271035, 91430213 and 11421101. Authors also like to thank Prof. Lin-bo Zhang for helpful discussions on implementation of H(div) basis.