title: On False Accuracy Verification of UMUSCL Scheme
authors: Hiroaki Nishikawa
text: 
  
  
  
   
  
  
   
  
  
  
  
  
   On False Accuracy Verification of UMUSCL Scheme
      
  Hiroaki NishikawaAssociate Research Fellow (hiro@nianet.org),
   100 Exploration Way, Hampton, VA 23666 USA.
  
    National Institute of Aerospace,
   Hampton, VA 23666, USA 
  
      December 30, 2023
  =============================================================================================================================================================================
  
  
  
  
   
  In this paper, we reveal a mechanism behind a false accuracy verification encountered with unstructured-grid schemes based on solution reconstruction 
  such as UMUSCL. Third- (or higher-) order of accuracy has been reported for the Euler equations in the literature, but UMUSCL is actually second-order 
  accurate at best for nonlinear equations. False high-order convergence occurs generally for a scheme that is high order for linear equations but second-order for nonlinear equations. It is caused by unexpected linearization of a target nonlinear equation due to too small of a perturbation added to an exact solution used for accuracy verification. To clarify the mechanism, we begin with a proof that the UMUSCL scheme is third-order accurate only for linear equations. Then, we derive a condition under which the third-order truncation error dominates the second-order error and demonstrate it numerically for Burgers' equation. Similar results are shown for the Euler equations, which disprove some accuracy verification results in the literature. To be genuinely third-order, UMUSCL must be implemented with flux reconstruction. 
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  
  This paper is a sequel to the two previous papers <cit.>, where we clarified the MUSCL and QUICK schemes towards the clarification of economical high-order unstructured-grid schemes for practical computational fluid dynamics (CFD) solvers, e.g., third-order UMUSCL with κ=1/2 <cit.>, κ=1/3 <cit.>, or κ=0 <cit.>. 
  In this paper, we will clarify one more confusion: the  false  accuracy verification of the UMUSCL scheme. 
  
  
  
  
  
  
  The UMUSCL scheme of Burg <cit.> is
   generally considered as an unstructured-grid extension 
  of Van Leer's κ-reconstruction scheme <cit.> and has been widely employed in practical CFD solvers <cit.> 
  with a confusion over the value of κ for giving third-order accuracy on regular or one-dimensional grids: κ=1/2 <cit.>, κ=1/3 <cit.>, or  
  κ=-1/6 <cit.>. 
  The confusion arises mainly from taking different combinations of numerical solution and discretization types  : numerical solutions stored as point-valued or cell-averaged
  solutions; discretizations of a differential form of a conservation law at a point (finite-difference) or of an integral form over a cell (finite-volume). 
  For example, Burg originally proposed the UMUSCL scheme as a finite-volume scheme with point-valued numerical solutions stored at nodes on an unstructured grid <cit.>.  Note that we know he used point-valued numerical solutions because that is the only way he could obtain third-order accuracy
  with κ=1/2 for a one-dimensional nonlinear system (see Refs.<cit.> for details). This is already confusing because the MUSCL scheme is based on cell-averaged numerical solutions, not point-valued solutions (see Ref.<cit.>). As clarified in the previous paper  <cit.>, a third-order finite-volume scheme with point-valued numerical solutions is nothing but the QUICK scheme and therefore the UMUSCL scheme should have been called the UQUICK scheme. In fact, third-order accuracy with κ=1/2, which is true for the QUICK scheme, has been confirmed for a one-dimensional steady conservation law by Burg <cit.> (see also Ref.<cit.>). Third-order accuracy demonstrated by Burg is genuine but only for one-dimensional problems; it cannot be third-order in multi-dimensions even for Cartesian grids unless the flux is integrated over a face by a high-order quadrature formula. 
  
   To be even more confusing, in many or perhaps all practical unstructured-grid codes, the UMUSCL scheme is implemented 
    not as a finite-volume scheme but as a point-wise scheme with the time derivative and source/forcing terms evaluated at a solution point (a node or a cell center)  <cit.>.
   
    In this paper, we will focus on this particular implementation. We will call it simply the UMUSCL scheme but it should not be confused with Burg's 
   finite-volume UMUSCL scheme. 
   Then, the fact that the scheme has been shown to achieve up to fourth-order accuracy with a single flux evaluation per face on Cartesian grids <cit.> indicates that the scheme is actually a finite-difference scheme, approximating the differential form of a target equation at a solution point.  Note that it does not matter how the discretization is derived; the resulting 
   discretization must be high-order as a finite-difference scheme, not as a finite-volume scheme because the time derivative and source/forcing terms are not integrated with high-order 
   quadrature over a cell. Therefore,  the UMUSCL scheme corresponds to neither the MUSCL scheme nor the QUICK scheme and does not achieve high-order accuracy in the same way as the MUSCL scheme does. As we will show, the UMUSCL scheme is third-order accurate for linear equations with κ=1/3, but only second-order accurate when applied to nonlinear conservation laws. 
   This feature is common to conservative finite-difference schemes with a flux evaluated with reconstructed solutions as in MUSCL (e.g., those in Refs.<cit.>). 
   Hence, high-order verification results reported in the literature for the UMUSCL scheme applied to the Euler equations  are misleading
  and/or misinterpreted. 
   
   
    To be genuinely third- or higher-order accurate, it is necessary to directly reconstruct the flux as pointed out for somewhat similar schemes in Ref.<cit.>. In other words, the UMUSCL scheme of Refs.<cit.> can be easily made third- or higher-order by direct flux reconstruction. 
   
   It is worth pointing out that a similar unstructured-grid MUSCL scheme had already been proposed in Refs.<cit.>. The scheme involves a parameter β (instead of κ) and is called the 
  β-scheme. Its basic form is equivalent to the UMUSCL scheme with κ=1-2β. In fact, the β-scheme is proposed as a finite-difference-like scheme 
  with point-valued solutions at nodes, not as a finite-volume scheme, and has been shown to achieve third-order accuracy on a regular 
  grid with β=1/3 (i.e., κ=1/3, not κ=1/2) for linear equations <cit.>. Therefore, all discussions in this paper will equally apply to 
  the β-scheme.  
  
  
  
  
  
   
  It should be noted also that the UMUSCL scheme has been shown to bring significant improvements to complex flow simulations and thus it is indeed useful   <cit.>. It is just that the improved resolution is largely due to reduced dissipation 
  by high-order solution reconstruction, not by high-order accuracy (see Refs.<cit.> for a relevant discussion). For scale-resolving turbulent-flow simulations, however, it would be strongly desired to achieve high-order accuracy since a high-order scheme is expected to be much more efficient than second-order methods on highly refined grid. 
  It is therefore important to  reveal the accuracy limitation of UMUSCL and develop a genuinely high-order version.
  
  
  
  
  
  
  
  To accomplish the task, we first show that the UMUSCL scheme,  in the form typically implemented in a practical code with the time derivative and source/forcing 
  terms evaluated at a solution point, is equivalent to the high-order conservative scheme of Shu and Osher <cit.> only for linear equations and thus cannot be high-order for nonlinear equations. Then, we reveal the mechanism behind the  false  third-order error convergence and derive a condition under which a third-order truncation error dominates the second-order one. Finally, we numerically demonstrate the  false  third-order accuracy of UMUSCL and genuine third-order accuracy of a flux-reconstruction version for the Burgers and Euler equations. Our focus is on third-order accuracy for smooth solutions, which is sufficient to illustrate the problem. A monotonicity property (which is important for discontinuous solutions) and higher-order accuracy are beyond the scope of the paper and will be discussed elsewhere.  Accuracy verification will be performed with exact solutions or by the method of 
  manufactured solutions <cit.>. See Ref.<cit.> for a deeper background on verification; but the issue discussed in this paper is not 
  addressed in Ref.<cit.>.
  
  
   
  The paper is organized as follows. In Section 2, we will explain why the UMUSCL scheme cannot be third-order accurate for nonlinear
  equations unless the flux is directly reconstructed. 
  In Section 3, we will describe the mechanism of how a target nonlinear equation is linearized by an exact solution. 
  In Section 4, we will present numerical results to confirm the  false  third-order error convergence of UMUSCL
  and genuine third-order accuracy with flux reconstruction.  
  In Section 5, we conclude the paper with remarks.
  
  
  
  
  
  § SECOND/THIRD-ORDER FOR LINEAR/NONLINEAR EQUATIONS
  
  
  
  
  
  
  
  
   §.§ UMUSCL is second-order accurate
  
   
  
  Consider a conservation law in two dimensions:
  
      ∂_t u +  ∂_x f  + ∂_y g = s(x,y),
  
  where u is a solution variable, f and g are fluxes, and s is a forcing function.
  Our target scheme is the UMUSCL scheme with point-valued numerical solutions stored 
  at nodes and point evaluations of the time derivative and forcing terms as widely used in practical unstructured-grid 
  solvers <cit.>. For a general unstructured
  grid, the UMUSCL scheme is given at a node j by 
  
      d u_j/dt   + 1/V_j∑_ k ∈{ k_j }ϕ_jk(u_L,u_R)  A_jk = s_j,
   
  where u_j is a point-valued solution at the node j, V_j is the measure of the  dual control volume around the node j, A_jk is the area (length in two dimensions) of the face between j and k,  { k_j } is a set of neighbor nodes of j, ϕ_jk is a numerical flux along a face normal 
  as a function of two states u_L and u_R reconstructed from the nodes j and k, respectively, by the so-called UMUSCL reconstruction scheme <cit.>, and s_j is a point evaluation of s(x,y) at j.  Note that we store numerical solutions as point values at nodes (there are no cell-averaged 
  solutions stored anywhere) and the time derivative ∂_t u is simply evaluated at the node as d u_j/dt. 
  
  As clarified in the previous paper <cit.> and will be discussed in detail later, the UMUSCL scheme corresponds to  a generalized version of the QUICKEST scheme of Leonard  <cit.> as well as the κ-family finite-difference scheme of Van Leer <cit.> or more generally the conservative finite-difference scheme of Shu and Osher <cit.> in one dimension or on Cartesian grids. 
  
  
  A major advantage of the UMUSCL scheme implemented in the form (<ref>) is its efficiency: it can achieve third-order accuracy on regular grids (at least for linear equations) whereas the original finite-volume version <cit.> requires high-order flux quadrature to achieve third-order accuracy even on regular grids and also requires a consistent treatment of the cell-averaged time derivative for unsteady problems <cit.>. A further discussion on efficiency will be given in a subsequent paper.
   
  In this paper, we focus on a Cartesian grid as shown in Figure <ref>, which is sufficient to illustrate the  false  accuracy problem, and denote the numerical solution at a node (i,j) by u_i,j. Then, the UMUSCL scheme reduces to
  
      d u_i,j/dt   +  F_i+1/2,j -   F_i-1/2,j/h +    G_i,j+1/2 -  G_i,j-1/2/h= s_i,j,
   
  where x_i+1,j - x_i,j = y_i,j+1 - y_i,j = h, s_i,j = s(x_i,j , y_i,j ), and F and G are numerical fluxes defined at a face,
  
      F(u_L,u_R) =  1/2[  f(u_L)    + f(u_R)     ]  - D/2  ( u_R - u_L),
  
  with the dissipation coefficient D= |∂ f/ ∂ u|, and u_L and u_R are solutions reconstructed with Van Leer's κ-reconstruction scheme <cit.> implemented in the form applicable to general unstructured grids 
  as suggested by Burg <cit.>: e.g., at a face (i+1/2,j), 
  
      u_L    =   κu_i +  u_i+1/2 +   (1-κ) [  u_j  + h/2(u_x)_i] ,     (u_x)_i  = u_i+1 -  u_i-1/2h, 
      
       [1ex]
      u_R    =   κu_i+1 +  u_i/2 +   (1-κ) [  u_i+1 - h/2(u_x)_i+1],    (u_x)_i+1  = u_i+2 -  u_i/2h,
  
  where the other subscript j has been omitted for brevity and κ is a parameter. 
  The other numerical flux G is computed in a similar manner in the y-direction. Refs.<cit.> add extra terms to the above reconstruction schemes in an attempt to achieve even higher-order accuracy. 
  As we will show, it is the flux evaluation with the reconstructed solution f(u_L), not the reconstruction scheme itself, that leads to accuracy deterioration for nonlinear equations. Therefore, such extra terms are irrelevant to the discussion here and thus not considered. 
  
   
    
  
  
  At first glance, the UMUSCL scheme (<ref>) looks like a finite-difference scheme approximating the differential form (<ref>) at the node (i,j). 
  As such, it is second-order accurate for any κ because it consists of the central difference approximations: 
  
      F_i+1/2,j -   F_i-1/2,j/h≈∂_x f   + O(h^2),     G_i,j+1/2 -  G_i,j-1/2/h≈∂_y g   + O(h^2),
  
  no matter how accurate the fluxes F and G are. 
   However, as mentioned in Introduction, third- or higher-order accuracy has been observed with κ=1/3 as
  reported in the literature  <cit.>. To clarify the apparent contradiction, we 
  will show that the scheme (<ref>) has a slightly different interpretation and does achieve third-order accuracy for linear conservation laws, at least. 
  
  
  
  
  
  
  
  
  
  
   
  
  
  
   §.§ UMUSCL can be third-order accurate for linear equations
  
   
  
  
  In 1989, Shu and Osher <cit.> proposed a high-order conservative finite-difference method based on the following 
  relationship (which is exact on uniform grids; see Ref.<cit.> for a discussion on its accuracy on non-uniform grids):
  
      ∂_x f   =   F̃_i+1/2,j -   F̃_i-1/2,j/h,
  
  where F̃ is a function whose cell-average in the x-direction is f(u_i,j), which can be easily derived as 
  
      ∂_x f   =  ∂_x (  1/h∫_x-h/2^x+h/2F̃ dx )  =  F̃_i+1/2,j -   F̃_i-1/2,j/h.
  
  Therefore, a high-order scheme can be constructed in the form: 
  
      d u_i,j/dt   + F̃_i+1/2,j -  F̃_i-1/2,j/h +  G̃_i,j+1/2 -  G̃_i,j-1/2/h= s_i,j,
   
  by accurately reconstructing F̃ and G̃. Then, since these fluxes are point values at faces and the data available at nodes are their cell-averages, 
  we have to reconstruct the point-valued flux, e.g., F̃, from its cell-averages f_i,j = f(u_i,j). 
  That is what Van Leer's κ-reconstruction scheme does exactly for a quadratic function with κ=1/3 <cit.>: e.g., again at a face (i+1/2,j),   
    
      f_L    =   κ f_i +  f_i+1/2 +   (1-κ) [  f_j  + h/2(f_x)_i] ,     (f_x)_i  = f_i+1 -  f_i-1/2h, 
      
       [1ex]
      f_R    =   κ f_i+1 +  f_i/2 +   (1-κ) [  f_i+1 - h/2(f_x)_i+1],    (f_x)_i+1  = f_i+2 -  f_i/2h,
  
  
  where we have denoted, for brevity, f_i-1  = f(u_i-1,j), f_i  = f(u_i,j), f_i+1  = f(u_i+1,j), f_i+2  = f(u_i+2,j).
   Note that the above expressions may not look like but are equivalent to the original κ-scheme (the terms 
  proportional to κ form a quadratic term; see Ref.<cit.>). 
  Then, with the following flux function, 
  
      F̃(f_L, f_R, u_L,u_R) =  1/2[  f_L    + f_R   ]  - D/2  ( u_R - u_L),
  
   where D= |∂ f/ ∂ u|, the conservative finite-difference scheme (<ref>) achieves third-order accuracy with κ=1/3, which has been verified analytically and numerically in the previous paper <cit.>.
  
  
  
  
  
  
  Notice now that the scheme (<ref>) is equivalent to the UMUSCL scheme, 
   which evaluates the left and right fluxes with reconstructed solutions,  f(u_L) and f(u_R),
   for a linear equation, i.e.,  f = a u with a constant a, 
  because we then have 
  
      f(u_L) = f_L ,      f(u_R) = f_R ,
  
  and thus
  
      F̃(f_L, f_R, u_L,u_R) =  F(u_L,u_R).
  
  That is, the solution reconstruction is equivalent to the flux reconstruction for a linear flux. 
  Therefore, the UMUSCL scheme achieves third-order accuracy with κ=1/3 for a linear equation on Cartesian grids. 
  
  
  
  
  
   §.§ UMUSCL cannot be third-order accurate for nonlinear equations
  
   
  
  For nonlinear fluxes,  the solution reconstruction is no longer equivalent to the flux reconstruction:
  
      f(u_L)    f_L  ,     f(u_R)   f_R,
  
  and thus
  
      F̃(f_L, f_R, u_L,u_R)   F(u_L,u_R).
  
  Therefore, the UMUSCL scheme is not equivalent to the conservative finite-difference scheme and simply reduces to 
  a second-order finite-difference scheme as in Equation (<ref>). 
  It implies also that any point-wise numerical scheme with F(u_L,u_R), e.g., those in Refs.<cit.>, is second-order accurate at best. 
  
  To see this more clearly, consider a quadratic function whose cell-average is the flux at a node j: 
  
      f̃(x)  =   f_i + (f_x)_i  (x - x_i ) + 1/2 (f_xx)_i[  (x - x_i )^2 -   h^2/12]
        ,    (f_xx)_i  = f_i+1 - 2 f_i  + f_i-1/h^2,
  
  
  
  which gives at a face i+1/2 
  
      f̃_i+1/2  =   f_i + h/2 (f_x)_i  + h^2/12 (f_xx)_i,
  
  which can be expanded, for Burgers' equation f=u^2/2, as 
  
      f̃_i+1/2  =  u_i^2/2  + 1/2 ( u_i ∂_x u ) h + 1/12[     (∂_x u)^2  + u_i ∂_xx u   ] h^2
        +O(h^3).
  
  This at least needs to be matched exactly by a numerical flux in order to generate a third-order scheme. For the direct flux reconstruction with κ=1/3, we find
  
      f_L    =   κf_i +  f_i+1/2 +   (1-κ) [  f_j  + 1/2(f_x)_i  h ]   
       [2ex]
         =   u_i^2/2 +   1/2  ( u_i ∂_x u ) h +  1/12[     (∂_x u)^2  + u_i ∂_xx u   ] h^2+ O(h^3),
    
  which matches the exact expansion (<ref>). 
   However, for the UMUSCL scheme, we find with κ=1/3, 
  
      f(u_L)  
           =   1/2{κu_i +  u_i+1/2 +   (1-κ) [  u_j  + h/2(u_x)_i] 
          }^2  
       [2ex]
         =   u_i^2/2 +   1/2  ( u_i ∂_x u )h  + 
       1/12[  3/2 (∂_x u)^2  + u_i ∂_xx u   ]   h^2 + O(h^3),
    
  which does not match the exact expansion (<ref>) in the quadratic term. Therefore, the UMUSCL scheme is only second-order accurate for nonlinear equations. To understand it more deeply, one has to understand that the flux at a face needs to be computed as a function whose cell-average is f_i and directly reconstructed from the cell-averages f_i, f_i+1, etc., whereas in the latter the flux is evaluated with the face value of a function whose cell-average is u_i, reconstructed from the cell averages u_i, u_i+1, etc.; it is well known that they differ by a second-order error <cit.> (the error in the above can be predicted by the theorem in Ref.<cit.>). See also Ref.<cit.> for discussions on the QUICKEST scheme, which is equivalent to the UMUSCL scheme. A further discussion will be given in a subsequent paper. 
  
  Nonetheless, it should also be clear by now that the UMUSCL scheme can be made third-order for nonlinear equations by performing the flux reconstruction to replace f(u_L) by f_L. 
   In the rest of the paper, the flux-reconstruction version of the UMUSCL scheme with the 
  numerical flux (<ref>) will be referred to as FSR (flux and solution reconstruction). The word `UMUSCL' is not used here because it is not MUSCL as mentioned earlier. We will come back to this scheme later in numerical experiments. Note that the necessity of flux reconstruction was discussed in Ref.<cit.> for a similar reconstruction-based scheme; but they did not rely on the exact relation (<ref>) and instead directly constructed a high-order approximation to the flux derivative. 
  
   It is confusing that third-order accuracy of the UMUSCL scheme has been observed for the nonlinear Euler equations 
   and such accuracy verification results have been used as a confirmation of third- and higher-order accuracy of the 
  UMUSCL scheme <cit.>. 
  Note that in Ref.<cit.>, the author observed third-order accuracy with a similar scheme but later realized that it was  false  for the reason discussed in this paper. Note also that Refs.<cit.> mention high-order flux quadrature but it 
   is not sufficient
  
   because they treat source terms as point values (thus not finite-volume). These results are not genuine and thus misleading (no other accuracy verifications are shown in Refs.<cit.>). 
  
  
  Here, the problem exists in the accuracy verification process as we will discuss in the next section.
  
  
  
  
  
  
  § UNEXPECTED LINEARIZATION BY EXACT SOLUTION
  
  
  
  
  
  
  
  
  
   §.§ Burgers equation
  
   
  
  The problem is illustrated for the steady Burgers equation, a representative of nonlinear conservation laws, in one dimension:
  
      ∂_x f = s(x),
  
  where f = u^2/2 and s(x) is defined such that a chosen function u_e(x) is made an exact solution: 
  s(x) = ∂_x  ( u_e^2 / 2). Consider a typical exact solution used in accuracy verification tests:
  
      u_e(x)   = u_∞ +  ϵ u_p(x),
  
  where u_∞ and ϵ are positive constants, and u_p(x) is an arbitrary smooth function (e.g., u_p(x)=sin(x)). 
  This is the type of solution that has been used to demonstrate high-order accuracy of the UMUSCL scheme <cit.>. 
  Consider the UMUSCL scheme for the steady Burgers equation:
  
      F_i+1/2 -   F_i-1/2/h = s_i,
  
  where the fluxes are computed as described in Section <ref>. 
  From the discussion in the previous section, we know this scheme cannot be third-order accurate because the Burgers equation is nonlinear. 
  However, it is possible that third-order accuracy is observed accidentally when it is tested with the exact solution of the form (<ref>) as
  we will explain in the next section.
  
  
  
  
   §.§ Linearized by exact solution
  
   
  
  
  
  
  
  False third-order error convergence is observed typically when the parameter ϵ is small, which effectively 
  linearizes the target nonlinear equation and allows a third-order error to dominate. 
  This is the mechanism behind the  false  third-order accuracy verification. 
  To see this, substitute the exact solution (<ref>) into the Burgers equation to get
  
      ϵ∂_x f' = s(x),
  
  where 
  
      f' =u_∞ u_p +   ϵ  u_p^2,
  
  which shows that the leading term of the flux is linear. Therefore, 
  the Burgers equation will behave like a linear equation as ϵ→ 0; and the UMUSCL scheme exhibits third-order accuracy. 
  A question arises, then, about how small ϵ should be in order to observe third-order error convergence, which is the subject for the next section.
  
  
  
  
   §.§ Dominated by third-order error
  
   
  
  To estimate how small ϵ should be for the  false  third-order error convergence to occur, we substitute a smooth exact solution 
  into the residual,
  
      Res_i =    F_i+1/2 -   F_i-1/2/h -  s_i,
  
  and expand it to obtain the truncation error TE (see Ref.<cit.>): 
  
      TE = 
       1/24[     f_uuu (u_x)^3 +  6 κ  f_uu u_x  u_xx   +  2 (  3 κ - 1   ) f_u  u_xxx] h^2    
      + 
      1/12[  u_i u_xxxx  +   u_x u_xxx] h^3  + 
       O(h^4).
  
  where the second- and third-order terms come from the averaged flux term and the dissipation term in the numerical flux (<ref>), respectively. 
  For simplicity, we have assumed u_i > 0 for all i. 
  Observe that the second-order error vanishes with κ=1/3 for linear equations, for which f_uuu= f_uu = 0.
  However, the Burgers equation gives f_uu=1 and thus the second-order term remains:   
  
      TE =   1/12 ( u_x )( u_xx  ) h^2 + 
      1/12[  u_i u_xxxx  +   u_x u_xxx] h^3  + O(h^4).
   
  If the smooth exact solution is given by Equation (<ref>) with u_p(x) =  sin( ω x), where ω is a constant, then
  
      TE =   ϵω^3  [  T_2     +     T_3   ] h^2 + O(h^4),
   
  where
  
      T_2 = - ϵsin( 2 ω x ) /24
      ,    T_3 =  -ω{ϵ ( 1 - 2  sin^2( ω x)  )  - u_∞sin(ω x)    }  h  /12.
   
  Our interest is to see how small ϵ should be in order for T_3 to dominate T_2, and thus the scheme 
  will become third-order accurate. It suffices then to compare the upper bounds (|  T_3  | reaches the maximum at x=π/2 ω): 
  
      |  T_2  |_max =   ϵ/24   , 
          |  T_3  |_max =   ω  ( u_∞ + ϵ  ) h  /12,
   
  from which we find 
  
      |  T_3  |_max >   |  T_2  |_max    ϵ/u_∞< 2 h ω/1- 2 h ω .
   
  Therefore, the scheme will behave as if it is a third-order scheme when ϵ satisfies the above condition.
  For example, if an accuracy verification study is performed with ω=2 π over a series of grids with the finest grid of 128 nodes in a unit 
  domain giving h =1/127, we find
  
      ϵ/u_∞ < 0.1098136158 …,
   
  which one would easily choose, since the parameter ϵ is meant to be a perturbation, if not aware of the problem. 
  It is important to note that the analysis shows that the order of accuracy changes based on the value of ϵ. This is quite different from a typical observation that 
  a lower-order scheme is more accurate on coarse grids, where the order of accuracy does not change. 
   Note also that the problem does not occur even on coarse grids if ϵ is sufficiently large. 
  One can estimate a critical value of h by solving |  T_3  |_max >   |  T_2  |_max for h:
  
      h > ϵ/u_∞/ 2 ω ( ϵ /u_∞ + 1 ).
  
  For example, for ω= 2 π, we find that third-order accuracy will be observed for grids coarser than the one with 
  h = 0.007234315595 ≈ 1/127 for ϵ/u_∞ = 0.1 and 
  h = 0.02652582384 ≈ 1/38 for ϵ/u_∞ = 0.5. 
  In other words, one must use grids finer than h=1/127 and h=1/38, in the former and latter cases, respectively, in order to avoid false third-order error convergence.
  
  
  
  
  
   §.§ Euler equations
  
   
  
  The same mechanism applies to the Euler equations:
  
      ∂_x  f  =  s(x) ,    f = ( ρ, ρ u^2 + p, ρ u H ),
  
  where  ρ is the density, u is the velocity, p is the pressure, and ρ H = γ p/ (γ-1) + ρ u^2/2 with γ=1.4.
  A similar unexpected linearization occurs when an exact solution for the velocity is defined as in Equation (<ref>) 
  while exact solutions for the density and the pressure can be defined by arbitrary smooth functions.
  
  
  
  
  
  Substitute the velocity (<ref>) into the Euler equations, 
  and obtain
  
      ∂_t  u + ∂_x  f_1 + ϵ∂_x  f_2  =  s(x) ,
  
  where
  
      f_1
       =
       [  [                             ρ  u_∞; [1ex]
                     ρ u_∞^2   + p;          [1ex]
                    u_∞ H ]]
      ,    f_2 =  
       [  [                                                                              ρ   u_p;                                    [1ex]
                   2   ρ u_p    u_∞  +  ϵρ  u_p^2; [1ex]
              u_p   [  H +     ρ  u_∞^2   ] +   ϵρ/2 (    3 u_∞ + ϵ  u_p    )  u_p^2 ]],
                       H = γ p /γ-1  + 1/2ρ u_∞^2  .
  
  The factor ϵ in the second flux derivative term indicates that the system will be linear as ϵ→ 0.
  Therefore, we expect to observe third-order error convergence for the Euler equations if the solution (<ref>)
  is used as an exact solution with a small ϵ. For an estimate of the largest ϵ, we do not attempt to derive it  
  for the Euler equations and instead demonstrate that the estimate (<ref>) serves rather well also for the Euler equations.
  
  Before proceeding, we remark that there are two other cases where the Euler equations are effectively linearized. 
  These case are somewhat trivial but have been used for accuracy verification in the literature, e.g., in Ref.<cit.>.
  One is a case where the velocity and the pressure are constant (i.e., the so-called entropy wave). In this case, only the density varies in space and the Euler equations 
  become linear since the flux f is a linear function of ρ. The other is a case where the velocity is constant. 
  In this case, only the density and the pressure vary in space and the Euler equations again become linear since the flux f is a linear function of ρ and p. 
  Therefore, third-order accuracy will be observed with the UMUSCL scheme. In Ref.<cit.>, high-order accuracy of a WENO-type scheme 
  based on solution reconstruction is verified for an entropy-wave solution in one dimension. Their scheme seems different from but should be similar to the UMUSCL scheme because the third-order part is equivalent to the UMUSCL scheme with κ=1/3 and also it is shown to produce high-order solutions on regular grids in two dimensions without high-order flux 
  quadrature over each face. Then, it must be a finite-difference scheme and thus flux reconstruction is required for achieve third- and higher-order accuracy. 
  
  
  
  
  
  § NUMERICAL RESULTS
  
  
  
   
  In this section, we present numerical results to demonstrate that third-order accuracy is indeed observed for the Burgers equation when the condition (<ref>) is satisfied,
  that the same is observed for the Euler equations in one dimension, 
  and finally that the same is observed in two-dimensional problems. The last part is intended for disproving high-order accuracy reported in the literature.  
  In all numerical experiments, the numerical solutions are point values stored at nodes (which can be interpreted as cell centers also). 
    
  
  
  
   §.§ Burgers equation in one dimension
  
  
  
  We consider the Burgers equation (<ref>) in a unit domain x ∈ [0,1] and the exact solution (<ref>) with 
  
      u_∞ = 0.3 , 
         
      u_p(x) = sin (2  π x),
  
  which corresponds to taking ω=2 π, and solve it over a series of grids with n = 16, 32, 64, 128, and  256 nodes. 
  To exclude effects of boundaries (which are beyond the scope of this paper), 
  we fix the exact solution at two boundary nodes and their neighbors also, and solve the system of nonlinear residual equations
  for those stored at the rest of the cells. The dissipation coefficient D is evaluated with the averaged solution at a face: (u_L + u_R )/2.
  An implicit solver based on the exact Jacobian of the first-order scheme is used to solve the residual equations,
  which is taken to be converged when the L_1 norm of the residual is reduced by nine orders of magnitude from an initial residual norm computed
  with the initial solution u_i = 1, n=3,4, ⋯, n-3, n-2.
  See Ref.<cit.>, for example, for details of the implicit solver for a one-dimensional finite-volume-type scheme. 
  
   
  To investigate the effect of the parameter ϵ, we consider varying ϵ as follows,
  
      ϵ  = c_ϵ u_∞,    c_ϵ = 0.05,     0.1,     0.3,     0.5,
  
  which correspond, respectively, to the critical mesh spacing (<ref>), denoted by h_cr, 
   
      h_cr =     0.000410  ,       0.00201,     0.00747,       0.0107.
   
  Comparing these with the actual grid spacings, 
   
      h  =  0.0666,   0.0322 ,   0.0159,   0.00787 ,      0.00392,
  
  we expect that third-order accuracy will be observed for c_ϵ = 0.05 and 0.1, and second-order accuracy will begin to 
  dominate somewhere between the last two grids for c_ϵ = 0.3 and between the third and fourth grids for c_ϵ = 0.5.  
  
  
  
  
  
  
   Also, we consider three values of κ, κ=0, 1/2, and 1/3, to remind that third-order accuracy 
  is possible only with κ=1/3.
  
  Error convergence results are shown in Figure <ref>, where the discretization error (i.e., solution error) is defined by
  
      L_∞(  E_p) = max_i  ∈{ I } | u_i - u_i^exact |,
        
      u_i^exact = u_e(x_i) = u_∞ + ϵsin (2 π x_i ).
  
  where { I } denotes the set of nodes in a given grid. 
   
  As expected, third-order error convergence is observed for κ=1/3 with c_ϵ = 0.05 and 0.1, as shown in the first two plots.
  Figure <ref> shows the result for c_ϵ = 0.3. Second-order convergence is observed in the last two grids as expected.
  Finally, Figure <ref> shows that c_ϵ = 0.5 leads to third-order error convergence begins to deteriorate somewhere
  between the third and fourth grids, again as expected.
  
  
  
  
  
  
  
  
    
  
  
  
  
  
  
   §.§ Euler equations in one dimension
  
   
  
  Next, we consider the Euler equations (<ref>) with the exact solution defined in terms of the primitive variables:
  
      w^exact (x)
      =
      [
      [ ρ^exact (x); u^exact (x); p^exact (x) ]]
      =
      [
      [   1 +  ϵ_ρsin( 2.3 π );     u_∞ +  ϵsin( 2 π ); 1 + ϵ_p sin( 2.5 π x ) ]],
  
  and the forcing term defined by s(x) = ∂_x  f( w^exact(x)). The amplitudes for the density and the pressure are set ϵ_ρ=ϵ_p=0.2 unless otherwise stated. For the velocity solution, we set u_∞ = 0.3 and define ϵ by Equation (<ref>) to see how well the formula derived for the Burgers equation works for the Euler equations.
   
  
  In the UMUSCL scheme, we perform the solution reconstruction in the primitive variables w=(ρ,u,p) instead of the conservative
  variables u=(ρ, ρ u, ρ H - p). Then, the residual is given by
  
      Res_i =   F_i+1/2 -  F_i-1/2/h -   s_i,
  
  where the numerical flux is given by the Roe flux <cit.>,
  
      F( w_L, w_R) =  1/2[    f( w_L)   +   f( w_R)   ]  - D̂/2[   u( w_R) -   u( w_L) ],
  
  D̂ is the Roe dissipation matrix evaluated with the Roe averages <cit.>, and the two states w_L and w_R are computed by
  the reconstruction formulas (<ref>) and (<ref>).  As before, we consider three values of κ: κ=0, 1/2, and 1/3.
  
  
  The resulting system of steady residual equations is solved by a pseudo-time integration with the three-stage SSP Runge-Kutta scheme <cit.> 
  with a local time step at CFL=0.99. The solver is taken to be converged when the residual is reduced by seven orders of magnitude in the L_1 norm, starting from
  the initial norm computed with the initial constant solution: w_i = (1, u_∞ , 1 ), i=3,4, ⋯, n-3, n-2. Again, the exact solutions are specified at i=1,2 , n-1, and n, to exclude boundary effects.
  For simplicity, the discretization error is defined as the maximum of the L_∞ norms among the three variables:
  
      L_∞(  E_p) = max_k=1,2,3   max_i  ∈{ I }   |  w_i(k) -  w^exact_i(k) |.
  
  
  Results are shown in Figure <ref>.  As expected, third-order accuracy is observed for c_ϵ=0.05 and 0.1, when κ=1/3 as shown in Figures <ref> and <ref>. For c_ϵ=0.3 and 0.5, the second-order error begins to dominate fairly accurately as predicted by the estimates h_cr, as can be seen in Figures <ref> and <ref>. 
  From the fact that the actual value of ϵ is 0.15 for the last case (c_ϵ=0.5), we would expect that 
  false third-order accuracy can be entirely avoided if ϵ is set equal to the corresponding coefficients for the other variables, i.e., ϵ = 0.2. 
   
  
  
  
  
  
  
  
  
  
  
     
     
     
  
  
  
  
  
    
  
  
  To demonstrate that the UMUSCL scheme can be made genuinely third-order with the flux reconstruction (i.e., FSR), we consider the following numerical flux,
  
      F( w_L, w_R,  f_L ,  f_R) =  1/2[    f_L   +   f_R  ]  - D̂/2[   u( w_R) -   u( w_L) ],
  
  where the left and right fluxes f_L and f_R are directly reconstructed as in Equations (<ref>) and (<ref>). 
  This is slightly more expensive because w_L and w_R still need to be computed as well, for example. However, the scheme is genuinely third-order 
  and thus third-order accuracy is obtained for any value of c_ϵ as can be seen in Figure <ref>. 
  
  
  
  
  
  
    
  
  
  For the sake of completeness, we present results for the following three cases: 
  (a)  ϵ_ρ=0.2 and ϵ=ϵ_p=0 (an entropy wave solution <cit.>),  
  (b)  ϵ_ρ=ϵ_p=0.2 and ϵ=0,
  (c)  ϵ_ρ=ϵ_p=ϵ=0.2. 
  As mentioned in the previous section, the Euler equations are linearized in the cases (a) and (b), where the velocity is constant. Therefore,
  the UMUSCL scheme may give third-order accuracy. However, the case (c) renders the Euler equations fully nonlinear and thus 
  is expected to reveal second-order accuracy of the UMUSCL scheme. For all the cases, UMUSCL is compared with FSR and κ=1/3 is used for both schemes.
  Results are shown in Figure <ref>. As expected, UMUSCL gives third-order accuracy for 
   the cases (a) and (b), but deteriorates to second-order in the case (c). 
  On the other hand, FSR gives third-order accuracy for all the cases. 
  
  
    
  
  
  
  
   
  
  
  
  
   §.§ Euler equations in two dimensions
  
   
  
  Finally, we consider two-dimensional cases for the Euler equations, which have been used to mistakenly verify high-order accuracy of the UMUSCL scheme 
  in Refs.<cit.> (and also of a relevant scheme in Ref.<cit.>). 
  We will demonstrate that third-order accuracy is due to the unexpected linearization and thus not genuine but that 
  genuine third-order accuracy can be achieved with the flux reconstruction.
  
  Consider the Euler equations in two dimensions, 
  
      ∂ u/∂ t +   F 
        =  s ,
           u =  [  [                        ρ; [1ex]
                     ρ v; [1ex]
                     ρ E ]],    F =  
       [  [                               ρ v^t; [1ex]
                     ρ v v  + p   I;        [1ex]
                     ρ v^t H ]] ,
  
  where  v is the velocity vector with Cartesian components v = (u,v)^t (which is a column vector; the superscript indicates the transpose),  denotes the dyadic product, I is the identity matrix, 
  and E = H-p / ρ  = (p/ρ) / (γ-1)  +  v^2 /2 is the specific total energy. 
  
  
   
    
  
  
  The UMUSCL scheme is implemented in the form applicable to a general unstructured grid <cit.>, i.e., at a node j as in 
  Figure <ref>, 
  
      d  u_j/dt   + 1/V_j∑_ k ∈{ k_j }Φ_jk(u_L,u_R)  A_jk=  s_j,
  
  with the Roe flux,
  
      F( w_L, w_R) =  1/2[    f( w_L)   +   f( w_R)   ]  - D̂_n/2[   u( w_R) -   u( w_L) ],
  
  where f =  F·n̂_jk is the flux projected along the face normal n_jk = n̂_jk | n_jk| = n̂_jk  A_jk, D̂_n = | ∂ f / ∂ u  | is the dissipation term evaluated with the Roe averages, the solution reconstruction is performed with the primitive variables 
  w = (ρ,u,v,p):
  
      w_L    =   κ w_j +   w_k/2 +   (1-κ) [   w_j  + 1/2∇ w_j · (   x_k -  x_j  )  ] ,  
      
       [1ex]
       w_R    =   κ w_k +  w_j/2 +   (1-κ) [   w_k  +  1/2∇ w_k · (   x_j -  x_k  )   ],
  
  where x_j and x_k denote the nodal coordinates of j and its neighbor k, respectively and the gradients are computed by a linear least-squares
  method. On Cartesian grids, the scheme (<ref>) reduces to the form (<ref>) with the κ-reconstruction scheme applied as a one-dimensional algorithm in each coordinate direction. Note that no high-order flux quadrature is used; the flux is computed only at the midpoint of each face. 
  
  For the flux reconstruction version, FSR, the numerical flux is computed as
  
      F( w_L, w_R,  f_L ,  f_R) =  1/2[    f_L  +   f_R   ]  - D̂_n/2[   u( w_R) -   u( w_L) ],
  
  where
  
      f_L    =   κ f( w_j) +  f( w_k)   /2 +   (1-κ) [  f( w_i)  + 1/2(  ∂ f/∂ w)_ j∇ w_j · (   x_k -  x_j  )  ] ,  
      
       [1ex]
       f_R    =   κ f( w_k) + f( w_j) /2 +   (1-κ) [   f( w_k) +  1/2(  ∂ f/∂ w)_ k∇ w_k · (   x_j -  x_k  )   ],
  
  
  
  
  
  
  
  
  
   
  
  and the flux Jacobian is given, with the unit face normal vector n̂_jk and v defined as column vectors and the notation u_n =  v· n_jk, as
  
      ∂ f/∂ w 
       =
        [  [                       u_n                  ρn̂_jk^t                         0;              [1ex]
      u_n  v    ρ ( u_n  I +  vn̂_jk )                     n̂_jk;         [1ex]
      u_n  v^2 /2 ρ(H  n̂_jk^t + u_n  v^t )              γ u_n /(γ-1) ]].
  
  Note that the flux gradients are efficiently computed from the solution gradients by the chain rule. Therefore, it is not necessary to compute the flux gradients directly by a least-squares method and thus does not also require extra storage for flux gradients. This efficient version of FSR is referred to as FSR-CR (chain rule) in the rest of the paper. 
  Ref.<cit.> seems to perform the flux reconstruction and therefore their method may be genuinely high-order accurate. However, their method requires the computation and storage of the the flux gradients; it can be very expensive in three dimensions. A further discussion including fourth- and fifth-order accurate schemes will be given in a subsequent paper.
  
  
  
   
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    §.§.§ Steady problem
  
   
  
  Let us begin with a steady case in a unit square, similar to the one used in Ref.<cit.>, with the exact solutions defined by
  
      ρ   =     1  + 0.2  sin[   π (  2.3  x + 2.3 y ) ], 
      
      u                       =     u_∞  + ϵsin[    π ( 2  x + 2 y )  ], 
       
      v                       =     v_∞  + ϵsin[    π ( 2   x + 2 y ) ], 
       
      p                      =     1  + 0.2 sin[     π ( 2.5  x + 2.5 y ) ],
  
  where u_∞  = 0.15, v_∞ = 0.02, and ϵ is given the following values:
  
      ϵ   = 0,    0.05,    0.1,    0.2.
  
  The forcing term s_j is numerically computed with the above exact solutions at a node j and added to the residual (see Ref.<cit.> for details). 
  To exclude boundary effects, we specify the exact solution at nodes on the boundary and at their neighbors. 
  The residual equations are solved using an 
  implicit defect-correction solver (see, e.g., Ref.<cit.>) until the residual norm is reduced
   by seven orders of magnitude for a series of regular quadrilateral grids: 49×49, 65×65, 81×81,
   97×97, 113×113, 129×129. Figure <ref> shows exact pressure contours on a
   coarser 17×17 grid for the sake of illustration. For this problem, we compare UMUSCL and FSR-CR with κ=1/3.
  
  
  
  
    
  
  
  
  
  Figure <ref> shows the L_∞ error convergence results for the pressure obtained with the UMUSCL scheme.
  As expected, the convergence rate deteriorates as ϵ increases towards 0.2. 
  Nearly third-order convergence is observed 
  for ϵ=0.05, but it eventually deteriorates to second-order at ϵ=0.2. 
   The exact solution used in Ref.<cit.> is similar to the above solution with ϵ=0.05; thus third-order accuracy reported in 
   Ref.<cit.> is not genuine (the author did not realize it at the time of submitting the paper). It is important to observe that the convergence rates are all close to third-order on coarse grids in Figure <ref>: accuracy verification must be conducted with sufficiently fine grids to observe the true order of accuracy. 
    On the other hand, if we perform the flux reconstruction, then third-order accuracy is achieved for all values of ϵ as 
   shown in Figure <ref>. 
  
  
  
  
  
  
  
    §.§.§ Unsteady inviscid vortex transport
  
   
  
  Finally, we consider an inviscid-vortex problem frequently used for accuracy verification
  <cit.>, which
  is an exact solution to the Euler equations without forcing terms. Following 
  Refs.<cit.>, we set the exact solutions as 
  
      u =  u_∞  -   K y/2 πexp(  1-r ^2 /2) ,    
      v  = v_∞   +  K x/2 πexp(  1-r ^2 /2),
  
  
  
  
  
  
  
  
  and
  
      T =  1  -   K^2 (γ-1)  /8 π^2 exp(  1-r ^2   ) ,   ρ = T^1/γ-1 ,   
      p = ρ ^γ/γ    ,
  
  where x = x  - u_∞ t, y = y  - v_∞ t, r ^2 =x^2 + y^2, and (u_∞,v_∞)=(0.2,0.0).  
  The initial solution at t=0 is shown in Figure <ref>. Here, the parameter K corresponds to ϵ in the previous cases; 
  
  we would expect that the unexpected linearization 
  occurs when K is small. Two cases are considered: (1)K=1 as in Refs.<cit.> (the one in Ref.<cit.> is much smaller) and (2)K=5 as in Ref.<cit.>. 
  As we will see, K=1 is small enough for the unexpected linearization to occur, and K=5 is large enough to reveal second-order accuracy of the UMUSCL scheme.
  
  
  
  
  For our purpose, it suffices to perform the calculation for a short time. Thus, 
  we compute the solution at the final time t_f=1.0 with the three-stage  
  SSP Runge-Kutta scheme <cit.> for the total of 1000 time steps with a constant time step Δ t = 0.001,
  which is so small that errors are dominated by the spatial discretization. To verify the spatial order of accuracy, we perform the computation
  over a series of n×n regular quadrilateral grids, where n=48, 64, 80, 96,112, 128, 144, 160, 176, 192, 208, 224, 240, 256. 
  The coarsest grid is shown with pressure contours in Figure <ref>.
   
  Error convergence results for K=1 are shown in Figure <ref>. 
  Here, for consistency with the results reported in Refs.<cit.>, we compute the L_2 error norm for the pressure. 
  As can be seen, the UMUSCL scheme exhibits nearly third-order accuracy with a slight deterioration on fine grids. Also plotted is error convergence
  obtained with the FSR-CR scheme. As can be clearly seen, it is genuinely third-order with no sign of deterioration. 
  
  Figure <ref> shows the results obtained for K=5. 
  It is clearly seen now that  the UMUSCL scheme quickly deteriorates to second-order. It is not genuinely third-order as expected. 
  On the other hand, the FSR-CR scheme maintains third-order accuracy.
  
  These results indicate that high-order error convergence reported in Refs.<cit.> 
  is not genuine, but also that the UMUSCL scheme of Refs.<cit.> can be made genuinely 
  high-order by flux reconstruction, which can be performed efficiently with the chain rule as suggested. 
  
  
  
  
  
    
  
  
  
  
  
   
  
  
  
  § CONCLUSIONS
  
  
  
    
  We have shown that the UMUSCL scheme of Burg <cit.>, which is 
  defined with point-valued numerical solutions, is second-order accurate at best for nonlinear equations 
  if used with point-valued time derivatives and source/forcing terms as in Refs.<cit.>. 
  As shown, the UMUSCL scheme is equivalent to the high-order conservative finite-difference scheme of Shu and Osher <cit.> for linear equations, but not for nonlinear equations. Thus, it cannot be high-order for nonlinear equations unless the flux is directly reconstructed. 
  Nevertheless, as we have shown, it can exhibit third-order error convergence for a manufactured solution with a small perturbation
  because a target nonlinear equation is effectively linearized by such a solution. An estimate of a critical value of the perturbation parameter has been
  derived for the Burgers equation and its validity has been confirmed by numerical experiments. Similar results have been obtained for the Euler equations
  in one dimension, which suggests the perturbation parameter for the velocity should be of the same order as those for other variables
  to avoid the  false  high-order error convergence. Finally, we have demonstrated  false  third-order accuracy of UMUSCL
  and genuine third-order accuracy of FSR (a flux and solution reconstruction scheme) for a steady problem and an unsteady inviscid vortex transport problem widely used for accuracy 
  verification. In conclusion, any point-wise numerical scheme based on fluxes evaluated with reconstructed solutions cannot be high-order for nonlinear equations 
  and accuracy verification must be carefully performed with exact solutions of a significant variation in order to avoid  false  high-order error convergence.
  
  Thus far, we have clarified third-order accuracy of the MUSCL scheme <cit.>, third-order accuracy of the 
  QUICK scheme <cit.>, and the  false  third-order accuracy of the UMUSCL scheme; we are now ready to clarify economical high-order 
  unstructured-grid schemes and identify one of the most efficient third- and higher-order schemes for practical simulations, which will be discussed in detail in a subsequent paper. 
  
  
  
  
  tocsectionAcknowledgments
  
  
  § ACKNOWLEDGMENTS
  
  
   
  The author gratefully acknowledges support from Software CRADLE, part of Hexagon, the U.S. Army Research Office 
  under the contract/grant number W911NF-19-1-0429 with Dr. Matthew Munson as the program manager, and 
   the Hypersonic Technology Project, through the Hypersonic Airbreathing Propulsion Branch of the NASA Langley
   Research Center, under Contract No. 80LARC17C0004. 
   
  
  
  tocsectionReferences
  
  aiaa
  
  
  
  
  
  
  
  
  
  
  
   
  