R Duvigneau

CAD-consistent adaptive refinement using a NURBS-based Discontinuous Galerkin method

This study concerns the development of a new method combining high-order CAD-consistent grids and adaptive refinement / coarsening strategies for efficient analysis of compressible flows. The proposed approach allows to use geometrical data from Computer-Aided Design (CAD) without any approximation. Thus, the simulations are based on the exact geometry, even for the coarsest discretizations. Combining this property with a local refinement method allows to start computations using very coarse grids and then rely on dynamic adaption to construct suitable computational domains. The resulting approach facilitates interactions between CAD and Computational Fluid Dynamics (CFD) solvers and focuses the computational effort on the capture of physical phenomena, since geometry is exactly taken into account. The proposed methodology is based on a Discontinuous Galerkin (DG) method for compressible Navier-Stokes equations, modified to use Non-Uniform Rational B-Spline (NURBS) representations. Local refinement and coarsening are introduced using intrinsic properties of NURBS associated to a local error indicator. A verification of the accuracy of the method is achieved and a set of applications are presented, ranging from viscous subsonic to inviscid trans-and supersonic flow problems.

, mainly concerning incompressible viscous flows. Indeed, most CFD methods for convection-dominant problems rely on FV or DG formulations, which are more suited to hyperbolic conservation laws.

.

. The demonstrated results are convincing, in terms of accuracy and convergence, but the proposed methodology may suffer from the complexity of the underlying mappings used for spatial integration.

To fully integrate CAD data in a compressible flow solver environnment, we developed in [29] a CAD-consistent DG method, able to preserve exactly CAD geometries defined using NURBS. The high-accuracy obtained was demonstrated on some academic test-cases. However, the application range was limited, because of the difficulty to generate NURBSbased grids, which are fully adapted to both geometry and solution. Therefore, an extension is presented in this work, which relies on adaptive refinement and coarsening techniques in the context of NURBS-based grids. With this adaptive approach, we obviously aim at improving the quality of the solution and the computational efficiency of the solver, but also at facilitating the initial grid generation. Indeed, we intend to intensively use local refinement procedures operating on a very coarse initial curved mesh, whose generation should be far easier than usually. The article is organized as follows: the underlying NURBS representations are described in section 2, while grid generation and refinement/corsening techniques are presented in section 3. The CAD-consistent DG method is described in section 4, followed by the adaptive strategy in section 5. Finally, a verification exercise is analysed in section 6 and some more complex applications are presented in section 7, ranging from a subsonic viscous flow to inviscid trans-and supersonic flows.

Introduction

Computation Fluid Dynamics (CFD) is now commonly employed in engineering design for several applications. This evolution is due to the increase of computational facilities, which has significantly reduced the computational burden related to CFD, and also to the maturity of underlying numerical methods. Nevertheless, the time necessary to set up CFD simulations is still significant because of the complexity of pre-processing tasks, like mesh generation from Computer-Aided Design (CAD) data. Despite the progress of CAD and grid generation software, the automation of geometrical processing is limited and, as consequence, time-consuming human work is still necessary.

