1.1.1. Linear correction

Corrections have been addressed by different means in the past with the dual weighted residual (DWR) method (Becker and Rannacher 1996; Giles and Suli 2002a; Venditti and Darmofal 2002; Jones et al. 2006; Leicht and Hartmann 2010; Fidkowski and Darmofal 2011), and more generally the goal-oriented method (Loseille et al. 2010b) and the error transport equations (ETEs) (Layton et al. 2002; Pierce and Giles 2004; Hay and Visonneau 2006; Derlaga and Park 2017). DWR/goal-oriented methods focus on the correction of a functional. In some cases, the corrected functional is a higher order approximation of it. ETE methods focus on solution fields. Both rely on the same principles. As concerns DWR/goal-oriented methods, the linearization of the functional reads

[1.4]

Since the projection Π_hW is unknown, we introduce a finer mesh h/2. Put directly, we have

In our standpoint, computing the fine-grid solution W_h/2 for evaluating a corrector is too costly. Conversely, knowing W_h, the evaluation of Ψ_h/2(Π_h/2W_h) is affordable, and we have

[1.5]

where is the Jacobian of Ψ_h/2. From relations [1.4] and [1.5], we obtain

where is the adjoint of the problem. In order to avoid the computation and the inversion of the large matrix A_h/2, the functional j and the numerical system are assumed to have a similar behavior on coarser and finer meshes so that is replaced by W_h, its analog on h:

Then, the RHS is a corrected approximation of j_h(W_h), which is obtained with the extra evaluation of the adjoint state and the finer grid residual .

Similarly, the ETE method computes the corrected solution field from the linearization of the problem as

where is a defect residual obtained by different reconstruction means to estimate the error introduced by the discretization. The linearization of the problem is a determining factor for the quality of the corrected solution (Yan and Gooch 2015; Yan and Ollivier-Gooch 2017). This is true for both approaches. In particular, an approximate Jacobian obtained by passing to a lower spatial order of accuracy (as used in the Navier–Stokes solver introduced in Chapter 1) should not be used as it will lead to a poor correction. Lastly, would the computation of the matrix even be exact, these two approaches are limited by the fact that the residual is linearized.

1.1.2. Nonlinear correction

Let us consider the ETE iteration with a RHS evaluated on a finer grid:

where is a transfer from finer grid h/2 to grid h. If the updating can be reiterated and converged, we can gain the following advantages:

– the implementation is simple when we already have a Jacobian;

– it does not require a more accurate approximate Jacobian than the one of the flow solver;

– it automatically takes into account all numerical features and nonlinearities of the solver (limiter, flux reconstruction, Riemann solver, etc.).

This chapter is organized as follows. We first discuss in section 1.2 two ways of defining correctors for the usual P₁ approximation of a Poisson problem. We then recall in section 1.3 the notations related to the RANS equations and to their discretization. In section 1.4, the nonlinear discrete corrector for that scheme is explained. In section 1.5, we show the application of the finite volume nonlinear corrector to an inviscid supersonic flow simulation for a Sonic Boom prediction workshop.

1.2.1. Notations

Let Ω being a sufficiently smooth computational domain of . The continuous PDE system is written as

[1.6]

where

and where is a positive, possibly discontinuous, scalar field in Further, we assume that the bilinear form a is coercive in space V , that is, there exists a positive α such that

This model exemplifies the pressure equation in some multi-phase incompressible flow formulation for which mesh adaptation is useful (e.g. Guégan et al. 2010). Let Ω_h =Ω for simplicity, τ_h a triangulation of Ω_h and V_h be the usual P₁ -continuous finite-element approximation space related to τ_h:

The finite-element discretization of [1.6] is written in variational and operational form:

[1.7]

in such a way that u_h is a linear function of f_h, which we denote .We denote by Π_h the usual interpolation operator:

Scalar correctors, that is, correctors for scalar outputs j(u_h) depending on the solution, for example, j(u_h)= (g, u_h) with g prescribed, have been defined by Giles and Pierce (1999). Our interest concerns the correction of the unknown field itself. Two options are now proposed.

