A comprehensive high-order solver verification methodology for free fluid flows
Farshad Navah Siva Nadarajah
physics.comp-ph physics.flu-dyn





















]Farshad Navahmycorrespondingauthor
farshad.navah@mail.mcgill.ca

]Siva Nadarajah
Mechanical Engineering Department, 
McGill University,
Montréal, Québec, Canada, H3A 0C3


[mycorrespondingauthor]Corresponding author













The aim of this article is to present a comprehensive methodology for the verification of computational fluid dynamics (CFD) solvers with a special attention to aspects pertinent to discretizations with orders of accuracy (OOAs) higher than two. The method of manufactured solutions (MMS) is adopted and a series of manufactured solutions (MSs) is introduced that examines various components of CFD solvers for free flows (not bounded by walls), including inviscid, laminar and  turbulent problems when the latter are modelled by the Reynolds-averaged Navier-Stokes (RANS) equations. The treatment of curved elements is also examined. These MSs are furthermore conceived with demonstrated suitability for the verification of OOAs up to the sixth. Each MS is as well utilized to discuss salient aspects useful to the code verification methodology such as the relative qualities of the most useful norms in measuring the discretization error, the sensitivity analysis of the verification process to forcing function terms, the relation between residual minimization and discretization error convergence in iterative solutions and finally the sensitivity of high-order discretizations to grid stretching and self-similarity. Furthermore, scripts and code are provided as accompanying material to assist the interested reader in reproducing the verification results of each manufactured solution (MS).



Code verificationHigh-order accuracy Method of manufactured solutions RANS Curved grids Flux reconstruction











§ INTRODUCTION


Code verification is a crucial step prior to the application of a scientific simulation software to the solution of practical problems as it aims at examining the soundness of the implementation of the governing equations in the numerical framework. With the increasing interest of the research community in the design and application of high-order-of-accuracy discretization methods for CFD problems, there is an imperative need to extend the verification methodology to this class of methods. Code verification is in fact even more critical for higher-order methods since it is the only means to provide assurance that the effort invested in their design and development is justified by the delivery of the expected higher performance in terms of accuracy per computational effort. We hence  present in this paper the fundamental aspects towards a comprehensive code verification methodology for CFD solvers with all orders of accuracy. 

To carry out the demonstration of the methodology and without loss of generality, we choose the numerical framework composed of the compressible RANS equations closed by the original and modified versions of the Spalart-Allmaras (SA) model of turbulence and discretized by the correction procedure via flux reconstruction (CPR) scheme.

The article is structured as this: in Section <ref>, the context and contributions of this work are introduced, followed by a comprehensive presentation of the theoretical background in verification and validation (V&V) in Section <ref>. The governing equations as well as the compact high-order numerical method are respectively exhibited in Sections <ref> and <ref>, including a precise description of all the employed boundary conditions. The application results of the verification and the discussion of the salient aspects of the methodology appear in Section  <ref> and the article ends with conclusions in the last Section.




 §.§ Contributions

A series of trigonometric manufactured solutions for the sequential verification of high-order RANS solvers is devised such that it demonstrably achieves all OOAs up to the sixth order on moderately fine isotropic grids, without being trivially reached on the coarsest ones. Attention is invested in ensuring that the MSs produce a fair balance between different terms of the governing equations. The sequence of MSs targets constitutive components of solvers in an isolated fashion and with incremental complexity such that systematic debugging is enabled and gathering cumulative evidence on the soundness of high-order CFD solver implementation is made possible. The MSs serve thus to examine the implementation of Euler, Navier-Stokes (NS) and RANS equations along with the original and also with the modified SA model, for free flows, i.e., flows that are not bounded by walls. The set of MSs is as well employed to explore the following concepts:


  * The comparative description of different norms and a demonstration of the importance of L_∞ norm in code verification;

  * The need for the inclusion of a relatively high order of accuracy in code verification;

  * The significance of the balancing of forcing function terms of the MMS and the sensitivity analysis of the verification process to terms with the lowest magnitude in the forcing functions;

  * The verification of both the original and the modified SA models of turbulence;

  * The relation of residual convergence level with regards to discretization error magnitude and insight on the necessary level of residual convergence;

  * The examination of the treatment of non-affine mapping of curved elements;

  * The effect of grid self-similarity and stretching on grid convergence of solutions with smooth gradients.


Accompanying (see <cit.>) IPython <cit.> notebook and C routine facilitate the application of the described verification methodology through the reproduction of the manufactured fields and forcing functions of the presented MSs.



§ THEORETICAL BACKGROUND


In this section, first the terminology involved in V&V is completed and defined further, the MMS is formalized then and finally, a short review of the previous works with a focus on verification via the MMS in CFD is presented.



 §.§ Terminology in V&V

Figure <ref> illustrates the relation between major concepts of interest under the three themes of simulation process, error sources, as well as verification and validation. 

As a scientific simulation process takes place, errors from various sources slip into its different steps, contaminating incrementally the outcome of the process. Abstractly, as the reality that we aim to capture cascades through a simulation, it diminishes at each step of the process. The role of V&V is hence to ensure that the amount of original reality captured by the simulation is sufficient for the purpose that the simulation is meant to serve, by ideally providing a dependable measurement of the discrepancies. In what follows, we describe more precisely these ideas with reference to Fig. <ref>.




Any scientific simulation process starts from a reality, a physical phenomenon that it aims at reproducing. Based on experimental data and previous theoretical knowledge, the relation between various quantities playing a role in the physical process is described by a conceptual model, i.e., a series of mathematical equations such as the partial differentials equations (PDEs) of Euler, NS or those of RANS-based turbulence models. Almost in all cases, these models are a mere, yet hopefully reliable, approximation of the inherent complexity of physics and as such they contain a modelling error. In this regard, model validation questions how well the phenomenon of interest is approached by the conceptual model. For complex problems such as those encountered in CFD, this is however only considerable once the simulation process has ended and a numerical solution is available. In order to solve the conceptual model, a numerical algorithm, a scheme, is often needed. The application of this numerical recipe to the conceptual model yields a numerical model such as the discretization of the RANS-SA system of PDEs by the CPR scheme. As the size of the discrete problem is increased for well-posed and smooth solutions, a numerical model is expected to tend towards the conceptual model with a rate known as the formal order of accuracy. A numerical model with suitable properties such as stability and efficiency is translated to a computer code. Considering the complexity of the numerical model, programming errors often occur at this step.  Indeed, according to an exhaustive analysis of the quality of scientific computing codes: "There were about 8 serious static faults per 1000 lines of executable lines in C, and about 12 serious faults per 1000 lines in Fortran"<cit.>.  Code verification has for purpose to identify and eliminate the mistakes affecting the correspondence between the scientific software and the conceptual model via the formal order of accuracy. Similarly to model validation, code verification relies on the numerical solution of specific problems. To solve a given problem, the spatial and temporal domains are discretized by a set of points, called degrees of freedom (DOFs), to which the discrete solution is associated. The discretization error is the difference between the continuous and discrete solutions and the solution verification is the estimation of this error for a given solution. The solution process refers to the minimization of discrete residual equations by iterative methods, mandatory for tackling non-linear systems, and by algorithms handling linear algebraic systems. The truncation of real values for representation on computer architectures, by double precision types for example, introduces a round-off error that affects the numerical solution by propagating through the discrete equations. On the other hand, a lack of sufficient minimization of the discrete residual equations results in an iterative convergence error that imposes a gap between the achieved numerical solution and the actual solution of the discrete problem. Both the round-off and iterative convergence errors need to be controlled during the solution process to ensure that these sources of error are minimized such that the discretization error is isolated as the major source of numerical error. This condition enables code and solution verifications to be carried out since they operate only on the discretization error and its rate of reduction for increasing DOFs.



  §.§.§ Code verification methods

The evaluation of OOAs in code verification requires the knowledge of the exact solution which could be provided by devising a problem for which such a solution exists. This approach is called the method of analytic solutions (MAS) which has the advantage of avoiding code modification and hence being non-intrusive. However, MAS is not always applicable as in the case of RANS equations for which analytic solutions of interest do not exist. Whenever applicable, MAS often involves simplified solutions that might not properly exercise all the governing equation terms. An example of an analytic solution for the NS equations is the Couette problem <cit.>. 

An alternative and more powerful approach is via the MMS which enables the verification of any arbitrary set of governing equations but which requires point-wise code modifications and hence is considered to be an intrusive method. Other advantages of the MMS are:


  * The possibility of debugging by systematic simplification and elimination of the governing equation terms;

  * The examination of arbitrary boundary conditions;

  * The verification of specific scheme properties by devising specific MSs.
 

For these reasons, we adopt the MMS as the verification method of choice.



 §.§ Verification via the MMS

Consider a linear differential operator, 𝒜, producing a homogeneous output, uniquely when applied to Q^ex, a time-invariant and smooth scalar field on a bounded spatial domain, viz.,

    𝒜 (Q^ex)= 0.

For a tessellation of the domain augmented by a polynomial expansion of degree P resulting in a representative resolution size, h, the equivalent discrete equation reads

    𝒜^h(Q^h) = ϵ≈0,

where 𝒜^h and Q^h  are respectively the discrete operator and solution and
ϵ is a finite precision representation of zero such as the machine precision. The truncation error is defined as

    Θ_𝒜^h = 𝒜^h (Q^ex) - 𝒜 (Q^ex) = 𝒜^h (Q^ex).



