{'fraction_non_alphanumeric': 0.04544675012862288, 'fraction_numerical': 0.01800720288115246, 'mean_word_length': 4.303547066848568, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Two-Level Finite Element Discretization of the Streamfunction Formulation of the Stationary Quasi-Geostrophic Equations of the Ocean

In this paper we proposed a two-level finite element discretization of the nonlinear stationary quasi-geostrophic equations, which model the wind driven large scale ocean circulation. Optimal error estimates for the two-level finite element discretization were derived. Numerical experiments for the two-level algorithm with the Argyris finite element were also carried out. The numerical results verified the theoretical error estimates and showed that, for the appropriate scaling between the coarse and fine mesh sizes, the two-level algorithm significantly decreases the computational time of the standard one-level algorithm.


Introduction
Two-level algorithms are computationally efficient approaches for finite element (FE) discretizations of nonlinear partial differential equations [6,7,17,30,44]. A two-level FE discretization aims to solve a particular nonlinear elliptic equation by first solving the nonlinear system on a coarse mesh and then using the coarse mesh solution to solve the linearized system on a fine mesh. The appeal of such a method is clear; one need only solve the nonlinear equations on a coarse mesh and then use this solution to solve on a fine mesh, thereby reducing computational time without sacrificing solution accuracy. The development of the two-level FE discretization was originally performed by Xu in [44]. Later algorithms were developed for the Navier-Stokes equations (NSE) by Layton [30] (see also [17,18,39,31,45,32]) and for the Boussinesq equations by Lenferink [34].
As computational power increases, complex models are becoming more and more popular for the numerical simulation of oceanic and atmospheric flows. Computational efficiency, however, remains an important consideration for geophysical flows in which long time integration is needed. Thus, simplified mathematical models are central to the numerical simulation of such flows. For example, the quasigeostrophic equations (QGE), a standard mathematical model for wind driven large scale oceanic and atmospheric flows [35,42], are often used in climate modeling [14].
Most FE discretizations of the QGE are for the streamfunction-vorticity formulation. The reason is that the streamfunction-vorticity formulation allows the use of low order (C 0 ) FEs, although one needs to discretize two flow variables, the potential vorticity, q, and the streamfunction, ψ. We note that the streamfunction-vorticity formulation is often used in the numerical discretization of the 2D NSE, to which the QGE are similar in form. Alternatively, one can, instead, use the pure streamfunction formulation of the QGE. The advantage lies in an equation that contains only one flow variable, the streamfunction, ψ, at the price of having to deal with a fourth-order partial differential equation. Thus, the numerical discretization of the pure streamfunction formulation of the QGE with conforming FEs requires the use of high-order (C 1 ) FEs, e.g., the Argyris element [1,9].
The streamfunction formulation of the QGE still suffers from having to solve a large nonlinear system of equations. This is usually done by using a nonlinear solver, such as Newton's method. These nonlinear solvers typically require solving large linear systems multiple times to obtain the solution to the nonlinear system. Solving these large linear systems multiple times can be time consuming. Thus, a two-level algorithm can significantly reduce computational time over the standard nonlinear solver, since we need only solve the nonlinear system on a coarse mesh and then use that solution to solve a linear system on a fine mesh.
In this paper, we propose a two-level algorithm for the FE discretization of the streamfunction formulation of the stationary QGE (SQGE). Just as in the NSE case [25,40,33], we regard the stationary QGE as a stepping-stone to the time-dependent QGE, which are the ultimate goal of the two-level algorithm put forth in this report. The conforming FE discretization is based on the Argyris element. Additionally, we present a rigorous error analysis for the two-level FE discretization. The theoretical error bounds as well as the increased computational efficiency are illustrated numerically for two test problems.
The rest of the paper is organized as follows: In Section 2 we present the SQGE. In Section 3 we present the weak formulation of the SQGE, including notation and functional spaces. Section 4 contains the presentation of the one-level FE discretization of the SQGE. In Section 5, we discuss both the twolevel algorithm and its application to the SQGE. Next, in Section 6 we provide rigorous error bounds for the two-level FE discretization of the SQGE and we discuss the scaling between the fine mesh size, h, and the coarse mesh size, H. Section 7 includes numerical results which both verify the theoretical error bounds presented in Section 6 and illustrate the computational efficiency of the two-level algorithm over the standard one-level method. Finally, in Section 8 we present our conclusions.

Streamfunction Formulation
The SQGE in a simply connected domain Ω is where are the Jacobian, Rossby number, Reynolds number, respectively, and β, A, U , and L are the coefficient in the beta plane approximation, the eddy viscosity, the characteristic velocity scale, and the characteristic length scale, respectively (see [21,38]). To completely specify (1), we need to impose boundary conditions (see [13,42,37] for a careful discussion of this issue). In this report, we consider where n represents the outward unit normal to Ω. These are also the boundary conditions used in [26,17,19] for the 2D NSE.

Weak Formulation
Now we can derive the weak formulation of the SQGE (1). To this end, we first introduce the appropriate functional setting. Let Multiplying (1) by a test function χ ∈ X and using the divergence theorem, we get, in the standard way (see [26]), the weak formulation of the SQGE: Find ψ ∈ X such that a(ψ, χ) + b(ψ; ψ, χ) + c(ψ, χ) = (χ), ∀χ ∈ X, where We note that in the space H 2 0 the semi-norm | · | 2 and the norm · 2 are equivalent (see (1.2.8) in [12]).
Lemma 1. Given ψ, ξ, ϕ ∈ H 2 0 (Ω) and F ∈ H −2 (Ω), the linear form , the bilinear forms a and c, and the trilinear form b are continuous: there exist Γ 1 > 0 and Γ 2 > 0 such that For a proof, see [11]. For small enough data, one can use the same type of arguments as those used in Chapter 6 in [33] (see also [24,25]) to prove that the SQGE (1) are well-posed [3,43]. In what follows, we will always assume that the small data condition involving Re, Ro and F , is satisfied and, thus, that there exists a unique solution ψ to (1).
The following stability estimate was proven in Theorem 1 in [21]: The solution ψ of (1) satisifies the following stability estimate:

Finite Element Formulation
Let T H denote a FE triangulation of Ω with meshsize (maximum triangle diameter) H. We consider a conforming FE discretization of (6), i.e., X H ⊂ X = H 2 0 (Ω). The FE discretization of the SQGE (6) reads: Find ψ H ∈ X H such that Using standard arguments [24,25], one can prove that, if the small data condition used in proving the well-posedness result for the continuous case holds, then (13) has a unique solution ψ H (see Theorem 2.1 and subsequent discussion in [11]). Furthermore, one can prove the following stability result for ψ H using the same arguments as those used in the proof of Lemma 2 for the continuous setting (see Theorem 2 in [21]).
Lemma 3. The solution ψ H of (13) satisfies the following stability estimate: As noted in Section 6.1 in [12] (see also Section 13.2 in [26], Section 3.1 in [29], and Theorem 5.2 in [8]), in order to develop a conforming FE for the SQGE (6), we are faced with the problem of constructing subspaces of H 2 0 (Ω). Since the standard, piecewise polynomial FE spaces are locally regular, this construction amounts in practice to finding FE spaces X H that satisfy the inclusion X H ⊂ C 1 (Ω), i.e., finding C 1 FEs. Several FEs meet this requirement (see, e.g., Section 6.1 in [12], Section 13.2 in [26], and Section 5 in [8]): the Argyris triangular element, the Bell triangular element, the Hsieh-Clough-Tocher triangular element (a macroelement), and the Bogner-Fox-Schmit rectangular element. We note that any C 1 FE can be used for this study. We chose the Argyris element, however, because it is a triangle, which allows for easy treatment of complex boundaries, and because of the recent development of a transformation that allows for calculations to be done on the reference element [15]. We also note that, in addition to the use of C 1 conforming FEs one might use the newly developed method of isogeometric analysis as in [2].

Two-Level Algorithm
In this section we propose a two-level FE discretization of the SQGE (6). We let X h , X H ⊂ X = H 2 0 (Ω) denote two conforming FE meshes with H > h. The two-level algorithm consists of two steps. In the first step, the nonlinear system is solved on a coarse mesh, with mesh size H. In the second step, the nonlinear system is linearized around the approximation found in the first step, and the resulting linear system is solved on the fine mesh, with mesh size h. This procedure is as follows:

Algorithm 1 Two-Level algorithm
Step 1: Solve the following nonlinear system on a coarse mesh for ψ H ∈ X H : Step 2: Solve the following linear system on a fine mesh for ψ h ∈ X h : The well-posedness of the nonlinear system (15) was proven in [11] (see also [21]). The following error estimate for the approximation in Step 1 of the two-level algorithm (Algorithm 1) was proven in Theorem 2 in [21]: Theorem 1. Assume that the following small data condition is satisfied: Let ψ be the solution of (6) and ψ H be the solution of (15). Then the following error estimate holds: where The following lemma proves the well-posedness of the linear system (16): Lemma 4. Given a solution ψ H of (15), the solution ψ h of (16) exists uniquely.
Proof. First, we introduce the bilinear form B : X h × X h → R given by Lemma 1 yields the following inequality: The stability estimate for ψ H in Lemma 3 and inequality (19) imply that B is continuous. Additionally, the fact that b(ψ H ; ψ h , ψ h ) = 0 and c(ψ h , ψ h ) = 0 for all ψ h ∈ X h combined with the Poincaré-Friedrichs inequality gives Thus, B is coercive. Therefore, by the Lax-Milgram lemma, ψ h exists and is unique.
In addition to the existence, uniqueness, and stability of the solution of the continuous linear system (Lemma 4) we also have a stability bound for the solution on the discrete fine mesh, h.
Lemma 5. The solution ψ h of (16) satisfies the following stability bound: Proof. Setting χ h = ψ h in (16), and noting that c(ψ h , where in the last inequality we used (11). Therefore, it follows that

Error Bounds
The main goal of this section is to develop a rigorous numerical analysis for the two-level algorithm (Algorithm 1). The proof for the error bounds follows the pattern used in [17].
We first introduce an improved bound on the trilinear form b(ζ; ξ, χ). To this end, we use the following inverse inequality, proven in Lemma 6.4 in [41] (see also page 122 in [17]): The following lemmas will be useful in determining the error bounds for Step 2 of the two-level algorithm. The first lemma, which corresponds to Lemma 5.1 in [17], follows from (21) and (9) and places error bounds on the trilinear form b(ψ; χ h , ξ): Lemma 6. For any χ h ∈ X h , the following inequality holds: The following lemma, which corresponds to Lemma 5.6 in [17], will be useful for proving the error bounds for Algorithm 1, by allowing one to permute the terms of the trilinear form: The following theorem gives the error bound after Step 2 of the two-level algorithm (Algorithm 1) and is the main result of this paper. The proof of this theorem is similar to the proof of Theorem 5.2 in [17].
In what follows, we consider both X h and X H Argyris FE spaces. We emphasize, however, that both Algorithm 1 and the error estimate in Theorem 2 remain valid for other conforming FE spaces, e.g., the Bell element, the Hsieh-Clough-Tocher element, or the Bogner-Fox-Schmit element.
For the Argyris triangle, we have the following inequalities, which follow from approximation theory [5] and Theorem 6.1.1 in [12]: where j = 0, 1, 2 and ψ, the solution of (6), is assumed to satisfy ψ ∈ H 6 (Ω) ∩ H 2 0 (Ω). Of course, if the solution ψ does not possess the assumed regularity, the convergence rates of the FE discretization can significantly degrade. This behavior is standard for high-order numerical methods. Corollary 1. Let X h , X H ⊂ H 2 0 (Ω) be Argyris finite elements. Then ψ h , the solution of the two-level algorithm (Algorithm 1) satisfies the following error estimate: Proof. This follows directly by substituting the inequalities (33) into (24).

Numerical Results
The goal of this section is twofold: first, we illustrate the computational efficiency of the two-level method, and second, we verify the theoretical rates of convergence developed in Section 6. To illustrate the computational efficiency of the two-level method, we compare solution times for the full nonlinear one-level method and for the two-level method applied to the SQGE. We choose coarse mesh/fine mesh pairs such that the ratio is 1 /2. To verify the theoretical rates of convergence, we compare the theoretical error estimates to the observed rates of convergence from our numerical tests. The calculations on the coarse mesh employ our original code that was benchmarked in [21].
The FE solver was written in MATLAB (the 2010b version), with the Argyris implementation written in C. The nonlinear systems were solved with Newton's method; at each step the resulting linear systems were solved with UMFPACK. The numerical tests were run on a Mac Pro with 16 gigabytes of RAM and two quad-core Intel Xeon processors.
The homogeneous boundary conditions are ψ = ∂ψ ∂ n = 0 and the forcing function F corresponds to the exact solution (35).

Practical Considerations
A key part of two-level algorithms is accessing a previous coarse mesh solution, i.e., finding the parent element given a child element. This step can negate any performance benefits if not implemented wisely. Indeed, let n be the number of elements in the FE discretization. For the unit square, a naïve search across every element takes O(n/2) operations. This procedure may be improved with a binary search, which is summarized in Algorithm 2.
We note that every element on the fine mesh corresponds to exactly one element on the coarse mesh. However, a coarse mesh element may correspond to multiple elements on the fine mesh.
Algorithm 2 Given an element on the fine mesh, determine the parent element on the coarse mesh. Before examining the fine mesh, sort the coarse mesh elements by their centroid values.
Step 1: Select an element on the fine mesh and compute its centroid.
Step 2: Perform a binary search across the coarse mesh elements until the difference between the x-values of the fine mesh centroid and coarse mesh centroids is less than H, the coarse mesh step size. There should be many elements that fit this condition; save them as a list.
Step 3: Search through this list until we find the correct coarse mesh element (that is, the centroid of the fine-mesh element is an interior point of the correct coarse mesh element).
For the considered unit square, the binary search will examine on average log(n) elements, while the linear search component involves at most √ n/2 elements. Therefore the search requires a O( √ n/2) number of element checks. Profiling results indicate that using Algorithm 2 to identify parent elements takes much less time than either setting up or solving the systems, so this approach is fast enough that lookup time does not contribute significantly to overall solution time.

Computational Efficiency
To illustrate the computational efficiency of the two-level method, we compare the simulation time for the standard one-level method (i.e., the full nonlinear system, without the two-level method) with the simulation time for the two-level method.
In Table 1, the L 2 -norm of the error (e L 2 ), the H 1 -norm of the error (e H 1 ), the H 2 -norm of the error (e H 2 ) and the simulation times are listed for various mesh sizes. For each fine mesh, we choose a coarse mesh that ensures the same order of magnitude for the errors in the one-level and two-level methods.
For small values of the fine mesh size, h, the two-level method was significantly faster than the one-level method. The errors in the H 2 -norm were nearly identical, while the error in the L 2 -norm were generally of the same order of magnitude. We also note that the tolerance in Newton's method seems to cause a plateau in the L 2 -norm of the error. The results in Table 1 are illustrated graphically in Figure 1. In this figure the simulation times of the one-level method (green) and of the two-level method (blue) are displayed for all the pairs (h, H) in Table 1. Figure 1 clearly shows that as the number of degrees of freedom (DoFs) increases, the computational efficiency of the two-level method increases as well.

Rates of Convergence
The goal of this subsection is to numerically verify the theoretical rates of convergence in estimate (34) of Corollary 1. Unlike the theoretical error estimates for the one-level method developed in [21], for the two-level method we must verify rates of convergence for two different mesh sizes: the fine mesh size, h, and the coarse mesh size, H.
To numerically verify the theoretical rate of convergence given in estimate (34) with respect to H, we fix h to a small value and we vary H. Thus, the total error in estimate (34) will be dominated by the H term, i.e., the total error will be of order O(H 5 ). In Table 2, we fix h = 0.0063 and we vary H. The error in the L 2 -norm (e L 2 ), the error in the H 1 -norm (e H 1 ), the error in the H 2 -norm (e H 2 ), and the rate of convergence with respect to H are listed in Table 2. The rate of convergence with respect to H of the error in the H 2 -norm follows the theoretical rate predicted in estimate (34) (i.e., fifth-order). For the last mesh pair, however, the rate of convergence appears to drop off. This occurs because, for small values of H, the total error in estimate (34) is not dominated anymore by the H term. We note that in Table 2 the number of DoF s for the fine mesh varies even though the mesh size h is constant. The reason for this is that for the coarse mesh a Delauney triangulation is performed around a grid spacing corresponding to H and then a red refinement is performed to obtain the fine mesh corresponding to     To numerically verify the theoretical rate of convergence given in estimate (34) with respect to h, we must proceed with caution. The reason is that a straightforward approach would fix H and let h go to zero. This approach, however, would fail, since the H term would eventually dominate the total error. To avoid this, we consider the following scaling between the mesh sizes: where C > 1. The scaling in (36) implies that the total error in estimate (34) is of order O(h 4 ). Indeed, the second term on the right hand side of estimate (34) now scales as follows: where in the last relation in (37) we used the fact that | ln(h)| h → 0 when h → 0 (which follows from l'Hospital's rule).

Remark 1.
We emphasize that the scaling in (36) is not needed in the two-level algorithm. We only use it in this subsection to monitor the convergence rate with respect to h.
In this subsection, we consider C = 2 in (36). We note, however, that any other constant C > 1 could be used in (36). With this choice, we are now ready to numerically verify the theoretical rate of convergence given in estimate (34) with respect to h, which, as shown in (37), will be of order O(h 4 ). In Table 3, for various mesh size pairs (H = 2h, h), we list the L 2 -norm of the error (e L 2 ), the H 1 -norm of the error (e H 1 ), the H 2 -norm of the error (e H 2 ), and the rate of convergence. The rate of convergence with respect to h of the error in the H 2 -norm follows the theoretical rate predicted in estimate (34) (i.e., fourth-order).
To investigate the performance of the two-level method in a more realistic setting, we consider the same setting as that used in Test 2 in [10]. The exact solution is given by where the computational domain is the [0, 3] × [0, 1] rectangle. The exact solution displays a sharp boundary layer, which is similar to the western boundary layer encountered in the large scale ocean circulation in the northern hemisphere [10,42,21]. Furthermore, we use the parameter values Re = 5 and Ro = 10 −4 , which are similar to the realistic values used in [4,22,23,21] for the Mediterranean Sea. In Table 4, for various mesh size pairs (H ≈ 2h, h), we list the L 2 -norm of the error (e L 2 ), the H 1 -norm of the error (e H 1 ), the H 2 -norm of the error (e H 2 ), and the rate of convergence. The rate of convergence with respect to h of the error in the H 2 -norm follows the theoretical rate predicted in estimate (34) Table 4: Two-level method: the L 2 -norm of the error (e L 2 ), the H 2 -norm of the error (e H 2 ), and the convergence rate with respect to h.

Conclusions
In this paper, we proposed a two-level FE discretization of the (nonlinear) stationary QGE in the pure streamfunction formulation. The two-level algorithm consists of two steps. In the first step, the nonlinear system is solved on a coarse mesh. In the second step, the nonlinear system is linearized around the approximation found in the first step, and the resulting linear system is solved on the fine mesh.
Rigorous error estimates for the two-level FE discretization were derived. These estimates are optimal in the following sense: for an appropriately chosen scaling between the coarse mesh size, H, and the fine mesh size, h, the error in the two-level method is of the same order as the error in the standard one-level method (i.e., solving the nonlinear system directly on the fine mesh).
Numerical experiments for the two-level algorithm with the Argyris element were also carried out. The numerical results verified the theoretical error estimates, both with respect to the coarse mesh size, H, and the fine mesh size, h. Furthermore, the numerical results showed that, for an appropriate scaling between the coarse and fine mesh sizes, the two-level method significantly decreases the computational time of the standard one-level method.
We plan to extend this study in several directions. We will treat the case of multiply connected domains [36,26,27,28] in order to allow the numerical investigation of more realistic computational domains (e.g., islands in the Mediterranean Sea and in the North Atlantic). We will also consider the time-dependent QGE and the two-layer QGE (which will allow the study of stratification effects). Finally, we plan to investigate various preconditioning techniques to improve the performance of the linear solvers used in this report (see [16] for the NSE in the primitive variable formulation and [20] for the NSE in the streamfunction-vorticity formulation).