1.2.2. A priori corrector for the PDE solution

A first rather simple a priori corrector can be derived from an analysis of the error RHS. We observe that:

Assuming that the solution u is continuous, we can apply to u the usual P₁ interpolator Π_h and get:

[1.8]

We call the implicit error. By implicit, we mean that it can be obtained through solving a discrete system. It differs from the approximation error by an interpolation error:

[1.9]

In order to find an approximate of the implicit error, we need to evaluate the RHS of [1.8] for any test function φ_h. The second term of RHS of [1.8] is easy to evaluate (we know f and f_h). The first term of RHS of [1.8] can be transformed as follows (assuming that is constant on each element):

Then, we get

[1.10]

where the sum is taken for all edges ∂T_ij separating triangles T_i and T_j of the triangulation. The unit vector n_ij normal to ∂T_ij is pointing outward T_i.Now, we do not know u but u_h. In order to evaluate the interpolation error, we first approximate the Hessian of u by an approximation H_h(u_h) in V_h of the Hessian of u_h. We apply the recovery method of Zienkiewicz-Zhu type (1992) described in section 5.4.3 of Volume 1. This produces the values at vertices of the approximate Hessian H_h(u_h). Then, Π_hu − u is approximated by the continuous interpolation error π_hu_h − u_h defined in section 5.4.6 of Volume 1. This P₂ by element function is zero at vertices and at mid-edges ij. We apply the same type of estimate to f − f_h, namely, In order to approximate the implicit error Π_hu − u_h, solution of [1.10], we define our a priori implicit corrector as follows:

[1.11]

This corrector is an approximation of Π_hu − u_h. Replacing again the unknown interpolation error by its evaluation on the discrete approximation, we define our a priori corrector as follows:

[1.12]

Assuming that the approximations made between [1.10] and [1.11] are small, the a priori corrector is then built in such a way that (formally)

This corrector is easy to compute, but we shall see in Chapter 6 of this volume that its accuracy is rather low.

1.2.3. Finer-grid DC corrector for the PDE solution

A second option relies on using a fictitious finer grid. Here, we write it in the linear case. Let us assume that the approximation is in its asymptotic mesh convergence phase for the mesh Ω_h under study of size h. Then, this will be also true for a strictly two times finer embedding mesh Ω_h/2 and, for our second-order accurate scheme applied to a sufficiently smooth problem, the heuristics of Richardson analysis is written as

[1.13]

where u_h and u_h/2 are, respectively, the solutions on Ω_h and Ω_h/2. Note that in the case of local singularities, second-order convergence [1.13] is not true for uniform meshes since global convergence degrades to first-order convergence. But Chapter 2 of Volume 1 tends to show that it holds almost everywhere for a sequence of adapted meshes (see also Loseille et al. 2007).

Then, at least formally, we have

which implies

Now,

where P_h→h/2 linearly interpolates coarse values on fine mesh. This allows to evaluate Π_hu − u_h but needs to solve a finer grid system. Let us introduce the residual transfer R_h/2→h, which accumulates on coarser grid vertices the values at finer vertices in neighboring coarse elements multiplied with barycentric weights. In order to reduce the computational cost to solving a coarser grid system, we approximate by

This motivates the definition of a finer grid DC corrector as follows:

[1.14]

The DC corrector approximates Π_hu − u_h instead of u − u_h and can be corrected as the previous one:

[1.15]

1.3.1. Vector form of the RANS system

The compressible Navier–Stokes equations for mass, momentum and energy conservation with Spalart–Allmaras turbulence model are those described in Chapter 1 of Volume 1. It is useful to re-write the RANS system in the following more compact vector form:

(+ boundary conditions), where (idem for S). W is the non-dimensionalized conservative variables vector:

are the convective (Euler) flux functions:

[1.16]

S(W) = (S₁(W),S₂(W), S₃(W)) are the viscous fluxes:

[1.17]