Besides, several recent studies have pointed out the limitations related to the use of piecewise linear grids and the necessity to use curved meshes to obtain the theoretical convergence rates. The use of a piecewise linear approximation of the geometry may also lead to spurious phenomena in the solution [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible navier-stokes equations[END_REF][START_REF] Silveira | Higher-order surface treatment for discontinuous Galerkin methods with applications to aerodynamics[END_REF][START_REF] Costa | High-order accurate finite volume scheme on curved boundaries for the two-dimensional steady-state convection-diffusion equation with dirichlet condition[END_REF]. This is especially true in the context of high-order schemes based on Discontinuous Galerkin (DG) or Finite-Volume (FV) methods, which are more and more employed to solve complex problems, for instance in Large Eddy Simulation (LES). As shown by some authors, the use of a high-order geometry description is beneficial, not only in terms of solution accuracy, but also in terms of computational efficiency, because it avoids performing excessive mesh refinement to describe curved boundaries [START_REF] Sevilla | NURBS-enhanced finite element method for euler equations[END_REF].

Therefore, the objective of the current work is to adress both issues: facilitate the transfer of geometrical data from CAD to analysis and use curved grids to maximize the accuracy and efficiency of CFD solvers. A similar idea is adopted in the so-called isogeometric analysis approach proposed ten years ago in the context of Finite-Element (FE)

NURBS representations

NURBS curves (and surfaces) are now considered as standard to define geometries in CAD [START_REF] Farin | Curves and Surfaces for Computer-Aided Geometric Design[END_REF]. In particular, they allow to represent exactly a broad class of geometric curves commonly encountered in engineering, like conic sections. Moreover, they permit a very intuitive definition and modification of geometrical objects through the handling of control points. In the following sections, we describe the most relevant features of such representations and explain how to construct a computational grid based on NURBS, yielding a domain geometrically exact with respect to CAD.

Basis functions

NURBS functions [START_REF] Piegl | The NURBS book[END_REF] are rational extensions of B-Spline functions [START_REF] Boor | A Practical Guide to Splines[END_REF]. Both are defined using a so-called knot vector Ξ = (ξ 1 , . . . , ξ l ) ∈ R l , which consists of nondecreasing l real numbers. This knot vector defines a discretization of the parametric domain [ξ 1 , ξ l ]. Open knot vectors, i.e. knot vectors with first and last knots of multiplicity p + 1, are usually used for representations of degree p to impose interpolation and tangency conditions at both extremities [START_REF] Boor | A Practical Guide to Splines[END_REF].

Therefore ξ 1 = . . . = ξ p+1 and ξ n+1 = . . . = ξ n+p+1 . B-Spline basis functions (N p i ) i=1,••• ,n are defined recursively as [START_REF] Boor | A Practical Guide to Splines[END_REF]:

N 0 i (ξ) = 1 if ξ i ≤ ξ < ξ i+1 0 otherwise (1)
N p i (ξ) = ξ -ξ i ξ i+p -ξ i N p-1 i (ξ) + ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 N p-1 i+1 (ξ). (2)
Note that the quotient 0/0 is assumed to be zero. According to this definition, B-Spline functions of degree 0 and 1 coincide with classical piecewise constant and linear FE basis functions. However, B-Spline functions of higher degree differ from classical Lagrange polynomials. In particular, they exhibit a higher regularity [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF]. We underline that this recursive definition is practically convenient to construct numerical schemes of arbitrary order and promotes the use of p-adaptive strategies. The degree of the functions p, the number of knots l and functions n are related by the relation l = n + p + 1 [START_REF] Piegl | The NURBS book[END_REF]. A set of such basis functions is illustrated in Fig. (2a).

Multi-dimensional NURBS basis functions of degree p are defined by using tensor products of B-Spline functions and by associating a weight to each function. We present below the two-dimensional case in detail, but the cases corresponding to other dimensions can be derived in a straightforward way. Thus, the two-dimensional NURBS basis function of index i 1 × i 2 is defined in the parameteric domain [ξ 1 , ξ l ] × [η 1 , η l ] (for the sake of simplicity we assume the number of knots and the degrees equal for the two parameters) by associating the weight w i1 i2 to the tensor product of B-Spline functions according to:

R p i1 i2 (ξ, η) = w i1 i2 N p i1 (ξ) N p i2 (η) n1 j1=1 n2 j2=1 w j1 j2 N p j1 (ξ)N p j2 (η) . (3)

NURBS patches

Once the NURBS basis is defined, planar NURBS surfaces, that will be considered then as building blocks of the computational domain and denoted as patches, can be easily defined by associating to the functions a set of geometrical degrees of freedom:

x(ξ, η) = n1 i1=1 n2 i2=1 R p i1 i2 (ξ, η)X i1i2 , (4)
where x = (x, y) are the cartesian coordinate values spanned by the patch. In this context, the coefficients associated to the basis functions (X i1i2 ) i1=1,••• ,n1 i2=1,••• ,n2 are denoted as control points. In the following, we simplify this expression by omitting the double index notation and the degree, which gives:

x(ξ, η) = n1×n2 i=1 R i (ξ, η)X i . ( 5
)
This construction is illustrated in Fig. [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible navier-stokes equations[END_REF], for a cubic patch (p = 3) with a 4×4 control point lattice (n 1 = n 2 = 4), in the particular case ξ 1 = ξ 2 = ξ 3 = ξ 4 = 0 and ξ 5 = ξ 6 = ξ 7 = ξ 8 = 1 (same knot vector for η). An important property is the fact that any NURBS curve or surface can be considered as the projection of a B-Spline curve or surface defined in a space of higher dimension using its weights [START_REF] Boor | A Practical Guide to Splines[END_REF]. For instance, a NURBS curve lying in the plane (x, y), characterized by a set of control points

(X i ) i=1,••• ,n = (x i , y i ) i=1,••• ,n and weights (ω i ) i=1,••• ,n
, can be described as the projection of a B-Spline curve lying in a 3D space and defined by the control points (x i ω i , y i ω i , ω i ) i=1,••• ,n . This representation of NURBS objects as B-Spline ones can be useful for practical manipulations.

ξ η 0 1 1 Ω • • • • • • • • • • • • • • • • X i1i2 Ω F x y Figure 1: NURBS representation
The transformation of the parametric domain Ω = [ξ 1 , ξ l ]×[η 1 , η l ] to the physical domain Ω can now be introduced:

F : Ω → Ω, (ξ, η) → x(ξ, η) = n1×n2 i=1 R i (ξ, η)X i . ( 6
)
This transformation is in general nonlinear. We make the asumption that it is injective and defines a map between the physical domain and the parametric one. This type of representation was used originally by T. Hughes and co-authors for the foundation of the isogeometric analysis paradigm, which is in short a Finite-Element (FE) approach relying on NURBS representations. The knot intervals in the parametric domain are considered as reference elements used in a Galerkin method. More sophisticated representations, such as T-Splines [START_REF] Bazilevs | Isogeometric analysis using T-splines[END_REF] and LR-Splines [START_REF] Dokken | Locally refined splines[END_REF] were proposed latter to achieve efficient local refinement. Nevertheless, in the context of a DG formulation, NURBS patches cannot be employed directly and an additional ingredient must be introduced to generate discontinuous solutions at element interfaces. This topic is addressed in section 3.

Knot insertion and Bézier extraction

An important property of NURBS representations is the capability to insert a new knot, and thus a new basis function, without altering the geometry [START_REF] Boor | A Practical Guide to Splines[END_REF]. It can be considered as a local h-refinement procedure [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF]. The one-dimensional case is described below, the extension to higher dimension is straightforward thanks to the tensor product structure.

Any NURBS curve defined by the knot vector Ξ = (ξ 1 , . . . , ξ l ) ∈ R l , n basis functions and control points, as described above, can be identically represented using the knot vector (ξ 1 , . . . , ξ q , ξ, ξ q+1 . . . , ξ l ) ∈ R l+1 , that includes an additional knot ξ inserted between ξ q and ξ q+1 . If we consider first the case of a B-Spline curve, the new representation can be written as:

x(ξ) = n i=1 N p i (ξ)X i = n+1 i=1 N p i (ξ) Xi , (7)
with a new set of control points defined as:

Xi = (1 -α i )X i-1 + α i X i α i =    1 if i ≤ q -p ξ-ξi ξi+p-ξi if q -p + 1 ≤ i ≤ q 0 if i ≥ q + 1 (8)
The extension to NURBS curves relies on the projection property already mentioned [START_REF] Piegl | The NURBS book[END_REF]. The extension for 2D or 3D patches simply consists in applying the above algorithm for each direction independently.

It is important to underline now a particular case: if a new knot is inserted at an existing knot location, the regularity of the curve is decreased. More generally, the curve at the knot ξ q has the regularity C p-r , where r is the multiplicity of the knot q [START_REF] Boor | A Practical Guide to Splines[END_REF]. Therefore, if one inserts p knots at an existing knot location, the geometry of the curve is preserved but the curve is now divided in two independent parts. If this multiple-knot insertion is achieved for all inner knots, the curve is decomposed into a set of independent Bézier curves. This procedure is named Bézier extraction and will be used to define an approximation space suitable to DG formulation, as explained below.

3 Generation and refinement of CAD-consistent domains

Initial NURBS-based domain

We consider now an arbitrary physical domain Ω ⊂ R 2 with the boundary ∂Ω, that will be used as computational domain. We make the asumption that the geometry of this boundary is defined thanks to NURBS representations and, consequently, that we can define a set of NURBS patches that span the physical domain exactly. Note that the practical construction of such patches is far from being straightforward. It is analogous to building structured meshes, or multi-block structured meshes, but in the context of NURBS representations. In particular, specific techniques should be employed to ensure injectivity [START_REF] Xu | Parametrization of computational domain in isogeometric analysis: methods and comparison[END_REF][START_REF] Xu | Constructing analysis-suitable parameterization of computational domain from cad boundary by variational harmonic method[END_REF]. The description of general methods to construct such domains is beyond the scope of the current study, but the reader can refer to [START_REF] Xu | Parametrization of computational domain in isogeometric analysis: methods and comparison[END_REF][START_REF] Xia | Continuity and convergence in rational triangular bézier spline based isogeometric analysis[END_REF][START_REF] Engvall | Isogeometric triangular bernstein-bézier discretizations: automatic mesh generation and geometrically exact finite-element analysis[END_REF][START_REF] Engvall | Isogeometric unstructured tetrahedral and mixed-element bernstein-bézier discretizations[END_REF][START_REF] Xu | Analysis-suitable volume parameterization of multi-block computational domain in isogeometric analysis[END_REF] to have some examples of admissible methods. The construction of a NURBS-based computational domain is certainly the main issue that isogeometric analysis faces currently. However, the approach proposed in this study aims at alleviating the difficulty by intensive use of local refinement techniques, starting from an extremely coarse initial computational domain, composed of a few elements only.

As explained above, we suppose therfore that the computational domain can be defined by a set of NURBS patches, that exactly matches the boundary of the physical domain. However, the underlying NURBS functions are not suitable to DG methods, since each support of function covers several knot intervals, as shown in Fig. (2a), which does not allow to define solution jumps at the interface between elements. Therefore, the patches have to be transformed into a set of fully disconnected elements, and associated basis functions, which can exhibit discontinuities at the interfaces. Obviously, this should be done without altering the geometry. Hopefully, this task can be easily achieved thanks to the Bézier extraction procedure, described in section 2.3. Indeed, as explained above, if p knots are inserted at existing inner knots of a patch of size n 1 × n 2 , the original NURBS patch is divided to a set of (n 1 -p) × (n 2 -p) rational Bézier elements of degree p. The geometry of the physical domain is unchanged, but each element is now defined according to its own basis of size (p + 1) × (p + 1), which enables the generation of discontinuities at the interfaces.

This procedure is illustrated in Fig. (2) for a one-dimensional case. Five NURBS functions of degree p = 2, which are non-zero on three selected knot intervals, are plotted in Fig. (2a). Their support is composed of p + 1 = 3 knot intervals. The knot insertion procedure applied twice for each existing inner knot yields nine quadratic rational Bézier functions, each of them with a support restricted to one knot interval, as illustrated by Fig.

(2b).

A two-dimensional illustration is provided by Fig. [START_REF] Costa | High-order accurate finite volume scheme on curved boundaries for the two-dimensional steady-state convection-diffusion equation with dirichlet condition[END_REF]. One considers the domain delimited by two quarters of cylinder, which can be described as a NURBS patch of degree p = 2 using a lattice of 4 × 4 control points, as depicted in Fig. (3a). The knot vector for this example is [0 0 0 1 2 2 2]. The basis functions span the whole parametric intervals, which prevent from generating discontinuities at interface ξ 4 (and η 4 ). As explained above, by inserting two additional knots at inner knot ξ 4 (and η 4 ), the domain is splitted in four rational Bézier elements, each of them being defined according to its own basis functions (see Fig. (3b)), yielding a DG-compliant computational domain whose geometry matches exactly the original CAD representation. Note that this Bézier extraction technique can also be useful for continuous elements, see for instance [START_REF] Hennig | Bézier extraction and adaptive refinement of truncated hierarchical NURBS[END_REF][START_REF] D'angella | Multi-level Bézier extraction for hierarchical local refinement of isogeometric analysis[END_REF]. As summary, the proposed approach to generate an initial CAD-consistent DG-compliant domain is composed of two steps: firstly, generate a set of NURBS patches matching CAD boundaries; secondly, apply Bézier extraction to transform the patches into DG-compliant rational Bézier elements. This approach will be illustrated on several test-cases in the application section.

ξ ξ 1 ξ 2 ξ 3 ξ 4 (

Local refinement

In the original isogeometric analysis paradigm proposed by T. Hughes, local refinement necessitated to introduce sophisticated approaches to avoid the propagation of refinement due to the underlying tensorial representation of NURBS patches, while maintaining a high-order regularity [START_REF] Cottrell | Studies of refinement and continuity in isogeometric analysis[END_REF]. In the present approach, local refinement is much simpler because the computational domain is composed of rational Bézier elements which are already disconnected. Therefore, it just consists in inserting p + 1 new knots in a given element to split it into two new ones. The price to pay is obviously the loss of regularity at the new interface, which is however natural in a DG framework. The refinement of each element can be achieved in an isotropic or anisotropic manner, as illustrated by Fig. [START_REF] Sevilla | NURBS-enhanced finite element method for euler equations[END_REF], by inserting the knots in one direction ξ or η only, or in both directions simultaneously. The procedure is illustrated for a real computational domain in Fig. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF], in the case the refinement is ruled explicitly by the user. The domain is initially defined by a single cubic NURBS patch, according to the CAD geometry of an airfoil. This domain is splitted into 8 × 6 rational Bézier elements (see Fig. (5a)), yielding a first very coarse DG-compliant grid, that matches exactly the airfoil. One refines then the area located around the airfoil by applying successive refinement in user-defined boxes (see Fig. (5b)). This approach generates a suitable computational domain that preserves the initial CAD geometry, but it supposes to know a priori where refinement should be achieved. This is not an easy task in general, especially for complex time-dependent problems. Therefore a fully automated local adaption procedure is developed, in the framework of the DG discretization described below.

Discontinuous Galerkin method

We describe in this section how a standard DG method can be modified to account for the proposed rational Bézier elements. The discretization is first detailed. Then, the issue of shock capturing is considered in the context of these

• • • • • • • • • • • • • • • • ξ ξ 1,2,3 ξ 4 ξ 5,6,7 (a) One NURBS patch • • • • • • • • • • • • • • • • • • • • • • • • • ξ ξ 1,2,3 ξ 4,5,6 ξ 7,8,9
(b) Four rational Bézier elements

Flow model

We consider the two-dimensional compressible Navier-Stokes equations, that can be written in the conservative form as follows:

∂W ∂t + ∇ • ( F -G) = 0, (9)
where W represents the conservative flow variables (ρ, ρu, ρv, ρe), with ρ the density, U = (u, v) the velocity vector and e the total energy per unit of mass. We denote c the sound speed. F = (F x (W), F y (W)) is the vector of the inviscid fluxes, whereas G = (G x (W, S), G y (W, S)) is the vector of the viscous fluxes, with S = ∇W the gradient fields. The pressure p is obtained from the perfect gas state equation:

p = ρ(γ -1)(e - 1 2 - → U 2 ) = ρ(γ -1)e i , (10)
where γ is the ratio of the specific heat coefficients and e i the internal energy. The inviscid fluxes are given by:

F x (W) =     ρu ρu 2 + p ρuv ρu(e + p ρ )     F y (W) =     ρv ρvu ρv 2 + p ρv(e + p ρ )     . (11)
The viscous fluxes write as:

G x (W) =     0 τ xx τ yx uτ xx + vτ yx -q x     G y (W) =     0 τ xy τ yy uτ xy + vτ yy -q y     , (12)
where the symmetric viscous stress tensor τ for Newtonian fluids is:

τ xx = 2 3 µ 2 ∂u ∂x - ∂v ∂y τ xy = µ ∂u ∂y + ∂v ∂x τ yy = 2 3 µ 2 ∂v ∂y - ∂u ∂x . ( 13
)
The heat flux q is based on Fourier's law:

(q x , q y ) = -λ µ P r ∂e i ∂x , ∂e i ∂y , (14)
where P r is the Prandtl number.

Discretization

To apply a DG method [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF], one first writes a weak formulation of the problem by multiplying Eq. ( 9) by an arbitrary test function ϕ and by integrating over an element Ω j . In the context of the CAD-consistent computational domain presented above, Ω j is defined as a rational Bézier element. Thus, the weak formulation writes:

Ωj ∂W ∂t ϕ dΩ + Ωj ∇ • ( F (W) -G(W, S)) ϕ dΩ = 0. (15)
After integration by parts, one obtains classically:

Ωj ∂W ∂t ϕ dΩ - Ωj ( F (W) -G(W, S)) • ∇ϕ dΩ + ∂Ωj ( F (W) -G(W, S)) • n ϕ dΓ = 0. ( 16
)
Since the solution is a priori discontinuous at the interfaces, the normal fluxes are evaluated via numerical flux functions F n (W + , W -) and G n (W + , W -, S + , S -), defined according to the values of the solution that prevail at each side of the interface. Several flux functions are classically used in DG methods [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF][START_REF] Cockburn | High-Order Methods for Computational Physics, chap. Discontinuous Galerkin Methods for Convection-Dominated Problems[END_REF]. In this work, the HLLC flux is employed [START_REF] Batten | On the choice of wavespeeds for the hllc Riemann solver[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] for the inviscid terms, whereas the LDG flux [START_REF] Cockburn | The local discontinuous galerkin method for time-dependent convection-diffusion systems[END_REF] is adopted for the viscous terms. The latter requires the evaluation of the gradient fields:

S = ∇W = ( ∂ρ ∂x ,
By multiplying this definition by a test function and integrating by parts, one obtains a DG formulation for the gradient fields:

Ωj S ϕ dΩ - Ωj H(W) • ∇ϕ dΩ + ∂Ωj H(W) • n ϕ dΓ = 0, (18)
where H = (H x (W), H y (W)) is the vector of the gradient fluxes, expressed as:

H x (W) =             ρ 0 ρu 0 ρv 0 ρe 0             H y (W) =             0 ρ 0 ρu 0 ρv 0 ρe             . ( 19
)
One introduces H n (W + , W -), the corresponding numerical flux computed at element interfaces according to the LDG method.

Finally, one has to choose a basis for the representation of the solution on each element Ω j , denoted W| j . Any basis defined on Ω j could be selected, but it has been found convenient to employ again the rational Bézier basis of degree p used to describe the geometry of the element. This choice was already tested in [START_REF] Duvigneau | Isogeometric analysis for compressible flows using a Discontinuous Galerkin method[END_REF] with satisfactory results. It was compared to a possible alternate choice consisting of Lagrange polynomials to define the solution associated to a NURBS-based geometry [START_REF] Michoski | Foundations of the blended isogeometric discontinuous Galerkin (BIDG) method[END_REF]. The accuracy of the two approaches was found very similar. Therefore, for the sake of simplicity, we choose to use a unique basis for both geometry and solution fields. In this sense, the proposed approach can be qualified as isogeometric because the same CAD-consistent representation is employed for the geometry and the solution. Therefore, conservative variables are expressed as:

W| j (ξ, η) = (p+1) 2 i=1 R i (ξ, η)W i . ( 20
)
The same basis is also employed for the test functions. Finally, by substitution and expressing integrals in the parametric domain, one obtains:

(p+1) 2 i=1 ∂W i ∂t Ωj R i R k |J Ω | d Ω = Ωj ( F (W) -G(W, S)) • ∇R k |J Ω | d Ω - ∂ Ωj (F n (W + , W -) -G n (W + , W -, S + , S -)) R k |J Γ | d Γ ( 21
) (p+1) 2 i=1 S i Ωj R i R k |J Ω | d Ω = Ωj H(W) • ∇R k |J Ω | d Ω - ∂ Ωj H n (W + , W -) R k |J Γ | d Γ. ( 22
)
|J Ω | represents the determinant of the Jacobian matrix of the geometric transformation F, defined according to Eq. ( 6) and |J Γ | is the lineal counterpart. Obviously, one identifies easily a mass matrix on the left-hand side of Eqs. [START_REF] Aigner | Swept volume parameterization for isogeometric analysis[END_REF][START_REF] Engvall | Isogeometric triangular bernstein-bézier discretizations: automatic mesh generation and geometrically exact finite-element analysis[END_REF], volumic residuals and interface fluxes on the right-hand side. The local mass matrix is inverted once and its inverse is stored, an explicit third-order four-step Strong Stability Preserving (SSP) Runge-Kutta method being used for time integration. We underline that the spatial integrals in Eq. [START_REF] Aigner | Swept volume parameterization for isogeometric analysis[END_REF][START_REF] Engvall | Isogeometric triangular bernstein-bézier discretizations: automatic mesh generation and geometrically exact finite-element analysis[END_REF] have to be evaluated carefully. Indeed, the nonlinearities in the flow variables and in the geometrical transformation necessitate high-order quadratures, depending on the basis degree. For the present work, classical Gauss-Legendre quadrature rules have been applied, but involving a large number of evaluation points when required. On most computations presented below, the use of p + 1 quadrature points per direction has been found satisfactory. However, if highly distorded elements are observed in the mesh, p + 2 quadrature points may be required to reach asymptotic convergence rate.

Initial and boundary conditions

For unsteady problems, the initial flow solution has to be expressed accordingly to the approximation space selected, i.e. in terms of discontinuous rational Bézier functions, as in Eq. ( 20). This task is achieved by solving a set of local least-squares fitting problems. If the initial solution is defined globally by w 0 , the initial solution in each element Ω j is obtained by:

W| j t=0 = argmin 1 2 Ωj   (p+1) 2 i=1 R i W i -w 0   2 dΩ. (23)
The optimality condition with respect to the degree of freedom W k , expressed in the parametric domain, yields a linear system involving the local mass matrix:

(p+1) 2 i=1 W i Ωj R i R k |J Ω | d Ω = Ωj w 0 R k |J Ω | d Ω. (24)
Boundary conditions are imposed weakly via the normal fluxes at boundary, which can be computed using different approaches [START_REF] Mengaldo | A guide to the implementation of boundary conditions in compact high-order methods for compressible aerodynamics[END_REF], either physical or numerical. For inlet and outlet boundaries, a numerical flux is evaluated using the interior state on one side and a reconstructed exterior state on the other side, computed from the free-stream conditions using Riemann invariants [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF][START_REF] Mengaldo | A guide to the implementation of boundary conditions in compact high-order methods for compressible aerodynamics[END_REF]. For wall boundaries, the physical flux is directly imposed, involving only the pressure terms for the inviscid part and the shear stress and heat flux for the viscous part.

Shock capturing

It is necessary to stabilize oscillations appearing in the solution in presence of shocks, either by introducing a solution limitation procedure or by adding artifical viscosity. In this work, the sub-cell capturing approach proposed by Persson and Peraire [START_REF] Persson | Sub-cell shock capturing for discontinuous galerkin methods[END_REF] is extended to rational Bézier representations.

An artifical local viscosity µ j , considered as constant over each element Ω j , is introduced in the model. As shown in [START_REF] Barter | Shock capturing with PDE-based artificial viscosity for DGFEM: Part i. formulation[END_REF], a more regular viscosity field could be beneficial and will be envisaged in a future work. In the case of Navier-Stokes equations, this local viscosity is just added to the physical diffusion coefficient µ. In the case of Euler equations, it necessitates to add the viscous terms in the discretization. As shown in [START_REF] Persson | Sub-cell shock capturing for discontinuous galerkin methods[END_REF], the use of the physical viscous flux offers some advantages with respect to a more simple laplacian-based smoothing operator. More precisely, the artifical viscosity µ j is defined by the following smooth function:

µ j =      0 if s j < s 0 -κ µ 0 j 2 1 + sin π(sj -s0) 2κ if s 0 -κ < s j < s 0 + κ µ 0 j if s j > s 0 + κ (25)
where µ 0 j is the maximum artificial viscosity admissible in the element Ω j , s j a sensor of the oscillations of the solution in the element and (s 0 , κ) global numerical parameters that control the artificial viscosity distribution. As analysed in [START_REF] Persson | Sub-cell shock capturing for discontinuous galerkin methods[END_REF], µ 0 should be defined according to local length and velocity scales. To capture the shock in a single cell, the length scale should be h/k, where h is a measure of the element length and k the scheme order. In the context of the proposed rational Bézier representation, the min distance between element corners is chosen as element length, i.e. h = min(|X 11 -

X 1n2 |, |X 11 -X n11 |, |X n11 -X n1n2 |, |X 1n2 -X n1n2 |)
, which is easy to estimate whatever the element shape. Regarding the velocity scale, the free-stream velocity u ∞ is chosen, yielding a maximum viscosity defined locally as:

µ 0 j = h p + 1 u ∞ (26)
The parameter s 0 is scaled according to a Fourrier analysis [START_REF] Persson | Sub-cell shock capturing for discontinuous galerkin methods[END_REF]: if the solution is smooth, the amplitude of highfrequencies should decay like O(1/k 2 ). Therefore, the switch s 0 is set as s0 /(p + 1) 2 , with s0 a problem dependent parameter. Its value controls the activation of the artificial viscosity. For some applications, it has to be tuned precisely to avoid affecting local extrema in the solution, as will be shown in section 7. κ controls the extent of the transition and is chosen equal to s0 /2. It remains to define the local sensor s j . In the original method [START_REF] Persson | Sub-cell shock capturing for discontinuous galerkin methods[END_REF], the sensor is based on the L 2 -norm of the highest term in a modal expansion of the solution. In this work, we propose to define the sensor by quantifying the oscillations in the control point lattice that represents the solution in each element. More precisely, one evaluates the following quantities, that compare the total variation of control points in each direction to the variation observed between element extremities (we suppose here that p > 1):

s i2 ξ = p i1=1 |W i1+1 i2 -W i1i2 | -|W p+1 i2 -W 1i2 | (27)
s i1 η = p i2=1 |W i1i2+1 -W i1i2 | -|W i1 p+1 -W i11 | (28)
These quantities vanish if the solution is monotonic over the element. If not, it is a measure of the oscillatory nature of the solution in each parametric direction, thanks to the convex hull property of NURBS curves [START_REF] Piegl | The NURBS book[END_REF]. The sensor for a given element is finally defined as:

s j = 1 2(p + 1) p+1 i2=1 s i2 ξ + p+1 i1=1 s i1 η 1 W ( 29
)
where W is the mean solution over the element. In the illustration presented below, this sensor is applied to the density field, but other choices could be envisaged, such as Mach number for instance.

Adaptive refinement

As explained in section 3.2, a procedure based on successive user-defined refinements is very limited in terms of applications, because it relies on the intuition of the user. A common approach in engineering is to refine the domain in the vicinity of walls for an accurate approximation of the geometry and boundary layers, in the wakes to capture vortex shedding, possibly in regions where shocks could appear, etc. Obviously, this becomes tedious when unsteady flows are considered. This task is even harder when using high-order schemes, because the required resolution level is less intuitive. The aim of the present work is to make adaption easier and more efficient. Accounting exactly for the CAD geometry, whatever the discretization level used for the resolution, is a first ingredient (exact-geometry property). It avoids refining locally the grid only to approximate the geometry and, thus, allows to focus refinement on the accuracy of the physical solution. The second ingredient is a fully automated, time-dependent, refinement and coarsening procedure, excluding user intuition to select areas to be refined. The joint use of both ingredients allows computations starting with an extremely coarse CAD-compliant domain, which is automatically and continuously adapted to the physical solution along time. The refinement procedure for the geometry has already been detailled in section 3.2. In the following sections, we describe the dynamic implementation of the computational domain, the refinement and coarsening criteria and the treatment of the solution.

Adaptive domain management

On the basis of the refinement procedure for rational Bézier elements described in section 3.2, a dynamic grid management is adopted using a tree-based approach. The root is composed of elements of level 0, defined by the user from CAD data (initial grid), as described in section 3.1. Then, the refinement of any element of level l yields four elements of level l + 1, an isotropic splitting being chosen for the sake of simplicity. These four elements are denoted as brothers. Any element of level l has the knowledge of its corresponding element at level l -1, denoted as father, and corresponding elements at level l + 1, denoted as sons, if they exist. Thus, the grid structure can be refined arbitrarily, with the constraint that two neighboring elements cannot have a difference of levels higher than one. This constraint is added to prevent an element from being faced to more than two refined elements in any direction, for mesh regularity reasons. A maximum refinement level l max is imposed to avoid an infinite series of refinements, that may occur in case of shock adaption for instance. Obviously, all existing elements are not used at each time step for the resolution. Therefore, elements are flagged as active if they are currently used for integration, or inactive. A very similar management is adopted for faces, which are stored using a tree-based approach. The refinement of a given face yields two faces, denoted as daughters, whereas the original face is denoted as mother. Faces can also be flagged as active or inactive.

The use of such a locally refined computational domain has a very limited impact on the implementation of the DG method, thanks to the flexibility of the DG scheme to handle non-conforming grids. The only one modification concerns the flux integration at the interface between two elements of different levels, as illustrated in Fig. [START_REF] Bazilevs | Isogeometric analysis: approximation, stability and error estimates for refined meshes[END_REF]. In this case, the numerical quadrature is achieved by using the integration points of the smallest element. Thus, local refinement has no impact on the scheme accuracy, as will be verified later. Obviously, the parametric representation of two neighboring elements could be different, which should be tacken into account for the evaluation of left and right quantities.

x y

• • • η = 0 η = 0.5 η = 1 η = 0 η = 1

Refinement and coarsening criteria

Several refinement criteria can be found in the literature to drive the adaption procedure in the context of DG discretizations, see for instance [START_REF] Leicht | Anisotropic mesh refinement for discontinuous galerkin methods in two-dimensional aerodynamic flow simulations[END_REF][START_REF] Papoutsakis | An efficient adaptive mesh refinement (AMR) algorithm for the discontinuous galerkin method: Applications for the computation of compressible two-phase flows[END_REF][START_REF] Sun | Anisotropic and dynamic mesh adaptation for discontinuous galerkin methods applied to reactive transport[END_REF]. In the present study, we choose an approach based on the measure of interface jumps, able to tackle a broad range of problems easily. More sophisticated approaches could be used as well and will be studied in forthcoming works. Nevertheless, as will be shown in application section, this indicator has the ability to identify regions exhibiting both discontinuities and under-resolved elements in case of smooth solutions. In practice, for each element Ω j , one evaluates the solution jump at the interface between Ω j and the neighboring elements:

j = k∈Nj Γ jk W| j -W| k dΓ, (30)
where N j represents the set of elements around Ω j and Γ jk the interface between Ω j and Ω k . Then, the element Ω j is flagged for refinement if the indicator j exceeds the indicator average ¯ by a user-defined ratio ref :

j ¯ > ref . (31)
Similarly, the element is flagged for coarsening if j is lower than ¯ , with a user-defined ratio coa :

j ¯ < coa . (32)
From geometrical point of view, the refinement of a selected element is achieved using the knot insertion procedure described in section 2.3. The coarsening can be performed only if all the four brother elements are flagged for coarsening. In this case, no new computation is required, it just consists in changing their status as inactive, while their father becomes active again. Faces are also updated accordingly.

Solution refinement and coarsening

When the computational domain is refined or coarsened, the current solution has to be updated. In the case of refinement, we benefit from the fact that geometry and solution have the same rational Bézier representation. Therefore, both are splitted simultaneously and identically using the same knot insertion procedure described above. As consequence, no loss of information occurs during refinement.

In the case of coarsening, the solution should be approximated in the father element, from data located in the four sons elements. This is achieved by using a least-squares projection. Similarly as the approximation of the initial condition, see Eqs [START_REF] Engvall | Isogeometric unstructured tetrahedral and mixed-element bernstein-bézier discretizations[END_REF][START_REF] Jaxon | isogeometric analysis on triangulations[END_REF], one first computes a solution estimate W| f LS in the father element Ω f from the solutions (W| s ) s=1,...,4 available in the son elements (Ω s ) s=1,...,4 by solving the following local linear system:

(p+1) 2 i=1 W i | f LS Ω f R i R k |J Ω | d Ω = 4 s=1 Ωs W| s R k |J Ω | d Ω. (33)
Then, this approximation is shifted to enforce conservation, yielding a solution in the father element defined as:

W| f = W| f LS + 1 A f 4 s=1 Ωs W| s |J Ω | d Ω - Ω f W| f LS |J Ω | d Ω , ( 34
)
where A is the area of a given element. The integration of the above equation yields the conservation of the mean solution during coarsening, up to quadrature accuracy, which is a necessary condition. Contrary to the refinement process, coarsening leads however to a loss of information during the least-squares fitting procedure.

6 Verification

Steady problem with curved geometry

We consider firstly the inviscid flow between two quarters of cylinders, which allows to verify the scheme accuracy in presence of curved wall boundaries. The geometry of the computational domain is exactly defined by a single NURBS patch, which is then uniformly refined as described in section 3.1. The analytical steady solution writes [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF]:

ρ = 1, u = -sin(θ) u θ , v = cos(θ) u θ , p = 1 + 1 75 2 r 2 2 -32 ln(r) - 128 r 2 , (35)
where u θ = 1 75 (-r + 16 r ), r = x 2 + y 2 and θ = tan -1 ( y x ). The analytical solution is provided as exterior state for inlet and outlet, which are located on straight boundaries (see Fig. (3)). Wall conditions are prescribed on the two cylinders. Time integration is achieved until convergence, based on the monitoring of the L 2 (Ω) norm of the error for all solution components. The evolution of the error for the energy field with respect to the number of degrees of freedom is shown on Fig. [START_REF] Cottrell | Isogeometric analysis : towards integration of CAD and FEA[END_REF], using both a geometrically exact grid based on NURBS (left figure) and an approximated grid based on a classical piecewise-linear representation (right figure). As can be seen, the proposed approach based on NURBS exhibits a convergence rate of value p + 1 (optimal value when using a classical DG method of degree p), whereas the use of a piecewise-linear mesh yields a second-order accuracy due to geometrical errors. More precisely, the examination of the solution on Fig. [START_REF] Bazilevs | Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly-enforced boundary conditions on unstretched meshes[END_REF] shows that spurious rarefaction waves are generated at the wall boundary, for each interface between elements, due to the discontinuity of the boundary normal when the grid is piecewise-linear. maximum width of the nozzle is 2 and its minimum 1.6. The flow is supposed to be symmetric with respect to y = 1 and, therefore, only half of the geometry is represented as computational domain. As inlet, a density of value 1.4 and a velocity of value 0.65 are weakly imposed, while a pressure of value 1 is weakly prescribed at outlet. The solution is initialized using these values, yielding a reference Mach number M 0 = 0.65. Since the geometry of the nozzle is described by a single NURBS curve, one could easily define an initial computational domain by a single NURBS patch. After Bézier extraction procedure, one obtains a set of 8 × 4 rational Bézier elements of degree p = 4 as initial grid, as illustrated in Fig. [START_REF] Hosseini | Isogeometric analysis of the navier-stokes equations with taylor-hood b-spline elements[END_REF]. The adaption procedure is achieved using refinement and coarsening criteria set to ref = 2 and coa = 0.5 respectively. The maximum refinement level is l max = 4. The parameter controling the artificial viscosity switch is set to s0 = 10. Some tests have shown that this parameters can vary significantly without impacting the results.

The time evolution of the energy field and the corresponding adapted grid are shown in Fig. [START_REF] Persson | Curved mesh generation and mesh refinement using lagrangian solid mechanics[END_REF]. As observed, some waves are generated at initial time at the bump location, which are then propagated and reflected. The flow finally converges towards a solution characterized by a steady shock. The adaption procedure allows an accurate description of the wave system and a shap capture of the shock. One can underline that the exact-geometry property, permitted by the use of a NURBS-based grid, allows to employ a very coarse grid at the bump when the flow is regular. Clearly, it would not be possible to use such a coarse mesh with classical piecewise linear elements. As illustrated in Fig. [START_REF] George | Construction of tetrahedral meshes of degree two[END_REF], the shock is located in a single element, as expected using the sub-cell capturing approach, and the artificial viscosity field is non-zero only in a very narrow region around the shock.

Subsonic viscous flow around a cylinder

As second test-case, we consider the subsonic flow around a cylinder, with a Reynolds number R ∞ = 200 and a Mach number M ∞ = 0.1 as free-stream conditions. The computational domain is the rectangle [-30, 40] × [-30, 30], the cylinder radius being equal to one. The initial grid is composed of only 15 rational Bézier elements of degree p = 3, as shown in Fig. [START_REF] Abgrall | A method for computing curved meshes via the linear elasticity analogy, application to fluid dynamics problems[END_REF]. It is constructed by defining first a set of four rational quadratic elements around the cylinder, as explained in section 3, which are then transformed into cubic elements by degree elevation [START_REF] Piegl | The NURBS book[END_REF], and finally connected to a set of surrounding cubic elements to cover the computational domain until the exterior boundary. This initial grid is extremely coarse, but very easy to define. We rely then on the CAD-consistent adaptive procedure to construct automatically a suitable computational grid. The adaption parameters are chosen as ref = 2 and coa = 0.5, and different maximum refinement levels are tested l max = 2, l max = 4 and l max = 6. Note that a small rotation of 0.5 degree is carried out to favor the trigger of the vortex shedding process.

The results are plotted in Fig. [START_REF] Geuzaine | The generation of valid curvilinear meshes[END_REF], at non-dimensional time t f in = 350, in terms of adapted grid and density field, with a zoom in the cylinder wake. As can be seen, the adaption algorithm targets for refinement both the vicinity of the cylinder and the wake. For l max = 2, the discretization of the cylinder is composed of only 16 faces, which remains very coarse. Nevertheless, the vortex shedding process is captured. However, the description of the vortices in the wake is very poor and discontinuities of the solution at element interfaces are signficant. For l max = 4, the flow description is far better. In particular, the density field in the vicinity of the cylinder looks smooth and vortices are now well represented, although some interface discontinuities can still be observed in the wake. One can notice that the refinement area follows correctly the oscillation of vortices in the wake. Finally, for l max = 6, a very accurate and smooth representation of the density field is obtained everywhere, thanks to the refinement that fits perfectly the characteristics of the wake.

Supersonic inviscid flow around a cylinder in a duct

Finally, the supersonic flow around a cylinder in a duct is considered. The flow is supposed to be inviscid and inlet boundary conditions correspond to M ∞ = 2. A cylinder of diameter D is located between two walls forming a duct of height H = 8D, as studied in [START_REF] Papoutsakis | An efficient adaptive mesh refinement (AMR) algorithm for the discontinuous galerkin method: Applications for the computation of compressible two-phase flows[END_REF]. An initial grid is constructed in the spirit of the one used in the previous section. It counts only six rational Bézier elements of degree p = 3, as depicted in Fig. [START_REF] Bazilevs | Isogeometric analysis using T-splines[END_REF]. The flow is characterized by the generation of a bow shock at the front of the cylinder, which is reflected symmetrically on the duct walls. In the cylinder wake, some vortices are generated, which interact with the reflected shocks to form complex vortical structures, as shown in Fig. [START_REF] Xu | Parametrization of computational domain in isogeometric analysis: methods and comparison[END_REF] for a fine mesh, at time t f in = 12.5 D/c. The objective of this test-case is therefore to study the capability of the method to accurately capture, by local refinement, both shocks and vortices, starting from a very coarse discretization.

In this perspective, the maximum number of refinement levels is set to the value l max = 6, with varying refinement and coarsening criteria. Additionally, a fine mesh is generated by applying the six refinement levels uniformly. Therefore, the smallest elements have the same size in the two cases. The grid parameters are summarized in table [START_REF] Silveira | Higher-order surface treatment for discontinuous Galerkin methods with applications to aerodynamics[END_REF]. The parameter controling the artificial viscosity switch is set to the value s0 = 5. The flow solution and the grid are shown at time t f in in Fig. [START_REF] Xu | Parametrization of computational domain in isogeometric analysis: methods and comparison[END_REF] for the fine mesh, and in Fig. [START_REF] Xu | Analysis-suitable volume parameterization of multi-block computational domain in isogeometric analysis[END_REF] for the adapted grids. For all cases with adaption, one can notice that the grid is quite coarse in the vicinity of the cylinder, which is not dommageable thanks to the exact-geometry property of the proposed approach. For the adapted grid number 1, characterized by large ref and coa values, the refinement is focused around the bow shock, which is accurately captured. However, the regions corresponding to the reflected shocks and the cylinder wake are below the refinement threshold for the last levels. As consequence, the grid remains too coarse in these regions to capture correctly these phenomena. For the adapted grid 2, a decrease of the refinement criteria yields a far better grid. In particular, the reflected shocks are well captured, as well as the interactions with the vortex shedding process. As can be seen in table [START_REF] Silveira | Higher-order surface treatment for discontinuous Galerkin methods with applications to aerodynamics[END_REF], this improvement is obtained for a significant increase of the number of elements, which remains however far lower than that of the uniform fine grid. Nevertheless, some discrepencies with the reference results can still be observed in the region where shocks and vortices interact. They can be explained by two reasons: firstly, the generation of the vortices is initiated in two different ways ; secondly, the use of the adaption procedure seems to promote the generation of larger vortices, due to successive refinement / coarsening steps during time evolution. Finally, the use of the adapted grid 3 yields results very close to the reference ones. Indeed, using lower refinement criteria prevents coarsening in the wake region. As consequence, the transport of small vortices is better represented. The comparison of the computational times, provided in table [START_REF] Silveira | Higher-order surface treatment for discontinuous Galerkin methods with applications to aerodynamics[END_REF], shows a signifiant increase of the computational cost as the refinement criteria are reduced. As consequence, the values of coa and ref can be adjusted to balance the accuracy and the computational effort, depending on the objective pursued. However, we underline that, if one compares the last configuration with adaption and the uniform fine mesh, a significant reduction of the computational time can be observed for very similar solutions, despite the additional cost related to the adaption procedure.

Conclusion

In this work, a CAD-consistent adaptive and refinement procedure was presented, in the framework of a DG method for compressible Euler and Navier-Stokes equations. In particular, we explained how intrinsic properties of NURBS representations can be used to generate a DG-compliant grid from CAD data, based on rational Bézier elements, and to define refinement and coarsening strategies. It should be underlined that the modifications of the DG method are limited and only concern the definition of the basis functions. If an existing DG code is written in a generic way, implementing the proposed approach may represent a small coding effort in practice. We also presented a modification of the sub-cell shock capturing method to account for the parametric representation employed. The accuracy of the method was verified for two problems with analytical solutions, one with curved boundaries and one including adaption, and convergence rates close to the value p + 1 were found. Finally, a set of flow problems were tested, to assess the capability of the method to capture accurately flow characteristics. As shown, the adaption strategy was able to consider a range of flow regimes, from subsonic to supersonic, and detect regions of interest like shocks and vortices, although the error indicator used was very simple. Obviously, the accuracy and the computational cost depend on the number of refinement steps allowed, as well as the refinement / coarsening criteria used. It was reported that a good choice of these parameters may be tedious when different phenomena are simultaneously present, like shocks and vortices in the case of the supersonic flow around a cylinder. The use of curved geometries associated to refinement / coarsening capabilities was especially interesting. Indeed, it was possible to use very coarse initial grids and rely on the adaption procedure to define a suitable computational domain. Moreover, we observed that, in some cases, refinement was not necessary at wall boundaries, provided that the geometry was exactly represented by coarse curved elements.

Efforts should now be focused on the achievement of more complex problems, for instance three-dimensional or turbulent flows, which will require to develop more efficient tools to manage CAD geometries in practice. Coupled problems, like fluid-structure interactions, may also benefit from the developed methodologies.

  a) Five NURBS functions defined on three knot intervals. Nine rational Bézier functions obtained by multiple knot insertions.

Figure 2 :

 2 Figure 2: Procedure to transform a NURBS basis to a set of DG-compliant rational Bézier basis.

Figure 3 :Figure 4 :

 34 Figure 3: 2D illustration of Bézier extraction.

  (a) Initial CAD-consistent grid. (b) Refined grid.

Figure 5 :

 5 Figure 5: Illustration of the local refinement procedure.

Figure 6 :

 6 Figure 6: Flux quadrature in case of local refinement.

  (a) Curved boundaries. (b) Piecewise linear boundaries.

Figure 7 :

 7 Figure 7: Error assessment for the cylinder case.

  (a) Curved boundaries. (b) Piecewise linear boundaries.

Figure 8 :

 8 Figure 8: Energy field using bases of degree five for the cylinder case.

  (a) Error versus number of degrees of freedom. (b) Error versus computational time

Figure 9 :

 9 Figure 9: Error assessment for the isentropic vortex case.

Figure 10 :

 10 Figure 10: Momentum-x field and adapted grid.

Figure 11 :

 11 Figure 11: Initial grid for the nozzle test-case.

Figure 12 :

 12 Figure 12: Transonic flow test case at selected time steps.

  (a) Energy field (b) Artificial viscosity field

Figure 13 :

 13 Figure 13: Zoom around shock location at final time.

Figure 14 :

 14 Figure 14: Initial grid for viscous cylinder test-case.

Figure 15 :

 15 Figure 15: Subsonic viscous test case for l max = 2 (top), l max = 4 (middle) and l max = 6 (bottom).

Figure 16 :

 16 Figure 16: Initial grid for supersonic test-case.

Table 2 :

 2 coa refDOF number (at t f in ) CPU (hRefinement parameters for the supersonic test-case.

Figure 17 :

 17 Figure 17: Reference results for the supersonic test-case.

Figure 18 :

 18 Figure 18: Results with adaption for the supersonic test-case for different adaption parameters

Unsteady problem with adaption

As second problem, we consider the transport of an isentropic vortex [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods[END_REF], whose analytical solution is provided by:

with r = (x -t -x 0 ) 2 + (y -y 0 ) 2 , x 0 = 5, y 0 = 0 and β = 5. The physical domain Ω is the square

The analytical solution is used to define the initial condition and the exterior state for boundary fluxes. Error is computed at final time t f in = 2. The assessment is first carried out using regular grids, for degrees ranging from one to four. Then, the proposed adaptive refinement and coarsening procedure is tested, using different maximum levels, for the same basis degrees. The parameters that drive the adaption are chosen as ref = 1 and coa = 0.5. The L 2 -norm of the error is computed by numerical quadrature, for the energy field. Figure (9a) shows the evolution of the error with respect to the number of degrees of freedom in log-log scale, whereas figure (9b) depicts the error versus the computational time. Table [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible navier-stokes equations[END_REF] provides, for the finest computations, the number of degrees of freedom, the error, an estimate of the convergence rate and the computational time. For the adaptive cases, the number of degrees of freedom corresponds to the final time and the maximum refinement level l max = 3. The convergence rate is estimated from the slope of the error curve in Figure (9a). Globally, one observes that the convergence rate obtained using uniform grids is slightly above the value p + 1. When using the adaptive procedure, the number of degrees of freedom is drastically reduced, while the error is maintained at a similar level. The convergence rate is tedious to evaluate in the case of non-uniform adaptive grids, nevetheless its estimation seems to be close to the value p + 1. Obviously, the adaption procedure yields an additional computational cost, due to the refinement of elements, update of connectivity tables, projection of the solution and evaluation of the inverse of the mass matrix. Nevertheless, the adaption procedure permits a significant reduction of the computational time for a given accuracy. Figure [START_REF] Nortoft | Isogeometric analysis of Navier-Stokes flow using locally refinable B-Splines[END_REF] illustrates the solution obtained using three refinement levels.

Uniform

Applications

Transonic inviscid flow in a nozzle

We consider as first application the inviscid flow in a nozzle at transonic regime. The geometry of the nozzle is defined using a single NURBS curve of degree four. Inlet and outlet are respectively located at x = 0 and x = 4. The