A Quadratic C Interior Penalty Method for the Quad-Curl Problem

In this paper we study the C interior penalty method for a quad-curl problem arising from magnetohydrodynamics model on bounded polygons or polyhedrons. We prove the well-posedness of the numerical scheme and then derive the optimal error estimates in a discrete energy norm. A post-processing procedure that can produce C approximations is also presented. The performance of the method is illustrated by numerical experiments.


Introduction
The magnetohydrodynamics (MHD) model [1,11] has wide range of applications in plasma physics, astrophysics, magnetospheric and thermonuclear fusion. It describes the motion of electrical conducting fluid in the presence of a

Copyright c 2020 The Author(s). Published by VGTU Press
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. magnetic field. In the MHD model, the behavior of a continuous fluid is governed by a simplified form of Maxwell's equation, together with Ohm's law and Navier-Stokes equations. The MHD system with hyper-resistivity is described as follows: where u is the fluid velocity, B is the magnetic field, ρ is the density, p is the static pressure, η 1 is the resistivity, η 2 is the hyper-resistivity, µ 0 is the magnetic permeability of free space, µ is the viscosity, and d i is the ion inertial length.
When simulating the MHD system, the velocity u and pressure p are often discretized by using the standard finite element method. However, for the magnetic field variable B, such approach will encounter difficulties because of the existence of the fourth order term (∇×) 4 B. In particular, an H 2 -conforming method would have 220 degrees of freedom per element to approximate the solutions. Hence it is more practical to use the nonconforming finite element method [26] or mixed finite element method [20] to solve the above quad-curl model problem. However, other problems may arise. More precisely, though the number of degree of freedom can be greatly reduced by using nonconforming and mixed finite element methods, the implementations of Morley-type elements [14,16] and high order Nédélec's edge elements [17,18] are more challenging than the simple standard P k Lagrange finite elements. Moreover, the error analysis for nonconforming and mixed methods is more complicated. The MHD model problem was also solved by the discontinuous Galerkin (DG) method [13] using weakly H(curl)-conforming elements. Recently, Brenner et al. [5] studied the P k Lagrange finite element methods for the quad-curl problem based on the Hodge decomposition for divergence-free vector fields. Their idea is to decompose the two-dimensional model vector problem into a sequence of second order elliptic scalar problems, which makes the two-dimensional problem easier to be solved by standard conforming FEM. Recently, Zhang et al. [24] developed an H 2 (curl)-conforming finite element in two dimensions and applied it to solve the quad-curl model problem. Chen et al. [12] proposed a mixed finite element method for a problem with a quad-curl term. In this work, they also used the standard P k Lagrange finite element to approximate the variables. However, the stability requires the inf-sup condition, which is a complicated process. This analysis was obtained by using the discrete Sobolev embedding inequalities for the piecewise polynomials.
In this work, we focus on the following three dimensional quad-curl model problem which is deduced from the magnetic induction equations in the MHD model [26]. Let Ω be a bounded polyhedral domain in R 3 and f ∈ [L 2 (Ω)] 3 . We consider the model problem in Ω, where (∇×) 4 = (∇ × ∇ × ∇ × ∇×), β and γ are nonnegative constants, and γ > 0 if Ω is not simply connected. The quad-curl model problem (1.1) also appears in the Maxwell's transmission problem [15] and the inverse electromagnetic scattering theory [9,10]. The weak problem of (1.1) is as follows: find u ∈ E such that in Ω, and n × v=0 on ∂Ω .
In this paper, we consider the C 0 interior penalty (C 0 -IP) method [2,6] to solve the quad-curl problem (1.2). One major merit of the C 0 -IP method is that it uses lower order P k polynomials, which is much simpler than C 1 -conforming and curl-curl conforming finite element methods. Therefore, both the construction of the finite element space and the error analysis of the C 0 -IP method are straightforward. Moreover, the positive definiteness of the continuous problem can also be preserved by the numerical scheme even for more complicated elliptic systems. On the other hand, for the C 0 -IP method, there is no need to include penalization terms involving the jumps of discrete functions. Hence the computational cost is considerably less than for the conventional DG methods. These advantages make C 0 -IP method more attractive for solving fourth order problems [8]. Furthermore, the C 0 -IP scheme for lower order elliptic problems can be exploited in the design of an effective smoother for fast solvers, such as multigrid algorithms for higher order problems [7,21].
The rest of this paper is organized as follows. In Section 2, we introduce the C 0 interior penalty method for (1.2) and prove the well-posedness of the numerical scheme. In Section 3, we derive optimal error estimates in a discrete energy norm. A post-processing procedure which further produces C 1 approximations is also studied. Several numerical experiments are carried out to illustrate the performance of the numerical scheme. The results are reported in Section 4.
For convenience, we will use the notation A B throughout the paper to represent the inequality A ≤ Constant×B, where the constant, unless otherwise specified, depends only on the shape regularity (or equivalently the minimum angle) of T h . The statement A ≈ B is equivalent to A B and B A.
2 The C 0 interior penalty method