where T_ij are the components of laminar stress tensor :

Q(W) are the source terms, that is, the diffusion, production and destruction terms from the Spalart–Allmaras turbulence model:

[1.18]

Note that Q = 0 in the case of the laminar Navier–Stokes equations, additional source terms are added (to take into account gravity, for instance).

1.3.2. Formal discretization

1.3.2.1. Finite element discretization

In the standard Galerkin (possibly stabilized) approach, the strong form of the equations over a given domain Ω is reshaped in a weak form:

with additional boundary conditions and where V (Ω) = = is the functional space in which the solution is sought. The problem is discretized by restricting the functional space to a subset V_h(Ω) of V (Ω). In the standard finite element approach, the domain Ω is split in a finite number of discrete elements K_i (triangles, tetrahedra, hexahedra, prisms, etc.) on which the restriction φ_h|Ki is a polynomial for each function φ_h ∈ V_h(Ω) and so that ∪_iK_i =Ω. The value of φ_h|Ki on each element is determined by a set of given nodal values and the functions φ_i – defined by setting the ith node to 1 and the other to 0 – form a basis of V_h(Ω). The residual vector of the discrete problem is defined by

1.3.2.2. Finite volume discretization

In the standard finite volume approach, the domain Ω is split in a finite number of cells ∪_iC_i = Ω and the residual is defined as

for any continuous field W, where Cⁱ is the cell associated with vertex P_i. The problem is discretized defining the vector W_i = ∫_C WdΩ, supposed to be the mean value of the field on each cell, and a reconstructionⁱ operator U that recovers a field W on each cell from the discrete values W_i. For flow problems, the advection terms can be written as a divergence; thus using the Green formulae, the integral of advection term is the integral of the advection fluxes over the boundaries of the cell:

where n is the inward normal and U a reconstruction operator that recovers a field W on each cell from the discrete values W_h.

1.3.3. Notations for discretization

The spatial discretization of the RANS model is described in Chapter 1 of Volume 1. We give here some details of the implementation of our in-house finite volume– finite element flow solver Wolf as this induces several constraints – as we will see – on the implementation of the corrector. In the sequel, we denote for a given vertex P_i:

– E(P_i) the set of edges joining P_i to another vertex P_j of the mesh;

– V¹(P_i) the set of vertices that are neighbors of P_i, that is, the set of vertices connected to P_i by an edge. V¹(P_i) is named the first-order ball of P_i.

– V²(P_i) the set of vertices that are neighbors of the neighbors of P_i. V²(P_i) is named the second-order ball of P_i. Similarly, V^k(P_i) is the kth iteration of this process.

The discretization of Euler fluxes are sufficiently described in Chapter 1 of Volume 1. Concerning the viscous terms using the P₁ finite element method (FEM), we use the following notations:

where ∂C_ij is the common interface between cells C_i and C_j . Let φ_i be the P₁ finite element basis function associated with vertex P_i, we have

and as the solution is represented on the P₁ basis, S(W_i) (which comes from a gradient) is assumed constant by parts on each element K, for example,

where μ|_K is the mean value of μ on the element to be constant. Then, we obtain

The effective computation of the previous integral then leads to the computation of integrals of the following form:

In this expression, ∇φ_i|_K is the constant gradient of basis function φ_i associated with vertex P_i. In the facts, these terms are computed element-wise and each contribution is added to each vertex of the element. Given an element K =(P_i,P_j , P_k, P_l), we have the partial flux associated with each vertex:

1.4.1. Finite volume nonlinear corrector

The central idea is to approximate Ψ(W_h) by Ψ_h/n(W_h). To accumulate on the h-mesh the residual computed on the h/n-mesh, we split the accumulation process depending on the possible configurations. If we consider a h-mesh vertex P_i, for the accumulation of the residual at P_i the different cases are as follows:

– the h/n-mesh vertex Q_k coincides with the h-mesh vertex P_i, then the full residual is accumulated: ε_h/n(Q_k) = ε_h/n(P_i);