According to the Lax equivalence theorem, if 𝒜^h is stable, i.e., its variations with regards to finite variations of Q^h are bounded, and if 𝒜^h is furthermore consistent, i.e., 
    lim_h→ 0Θ_𝒜^h = 𝒜 (Q^ex)=0,
 sufficient conditions are fulfilled for the discrete solution to be convergent, i.e., 

    lim_h→ 0 Q^h = Q^ex,

and hence the discretization error, ℰ_Q= (Q^h - Q^ex ),   tends to zero with ℰ∼𝒪(h^o) where · is a measurement in an appropriate norm and o is the formal rate of convergence (order of accuracy).

Now, let's apply the operator 𝒜 to an arbitrary smooth scalar field, Q^MS≠ Q^ex. Obviously,

    𝒜 (Q^MS) ≠ 0.

But, Eq. (<ref>) can be modified by the addition of a suitable forcing function, ℱ(Q^MS), such that

    ℛ(Q^MS)  = 𝒜 (Q^MS)  - ℱ(Q^MS)= 0.

where ℛ designates the modified operator. An immediate choice of forcing function is ℱ(Q^MS)= 𝒜 (Q^MS). The modified counterparts of the discrete equation (<ref>) and the truncation error  (<ref>) respectively are 

    ℛ^h(Q^h) = 𝒜^h(Q^h) - ℱ(Q^MS)  = ϵ≈0,

and

    Θ_ℛ^h = ℛ^h (Q^MS) - ℛ (Q^MS) = ℛ^h (Q^MS),

via Eq. (<ref>).

The stability and consistency of ℛ^h follow from those of 𝒜^h and hence the convergence to the manufactured solution is established as

    lim_h→ 0 Q^h = Q^MS,

and ℰ_Q = Q^h - Q^MS∼𝒪(h^o) is expected. The reader is referred to the literature for further information on the value of the formal OOAs in finite difference and finite volume frameworks. As for continuous finite element and discontinuous compact variational methods, we refer to the comprehensive analyses in <cit.>. A generic result for such methods applied to linear (and by extension to non-linear) advection-diffusion problems is that the expected order in L_1, L_2 and L_∞ norms is o=P+1, i.e., ℰ∼𝒪(h^P+1), whereas it rather is o=P in H_1 norm and semi-norm. The definition of each norm is provided in  <ref>.

Let's note that the focus of this work is on the verification of the steady-state RANS equations since a semi-discrete scheme is considered. In other words, the integration in time is segregated from the spatial discretization and could hence be verified in an independent fashion.



 §.§ Review of the literature

The V&V in CFD originates from pioneering and comprehensive contributions of Roache <cit.>, Oberkampf et al. <cit.>, Salari and Knupp <cit.> as well as Roy et al. <cit.>. The applications of the MMS via trigonometric manufactured solutions are reported for Euler and NS solvers <cit.>, for RANS solvers <cit.>,  for fluid-structure interaction <cit.> and more recently for multiphase flows <cit.>. Physically realistic manufactured solutions for the verification of RANS-based turbulence models, that mimic a turbulent boundary layer in the vicinity of no-slip wall, are introduced by Eça et al. <cit.> as well as by Oliver and colleagues <cit.>. Finally, a library for verification via MMS is presented in <cit.>. The mentioned applications involve low-order finite volume and finite element codes. Applications of the MMS to the verification of high-order spatial discretizations are found in <cit.> for DG/FR/CPR methods and in <cit.> for mixed finite differences methods. However, these applications are of limited scope as each targets a few specific solver features at a time. 





§ GOVERNING EQUATIONS





The governing equations are the steady-state, compressible RANS equations assuming a Newtonian and calorically perfect gas along with the conservative form of the original and revised SA models of turbulence <cit.>. The SA model is chosen due to its simplicity and proven effectiveness for the simulation of aerodynamic flows. The revised SA model supplements the original model by conditional variations of some of its terms in order to provide enhanced numerical stability despite the occurrence of negative values of the model's working variable, especially in presence of coarse spatial discretizations of the boundary layer in high-order solutions.

We cast the governing PDEs into the following compact form for general unsteady advection-diffusion models:

    ∂_t(Q_k) +∂_j(F^inv_kj) - ∂_j(F^vis_kj) = S_k,

where Q_k is a state variable and the solution of the k^th partial differential equation (PDE) with k ∈ [1 .. N_eq] where  N_eq=N_d+3 is the number of equations based on  N_d  space dimensions; F^inv_kj and F^vis_kj respectively express the k^th inviscid (advective) and viscous (diffusive) fluxes for j∈ [1 .. N_d],  S_k  refers to the source term for equation k and we use the Einstein's summation convention for repeated indices. The expanded form of the conservation laws is retrieved by substituting the following expressions and values in Eq. (<ref>):
















  - Conservation of mass (k=1)

    Q_k=ρ,         F^inv_kj=ρ u_j,     F^vis_kj=0,     S_k=0;


  - Conservation of momentum (k ∈ [2 .. N_d+1] and i=k-1)

    Q_k=ρ u_i,     F^inv_kj=ρ u_ju_i +p δ_ij,     F^vis_kj=τ_ij,     S_k=0;


  - Conservation of energy (k=N_d+2)

    Q_k=ρ E,       F^inv_kj=ρ u_j H,       F^vis_kj=u_i τ_ij +ω_j,       S_k=0;


  - Transport of the turbulent working variable (k=N_d+3)

    Q_k=ρν̃,      F^inv_kj=ρ u_j ν̃,     F^vis_kj=1/σ(μ+ρν̃f_n) ∂_j ν̃,     
    S_k=  ρ 𝒫 - ρ 𝒟 + ρ 𝒯
    +c_b2/σρ ∂_j ν̃ ∂_j ν̃ 
     - 1/σ(ν+ν̃f_n) ∂_j (ρ ∂_j ν̃) .











The quantities appearing in equations (<ref>) through (<ref>) are defined as follows: ρ is the density, u = e_i u_i is the velocity vector with e_i being the i^th orthonormal basis vector of the Euclidean spatial system, E is the total energy per mass defined as E=e+1/2 (u_iu_i) where e is the internal energy which for a calorically perfect gas is defined as e=R/γ-1T where R is the gas constant and T is the temperature. The total enthalpy is defined as H=E + p/ρ with p denoting the pressure related to the energy via the ideal gas law as

    p=(γ-1)ρ(E-1/2(u_iu_i)),

where γ is the specific heat ratio (γ=1.4 for air).

In Eq. (<ref>), τ_ij are the components of the viscous stress tensor, τ, which for compressible Newtonian fluids read

    τ_ij=2  μ_eff  S_ij,
      with  
    S_ij = 1/2(∂_i u_j + ∂_j u_i)-1/3∂_k u_k δ_ij,

where μ_eff is the effective viscosity, defined as the sum of the dynamic viscosity, μ, and the eddy viscosity, μ_t, viz., μ_eff= μ+μ_t; and δ_ij represents the Kronecker delta. Note that the dynamic viscosity is assumed to be spatially constant throughout this work.

In Eq. (<ref>), ω_j =λ_eff ∂_j T is the j^th component of the heat flux vector where λ_eff is the effective thermal conductivity defined as λ_eff = λ + λ_t, with λ=γ R/(γ-1)μ/Pr, the molecular conductivity,  and λ_t=γ R/(γ-1)μ_t/Pr_t, the eddy conductivity. The laminar and turbulent Prandtl numbers are respectively set to Pr=0.7 and Pr_t=0.9 unless specified.







In the SA model, μ_t, the turbulent (eddy) viscosity, is expressed by



[left=μ_t= ρν_t =]alignat=2
   ρν̃ f_v1        ν̃≥0, 

   0                              ν̃< 0,    


where

    f_v1=χ^3/χ^3+c_v1^3,
      χ=ν̃/ν,
      
    c_v1=7.1,

and ν̃ is the working variable of the SA model that represents a turbulent kinematic viscosity.

In Eq. (<ref>), the production term, 𝒫, is defined as


[left=𝒫=]alignat=2
    c_b1(1-f_t2)s̃ ν̃         ν̃≥0, 

    c_b1(1-c_t3)       s ν̃         ν̃< 0,    


where c_b1 = 0.1355, f_t2=c_t3 exp(-c_t4 χ^2) with c_t3=1.2 and c_t4=0.5 is the laminar suppression term, s=|ε_ijk∂_j u_k | is the vorticity magnitude with ε_ijk standing for the Levi-Civita symbol of permutation, and s̃ is the modified vorticity defined as


[left=s̃=]alignat=2
    s+s̅               s̅ ≥-c_v2 s, 

    s+s(c_v2^2s+c_v3s̅)/(c_v3-2c_v2)s-s̅          s̅ < -c_v2s, 


where

    s̅=ν̃ f_ν 2/κ^2 d_w^2,   f_v2 = 1 - χ/1+χ f_v1,
       c_v2=0.7,    c_v3=0.9,   κ = 0.41,

and d_w being the distance to the closest wall. In Eq. (<ref>), the destruction term, 𝒟, is defined as


