Data assimilation finite element method for the linearized Navier-Stokes equations with higher order polynomial approximation

In this article, we design and analyze an arbitrary-order stabilized finite element method to approximate the unique continuation problem for laminar steady flow described by the linearized incompressible Navier--Stokes equation. We derive quantitative local error estimates for the velocity, which account for noise level and polynomial degree, using the stability of the continuous problem in the form of a conditional stability estimate. Numerical examples illustrate the performances of the method with respect to the polynomial order and perturbations in the data. We observe that the higher order polynomials may be efficient for ill-posed problems, but are also more sensitive for problems with poor stability due to the ill-conditioning of the system.


Introduction
The question of how to assimilate measured data into large-scale computations of flow problems is receiving increasing attention from the computational mathematics community [3,7,22,26,28,36].There are several different situations where such data assimilation problems as can be seen in the above examples.One situation is when the data necessary to make the flow problem well-posed is lacking, for instance, when the data on the boundary of the domain is unknown; instead, measurements are available in some subset of the bulk domain or boundary to make up for this shortfall.In such a case, the problem is typically ill-posed, and numerical simulations are significantly more challenging to perform than when handling well-posed flow problems.Ill-posed problems usually come up in inverse problems and data assimilation.Traditionally, these ill-posed problems have been solved by regularizing at the continuous level, using e.g.Tikhonov regularization [38] or quasi-reversibility [34].The regularized problem is well-posed and may be discretized using any appropriate numerical technique.Then, the regularization parameter must be tuned to the optimal value for the noise in the data.There is considerable literature of research on Tikhonov regularization and inverse problems, and we suggest the reader to [32] and its references for an overview of computational approaches employing this strategy.The quasireversibility methods relevant to the current study may be found in [10][11][12]21].
The goal of the current contribution is to develop a finite element approach directly applied to the ill-posed variational data assimilation form.Regularization is then introduced at the discrete level utilizing stabilized finite element methods that allow for a comprehensive analysis employing conditional stability estimates.The idea is presented in [13] for standard  1 -conforming finite element methods.Ill-posed problems are analyzed in [14], and in [16], the technique is extended to nonconforming approximations.In both cases, low-order approximation spaces are considered.The error analysis requires the availability of sharp conditional stability estimates for the continuous problem.The estimates are conditional in the sense that a particular a priori bound must be assumed to hold for the solution, and the continuity provided in this bound is often merely Hölder [33].In the literature, these estimates are referred to as quantitative uniqueness results and employ theoretical methods such as Carleman estimates or three-ball estimates [1,31].Error bounds derived using conditional stability estimates can be optimal because they reflect the approximation order of the finite element space and the stability of the ill-posed problem.In particular, when applied to a well-posed problem, the finite element method recovers optimal convergence.
The ill-posed problem that we consider here is the unique continuation problem.The unique continuation problem for the Stokes equations was initially studied in [25].The analysis of the stability properties of ill-posed problems based on the Navier-Stokes equations is a very active field of research, and we refer to the works [4-6, 8, 29, 30, 35] for recent results.
This study aims to determine whether using high-order methods in the primal-dual stabilized Galerkin methods is as helpful in the ill-posed case as in the well-posed situation.Inspired by the approach proposed in [9] for the lowest-order finite element discretization of the unique continuation problem subject to the Navier-Stokes equations, here we generalize the method to arbitrary polynomial orders and investigate the benefits of using higher-order polynomials in numerical experiments.
The rest of the paper is organized as follows.In Section 2, we introduce the considered inverse problem and some related stability estimates.In Section 3, we describe the proposed stabilized finite element approximation of the data assimilation problem and state the local error estimate.The numerical analysis of the method is carried out in Section 4. Finally, Section 5 presents a series of numerical examples which illustrate the performance of the proposed method.