– the h/n-mesh vertex Q_k belongs to a h-mesh edge e_j containing P_i, then we compute the barycentric coordinates β_i(e_j ) of Q_k w.r.t. P_i on edge e_j (Alauzet 2016). The accumulated residual is β_i(e_j )ε_h/n(Q_k);

– the h/n-mesh vertex Q_k belongs to a h-mesh face f_j containing P_i, then we compute the barycentric coordinates β_i(f_j ) of Q_k with respect to P_i on face f_j . The accumulated residual is β_i(f_j )ε_h/n(Q_k );

– the h/n-mesh vertex Q_k is inside a h-mesh element K_j containing P_i, then we compute the barycentric coordinates β_i(K_j ) of Q_k with respect to P_i in element K_j . The accumulated residual is β_i(K_j )ε_h/n(Q_k).

For the h-mesh vertex P_i, the accumulation process is expressed as

[1.19]

For n = 2, this gives Algorithm 1.1. A 2D example of this accumulation process at the element level is illustrated in Figure 1.1.

Applying the NS solver to compute (1) the initial solution, (2) the residual on the h/2-mesh and (3) the corrected solution ensure that the corrector encompasses all the features of the solver (extrapolation, Riemann solver, limiters etc.). The number of iterations required to compute the corrected solution depends on the considered problem (essentially its nonlinearity), the simplified Jacobian and the linear solver used (preconditioner and fixed point iteration) in case of incomplete linear convergence. As the initial and the corrected solution are relatively close to each other, the number of iterations required to compute the corrected solution is typically about 1/10 to 1/100 of the number of iterations required to compute the initial solution. Finer subdivision (h/4, h/8, etc.) of the initial mesh can also be used to improve the predictions. It is very important to note that linearly interpolating the discrete solution W_h on the h/2-mesh (i.e. using Π_h/2(W_h)) is nonetheless convenient but also essential to the precision of the correction. Replacing Π_h/2(W_h) by a better polynomial reconstruction tends to reduce the efficiency of the correction. Indeed, as we want to quantify the error of piecewise linear solution W_h with respect to the h/2-mesh; it is counter productive to try to improve Π_h/2(W_h). More numerical implementation details are given in Frazza et al. (2019).

1.4.2. Finite element corrector

We transpose the variational derivation of section 1.2.2:

In particular, the exact solution verifies

but generally

We have

and by linearizing the fluxes in Π_hW and W_h, we have

and

From these three relations above, we deduce for all φ_h ∈ V_h:

[1.20]

[1.21]

This expression relies on the fact that W_h is relatively close to Π_hW, so that the derivatives are approximatively the same. It has two major advantages. First, it does not require to compute the derivative , which saves a lot of trouble: as we mentioned previously, in linear approaches, the quality of the correction depends on the choice and the precision of these derivatives. Second, it automatically takes into well of the residual. As long as the higher derivatives the nonlinearitiesare relatively close, which is a reasonable assumption as the correction is a purely local correction, the Taylor developments mentioned previously can be expanded.

In relation [1.21], we still need to estimate W − Π_hW that represents a loss of smoothness due to the discretization of the solution; it does not depend on the absolute value of the solution but on its variations. The term W+ W − Π_hW can be thus estimated by reconstructing a smooth solution (Clément 1975; Zienkiewicz and Zhu 1992; Bank and Smith 1993) from the actual numerical solution W_h as shown in Figure 1.2. Finally, the source term can be expressed as

[1.22]

where U^k+1 holds for the (k +1)^th order reconstruction if the solution is of order k. We propose Algorithm 1.2. to compute the corrector. Note that there is no need for a finer mesh in this case as the integrals are computed on the same elements, but more quadrature points are required to have a consistent computation of the source term given by relation [1.22].

REMARK 1.1.– For both nonlinear corrector methods, we take great care to propose methods that do not require any additional memory costs. With regard to the computational cost, the computation of the source term S_h is approximatively the cost of one flow solver iteration, which is negligible.