domain: spi phys

text:
Hubert M J Cantalloube

Philippe Martineau

Joseph Martinot-Lagarde

Nicolas Castet

Maritime circular SAR processing issues: Antenna pattern verification, sea-clutter analysis & ship image autofocus

During two SAR acquisition campaigns in 2015 and 2017, several circular SAR acquisitions have been performed on sea surface and ships. The 2015 set is full polarimetric acquisition at both X & L band (8 channels with 300 & 150 MHz bandwidth each). The 2017 set is vertical polarisation at X band (two channels with 2 GHz bandwidth each). Several processing issues had to be addressed with these experiments: -Simultaneous antenna pattern & sea clutter backscatter table evaluation was required because of antenna pattern compensation abnomalies at X band. -The autofocus used above ground was useless on the moving sea surface, hence the ships images required autofocus with tiepoints on the ship itself to compensate her own motion. During the 2017 campaign, low resolution images were synthesised on board to assess correct target acquisition during flight.

Context

In 2015, ONERA participated in an oil-spill detection experiment in Norway [START_REF] Viallefont-Robinet | Hyperspectral and radar airborne imagery over controlled release of oil at sea[END_REF], for which it performed stripmap polarimetric SAR acquisitions at X and L band. During that experiment (and especially its preparation test-flight) several circular acquisitions have been performed. The motive for these, was both science (study of the evolution of backscatter with incident angle & heading angle w.r.t. the wind), and helping to improve incorrect antenna pattern compensation at X-band (as sea is a very uniform landscape). Antenna pattern determination from acquisition signal have been done in 2001 for the P-band antenna, because the pattern measured in an-echoic chamber was significantly altered by the surrounding aircraft cell. Acquisitions on an extended forested area [START_REF] Pastoré | Imagerie radar par synthèse d'ouverture en basse fréquence[END_REF] was used for that purpose assuming a Lambertian behaviour of forest at this wavelength. However, The method needed adapting to the highly varying and unknown sea surface RCS (backscattering) namely, the sea surface RCS increases at low incidence angle and around the upwind/downwind heading. In 2017, a spare flight in campaign of high resolution imaging, was used for testing both the system on maritime targets and the onboard processing. The campaign was aimed at 2 GHz bandwidth acquisition of ground targets. Our system has the capacity of covering the full (8 to 12 GHz) X-band, with two channels spanning upper and lower X-band (8-10 Ghz & 10-12 GHz resp.) with three interleaved 720 MHz pulses (there is a 50 MHz overlap between adjacent pulses in order to fine tune clock drifts between local oscillators and/or between the two channels). Technically, the lower channel is downmodulated using the higher side band, while the upper channel is down-modulated using the lower side band (thus its signal is conjugated w.r.t. the signal of the lower channel). Since the required 2 GHz bandwidth corresponded to one pulse of the lower channel & two pulses of the upper channel we needed to digitise both channels see [START_REF] Du | Ramses-ng x-band uhr capability[END_REF] for details on the current RAMSES system capabilities).

simultaneous antenna pattern and sea surface back-scatter curve evaluation

The principle of simultaneous antenna pattern and sea surface back scattering coefficient evaluation is to start from a serie of SAR images in range Doppler coordinates each obtained from a fixed aperture as Figure 1 illustrate. In order to differentiate from stripmap SAR images (for which the aperture drifts from image row to image row), these images are further denoted flashlook SAR images. Such an image geometry is used for each burst during satellite SAR-SCAN processing (and the burst for the same range span are later seamed), and is also sometimes used for spotlight SAR acquisitions. Circular SAR images are generally computed in this geometry (though section 4 is a counterexample) since the stripmap geometry is not sound for imaging with a swath extending both near and far side of the trajectory centre (the successive image row would progress in the aircraft flight direction in the near range and in the opposite direction beyond the centre!).

Assuming short aperture hence no significant attitude changes during integration, illumination of each image is close to the illumination at the middle of the aperture.

Hence if we map each image to the antenna axes, assuming a "flat" ocean surface (as in Figure 2), we have an image sequence of a moving sea surface shaded by a fixed antenna pattern. If the imaged landscape were an homogeneous perfectly Lambertian surface, compensating the cosine of the incidence angle would directly provide a measure of the antenna pattern. Tropical rain forest provides extended areas which are very close to this ideal case, but (unlike SAR satellites) most aircraft SAR systems do not have frequent access to such areas, while sea or great lakes surfaces are more easily accessible. However, water surface is far from Lambertian: though quite homogeneous, its RCS increases dramatically at low incidence angle, and the anisotropy of waves induces two RCS peaks in heading at the downwind and upwind directions. (In case of wave reflections close to the shore, or in case of recent change in wind directions, more secondary RCS peaks in heading could be observed).

Since the antenna pattern is by definition constant in the antenna coordinates, the fluctuations in image level at a given pixel correlated with the incident angle (& heading) fluctuations provides a measure of the derivatives of the back-scattering curve at a given incidence (& heading) angle. This derivative can be integrated yielding a back scattering curve (Figure 3). Note, this curve, obtained by integration, is offset by an arbitrary constant.

Of course, for our purpose, the antenna beam should not be steered, hence aircraft roll & yaw fluctuations induce changes in incidence (& optionally heading) at any pixel location in the mapped sequence. Once the back scattering curve is known, the incidence angle variation can be compensated for in the sequence, and the average of all images in the sequence after incidence angle compensation gives an estimate of the antenna pattern (again offset by an arbitrary constant Figure 4). On this figure, we observe a bias in elevation of reconstructed the antenna pattern with respect to the anechoic chamber measures. Further inquiries have shown that a rotary sensor of the antenna was mis-aligned during the 2015 campaign due to vibrations inducing a systematic error in the antenna elevation, and introduction of this bias in the antenna pattern compensation significantly improved the antenna pattern compensation quality.