The linearized Navier-Stokes problem
Let Ω be an open polygonal (polyhedral) domain in R  ,  = 2, 3. Let (,  ) be the solution of the stationary incompressible Navier-Stokes equations and consider some perturbation (, ) of this base flow.If the quadratic term is ignored, the linearized Navier-Stokes equations for (, ) can be written where Here,  is a diffusion coefficient.We assume that  belongs to [ 1,∞ (Ω)]  and that (, ) satisfies the regularity For this problem, we assume that measurements on  are available in some subdomain   ⊂ Ω having a nonempty interior and our purpose is to reconstruct a fluid flow perturbation of  for system (1)-(2) based on the measurements of velocity.Now, we will present some useful notations.Consider the following spaces: where  2 0 (Ω) = { ∈  2 (Ω) : ∫︀ Ω  = 0}.We also define the norms, for  = 1 or , Observe that in the definitions, we employ the same notation for  = 1 and  = .For any subdomain  ⊂ Ω, we set Next, define the bilinear forms as: for all (, ) where  :  := ∑︀  ,=1  ,  , and, for all (, ) The weak form of the inverse problem can be expressed as: and Here,  ∈ [ 1 (  )]  corresponds to the exact fluid velocity on   , i.e.  is a solution to the linearized Navier-Stokes' equations in   and has an extension  to all of Ω. Below in the finite element method we will assume that we do not have access to , but only some measured velocities   =  + .So   corresponds to the exact velocity polluted by a small noise  ∈ [ 2 (  )]  .Consider the linearized Navier-Stokes problem with a non-zero velocity divergence We assume that if the boundary conditions of system ( 8)-( 9) are homogeneous Dirichlet boundary conditions, then it is well-posed.More precisely, we make the following assumption: Assumption 2.1.For all  ∈  ′ 0 and  ∈  0 we assume that system (8)-( 9) admits a unique weak solution (, ) ∈  0 ×  0 and that there exists a constant   > 0 depending only on ,  and Ω such that × is small enough, then the Lax-Milgram lemma implies that Assumption 2.1 holds.The assumption of smallness on ∇ is a sufficient condition, there are reasons to believe that Assumption 2.1 holds in more general cases.In the homogeneous case (which corresponds to  = 0 in (1)- (2) or to  = 0 and  = 0 in ( 8)-( 9)), a solution (, ) satisfies a three-balls inequality which only involves the  2 norm of the velocity.This three-balls inequality result is stated in [35] (with their notations,  corresponds to  and  to ∇ ).Theorem 2.2.(Conditional stability for the linearized Navier-Stokes problem).Let  ∈  ′ 0 ,   ⊂ Ω and  ∈  be given.For all  ⊂⊂ Ω, there exist  > 0 and 0 <  < 1 such that for all (, ) ∈ [ 1 (Ω)]  ×  1 (Ω) solution of ( 8)- (9).
Proof.For the proof we refer the reader to [9, Appendix A].Theorem 2.2 provides a conditional stability result for ill-posed problems [1] in the sense that, for this estimate to be helpful, it must be accompanied by an a priori bound on the solution on the global domain (due to the presence of ‖‖  on the right-hand side).Specifically, Theorem 2.2 implies that a solution (, ) in [ 1 (Ω)]  ×  1 (Ω) of problem ( 6) and (7), must be unique.For the pressure uniqueness holds up to a constant.Moreover, in inequality (11), the exponent  depends on the dimension , the size of the measure domain   and the distance between the target domain  and the boundary of the computational domain Ω.