[left=𝒟=]alignat=2
    (c_w1 f_w -c_b1/κ^2f_t2) ν̃^2 /d_w^2         ν̃≥0, 

    - c_w1  ν̃^2 /d_w^2            ν̃< 0, 


where

    c_w1 = c_b1/κ^2+1+c_b2/σ,1/1/1^2/1^1/2
    c_b2 = 0.622,   σ =2/3,   and   f_w = g(1+c_w3^6/g^6+c_w3^6)^1/6.


We refer the reader to <cit.> for the full definition of the trip term, 𝒯, in Eq. (<ref>) which serves to mimic the effect of a forced transition. A value of 𝒯=0 is considered throughout this work. The remaining closure functions and constants of the SA model are


[left=f_n=]alignat=2
    1                                           ν̃≥0, 

    c_n1+χ^3/c_n1-χ^3          ν̃< 0, 













    g = r + c_w2(r^6-r),
      
    r = min(ν̃/s̃κ^2 d_w^2, r_lim),
      
    r_lim=10,
      
    c_n1=16,
      
    c_w2 = 0.3,
      and  
    c_w3 = 2.


The original SA model is represented by equations <ref>, <ref>, <ref>, <ref>, <ref> and <ref> whereas the modified SA model corresponds to equations <ref>, <ref>, <ref>, <ref>, <ref> and <ref>.








§ COMPACT HIGH-ORDER NUMERICAL METHOD


The high-order discretization scheme adopted in this work is the CPR <cit.>, that is an extension of the flux reconstruction (FR) <cit.> scheme to simplices by a lifting collocation penalty technique. Variations of the correction function at the core of this scheme allow to recover a number of prominent compact high-order methods such as the spectral difference, the spectral volume and the discontinuous Galerkin (DG). Hence, the FR/CPR method is regarded as a unifying compact high-order scheme and a representative member of this class.



 §.§ CPR Scheme



To discuss the formulation of the CPR scheme, let's consider a scalar hyperbolic conservation law,

    ∂_t Q + ∂_j F_j=0,

where Q is the state variable and F=e_q F_q(Q)∈ IR^N_d is a generic flux vector which could be or not function of state derivatives.
Multiplying by an arbitrary test function, ϕ, integrating and applying the divergence theorem, the variational formulation of Eq. (<ref>) in Green's form is obtained

    ∫_Ω∂_t(Q)  ϕ  dΩ - ∫_Ω F_j ∂_j(ϕ)  dΩ + ∫_Γ F_j  n_j ϕ  dΓ=0,

where Ω is a bounded spatial domain in IR^N_d, Γ is its frontier and n= e_q n_q is the local unit outward normal. Considering a tessellation, Ξ, of Ω, into a set of continuous and non-overlapping elements, e_i∈Ξ, the variational formulation reads

    ∑_e_i(∫_Ω_e_i∂_t(Q)   ϕ  dΩ 
    - ∫_Ω_e_i F_j ∂_j(ϕ)  dΩ  + ∫_Γ_e_i F_j n_j ϕ  dΓ)=0.

Taking into account the inner, Q^-, and the outer, Q^+, states defined with regards to the direction of n at a given element boundary, the inter-element coupling of discontinuous solutions can be achieved via a common Riemann flux, F̂^n(Q^-,Q^+,n^-,n^+),  such that

    F_jn_j ≡ F^n ≈F̂^n(Q^-,Q^+,n),

since for a well-constructed mesh, n≡n^-=-n^+ where - and + exponents respectively refer to the unit outward normal from the element e_i and from its neighbour at the same point on the interface.

Integrating by parts and applying the divergence theorem once again yields the variational formulation in the divergence form,

    ∑_e_i(∫_Ω_e_i∂_t(Q)   ϕ dΩ + ∫_Ω_e_i∂_j(F_j) ϕ dΩ  + ∫_Γ_e_i (F̂^n-F^n)  ϕ  dΓ)=0.

The boundary term can be projected into the element. This is achieved by the following lifting operator that provides a correction field, ψ, as output:

    ∫_Ω_e_i ψ ϕ  dΩ=∫_Γ_e_i (F̂^n-F^n) ϕ  dΓ.

We are interested in a correction that belongs to the space of polynomials of degree N_P or less, viz., ψ∈ IP^N_P since
the solution and the correction on element e_i are discretized by a set of interpolation polynomials φ_l with the property that φ_l ∈ IP^N_P on element e_i and identically 0 on the others, such that

    Q_e_i =  Q_l   φ_l,

and

    ψ_e_i = ψ_l  φ_l,

where Q_l and ψ_l are respectively the state and correction values at the solution point l ∈ [1 .. N_node] where N_node is the number of solution points per element (N_node=(N_P+1)^N_d for a tensor-product element). Throughout this article, the notation P is used equivalently to N_P. The short notation, P3, P4, etc. will also signify P=3, P=4, etc. Remark that integral operators in Eq. (<ref>) could be computed via quadratures on a reference element once and stored for use on physical elements, assuming that a valid geometrical mapping exists.

In the case where the flux is a non-linear function of the state, the flux divergence will not necessarily belong to the space IP^N_P. We therefore consider the following projection of the flux divergence onto the space IP^N_P:

    ∫_Ω_e_i∏_e_i (∂_j F_j)  ϕ  d Ω = ∫_Ω_e_i∂_j F_j ϕ  d Ω,

where ∏_e_i (∂_j F_j) is the projected flux divergence.
The resulting discrete formulation is

    ∑_e_i∫_Ω_e_i(∂_t(Q_e_i) + ∏_e_i (∂_j F_j) +  ψ_e_i) ϕ  d Ω=0,

where all terms are of degree N_P or less. Choosing a proper test space that guarantees solution uniqueness yields

    ∑_e_i(∂_t(Q_e_i) + ∏_e_i (∂_j F_j) +  ψ_e_i)=0,

which is the variational formulation in the divergence form of Eq. (<ref>) and which can be cast into a purely differential scheme, thus avoiding costly explicit quadratures. In fact, the direct projection of the divergence term via Eq. (<ref>) still requires the use of quadratures. Two more efficient alternatives could be considered <cit.>: the Lagrange polynomial (LP) and the chain rule (CR) approaches. 

The LP method consists of interpolating the flux before the application of the divergence operator, as

    (F_j)_e_i = F_j_l  φ_l,

where F_j_l is the j^th spatial component of the flux evaluated at the solution point l. The projection of the divergence can then be computed as

    ∏_e_i (∂_j F_j) = ∂_j (F_j)_e_i =  ∂_j(φ_l)  F_j_l.


The CR approach employs the flux Jacobian, A_j = ∂F_j/∂ Q, and the spatial derivatives of the state variable evaluated via

    ∂_j Q_e_i=∂_j(φ_lQ_l) = ∂_j(φ_l) Q_l,

for a projection onto the space IP^N_P that reads

    ∏_e_i (∂_j F_j) = A_j ∂_j Q_e_i.


Numerical experiments have shown <cit.> that while the LP is fully conservative, the CR is more accurate at the expense of a slight loss in strict conservation. A fix to retrieve the full conservation of the CR has been proposed in <cit.>. In this work, the CR and LP methods are used for the evaluation of the divergence of inviscid and viscous fluxes respectively. Furthermore, we choose the Lagrange polynomials as both test functions (ϕ in Eq. (<ref>)) and interpolation functions (φ in Eqs. (<ref>) and  (<ref>)), thus retrieving a pre-integrated version of the discontinuous Galerkin method. Finally, we utilize the Gauss-Legendre-Lobatto (GLL) set, also known as Gauss-Lobatto, as solution nodes on quadrangular elements.




 §.§ Numerical fluxes

For finite number of DOFs, the scheme described in the precedent section yields discontinuous solutions at element interfaces. The common interface flux of Eq. (<ref>) couples these solutions by taking the following form for the considered governing equations:

    F̂^n=ℋ^inv(Q^-,Q^+,n) + ℋ^vis(Q^-,Q^+,n,e_q(∂_q Q)^-,e_q(∂_q Q)^+),

where ℋ^inv and ℋ^vis are proper numerical fluxes for respectively the inviscid and viscous components of the equations and ·^- and ·^+ refer to the traces of a given quantity from the inner and outer sides of the interface respectively. The choices of numerical fluxes are further described in what follows.


  §.§.§ Inviscid flux

We opt for the Roe's approximate Riemann solver <cit.> as ℋ^inv in Eq. (<ref>), expressed as

    ℋ^inv = 1/2(F^inv_j(Q^-) n_j+F^inv_j(Q^+) n_j +D),

where D is the dissipation added to the central flux, stabilizing it through upwinding.

As we solve the RANS and the SA equations in a fully coupled fashion, the Roe numerical flux is re-derived following Appendix D of <cit.> to account for this coupling. The resulting upwinding contrasts with that of decoupled approaches in which only the normal velocity at the surface is considered in the dissipation term of the SA common flux. In the coupled upwinding, the dissipation term accounts as well for sound speed and pressure. The coupling is reported <cit.> to enhance the smoothness and robustness of the RANS-SA solutions.



  §.§.§ Viscous flux


An intuitive approach for deriving a viscous numerical flux would be to consider a central flux based on the element-wise approximation of  ∂_q Q, computed as ∂_q(φ_l)Q_l, such as

    ℋ^vis  = 
    1/2( F^vis_j(Q^-,e_q(∂_q Q)^-) +  F^vis_j(Q^+,e_q(∂_q Q)^+)   )  n_j.