The same approach also helped in detecting a mecanical damage to one of the parabola during expedition of the X-band component for the 2017 campaign. This method could be used on linear acquisitions (as a matter of fact in [START_REF] Pastoré | Imagerie radar par synthèse d'ouverture en basse fréquence[END_REF] linear acquisitions on a conifer forest was used for antenna pattern retrieval) provided aircraft attitude fluctuation are significant, but the heading dependency of the RCS (much more significant for sea surface than for forests) makes the method less robust than with circular acquisition (the technical reason is that the heading is correlated with the left-right position in the antenna coordinates in linear case, and not on circular case). The drawback of circular acquisition is that the RCS variation of sea surface should be estimated for both incidence and heading angles.

ship autofocusing

Ships are focused using the algorithm described in [START_REF] Cantalloube | Sar retrieval of a ship vertical profile from her roll and pitch motion[END_REF] which is derived from an autofocus technique described in [START_REF] Hubert | Multiscale local map-drift-driven multilateration sar autofocus using fast polar format image synthesis[END_REF] but with tie-points on the moving target. In order to have an image resolved enough to isolate at least three tie-points, the mean velocity of the ship (if running) is evaluated e.g. by contrast optimisation autofocus Once the sensor trajectory relative to the ship frame is determined, the focused ship images can be synthesised. Note that the relative trajectory is non linear even in the case of a linear acquisition (in fact, figure 7 and 8 are from [START_REF] Cantalloube | Sar retrieval of a ship vertical profile from her roll and pitch motion[END_REF] and acquired in a previous campaign along a linear trajectory. In the definitive version of the paper, it will be replaced by images obtained on the circular acquisition illustrating next section). The strong non-linearity caused by ship roll & pitch and to a lesser extend by aircraft trajectory circle curvature may impose to use 3D processing to mitigate the focus dependency to height above sea surface [START_REF] Cantalloube | Three-dimensional processing for (circular) high resolution sar imaging of non planar targets[END_REF].

onboard processing

For the 2017 campaign, as the radar hardware was redesigned, room and electrical power become available to accommodate a third operator seat with computer workstation equipped with a Kepler GPU, linked to the signal recorder by an optical fibre, and recording the inertial navigation unit (INU) messages from the CAN bus.

Between the acquisitions, the processing computer reads back part of the signal from the recorder and synthesise images. It is possible to read back the full bandwidth (2 channels × 3 pulses of 720 MHz) in principle. However, the GPU is just powerfull enough to process two the subbandwidths (1.4 GHz bandwidth) and furthermore, the real-time trajectory obtained from the INU is not accurate enough to yield 10 cm resolution in azimuth. Since the autofocus application was not paralleled in June 2017 (and has not yet been ported to GPU to date) it was too slow to use in the small time slot between two successive acquisitions. Hence, we only processed one pulse of one channel (thus a 720 MHz bandwidth) yielding images with 20 cm resolution in a few tens of seconds after the signal is read back (Figure 9). On the ground, signal transfer to the processor (limited to 1/6 of the bandwidth) during acquisition was tested but not experimented in flight. The transfer rate to the processor is limited to 25% of that of the recorder because the link with the recorder is just one optical fibre, while the link between radar backend and the recorder uses four fibers in parallel.

Only the stripmap processing is ported to GPU (it implements an ω -κ frequency domain algorithm with space varying motion compensation). As mentioned above, offline processing of circular acquisitions is generally done using the flash-look processor (a two stage factored back-projection time-domain algorithm) or sometimes with polar format algorithm (a frequency domain processor limited in image size). None of these algorithm are ported to GPU yet, hence for processing circular acquisition onboard, we process a small arc of the circular acquisition as a linear acquisition with a curved trajectory (motion compensation can manage it as long as the arc is short, but for 20 cm resolution the integration is 4 • , hence an arc of 10 • gives a decent image). Figure 10 and 11 show images thus obtained for a ship not running and a running ship (that of Figure 5). The ship velocity is estimated and compensated for using the interactive tool illustrated on Figure 5 right (which applies a simple refocusing on a small part of the image, such as altitude, velocity, nearrange error).

Perspective

On board processing capabilities are demonstrated here, however, true real time processing (i.e. processing starting before the end of acquisition) with simultaneous recording of the raw signal remains to be proven on flight. Further work could be to transmit the processor station screen to the other operators and to the ground using a RF datalink. With a datalink and a remote control protocol such as VNC, it could be even envision to have the processing operator on the ground that can zoom to details on the computed image (transmitting the full images at native resolution would require an prohibitory data rate).

Figure 1 :

 1 Figure 1: Flashlook image (Vv at X-band component) during a 2015 circular acquisition. Images in "ground" coordinates. Illumination is upwind (from top right corner).

Figure 2 :

 2 Figure 2: Flashlook image in antenna coordinates. Antenna pattern is fixed in this frame. This geometry is similar to what would see an optical camera tied to the antenna.

Figure 3 :

 3 Figure3: Sea RCS as a function of incident angle (left) and heading (right). The downwind RCS peak is higher and narrower than the upwind peak. The asymmetry of both peak is possibly due to a change in direction between the gust and steady winds.

Figure 4 :

 4 Figure 4: Resulting antenna pattern in antenna coordinates (left). Comparison of the azimuth (middle) and elevation (right) slice of the result and anechoic chamber measures.

Fig- ure 5 .

 5 Note on this figure, the focus blur on the ship due to the along track component of its velocity (it is visible, though only 300 MHz of the bandwidth is processed on this quick-look image), and the cross-range position aberration (the ship appears offset to the right of her wake). The red box on the full image on the background indicates the true position of the ship, the lower foreground window shows the refocused ship image (by adjusting the across range velocity component). On this image, the velocity vector is overlaid in red (and this vector is aligned to the hull axis by adjusting the along range velocity component, the alignment of the true target position with her wake further confirms this velocity component).

Figure 5 :

 5 Figure 5: Raw flashlook quick-look image (no pattern compensation) with a running ship (left). Coarse estimation of the ship velocity and heading (right).

Figure 6 :

 6 Figure 6: Raw flashlook images after coarse ship velocity compensation at the beginning of the circular acquisition (left) and after 20 s course (right).

Figure 7 :

 7 Figure 7: Unfocused ship image with the 4 details focused and tracked to refine sensor trajectory relative to ship hull.

Figure 8 :

 8 Figure 8: Refocused ship image (colour indicates the focusing height above see surface).

Figure 9 :

 9 Figure 9: Stripmap SAR image with 20 cm resolution obtained from a linear acquisition on the airport just after takeoff in the aircraft.

Figure 10 :

 10 Figure10: Circular SAR image of a ship not running (adrift). The image was also computed with the GPU version of the stripmap processor, by feeding it with the signal from a small arc of the full circle.


text:
Pierre Schickele
email: pierre.schickele@onera.fr

Xavier Ferrieres
email: xavier.ferrieres@onera.fr

Jean-Philippe Parmantier
email: jean-philippe.parmantier@onera.fr

FEM-MTLN Hybridization Technique to Evaluate Electrical Current on Multiconductor Cables inside Enclosures Illuminated by a Plane Wave

Keywords: EM modeling, Maxwell's equations, multiconductor cables, hybridization, multiconductor transmission lines, FEM, field-to-transmission-line, EMC

In this paper, we introduce a strategy to evaluate current on multiconductor electrical cables installed inside enclosures. This strategy consists of the hybridization of a FEM and a MTLN approach. The motivation of this strategy relies on the fact that the high order FEM scheme has less numerical dispersive effect than a usual first order scheme such as an FDTD scheme. More accuracy is thereby obtained in the MTLN response. This property is demonstrated on a generic test case configuration.

I. INTRODUCTION

Marketing industrial vehicles requires compliance with Electromagnetic Compatibility (EMC) standards. In this context, we are interested in the problem of currents induced on electrical cables installed in a vehicle illuminated by a electromagnetic (EM) plane wave. A commonly adopted approach to solve this problem consists in solving Maxwell's equations coupled with Multiconductor Transmission Line Network (MTLN) equations to model respectively the EM fields in the 3D structure and the electrical currents on complex multiconductor cables. When reaction on the induced currents on the incident fields (taken as sources to drive those currents) are neglected, this approach is known under the name "Fieldto-Transmission-Line" (FTL) and allows modeling of large wiring systems of industrial complexity because of its capacity to decouple 3D and MTLN problems as well as to account for MTLN models (whereas cable models embedded in 3D models are usually limited to thin single-wire models, which is a significant simplification). In EMC studies, the Yee scheme applied on a Cartesian mesh to evaluate Maxwell's equations is a widely used approach because of its numerous advantages. However, this technique is known to introduce a numerical dispersion error on the EM fields driving themselves an error on the induced currents on the cable models. To avoid this problem, one solution is to reduce the size of the mesh cells, but this strongly increases the computational costs. Therefore an other idea consists in applying a higher order scheme such as the Finite Element Method (FEM) to solve both Maxwell's equations, still on a Cartesian grid and MTLN equations. This approach gives more accuracy with fewer cells in the mesh and maintains reasonable computational costs. Moreover, by taking into account a higher order of approximation for the solution, calculations can be done at higher frequencies using the same cell size in the mesh. In the first and second sections of this paper, we describe the FEM method for Maxwell's equations and the MTLN equations. In the third section, we describe the hybridization principle between the two methods. Finally, in the last section, we present preliminary results obtained for a wire in a generic enclosure.

II. HIGH SPATIAL ORDER FEM METHOD TO SOLVE MAXWELL'S EQUATIONS

To evaluate EM fields (E, H) on Ω × [0, T ], we solve Maxwell's equations, on a computational domain Ω approximated by a Cartesian mesh with a high order (r) Finite Element Method (FEM) [START_REF] Deymier | New High order FDTD method to solve EMC problem[END_REF].

As an illustration, (1) shows the approximation used for fields. In the FEM scheme, for E x and H z components :

E x (t, x, y, z) = r i,j,k=1 E x ijk (t)L G i (x)L GL j (y)L GL k (z) H z (t, x, y, z) = r i,j,k=1 H z ijk (t)L G i (x)L G j (y)L GL k (z) (1)
where L G i and L GL i are Lagrange polynomials defined on [0, 1] respectively by using r Gauss ((ξ g k ) k=1...r ) and r + 1 Gauss-Lobatto ((ξ gl k ) k=1...r+1 ) points. These points are also the quadrature points used to evaluate the integrals in the FEM scheme. By considering this approximation, we obtain a numerical scheme which provides an efficient and accurate solution for Maxwell's equations on EM fields.

III. MTLN FORMALISM FOR COMPLEX HARNESS

The organization of cables inside vehicles define networks of Multiconductor Transmission Lines (MTL) connected to nodes (also called "junctions") as presented in Fig. 1. On each MTL, we define unknown currents I i and voltages V i distributed along the conductors of the MTL that verify (2). C i , L i , R i and G i are the capacitance, inductance, resistance and conductance matrices of the MTL and E i the incident field. To close (2), we introduce on each conductor some limit conditions at the ends of the MTLs using several models of junctions that represent electrical circuits as illustrated in Fig. 2. To approximate (2), we decompose each MTL into several segments where the unknowns I and V are located respectively at Gauss and Gauss-Lobatto quadrature points. By using the Lagrange polynoms previously defined for the EM fields, we have :

   C i ∂ t V i + G i V i = -∂ l I i L i ∂ t I i + RI i = -∂ l V i + E i V i (t = 0) = I i (t = 0) = 0 (2)
I(t, x) = r i=1 I i (t)L G i (x) and V (t, x) = r+1 i=1 V i (t)L GL (x)
(3) As for Maxwell's equations on EM fields, this high order FEM to solve MTLN equations gives an efficient and accurate solution with low dispersion.

IV. HYBRIDIZATION PROCESS

The coupling terms between Maxwell's equations and the MTLN equations are given by the electric fields along the MTLs for the MTLN equations and the currents on the wires for Maxwell's equations [START_REF] Muot | A combined FDTD/TLM time domain method to solve efficiently electromagnetic problems[END_REF]. In this coupling process, the electric fields evaluated along the direction of the MTLs are taken as distributed generators of the MTLs and the J source term in Maxwell's equations is approximated by the sum of currents on the conductors of the MTLs, for the unknown fields located in the cell where the MTLs are located.

V. TEST CASE

In this preliminary version of the paper, we consider a perfectly metallic enclosure with an aperture. An electric wire connected to 50 Ω loads to the ground is routed inside the enclosure. The enclosure is illuminated by a Gaussian plane wave (of mean 3e-10 s and standard deviation 1e-10 s) at the level of the aperture, as shown in Fig. 3. The electrical current has been calculated in the middle of the wire using "full" FDTD and 3-order FEM models (for which the wire model is imbedded in the solver) on the one side (Fig. 4) and the full FEM model and a FEM/MTLN hybridized model on the other side (Fig. 5). FDTD and FEM models have a 0.5 cm cell-size and a 1.5 cm cell-size meshes respectively. Fig. 4 clearly shows a shift in time on the FDTD solution due to FDTD numerical dispersion. Fig. 5 shows that the hybridization technique using FEM provides almost the same results as the FEM full reference results. The slight amplitude difference is due to the fact that the hybridization used here did not account for the reaction of the wire current on the field (the FTL approximation which does not reproduce EM radiation phenomenon). The final version of the paper will provide other results on a more exhaustive test case made of several enclosures and a multiconductor cable network.

Fig. 1 .

 1 Fig. 1. Example of a multiconductor cable network modeling.

Fig. 2 .

 2 Fig. 2. Examples of electrical circuits in junctions.

Fig. 3 .

 3 Fig. 3. Wire in an enclosure illuminated by an EM plane wave.

Fig. 4 .

 4 Fig. 4. Current obtained in the middle of the wire by full FDTD and FEM methods.

Fig. 5 .

 5 Fig. 5. Current obtained in the middle of the wire by full FEM method and the hybridization of MTLN with FEM method.


text:
Fabio Naddei
email: fabio.naddei@onera.fr

Marta De La

Llave Plata

Vincent Couaillier
email: vincent.couaillier@onera.fr

Ph D Student

A comparison of refinement indicators for p-adaptive discontinuous Galerkin methods for the Euler and Navier-Stokes equations

This paper presents an analysis of refinement indicators for the numerical solution of the Euler and the Navier-Stokes equations using a p-adaptive discontinuous Galerkin (DG) method. Residual error, discretization error and feature-based indicators are studied in the context of p-adaptive DG simulations of two inviscid flow configurations, namely, the flow over a Gaussian bump and the flow past a circular cylinder. As expected, the superiority of local p-refinement as compared to uniform p-refinement is clearly observed. The results obtained highlight the efficiency of the considered refinement indicators with respect to an indicator based on the local entropy error. In particular, they identify regions where numerical errors are being produced and not those where errors are being convected. Finally the p-adaptive algorithm is applied to the simulation of the laminar flow past a circular cylinder at Re = 40. The great potential of the spectral decay 1 and VMS 2 indicators, which combine accuracy and computational efficiency, is also demonstrated.

I. Introduction

High-order discontinuous Galerkin (DG) methods [START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF][START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework[END_REF][START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF] have been gaining growing interest among the scientific community in recent years. [START_REF] Wang | High-order CFD methods: current status and perspective[END_REF][START_REF] Kroll | ADIGMA-A European Initiative on the development of adaptive higher-order variational methods for aerospace applications: results of a collaborative research project funded by the European Union[END_REF] The DG methods are based on the variational projection of the Navier-Stokes equations and combine some features of Finite Volume (FV) and Finite Element (FE) methods. Remarkable properties of DG methods include the possibility of achieving arbitrarily high-order of accuracy, even at physical boundaries and on unstructured meshes, and high parallel efficiency on multi-core architectures, thanks to the compact stencil of the numerical scheme, each element requiring only information from direct neighbours. [START_REF] Biswas | Parallel adaptive finite element methods for conservation laws[END_REF] These properties play a fundamental role in the simulation of turbulent flows for which the finest scales of turbulence must be accurately resolved over long integration periods thus requiring efficient numerical methods with low dispersion and dissipation errors.

A large reduction in the computational cost of such simulations could be achieved by exploiting the possibility of increasing the local spatial resolution of DG methods by local mesh refinement (h-refinement), and by locally increasing the polynomial order within the elements (p-refinement). Local p-refinement can reduce dissipation and dispersion errors in regions where the solution is smooth by employing high polynomial degrees, thus allowing for the accurate resolution of unsteady turbulent phenomena with a lower number of degrees of freedom (DOF) as compared to Finite Difference (FD) and FV methods. [START_REF] Chapelier | Evaluation of a high-order discontinuous Galerkin method for the DNS of turbulent flows[END_REF][START_REF] De La Llave Plata | On the use of a high-order discontinuous Galerkin method for DNS and of wall-bounded turbulence[END_REF] On the other hand, h-refinement can isolate regions characterized by geometrical and physical discontinuities.

hp-adaptive techniques are based on the definition of refinement indicators and an hp-decision algorithm. Refinement indicators, usually relying on a posteriori error estimates, identify regions where an increased spatial resolution is required. The hp-decision procedure, on the other hand, measures the level of smoothness of the solution and selects the refinement strategy (h or p). Considerable efforts have been dedicated to the development of robust a posteriori error estimators for partial differential equations. We refer to the works by Roy, [START_REF] Roy | Review of discretization error estimators in scientific computing[END_REF] Houston and Süli, [START_REF] Houston | A note on the design of hp-adaptive finite element methods for elliptic partial differential equations[END_REF] and Mitchell and McClain, [START_REF] Mitchell | A comparison of hp-adaptive strategies for elliptic partial differential equations[END_REF] for a review of common error estimation strategies and their use for hp-adaptive algorithms.

In this context three main classes of refinement indicators can be identified:

Feature based -this class of refinement indicators is based on physical or theoretical properties or on the observation that some phenomena must be fully resolved in order to obtain an accurate representation of the flow Goal Oriented -Goal Oriented techniques have been shown to be optimal for the definition of refinement indicators when a target quantity is the goal of a CFD simulation (e. g. the drag or the lift coefficient). That is to say, they provide the lowest computational cost to achieve a prescribed level of accuracy in the measure of the target functional. [START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations[END_REF][START_REF] Fidkowski | Review of output-based error estimation and mesh adaptation in computational fluid dynamics[END_REF][START_REF] Hartmann | Multitarget error estimation and adaptivity in aerodynamic flow simulations[END_REF] Their optimality reflects their ability to consider the inherent mechanism of error propagation in hyperbolic and nearly-hyperbolic problems. Their derivation requires the computation of element-wise residuals associated with the discrete solution. These residuals must be weighted by the solution of an adjoint problem derived from the discretization employed and the target quantity. One possible drawback is the sensitivity of such methods to the target selected. Moreover the solution of the adjoint problem must be performed in a refined polynomial space and requires the backward integration in time for unsteady problems, thus requiring the storage of the solution at each previous time step. Given the high computational cost of these operations, adjoint methods are not commonly used for unsteady configurations.

There exists a vast literature dedicated to the analysis of error indicators and refinement techniques. However, given the large variety of indicators introduced to date, the selection of such an indicator in the implementation of a new adaptive algorithm can be a daunting task. The aim of this research is to provide guidelines for the design of adaptive refinement algorithms with a particular focus on their suitability for unsteady turbulent flow simulations.

As previously mentioned, goal-oriented refinement indicators have been shown to be vastly superior to other error estimators for steady problems. Unfortunately, their high computational cost impedes their application to LES and DNS of turbulent flows. For this reason, adjoint-based methods are not considered in this work. On the other hand DE and feature-based indicators have already been successfully applied to LES. [START_REF] Abbà | Dynamic models for Large Eddy Simulation of compressible flows with a high order DG method[END_REF][START_REF] Mitran | A comparison of adaptive mesh refinement approaches for Large-Eddy Simulation[END_REF][START_REF] Antepara | Parallel adaptive mesh refinement for Large-Eddy Simulations of turbulent flows[END_REF][START_REF] Tugnoli | A locally p-adaptive approach for Large Eddy Simulation of compressible flows in a DG framework[END_REF] However, to the best of the authors' knowledge, there have been no attempts to extend and assess the efficiency of RE-based indicators for LES/DNS. This analysis marks the beginning of an extended comparison of different refinement indicators. Here, we will assess the performance of DE, RE and feature-based refinement indicators in guiding adaptive p-refinement for the solution of the steady Euler and the Navier-Stokes equations using a DG method based on a modal approach. We focus on their ability to reduce the computational cost necessary to achieve a prescribed level of accuracy as well as potential implementation issues. We do not compare the accuracy of different refinement indicators on the evaluation of the exact error (usually measured by the effectivity index i eff = error estimate exact error ) as different error norms are employed for different indicators.

This paper is organised as follows. In section II the considered equations and the DG method are briefly introduced. In section III a brief review of the considered refinement indicators is presented. A short introduction is provided along with their relationship with other refinement indicators employed in the literature. In section IV the direct comparison of refinement indicators on inviscid and viscous flows is provided. In section V we present our observations about possible implementation issues and the computational cost of the considered indicators. While the latter is of limited interest when considering steady flows, it will be of primary importance when using these indicators for unsteady flows computations. The main conclusions from this research are presented in section VI.

II. Model problem and discretization

The simulations presented in this work have been performed using the DG unstructured solver Aghora developed at ONERA 32 solving the compressible Navier-Stokes equations. The physical model and the DG method used in this work are recalled in this section.

II.A. Physical model

Let Ω ∈ R d be a bounded domain, where d is the space dimension, the compressible Navier-Stokes equations take the form :

∂ t u + ∇ • F c (u) -∇ • F v (u, ∇u) = 0 (1)
with associated initial conditions u(•, 0) = u 0 (•) and appropriate boundary conditions prescribed on ∂ Ω. The vector u = (ρ, ρv, ρE) T represents the conservative state variables, with ρ being the density, v = (U 1 ,U 2 ,U 3 ) T the velocity vector and E = p (γ-1)ρ + v•v 2 the total specific energy. Here, p is the static pressure and γ = C p C v > 1 the ratio of specific heats. The nonlinear convective fluxes and the diffusive fluxes in Eq. ( 1) are defined as:

F c,i (u) =        ρU i ρU 1 U i + pδ i1 ρU 2 U i + pδ i2 ρU 3 U i + pδ i3 ρU i E        F v,i (u, ∇u) =        0 τ i1 τ i2 τ i3 τ ik U k -q i        (2) with τ i j = µ ∂U j ∂ x i + ∂U i ∂ x j - 2 3
∂U k ∂ x k δ i j (3)
q i = -λ ∂ T ∂ x i (4)
T = 1 C v E - 1 2 U i U i ( 5
)
where µ is the dynamic viscosity, described by the Sutherland law, T denotes the temperature, λ = µ

C p
Pr is the thermal conductivity and Pr refers to the Prandtl number.

II.B. The DG method

The domain Ω is partitioned in a shape-regular mesh T k consisting of N non-overlapping and non-empty elements K of characteristic size h K . We also define the sets E i and E b of interior and boundary faces in

T k such that E h = E i ∪ E b . Let S p h = φ ∈ L 2 (Ω) : φ | K ∈ P p (K)
, ∀K ∈ T k be the functional space formed by piecewise polynomials of either total or partial degree at most p defined in the element, and (φ 1 K , . . . φ N p K ) ∈ P p (K) a hierarchical and orthonormal basis of P p (K) of dimension N p . The polynomial order p will in general depend on the element and will be indicated as p K when necessary. For the DG method used in this research, the orthonormal basis is obtained by applying a modified Gram-Schmidt procedure to a tensor product of monomials of degree at most p. This produces a diagonal mass matrix even on curved elements. [START_REF] Botti | Influence of reference-to-physical frame mappings on approximation properties of discontinuous piecewise polynomial spaces[END_REF] The solution in each element is thus expressed as

u h (x,t) = N p ∑ l=1 φ l K (x)u l K (t), ∀x ∈ K, K ∈ T k , ∀t > 0 (6)
in which the polynomial coefficients (u l k ) 1≤l≤N p represent the DOFs of the discrete problem in element K. The semidiscrete variational form of Eq. (1) therefore reads: find u h in S p h such that ∀φ h ∈ S p h

T K φ h ∂ t u h dV + L c (u h , φ h ) + L v (u h , φ h ) = 0 (7)
In Eq. ( 7) L c and L v represent the variational projection of the convective and viscous terms onto the functional space S p h . For a given interface e in E i we define the average operator { {u} } = (u + + u -) /2 and the jump operator

[[u]] = u + ⊗ n -u -⊗ n,
where n is the local outward unit normal. The DG discretization of the convective terms then reads

L c (u h , φ h ) = - T K F c (u h ) • ∇ h φ h dV + E i [[φ h ]]H c (u + h , u - h , n)dS + E b φ + h F c (u b ) • ndS (8)
The boundary values u b = u b (u + h , u ext , n), with u ext a reference external state, are computed in order to impose the appropriate boundary condition on each boundary in E b .

The function H c is a numerical flux that approximates the convective flux on an element face and must be chosen to be conservative and consistent. In the simulations presented in this work the local Lax-Friedrichs (LLF) flux and the Roe flux have been employed. For the discretization of the viscous term in Eq. ( 7) the BR2 method of Bassi and coworkers [START_REF] Bassi | A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows[END_REF][START_REF] Bassi | Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k-ω turbulence model equations[END_REF] was employed, namely:

L v (u h , φ h ) = T K F v (u h , ∇ h u h + L L L h ) • ∇φ h dV - E i [[φ h ]]{ {F v (u h , ∇ h u h + η e l l l e h )} } • ndS - E b φ + h F v (u b , ∇u b + η e l l l e h ) • ndS (9)
where η e is a numerical parameter that ensures the stability of the method. The global lifting operator L L L h is defined as the sum of the local lifting operators l l l e h : L L L h = ∑ e∈∂ K l l l e h where l l l e h has support on the elements adjacent to e ∈ E i and is obtained from the solution of the following problem on the internal faces:

K + ∪K - φ h l l l e h dV = - e { {φ h } }[[u h ]]dS ∀φ h ∈ S p h (10)
A similar expression for the local lifting coefficients, consistent with the boundary conditions, is obtained on E b . The semi-discrete equation Eq. ( 7) is discretized in time by means of the implicit backward-Euler scheme relying on the GMRES method and ILU(0) preconditioning for the solution of the corresponding linear system. The integrals in Eq. ( 7) are computed by means of the Gauss-Legendre quadrature with q = p + 1 + m points in each space direction, where m denotes the number of overintegration points used for dealiasing purposes.

III. Refinement and marking procedure

An adaptive p-refinement algorithm defines locally the polynomial degree p K used to approximate the solution within a given element. As mentioned in section I, in this work only stationary configurations are considered and the refinement procedure acts as follows. First a converged steady solution is obtained for a given refinement level. A global refinement indicator η glob is then computed in order to estimate the accuracy of the obtained solution, as an example the total entropy variation in an inviscid uniform flow. This quantity is then compared to a user-defined threshold ω specifying the accuracy level desired. Refinement is therefore deemed necessary if η glob > ω. In this case, element-wise refinement indicators η K are computed and used to identify elements where spatial resolution must be increased.

In this work, we use the maximum marking criterion. Namely, an element K is marked for refinement if :

η K > θ max L∈T K η L (11)
where θ ∈ [0, 1] is a user-defined marking parameter. In this work we set θ = 0.5. Moreover, we have decided to limit to 1 the maximum change in polynomial degree between two neighbouring elements. For this reason additional elements may be marked for refinement.

Once the marking is completed, a new set of p values is computed by increasing by one the polynomial degree in marked elements. A maximum local polynomial degree is enforced to avoid instability issues and the need of excessive memory storage. The simulation is then restarted from the previous solution based on this the new set of p values.

III.A. Refinement Indicators

In this section we present the refinement indicators considered in this research. First, three refinement indicators that can be classified as DE-or feature-based will be presented: the VMS indicator, the spectral decay indicator and the nonconformity error indicator. Then, two more indicators will be introduced: the residual-based and residuum-ncf based indicator. These can be described as RE-based indicators.

III.A.1. The VMS indicator

The first refinement indicator that we consider in this work was developed by Kuru et al. [START_REF] Kuru | An adaptive Variational Multiscale Discontinuous Galerkin Method For Large Eddy Simulation[END_REF] in the context of the variational multiscale simulation (VMS) approach. [START_REF] Hughes | Large Eddy Simulation and the variational multiscale method[END_REF] In the VMS framework the numerical solution vector u h is expressed as:

u h = u h,L + u h,S (12)
where u h,L and u h,S are the contributions to the numerical solution corresponding to the large-resolved and the smallresolved scales. In their work, the authors define u h,L ∈ S p L h and u h,S ∈ S p h \ S p L h being 0 ≤ p L < p the polynomial degree that defines the cut-off between large and small scales.

The refinement indicator proposed by Kuru et al. is then defined as the L 2 -norm of the momentum density of the small-resolved scales:

η K . . = K (ρv) h,S 2 dV 1 2 η glob . . = ∑ K η 2 K 1 2 (13)
This refinement indicator was therefore interpreted as an approximation of the DE associated to the coarse-scale solution measured by the L 2 -norm of the momentum density

(ρv) h,L -(ρv) L 2 (Ω) .
We define the VMS refinement indicator as:

η K . . = K (ρv) h,p -(ρv) h,p-1 2 dV 1 2 (14)
where (ρv) h,p-1 is obtained from the projection of the numerical solution on S p-1 h . This refinement indicator acts as a DE estimate for u h,p employing a lower-order solution. It can also be interpreted as a feature-based refinement indicator, measuring the 'energy' of the highest-order modes.

This estimator bears strong similarities with the error indicator developed by Mavriplis [START_REF] Mavriplis | Adaptive mesh strategies for the spectral element method[END_REF][START_REF] Mavriplis | A posteriori error estimators for adaptive spectral element techniques[END_REF] where an additional term depending on the decay rate of polynomial coefficients is also included. This decay rate is helpful in identifying underresolved regions of the flow. However, the corresponding DE estimate has been shown to produce relatively low effectivity index values when considering two-and three-dimensional problems. [START_REF] Mavriplis | Nonconforming discretizations and a posteriori error estimators for adaptive spectral element techniques[END_REF] For this reason, the refinement indicator derived by Mavriplis is not considered in this work.

III.A.2. The spectral decay indicator

The so-called spectral decay indicator was first introduced as a discontinuity sensor by Persson and Peraire [START_REF] Persson | Sub-cell shock capturing for discontinuous Galerkin methods[END_REF] for shock capturing in order to stabilize numerical simulations in the presence of discontinuities. The spectral decay indicator, also known as smoothness indicator, is defined as:

S e . . = f (u h,p ) -f (u h,p-1 ) L 2 (K) f (u h,p ) L 2 (K) (15)
where again u h,p-1 represents the projection of the numerical solution on S p-1 h and f is a function that depends on the solution. Persson and Peraire employed this smoothness indicator as a shock capturing parameter based on either the entropy or the enthalpy for f (u). Gassner et al. [START_REF] Gassner | Explicit discontinuous Galerkin schemes with adaptation in space and time[END_REF] considered Eq. ( 15) based on the momentum density in one direction and employed this estimator as a refinement indicator. Tumolo et al. [START_REF] Abbà | Dynamic models for Large Eddy Simulation of compressible flows with a high order DG method[END_REF][START_REF] Tugnoli | A locally p-adaptive approach for Large Eddy Simulation of compressible flows in a DG framework[END_REF][START_REF] Tumolo | A semi-implicit, semi-Lagrangian, p-adaptive discontinuous Galerkin method for the shallow water equations[END_REF][START_REF] Tumolo | A semi-implicit, semi-Lagrangian discontinuous Galerkin framework for numerical weather prediction[END_REF] used this refinement indicator for the simulation of various CFD configurations. Following their approach, we define the spectral decay indicator as:

η K . . = (ρv) h,p -(ρv) h,p-1 L 2 (K) (ρv) h,p L 2 (K) (16)
Equation ( 16) corresponds to the error estimator provided in Eq. ( 14) normalized by the total 'energy' of the cell defined by (ρv) h,p L 2 (K) .

III.A.3. The nonconformity error indicator

The third refinement indicator considered in this work is the so-called nonconformity (NCF) error. This indicator is based on the assumption that, unless a physical discontinuity is present, the analytical solution is continuous across element faces. This implies that a jump in the numerical solution can be considered as a a measure of the error. Krivodonova et al. [START_REF] Krivodonova | Error estimation for discontinuous Galerkin solutions of two-dimensional hyperbolic problems[END_REF][START_REF] Krivodonova | Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws[END_REF] introduced the following discontinuity sensor:

η K . . = max E ∈∂ K max j [[u h (x j )]] 2{ {u h (x j )} } (17)
where u h is a variable like density or pressure and x j are the quadrature points on the element faces. This error estimator can also be considered as a local DE estimator thanks to the superconvergent phenomena at outflow boundaries for some DG methods. [START_REF] Krivodonova | Error estimation for discontinuous Galerkin solutions of two-dimensional hyperbolic problems[END_REF][START_REF] Krivodonova | Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws[END_REF] Gassner and Altmann 40 employed this quantity as a refinement indicator. Following their approach, we indicate as nonconformity error indicator:

η K . . = max E ∈∂ K max j (ρv h ) + (x j ) -(ρv h ) -(x j ) 2 { {(ρv h )(x j )} } (18)

III.A.4. The residual-based indicator

The residual-based refinement indicator that will be described here follows the classical derivation of residual-based error estimators for FE method [START_REF] Ainsworth | A posteriori error estimation for Discontinuous Galerkin finite element approximation[END_REF] and the derivation by Hartmann and Houston. [START_REF] Hartmann | Error estimation and adaptive mesh refinement for aerodynamic flows[END_REF] Let us write the semidiscrete form of Eq. ( 7) as:

N h (u h , φ ) = 0 ∀φ ∈ S p h ( 19
)
and let J(u) be a nonlinear target functional. Provided that the target functional satisfies a compatibility condition and the discretized operator is adjoint consistent, [START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations[END_REF] it is possible to show that:

J(u) -J(u h ) = -N h (u h , z) (20)
where z is the solution of an adjoint problem derived from the target functional J(•). The derivation of this result corresponds to the general derivation of duality-based a posteriori error estimators. Manipulating the semidiscrete operator, it is possible to obtain:

N h (u h , φ ) = ∑ K∈T K   K R(u h )φ dV + ∂ K∩E i r(u h )φ + dS + ∂ K∩E b r b (u h )φ + dS   (21)
in which the terms on the right-hand side are defined as:

R(u h ) = ∂ u h ∂t + ∇ • F c (u h ) -∇ • F v (u h , ∇ h u h + L L L h ) (22) r(u h ) = H c (u + h , u - h , n) -F c (u + h ) • n -{ {F v (u h , ∇ h u h + η e r e h )} } -F v (u + h , ∇ h u + h + L L L h ) • n (23) r b (u h ) = F c (u b ) -F c (u + h ) • n -F v (u b , ∇ h u b + η e r e h ) -F v (u + h , ∇ h u + h + L L L h ) • n (24)
Where the internal residual Eq. ( 22) measures how accurately the numerical solution within an element satisfies the continuous mathematical model, the interelement residuals Eq. ( 23) measures the jump in normal fluxes at internal faces and the boundary residual Eq. ( 24) measures the error in the imposition of boundary conditions. Provided that the adjoint solution is z| K ∈ H s (K) with 2 ≤ s ≤ p K + 1, ∀K ∈ T K , it can be shown [START_REF] Hartmann | Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations[END_REF] that the error in the target functional can be bounded as

|J(u) -J(u h )| ≤ C int ∑ K∈T K (η K ) 2 1 2 (25)
where C int is a positive constant which depends on the shape-regularity of the mesh, the polynomial degree p and the regularity of the dual solution, and

η K . . = h s K R(u h ) L 2 (K) + h s-1 2 K r(u h ) L 2 (∂ K∩E i ) + h s-1 2 K r b (u h ) L 2 (∂ K∩E b ) (26)
On the basis of numerical experiments, Hartmann and Houston employed therefore the following residual-based inticator:

η K . . = h K R(u h ) L 2 (K) + h 1 2 K r(u h ) L 2 (∂ K∩E i ) + h 1 2 K r b (u h ) L 2 (∂ K∩E b ) η glob . . = ∑ K∈T K (η K ) 2 1 2 (27)
Rivière and Wheeler 26 derived a similar error estimator for an elliptic problem, with linear diffusion and production, by employing the energy norm of the operator in place of the L 2 -norm.

III.A.5. The residuum-NCF based indicator

Dolejší [START_REF] Dolejší | hp-DGFEM for nonlinear convection-diffusion problems[END_REF] proposed a refinement indicator based on the following global error definition:

E . . = R(u h ) 2 + NCF(u h ) 2 1 2 (28)
where R(u h ) is called residuum error and NCF(u h ) is a measure of the nonconformity error. Denoting by X . . = H 2 (Ω, T K ) = φ : φ ∈ H 2 (K), ∀K ∈ T K the broken Sobolev space, these quantities are defined as:

R(u h ) . . = N h (u h , •) X = sup φ ∈X\{0} N h (u h , φ ) φ X ( 29
) NCF K (u h ) . . = ∂ K∈E i h -1 ∂ K [[u h ]] 2 dS + ∂ K∈E b h -1 ∂ K (u h -u b ) 2 dS NCF(u h ) . . = ∑ K∈T K NCF K (u h ) 2 1 2 (30)
For the analytical solution we have that N h (u, φ ) = 0, ∀φ ∈ X, while for the numerical solution this is true only for φ ∈ S p h , thus the residuum error is a measure of the residual in the weak form of Eq. ( 1) . Using Eq. ( 29) and Eq. ( 20), it can be readily shown that

|J(u) -J(u h )| = |N h (u h , z)| ≤ R(u h ) z X (31)
Therefore, limiting the residuum error has a similar effect to limiting the residual-based indicator defined by Eq. ( 27) in the evaluation of the error in a target functional. The residuum error cannot be exactly computed, as it is defined on an infinite-dimensional space. Thus, the residuum error estimator and the local or element residuum error estimator are used in practice and are defined as:

ρ h (u h ) . . = sup φ ∈S p+1 h \{0} N h (u h , φ ) φ X (32)
ρ h,K (u h ) . . = sup φ ∈S p+1 K \{0} N h (u h , φ ) φ X (33)
The spaces S p+1 h and S p+1 K correspond respectively to the spaces of piecewise polynomials of maximum degree p + 1 and its restriction to the considered element for the evaluation of the element residuum estimator, defined by:

S p+1 K . . = φ ∈ L 2 (K) : φ | K ∈ P p+1 (K), φ | K = 0, ∀K, K ∈ T k , K = K (34)
In this work, following Dolejší, [START_REF] Dolejší | hp-DGFEM for nonlinear convection-diffusion problems[END_REF] we consider the norm

• X . . = • 2 L 2 (Ω) + 1 Re | • | 2 H 1 (Ω,T K ) 1 2 ( 35
)
where Re is the Reynolds number. Other norms have also been employed in the literature. [START_REF] Dolejší | Residual based error estimates for the space-time discontinuous Galerkin method applied to nonlinear hyperbolic equations[END_REF][START_REF] Dolejší | A design of residual error estimates for a high order BDF-DGFE method applied to compressible flows[END_REF] It can be shown [START_REF] Dolejší | Residual based error estimates for the space-time discontinuous Galerkin method applied to the compressible flows[END_REF] that, based on this choice, the residuum error estimator can be directly computed from the element-wise values:

ρ h (u h ) 2 = ∑ K∈T K ρ h,K (u h ) 2 (36)
The local residuum estimator ρ h,K (u h ) is computed by seeking a solution ψ = ∑ N i=1 ξ i φ i , where {φ i } is a basis for S p+1 K and ξ i are the unknown coefficients, such that |N h (u h , ψ)| is maximum and ψ X = 1. This leads to a constrained optimisation problem which can be solved via the Lagrange multipliers techniques, which requires the computation of the residuals N h (u h , φ i ).

It should be noted that N h (u h , ψ), ∀φ i ∈ S p K should be identically null by virtue of the Galerkin orthogonality property. In practice this is not the case due to the presence of an algebraic error which is the combination of the aliasing and iterative errors. Thus, in order to evaluate this contribution, we compute the quantity N h (u h , ψ) by increasing the number of overintegration points.

In conclusion, we indicate as residuum-NCF refinement indicator :

η K . . = ρ h,K (u h ) 2 + NCF K (u h ) 2 1 2 (37)

IV. Evaluation of the performances of refinement indicators

In this section the results from the application of the adaptive p-refinement methodology based on the indicators described above are presented. Firstly, we consider the inviscid flow over a Gaussian bump at Mach number M = 0.5 and the inviscid flow past a cylinder at M = 0.3. Secondly, we analyse the laminar flow past a cylinder at M = 0.1 and Reynolds number Re = 40, based on the upstream velocity and the cylinder diameter. These configurations are subsonic and are characterised by solid boundaries free of geometrical singularities. For this reason, they are suitable for adaptive p-refinement. In fact, configurations presenting geometrical singularities or shock discontinuities are best treated by a combination of h and p refinement.

IV.A. Steady Euler flow over a Gaussian Bump

The first test case considered is the inviscid flow over a Gaussian bump at M = 0.5. This test case has been widely studied at the series of International Workshops of High-Order CFD Methods. [START_REF] Wang | High-order CFD methods: current status and perspective[END_REF] The configuration is 2D with a subsonic inlet boundary condition, a prescribed static pressure at the outlet boundary and slip conditions imposed at the upper and bottom walls. The mesh employed is a 4 th -order mesh which is generated based on the analytical expression of the bump that folows a Gaussian profile

The performance of the refinement indicators will be assessed by comparing the entropy error on the full domain, namely:

e 2 L 2 (Ω) = Ω p p ∞ ρ ∞ ρ γ -1 2 dΩ ( 38
)
where p ∞ and ρ ∞ are the static pressure and density corresponding to inlet conditions. Figure 1 shows the variation of this error measure versus the number of DOFs when uniform p-refinement is performed (solid black line with circles) as well as for the locally p-adapted simulations. For the latter, we report the results based on the refinement indicators introduced in the previous section. All the adaptive simulations have been performed starting from uniform p = 1 and limiting the maximum local degree to p = 4.

The global refinement indicator η glob is here defined as the entropy error given in Eq. ( 38) and the threshold value is set equal to a value slightly higher than that corresponding to uniform p = 4. Figure 2 illustrates the refinement levels obtained at the end of the refinement procedure.

As reference, we also report the results obtained by using the local entropy error as refinement indicator (light blue dots in Fig. 1). For inviscid subsonic flows with uniform inlet conditions, this quantity can be considered as an exact measure of the discretization error. Therefore, these results illustrate one of the possible drawbacks of DEbased refinement indicators as mentioned in section I. In fact, for this configuration it is expected that a higher spatial resolution will be required in the proximity of the bump. In underresolved simulations an error in the entropy will be produced in this region and it will be convected downstream. For this reason, as can be seen in Fig. 2, this refinement indicator will effectively select for further refinement the region behind the bump. This will cause the suboptimal performance demonstrated in Fig. 1.

With regards to the other refinement indicators,we can see from Fig. 1 that they all show similar results, in spite of their fairly different derivation, implementation, and computational cost. The residual-based, residuum-NCF and NCF refinement indicators show slightly better results after the initial refinement step with the residual-based being the best performing indicator. The residual-based indicator reduces the number of DOFs necessary to achieve the same convergence level as uniform p = 4 by ~63% while the VMS and spectral decay indicator lead to a reduction of ~58%. 10 1.75 10 1.8 10 1.85 10 1.9 10 1.95 10 2 10 2.05 10 2.1 10 2. It is also interesting to note that the VMS and NCF indicators, despite being DE estimators, do not suffer from the drawbacks illustrated for the local entropy error indicator.

Coming back to Fig. 2 it can be observed that the residual-based, NCF and residuum-NCF indicators exhibit fairly similar p-refinement patterns while those of the VMS and spectral decay indicators are nearly identical. As expected, the residuum-NCF indicator produces a refinement pattern which is in between the ones produced by the residual-based and the NCF-based estimators. Finally we point out that the NCF indicator generates a checkerboard-like pattern. This is despite the smooth distribution of the error indicator computed from the starting p-uniform solution.

IV.B. Steady Euler flow past a cylinder

The second test case considered is the inviscid flow past a cylinder at Mach number 0.3. Slip boundary conditions are imposed on the wall and non-reflecting boundary conditions are specified on the external boundary. The mesh employed is a O-type 4 th -order mesh generated using Gmsh. [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] This mesh has a radius R = 400D, with a number of elements N r = 20 and N θ = 26 in the radial and azimuthal directions respectively. In the radial direction, the element size changes following a geometric progression of ratio 1.2 up to r = 4D, and 1.6 up to the domain boundary. As in the previous test case the p-adaptive algorithm is applied by starting from uniform p = 1 and limiting the maximum local degree to 4.

Due to the fact that a geometric progression is employed in the radial direction, there is a significant variation in element size for the employed mesh.

In this case, the application of the VMS and residual-based indicators leads to extremely poor results, selecting for refinement the largest elements even in the presence of low average values of the error. This is possibly caused by the direct dependency to the element size in Eq. ( 14) and Eq. ( 27). This is not the case for the other considered indicators which all employ some form of normalization. This is the reason why we consider in the following the normalized version of the VMS indicator given in Eq. ( 14) :

η K . . = K (ρv) h,p -(ρv) h,p-1 2 dV |K| 1 2 ( 39
)
where |K| is the volume of the element K and the residual-based indicator is also normalized by dividing Eq. ( 27) by the characteristic size of the element h K :

η K . . = R(u h ) L 2 (K) + h -1 2 K r(u h ) L 2 (∂ K∩E i ) + h -1 2 K r b (u h ) L 2 (∂ K∩E b ) (40)
Figure 3 reports the entropy error corresponding to uniform and adaptive p-refinement. As in the previous test case, all refinement indicators provide improved performances with respect to uniform refinement. The number of DOFs necessary to achieve the same accuracy as uniform p = 4 is reduced by ~60%. Differences can be observed for the last steps of the refinement process, with the NCF indicator yielding the poorest results. Figure 4 shows the distribution of the local polynomial degree across the domain obtained using the different refinement indicators. As in the previous section, we observe that the VMS and the spectral decay indicators provide very similar refinement patterns. As regards the NCF indicator, the checkerboard-like pattern observed in the previous test case (Fig. 2) is also visible here. Interestingly compared to the residual-based estimator, the residuum-NCF indicator refines over a larger region around the cylinder before increasing the local p values in its proximity. This is possibly the cause of the slight loss in performance observed in Fig. 3 at high number of DOFs.

IV.C. Steady laminar flow past a cylinder at Re = 40

Finally, we consider the laminar flow past a circular cylinder at Re = 40 and M = 0.1. For the considered Reynolds number the flow field is two-dimensional, symmetric and steady, thus allowing for computationally inexpensive comparisons on a laminar flow configuration. The flow pattern presents separation leading to the formation of two recirculation bubbles in the near wake that need to be correctly resolved in order to compute the drag coefficient accurately. This is the reason why this configuration is often employed for the validation of CFD codes.

Two-dimensional simulations are therefore performed by imposing no-slip boundary conditions at the wall and non-reflecting boundary conditions on the external boundary. The accuracy of the solution will mainly depend on three parameters: the spatial resolution, the extension of the computational domain and the accuracy in the representation of the cylinder geometry.

It has been shown that the size of the computational domain for unbounded flows is of critical importance as it is responsible for blockage effects and can distort phenomena generated in the internal domain convected through the artificial boundary. [START_REF] Gautier | A reference solution of the flow over a circular cylinder at Re= 40[END_REF][START_REF] Posdziech | A systematic approach to the numerical calculation of fundamental quantities of the two-dimensional flow over a circular cylinder[END_REF] A computational domain of radius R = 2000D is therefore employed here and has been verified to have an effect on the drag coefficient in the order of O(10 -5 ).

The mesh employed for the following simulations is a O-type 4 th -order mesh with N r = 18 and N θ = 28 elements in the radial and azimuthal direction respectively. In the radial direction the element size changes following a geometric progression of ratio 1.25 up to r = 6D and 1.65 up to the external boundary.

In order to reduce the effect of quadrature errors, all the simulations where performed by using a number of integration points q = p + 5 in each spatial direction. We have verified that, increasing this value has negligible effects on the obtained solutions.

For the simulations presented in this section, the adaptive p-refinement has been performed by starting from uniform polynomial degree p = 2 and limiting the maximum local polynomial degree to p = 6. By limiting the maximum polynomial degree to p = 6 we can expect the solution of adaptive refinement to reach a level of accuracy similar to that obtained from uniformly refined simulation at p = 6, albeit for a reduced number of DOFs. Figures 5 and6 show the evolution of the error in the drag coefficient and in the L 2 -norm using uniform and adaptive p-refinement with respect to a reference solution obtained with uniform p = 7 (C D ref = 1.50219). As the uniform p = 6 solution produces an error in the drag coefficient of the order of 10 -5 , any further reduction of this error with p-adaptive refinement can be explained as small oscillations around the C D value obtained from the uniform p = 6 solution.

All refinement indicators are able to reduce by ~75% the number of DOFs necessary to achieve the same accuracy as that provided by uniform p = 6.

Only minor differences can be observed between different refinement indicators in Fig. 5 and6. If we consider the error in the drag coefficient, we cab observe from Fig. 5b that the residual-based and residuum-NCF refinement indicators show slightly better performance. As regards the discretization error in Fig. 6, the residuum-NCF refinement indicator and the VMS estimator appear to lead to the fastest convergence, with the residuum-NCF indicator performing best. However, after the initial steps, similar results are produced overall by the other refinement indicators.

Figures 7 and8 show the distribution of local p-refinement levels corresponding to the last iteration of the adaptive refinement procedure for different estimators. Note that, as the maximum marking criterion was employed, the number of additional DOFs introduced at each refinement iteration could not be controlled. Therefore, the results presented correspond to a similar but not exactly the same number of DOFs and number of iterations. While some differences can be observed, the distributions of local polynomial degree are fairly similar for all the considered refinement indicators. All the indicators appear to select for refinement an approximately circular region around the cylinder and the shear region extending up to the outer boundary (see Fig. 8). Moreover, all refinement indicators select for further refinement the regions right in front of the cylinder and those located at an angle of 90 • , before the separation of the flow which is located at 125 • . [START_REF] Gautier | A reference solution of the flow over a circular cylinder at Re= 40[END_REF] One of the main differences that can be observed is the disparity in the refinement levels obtained in the recirculation region and in proximity of the cylinder base. Among all refinement indicators, the VMS and spectral decay 10 indicators lead respectively to the lowest and highest refinement levels in this region. This difference in behaviour between the VMS and spectral decay indicators is somewhat expected. Indeed the spectral decay indicator normalizes the estimate of DE by the total 'energy' of the flow which means that it will tend to select for further refinement regions that are characterized at the same time by lower values of the error but also lower 'energy' (e. g. recirculation regions or stagnation points). This difference could be the cause of the slower convergence rate demonstrated by the spectral decay indicator in terms of the L 2 -norm of the error in the momentum density. In Fig. 8 the same refinement levels are shown for the far-wake region. As already observed, all refinement indicators select for refinement the far-wake region extending to the outer boundary. This is possibly because the large size of the elements requires an increase in polynomial degree in order to correctly capture the behaviour of the flow in this region. Therefore, this refinement reduces the L 2 -norm of the error in the momentum density even if it has a limited effect on the convergence of the drag coefficient.

Finally we note that the checkerboard-like pattern shown for the NCF indicator, in the Euler case, is not observed in Fig. 7 although it has been observed in the intermediate steps of the refinement process.

V. Implementation issues and computational cost

In section IV, the performances of the different refinement indicators have been analysed in terms of the reduction in the number of DOFs yielded by each estimator. It is indeed difficult to compare the computational time required by each approach due to the presence of several parameters that might affect the results obtained and therefore the derived conclusions. Firstly, we have observed that the choice of the marking parameter θ has a significant effect on the total computational time, while limited changes where observed in terms of the evolution of the error versus the number of DOFs of the adaptive refinement procedure. This is caused by the direct link between the computational time and the number of refinement steps which in turn depends on the marking strategy employed. Moreover, the computational cost of the refinement indicators strongly depends on the details of the implementation. Finally, direct comparison of the computational times for distributed memory computations will also be influenced by the possibility of achieving a good load balance by taking into account the local polynomial degree p.

In this section, we make some general observations that might be useful even when considering a different implementation from the one employed in our work. In particular, we consider the locality properties of the refinement indicators and their possible effect on parallel efficiency, the ease of implementation in existing CFD codes and the number of operations required by each indicator.

As regards locality, we can see from the definition of the different refinement indicators provided in section III that, the only two indicators which are fully local are the VMS and spectral decay indicators as their computation (see. Eq. ( 14) and Eq. ( 16)) only requires the knowledge of the solution inside the element. On the other hand, the 10 1.9 10 1.95 10 2 10 2.05 10 2.1 10 2.15 10 2.2 10 2.25 10 2. computation of the NCF indicator (given in Eq. ( 18)) requires knowledge of the projection at the interface of the solution from the direct neighbours, while the residual-based and the residuum-NCF indicators will require additional information to be exchanged for the computation of the numerical fluxes. In the framework of distributed memory computations this lack of locality will involve blocking message passing operations which may lead to a reduced parallel efficiency.

Considering now the number of operations involved in the computation of each indicator, we observe that the VMS and spectral decay indicators require the projection of the solution on the reduced order space S p-1 h and the computation of the integrals to compute the L 2 -norm. In the case of a modal DG approach based on hierarchical basis functions, the scale separation is readily available and the required integrals can be efficiently computed by employing information contained in the mass matrix (see Tumolo et al. [START_REF] Tumolo | A semi-implicit, semi-Lagrangian, p-adaptive discontinuous Galerkin method for the shallow water equations[END_REF] ). On the other hand, if a non-hierarchical basis is employed (e. g. nodal approach), the additional projection operation and the integration by quadrature will lead to higher computational cost. In this particular case, the NCF indicator would be more computationally efficient as it only requires the interpolation of the solution at the gauss points on each interface. As regards the residual-based and residuum-NCF indicators, their computation is more expensive as compared to the other indicators. Indeed, the evaluation of the former requires the computation and integration over the elements and faces of the element residuals, numerical fluxes and the convective and diffusive fluxes obtained from the internal reconstruction F (u + h , ∇u + h ). Finally, we expect the residuum-NCF to be the most expensive among the considered indicators, regardless of the particular CFD solver considered. This is because, the contribution of the nonconformity error is more expensive to compute than the NCF indicator alone as it requires the integration of the jump terms that appear in the expression of the NCF indicator. The evaluation of the residuum-error contribution, on the other hand, requires the computation of the discrete residuals N h (u h , φ i ) with φ i being the basis functions corresponding to the space S p+1 h . To this we have to add the cost of the Lagrange multiplier method involved in the computation of this indicator which requires the solution of a linear system of size (p k + 2) d for each element.

Considering now the ease of implementation, it is interesting to note that the complexity of computing the residuals corresponding to the refined polynomial space S p+1 h will depend on the numerical implementation of the CFD method. In the case of the Aghora solver, [START_REF] Renac | Aghora: a high-order DG solver for turbulent flow simulations[END_REF] this required the computation and storage of the basis functions for S p+1 h as well as their derivatives on the corresponding set of quadrature points. In general, important modifications of the employed CFD solver might be required to achieve a computationally efficient implementation of this method. This is not the case for the VMS, spectral decay and NCF indicators which can be easily implemented even as external post-processing tools. On the other hand, the definition of the different terms involved in the definition of the residual-based indicator (Eq. ( 22), Eq. ( 23) and Eq. ( 24)) depends on the set of equations being solved and the discretization employed.

These observations might be of limited interest in the context of steady-flow simulations, since the computation of a steady solution is likely to be more time consuming than the refinement procedure. However, these considerations will be of primary importance when extending these techniques to unsteady-flow simulations for which the refinement procedure is carried out over the course of the simulation.

VI. Conclusion

In this article we have assessed the performance of various refinement indicators for the implementation of an adaptive p-refinement algorithm in the context of a DG method. Five refinement indicators of different types, namely, feature-based, discretization and residual error-based have been employed for the simulation of inviscid and laminar flow configurations. The choice of the error estimators considered in this work has been made based on their potential suitability for the simulation of turbulent flows by LES. The goal-oriented approach to error estimation has therefore been discarded due to its high computational cost when dealing with unsteady flow simulations.

Three flow configurations have been considered: the inviscid flow over a Gaussian bump at M = 0.5, the inviscid flow past a circular cylinder at M = 0.3 and the laminar flow past a circular cylinder at M = 0.1 and Re = 40. The advantage of adaptive refinement as compared to uniform refinement in terms of number of DOFs has been clearly observed. Indeed, for all the configurations considered, the number of degrees of freedom necessary to achieve a prescribed level of accuracy was reduced by 60% to 75% thanks to the use of the p-adaptive algorithm.

Minor advantages with respect to other estimators have been observed from the application of the residual-based indicator to Euler flow simulations and from the application of the residuum-NCF indicator to the laminar flow simulations. However, similar results have been observed for the VMS, spectral decay and NCF indicators both in terms of the refinement pattern and convergence history. This is interesting if we consider the lack in their formulation of any specific terms related to the set of equations being solved and the discretization being employed. Moreover these estimators can be easily implemented as post-processing tools and require a reduced computational cost as we have seen from section V. In particular, the VMS and spectral decay indicators can be evaluated with very limited computational cost when a modal DG method based on hierarchical basis functions is employed. In contrast the computation of the residual-based and residuum-NCF indicators might require a considerable effort in order to obtain an efficient implementation well adapted to the particular solver at hand.

With regards to the VMS and spectral decay indicators, the former demonstrated similar if not superior performances as compared to the latter in each case considered. An advantage of the VMS over the spectral decay indicator is its lower computational cost regardless of the implementation as the normalization term at the denominator of Eq. ( 16) does not need to be computed.

Finally we have observed that the NCF indicator generates a checkerboard-like pattern which appears to depend on the refinement history. This pattern could potentially damage the convergence of the method as it introduces an irregularity in the spatial discretization.

All these results are of interest as they provide guidelines for the application of p-refinement to unsteady flow simulations. Indeed, for such applications, the computational cost of the refinement indicators will be a fundamental parameter in the choice of the adaptive refinement strategy. Other important issues that must be considered when designing adaptive refinement methods for unsteady flow computations will include the definition of a suitable marking criterion, refinement and derefinement thresholds and the possibility for the user to define relevant values for these a priori. The extension of this comparison to h and hp-refinement is the subject of current research.

15

 15

Figure 1 :

 1 Figure 1: Euler flow over a Gaussian bump at M = 0.5 : evolution of global entropy error Eq. (38) under uniform and adaptive refinement.

Figure 2 :Figure 3 :

 23 Figure 2: Euler flow over a Gaussian bump at M = 0.5 : map of local polynomial degrees obtained based on different refinement indicators.

Figure 4 :

 4 Figure 4: Euler flow past a cylinder at M = 0.3 : map of local polynomial degrees obtained based on different refinement indicators .

Figure 5 :

 5 Figure 5: Evolution of the drag coefficient and corresponding error under uniform and adaptive p-refinement for the flow past a cylinder at Re = 40 and M = 0.1 .

3 10Figure 6 :

 36 Figure 6: L 2 -norm of the error in the momentum density under uniform and adaptive p-refinement for the flow past a cylinder at Re = 40 and M = 0.1 .

Figure 7 :

 7 Figure 7: Local polynomial orders obtained with different refinement indicators for the laminar flow past a cylinder at Re = 40 and M = 0.1 .

  (a) VMS indicator (b) Spectral decay indicator (c) NCF indicator (d) Residual-based indicator (e) Residuum-NCF indicator

Figure 8 :

 8 Figure 8: Local polynomial orders obtained with different refinement indicators for the laminar flow past a cylinder at Re = 40 and M = 0.1 .

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675008. This work was performed using HPC resources from GENCI (Grant 2017-A0022A10129). We also thank Mathieu Lorteau, Emeric Martin, Florent Renac and Marie-Claire Le Pape for the valuable discussions.