Preliminaries
Let T h be a regular tetrahedron mesh of Ω ∈ R 3 . The diameter of a tetrahedron The construction of the quadratic C 0 interior penalty method requires the concepts of jumps and average of functions across the element interfaces F h , which are defined below.
Let F ∈ F i h be an interior face shared by two tetrahedrons T − and T + , and n e be the unit normal of F pointing from T − to T + . We define on F we take n e to be the unit normal of F that points towards the outside of Ω and define Remark 1. Note that in R 2 , ∇ × v is a scalar function. Hence the average terms are defined as the same formulation as in three dimension. However, the jump terms are defined to be vector functions as below.
where t e is the unit tangent vector which is obtained by rotating n e by a right-angle clockwisely.

The discrete problem
be the P 2 -vector Lagrange finite element space associated with T h , such that the tangential components of its members vanish on all nodal points on the boundary of Ω. More precisely, n × v T = 0 on the vertices and midpoints of each edge of F ∈ F b h }.
Let u be the solution of the quad-curl problem (1.1) and v ∈ [P 2 (T )] 3 . We assume that u ∈ [H 2+α (Ω)] 3 so that the construction of the discrete scheme makes sense. For any T ∈ T h , by using Green's formula we have By summing up (2.1) over all T ∈ T h and considering the definition V h , we find where |e| is the diameter of the face F , and σ > 0 is the penalty parameter to be chosen later. The C 0 interior penalty method for the quad-curl problem (1.2) is defined as follows: find u h ∈ V h such that Remark 2. The first four terms on the right-hand side of (2.4) come from integration by parts. The fifth term is added for symmetry. The sixth term is the penalization term (with penalty parameter σ > 0). In the analysis of wellposedness, the penalty parameter σ is chosen to be large enough in order to derive the coercivity of the discrete bilinear form a h (·, ·). The last term is added to the numerical scheme to (weakly) enforce the divergence-free condition for the discrete solution. The inclusion of this term provides a good approximation of divergence free condition. We can also remove this term from the discrete scheme by adding a locally divergence free condition in the discrete space V h . However, such approach will make the construction of V h more complicated.

Well-posedness of the discrete problem
It follows from (2.3) and Lemma 1 that the following Galerkin orthogonality holds Define the discrete energy norm as

Lemma 2 [Continuity]
. It is easy to show that the bilinear form a h (·, ·) is bounded in the sense that There exists a positive constant C 1 depending only on the shape regularity of T h such that where the seminorm ||| · ||| H 2 (Ω,T h ) is defined by Proof. The first inequality in (2.8) is obvious by (2.6) and (2.9). By the trace theorem(with scaling) and standard inverse estimates, we have Moreover, it follows from Poincare-Friedrichs inequality (cf. [4]) for piecewise H 1 functions that The second inequality in (2.8) now follows from (2.6), (2.9), (2.10) and (2.11).