However, this choice suffers from instabilities for under-resolved elliptic problems  due to the singularity of the resulting discrete viscous operator <cit.>. Different approaches have been devised to alleviate this problem that mostly rely on the basic idea of introducing a correction on the solution derivatives that takes into account the data from neighbouring elements and penalizes the derivatives for solution jumps at the interface. We consider a representative method of this class that is the second flux of Bassi and Rebay <cit.>, referred to as BR2, which employs an auxiliary vector function, θ = e_qθ_q,  such that

    θ_q≈∂_q Q.


Applying the approach described in section <ref> to Eq. (<ref>) yields its variational formulation in the divergence form,

    ∫_Ω_e_iθ_q ϕ  dΩ - ∫_Ω_e_i∂_q (Q) ϕ  dΩ=∫_Γ_e_i (Q̂-Q)  ϕ  n_q  dΓ,

where Q̂ =1/2(Q^-+Q^+) is the common interface value of the state variable. 

By considering a lifting operation such that

    ∫_Ω_e_i R_q ϕ  dΩ =∫_Γ_e_i (Q̂-Q)  ϕ  n_q  dΓ,

where R_q is a correction field accounting for the data from the neighbouring element, the term ∂_q Q≡θ_q=∂_q Q+R_q can be interpreted as a corrected divergence that serves in the computation of the viscous flux, F^vis = e_jF^vis_j(Q,e_q∂_q Q) for the discontinuous interface flux value in Eq. (<ref>) and for flux divergence value in Eq. (<ref>). Note that lim_h→0 R_q = 0 is expected and θ_q = ∂_q Q is recovered in the continuum limit. As for the viscous numerical flux in Eq. (<ref>), such a correction for the divergence results in an extended stencil since it couples the elemental solution to the data from neighbours of its immediate neighbours. Nevertheless,  considering a segmentation of Γ_e_i into N_f non-overlapping and continuous faces, the viscous numerical flux on face (Γ_e_i)_f, with f ∈ [1  ..  N_f], can be evaluated as 


    ℋ^vis =
    1/2( F^vis_j(Q^-,e_q(∂_q Q)^-) +  F^vis_j(Q^+,e_q(∂_q Q)^+)   )  n_j,

where ∂_q Q=∂_q Q+C_r_q r_q designates a partially corrected derivative. The constant C_r_q ensures the stability of the scheme for highly diffusive problems and is set to the value of 2 in this work. Finally, the partial correction, r_q, is the result of the following lifting operation:

    ∫_Ω_e_i r_q ϕ  dΩ =∫_(Γ_e_i)_f (Q̂-Q)  ϕ  n_q  dΓ.

To conserve the efficiency of tensor product operators, Eq. (<ref>) is applied in a 1D manner in the context of FR/CPR schemes on quadrangular and hexahedral elements <cit.>.




 §.§ Boundary conditions


The soundness of a CFD implementation can not be established without the verification of its treatment of boundary conditions  (BCs). Therefore, in this section, we precisely define the expression of the BCs that will be examined in this work.



  §.§.§ Riemann BC

The inviscid boundary conditions for free flows are imposed in the sense of weak-Riemann of <cit.> via the numerical flux of Eq. (<ref>) that takes into account the incoming and outgoing characteristics arising from the inner, Q_k^- and the outer, Q_k^+=Q_k^BC, states where Q_k^BC is the k^th boundary state defined at ghost nodes. 

In the case of a manufactured solution we set

    Q_k^BC = Q_k^MS|_Γ_BC,

where Q_k^MS|_Γ_BC designates the value of the manufactured solution evaluated at the point of the boundary where the BC is imposed.

It is important to note that for verification purposes, the treatment of the boundary conditions should be the same in the manufactured case as in the general flow problem. Hence, for the general flow problem, one would simply provide the desired boundary state at the ghost node instead of Q_k^MS|_Γ_BC. For a farfield BC for example, the free-stream values, Q_k^∞, are imposed by setting Q_k^BC=Q_k^∞.



  §.§.§ Viscous BC

The boundary condition for the viscous terms is enforced through the numerical flux of Eq. (<ref>), by setting the desired values of the state variables and their derivatives at the ghost nodes via Q_k^+=Q_k^BC and (∂_q Q_k)^+=(∂_q Q_k)^BC, and by taking the common state in Eqs. (<ref>) and (<ref>) to be Q̂ = 1/2(Q^-+Q^BC). 

In the case of manufactured solutions, the boundary value of the state is set via Eq. (<ref>) and its derivatives read

    (∂_q Q_k)^BC = (∂_qQ_k)^MS|_Γ_BC.

In the case of a farfield boundary employed in an application problem, the boundary values are set to Q_k^BC=Q_k^- and (∂_q Q_k)^BC=(∂_qQ_k)^-, such that ℋ^vis =
 F^vis_j(Q^-,e_q(∂_q Q)^-)   n_j at the farfield boundary. 




























 §.§ Treatment of curved elements

We refer the reader to Eq. (37) of Section 2.3 of <cit.> for a detailed description of the approach adopted in this work to treat curved elements. Furthermore, all mappings from computational element to physical element are isoparametric, i.e., they are based on polynomial expansions of same degree as the solution.





 §.§ Solution process

Starting from an initialization by exact values, an implicit Euler scheme, i.e., the relaxed Newton's method with pseudo-time integration, is employed to smoothly converge the non-linear residuals to the vicinity of the final solution by gradually decreasing the relaxation (increasing the pseudo-time step) down to the final stage where the full Newton's method yields a quadratic convergence to the discrete solution. The linearization of the residual equation with regards to the solution (the Jacobian matrix) is required by the Newton's method to determine the direction of descent. In order to ensure optimal residual minimization (to machine precision), an analytic Jacobian is implemented for the fully coupled system of RANS-SA including the original and modified portions of the SA equation. This linearization is then verified to provide significant digits of double precision when compared to a linearization by complex step via operator overloading <cit.>. The operator overloading automatically ensures that the Jacobian matrix is consistent with the residual equation such that any modification in the latter is inherently accounted for. The effort involved in the achievement of an exact Jacobian  is justified since it enables the full operation of the Newton's method in the minimization of residuals. In fact, as iterative and round-off errors scale with residuals <cit.>, their minimization permits to safely direct the focus towards the discretization error. The linear system is solved by a generalized minimal residual (GMRES) method along with a global block-Jacobi and local ILU(0) preconditionners via the PETSc package <cit.>.







§ VERIFICATION CASES AND METHODOLOGY






The five manufactured solutions considered in this work, i.e., MS-1 through MS-5 (see Table <ref>) are defined by a combination of trigonometric functions with terms depending exclusively on x or on y (mono terms) along with cross terms depending on both independent variables simultaneously. The cross terms are necessary for the verification of the viscous terms of the governing equations which feature mixed second order partial derivatives  whereas the mono terms are included to maintain the spatial dependency in the case where the cross terms are turned off for debugging purposes by setting their coefficients to zero. The generic form of the MS is

    ρ^MS   = ρ_0   + ρ_x sin(a_ρ_xπ   x / L)  + ρ_y cos(a_ρ_yπ   y / L)  + ρ_xycos(a_ρ_x yπ   x / L)  cos(a_ρ_x yπ   y / L),   
    
    u^MS   = u_0     + u_x   sin(a_u_xπ   x / L)  + u_y    cos(a_u_yπ   y / L)   + u_xycos(a_u_xyπ   x / L)  cos(a_u_xyπ   y / L),   
    
    v^MS   = v_0     + v_x   cos(a_v_xπ   x / L)  + v_y  sin(a_v_yπ   y / L)   + v_xycos(a_v_xyπ   x / L) cos(a_v_xyπ   y / L),   
    
    p^MS   = p_0     + p_x  cos(a_p_xπ   x / L)  + p_y       sin(a_p_yπ   y / L)   + p_xycos(a_p_xyπ   x / L)  cos(a_p_xyπ   y / L),   
    ν̃^MS   = ν̃_0     + ν̃_x   cos(a_ν̃_xπ   x / L)  + ν̃_y       cos(a_ν̃_yπ   y / L)   + ν̃_xycos(a_ν̃_xyπ   x / L)  cos(a_ν̃_xyπ   y / L),

where u and v respectively refer to the first and second velocity components, L=1.0 is a reference length and the manufactured total energy, E^MS, is obtained via Eq. (<ref>) with γ=1.4. 


We choose to define the MS in terms of the primitive variables since this procures a few advantages: firstly, since the Mach (Ma) number features the primitive variables directly, tuning the primitive variables of the MS allows for an easier control over the flow regime (subsonic vs. supersonic, compressible vs. incompressible). Secondly, the derivatives appearing in the forcing functions are those of the primitive variables and the MS definition in terms of the latter spares futile chain-rule operations. Finally, the expression of MS in primitive variables is forthwith compatible with the non-conservative formulation of incompressible solvers.

The constants appearing in Eqs. (<ref>) serve to tune the magnitudes and frequencies of the trigonometric functions of the manufactured primitive variables such that the manufactured solution verifies all the terms of the PDE in the desired flow regime and furthermore features a numerically desirable behaviour. For example, a suitable manufactured solution would yield an asymptotic (monotonic) range that is attainable on reasonably fine grids and that is nevertheless not trivially obtained on the coarsest grids. Finally, let's recall that the manufactured solution is not bounded by strict fidelity to physics since the verification is rather a mathematical exercise; it is however desirable to produce quantities that remain within the range of physically valid variations such that density and absolute temperature are defined positive for example.