Stabilized finite element approximation
In this section, we first introduce a discretization of problem (13) using a standard finite element method.Then, the discrete inverse problem is reformulated as a constrained minimization problem in the discrete space where the regularization of the cost functional is achieved through stabilization terms.Finally, the estimation of the error between the exact continuous solution and the discrete solution of our minimization problem is stated in Theorem 4.12 which corresponds to our main theoretical result.
Let {T ℎ } ℎ be a family of affine, simplicial meshes of Ω.For simplicity, the family {T ℎ } ℎ is supposed to be quasi-uniform.Mesh faces are collected in the set F ℎ which is split into the set of interior faces, F int ℎ , and of boundary faces, F ext ℎ .For a smooth enough function  that is possibly double-valued at  ∈ F int ℎ with  =  − ∩  + , we define its jump at  as [] =:   − −   + , and we fix the unit normal vector to  , denoted by   , as pointing from  − to  + .The arbitrariness in the sign of [𝑣] is irrelevant in what follows.
We next define a piecewise polynomial space as where P  ( ),  ≥ 0, is the space of polynomials of degree at most  over the element  .Further, define a conforming finite element space as For the analysis below the polynomial degrees of the above spaces may be chosen as  ≥ 1,  1 ≥ 1,  2 ∈ {max{1,  − 1}, } and  3 ≥ 1 and the convergence order will be given in terms of .To make the notation more compact we introduce the composite spaces V ℎ :=   ℎ ×  0 ℎ and W ℎ :=  ℎ ×  ℎ .We may then write the finite element approximation of (13): Find ( ℎ ,  ℎ ) ∈ V ℎ such that (( ℎ ,  ℎ ), ( ℎ ,  ℎ )) = ( ℎ ), (14) for all ( ℎ ,  ℎ ) ∈ W ℎ .Let us introduce the measurement bilinear form to take into account the measurements on   given by (6).
where  = max(, ‖ ‖ [ ∞ (Ω)] × ℎ) and   > 0 will correspond to a free parameter representing the relative confidence in the measurements.The objective is then to minimize the functional under the constraint that ( ℎ ,  ℎ ) satisfies (14).
We now introduce the following discrete Lagrangian for ( If we differentiate with respect to ( ℎ ,  ℎ ) and ( ℎ ,  ℎ ), we get the following optimality system: for all ( ℎ ,  ℎ ) ∈ V ℎ and ( ℎ ,  ℎ ) ∈ W ℎ .However, the discrete Lagrangian associated to this problem leads to an optimality system which is ill-posed.To regularize it, we introduce stabilization operators that will convexify the problem with respect to the direct variables  ℎ ,  ℎ and the adjoint variables  ℎ ,  ℎ .We introduce The choice of stabilization terms will be discussed later.For compactness, we introduce the primal and dual stabilizers: for all ( ℎ ,  ℎ ), ( where and   , ,   , and  div are positive user-defined parameters.And for all ( ℎ ,  ℎ ), ( where  *  and  *  are positive user-defined parameters.Let us make some comments on these stabilization terms.The stabilization of the direct velocity acts on fluctuations of the discrete solution through a penalty on the jump of the solution gradient over element faces and has no equivalent on the continuous level.The form   (•, •) is a Galerkin least squares stabilization.Let us mention that there is some freedom in the choice of dual stabilization, e.g.set  *  ( ℎ ,  ℎ ) =  *  ∫︀ Ω ∇ ℎ ∇ ℎ x.We will only detail the analysis for the first choice (23) below.We refer the reader to [13,15] for a more general discussion of the possible stabilization operators.

Stability and Error Analysis
To prove the stability of our formulations, we need the following result.
Lemma 4.1.There exists   such that for all  ℎ ∈  ℎ there holds Proof.The following Poincaré inequality is well known Lemma B.63 of [24].If  :  1 (Ω) → R is a linear functional that is non-zero for constant functions then For instance, we may take As an immediate consequence we have the bound (27).
Let us prove that the discrete problem is well-posed.We can write the discrete formulation in a more compact form.
Using the above bounds to the componentwise extension of  ℎ to vectorial functions, we deduce the following approximation bound.
In this section, we will present and prove several technical results.First observe that the formulation ( 25)-( 26) is weakly consistent in the sense that we have a modified Galerkin orthogonality relation with respect to the scalar product associated to : Lemma 4.7.(Consistency).Let (, ) satisfy ( 1) and ( ℎ ,  ℎ ) be a solution of ( 25)- (26).Then there holds Proof.The result follows by taking the difference between ( 13) and (25).
Proof.The approximation bounds can be deduced using the component-wise extension of  ℎ to vector functions.
Proof.Let us derive the estimate (41).Using the definition of (•, •) : Consider the first term on the right hand side of (42).Using the Cauchy-Schwarz inequality and Poincaré inequality, The second term of (42) can be handled as: The last term of (42) can be handled as: Finally, the result follows by combining all the above estimates.
Lemma 4.10.We assume that the solution (, ) ∈ [ +1 (Ω)]  ×  2 0 ∩   (Ω) and we consider ( ℎ ,  ℎ ) ∈ V ℎ and ( ℎ ,  ℎ ) ∈ W ℎ the discrete solution of ( 25)- (26).Then there holds, Proof.We introduce the discrete errors The first term of (44) can be handled by using the Lemma 4.8.Consider the second term of ( 44) To estimate the right-hand side, we notice that, using the second equation of ( 26) with ( ℎ ,  ℎ ) = ( ℎ ,  ℎ ) Using Lemma 4.7, we obtained Adding ( 45) and (46), We bound the terms (1)-( 3) term by term.The first term is handled by using Lemma 4.9 Consider the second term on the right hand side of (47) We now estimate the terms on the right hand side of (49).Using the  1 -stability of  ℎ , the first term of (49) can be handled as: .
The last term can be handled as: Put together (47) leads to The combination of the above estimates concludes the claim.
Corollary 4.11.Under the same assumptions as for Lemma 4.10, there holds and Proof.Using Lemma 4.10, we see that The estimate (51) is immediate by using the triangle inequality and the estimate (50).
The following theorem is the main theoretical result of the paper and states an error estimate for this method.
We introduce  ℎ and  ℎ being fixed.The linear forms   and   on  0 and  respectively defined by: For all  ∈  0 and  ∈  ⟨  , ⟩  ′ 0 , + (  , )  := (, )  − (( ℎ ,  ℎ ), (, )). (53) It follows that (, ) is the solution of ( 8)-( 9) with  and  in the right hand sides replaced respectively by   and   .Applying now Corollary 2.2, we directly get Using (25), we can write the residuals: for all ( ℎ ,  ℎ ) We take  ℎ =  ℎ  and  ℎ =  ℎ  in (55).Now, let us estimate the terms on the right hand side of (55).Consider the first two terms of (55) Applying an integration by parts to the first two terms of (56) and using Lemma 4.10, The last term is handled using the Cauchy-Schwarz inequality and Lemma 4.10, Applying the above bounds in (56) leads to The last term of (55) is handled using Lemma 4.10 As a consequence we can bound the quantity defined in (55) by Since this bound holds for all  ∈  0 and  ∈ , we conclude that Using the Poincaré inequality ( 27), we have the bound Thus, we can bound the terms in the right-hand side of (54) in the following way: and Using these two bounds in (54), we conclude that which completes the proof of the theorem.

Numerical simulations
In this section, we use several two-dimensional numerical examples to apply the methodology described in Section 3.All experiments have been implemented using the open-source computing platform FEniCSx [2,37].A docker image to reproduce the numerical results is available at https://doi.org/10.5281/zenodo.7442458.The free parameters in ( 25)-( 26) are set to In all the numerical examples.In the first example we will verify the convergence orders for different polynomial orders using equal order interpolation  for all variables.Then we consider the same test case using the minimal polynomial order that results in the same error bounds,  1 = 1,  2 = max{ − 1, 1} and  3 = 1.Finally, we study the robustness of the error estimate with respect to the viscosity for a configuration where the target subdomain  is strictly downwind the data subdomain   , so that every point  is on a streamline intersecting   .

Convergence study: Stokes example
To demonstrate the convergence behaviour of the method introduced in Section 3, we take the test case for the Stokes problem from [17].Let Ω = [0, 1] 2 be the unit square.We consider the velocity and pressure fields given by It is simple to demonstrate that (, ) is a solution to the homogeneous Stokes problem with  = 1, corresponding to the system (1)-( 2) with  = 0 and  = 0.As a result, we consider ( 25)-( 26) with  = 0 and  = 0. Two different geometric settings are considered: one in which the data is continued in the convex geometry, inside the convex hull of   , and one in which the solution is continued in the non-convex geometry, outside the convex hull of   .The convex geometry represented by Figure 1a is given as: We begin by performing the computation using unperturbed data.The relative  2 -norm errors ‖ −  ℎ ‖ [ 2 ()]  / ‖‖ [ 2 ()]  are computed in the subdomain .In addition, we present the history of convergence of the residual quantity for velocity stabilization: .
Figure 2 displays the velocity, pressure errors and residual quantity in the convex and non-convex geometry.Filled squares, circles and triangles represent the velocity errors; dashed lines represent the pressure error, and the plain thin lines represent the residual.The expected order of convergence is observed for the residual in Lemma 4.10.The local velocity error behaves consistently with the convergence rates obtained in Theorem 4.12.We can also see in Figure 2 that the higher order polynomials are more satisfactory for ill-posed problems.Next, we proceed with the numerical verification of the above method with data perturbation.Consider the perturbed data   = |   + , with random perturbations for some  ∈ N 0 available for implementing our method.According to Theorem 4.12, we have the estimate  consequently, convergence requires the condition  −  > 0. Figures 3-4 present the convergence history of the velocity, pressure and residual quantities with the data perturbation in the convex and non-convex geometry, respectively.The effect of different values of  for relative  2 -error are studied in Figures 3-4.The relative error for  = 0 is displayed in Figures 3a and 4a.In both cases, the results are in agreement with the Theorem 4.12.
As stated in (65), the  = 1 polynomial approximation may diverge for  = 1, which is confirmed by Figure 3b.
Next, the method  = 2 converges linearly, whereas  = 3 still manages to converge, albeit at a slower rate.As shown in the Figure 3c, this result is consistent with the fact that for  = 2, convergence is no longer observed for any  ≤ 3. Similar convergence results are observed in the non-convex domain as shown in Figure 4.The results of Figures 3-4 indicate that for the convex geometry  ≈ 1 and for the non-convex geometry  ≈ 2 3 .In Figures 5-7 the same results are presented for the case where the minimal polynomial order is considered that is  1 = 1,  2 = max{ − 1, 1},  3 = 1.The results are very similar and we conclude that for these numerical examples there is no disadvantage in taking the smallest possible dual space.

Convergence study with varying viscosity
In this subsection, we consider the flow of a viscous Newtonian fluid between two solid boundaries at  = , − driven by a constant pressure gradient.The source term  is chosen such that the solution of the plane (a) As in the previous section, we have examined the convergence of the method by performing numerical tests on both unperturbed and perturbed data.We vary the viscosity between  = 1 and  = 0. Observe that since no boundary conditions are imposed nothing needs to be changed in the formulation in the singular limit.Also note that the choice of   and   in ( 20)-( 21) mimicks the choice for the stabilized method for (the well-posed) Oseen's problem used to improve robustness in the high Reynolds limit.Also with reference to high Reynolds computations for the well-posed case we here consider equal order interpolation for all fields.We wish to explore if the results on stability for the unique continuation for convection-diffusion equations in the limit of small diffusivity [19,20] carry over to the case of incompressible flow.The key observation there was that for smooth solutions to the convection-diffusion equation the method had Hölder stable error estimates when diffusion dominates, similar to the analysis above, but in the convection dominated regime the stability in a subdomain slightly smaller than that spanned by the characteristics intersecting the data zone is Lipschitz.In that zone the convergence for the ill-posed problem coincides with that of the well-posed problem for piecewise affine approximation.As a means to study the effect of incompressibility we compare with the case where in addition to  in   ,  is also provided as data in Ω.The proposed method can be modified to accommodate this case by including 1  2 ‖ ℎ − ‖ Ω , as an additional term in the Lagrangian (17).Note that when the pressure is added the velocity pressure coupling is strongly reduced.The relative  2 -errors for running the same problem as above are displayed in Figure 9. Left side plots of Figure 9 show the results without adding any additional pressure term, and right side plots of Figure 9 display the results by including the pressure data.We observe that the results with pressure information are consistently better than those without.In particular for high order polynomials and high Reynolds number the information on the pressure appears to provide a very strong enhancement of the stability.Further, the effect of data perturbations for different values of the viscosity coefficient is studied with and without the pressure augmentation, see Figures 10-11.We observe that if a priori information on the pressure is added and viscosity is reduced the convergence order for the relative  2 -error increases.This is consistent with the results of [19,20].If the pressure is not added however we do not observe this effect and it appears from these computational examples that we can not expect the result from [20] to hold for linearized incompressible flow.
In Figure 10-11 we present the results under perturbations of data.These results show that the robustness under perturbations is also substantially enhanced if the pressure is known, indicating that the pressure velocity coupling introduces a strong sensitivity to perturbations.

Conclusions
We have introduced a finite element data assimilation method for the linearized Navier-Stokes' equation.We proved the natural extension of the error estimates of [9] valid for piecewise affine approximation to the case of arbitrary polynomial orders.The expected increase in convergence rate was obtained, but the estimates also show that the sensitivity of the system to perturbations in data increase.The theoretical results were validated on some academic test cases.The main observations are that high order approximation for the illposed linearized Navier-Stokes' equations pays off, at least for sufficiently clean data.The spaces for the dual variables on the other hand can be chosen with piecewise affine approximation without loss of accuracy of the approximation.A study where the viscosity was varied showed that the incompressibility condition and the associated velocity-pressure coupling severely compromise the convective Lipschitz stability that is known to hold in the zone in the domain defined by points on the characteristics intersecting the data zone.If additional  data in the form of global pressure measurements were added the results improved and were similar to the those of the scalar convection-diffusion equation.Future work will focus on the nonlinear case and the possibility of enhancing stability by adding knowledge of some other variable than the pressure, such as for example a passive tracer as in scalar image velocimetry [18].

Figure 1 .
Figure 1.Sketch of the domains used for computations in Section 5.1.(a) Convex geometry.(b) Non-convex geometry.

Figure 2 .
Figure 2. Relative error for geometrical setup displayed in Figure 1.(a) Errors for convex geometry in Figure 1a.(b) Errors for non-convex geometry in Figure 1b.

Figure 5 .
Figure 5. Relative error with using the minimal polynomial order for geometrical setup displayed in Figure 1.(a) Errors for convex geometry in Figure 1a.(b) Errors for non-convex geometry in Figure 1b.

Figure 6 .
Figure 6.Relative error with using the minimal polynomial order for geometrical setup Figure 1a in terms of the strength of the data perturbation.(a)  = 0. (b)  = 1.(c)  = 2.

Figure 7 .
Figure 7. Relative error with using the minimal polynomial order for geometrical setup Figure 1b in terms of the strength of the data perturbation.(a)  = 0. (b)  = 1.(c)  = 2.

Figure 8 .
Figure 8. Data set   and error measurement regions (B).