Lemma 4 [Coercivity]
. When the penalty parameter σ is chosen to be sufficiently large, there exists a constant C * such that Proof. By Cauchy-Schwarz inequality and (2.10), we have for all v ∈ V h , where > 0 is arbitrary and C 2 > 0 depends only on the shape regularity of T h . By choosing = 1 2C2 , σ = 1 2 + 1 , and in view of (2.4), (2.8) and (2.9), we get Combining Lemma 1, Lemma 2 and Lemma 4, the discrete problem (2.3) is well-posed under the assumptions of Lemma 4.

Enriching operator
In this subsection, we introduce an enriching operator which will be used later to construct a post-processing procedure that produces a C 1 approximation solution. The enriching operator measures the distance between V h and the Sobolev space 3 . The construction of E h uses the C 1 -Ženíšek finite element space (cf. [22]).
We briefly introduce the idea as follows. More details can be found in [6]. Let P i , P j , P k , P l be the four vertices of the tetrahedron T , Q (1,s) jk , ..., Q (s,s) jk be the points dividing the edge P j P k into s + 1 equal parts, Q i be the triangular face which lies opposite to the vertex P i , P 0 be the center of gravity of the tetrahedron T , as shown in Figure 1. For any v ∈ V h , we define w = E h v by averaging. More precisely, at the interior vertices P of T h , we take (2.14) where T p is the set of tetrahedrons in T h whose closures share the nodal point P .
At the points on edges, we take jk ), |β| = s, r = 1, ..., s, . s = 1, 2, m = 1, 2, 3. At the center of gravity Q of the triangular face, we take (2.16) At the center of gravity P 0 of the tetrahedron, we take Proof. Let T ∈ T h be arbitrary and N(T ) be the set of nodal variables of the C 1 -Ženíšek finite element. It follows from (2.13), (2.16), (2.17) and scaling that Here V T is the set of the four vertices of T and M T is the set of the centers of gravity of the triangular faces. From (2.15) we have It then follows from scaling that Similarly, by using (2.14) we have It follows from (2.22) that where F p is the set of the faces in F h that share the common endpoint P . Similarly we have

Corollary 1. We have
and v have identical dofs, it directly follows from the construction of the enriching operator that the following property holds The following result shows that the composition E h • Π h behaves like a quasi-interpolation operator (cf. [6]). Lemma 6. For any s ≥ 2, there exists a positive constant C depending only on s and the shape regularity of T h such that where S T is the union of the tetrahedrons in T h that share a common vertex with T .

Error estimates in energy norm
In this section, we derive discretization error estimates for the C 0 interior penalty method introduced in section 2. The next lemma, which gives the abstract error estimate, is a direct consequence of the (2.5), (2.7) and (2.12). Lemma 7. Let u and u h be the solutions of (1.2) and (2.3), respectively. Then the following estimate holds: By the definition of Π h , we have the standard interpolation error estimate (cf. [6]): Lemma 8. Let u ∈ E be the solution of problem (1.2), and assume u ∈ [H 2+α (Ω)] 3 , 1 2 < α ≤ 1. Then the following interpolation error estimate hold: Proof. According to (2.6), we have It follows from (3.1) that By using trace theory with scaling and standard interpolation error estimate, we have It follows from (2.10) and (3.4) that Finally, the interpolation error estimate (3.2) holds by combining (3.3)-(3.5).
The error estimate for the C 0 -IP method is given in the following theorem, which directly follows from Lemma 7 and Lemma 8.

Error estimate in a lower order norm
Now we investigate the discretization error in lower order norms for C 0 interior penalty method.
Also, we have the duality formula The estimate (3.7) follows from (3.9) and (3.10). 3 be the enriching operator introduced in subsection 2.4. The C 1 finite element function E h u h obtained by post-processing provides an (smooth) approximation of u in the space W h . Recall that S T is the union of polyhedrons in T h that share a common vertex with T , and note that u ∈ H 2+α(S T ) (S T ), where α(S T ) = min T ∈S T α(T ).