The spatial domain for this set of manufactured solutions is Ω = [0,1]^2 and unless specified, all the grids are formed by isotropic quadrangular elements and refined by halving each element in each direction. The typical element size, h, is defined as 

    h=√(N_DOF),

where N_DOF=∑_e_i(N_node)_e_i is the total number of DOFs per equation.

In the next sections, we present manufactured solutions for free (unbounded) flows that target a specific set of solver capabilities. Each manufactured solution serves as well to exemplify and discuss important aspects of the verification process in general and in particular with attention to what is pertinent for high-order solvers. Extensive results are presented in <ref> as reference to facilitate comparison with other high-order codes.



 §.§ Inviscid flows in subsonic regime - MS-1


MS-1 is devised to verify both the incompressible (0.13≤Ma≤ 0.3) and  compressible (0.3 ≤Ma≤ 0.49) portions of the subsonic regime. The parameters of Eqs. <ref> that define this MS are given in Table <ref> and MS-1 fields are illustrated in Fig. <ref>. The results of the grid refinement study for polynomial degrees P1 to P5 are presented in Figs. <ref> in terms of discretization error in L_1, L_2 and L_∞ norms versus mesh size for all conservative variables. For each norm and P combination, the slope of the observed curve is compared to that of the expected OOAs, i.e., 𝒪(h^P+1).  It is as well notable that by increasing the degree of the polynomial discretization, P, the error level decreases for the same number of DOFs as in the previous case. This, in fact, is a demonstration of the higher accuracy/effort performance of high-order methods that constitutes their core interest.



For the same spatial distribution of the discretization error, each norm in Fig. <ref> outputs a different value. A general trend can be identified, with the L_∞ yielding the highest value, followed by the L_2 and L_1 norms that are closer together. This is expected since the L_∞ reports the maximal value of the error throughout the domain. We will exhaustively dwell on the importance of this norm  in the following sections. The L_1 norm measures the magnitude of the error, averaged over the domain, whereas the L_2 norm can be regarded as the magnitude of the error, weighed by itself first, and  averaged then. Hence the L_2 norm is more sensitive to the error deviations from the average value. Consequently, the disparity between the L_1 and L_2 norms reflects the irregularity of the spatial distribution of the error.

In terms of verification, the crucial verdict on grid convergence is declared as soon as the observed slope steadily recovers the expected order on at least three consecutive grids <cit.>. However, judging the satisfaction of this criterion merely from the error plots, such as the ones in Fig. <ref>, is an arduous task. It is therefore recommended to plot the evolution of the OOAs (the observed slope of the error) versus mesh refinement such as in Fig. <ref>, from which (see plot (c)), one can easily detect for example that in the case of P4, the ρ v variable has still not fully grid-converged on the fourth finest grid (h ≃ 6 × 10^-3). One can nevertheless acknowledge that this assessment is much harder to achieve by visually comparing the slopes in the plot (c) of Fig. <ref>. 




  §.§.§ Solution initialization and residual convergence


In order to study the evolution of the discretization error with mesh refinement in an isolated fashion, other sources of numerical error, such as iterative and floating point errors should be reduced as much as possible, by driving the norm of the residual equations as close to machine precision as possible. This is however a challenging goal in practical situations such as those encountered in external aerodynamics where solution non-linearity, mesh anisotropy and inadequate solution initialization contribute to hinder the convergence of iterative methods such as the Newton's. In the context of verification via the MMS, one could track the evolution of the exact discretization error during the convergence process. To provide insight on the adequate level of iterative convergence, we consider two scenarios: initialization with the exact[Let's note that for finite number of DOFs, the exact solution of the continuous equation is not the exact solution of the discrete equation.] (manufactured) solution, (Q^0=Q^MS), and initialization with an inexact solution, (Q^0=α  Q^MS  with α≠ 1 ).  Figure <ref> compares, in either of these cases, the evolution of the discretization error of ρ u as well as that of the absolute and relative residuals of the x-momentum equation in L_2 and L_∞ norms versus the number of iterations, n. For this particular combination of MS and spatial discretization (P5 and 32 × 32 elements), an absolute residual norm drop to values below ℛ_ρ u ≃ 10^-11 is sufficient for the discretization error to converge to its final value (ℰ_ρ u≃10^-10). In terms of relative residuals, the convergence of the discretization error is achieved as soon as the relative residuals are dropped below ℛ_ρ u ^n/ℛ_ρ u ^0 ≃ 10^-10 (by 10 orders of magnitude in other words) in the case of inexact initialization and by ℛ_ρ u ^n/ℛ_ρ u ^0 ≃ 10^-5 in the case of exact initialization. This indicates that a relative reduction in terms of orders of magnitude of the residual can not serve as a universal convergence criterion in verification exercises due to the inherent dependency of such criterion on the initial residual magnitude. We will dwell further on the necessary level of residual convergence for code verification while discussing laminar cases. 

The residual drop is understandably much faster achieved for the case of exact initialization that constitutes a sufficiently fine guess to ensure the stability and the quadratic convergence of the full Newton's method. Indeed, ensuring the iterative convergence in the case of inexact initialization required the relaxation of the Newton's method by the introduction of pseudo-time integration via the implicit Euler scheme. 
We hence choose to employ exact initialization throughout the rest of this work and consistently reduce discrete residuals to reach a stable minimum.


















 §.§ Inviscid flows in supersonic regime - MS-2


For MS-2, which verifies the implementation of the Euler equations and the farfield boundary conditions in the supersonic regime and on curved domains, we found that the values of parameters in Table <ref> result in grid convergence on reasonably fine grids.



The curved domain is defined via a deformation of the initial domain of (𝒳, 𝒴) ∈Ω_0 = [0,1]^2 (where isotropic Cartesian grids are created) to the target domain of (x,y) ∈Ω via

    x    = 𝒳 + 0.1 sin(π𝒳+π𝒴),
    
    y    = 𝒳 + 0.1 sin(π𝒳+π𝒴),

along with an isoparametric elemental mapping. Figure <ref> illustrates the domain and a typical grid of MS-2.



MS-2 fields are plotted in Fig. <ref>. The evolution of the errors and the orders versus mesh refinement are given in Figs. <ref> and <ref> respectively. The observed OOAs recover the theoretical values for all P and all norms considered except for L_∞ norm of P5 that is affected by the vicinity of the machine precision to the discretization error on the finest grid.

The OOAs of ρ u in L norms for the curved domain are compared to the ones for the original domain in Fig. <ref>, showing the occurrence of the asymptotic convergence on coarser grids in the latter case. Two factors are at the source of this difference: first, the domain of definition of the manufactured fields are different since the two domains do not recover each other; second, the curved grid does not exhibit the same isotropy with regards to the solution  as the un-deformed grid. The effect of grid anisotropy is dwelt on at length in Section <ref>. 

In order to avoid delays in the onset of asymptotic convergence with the goal of keeping the required grid sizes reasonable and since the treatment of the curved elements is verified by MS-2, the next manufactured cases will be considered on un-deformed domains along with Cartesian grids.






  §.§.§ Importance of L_∞ norm

The importance of the L_∞ norm in code verification has previously been discussed in the literature (see <cit.> and <cit.>) as to be a valuable metric for the detection of localized inconsistencies such as coding errors in the implementation of the boundary conditions. But the discussion has mostly stemmed from theoretical and practical appreciations of the qualities of this norm.
Here, we demonstrate the importance of the L_∞ norm by rather providing a concrete example: a bug affecting the implementation of the boundary condition at a single point is introduced by changing the boundary value of the ρ u variable of MS-2 at (x,y)=(0,0) such that (ρ u )^BC=1.000001 ×(ρ^MS u^MS) |_BC, while keeping all other boundary conditions at all other points intact in Eq. (<ref>). Fig. <ref> shows the evolution of the resulting discretization errors and OOAs versus mesh refinement for the polynomial degrees P1 to P3 on the un-deformed domain of MS-2. For P1, none of the error norms has yet been affected by the punctual boundary error on the finest mesh. Indeed, the discretization error of the P1 is rather large and drops slowly with mesh refinement (ℰ∼𝒪(h^2)) when compared to higher orders.  Hence, for the boundary inconsistency to be manifested in a P1 polynomial discretization, one needs to drive the mesh refinement process up to very fine grids for the level of discretization error to be lower than that caused locally by the boundary bug that is around ℰ_ρ u(0,0)=|ρ u(0,0)-ρ u_ex(0,0)|≃ 2× 10^-5, see Fig. <ref> (a). For P2 however, the abnormality is detected on the finest mesh but only by the the L_∞ norm. As the degree of the polynomial discretization is increased, the punctual source of stagnating error  dominates the norms for lower mesh sizes since the discretization error diminishes more rapidly on coarser grids for higher P. As such, the L_∞ norm of P3 is affected by the spurious boundary condition as off the third grid with a mesh size of h≃ 1.5× 10^-2 which is the coarsest grid to manifest the bug amongst all the polynomial degrees considered. But still, it takes another level of mesh refinement in P3 for the L_1 and L_2 norms to pick the error up. The effect of the mesh size on the discretization error is better illustrated in Fig. <ref> that compares for P3, the higher levels of discretization error, distributed regularly on a 16×16 grid, to the lower levels of discretization error on a 32×32 grid where  the manufactured bug at (x,y)=(0,0) is exhibited as the global maximum error on the domain. These observations support the idea that the L_∞ norm is a crucial metric to ensure the early detection of localized sources of error on coarse grids. Finally, it is interesting to  note in Fig. <ref> (b), that the error is being convected in the stream-wise direction from the punctual source.



  §.§.§ Importance of the inclusion of a high-order discretization

