{'fraction_non_alphanumeric': 0.04842431808723944, 'fraction_numerical': 0.0169485113305338, 'mean_word_length': 4.00738776283387, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Adaptive Flux-Only Least-Squares Finite Element Methods for Linear Transport Equations

In this paper, two flux-only least-squares finite element methods (LSFEM) for the linear hyperbolic transport problem are developed. The transport equation often has discontinuous solutions and discontinuous inflow boundary conditions, but the normal component of the flux across the mesh interfaces is continuous. In traditional LSFEMs, the continuous finite element space is used to approximate the solution. This will cause unnecessary error around the discontinuity and serious overshooting. In arXiv:1807.01524 [math.NA], we reformulate the equation by introducing a new flux variable to separate the continuity requirements of the flux and the solution. Realizing that the Raviart-Thomas mixed element space has enough degrees of freedom to approximate both the flux and its divergence, we eliminate the solution from the system and get two flux-only formulations, and develop corresponding LSFEMs. The solution then is recovered by simple post-processing methods using its relation with the flux. These two versions of flux-only LSFEMs use less DOFs than the method we developed in arXiv:1807.01524 [math.NA]. Similar to the LSFEM developed in arXiv:1807.01524 [math.NA], both flux-only LSFEMs can handle discontinuous solutions better than the traditional continuous polynomial approximations. We show the existence, uniqueness, a priori and a posteriori error estimates of the proposed methods. With adaptive mesh refinements driven by the least-squares a posteriori error estimators, the solution can be accurately approximated even when the mesh is not aligned with discontinuity. The overshooting phenomenon is very mild if a piecewise constant reconstruction of the solution is used. Extensive numerical tests are done to show the effectiveness of the methods developed in the paper.


Introduction
In this paper, we consider the following linear transport equation in the conservative form: with β an advection field and Γ − the inflow boundary. Detailed descriptions of the equation can be found in Section 2.
It is crucial to realize that unlike the elliptic or parabolic equations whose solutions are generally smooth, hyperbolic equations commonly have discontinuous solutions. For the transport equation, the solution space is W = {v ∈ L 2 (Ω) : ∇ · (βv) ∈ L 2 (Ω)}.
For the space W, there is no simple finite element subspace of it. The usual continuous finite element space requires too much continuity while the discontinuous finite element space lacks the continuity of the flux βu. To separate the continuity requirement, we realize that for a true solution u ∈ W, the condition ∇ · (βu) ∈ L 2 (Ω) essentially means that u ∈ L 2 (Ω) and βu ∈ H(div; Ω).
In [13], we propose new least-squares variational formulations with a flux reformulation. Introducing the flux σ = βu ∈ H(div; Ω), we have a first order system (with appropriate boundary conditions): σ − βu = 0 and ∇ · σ + γu = f. (1. 2) The space requirement of u is reduced to L 2 (Ω) only. In [13], we use H(div; Ω)-conforming Raviart-Thomas finite element space RT k and discontinuous piecewise polynomial space P k to approximate σ ∈ H(div; Ω) and u ∈ L 2 (Ω) separately. Notice that the space RT k has enough degrees of freedom to approximate both σ and the divergence of σ in the same order k, this suggests us that it is possible to remove the u-variable from the system and get a system with the flux only to reduce the degrees of freedom of the discrete system. In this paper, two reformulations are suggested. In the first reformulation, under the assumption that |β| 0, we use σ = βu to get u = σ · β/|β| 2 , then a single equation for σ is obtained: Since continuous functions are used to approximate discontinuous solutions, it will introduce unnecessary error. The error indicator will always indicate unnecessary big errors for those elements on the discontinuity if an adaptive method is used. Compared with a P 0 discontinuous approximations for u which is easily available for the methods in this paper, Gibbs phenomenon like spurious over-shooting is more intense near the discontinuity in these continuous approximations, see also discussions in [17] on overshooting with adaptive finite elements.
The nonconforming LSFEM in [7] and the similar method in [14] use discontinuous approximations with a penalty term to enforce the continuity of the normal component of the flux. The H(div) flux-only methods use natural finite element spaces without any penalty terms. In addition, these two methods have to use at least piecewise linear functions, which can lead to non-trivial overshooting, see discussions in [17]. The methods developed in this paper can simply recover the solution in a piecewise constant space to avoid overshooting.
The paper is organized as follows. Section 2 describes the model linear hyperbolic transport problem. Based on two flux-only reformulations, two least-squares variational problem are presented in Section 3. Corresponding LSFEMs are developed in Section 4, a priori and a posteriori error estimates are established. Section 5 provides numerical results for many test problems. In Section 6, we make some concluding remarks.

Model Linear Hyperbolic Transport Equation
Let Ω be a bounded polyhedral domain in d (d = 2 or 3) with a Lipschitz boundary. Assume that (2.1) Define the inflow and outflow parts of ∂Ω in the usual fashion: where n ∂Ω (x) denotes the unit outward normal vector to ∂Ω at x ∈ ∂Ω. In order to specify the boundary condition in (1.1), define In this paper, we assume that the coefficients are nice enough to guarantee the existence and uniqueness of the solution.
Assumption 2.1. (Assumption on the data β and γ) Assume that for g ∈ L 2 (|β · n|; Γ − ), the linear transport equation (1.1) has a unique solution in W for the given data β and γ.
Several known conditions that guarantee the existence and uniqueness of the solution are: (i) There exists a positive γ 0 , such that The proof of the uniqueness and existence with the condition (i) can be founded in [7] and Chapter 2 of [9]. For the case with the condition (ii) (and more general cases), the proof is based on the standard ODE theory, and can be founded in [7,8].
2. An equivalent non-conservative reformulation is We can apply the methods developed in this paper to the non-conservative version by changing it to the conservative formulation.

Flux-Only Least-Squares Variational Problems
In this section, two flux only least-squares variational problems are introduced. The existence and uniqueness of the formulations are discussed.

Reformulation with both flux σ and solution u
In [13], a flux σ = βu is introduced to separate the continuity requirements, then the divergence of the flux is in L 2 (Ω), ∇ · σ = f − γu ∈ L 2 (Ω). Thus the flux σ ∈ H(div; Ω). Since the inflow condition is essentially a condition of the flux, we write the condition u = g on Γ − as Thus we get the following system of σ and u: in Ω, Define the following spaces: H 0,− (div; Ω) := {τ ∈ H(div; Ω) : τ · n = 0 on Γ − }.

Flux-only least-squares variational problem: first reformulation
In this reformulation, we assume: From the equation σ = βu, we have Then we obtain a single equation of σ: Such σ is not unique, since a single scalar equation is used to compute a vector function σ. For example, when γ = 0, the only requirement is ∇ · σ = f . Thus some extra conditions are needed to ensure the uniqueness of σ. Notice that since σ = βu, we should also have is an orthonormal family of vectors satisfying For d = 2, we can simply choose β ⊥ (x) = (−β 2 (x), β 1 (x))/|β(x)|. For d = 3, β (1) ⊥ (x) and β (2) ⊥ (x) can be obtained by first finding two linearly independent vector functions orthogonal to β(x), then performing a Gramm-Schmidt orthonormalization procedure.

Remark 3.2. It is also clear that
is a semi-norm on an element K ∈ T .

Flux-only least-squares variational problem: second reformulation
If we assume that γ 0, or more precisely, we assume that by σ = βu. Notice this is a vector equation, thus this reformulation does not need the extra orthogonal condition to guarantee the uniqueness of σ. We can also reformulate it as We get the following linear system with respect to σ for γ 0: The least-squares variational problem of (3.11) is: we seek σ ∈ H g,− (div; Ω), such that with the least-squares functional J 2 defined as Its corresponding Euler-Lagrange formulation is: find σ ∈ H g,− (div; Ω), such that where the bilinear form a 2 is defined as for all τ, ρ ∈ H(div; Ω), and Remark 3.4. Note that f 1 (τ) and f 2 (τ) are essentially the same if we weight the system with |β| 2 (or simply normalize β to make |β| = 1), but the a 1 and a 2 are not the same when γ 0 since (β · τ, β · ρ) (|β| 2 τ, ρ).
Thus these two reformulations are indeed different.
Lemma 3.5. Assuming that the coefficients β and γ 0 satisfy the assumption of the data ensuring the existence and uniqueness of original equation (1.1), the following defines a norm for τ ∈ H 0,− (div; Ω): Proof. We only need to check that This yields ∇ · (βv) + γv = 0. Similar to Lemma 3.1, we have v = 0 on Γ − , and thus v = 0 in Ω and τ = βv = 0. The lemma is proved.
Remark 3.7. In traditional least-squares methods, a norm equivalence is often seeked to prove to the existence and uniqueness. As we will see, it is not possible in our case. Thus an indirect proof is used. Another example of such least-squares method is its application in non-divergence equation, see [16].

Flux-only Least-Squares Finite Element Methods
In this section, least-squares finite element methods based on the flux-only least-squares variational problems are developed. A priori and a posteriori error estimates are derived.
Let T = {K} be a triangulation of Ω using simplicial elements. The mesh T is assumed to be regular. We denote the set of edges/faces of the triangulation T on inflow boundary Γ − by E − . For an element K ∈ T and an integer k ≥ 0, let P k (K) be the space of polynomials with degrees less than or equal to k. The space P k (F) is defined similarly on an edge/face F. Define the finite element spaces RT k and P k as follows: Assumption 4.1. (Assumption on the boundary data) For simplicity, we assume (β · n)g on Γ − can be approximated exactly by the trace of RT k space on Γ − , i.e., g| F ∈ P k (F), for all faces/edges F ∈ E − .
Note that this assumption still allows a discontinuous boundary condition, but it does require that the boundary mesh is aligned with the discontinuity. For an arbitrary g, we need to first interpolate or project (β · n)g to the piecewise polynomial space.

Interpolations and their properties
In order to derive a priori error estimates, we introduce some interpolations and their properties. Note that all properties here are local with respect to the element K.
For s > 0, denote by I rt k : H(div; Ω)∩[H s (Ω)] d → RT k the standard RT interpolation operator [5]. It satisfies the following approximation property: In addition, another approximation property holds: for τ ∈ H s (K) and ∇ · τ ∈ H s (K), s > 0, The discussion of the above two properties can be found in [13].

Least-squares finite element methods
Or equivalently, find σ h ∈ RT k,g,− , such that We call LSFEM1 and LSFEM2 for least-squares finite element methods with i = 1 and 2 separately. Note that we need the assumption |β| > 0 for LSFEM1, and the assumption γ 0 is needed for LSFEM2.

Post-processing for u h
Once the numerical flux σ h is known, since we have two equations to relate u and σ, there are two simple ways to recover a discrete approximation u h : The second recovery Note that both post-processing procedures can be used for both LSFEMs as long as the assumption of the coefficients holds. We can also use a simple element-wise L 2 -projection to get an approximation u h ∈ P k if needed. To ensure the overshoot effect is mild, we should project u h onto P 0 as suggested in [17].

A priori error estimates
The following best approximation property is simple and straight. satisfying the corresponding condition on β and γ and the assumption on the boundary data, the following best approximation result holds: Proof. For an arbitrary τ h ∈ RT k,0,− , the following error equation holds: From the definition of the norm | | | · | | | i , the error equation, and the Cauchy-Schwarz inequality, we have Since τ h is chosen arbitrarily, the theorem is proved.
Define the following piecewise function spaces on the triangulation T , with s a piecewisely defined function, s| K = s K > 0.
Remark 4.5. For both | | | · | | | i norms, we have the following simple upper bound w.r.t. the standard H(div)-norm: with C depending on β and γ only.
On the other hand, our numerical test will also disprove the possibility of a coercivity w.r.t. the standard H(div)-norm: or a weak discrete version with an h-independent C > 0, Because if one of such coercivity results holds, we will have For a uniform mesh and an RT 0 approximation of σ h , we should have for piecewise smooth solutions with a discontinuity aligned mesh or a global smooth solution with a Peterson's mesh, which is not the case in our numerical examples 5.2, 5.3, and 5.4.

A posteriori error estimation
The least-squares functionals can be used to define the following fully computable a posteriori local indicators and global error estimators: Theorem 4.6. The a posteriori error estimator η i is exact with respect to the least-squares norm | | | · | | | i : The following local efficiency bound is also true with C > 0 independent of the mesh size h: Proof. For LSFEM1, using the facts f = ∇ · σ +γ(β · σ) and σ · β (i) ⊥ = 0, ∀i ∈ I, for the exact solution σ, we get The proofs of LSFEM2 and the local exactness are identical.
With local bounds (4.9) and (4.10), the local efficiency bound for the standard H(div) norm can be easily proved.
Remark 4.7. Due to the fact that the least-squares functional norm is not equivalent to the standard H(div) norm, it is impossible to get the corresponding reliability result w.r.t. the H(div) norm. 11

Computational Examples
In all our numerical examples, the lowest order element RT 0 is used to approximate the flux σ. In the adaptive mesh refinement algorithm, the Döfler's bulk marking strategy with θ = 0.5 is used and the algorithm is stopped when the total number of DOFs reaches 10 6 . All refinements are based on the longest edge bisection algorithm. For all the numerical examples with domain (0, 1) 2 , the mesh shown in Fig. 1 is used as an initial mesh. Since we have two formulations for σ h and two recoveries for u h , there are four variants to find the numerical solution pair (u h , σ h ). We use LSFEMi-j to denote them, with i denoting the method to get σ h and j denoting the method to recover u h . For example, LSFEM1-2 means we use LSFEM1 to get σ h and use the second recovery to get u h . For almost all our numerical tests, we find that the results of the four combinations are identical. We only show the full 4 combinations for one example, and show only one option for other examples.

5.
1. An example with a constant advection field and a piecewise constant solution, on a matching grid In this example, we only need to do a thought experiment, although the actual computation does confirm our result.
Consider the following problem: Ω = (0, 1) 2 with β = (1/ √ 2, 1/ √ 2) T . The inflow boundary is {x = 0, y ∈ (0, 1)} ∪ {x ∈ (0, 1), y = 0}, i.e., the west and south boundaries of the domain. Let γ = 1 and choose g and f such that the exact solution u is If we choose the mesh aligned with the discontinuity, for example, any refinements of the mesh in Fig. 1. Note that the exact σ ∈ RT 0 . By the best approximation property Theorem 4.2, the numerical solution σ h is identical to the exact solution. The recovered u h is also identical to u. So no further refinements are needed. This is not true when C 0 finite elements are used to approximate the discontinuous u as in [1,8,3,4], many unnecessary refinements are needed.

An example with a global smooth solution
Consider the following simple problem: Ω = (0, 1) 2 with β = (1, 1) T . The inflow boundary is {x = 0, y ∈ (0, 1)} ∪ {x ∈ (0, 1), y = 0}, i.e., the west and south boundaries of the domain. Let γ = 1. Choose f and g such that the exact solution is u = sin(x + y). In Fig. 2, the convergence histories on uniformly refined meshes are shown. For all 4 methods, errors measured in the least-squares energy norms, σ − σ h H(div) , and u − u h 0 are all of order 1.
In Fig. 2, we also notice that all combinations of the methods have almost identical numerical results.
The exact solution is u = x. The mesh is chosen to be in the pattern on the left of Fig. 3. Since γ = 0, we only test this example with the method LSFEM1-1. We compute a series of solutions by LSFEM1-1 on meshes with h from 1/6, 1/12, · · · , to 1/768. The convergence result is plotted on the right of Fig. 3. It is observed that the error in LS norm still converges in the order of 1, but the L 2 -norm of u − u h and H(div)-norm of σ − σ h only converges in the order of 3/4.
We choose an initial mesh that matches the discontinuity (Fig. 1) and uniformly refine it for 8 times. In Fig. 4, we show the convergence histories. The convergence order of the error in leastsquares norms is 1, which matches the optimal convergence theory. The orders of σ − σ h H(div) and u − u h 0 are less than 1 (about 0.6 at late stages). This one and example 5.3 suggest that the norm equivalence (or in discrete sub-spaces): should not be true.

An example with a piecewise smooth solution, on non-matching grid
In this example, we discuss the over/undershooting of the solution when the mesh is not matched with discontinuity.
Consider the problem: Ω = (0, 2) × (0, 1) with β = (0, 1) T . The inflow boundary is {x ∈ (0, 2), y = 0}, i.e., the south boundary of the domain. Let γ = 1 and f = 1. Choose the inflow boundary condition We set the initial mesh to be as shown in Fig. 5. The bottom central node is (π/3, 0) and the top central node is (1, 1). So the inflow boundary mesh is matched with the inflow boundary discontinuity while the mesh is not aligned with the discontinuity in general and will never match it if a bisection mesh refinement is used. On the left of Fig. 6, we show the solution u h on a mesh after 8 uniform refinements of the initial mesh. In order to study the overshooting phenomenon on the outflow boundary, we only draw the graph of u h on y = 1, that is, we plot the u h value at the midpoint of x-axis of each elements with edges on y = 1. The solutions are almost identical for all 4 combinations. Small under/overshooting can be found. The maximum of numerical u h on y = 1 is 1.0430 with the exact solution u = 1, and the minimum of numerical u h is 0.6210 with the exact u = 0.6321.
On the right of Fig. 6, we plot the convergence results of uniform refinements. The order of convergence of errors measured in least-squares norms is still about 1, the optimal order. The order of u − u h 0 is smaller than 1/2. We then test the problem with adaptive algorithm. In Fig. 7, adaptively refined meshes after several iterations are shown. Clearly, more refinements are required along the discontinuity. We do find that for both two LSFEMs, the meshes are almost identical. In Fig. 8, the convergence histories for the adaptive meshes are shown. The convergence order of the error measured in the least-squares norm is optimal with order 1, while the order of u − u h 0 is about 1/2, which is the same order (also the best possible order) as uniform refinements for discontinuous solution on an aligned mesh, Example 5.4.
On the left of Fig. 9, we show the decreasing of the overshooting values by adaptive mesh refinements. Here, the overshooting value is defined as max(max(u h − 1), min(1 − e −1 − u h )) on y = 1, still along the outflow boundary. We clearly see that the overshooting value begins to decrease after the mesh is reasonably fine. When the mesh is very coarse, the overshooting is actually not very severe since u is approximated by P 0 .
On the right of Fig. 9, we show the numerical solution u h on the outflow boundary for the final mesh. It is clear that when with the adaptive mesh and P 0 recovery, the overshooting phenomenon is almost neglectable. 16 Note that with an initial mesh as in Fig. 1, any refinement of it will never match the discontinuity.
We show the convergence results after uniform refinements in Fig. 10. The convergence order in the least-square norm is about 0.8. Similar to example 5.5, it is worse than order 1 but better than order 1/2. The convergence order for u − u h 0 is about 0.3.
On the left of Fig. 11, adaptively refined meshes after several iterations are shown. Many refinements are generated near the discontinuity. On the right of Fig. 11, convergence histories of the adaptive LSFEMs are shown. The rate of convergence of the error in least-squares norms is about order 1, and the order of u − u h 0 is about 0.5.

Curved transport example
We consider an example similar to Example in 4.4.2 of [10]. Consider the problem on the half disk Ω = {(x, y) : x 2 + y 2 < 1; y > 0}. Let the inflow boundary be {−1 < x < 0; y = 0}. Choose the advection field β = (sin θ, − cos θ) T = (y/ x 2 + y 2 , −x/ x 2 + y 2 ) T , with θ the polar angle. Let γ = 0, f = 0, and the inflow condition and the exact solution be We use an initial mesh as shown on the left of Fig. 12. We choose the bottom central node to be (0, 0) and the node left of it to be (−0.5, 0). So the inflow boundary mesh is alighed with the inflow boundary discontinuity. Since the advection field is curved and so is the discontinuity, the mesh will never be aligned with the discontinuity even after refinements. Since the boundary is a half circle, when we refine the mesh, we take an extra step to map those boundary nodes to the right positions on the circle.
Since γ = 0 in this example, so we only use the LSFEM1-1 to compute it. We show the numerical recovered solution u h computed by LSFEM1-1 on a mesh after 8 uniform refinements of the initial mesh on the right of Fig. 12. Small overshootings can be observed near the discontinuity. Along the radius, the solution is essentially one dimensional, we project the graph of the solution onto the radius, see the left of Fig. 13. We do see the under and overshooting. The maximum and minimum values of numerical solution u h are 1.0404 and −0.0383, respectively.
With uniform refinements, the convergence rate of the error in the least-squares norm is about 0.75 and the rates of u−u h 0 and σ−σ h H(div) are about 0.28, see the right of Fig. 13. Since the mesh is not aligned with the discontinuities, the convergence order of the least-sqaures energy norm is smaller than 1. On the left of Fig. 14, we show the adaptive mesh generated by LSFEM1-1 after several iterations. Along the discontinuity, we see many refinements. Also, almost uniform refinements can be found in the half ring where u = 1. The reason is that even u is a constant 1, the flux σ = β is not a constant vector and has approximation errors. On the other hand, in the region where u = 0, the flux is also a zero vector and can be exactly computed. So no refinement is needed in the inner half ring. 19 On the right of Fig. 14, we show the convergence history of the adaptive method. With the adaptive least-squares finite element method, the convergence order of the error in the leastsquares norm is about 1 and is optimal, and the rates of u − u h 0 and σ − σ h H(div) are about 0.5. On the left of Fig. 15, we show the reduction of overshooting values obtained by LSFEM1-1. After the initial stages, the overshooting values are decreasing with refinements along the discontinuities. On the right of Fig. 15, the projected solution on the radius is shown on the final mesh. We can see that the overshooting is very small compared with the uniform refinements. The Gibbs phenomena is not observed.
20 When = 0.01, the layer can be fully resolved. When = 10 −10 , the layer is never fully resolved in our experiments and can be viewed as discontinuous.  Figs. 16 and 18. The behaviors of the methods are very similar to the piecewise smooth solution example with a non-matching grid, the example 5.6. The order of convergence of u − u h 0 is about 0.12. The contour of the solution on the right of Fig. 18 shows that the overshooting is neglectable when the mesh is fine enough.
Remark 5.1. The numerical tests here show that the convergence rate of the flux-only LSFEMs is the same as the original flux-solution LSFEMs that we suggested in [13].
As mentioned in [13], the common rumor that of the least-squares method is it tends to have a strong smearing effect which is actually often the effect of Galerkin least-squares or stabilized methods where some least-squares terms are added to variational problems, see for example  [11]. We do not observe smearing effect for the bona fide least-squares methods developed in this paper.

Concluding Remarks
In this paper, two flux-only least-squares finite element methods (LSFEM) for the linear hyperbolic transport problem are developed. We first reformulate the linear transport equation into a flux-solution system, then eliminate the solution from the system and get two flux-only formulations, and develop corresponding LSFEMs. The solution then is recovered by simple post-processing methods using its relation with the flux. These two versions of flux-only LS-FEMs use less DOFs than the method we developed in [13].
In this paper, we only discuss the versions of LSFEMs with boundary condition strongly enforced. We can also weakly enforce the condition as we did in [13].
Similar to the LSFEM developed in [13], both flux-only LSFEMs can handle discontinuous solutions better than the traditional continuous polynomial approximations. We show the existence, uniqueness, a priori and a posteriori error estimates of the proposed methods. With adaptive mesh refinements driven by the least-squares a posteriori error estimators, the solution can be accurately approximated even when the mesh is not aligned with discontinuity. The overshooting phenomenon is very mild if a piecewise constant reconstruction of the solution is used. Extensive numerical tests are done to show the effectiveness of the methods developed in the paper.