Post-processing
The error estimate for E h u h is given in the next theorem.
where α is the index of elliptic regularity and the positive constant C depends only on the shape regularity of T h .

Numerical results
In this section we report the results of numerical experiments in two dimensions. All three numerical experiments are carried out on quasi-uniform triangulations. We take α and β both to be 1. Besides the errors in the energy norm and L ∞ norm, we also include the errors in the norm · H(curl,Ω) which is defined by u H(curl,Ω) = ∇ × u L 2 (Ω) . Note that the choice of φ ensures that u ∈ E. The domain for this experiment is [0, 1] 2 . The relative errors are reported in Table 1. The results show that the In [23], we also studied a mixed finite element method for this example. In order to provide a comparison in terms of efficiency, we plot the convergent error of u − u h H(curl,Ω) / u H(curl,Ω) with computational time for the C 0 -IP method and the mixed method in Figure 2. We observe that the C 0 -IP method takes less time than the mixed finite element method in order to achieve the same accuracy. We notice that this example was also studied by Brenner et al. in [5] by using Hodge decomposition method. We give a comparison of accuracy for u − u h L 2 (Ω) with degree of freedoms (DOFs) in Table 2. We see that C 0 -IP method method takes less DOFs to achieve a higher accuracy for the P 1 case. The accuracy for Hodge decomposition using P 2 Lagrange finite element method is better than C 0 -IP method at the cost of more DOFs.
to be the right-hand side function of (1.2). We take Ω = [0, 1] 2 in this experiment. Since the exact solution is not available in this case, we use the differences of numerical solutions on consecutive levels to measure the errors. The relative errors are reported in Table 3. The errors converge with the first order rate in both the energy norm and the · H(curl,Ω) norm, which agree with the theoretical results.  1, 3)]. The exact solution is taken to be u = ∇ × φ, where φ(x) = sin 3 (4πx 1 ) sin 3 (4πx 2 ). In this case, Then we have and ∇ × u = 0 on ∂Ω. It is obviously that u ∈ E since ∇ · u = 0, and n × u = 0 on ∂Ω. We depict the second components of the (enriched) numerical solution E h u h on consecutive levels in Figure 3 and that of the exact solution on final level in Figure 5, i.e., the nodal values of exact solution on the mesh with h = 1/64 are depicted. We also present the second components of numerical solution u h obtained by C 0 -IP method in Figure 4. It can be observed that the post-processed numerical solutions are C 1 smooth and quickly converge to the true solution.  In [21], we studied multigrid fast solvers based on the C 0 -IP method. In that paper, we test W -cycle multigrid algorithm for the same numerical example. We present the errors and computational time on each level in Table 4. Here u M G h is the numerical solution obtained by multigrid method. We see that   the computation time is greatly reduced by using multigrid method on higher levels, while the errors remain the same order of accuracy.

Conclusions
In this paper, we developed the C 0 interior penalty method for solving the quad-curl problem which arises in the MHD model. We established the wellposedness of the numerical scheme, in which we weakly enforce the divergencefree condition by adding an extra penalty term. Optimal convergence rates are proved on convex polyhedral/polygonal domains based on uniform meshes. We further constructed an enriching operator for the numerical solution in three dimensions. It provides an approximation of u in the C 1 -Ženíšek finite element space and also preserves the convergence property of C 0 interior penalty method.
We would be able to further improve the convergence rate of the new C 0 -IP method on nonconvex domains by using graded meshes [3]. The grading strategy will be closely related to the elliptic regularity of the fourth-order model problem. Moreover, we would like to study the a posteriori error estimate for the C 0 -IP method and its fast solvers. Those will all be carried out in our future work.