One other essential conclusion can be drawn from Fig. <ref> pertaining to the verification of high-order codes: it is desirable to include  a fairly high polynomial degree in the verification process to ensure the early detection of bugs causing errors with small magnitudes. Indeed, since the higher polynomials reduce the discretization error at a higher rate, the presence of the slightest bug can be revealed on coarser grids than if lower orders were used, thus providing a greater confidence in the absence of minor bugs in the code.







Since the Euler equations are verified for the supersonic regime as well, we can focus on the verification of the NS and RANS-SA equations in the subsonic regime. This relies on the fact that flow regime in terms of Ma directly affects the treatment of the inviscid portion of the equations and its corresponding boundary conditions whereas the treatment of the viscous terms is independent of the Ma number.



 §.§ Laminar flows - MS-3







In a first attempt, MS-1, that served to verify the Euler equations in subsonic and simultaneously compressible and incompressible regimes, was extended to viscous mode by adding the diffusion operator of the NS equations to the forcing functions. Although attaining the asymptotic range on reasonably fine grids for both lightly (μ = 1 × 10^-4) and highly (μ = 1 × 10^+2) viscous flows, MS-1 resulted in very delayed asymptotic ranges for intermediary values of dynamic viscosity (μ = 1 × 10^-1) as shown in Fig. <ref> (a) for errors in ρ. In this moderately viscous mode, the inviscid and viscous terms of the forcing functions of the momentum and energy equations take similar magnitudes in some regions of the domain whereas for more extreme values of μ, the forcing functions are dominated by either. The superposition of first and second order derivatives of the primitive variables, emanating respectively from the  inviscid and viscous terms, translates to the combination of a large number of sinusoidal functions that could explain the delay in the manifestation of the monotonic convergence range. In fact, to capture the multi-modal forcing function, finer grids are necessary, especially for high-order polynomial discretizations that exhibit a higher sensitivity to sub-cell spatial variations, even on coarse grids. Nevertheless, by tuning the trigonometric MS parameters, one could devise a numerically benign MS that yields a reasonably fast asymptotic range for largely distinct values of μ. The search for such properties resulted in the elaboration of MS-3, that verifies the NS equations in the compressible portion of the subsonic regime. MS-3 is defined by Eqs. <ref> along with the parameters of Table <ref>. The spatial distributions of its primitive variables and Ma number are illustrated in Fig. <ref>. For this MS, the results of the grid convergence study in terms of the error and order curves are respectively presented in Figs. <ref> and <ref>. Figure <ref> compares the evolution of OOAs for ρ versus mesh size for MS-1 and MS-3. One can appreciate the faster convergence (in h) of MS-3, that is the result of the proper tuning of its parameters. It is interesting to note that although the convergence of the higher orders is delayed for MS-1, its observed orders for P1 converge to the value of 2.0 on rather coarse grids, demonstrating the particular difficulty in deriving numerically benign manufactured solutions for the verification of high OOAs. 










  §.§.§ Balancing of the forcing function terms

In realistic flow problems, the Reynolds number (Re) provides a measure of the relative magnitude of advective over the diffusive forces, that is also related to the ratio of the inviscid over the viscous terms of the NS equations. The magnitude of the viscous terms depends on both the value of the viscosity and the strain rate which involves the velocity gradients. The latter are intrinsically adjusted by flow features such as the boundary layer thickness, the free-stream  velocity, the characteristic length, etc. that are themselves correlated with Re (and consequently with μ). However, since a manufactured solution is not necessarily realistic, it might be absent of these features and the velocity gradients could take arbitrary values independently of the value of viscosity. Hence, to generalize the concept of Reynolds number to manufactured cases, we directly look into S^rel(𝐱) = S^inv(𝐱)/S^vis(𝐱), the ratio of the inviscid over the viscous terms of the manufactured forcing functions in the domain as in <cit.>. Figure <ref> compares the S^rel distributions for the x-momentum equation of MS-3 and three different non-dimensional values of dynamic viscosity, i.e.,  μ = 1 × 10^-4, μ = 1 × 10^-1 and μ = 1 × 10^+2. Using the same MS fixes the velocity gradients while  μ is changed. To ensure the verification of both the inviscid and viscous terms by the MS, the S^rel(𝐱) distribution should simultaneously include regions of viscous preponderance (|S^rel(𝐱)|≪ 1) and inviscid prevalence (1 ≪|S^rel(𝐱)|) along with non-negligible regions of equivalence (|S^rel(𝐱)|≈1). This is an important criterion, since when one mode outweighs the other by much, the verification process becomes biased and the detection of minor bugs by the MMS in the outweighed mode is jeopardized. Figure  <ref> suggests that for MS-3, a value of μ = 1 × 10^-1 guarantees a balanced leverage since for this value, non-negligible regions of inviscid and viscous dominance exist in different regions of the domain simultaneously.



To further demonstrate the importance of balancing different terms of the manufactured forcing functions in the verification process, we have introduced a bug by modifying the heat flux in Eq. (<ref>), such that ω_j =-αλ_eff ∂_j T, where the value of the factor α is changed from its original value of α=1.0 to α=1.0001. This bug can go undetected for μ=1× 10^-4 and only revealed when the viscosity is increased to μ=1× 10^-1 as shown in Fig <ref> for the error in ρ E of MS-3, discretized by P3.





  §.§.§ Residual convergence

The residual convergence is affected by the inviscid/viscous ratio in two fashions: firstly, as μ is increased for the same MS, the convergence of iterative linear system solvers such as GMRES is hindered since the prevalence of the viscous terms reduces the diagonal dominance of the system and consequently the system conditioning degrades. This translates to an increase in the number of GMRES iterations needed to satisfy its internal convergence criterion. Secondly, for an exact initialization, the initial residual takes larger values as μ is increased. This effect is shown in Fig. <ref> that compares the residual and discretization error convergences of MS-3 for viscosity values of μ = 1 × 10^-4, μ = 1 × 10^-1, μ = 1 × 10^+2 and μ = 1 × 10^+4 discretized by P5 and a grid of 32 × 32 elements. The increase in the initial residual with μ reflects the magnifying effect of the viscosity on the truncation error of the viscous terms since any norm of the discrete residual computed using the exact solution is in fact a measure of the truncation error (see Eq. <ref>). 

Once the viscous terms dominate the residual equation, increasing the viscosity any further just scales the level of the residual norm as it ensues from the comparison of residuals for μ = 1 × 10^+2 and μ = 1 × 10^+4 in Fig. <ref>. It is noteworthy, in the same figure, that although the residual level is sensitive, the norms of the discretization error are quasi invariable with regards to viscosity. This suggests that the level of the discretization error has a higher sensitivity to the construction of the manufactured variables and their derivatives than to the terms of the governing equations considered in this study.



Another important observation can be made form Fig. <ref> with regards to the required level of iterative convergence. It indeed is argued in <cit.> that the iterative errors scale closely with the level of the discrete residual in a wide range of flow problems and in <cit.> that the iterative errors should be driven down to two to three orders of magnitude lower than the discretization error with the goal of ensuring that the iterative and round-off errors are negligible and do not affect the verification process. Hence, from Fig. <ref>, it is obvious that the final value of the residual at convergence is a function of μ and could be higher than the discretization error. Consequently, it seems that a convergence criterion requiring the residual norm to be lower than the discretization error is overly conservative and in fact unrealizable in some cases. Based on our experience, for numerically benign problems such as the free flows considered in this work, reaching the minimum residual possible, by a fully implicit method via exact initialization, is sufficient for the realization of the expected OOAs. For stiff problems and in the case where the OOAs can not be achieved, one could rather verify their sensitivity to the residual convergence by under-converging a few cases by one or two orders of residual magnitude and recomputing the OOAs. If the latter change significantly, it could be an indication of a need for further residual convergence.



  §.§.§ H_1 semi-norms

In addition to the previously mentioned metrics, the verification of the viscous terms can be complemented by tracking the order of convergence in the H_1 semi-norm that monitors the discretization errors in state derivatives. For MS-3 with μ = 1 × 10^-1, the evolution of the errors in this metric versus mesh refinement  is provided in Fig. <ref> for both the uncorrected (∂_q Q_k) and fully corrected (∂_q Q_k)  derivatives, for which the semi-norms are respectively labelled H_1 and H_1. One can firstly appreciate that whenever there is a significant difference between them, the fully corrected derivatives are more accurate than the uncorrected derivatives. Secondly, the error in this semi-norm recovers convergence slope of 𝒪(h^P) as expected. Note that the regularity of its convergence to the theoretical OOAs can be better assessed in Fig. <ref>. 

In order to evaluate the relative importance of H_1 semi-norms versus the L norms, we have investigated the ability of each in detecting various types of spurious alterations to the implementation. In all cases, the performances of these types of norms were comparable and as such H_1 semi-norms showed no particular prevalence in exhibiting the presence of the bugs.

Figure <ref> shows that as the viscosity is lowered from μ = 1×10^+2 to μ = 1 × 10^-4, and consequently the diffusive terms loose their relative magnitude compared to the advective terms, the corrected derivatives undergo a significant delay in the attainment of the asymptotic range for P3 and P5 discretizations, whereas the OOAs in L norms, illustrated in Fig. <ref>, are not noticeably affected. Note in the same Fig. that there is a drop in the L norms of the finest P5 solutions which is due to the proximity of the error level to machine precision. The presence of large discrepancies between the OOAs of corrected and uncorrected derivatives in H_1 semi-norm seems correlated to the insignificance of the relative magnitude of viscous terms compared to inviscid terms and hence suggests that an examination of the balancing of the forcing function terms is required.























 §.§ Turbulent flows - MS-4 (original SA) and MS-5 (modified SA)





  
MS-4



We consider an extension of MS-3 to ensemble-averaged turbulent flows by activating the ν̃ field, first such that the original portion of the SA model (Eqs. <ref>, <ref>, <ref>, <ref>, <ref> and <ref>) is verified for positive ν̃ values. This manufactured solution is labelled MS-4 and is determined by parameters presented in Table <ref> along with Eqs. (<ref>) and a value of dynamic viscosity of μ = 1 × 10^-3 chosen such that the eddy viscosity be preponderant in comparison. The distance to the wall is set to d_w=y+1, for y ∈ [0,1], to avoid large values in SA source terms with dependencies on 1/d_w^2, that could otherwise saturate the entire forcing function. 

The balancing of different SA terms is discussed at length further in this section. The manufactured ν̃ field of MS-4 is presented in Fig. <ref>. The other fields are the same as those of MS-3 (Fig. <ref>). The errors and OOAs are respectively presented in Figs. <ref> and <ref> for L norms and in Figs. <ref> and <ref> for H_1 semi-norm. 




  §.§.§ Balancing of the forcing functions

As discussed previously for laminar flows, a balance between different terms of the forcing function is desirable with the goal of enabling the detection of slightest bugs on coarse grids, thus reducing the risk of a false pass verdict. We hence extend the budget analysis that involves the relative magnitudes of inviscid versus viscous terms in laminar flows, to include all the terms of the SA source function. In the laminar case, achieving a balance between the advective and diffusive terms is a relatively easy task since the value of μ directly scales their ratio. For turbulent flows however, the presence of source terms within the RANS model, involving  a range of quantities, makes this balancing a more challenging task. For the assessment of the sought balance, we propose considering a relative absolute magnitude, S^rel_term = |S_term| / S^sum_SA, as metric of the sensitivity of the verification process to each SA term where |S_term| is the absolute magnitude of the SA term and

    S^sum_SA = ∑_term|S_term| =     []∂_j (ρ u_j ν̃)_Advection  +
    []- ∂_j (1/σ(μ+ρν̃)∂_j ν̃)_Diffusion+
    []ρ  c_b1(1-f_t2)  s̃ ν̃_Production +  
       []-ρ(c_w1 f_w -c_b1/κ^2f_t2)ν̃^2 /d_w^2_Destruction + 
    []c_b2/σρ ∂_j ν̃ ∂_j ν̃_Distribution + 
    []1/σ(ν+ν̃f_n)  ∂_j (ρ ∂_j ν̃)_Conservation.


The six terms of the SA forcing function appearing in Eq. (<ref>) are those of the original portion of the SA model, Eqs. <ref>, <ref>, <ref>, <ref>, <ref> and <ref>.
For each of these terms, the spatial distribution of the S^rel_term is provided in Fig. <ref> showing that the advection, diffusion and destruction terms exhibit  spatial maxima of max(S^rel_term) ≈50% on the domain whereas the sensitivity to production, distribution, and conservation terms only accounts to  max(S^rel_term)≈ 5%. To establish whether such a low sensitivity is acceptable and how it relates to actual sensitivity to a given bug, we introduce an error in the distribution term of the SA model, (1+dα) c_b2/σρ ∂_j ν̃ ∂_j ν̃, such that the original term (dα=0) is modified sequentially by increasing dα from dα=1 × 10^-15 by an order of magnitude (× 10) in each step and the OOAs are recomputed until the bug is detected. Consequently, a significant degradation in the orders of the P5 discretization was detected for a value of dα=1 × 10^-7 (see Fig. <ref>) demonstrating that considerably small inconsistencies in the implementation of the distribution term can be detected and hence max(S^rel_term)≈5%  is a sufficient criterion to ensure that all terms are correctly verified. From Fig. <ref>, it can again be  emphasized that including a reasonably high (such as P5) polynomial degree in the verification process is crucial for the detection of minor inconsistencies in the implementation. 





Based on the sensitivity analysis presented and the results obtained for OOAs in Figs. <ref> (L norms) and <ref> (H_1 semi-norm), it is possible to conclude that the considered implementation of the original portion of the SA model is verified via  MS-4 and hence the focus can be moved on the modified portion of the SA model. 



  
MS-5
We propose MS-5 for the verification of the modified portion of the SA model; recalling that the revised model is intended to tackle the instabilities generated when negative values of the SA working variable occur in the numerical solution. Like MS-4, MS-5 is an extension of MS-3 by that it inherits the latter's manufactured fields which it complements by a manufactured solution for ν̃.



Table <ref> presents the parameters that define MS-5. The dynamic viscosity is set to μ=1 × 10^-1. For negative ν̃ values, the RANS equations are decoupled from the SA equation by eliminating the effect of the latter on the former via μ_t = 0. However, the SA equation remains unilaterally coupled to the RANS equations. Therefore, MS-5 only verifies the implementation of the SA model and regarding the state variables of all other equations, it should result in the same exact error and orders as observed for MS-3 in Figs. <ref>, <ref> for L norms and in Figs. <ref> and <ref> for H_1 semi-norm.  As for the SA equation, its manufactured field as well as errors and orders in L and H_1 norms are respectively presented in Figs. <ref>, <ref> and <ref>. 



A similar budget analysis of the SA forcing function to the one for MS-4 is conducted here by rather considering the modified portion of the SA model (Eqs. <ref>, <ref>, <ref>, <ref>, <ref> and <ref>) such that:

    S^sum_SA = ∑_term|S_term| =     []∂_j (ρ u_j ν̃)_Advection  +
    []- ∂_j (1/σ(μ+ρν̃c_n1+χ^3/c_n1-χ^3)∂_j ν̃)_Diffusion+
    []ρ  c_b1(1-c_t3)   s  ν̃_Production +  
       [] c_w1ν̃^2 /d_w^2_Destruction + 
    []c_b2/σρ ∂_j ν̃ ∂_j ν̃_Distribution + 
    []1/σ(ν+ν̃c_n1+χ^3/c_n1-χ^3)  ∂_j (ρ ∂_j ν̃)_Conservation.


This analysis resulted in the spatial sensitivity distributions of Fig. <ref> showing that in descending order of magnitude, the SA forcing function depends on the diffusion, destruction,  advection, conservation, distribution and production terms. The sensitivities to the latter term are very low, max(S^rel_term)≈0.1%,  and are particularly under the criterion of max(S^rel_term)≈5%, established previously. It is hence prudent to conduct an actual sensitivity analysis by introducing an inconsistency in the implementation of the production term of the SA model as (1+dα) (ρ  c_b1(1-c_t3)   s  ν̃) and computing the OOAs for an initially small and increasingly larger values of dα until a difference in the orders is observed. This exercise resulted in the detection of the bug as soon as dα = 1 × 10^-6. Therefore, the production term is verified by MS-5 even for minor bugs and MS-5 can be considered for the verification of the modified SA equation. 

The only term of the governing equations that remains unverified throughout the manufactured cases presented so far is the modified vorticity magnitude, s̃ of Eq. <ref>. In fact, this is the only term of the modified SA that is not activated by a condition on the sign of ν̃ values but rather by a criterion on the vorticity magnitude, s, that is often satisfied in the vicinity of solid walls. The verification of this term is hence better suited for wall-bounded manufactured cases.





  §.§.§ Effect of grid stretching and self-similarity


Manufactured solutions based on trigonometric functions do not feature sharp gradients and thereby constitute an ideal candidate for code verification on isotropic grids. Nevertheless, mesh stretching (element anisotropy and non-uniformity of element size distribution) is often mandatory in the case of realistic flows such as  boundary layers, jets or shocks in order to tackle sharp variations while ensuring residual convergence and solution accuracy. Determining the necessary and sufficient level of stretching  in every region of the domain is a complex subject. For RANS, one could adopt a heuristic approach by following general guidelines on the wall element size such as y^+ ≤5 where y^+ is a dimensionless quantity reflecting the flow intensity in the vicinity of solid wall. Although sufficient to ensure the validity of the numerical solution with regards to the turbulence model, this criterion does not guarantee that the grid is appropriate to realize the optimal OOAs for integrated quantities at the wall or for the smoother solution in outer boundary layer region. Furthermore, generating strictly self-similar grid sets for grid convergence study and uncertainty analysis is not always possible for complex geometries. Self-similarity of a grid set refers to the property of featuring a constant refinement ratio of the grid resolution (typical element size), h, throughout the domain, for any two consecutive grids of the set. Note that the effect of various common measures of h on the appearance of the asymptotic convergence rate and the achievement of optimal OOAs is discussed in <cit.> for low-order discretizations and self-similar and dissimilar grid sets. We choose the definition (<ref>) that is shown in <cit.> to produce the optimal rates for self-similar grid sets. In this section, we investigate the effect of grid stretching  as well as the impact of self-similarity of a grid set on the OOAs of trigonometric solutions. 

The effect of stretching and self-similarity of a sequence of grids  on the OOAs produced by a grid convergence study is analyzed in <cit.> for a low-order solver and L norms. We extend the analysis to higher OOAs and H_1 semi-norm by considering:


  * the effect of self-similarity on uniform grid distributions by applying two random perturbations of element sizes, with increasing amplitudes, around a uniform element size distribution;

  * the effect of grid anisotropy and non-uniformity on trigonometric manufactured cases by generating an expanded element size distribution via a geometric series with fixed size ratio of r=1.015;

  * the effect of self-similarity on stretched grids by applying two random perturbations of element sizes, with increasing amplitudes, around the expanded element size distribution via a geometric series with fixed size ratio of r=1.015.


In all these cases, the y coordinates are produced by replicating the x coordinates such that the resulting grid is symmetrical. Coarser grids are generated from the finest grid of the sequence by removing every other grid line repeatedly until the desired level of coarsening is reached. As a result, the size distribution pattern of the coarser grid is  similar to the finer grid, strictly if no random perturbation is applied. The random perturbation is applied to the  element e_i with size Δ x_i as Δ x'_i = α_i  Δ x_i
where Δ x'_i is the perturbed size and α_i is a random number bounded by 1-dα_m≤α_i≤1+dα_m. We consider three values of perturbation intensity, dα_m:


  * dα_0 = 0, the unperturbed grid set;

  * dα_1 = 0.25, moderately perturbed grid set;

  * dα_2 = 0.5, most intensely perturbed grid set.


Figure <ref> shows  the resulting six element size distributions corresponding to the finest grid of each set (uniform and expanded with three perturbation intensities each). MS-5 is solved on each grid set and the normalized discrepancy between observed, 𝒪(h^o), and  theoretical, 𝒪(h^P+1), OOAs for P1–P4 and ρν̃ is measured by |(P+1)-o|/(P+1) and tracked as the mesh is refined. The results are presented in Figs. <ref>–<ref> and leading to the following remarks: 


  * The observed order discrepancy in L_2 norm for the unperturbed uniform set (dα_0 in Fig. <ref>) does not reduce as regularly in the present two-dimensional case as in a one-dimensional case <cit.>. This suggests that the discretization error in multiple spatial dimensions does not scale as smoothly with mesh refinement as in a single dimension.



  * As the uniform grid is perturbed, the order discrepancy in L (Fig. <ref>) and H_1 (Fig. <ref>) norms often increases rather than reducing with mesh refinement. This trend is more pronounced as the perturbation intensity, dα, is increased. In other words, the order departs from the asymptotic range with mesh refinement under the effect of grid dissimilarity.

  * The perturbations of the uniform grid set often affect the order discrepancy in L_∞, L_2 and L_1 norms most adversely in the descending order (Fig. <ref>). This demonstrates the higher sensitivity of L_∞ norm followed by L_2 norm with regards to grid self-dissimilarity that is manifested here through the irregularity of discretization error reduction rate.

  * The H_1 norm of corrected derivatives (∂_q Q_k) undergoes larger order discrepancies than that of uncorrected derivatives (∂_q Q_k) for both unperturbed uniform (Fig. <ref>) and unperturbed expanded (Fig. <ref>) sets. Hence, the corrected derivatives reach the asymptotic range with a delay with regards to the uncorrected derivatives.

  * Although larger in magnitude compared to the unperturbed uniform grid set (dα_0 in Figs. <ref> and <ref>), the order discrepancy of the unperturbed expanded grid set (dα_0 in Figs. <ref> and  <ref>) diminishes with mesh refinement for all L and H_1 norms (except for the L_∞ norm of the P2 discretization). This illustrates that the unnecessary application of grid anisotropy delays the asymptotic range but in contrast with grid self-dissimilarity, it does not disrupt it.

  * As in the case of the uniform set, the random perturbations of the expanded set disrupt the asymptotic convergence for both the L (Fig. <ref>) and H_1 norms (Fig. <ref>) and more so as the perturbation increases from dα_1 to dα_2.


In the light of these observations, it can be concluded that for higher OOAs, consistently with low-order discretizations, the inadequate application of grid stretching (anisotropy and non-uniformity), to regions where solution gradients are rather smooth, only delays the occurrence of the asymptotic range, whereas if the grid set is not self-similar, the monotonic convergence of the OOAs is compromised and hence grid self-similarity needs to be enforced in a grid convergence study. 






§ CONCLUSIONS

The aim of this article is to extend the code verification methodology to high-order accurate  CFD solvers. To this end, a series of trigonometric manufactured solutions is introduced which enables the incremental examination of CFD solver components for the simulation of free flows in inviscid, laminar and turbulent (via the original and modified SA models) regimes. The capability of each of these MSs for the verification of OOAs from the first up to the sixth is demonstrated on a pre-integrated DG implementation via the FR/CPR scheme. These MSs served as well to reach the following conclusions:


  * The L_2 norm is shown to be more sensitive to deviations from the mean than the L_1 norm whereas the L_∞ norm is demonstrated to be the most useful in detecting localized inconsistencies and hence its inclusion in code verification is recommended;

  * The H_1 semi-norm of uncorrected and fully corrected (via BR2 scheme) derivatives is considered and the latter are found to be often more accurate than the former. Note however that among many attempts via introduction of various types of bugs,  no case showed a significant advantage  of H_1 semi-norms over L norms in their detection;

  * The inclusion of a relatively high OOA (such as the sixth) in code verification is recommended since its utility is demonstrated for the detection of minor bugs  which affect higher orders first on coarser grids;

  * A supersonic inviscid MS on deformed domain is introduced to verify the treatment of curved elements;

  * Manufactured cases are presented to verify the original and modified SA models coupled to compressible RANS equations. The importance of the balancing of forcing function terms in assuring the reliable examination of all terms of the governing equations is illustrated via an example. A sensitivity analysis method is introduced to enforce the conclusions of the verification exercise in cases where some terms are outweighed by others, thus providing assurance in the outcomes of the verification campaign. Let's note that the only terms of the governing equations left unverified through this study are those of the modified vorticity (Eq. (<ref>)) of the modified SA model;

  * The necessary level of residual minimization for the convergence of the discretization error is investigated and it is shown that for high μ values in a laminar manufactured case, the residual magnitude is larger than that of the discretization error, indicating that it is not always necessary to drive the residual level lower than the discretization error for the iterative and round-off errors to be negligible;

  * The effect of stretching and self-similarity of a grid set on grid convergence of solutions with smooth gradients is studied numerically for a turbulent MS and it is concluded that unnecessary grid stretching delays the occurrence of the asymptotic convergence whereas the loss of self-similarity disrupts it.


The discussion on the extension of verification methodology of high-order solvers on wall-bounded curved domains and realistic turbulent flows is intended to be pursued in a future work.







§ ACKNOWLEDGEMENTS

The authors gratefully acknowledge the generous support from the Fonds de recherche du Québec – Nature et technologies (FRQNT), the Natural Sciences and Engineering Research Council (NSERC), and the Department of Mechanical Engineering of McGill University.
























§ DEFINITION OF THE GRID CONVERGENCE METRICS


Here are the definitions of the norms employed throughout this work to measure the discretization error. The integrals are computed by GLL quadratures.




  * L_∞ norm:

    ℰ_Q _∞ =
    max| Q_i - Q_i^ex|  for  i ∈ [ 1 .. N_DOF]



  * L_1 norm:

    ℰ_Q _L_1 =
    ∫_Ω | Q - Q^ex|  dΩ/∫_Ω dΩ



  * L_2 norm:

    ℰ_Q _L_2 = ( ∫_Ω ( Q - Q^ex)^2dΩ/∫_Ω dΩ)^1/2


  * H_1 norm:

    ℰ_Q _H_1 = ( ∫_Ω ( Q - Q^ex)^2  dΩ +  ∫_Ω ∑_q=1^N_d(∂_q Q - (∂_q Q)^ex)^2   dΩ/∫_Ω dΩ)^1/2


  * H_1 semi-norm of uncorrected gradients:

    []ℰ_Q _H_1 =(∫_Ω ∑_q=1^N_d(∂_q Q - (∂_q Q)^ex)^2   dΩ/∫_Ω dΩ)^1/2


  * H_1 semi-norm of fully corrected gradients:

    []ℰ_Q _H_1 =(∫_Ω ∑_q=1^N_d(∂_q Q - (∂_q Q)^ex)^2   dΩ/∫_Ω dΩ)^1/2










§ MANUFACTURED SOLUTION RESULTS







 §.§ MS-1












 §.§ MS-2











 §.§ MS-3


















 §.§ MS-4

















 §.§ MS-5














§ GRID SENSITIVITY RESULTS























