domain: spi phys

text:
Milan Rozel
email: milan.rozel@onera.fr

Nicolas Gonçalves

Alec Reygrobellet

Pierre Bruneel

Construction of a Radar Cross-Section Database Using ADS-B Data from the OpenSky Network †

Keywords: database: radar, ADS-B, radar cross section

Building a labeled database usually requires extensive work when the labeling is not automated. This paper presents a methodology to build a labeled database of aircraft radar signature (specifically radar cross section) using the OpenSky Network data. This method relies on a radar system steered by the live API to automatically select, track and measure target aircraft according to custom user rules. The OpenSky ADS-B database allows for an efficient processing of the radar data, as it provides ranging and speed information to locate the region of interest inside the radar data. In addition, the ADS-B and Mode-S databases give an estimate of the attitude of the aircraft, which is paramount for signature measurement. We also compare the ADS-B data with the radar data.

Introduction

Signature is a set of elements that makes it possible to characterize the presence, the type or the identity of an object by a sensor. Radar signature refers to radio-frequency scattering mechanisms on targets. In the military field, the lower the signature, the stealthier the aircraft. Radar cross section (RCS) of a target quantitatively represents the amount of energy scattered back to a radar, and depends on the aspect angle of the target as seen by the radar. In low resolution, the target is considered a single-point object, whereas high resolution allow target imaging and the localization of targets' hot spots [START_REF] Titin-Schnaider | Full-polarimetric analysis of MERIC air targets data[END_REF]. Gathering a large amount of signature data inside a database is necessary for automatic target recognition using artificial intelligence or deep learning [START_REF] Oh | Micro-Doppler Mini-UAV Classification Using Empirical-Mode Decomposition Features[END_REF][START_REF] Gusland | Open Radar Initiative: Large Scale Dataset for Benchmarking of micro-Doppler Recognition Algorithms[END_REF]. This database might also help to benchmark radar signature simulation software, since few measured signature data are available. In fact, indoor aircraft-signature measurements [START_REF] Gürel | Validation through comparison: Measurement and calculation of the bistatic radar cross section of a stealth target[END_REF] require either great dimension anechoic chambers (for scale 1 aircrafts) or measurements at very high frequencies (typically above 100 GHz for reduced scale aircraft models), which are challenging for model positioning.

We propose a method to build a database of aircraft radar signature from outdoor measurements, using a radar developped by ONERA in Palaiseau called HYCAM [START_REF] Brouard | A new S band surface radar testbed[END_REF] and ADS-B data from the OpenSky network [START_REF] Schäfer | Bringing Up OpenSky: A Large-scale ADS-B Sensor Network for Research[END_REF]. The database gathers both low-resolution and high-resolution measurement data; in this article, we only present low-resolution measurements. The radar works in designation mode: the tracking feature uses GPS position communicated by the aircraft via ADS-B (a set of information broadcasted by most commercial aircraft to indicate their position) in real time, thanks to the OpenSky Network real-time API (an API is a communication framework for to interface different comuters). Moreover, the retrieved data allow the database documentation and labeling (ICAO24 identification number, aircraft type, engine type, etc.). Finally, the ADS-B data provide information about the groundtruth of the target (time, speed, position, and heading), which is necessary for radar signature analysis. The method provides a manual operation mode and an automatic operation mode. In manual mode, radar operators carry out the measurements, i.e., they manually choose which plane to measure. In automatic mode, the software can select the targets according to various selection rules, which can help to enrich the database according to custom user goals (diversification or consolidation).

We validate the proposed methodology by comparing the ADS-B data with actual radar measurement and assess that the quality of speed and range data yielded by ADS-B is good enough to lift ambiguity and to automatically focus on the region of interest inside radar data. The steering of the radar and the inertial data of the aircraft allow assessing the aircraft attitude. We present the methodology for the automatic measurement and the building of the database, as well as several radar measurements of planes labeled by the OpenSky API.

Real-Time Radar Architecture

The HYCAM radar is an S-band radar with a center frequency around 3 GHz, designed in 2008 for radar signature measurements purposes and to carry out software-defined radar studies. The antennas are mounted on a 360°az / 45°elevation positioner, allowing target tracking within this field of regard. It has a full polarimetric capacity, and can use arbitrary waveforms (with either narrow and wide bandwidth). In this article, we use pulsed chirps, as described in Section 3.2. Further information about HYCAM can be found in [START_REF] Brouard | A new S band surface radar testbed[END_REF]. The GPS coordinates of HYCAM are 48°42 47 N, 2°13 51 E, 175 m above sea level. The distances from the planes to the radar and the orientation of the turret use these coordinates as we detail in Section 3.3.

Target Tracking and Radar Workflow

As aforementioned, OpenSky real-time data are treated to send instructions to the HYCAM positioner and track aircrafts. In this article, the tracking algorithm linearly extrapolates the ADS-B position of the planes to fill the gaps between successive packets of data. The turret receives azimuth and elevation instructions to orient toward the plane. A radar with an automatic tracking system may autonomously perform the pursuit of the target after this choice, but this is not the case in this work. As the aperture of the radar is relatively important (10°× 4°), the positioning errors should not impact the observation of the plane. A sanity check of the radar data and turret log can showcase whether there are missing data due to poor estimation of the position. The start time (GPS time) of radar acquisitions are recorded. Meanwhile, OpenSky data are also timestamped using Unix time. Both are synchronized during processing.

The data from the OpenSky API are used not only to perform tracking of the target but also for further processing of the data. In particular, they are recorded jointly with the acquisition and may be of use to estimate the aircraft position and to automatically document the entries of the database. In the following, "ADS-B quantity" refers to the value of a quantity given or computed from ADS-B data. It can be the distance from the target to the radar (ADS-B range) or the radial speed from the target to the radar (ADS-B speed). On the contrary, a "radar quantity" is estimated from the radar data.

Target Selection

The OpenSky API yields the current state of the aerial traffic, and this helps to choose the plane to be measured according to a set of criteria. Using both inertial properties (position, heading, etc.) and metadata (model, ICAO24 unique identification number of the airframe) of each plane, we build up a methodology to optimally choose a plane to be measured by the radar. We denote by (r p , f p ) respectively the inertial properties and the metadata about a plane p, and choose a utility function U(r p , f p ). We propose a basic algorithm to optimally choose the plane p to be measured:

•

At each time step, compute the utilities U(r p , f p ) for all planes p in a given area around the radar.

• If no plane is currently measured, then assess whether the maximum utility exceeds a certain threshold U 0 . If this is true, then the plane with the maximal utility becomes tracked by the radar. • If a plane is under observation and its utility goes lower than a threshold U < , the measurement is stopped. We introduce a hysteresis U < < U 0 to avoid that the acquisition immediately stops in case of utility low stability. • If a plane p 1 with utility U 1 is under observation and the plane p 2 with utility U 2 has an utility which is substantially greater than U 1 , we allow for the system to switch from p 1 to p 2 . Typically U 2 > U 1 + U > , with U > > 0 to also introduce a hysteresis, to avoid that the system indefinitely switches between two planes with a close utility.

•

This algorithm is repeated until a stop command related to the size of the data gathered during the campaign is reached in order to avoid an overflow of the storage device.

This framework is very generic and can adapt to different sensor technologies. The utility function can implement exclusions criteria or hard constraints to match system limitations (turret mechanical stop). It can favor specific models or specific plane attitudes. This utility function can even depend on the state of the database: one can seek new plane models, or consolidate previous measurements by choosing similar conditions. For now on, the utility function we use is decreasing with distance to guarantee a good signal-to-noise ratio (SNR). The choices of the utility function and of the policy of selection can also add priority mechanisms, fuzzy logic or other advanced concepts such as those mentioned in [START_REF] Charlish | Autonomous Agents for Multi-Function Radar Resource Management[END_REF]. We use mechanical steering in our work; the utility and policy can be adapted to harness the assets of electronic beamsteering, such as the ability to dynamically switch between different areas.

Processing Workflow

The ADS-B data simplify the radar processing, and provide an estimation of the attitude of the measured aircraft to correlate it with the radar measurements.

Radar Equation

The radar equation [START_REF] Skolnik | Radar Handbook[END_REF] gives the power of the electromagnetic wave, which is reflected back by a target and received by the radar:

P r = P t G t 1 R 2 σ 1 R 2 G r λ 2 (4π) 3 (1)
The different terms of this equation represent all the phenomena from emission to reception of the electromagnetic wave:

• P r is the received power, and P t the power emitted by the radar. This equation is a scattering equation and states that the power reflected is proportional to the power emitted, • R the distance between the target and the radar. The two 1 R 2 factors are two spherical wave decays due to the two paths between the radar and the target (back and forth); • G t , G r are the gains of the antenna in emission and reception. They can be different if the reception antenna is different from the emission antenna; • σ is called the radar cross section (RCS) of the target, and represents the effective area of the target as seen by the radar; • λ is the wavelength of the transmitted waveform. The constants in this expression (gain, emitted power, and wavelength) can be compensated by a calibration procedure that we neither detail nor apply in the data shown in this work. In this paper, we focus on the time evolution of the RCS of the target and compare the ADS-B range and speed with the radar range and speed. We compensate for the 1/R 4 term to have consistent RCS values over a single acquisition, following the procedure detailed in Section 3.3.

The radar equation actually applies to a single point target, whereas real targets are composed of several scattering centers with different radar range and radial speeds. Depending on the radar range and speed resolutions, the scattering centers might be separated by radar processing. We call aircraft range extension the maximum radar range difference between two scattering centers of the same target. The aircraft range extension is of the order of the largest dimension of the plane. Speed extension refers to the maximum radial speed difference scattered back from the aircraft to the radar.

Radar Parameters

In this work, the radar periodically sends chirp pulses [START_REF] Klauder | The theory and design of chirp radars[END_REF], with a period equal to the pulse repetition interval (PRI). Distance aliasing happens when the distance D from the target to the radar is greater than the so-called ambiguous range D amb = cPRI 2 (where c is the speed of light). The target is thought to be at a distance equal to D mod D amb , without additional processing. When the radar records and simultaneously transmits its pulses, a portion of the transmitted power is received directly from the transmit antenna to the receive antenna and can corrupt the signature measurement so that the data in this time domain are often discarded. This is called coupling and is a blind spot of the radar and is related to the duration τ of the pulse sent d blind = cτ 2 . The distance resolution of the radar is given by the bandwidth B of the waveform δd = c 2B . The radial velocity of the target is measured by the Doppler effect experienced at the reflection of the pulse on the moving target. It is determined by a Fourier transform of the signal along the number of pulse axis. There is also a speed ambiguity due to a stroboscopic effect; the Doppler frequency is known modulo the pulse repetition frequency PRF

Ambiguity Lifting

Knowing the actual radar range is necessary to estimate the RCS of aircraft in order to compensate the 1 R 4 term of Equation (1). However, radar range and speed measurements can be ambiguous as it is sketched in Section 3.2. Techniques based on the Chinese remainder theorem can lift range and speed ambiguities [START_REF] Li | A robust Chinese remainder theorem with its applications in moving target Doppler estimation[END_REF], but they require specific waveform schemes. ADS-B positioning (latitude, longitude and altitude) and speed data (ground speed and vertical rate) allow the computation of ADS-B range and ADS-B speed to help lifting ambiguities. Let O be the radar location, and A the aircraft location. ADS-B range is equal to OA . The track angle, ground speed and vertical rate provided by OpenSky data are used to estimate the speed vector v of the aircraft. Then, the ADS-B radial speed is given by OA OA . v. A positive speed means that the aircraft flies away from the radar.

Figure 1 shows a Doppler-range map, representing the amplitude of the radar signal as a function of the radial velocity of the target along the x-axis and the range of the target in the y-axis. The horizontal line near the top and bottom of the figure are due to the coupling with the emission of the radar. The vertical line in the center of the x-axis represents the return of the ground (clutter), which is indeed at zero speed. The cloud of nonzero data at around -5 m.s -1 speed and from 0 to 15 km range is the radar return of the rain captured during the measurement. If the signature of the target dwells in the coupling, the clutter or the rain cells, the data can be corrupted. In this map, the range and speed of the aircraft is ambiguous: the radar measurements state that the aircraft is located at 4.95 km from the radar, at a radar speed of -105 m.s -1 , whereas the ADS-B range is 35.0 km, and the ADS-B speed is 135 m.s -1 . ADS-B range and speed make ambiguity lifting straightforward.

Data Reduction

The radar continuously records data, of which the aircraft signature only occupies a fraction, which we call the region of interest of radar data, as shown by the inset in Figure 1. The signature data are indeed located in a red area around the ADS-B position of the aircraft in Doppler-range maps, which is only a small portion of the total range swath (D amb = 30 km) and speed swath (v amb = 240 m.s -1 ). Hence, the signature data are cropped thanks to OpenSky data. An interpolation step is performed around the maximum value in the region of interest to compute the radar range and speed with a sub-cell accuracy.

Figure 2a,b show the difference between radar and ADS-B range and speed, for the aircraft with ICAO24 4b1a2a flying the 7 October 2022 between 08:26:36 UTC and 08:27:46 UTC, computed with the method described in Section 3.3. In this example, the standard deviation of the difference between radar and ADS-B range is 8.5 m. The difference between radar and ADS-B speed remains between ±0.25 m.s -1 . We observe steps on the speed difference chart, which can be explained by the quantization of the ground speed (and vertical rate) data provided by OpenSky, as illustrated in Figure 2b,c.

In addition to this quantization, the remaining differences in range and speed may be due to a limited accuracy of radar range and speed measurements in case of low SNR. Moreover, ADS-B provides the GPS position and speed of a fixed point of the aircraft, whereas the radar measures the range and speed of scattering centers whose location can vary with aspect and differ from the GPS point.

Yet, the remainder between radar measurement and OpenSky data is small, and makes it possible to automatically keep the relevant data for RCS measurement and analysis. To protect from the remaining errors, the method described in this article keeps a range swath of ±300 m around the ADS-B range, and a speed swath of ±4 m.s -1 around the ADS-B speed. In Figure 1, the swath of interest is inside the red frame. This significantly reduces the amount of data stored in the database.

RCS Computation

During the data reduction step presented in Section 3.4, the maximum amplitude in the signature data is measured, which is corrected by the unambiguous radar range (following Equation ( 1)), leaving only constant factors multiplying the RCS. A calibration process remains mandatory to yield actual RCS values. Figure 3 shows the computed RCS as a function of time for the plane with ICAO24 46316e flying the 12 October 2022 between 12:05:59 UTC and 12:07:26 UTC. As a reminder, the signatures shown in this article are not calibrated-here, the plotted RCS is centered around its mean value. The signature of the aircraft fluctuates over almost 35 dB, due to changes of the aspect angle of the aircraft. The RCS is maximal in the beginning, when the speed vector of the aircraft and the radar steering vector are almost orthogonal. Such measurements are of great interest to study and build radar signature statistical models.

Aspect Angle Computation

A surveillance radar shall detect aircraft regardless of their relative orientation to the radar, and this orientation continuously evolves along the trajectory. As the RCS of an aircraft strongly varies with aspect angles (see Figure 3), the specifications of a radar must account for this fluctuation to guarantee performance of detection. Thus, signature analysis studies radar signatures as a function of aspect angle. Estimating the attitude of the aircraft via ADS-B data tags RCS data with aspect estimates in order to map the RCS vs. aspect angle function.

Figure 4 illustrates the two orthonormal bases we use in this work: one for the radar B u = ( u x , u y , u z ) and one for the aircraft B v = ( v x , v y , v z ). We denote by R u,v the transfer matrix corresponding to basis B u,v . The coordinates of the vectors are computed in a common xyz basis using the GPS and ECEF (Earth Centered Earth Fixed, a cartesian coordinate system for locating objects) coordinates [START_REF] Zhou | Sensor alignment with earth-centered earth-fixed (ECEF) coordinate system[END_REF]. u x is the radar steering vector, u y the vector along the horizontal of the antenna and u z the vector oriented along the vertical of the antenna. The aircraft is described by three principal axes, defining three rotation angles: roll, pitch and yaw. We compute the basis B v using the position of the aircraft from ADS-B data and assuming that roll, pitch and yaw values are known. Any vector v in B v can be written in B u according to Equation [START_REF] Oh | Micro-Doppler Mini-UAV Classification Using Empirical-Mode Decomposition Features[END_REF].

u = (R u • R v ) • v = M • v (2)
M can be decomposed as the product of three rotation matrices defined by the three Euler angles, yaw (α), pitch (β) and roll (γ):

M = M z (α) • M y (β) • M x (γ)
, where • is the matrix product. These are the apparent aspect angles of the aircraft seen from the radar. Actually, the yaw, pitch and roll of the aircraft are not completely determined from OpenSky data. We assume that v x is co-linear to the speed vector of the aircraft, meaning that the pitch is the steep of the trajectory, the yaw is the track of the trajectory and that the roll is zero. In fact, the three angles differ in windy conditions or during maneuvers. This estimation can be further refined by using external data from extended mode S as discussed in Section 4.2. The estimated aspect angles are stored in the database for the sake of signature analysis.

Database Presentation

This part gives an overview of the two OpenSky Network databases that allows retrieving the flight data, and presents the data stored in our radar signature database.

ADS-B Sources

The methodology of this article relies on ADS-B data for real time tasks (aircraft selection, tracking, labels) and post-processing tasks (signal processing and RCS computation). Real-time ADS-B data are yielded by the real-time OpenSky API [START_REF] Schäfer | Bringing Up OpenSky: A Large-scale ADS-B Sensor Network for Research[END_REF], which is recorded live. It is compared with actual turret positions to assess whether the tracking feature worked correctly. If not, the data are not stored in the database.

The Impala state_vectors_data4 database [START_REF]A Quick Guide To OpenSky's Impala Shell[END_REF] is the source for the processing proposed in Section 3. It provides aircraft positioning data sampled at around 1 Hz, ten times quicker than our requests to the real-time API. In the case of missing Impala data, positions from the real-time database remain available. These additional position data allow to estimate the error between the extrapolated and actual positions.

Mode-S Data

The OpenSky Network also provides access to the Impala rollcall_replies_data4 database, which contains the Mode-S messages. Some of these messages are known as Comm-B Data Selector (abbreviated as BDS) and contain information on the attitude of the aircraft. BDS 50 messages contain the roll information and BDS 60 messages contain the heading information, which allows for better knowledge of the attitude of the aircraft. The type of the message must be inferred and parsed accordingly to their type, with a method presented in [START_REF] Sun | Megahertz Riddle: A Guide to Decoding Mode S and ADS-B Signals[END_REF] implemented in the pyModeS python library [START_REF] Sun | Decoding mode-s surveillance data for open air transportation research[END_REF]. If present, these data improve the knowledge of the attitude of the plane.

Information Stored in the Database

The database contains measurement descriptors used to perform queries and to facilitate signature analysis. It provides direct aircraft signature information as well as Doppler-range maps so that users can build their own processing chain. Therefore, the database stores low level data: From this low-level data, custom treatments allow users to build signature data. This makes signature analysis as a function of time or aspect angle possible using the estimated attitude. The database can be queried by the ICAO24 identification number, aircraft type (when available) and date of the measurement.

Limitations Due to the Data Sources

ADS-B data can present several sources of errors and are not entirely reliable [START_REF] Schäfer | Assessing the integrity of crowdsourced mode S and ADS-B data[END_REF]. In addition, the extrapolation algorithm used in this work poorly behaves when the plane is turning such that turret steering can be improved. Several sanity checks can be performed to gauge the quality of the gathered data. As sketched in Section 4.1, the forecast plane position is compared with the consolidated ADS-B logs of the Impala database. This allows gauging the error in the extrapolating scheme, and can be useful to test tracking algorithms for planes. A bad steering of the turret toward the target induces a total loss of signal. Good signature data feature high SNR, which is impacted by the range and position inside the lobe.

Uncertainties on the inertial quantities also affect the estimation of the attitude of the plane. The true angle of observation can deviate from the angle of observation of the turret due to a bad estimation of the position. Yaw and pitch can only be approximated by the direction of the velocity of the plane as computed from ADS-B position data. Roll data and consolidated yaw data can be available a posteriori with the mode-S messages. In the absence of these data, uncertainty on the axis of the plane might be important, especially during turns, but can be bounded. A detailed computation of these uncertainties is beyond the scope of this work, but can entirely be performed offline.

Discussion

This article highlights the possibility to automatically build a labeled RCS database of commercial aircraft by using the OpenSky live API. Using the plane coordinates transmitted by ADS-B, an extrapolation of the trajectories allows tracking planes by the mechanical turret of the radar to measure their radar signature. The OpenSky data are also used to locate the region of interest inside the radar data and to deal with ambiguities and distance correction. The ADS-B data from OpenSky match the required accuracy to track planes with our radar, as shown by a comparison between ADS-B and radar quantities. The state of aerial traffic yielded by OpenSky is also of use to build an algorithm able to automatically choose a plane to be measured, using a framework of utility coupled to simple policies to choose planes matching an optimal criterion. This methodology is a first step for building a database of in situ radar signatures. An estimation of the attitude of the plane allows to map the signature with an estimation of the orientation of the plane with respect to the antenna.

This methodology is versatile and can adapt to other systems. The building blocks are a positioner (mechanical turret and beam steering) and a sensor (radar, infrared, optical, and camera). This framework allows automatically gathering measurements to build plane-related databases and can support the design of surveillance and tracking algorithms. For databases, automatic measurements with labels are necessary for many advanced dataprocessing applications (machine learning, recognition, and statistical analysis). Further work around the automatic decision framework might help to specifically tailor the measurements to put in the database in order to answer a specific question (focus on specific models or aspect angles. . . ). A comparison of extrapolation algorithms can be based on the secondary data of trajectory yielded by OpenSky, and background checks of the difference between steering data and ADS-B data can help tuning parameters of a control process. A ground truth might help to design surveillance algorithms using a steering beam to optimize the share of sensor time between search and track. The offline plane trajectory database also helps to specify realist situations.

  = 1 PRI . The Doppler frequency, f D , relates to the velocity of the target v through f D / f c = 2v/c, with f c , the central frequency of the signal. Thus the radial velocity of the target is known modulo the speed ambiguity v amb = cPRF 2 f c . In the following, we use PRI = 200 µs, τ = 10 µs, f c = 3.15 GHz and B = 2 MHz, giving a distance ambiguity of D amb = 30 km, a blind zone of d blind = 1.5 km, a range resolution of δd = 75 m and an ambiguous speed v amb = 240 m.s -1 .

Figure 1 .

 1 Figure 1. Doppler-range map (ICAO identifier: 406c43, date: 14 October 2022 08h41:01 UTC).

Figure 2 .

 2 Figure 2. (a) Difference between radar range and ADS-B range, (b) difference between radar speed and ADS-B speed, (c) ground speed evolution.

Figure 3 .

 3 Figure 3. RCS evolution of aircaft with ICAO 46316e the 12 October 2022 between 12:05:59 UTC and 12:07:26 UTC.

Figure 4 .

 4 Figure 4. Radar and aircraft referentials.

Funding: This research was funded by ONERA.

Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.

Data Availability Statement: Data concerning the positioning of the planes mentioned can be found in the Impala database https://opensky-network.org/data/impala, accessed on 19 October 2022.

Conflicts of Interest:

The authors declare no conflict of interest.


domain: phys spi

text:
Valentin Carlier
email: vcarlier@ipp.mpg.de

Florent Renac
email: florent.renac@onera.fr

INVARIANT DOMAIN PRESERVING HIGH-ORDER SPECTRAL DISCONTINUOUS APPROXIMATIONS OF HYPERBOLIC SYSTEMS

Keywords: Hyperbolic systems, Convex invariant set, Limiting, Spectral discontinuous method, Discontinuous Galerkin method AMS subject classifications. 65M12, 65M70, 35L65

come

We are interested here in high-order numerical solutions to hyperbolic systems of conservation laws (1.1)

∂ t u + ∇ • f (u) = 0, in D × (0, ∞), u(x, 0) = u 0 (x), in D,
where u(x, t) represents the vector of conserved variables with values in the set of states Ω a ⊂ R m which is assumed to be convex. The flux tensor f = (f 1 , . . . , f d ) :

Ω a ∋ u → f (u) ∈ R m×d is assumed to be smooth.

Solutions to (1.1) may develop discontinuities in finite time even if the initial data is smooth, therefore the equations have to be understood in the sense of distributions.

Nevertheless, in this setting we lose uniqueness of the solution and (1.1) must be supplemented with further admissibility criteria. We here focus on entropy inequalities on some twice differentiable strictly convex function η : Ω a → R associated with a smooth entropy flux q : Ω a → R d satisfying

(1.2) η ′ (u) ⊤ f ′ i (u) = q ′ i (u) ⊤ ∀u ∈ Ω a , 1 ≤ i ≤ d.
A weak solution to (1.1) is called an entropy weak solution if for every entropy pair of (1.1) we have

(1.3) ∂η(u) ∂t + ∇ • q(u) ≤ 0,
in the sense of distributions. Classical (smooth) solutions respect this condition with an equality. The inequality for discontinuous solutions may be derived from a vanishing viscosity or semilinear relaxation approximations to (1.1), the regularizing effect allowing to have a smooth unique solution in this case. Then, under some structural assumptions on (1.1), regularized approximations converge to an entropy measure valued solution to (1.1) and (1.3) [?, ?]. These results are in particular based on the existence of convex invariant domains B ⊂ Ω a for (1.1): if u is in B, then it remains in B almost everywhere in D × (0, ∞) [?, ?, ?, ?]. This property generalizes the notion of maximum principle for scalar equations. Numerical methods keeping this property at the discrete level are called invariant domain preserving (IDP).

We are here interested in the approximation of (1.1) using high-order discontinuous spectral methods (see, e.g., [?, ?, ?, ?, ?] and references therein) where the solution to ( The main contribution of this work is the derivation of a CFL condition on the time step (in Theorem 3.5) for a large class of discretization methods on general unstructured meshes with explicit time stepping, that guarantees that each cell-averaged solution is a convex combination of states lying in invariant domains, provided those domains are preserved by the numerical flux at mesh interfaces. As a consequence, the cell-averaged scheme inherits some robustness and stability properties of the firstorder finite volume scheme with the same numerical flux. The CFL condition is general in the sense that it only requires the scheme to be conservative, to satisfy some geometric conservation law, and to allow the existence of a quadrature rule to evaluate the cell-averaged solution that includes the traces of the numerical solution used to evaluate the numerical fluxes in the scheme. This CFL condition relies on two other contributions of this work. The first one concerns the existence of a state in the invariant domains which satisfies a trivial flux balance over the faces of the each element (Lemma 3.4). We call this state the pseudo-equilibrium state and then use tricks from [?] to expand the cell-averaged approximate solution as a convex combination of states lying in the invariant domains. The second one concerns the construction of the quadrature rule which generalizes [?, Eq. (3.5)] on triangular grids to general curved polygonal and polyhedral elements.

The interest in this CFL condition is that it allows the use of limiters to make the full high-order approximate solution within the cell IDP. Linear scaling limiters [?, ?, ?] represent an important contribution in this direction. These limiters scale the cellwise approximate solution around its cell average to impose bounds on the solution such as positivity (i.e., with B = Ω a ), maximum and minimum principles, etc., while preserving conservation of the method. Such limiters therefore require the cell-averaged solution to be IDP and were limited to Cartesian and simplicial meshes for this reason. We also refer to the works in [?, ?] which generalize the CFL condition to unstructured quadrangular meshes with straight-sided elements. In contrast, the present CFL condition holds for a broad class of time explicit spectral discontinuous methods on general unstructured multidimensional meshes with possibly curved elements. For application purposes, we use the IDP states appearing in the convex combination defining the cell-averaged solution as local bounds which we impose by

This manuscript is for review purposes only.

limiting the solution in modal and nodal discontinuous Galerkin (DG) schemes. This strategy is closely related to convex limiting [?] based on first-order IDP approximations defining local bounds and then forcing the high-order approximation to satisfy these bounds through flux limiting [?, ?]. This latter approach has been applied to finite element approximations in [?] and DG spectral element method in [?] among others.

The paper is organized as follows. In section 2 we introduce the notion of invariant domain preserving approximate Riemann solvers (ARS). In section 3 we state and prove our theorem on the existence of the pseudo-equilibrium state and the CFL condition for high-order schemes, and then introduce an algorithm for their evaluations. in Ω a , it is convenient for the present analysis to consider the Riemann problem in the direction n:

(2.1)

∂ t u + ∂ x f (u) • n = 0, in R × (0, ∞), u(x, 0) = u L , if x < 0, u R , if x > 0.
We will integrate (2.1) over the space time slab [-h 2 , h 2 ] × [0, ∆t] with h > 0 and ∆t > 0 the space and time steps. We suppose here that all the Riemann problems we consider have a self-similar entropy weak solution W( x t ; u L , u R , n). Let introduce the self-similar variable ξ = x ∆t and assume that there exist σ L , σ R such that:

W(ξ; u L , u R , n) = u L for ξ < σ L and W(ξ; u L , u R , n) = u R for ξ > σ R .
We then define the maximum wave speed in (2.1) by

(2.2) |λ|(u L , u R , n) = max (|σ L |, |σ R |),
and for ∆t h |λ|(u L , u R , n) ≤ 1 2 , we define the average over the Riemann fan [?, ?]

(2.3) ū(u L , u R , n, ∆t h ) := 1 h h 2 -h 2 W x ∆t ; u L , u R , n dx = u L + u R 2 - ∆t h f (u R )-f (u L ) •n, which satisfies (2.4) η ū(u L , u R , n, ∆t h ) ≤ η(u L ) + η(u R ) 2 - ∆t h q(u R ) -q(u L ) • n.
We will also use the definition of invariant domain from [?]: a convex set B ⊂ Ω a is an invariant domain for (1.1) if for all u L and u R in B, we have

ū u L , u R , n, ∆t h ∈ B ∀ ∆t h |λ|(u L , u R , n) ≤ 1 2 .
2.2. Two-point numerical fluxes and approximate Riemann solvers.

The discretization of (1.1) will rely on two-point numerical fluxes [?, ?] for the approximation of f (u) • n and we assume them to be consistent and conservative:

(2.5) h(u, u, n) = f (u) • n, h(u L , u R , n) = -h(u R , u L , -n) ∀u, u L , u R ∈ Ω a ,
and Lipschitz continuous. We also define the notions of IDP and entropy stable twopoint flux in the following definitions.

Definition 2.1. A two-point flux is said to be invariant domain preserving (IDP)

for B an invariant domain if we have

u - ∆t h h(u, u R , n) -h(u L , u, n) ∈ B ∀ u L , u, u R ∈ B,
under the half CFL condition

(2.6) ∆t h max |λ|(u L , u, n), |λ|(u, u R , n) ≤ 1 2 .
Definition 2.2. A two-point flux is said to be entropy stable for the pair (η, q)

in (1.3) if, under the half CFL condition (2.6) and for all u L , u, u R in Ω a , we have

η u - ∆t h h(u, u R , n) -h(u L , u, n) ≤ η(u) - ∆t h Q(u, u R , n) -Q(u L , u, n)
for some consistent and conservative numerical flux Q(•, •, •).

We now introduce the notion of ARS and IDP ARS, which will be used to derive IDP two-point fluxes.

(2.7) 1 h h 2 -h 2 W a x ∆t ; u L , u R , n dx = u L + u R 2 - ∆t h f (u R ) -f (u L ) • n.
Using consistency in (2.5), and setting ξ = x ∆t , one defines a two-point flux from an ARS as

h W a (u L , u R , n) = f (u L ) • n - 0 -λ W a (ξ; u L , u R , n) -u L dξ (2.8a) = f (u R ) • n + λ 0 W a (ξ; u L , u R , n) -u R dξ, (2.8b)
where λ = |λ|(u L , u R , n). Both definitions are equivalent due to (2.7). We can now define the notion of IDP ARS.

Definition 2.4. An ARS W a (ξ; u L , u R , n) is IDP for an invariant domain B if we have 1 λ 0 -λ W a (ξ; u L , u R , n)dξ ∈ B, 1 λ λ 0 W a (ξ; u L , u R , n)dξ ∈ B ∀u L , u R ∈ B.
This manuscript is for review purposes only.

As a consequence

|λ|(u L , u R , n) ≤ 1, we have u L - ∆t h h W a (u L , u R , n) -f (u L ) • n ∈ B, (2.9a) u R - ∆t h f (u R ) • n -h W a (u L , u R , n) ∈ B. (2.9b)
Proof. Let consider (2.9a), a similar argument holds for (2.9b). From (2.8a), we have

u L - ∆t h h W a (u L , u R , n)-f (u L )•n = 1- ∆t h λ u L + ∆t h λ 1 λ 0 -λ W a (ξ, u L , u R , n)dξ
(2.10) u - ∆t h h W a (u, u R , n) -h W b (u L , u, n) ∈ B ∀u L , u, u R ∈ B,
under the half CFL condition (2.6).

Proof. We rewrite (2.10) as

1 2 u -2 ∆t h h W a (u, u R , n) -f (u) • n + 1 2 u -2 ∆t h f (u) • n -h W b (u, u R , n)
and apply Lemma 2.5 to both terms with 2 ∆t h λ ≤ 1.

The previous lemma with W b = W a also proves that the three-point scheme built from an IDP ARS is also IDP [?].

3. A general CFL condition. We here state and prove our main results on the existence of an explicit condition on the time step to ensure that the cell-averaged solution from a high-order spectral discontinuous scheme is IDP. In subsection 3.1

we clarify the schemes we are considering in this work. Our results are based on the existence of a pseudo-equilibrium state allowing a balance of the numerical fluxes at the faces of each element which is introduced in subsection 3.2 where we prove its existence. The main result (Theorem 3.5) on the CFL condition on the time step is given in subsection 3.3, while subsection 3.4 introduces a fast algorithm to evaluate the pseudo-equilibrium state, and so the time step. For application purposes a limiting strategy based on convex bounds is described in subsection 3.5.

3.1. Cell-averaged fully discrete scheme. We now describe the main properties of the numerical methods we are considering in this work. We consider here discretely conservative high-order approximations of (1.1). Without loss of generality, we use an explicit forward Euler discretization in time. High-order time integration

This manuscript is for review purposes only.

is then performed using strong-stability preserving Runge-Kutta methods [?] that are convex combinations of explicit first-order schemes in time and thus keep their stability properties. For the spatial discretization, the approximate solution u h (x, t) is defined locally over each element κ of the partition Ω h of the domain D in a local function cell space V p h (κ). By u

(n+1) h (•) = u h (•, t (n+1)
) we denote the solution at time t (n+1) = t (n) + ∆t (n) with t (0) = 0 and ∆t (n) > 0 the time step. The approximate solution is assumed to satisfy the following relation for the cell-averaged solution

⟨u h ⟩ κ : (3.1) ⟨u (n+1) h ⟩ κ = ⟨u (n) h ⟩ κ -∆t (n) N f k=1 s κ k h u - h (x κ k , t (n) ), u + h (x κ k , t (n) ), n κ k ∀κ ∈ Ω h ,
where the x κ k are some points on the faces f in ∂κ and u

± h (x κ k , t (n) ) = lim ϵ→0 + u h (x κ k ± ϵn κ k , t (n)
) denote evaluations of the traces of the solutions at x κ k (see Fig. 1). The

S κ := N f k=1 s κ k . u - h u + h n κ k x κ k y κ i • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Fig. 1: Notations for d = 2 on a quadrangle: definitions of the unit outward normal vector n κ k , element quadrature nodes y κ i (black bullets), surface quadrature node x κ k (gray bullets), and inner and outer traces u ± h at x κ k .

The geometrical quantities depend on the numerical method under consideration and examples will be given in section 4. By h we denote a consistent, conservative h ⟩ κ is supposed to be evaluated through a suitable quadrature rule, that includes N v volume quadrature points y κ i in κ together with the N f surface points x κ i on ∂κ introduced in (3.1) (see Fig. 1):

(3.3) ⟨u (n) h ⟩ κ = Nv i=1 ν κ i u h (y κ i , t (n) ) + N f i=1 β κ i u - h (x κ i , t (n) ),
where the weights ν κ 1≤i≤Nv ≥ 0 and β κ 1≤i≤N f > 0 are assumed to satisfy

(3.4) Nv i=1 ν κ i + N f i=1 β κ i = 1.
This manuscript is for review purposes only.

Lemma 3.1 and Corollary 3.3 give a theoretical basis on the existence of such a quadrature on polygonal or polyhedral mesh elements. An explicit example is also

given in [?] in the case of triangles and in section 4.

Lemma 3.1. Let κ be a compact subset of R d and P κ a finite dimensional subspace of C 0 (κ, R), that contains the constant function f ≡ 1 κ . Suppose that there exists a quadrature (ϖ κ i , y κ i ) 1≤i≤Nv with positives weights ϖ κ i > 0 and points y κ i in κ that integrates exactly products of functions in P κ :

(3.5) 1 |κ| κ f (y)g(y)dV = Nv i=1 ϖ κ i f (y κ i )g(y κ i ) ∀f, g ∈ P κ .
Then, for given points

(x κ i ) 1≤i≤N f in κ there exist nonnegative (ν κ i ) 1≤i≤Nv and positive (β κ i ) 1≤i≤N f coefficients such that ⟨f ⟩ κ = Nv i=1 ν κ i f (y κ i ) + N f i=1 β κ i f (x κ i ) ∀f ∈ P κ . Proof. It is known that (f, g) → 1 |κ| κ f (y)g(y)
dV defines a scalar product on C 0 (κ, R) and therefore on P κ . Then using the Riesz representation theorem, every linear form φ : P κ → R can be represented using this scalar product: there exists f φ ∈ P κ such that for every

g ∈ P κ , φ(g) = 1 |κ| κ f φ (y)g(y)dV , then let α φ i = ϖ κ i f φ (y κ i ),
we obtain for every

g ∈ P κ , φ(g) = Nv i=1 α φ i g(y κ i ). Now, since f → N f i=1 s κ k f (x κ i )
defines a linear form on P κ , for any s κ k > 0, it can be represented this way: there exist

(α κ i ) 1≤i≤Nv such that (3.6) N f i=1 s κ i f (x κ i ) = Nv i=1 α κ i f (y κ i ) ∀f ∈ P κ .
As the constant function is in P κ , we have for

f in P κ , ⟨f ⟩ κ = Nv i=1 ϖ κ i f (y κ i ), so for ε κ > 0 ⟨f ⟩ κ = Nv i=1 ϖ κ i f (y κ i ) = Nv i=1 ϖ κ i f (y κ i ) -ε κ N f i=1 s κ i f (x κ i ) + ε κ N f i=1 s κ i f (x κ i ) = Nv i=1 ϖ κ i f (y κ i ) -ε κ Nv i=1 α κ i f (y κ i ) + ε κ N f i=1 s κ i f (x κ i ) = Nv i=1 (ϖ κ i -ε κ α κ i )f (y κ i ) + ε κ N f i=1 s κ i f (x κ i ).
Since the ϖ κ i are positive, for ε κ = min {i:α κ i >0} (

ϖ κ i α κ i ) > 0, the (ϖ κ i -ε κ α κ i ) are nonnegative, then (3.3) holds with (3.7) ν κ i = ϖ κ i -ε κ α κ i ≥ 0 ∀1 ≤ i ≤ N v , β κ i = ε κ s κ i > 0 ∀1 ≤ i ≤ N f .
Remark 3.2. Note that κ is usually a polyhedron and P κ a polynomial space, we can subdivide κ into simplices and since there are quadrature rules integrating

This manuscript is for review purposes only.

exactly arbitrary order polynomials on simplices, the previous lemma can be applied.

Likewise, the existence of the quadrature in Lemma 3.1 is required only for modal methods or when the DOFs are not defined at the x κ k in (3.1), so the present framework also holds for non polynomial nodal approximations with DOFs at the faces as in [?].

The following corollary allows to explicitly define the quadrature (3.3) for modal polynomial based methods.

Corollary 3.3. Given a basis (ϕ j ) 1≤j≤Np of P κ , of dimension N p , which is orthonormal with respect to the inner product (i.e., κ ϕ i (x)ϕ j (x)dV = δ ij where δ ij is the Kronecker symbol), then the α κ i in (3.7) read

(3.8) α κ i = ϖ κ i Np j=1 N f k=1 s κ k ϕ j (x κ k )ϕ j (y κ i ).
Proof. We know that there exists g ∈ P κ such that (see proof of Lemma 3.1) (3.9)

N f i=1 s κ i f (x κ i ) = Nv i=1 ϖ κ i f (y κ i )g(y κ i ) ∀f ∈ P κ . Expanding g ≡ Np k=1
g k ϕ k in the orthonormal basis and using (3.9) with f ≡ ϕ j we get

N f i=1 s κ i ϕ j (x κ i ) = Nv i=1 ϖ κ i ϕ j (y κ i ) Np k=1 g k ϕ k (y κ i ) = g j ,
for all 1 ≤ j ≤ N p , by orthonormality of the basis. Substituting g in (3.9) by its expansion in the basis gives

N f i=1 s κ i f (x κ i ) = Nv i=1 ϖ κ i f (y κ i ) Np j=1 N f k=1 s κ k ϕ j (x κ k )ϕ j (y κ i ) ∀f ∈ P κ ,
and we conclude by comparing this result with (3.6).

Finally, the scheme is assumed to preserve uniform states in the following sense:

(3.10)

N f k=1 s κ k n κ k = 0,
which is a discrete version of

h λ (u L , u R , n) = f (u L ) • n + f (u R ) • n 2 - λ 2 (u R -u L ),
This manuscript is for review purposes only.

for λ ≥ |λ|(u L , u R , n) defined in (2.
2). Note that the Rusanov flux is derived from the following ARS:

W λ (ξ, u L , u R , n) =    u L , ξ < -λ, u L +u R 2 -1 2λ (f (u R ) • n -f (u L ) • n), -λ < ξ < λ, u R , λ < ξ, so from (2.
3) we know that it is IDP, see also [?]. It is also entropy stable in the sense of Definition 2.2 with Q(•, •, •) defined by

(3.12) Q λ (u L , u R , n) = q(u L ) • n + q(u R ) • n 2 - λ 2 η(u R ) -η(u L ) .
We now state a result that will allow us to rewrite (3.1) with updates of threepoint schemes.

Lemma 3.4 (pseudo-equilibrium state). Suppose that the numerical scheme satisfies (3.1), (3.3) and (3.10). Let B be an invariant domain for (1.1) and suppose that the internal traces u

- h (x κ k , t (n) ) 1≤k≤N f are in B, then there exist u ⋆ κ = u ⋆ κ (t (n) ) in B and λ ⋆ κ = λ ⋆ κ (t (n) ) > 0 finite such that (3.13) N f k=1 s κ k h λ ⋆ κ u ⋆ κ , u - h (x κ k , t (n) ), n κ k = 0, with h λ ⋆ κ defined in (3.11) and λ = λ ⋆ κ where (3.14) λ ⋆ κ ≥ max 1≤k≤N f |λ|(u ⋆ κ , u - h (x κ k , t (n) ), n κ k ) ,
and the pseudo-equilibrium state is defined by

(3.15) u ⋆ κ = N f k=1 γκ k u - h (x κ k , t (n) ) - f u - h (x κ k , t (n) ) • n κ k λ ⋆ κ , γκ k := s κ k S κ .
Proof. For the sake of clarity, we remove the time dependence of u h since all evaluations are done at t (n) and write u

- h (x κ k ) for u - h (x κ k , t (n)
). We first remark that from (3.2) and (3.15), we have

(3.16) N f k=1 γκ k = 1.
We introduce the following two sequences:

(3.17)

                 u ⋆ 0 = N f k=1 γκ k u - h (x κ k ), λ 0 = max 1≤k≤N f |λ|(u ⋆ 0 , u - h (x κ k ), n κ k ) , λ p+1 = max λ p , 1 d(u ⋆ p ,∂B) + max 1≤k≤N f |λ|(u ⋆ p , u - h (x κ k ), n κ k ) , p ≥ 0, u ⋆ p+1 = N f k=1 γκ k u - h (x κ k ) + u ⋆ p 2 - f u - h (x κ k ) • n κ k -f (u ⋆ p ) • n κ k 2λ p+1 , p ≥ 0,
where d(u, ∂B) = inf v∈∂B ∥u-v∥ is the distance from u to the boundary of B. We will show that both sequences converge and the limits satisfy (3.13) and (3.14). We first

This manuscript is for review purposes only.

remark that (λ p ) is non decreasing and therefore converges to some λ ⋆ κ in R ∪ {+∞} and that for all p ≥ 0,

u ⋆ p is in B since λ p+1 ≥ |λ|(u - h (x κ k ), u ⋆ p , n κ k ) so
u ⋆ p+1 = N f k=1 γκ k u - h (x κ k ) + u ⋆ p 2 - f (u - h (x κ k )) • n κ k 2λ p+1 = u ⋆ p 2 + 1 2 N f k=1 γκ k u - h (x κ k ) - 1 2λ p+1 N f k=1 γκ k f (u - h (x κ k )) • n κ k .
Then the first step of the sequence (3.17) for u ⋆ κ reads

u ⋆ 1 = N f k=1 γκ k u - h (x κ k ) - 1 2λ 1 N f k=1 γκ k f (u - h (x κ k )) • n κ k ,
so applying the recurrence relation p more times, we obtain

(3.18) u ⋆ p+1 = N f k=1 γκ k u - h (x κ k ) - p+1 i=1 1 2 i λ (p+2-i) × N f k=1 γκ k f (u - h (x κ k )) • n κ k .
Let show that p i=1

1 2 i λ (p+1-i) above converges to 1 λ ⋆ κ . Assume λ ⋆ κ is finite and let ϵ > 0, since λ p converges to λ ⋆ κ , there exists p 0 such that | 1 λp -1 λ ⋆ κ | < ϵ for all p > p 0 .
Then for p > p 0 , we set

p i=1 1 2 i λ (p+1-i) - 1 λ ⋆ κ = p k=1 1 2 (p+1-k) λ k - 1 λ ⋆ κ = p0 k=1 1 2 (p+1-k) λ k + p k=p0+1 1 2 (p+1-k) λ k - 1 λ ⋆ κ = 1 2 p+1 p0 k=1 2 k λ k + p k=p0+1 2 k 2 p+1 1 λ k - 1 λ ⋆ κ - 1 λ ⋆ κ 1- p k=p0+1 2 k 2 p+1
and using a triangle inequality we obtain This manuscript is for review purposes only.

p i=1 1 2 i λ (p+1-i) - 1 λ ⋆ κ ≤ 1 2 (p+1) p0 k=1 2 k λ k + p k=p0+1 2 k 2 p+1 1 λ k - 1 λ ⋆ κ + 1 λ ⋆ κ 1 - p k=p0+1 2 k 2 p+1 ≤ 1 2 (p+1) p0 k=1 2 k λ k + p-p0 i=1 1 2 i ϵ + 1 λ ⋆ κ 1 - p-p0 i=1 1 2 i = 1 2 p+1 p0 k=1 2 k λ k + 1 - 1 2 p-p0 ϵ + 1 2 p-p0
Let now prove that λ ⋆ κ is finite by contradiction. Suppose that λ p → +∞, then

u ⋆ p → N f k=1 γκ k u - h (x κ k ), but the application u → max 1≤k≤N f (|λ|(u, u - h (x κ k ), n κ k )) + 1 d(u, ∂B)
is continuous in the interior of B and is hence locally bounded around

N f k=1 γκ k u - h (x κ k ),
implying that λ p is bounded and λ ⋆ κ is finite which is a contradiction. By (3.17), we

obviously have λ p+1 ≥ max 1≤k≤N f |λ|(u ⋆ p , u - h (x κ k ), n κ k
) and passing the inequality to the limit we obtain (3.14).

Let now prove that u ⋆ κ is in B, since for all p ≥ 0, u ⋆ p is in B, we already know that u ⋆ κ is in the closure B. Now by contradiction, assuming that u ⋆ κ is not in B, we necessarily have u ⋆ κ in ∂B, so d(u ⋆ p , ∂B) → 0 inducing λ p → +∞ by (3.17) which a contradiction. It remains to prove (3.13). Using (3.10), we add

1 λ ⋆ κ f (u ⋆ κ )• N f k=1 γκ k n κ k =
0 to (3.15) and get

u ⋆ κ = N f k=1 γκ k u - h (x κ k ) - f (u - h (x κ k )) • n κ k + f (u ⋆ κ ) • n κ k λ ⋆ κ .
Moving u ⋆ κ to the right-hand side and using (3.16) we obtain

N f k=1 γκ k u - h (x κ k ) -u ⋆ κ - f (u - h (x κ k )) • n κ k + f (u ⋆ κ ) • n κ k λ ⋆ κ = 0,
and multiplying the above quantity by -

λ ⋆ κ 2 N f i=1 s κ i = - λ ⋆ κ
2 S κ we finally get

N f k=1 s κ k f (u - h (x κ k )) • n κ k + f (u ⋆ κ ) • n κ k 2 -λ ⋆ κ (u - h (x κ k ) -u ⋆ κ ) 2 =
± h (x κ k , t (n) ) 1≤k≤N f and u h (y κ i , t (n) ) 1≤i≤Nv
are in B an invariant domain for (1.1), then under the following condition on the time step

(3.19) ∆t (n) max κ∈Ω h max 1≤k≤N f s κ k β κ k max λ ⋆ κ , |λ| u - h (x κ k , t (n) ), u + h (x κ k , t (n) ), n κ k ≤ 1 2 ,
where λ ⋆ κ is defined by (3.14), ⟨u

(n+1) h ⟩ κ is also in B.
Proof. Once again we remove the time dependence of u h when possible for the sake of clarity, except when explicitly needed in the evaluation of ⟨u

⟨u (n+1) h ⟩ κ =⟨u (n) h ⟩ κ -∆t (n) N f k=1 s κ k h(u - h (x κ k ), u + h (x κ k ), n κ k ) -h λ ⋆ κ (u ⋆ κ , u - h (x κ k ), n κ k )
This manuscript is for review purposes only.

= Nv i=1 ν κ i u h (y κ i ) + N f k=1 β κ k u - h (x κ k ) -∆t (n) s κ k h(u - h (x κ k ), u + h (x κ k ), n κ k ) -h λ ⋆ κ (u ⋆ κ , u - h (x κ k ), n κ k ) = Nv i=1 ν κ i u h (y κ i , t (n) ) + N f k=1 β κ k U κ,n k , (3.20)
where, from Lemma 2.6 and the condition (3.19), the updates

(3.21) U κ,n k := u - h (x κ k ) - ∆t (n) s κ k β κ k h u - h (x κ k ), u + h (x κ k ), n κ k -h λ ⋆ κ u ⋆ κ , u - h (x κ k ), n κ k , are in B. Then by (3.4) ⟨u (n+1) h
⟩ κ is a convex combination of quantities in B, which concludes the proof.

Note that when explicitly using the quadrature (3.3) defined by Lemma 3.1

(see subsection 4.1 for an example), from (3.7) and (3.8) and the definition ε κ = min {i:α κ i >0} (

ϖ κ i α κ i
), we have

(3.22) s κ k β κ k = 1 ε κ = max {i:α κ i >0} ( α κ i ϖ κ i ) = max 1≤i≤Nv Np j=1 N f l=1 s κ l ϕ j (x κ l )ϕ j (y κ i ) ∀1 ≤ k ≤ N f ,
and the CFL condition (3.19) now reads

(3.23) ∆t (n) max κ∈Ω h 1 ε κ max λ ⋆ κ , max 1≤k≤N f |λ| u - h (x κ k , t (n) ), u + h (x κ k , t (n) ), n κ k ≤ 1 2 .
Remark 3.6. Since the set of states Ω a is in general a convex invariant domain, Theorem 3.5 can then be applied to B = Ω a to ensure robustness of the scheme. Likewise, using the convex combination (3.20) with (3.4), other properties such as discrete minimum and maximum principles can be proved according to the conservation laws under consideration.

We end this section with a result on entropy stability that directly follows from Theorem 3.5.

(n+1) h ⟩ κ ) ≤ ⟨η(u (n) h )⟩ κ -∆t (n) N f k=1 s κ k Q u - h (x κ k , t (n) ), u + h (x κ k , t (n) ), n κ k .
Proof. Again we remove the time dependence when possible. By (3.1), (3.10), and (3.15), we decompose u ⋆ κ into

u ⋆ κ 2 + u ⋆ κ 2 + 1 2λ ⋆ κ f (u ⋆ κ )• N f k=1 γκ k n κ k = N f k=1 γκ k u ⋆ κ + u - h (x κ k ) 2 - f (u - h (x κ k )) -f (u ⋆ κ ) • n κ k 2λ ⋆ κ ,
and then using (2.3) and (2.4)

with ∆t h = 1 2λ ⋆ κ we obtain η(u ⋆ κ ) ≤ N f k=1 γκ k η(u ⋆ κ ) + η(u - h (x κ k )) 2 - q(u - h (x κ k )) • n κ k -q(u ⋆ κ ) • n κ k 2λ ⋆ κ
This manuscript is for review purposes only.

= N f k=1 γκ k η(u ⋆ κ ) + η(u - h (x κ k )) 2 - q(u - h (x κ k )) • n κ k + q(u ⋆ κ ) • n κ k 2λ ⋆ κ
where we have used (3.10) to subtract

q(u ⋆ κ ) λ ⋆ κ • N f k=1 γκ k n κ k = 0 in the last step.
Moving everything to the left-hand side, we deduce that (3.12) satisfies (3.25)

N f k=1 s κ k Q λ ⋆ κ u ⋆ κ , u - h (x κ k ), n κ k = S κ λ ⋆ κ N f k=1 γκ k q(u - h (x κ k )) • n κ k + q(u ⋆ κ ) • n κ k 2 - λ ⋆ κ 2 η(u - h (x κ k )) -η(u ⋆ κ ) ≤ 0.
Finally, by (3.4) and (3.20) and then (3.21) and the entropy stability of h and h λ ⋆ κ in the sense of Definition 2.2, we obtain η(⟨u

(n+1) h ⟩ κ ) ≤ Nv i=1 ν κ i η u h (y κ i ) + N f k=1 β κ k η(U κ,n k ) ≤ Nv i=1 ν κ i η u h (y κ i ) + N f k=1 β κ k η u - h (x κ k ) - ∆t (n) s κ k β κ k Q u - h (x κ k ), u + h (x κ k ), n κ k -Q λ ⋆ κ u ⋆ κ , u - h (x κ k ),
               λ 0 = 1 ϑ max 1≤k≤N f |λ| u - h (x κ k , t (n) ), u + h (x κ k , t (n) ), n κ k , u ⋆ 0 = N f k=1 γκ k u - h (x κ k , t (n) ), λ p+1 = max λ p , 1 ϑ max 1≤k≤N f |λ|(u ⋆ p , u - h (x κ k , t (n) ), n κ k ) , p ≥ 0, u ⋆ p+1 = N f k=1 γκ k u - h (x κ k , t (n) ) - f u - h (x κ k ,t (n) ) •n κ k λp+1 , p ≥ 0,
where

ϑ = max 0 ≤ t ≤ 1 : u ⋆ 0 - t λ1 N f k=1 γκ k f u - h (x κ k , t (n) ) • n κ k ∈ B ,
This manuscript is for review purposes only.

λ1 = max 1≤k≤N f |λ| u - h (x κ k , t (n) ), u + h (x κ k , t (n) ), n κ k , |λ| u ⋆ 0 , u - h (x κ k , t (n) ), n κ k .
Since u ⋆ 0 ∈ B, we have 0 < ϑ ≤ 1 and then u ⋆ 1 ∈ B by definition of ϑ. Now λ p is non-decreasing and

u ⋆ p+1 ∈ [u ⋆ 0 , u ⋆ p ] ⊂ [u ⋆ 0 , u ⋆ 1 ]
for all p ≥ 1 which is enough to prove convergence of both sequences and ensure that u ⋆ p ∈ B. Since [u ⋆ 0 , u ⋆ 1 ] is closed, the limit is also in B (no need to add the distance term as in (3.17)). Obviously, the limits satisfy (3.13), (3.14) and (3.15). But this time, if λ p+1 = λ p for p ≥ 1, then u ⋆ p+1 = u ⋆ p and from this point both sequences are stationary. Now, the evaluation of u ⋆ p is really fast and we only need to evaluate max

1≤k≤N f (|λ|(u ⋆ p , u - h (x κ k , t (n) ), n κ k )).
In the supplementary materials, we investigate the possibility of using local wave speed estimates λ ⋆k κ at each vertex x k κ in ∂κ instead of the global estimate λ ⋆ κ . Numerical experiments however show that there is no improvement in the solution, while the algorithm for computing λ ⋆k κ and u ⋆ κ is slower to converge.

λ i = 1, we have ψ( i λ i u i ) ≥ min i ψ(u i ) for all u i in B.
From Theorem 3.5 we see that for any quasiconcave function ψ we have

(3.27) ψ(⟨u (n+1) h ⟩ κ ) ≥ m ψ κ := min ψ u h (y κ i , t (n) ) 1≤i≤Nv , ψ U κ,n k 1≤k≤N f
, where the updates U κ,n k are defined in (3.21). We now limit the solution around its cell-average ⟨u

⟩ κ so that it satisfies the same bounds and we rely on scaling limiters introduced in [?, ?] to enforce the bounds from quasiconcave functions to points where u h needs to be evaluated. The limited solution is thus defined as

(3.28) ũ(n+1) h (x) ≡ (1 -θ κ )u (n+1) h (x) + θ κ ⟨u (n+1) h ⟩ κ ∀κ ∈ Ω h , where (3.29) θ κ = min z∈(y κ 1≤i≤Nv )∪(x κ 1≤i≤N f ) max 0 ≤ t ≤ 1 : ψ (1-t)u h (z, t (n+1) )+t⟨u (n+1) h ⟩ κ ≥ m ψ κ .
This strategy may be applied to a finite family (ψ i ) 1≤i≤nc of n c quasiconcave functions by using the minimum value θ κ = min{θ κ (ψ i ) : 1 ≤ i ≤ n c }. The limiter (3.28) is then applied locally to each cell κ. The cell-average is not modified,

⟨ũ (n+1) h ⟩ κ = ⟨u (n+1) h ⟩ κ ,
u h (x, t) = Np k=1 ϕ κ k (x)U κ k (t) ∀x ∈ κ, κ ∈ Ω h , ∀t ≥ 0,
This manuscript is for review purposes only.

Algorithm 3.1 Algorithm flowchart of the IDP time integration with limiting. In 3, R h denotes the space residuals and depends on the chosen discretization. where the basis functions ϕ κ k span the function space V p h (κ) restricted onto κ and N p is its dimension. For simplicity we consider a discrete scheme of the form

(4.2) M κ k U κ k,n+1 -U κ k,n ∆t (n) + R κ k (u (n) h ) = 0 ∀κ ∈ Ω h , 1 ≤ k ≤ N p , n ≥ 0, where U κ k,n = U κ k (t (n)
), the M κ k are the entries of the diagonal mass matrix, and

R κ k (u h ) represent
V p h = κ∈Ω h V p h (κ) = {ϕ ∈ L 2 (D) : ϕ| κ • x κ ∈ P p ( K) ∀κ ∈ Ω h },
where P p ( K) is a polynomial space over a reference element K of polynomials of total degree lower or equal to p and N p = (p+d)! d!p! is its dimension. Each physical element κ is the image of K through the mapping x = x κ (ξ) with ξ = (ξ 1 , . . . , ξ d ).

Likewise, each face f in F h is the image of a reference face F through the mapping

x = x f (ξ 1 , . . . , ξ d-1
). We further define the Jacobians of the transformations by

J κ (x) = |x ′ κ (ξ)| and J f (x) = |x ′ f (ξ 1 , . . . , ξ d-1 )|.
We consider a quadrature rule (ξ V i , ω V i ) 1≤i≤Nv on K and denote y κ i = x κ (ξ i ) (see Fig. 1). Similarly, we consider a quadrature (ξ f k , ω f k ) 1≤k≤n f on F and denote the faces of κ (f j κ ) 1≤j≤Nκ where N κ is the number of faces of κ and the f j κ ∈ F are distinct.

We define the x κ i in (3.1) by x κ (j-1)n f +k = x f j κ (ξ

f j κ k ) for 1 ≤ j ≤ N κ and 1 ≤ k ≤ n f and N f = N κ × n f .
The DG method consists in defining a discrete weak formulation of (1.1) by multiplying it with test functions ϕ κ k spanning V p h and integrating over κ, using integration by parts and approximating f (u h ) • n by two-point numerical fluxes and the integrals by the quadrature rules. The space discretization in (4.

2) reads

This manuscript is for review purposes only.

V. CARLIER AND F. RENAC

(4.3) R κ k (u h ) = - Nv i=1 ω V i J κ (y κ i )f u h (y κ i , t (n) ) • ∇ϕ κ k (y κ i ) + N f i=1 ω f i J f (x κ i )h u - h (x κ i , t (n) ), u + h (x κ i , t (n) ), n κ i ϕ κ k (x κ i ).
We build an orthonormal basis with respect to the discrete inner product on κ by using a Gram-Schmidt algorithm, hence

M k κ = |κ| := Nv i=1 ω V i J κ (y κ i ).
Setting

ϕ κ k = 1 κ
s κ k = ω f k J f (x κ k ) |κ| and β κ k = ε κ s κ k in (3.
3) where ε κ is evaluated from (3.22). Then by (3.23), the scheme (4.2) and (4.3) is IDP under the condition

(4.4) ∆t (n) max κ∈Ω h 1 ε κ max λ ⋆ κ , max 1≤k≤N f |λ| u - h (x κ k , t (n) ), u + h (x κ k , t (n) ), n κ k ≤ 1 2 .

Discontinuous Galerkin

k = k(i 1 , . . . , i d ) := 1 + d j=1 i j (p + 1) j-1 0 ≤ i 1 , . . . , i d ≤ p.
We define a basis (ϕ κ k ) 1≤k≤Np of V p h (κ) by using tensor products:

ϕ κ k (x) = ϕ κ k (x κ (ξ)) = Π d j=1 ℓ ij (ξ j )
, where ℓ 0≤i≤p denote the ith Lagrange interpolation polynomial associated to ζ i the ith Gauss-Lobatto quadrature node with 2) and by ω i we denote the associated weight. We therefore have the following cardinality relation at quadrature points

ζ 0 = -1 < ζ 1 < • • • < ζ p = 1 (see Fig.
ξ k ′ = (ξ i ′ 1 , . . . , ξ i ′ d ) in K: ϕ κ k (x κ k ′ ) = ϕ κ k (x κ (ξ k ′ )) = δ i1,i ′ 1 . . . δ i d ,i ′ d with δ i,i ′ the
Kronecker symbol, so the DOFs correspond to the point values of the solution:

U κ k (t) = u h (y κ k , t) with y κ k = x κ (ξ k )
and interpolation and quadrature points are collocated, hence

N v = N p .
Let introduce the discrete derivative matrix with entries

D ij = ℓ ′ j (ζ i ) with 0 ≤ i, j ≤ p. The DGSEM discretization takes the form (4.2) with M κ k = ω V k J κ (y κ k ), ω V k = Π d j=1 ω ij for k = k(i 1 , . . . , i d ) and (4.5) R κ k (u h ) = 2ω V k d j=1 p l=0 D ij l h sym U κ k , U κ k ′ j , {J κ ∇ξ j } (k,k ′ j ) + N f i=1 ϕ κ k (x κ i )ω f i J f (x κ i ) h U κ k,n , u + h (x κ i , t), n κ i -f U κ k,n • n κ i where for d = 3 k ′ 1 = k(l, i 2 , i 3 ), k ′ 2 = k(i 1 , l, i 3 ) and k ′ 3 = k(i 1 , i 2 , l), hence k ′ j = k + (i ′ j -i j )(p + 1) j-1 , N f = 2d(p + 1) d-1 , ω f i = Π d-1
j=1 ω ij , and the

{J κ ∇ξ} (k,k ′ j ) = 1 2 J κ (y κ k )∇ξ j (ξ k ) + J κ (y κ k ′ j )∇ξ j (ξ k ′ j ) ,
This manuscript is for review purposes only.

κ = κ - κ + f u - h u + h n κ k x κ k y κ i • • • • • • • • • • • • • • • •
Fig. 2: Notations for the DGSEM and d = 2: inner and outer elements κ -and κ + ; traces u ± h on the interface f and of the unit outward normal vector n κ k ; element quadrature nodes y κ i and surface quadrature node x κ k that are also included in the y κ i .

have been introduced to keep conservation of the scheme [?]. By h sym we denote a twopoint flux supposed to be symmetric in the sense that

h sym (u, v, n) = h sym (v, u, n). Note that ϕ κ k (x κ i ) = 1 if x κ i = y κ k and ϕ κ k (x κ i ) = 0 else.
The choice of the y κ k in the quadrature (3.3) is not unique and we here use the following decomposition

⟨u (n) h ⟩ κ = 1 |κ| y κ i ∈κ ω V i J κ (y κ i )U κ i,n = 1 |κ| y κ i ∈int(κ) ω V i J κ (y κ i )U κ i,n + 1 |κ| N f i=1 1 α κ i ω V i J κ (x κ i )u h (x κ i , t (n) ),
where int(κ) denotes the interior of κ, while

α κ i = d if x κ i = y κ j is a vertex of the d-dimensional hexahedron, α κ i = d -1 if x κ i = y κ j is on some edge, α κ i = 1 else. Summing (4.2) over 1 ≤ k ≤ N v gives for the cell-averaged solution ⟨u (n+1) h ⟩ κ = ⟨u (n) h ⟩ κ - ∆t (n) |κ| Nv k=1 R κ k (u (n) h ) = ⟨u (n) h ⟩ κ - ∆t (n) |κ| N f i=1 ω f i J f (x κ i )h u - h (x κ i , t (n) ), u + h (x κ i , t (n) ), n κ i ,
by conservation of the DGSEM [?, ?] and providing that the so-called metric identities are satisfied at the discrete level [?]. This relation can be identified with (3.1) with

s κ k = ω f k J f (x κ k ) |κ|
. Then we can apply Theorem 3.5 with

β κ k = ω V k Jκ(x κ k ) α κ k |κ|
and the DGSEM scheme (4.2) is IDP under the condition

(4.6) ∆t (n) max κ∈Ω h max 1≤k≤N f α κ k ω f k J f (x κ k ) ω V k J κ (x κ k ) max λ ⋆ κ , |λ| u - h (x κ k , t (n) ), u + h (x κ k , t (n) ) ≤ 1 2 .
Remark 4.1. The CFL conditions (4.4) and (4.6) may be compared with results in the literature. First, bounding the wave speeds by λ ∞ and assuming a first-order scheme, p = 0, and straight-sided elements we have ε κ = |κ|/|∂κ| from (3.22), so (4.4)

reduces to ∆t (n) max κ∈Ω h |∂κ| |κ| λ ∞ ≤ 1 2 which is the estimate obtained in [?, Th. 3].
Then, assuming a Cartesian grid with elements of size h and p ≥ 1, using This manuscript is for review purposes only. We aim to test the robustness and efficiency of the present approach on simulations of (1.1) and (5.1) with discontinuous solutions on one-dimensional and unstructured two-dimensional grids. We use the modal DG method in subsection 4.1 and the DGSEM in subsection 4.2, where the time step is respectively evaluated from conditions (4.4) and (4.6) with the wave estimate λ ⋆ κ from algorithm (3.26). The cell-averaged solution is thus guaranteed to be IDP and we apply the limiter (3.28) to further impose the high-order solution to be IDP. We will compare the following as expected because the solution remains all the time positive, so does the IDP scheme with the smoothness indicator. Removing the smoothness indicator (IDP (no TCI)) affects the accuracy and order of convergence, especially on the finest mesh. We also tested the bound relaxation (IDP (Rel)), see [?, Sec. 7.6], which allows to recover the expected convergence order. These results are in agreement with results observed in the literature [?, ?, ?, ?, ?, ?, ?]. We conclude that the limiter may alter the formal accuracy of the scheme, but there exist ways to fix this issue.

α κ k = d, J f (x κ k )/J κ (x κ k ) = 2 h and min ω ij = 2 p(p+1) , (4.6) becomes ∆t (n) λ∞ h ≤ 1 2dp (
∥ L 1 (Ω) O 1 ∥e h ∥ L 2 (Ω) O 2 ∥e h ∥ L ∞ (Ω) O ∞ POS 1/8

Riemann problems.

We here consider Riemann problems (2.1) with initial data given in Tab. 3. Results are shown in Figs. 3 to 5. The flow features are well captured by the schemes with the four present strategies. We observe some spurious This manuscript is for review purposes only.

oscillations in the DGSEM solution with POS limiter that are slightly damped with the IDP limiter. The POS limiter only ensures that the solution remains in the set of states Ω a and does not control spurious oscillations, while the IDP strategy further reduces the oscillations. In contrast, the oscillations almost vanish with the modal DG scheme and there is no evident benefice of using the IDP limiter compared to the POS limiter. Imposing the IDP property at all quadrature points may result in stronger limitation of the solution with the modal DG scheme as these are more quadrature points (compare Figs. 1 and2). Note that all the computations (with either DGSEM, or modal DG scheme) require to apply one of the limiting strategy to avoid non-physical solution.

5.3. Two-dimensional problems. We now consider two-dimensional problems and first consider the reflection of a Mach 10 shock over a 30 • wedge [?]. Ahead of the shock, the gas is at rest and has a density of 1.4 and pressure of 1. Inflow and outflow conditions are applied at the left and bottom boundaries, while a symmetry condition is applied at the top boundary. Initially, the shock is located at x = 0 corresponding to the beginning of the wedge. We use an unstructured mesh with 164850 quadrangles to solve the horizontally moving shock interacting with the inclined wall where slip conditions are applied. In Fig. 6 we plot the density contours and post-process the solution by splitting each element into four subelements to evaluate the resolution within elements, but also oscillations of the polynomial solution. We observe qualitatively similar results for all the tests. Some spurious oscillations remain after limiting, but the computations proved to be robust.

We then consider the diffraction of a shock wave at a corner. The computational domain is the union of [0, 1] × [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF][START_REF] Crean | Entropy-stable summation-by-parts discretization of the Euler equations on general curved elements[END_REF] and [START_REF] Avriel | r-convex functions[END_REF][START_REF] Don | Numerical study of pseudospectral methods in shock wave applications[END_REF] 6. Conclusions. We here investigate robustness and stability properties of highorder spectral discontinuous methods on general unstructured meshes with explicit time stepping for the approximation of hyperbolic systems of conservation laws. We derive a condition on the time step to guarantee that the cell-averaged approximate solution is a convex combination of DOFs at preceding time step and updates of IDP and entropy stable three-point schemes. As a consequence, the cell-averaged scheme preserves the invariant domains of the equations. For the sake of illustration, we apply

This manuscript is for review purposes only. balance over each mesh element, and the existence of a quadrature rule including the traces to evaluate the cell-averaged solution. We here prove their existence in the general case and provide an iterative algorithm to evaluate the pseudo-equilibrium state. We illustrate these results with the classical modal DG and DGSEM schemes.

Numerical experiments in one and two space dimensions are provided to illustrate the robustness and stability of the present approach. The extension of this framework to parabolic systems of conservation laws, e.g., the compressible Navier-Stokes, is a

This manuscript is for review purposes only.

1 .

 1 Introduction. Let D ⊂ R d be an open domain with d the space dimension.

  s κ k > 0 are local contributions to |f | |κ| with |f | and |κ| approximations of the face surface and element volume, and we introduce (3.2)

( 2 . 5 )

 25 , and IDP (see Definition 2.1) two-point flux derived from an ARS. The cellaveraged solution ⟨u (n)

  which clearly converges to ϵ when p → ∞, proving our statement and (3.15).

  0, which is exactly (3.13) from the definition of the Rusanov flux (3.11).3.3. Invariant domain preserving schemes. Using Lemma 3.4 we now state and prove the main result of this work in the theorem below. Theorem 3.5 (Time step condition). Assume that the numerical scheme satisfies (3.1), (3.3) and (3.10) and assume that u

  add the trivial quantity ∆t (n) × (3.13) to (3.1) and use (3.3) to get

Corollary 3 . 7 .

 37 Under the assumptions of Theorem 3.5 and further assuming that the two-point flux in (3.1) is entropy stable in the sense of Definition 2.2, we have (3.24) η(⟨u

3. 5 .

 5 Application to convex limiting. Convex limiting enforces the numerical solution to preserve invariant domains [?] through quasiconcave constraints [?]. Let recall that a function ψ : B → R is quasiconcave iff. for every family of convex coefficients (λ i ) ≥ 0, with

4 .

 4 and cellwise discrete conservation (3.1) still holds. Algorithm 3.1 describes the procedure of the IDP time integration with the above limiter at each time step with an explicit time integration. Examples of high-order spectral discontinuous methods. We here review some high-order spectral discontinuous approximations of (1.1) which satisfy the assumptions of discrete conservation (3.1), existence of a quadrature rule (3.3) and preservation of uniform states (3.10). Then, there exists pseudo-equilibrium states u ⋆ κ such that Lemma 3.4 and Theorem 3.5 hold, and the limiter (3.28) can be applied.In the following, we consider a partition Ω h of D ⊂ R d , composed of non-overlapping and non-empty elements κ, and by F h we denote the set of faces in the partition. The approximate solution is sought under the form(4.1)

  the indicator function of κ, the first sum vanishes and we obtain (3.1) with

  Spectral Element Method. In the DG spectral element method (DGSEM), the reference element is an hypercube : K = I d := {ξ = (ξ 1 , . . . , ξ d ) : -1 ≤ ξ j ≤ 1} and the polynomial space P p (I d ) is formed by tensor products of polynomials of degree at most p in each direction. The approximate solution is sought under the form (4.1) where (U κ k ) 1≤k≤Np are the N p = (p + 1) d DOFs in the element κ with indexing

2 . 4 . 3 .

 243 p+1) , which results in a time step reduced by a factor 2 when compared with [?,Th. 3.1]. This latter result however assumes Cartesian meshes to rewrite the numerical flux balance of opposite faces as the space residuals of a three-point scheme, thus gaining this factor Other methods. Properties (3.1),(3.3), and (3.10) also hold for other discretely conservative spectral difference methods on general curved elements provided the discretization operators satisfy the metric identities at the discrete level, which imposes some limits on the order of approximation of the mesh elements compared to the approximation order of the solution [?]. The limiter (3.28) may hence be applied to make these methods invariant domain preserving. We list some examples below. The skew-symmetric entropy stable modal DG methods [?] uses skew-hybridized summation-by-parts (SBP) operators allowing conservation and free-stream preservation under standard accuracy requirements of volume and surface quadratures. In [?] multidimensional discretization schemes based on SBP operators on general curved elements generalize the staggered finite differences from [?]. The discretizations with curved elements remain accurate, conservative, and entropy stable. Spectral differences [?] and staggered Chebyshev [?] methods on curved elements belong to a same family of conservative approximations satisfying the discrete metric identities. The methods use two sets of interpolation points for the solution and fluxes and impose the discrete residuals to be satisfied at solution points. Flux points contain points at faces of the elements where two-point numerical fluxes are used. Then, the space derivatives at solution points are evaluated by differencing the polynomials interpolating the fluxes. 5. Numerical experiments. Let consider the compressible Euler equations of gas dynamics. The conservative variables and fluxes in (1.ρvv ⊤ + pI d (ρE + p)v ⊤   , where ρ, v, and E denote the density, velocity vector, and specific total energy, respectively. The system is closed by defining the equation of state p = p( 1 ρ , e) with e = E -1 2 v • v the specific internal energy and the system is hyperbolic over the set of states Ω a = {u ∈ R d+2 : ρ > 0, v ∈ R d , e > 0}. We focus here on the polytropic ideal gas law p = (γ -1)ρe where γ = Cp Cv = 7 5 is the ratio of specific heats. The compressible Euler equations (1.1) and (5.1) possess the natural entropy -entropy flux pair η = -ρs, q = -ρsv, s = C v ln p ρ γ , and B = {u ∈ Ω a : s(u) ≥ s 0 }, with s 0 in R, is an invariant domain. We use the convex limiting strategy with the quasiconcave functions ψ 1 ≡ ρ and ψ 2 ≡ ρe.

5 . 1 .

 51 limiting strategies: a positivity limiter (POS) which imposes the solution to remain in B = Ω a thus extending [?] to unstructured grids; an IDP limiter which imposes the solution to remain in the convex hull of the states in (3.20) on density and internal energy. Imposing the IDP property may result in over-limiting of the solution and some strategies are usually applied such as bound relaxation [?], or subcell smoothness indicator [?]. We here follow the second strategy which relies on the smoothness indicator from [?] (see [?, Sec. 4.4] for details). Finally, we use the Suliciu pressure relaxation based numerical flux from [?, Sec. 2.4.6] at interfaces, while for h sym in the DGSEM scheme (4.5) we use the Kennedy and Grubber splitting from [?]. We systematically evaluated the convergence efficiency of algorithm (3.26) in terms of iterations to converge. Tab. 1 indicates the global average number of iterations evaluated over the whole computations for all schemes and different approximation orders. We observe that algorithm (3.26) converges very fast with a global average always lower than 1.2 iterations, thus meaning that most of the time u ⋆ κ = u ⋆ 1 given in (3.26) satisfies the requirement (3.14) so no more step is needed. Convection of an isentropic vortex. We first evaluate the accuracy of the fourth-order DGSEM scheme together with a limiter on the convection of an isentropic vortex in a uniform subsonic flow. Results with the modal DG method and with different approximation orders lead to the same conclusions. Error norms are gathered in Tab. 2 for different limiters. The POS scheme keeps high-order accuracyThis manuscript is for review purposes only.

  × [0, 11]. The initial condition is a pure right-moving shock of Mach 5.09, initially located at x = 0.5 and 6 ≤ y ≤ 11, moving into undisturbed air ahead of the shock with a density of 1.4 and pressure of 1. The boundary conditions are inflow at x = 0, 6 ≤ y ≤ 11, outflow at x = 13, 0 ≤ y ≤ 11, and 1 ≤ x ≤ 13, y = 0, and reflective at the walls 0 ≤ x ≤ 1, y = 6 and x = 1, 0 ≤ y ≤ 6, and top boundary 0 ≤ x ≤ 13, y = 11. Figure 7 shows the density contours (post-processed into four subelements) obtained on an unstructured mesh with 110242 quadrangles and we observe similar results as in the previous test. We finally consider the problem of a two-dimensional Kelvin-Helmholtz instability [?] in a square domain [-1, 1] 2 with periodic conditions and initial data as ρ 0 (x) = 1 2 + 3 4 B(x), u 0 (x) = 1 2 (B(x) -1), v 0 (x) = 1 10 sin(2πx) and p 0 (x) = 1 where B(x) = tanh(15y + 7.5) -tanh(15y -7.5). Though these are no shocks, the strong variations in the initial density and velocity fields make this test challenging the robustness and resolution capabilities of the methods. Figure 8 displays snapshots of density fields obtained on a Cartesian mesh where we observe a wide range of scales that slightly differ between different methods due to the sensitivity of the flow to small perturbations. All the computations where stable and robust with a resolution similar to computations with the same accuracy [?].

Fig. 3 :

 3 Fig. 3: Sod problem at t = 0.2 with p = 3 and N = 100 elements.

Fig. 4 :

 4 Fig. 4: Lax problem at t = 0.13 with p = 3 and N = 100 elements.

Fig. 5 :

 5 Fig. 5: Toro 4 problem at t = 0.035 with p = 3 and N = 100 elements.

Fig. 6 :

 6 Fig. 6: Double Mach reflection problem at t = 0.2: 41 equispaced density contours between 2.1 and 22.1.

Fig. 7 :

 7 Fig. 7: Shock diffraction problem at t = 2.3: 28 equispaced density contours between 0.5 and 7.5.

Fig. 8 :

 8 Fig. 8: Two-dimensional Kelvin-Helmholtz instability: density field at t = 10 with p = 3 and N = 64 × 64 elements.

  1.1) is sought under the form of discontinuous piecewise truncated series of analytic functions over a partition of the domain D. Such methods have been applied to a wide range of applications [?, ?], and have the potential to achieve high-order accuracy efficiently on modern parallel architectures [?, ?]. Unfortunately these approximations suffer from spurious oscillations around discontinuities of the exact solution due to

Gibbs phenomenon [?, ?] that may cause the approximate solution to become locally nonphysical, leading to robustness issues. A large body of research has been proposed to address such issues with, e.g., solution and flux limiters [?, ?, ?, ?, ?], entropy conservative subcell flux differencing [?, ?, ?], artificial viscosity [?, ?], shock-capturing terms [?, ?], entropy projection [?].

  W a (ξ, u L , u R , n)dξ is also in B. We have the following results linking IDP ARS and two-point numerical flux.Lemma 2.5 (Interface invariant domain preservation [?]). The ARS W a is IDP for B iff. for all u L , u R in B, and ∆t h

	1	λ
	2λ	-λ

  Let h W a and h W b be two-point fluxes from two different ARS that are IDP for B. Then, we have

	,
	and since 0 < ∆t h λ ≤ 1 this is a convex combination of states in B and therefore is in B. Conversely, taking ∆t h λ = 1 the above integrals are B iff. (2.9a) and (2.9b) hold
	and we conclude from Definition 2.4.
	This allows us to state the following result.
	Lemma 2.6.

Table 1 :

 1 Global average number of iterations required for algorithm (3.26) to converge for the different numerical experiments in section 5 (DMR: double Mach reflection; SDP: shock diffraction problem; KHI: Kelvin-Helmholtz instability) and different polynomial degrees 1 ≤ p ≤ 3.

	p	1	2	3	1	2	3	1	2	3
	test case		Sod			Lax			Toro 4
	POS DGSEM 1.02 1.03 1.06 1.00 1.05 1.13 1.02 1.08 1.09
	IDP DGSEM 1.02 1.03 1.06 1.00 1.05 1.12 1.02 1.10 1.10
	POS DGM	1.08 1.15 1.11 1.18 1.20 1.18 1.08 1.14 1.14
	IDP DGM	1.08 1.15 1.11 1.18 1.20 1.18 1.08 1.14 1.14
	test case		DMR			SDP			KHI
	POS DGSEM 1.02 1.14 1.18 1.12 1.18 1.18 1.04 1.11 1.15
	IDP DGSEM 1.02 1.14 1.18 1.12 1.18 1.18 1.04 1.11 1.15
	POS DGM	1.04 1.18 1.17 1.14 1.17 1.18 1.06 1.15 1.19
	IDP DGM	1.04 1.18 1.17 1.14 1.17 1.18 1.06 1.15 1.16

Table 2 :

 2 Convection of an isentropic vortex in a subsonic flow at Mach number M ∞ = 0.5. Domain [0, 1] 2 discretized with a Cartesian mesh with a 4th-order DGSEM (p = 3), 1/h×1/h elements and periodic boundary conditions. The vortex is isentropic with an exponential decay type profile and a radius of 0.05. Different norms of the error on the pressure field after 5 convection time scales, i.e., t = 5 × 1/U ∞ .

	limiter	h	∥e h

Table 3 :

 3 Initial conditions of Riemann problems (2.1) where x 0 indicates the abscissa separating the initial states.problem left state (ρ L , u L , p L ) ⊤ right state (ρ R , u R , p R ) ⊤ x 0

			1.3932e+01	-	3.6193e+01	-	5.6936e+02	-
		1/16 4.3118e+00 1.69 9.7078e+00 1.90 1.6683e+02	1.77
		1/32 1.1530e-01 5.22 3.9447e-01 4.62 9.2519e+00	4.17
		1/64 5.1484e-03 4.49 2.3328e-02 4.08 5.8875e-01	3.97
		1/8	1.3932e+01	-	3.6193e+01	-	5.6936e+02	-
	IDP	1/16 4.3118e+00 1.69 9.7078e+00 1.90 1.6683e+02 1/32 1.1530e-01 5.22 3.9447e-01 4.62 9.2519e+00	1.77 4.17
		1/64 5.1484e-03 4.49 2.3328e-02 4.08 5.8875e-01	3.97
		1/8	1.3559e+01	-	3.6117e+01	-	5.7004e+02	-
	IDP (no	1/16 4.3035e+00 1.66 1.0000e+00 1.85 1.7109e+02	1.74
	TCI)	1/32 1.1568e-01 5.22 3.9849e-01 4.65 8.6409e+00	4.31
		1/64 2.8711e-02 2.01 2.2270e-01 0.83 1.0147e+01 -0.23
		1/8	1.3932e+01	-	3.6193e+01	-	5.6936e+02	-
	IDP	1/16 4.3118e+00 1.69 9.7078e+00 1.90 1.6683e+02	1.77
	(Rel)	1/32 1.1530e-01 5.22 3.9447e-01 4.62 9.2519e+00	4.17
		1/64 5.1484e-03 4.49 2.3328e-01 4.08 5.8875e-01	3.97
	Sod		(1, 0, 1) ⊤		(0.125, 0, 0.1) ⊤	0
	Lax		(0.445, 0.698, 3.528) ⊤	(0.5, 0, 0.571) ⊤	0
	Toro 4 [?] (5.99924, 19.5975, 460.894) ⊤ (5.99242, -6.19633, 46.0950) ⊤ -0.1

This manuscript is for review purposes only.


domain: scco info

text:
Marcel Francis Hinss
email: marcel.hinss@isae-supaero.fr

Emilie S Jahanpour

Anke M Brock

Raphaëlle N Roy

Marcel F Hinss¹²

Emilie S Jahanpour¹

Anke M Brock²

Raphaëlle N Roy¹

Labeling mental fatigue for passive BCI applications: Accuracy vs applicability tradeoff

Introduction: Within BCI research, particularly passive-BCI (pBCI) research, mental fatigue estimation is of great popularity, as mental fatigue can be attributed to catastrophic events in risky work settings (e.g. aviation and nuclear plant domains). Mental fatigue is often understood as a psychobiological mental state, caused by prolonged cognitive activity resulting in an increased probability of a performance decrement [START_REF] Hinss | Cognitive effects of prolonged continuous human-machine interaction: The case for mental statebased adaptive interfaces[END_REF]. The most popular method to label brain activity data is time-on-task (TOT). In practice, data are often labeled using the first and last blocks of an experimental procedure for binary classification between non-fatigued and fatigued states [START_REF] Trejo | EEG-Based Estimation and Classification of Mental Fatigue[END_REF]. While these approaches often yield high classification accuracies, their usefulness in practical applications is debatable, as the classification is not directly linked to the behavior of the user. A more applicable approach may be to move towards performance estimation in the context of mentally fatiguing environments. Another aspect to consider is that using the entire frequency spectrum of the EEG signal may inflate classification accuracies, due for instance to different motor activity between conditions which generates different amounts of artifacts within epochs and trials. This study has been designed to directly assess the impact of the type of mental fatigue labeling on estimation accuracy using a surveillance task dataset. In addition, the impact of the range of frequencies used for estimation is also considered. We hypothesize that TOT labeling as well as using the entire EEG frequency spectrum will result in a more accurate classification as compared to using performance estimation and a restricted EEG frequency spectrum. Methods: Data from 13 participants performing a 60-minute visual search task was recorded using a 64-electrode EEG setup. The data was preprocessed, either filtering out any motor-related activity (>20 Hz) or keeping the entire spectrum (only a 50 Hz notch filter). The data was labeled using one of two approaches. The TOT approach used the first and third (out of 4) blocks. The second approach used behavioral performance to label the epochs. Data were then classified using a Riemannian minimum-distance-to-mean classifier. Accuracies were compared using a 2x2 ANOVA.

Results: Classification accuracies obtained using TOT labeling were significantly higher as compared to performance labeling (F(2,48)=89.87 p<0.001). While the difference between the two spectra used in this approach was not significant, using the entire frequency spectrum still resulted in an almost 6% increase (79.9% vs 73.9%) in accuracy. Discussion: As shown by these results, performance estimation is significantly more challenging as compared to TOT estimation. However, for real-world applications, TOT labeling has no real added value compared to actually knowing how much time has been spent on a task, whereas performance labeling presents more use cases for instance in order to predict and mitigate degraded attentional states via adaptive interfaces based on a pBCI [START_REF] Hinss | Cognitive effects of prolonged continuous human-machine interaction: The case for mental statebased adaptive interfaces[END_REF].

Figure 1 (

 1 Figure 1 (from left to right) Screenshot of the visual search task in which participants had to detect humans to launch a rescue. | Comparison of the two labeling approaches | Results of the analysis, showing the upper limit of the confidence interval for chance level (57.5%) according to [3].


Nesrine Bouhamed
email: bouhamednesrine34@gmail.com

Slim Souissi

Pierre Marechal

Mohamed Ben Amar

Olivier Lenoir

Physicochemical, Morphological and Thermal Characterization of Composites Based on Olive Wood Flour

Keywords: Physicochemical, Biocomposite, Coupling agent

This study deals with the evaluation of the physicochemical, morphological and thermal properties of a biocomposite material based on polypropylene (PP) reinforced with vegetable fibers olive wood flour (OWF). Initially, the olive wood flour was chemically treated with a coupling agent (amino-silane) and then characterized to provide the modification influx. In a second step, various PP/OWF bio-composites with variable contents of OWF (0, 10, 20 and 30% by weight) were implemented by twin-screw extrusion and then shaped by injection molding. The effect of the OWF content on the one hand and the chemical modification of the latter by adding silane as coupling agent on the other hand, are the two parameters that are the subject of the present sensitivity study. Finally, the characterization results are compared and discussed. The significant influence of the treatment of OWF with the silane coupling agent is demonstrated, clearly demonstrating that the silane promotes the adhesion of the plant-based filler that is OWF with the PP polymer matrix.

Introduction

The development of a composite material makes it possible to obtain properties that are better than those of the constituents that compose it. Thus, it is imperative to know the composition and the properties of the constituents, as well as the stages of implementation to understand the interaction mechanisms between the constituents. This study aims to focus on this approach and allowed us to determine the physicochemical, thermal and morphological characteristics of the fibers, in this case olive wood flour (OWF), the matrix, polypropylene (PP), and elaborate composite (PP/OWF).

2 Materials and Method

Materials

The polypropylene (PP) used in this study is a standard homopolymer PP H9069 (Braskem, Brazil). FBO is recovered from wood waste generated during agricultural work. It is sieved and stored to protect it from moisture. Two types of FBO were used as reinforcing fillers, untreated (OWF) and treated (OWF) at 3% by weight with a coupling agent [START_REF] Arbelaiz | Mechanical properties of flax fibre/polypropylene composites[END_REF], amino-propyl-triethoxy-silane (OWFT).

Filler Characterization

The biomass used in this study consisted of olive wood flour and sieved to obtain finer particles, from 84.6 to 254.9 lm in size. The chemical composition of the dried fibers was determined according to [START_REF] Jayaramudu | Characterization of new natural cellulosic fabric Grewia tilifolia[END_REF]. It was determined by carrying out a series of extractions to isolate the various substances contained in the fibers such as: cellulose, hemicellulose and lignin. The length and aspect ratio of OWF filler was determined by Morfi analyses using JEOL JSM-5400 Morfi equipment

Elaboration of Bio-composite

In order to eliminate all the moisture absorbed and to avoid the formation of agglomerates, the FBOs were pre-dried at 105 °C for 24 h before extrusion. The PP polymer and the FBO fibers were mixed in a co-rotating twin-screw extruder (DSM Xplore Netherlands) with 3 heating zones. The 10 mm diameter screw was set at a rotational speed of 100 rpm. The mixture is extruded through a cylindrical matrix of diameter 1 mm at a rate of 200 g/h. The rod is cooled, then cut into granules which are in turn dried for 24 h at 60 °C. After extrusion, the dried granules are injection molded using a molding machine (Demag Plastic Group Extra 50-200). The injection is produced at a temperature of 180 °C and a pressure of 1000 bar to produce the biocomposite samples required for the tests, in the form of plates and test pieces. Different charge rates of 0, 20 or 30% have been developed.

Scanning Electron Microscopy (SEM)

Scanning electron microscopy (SEM) was performed with a JEOL JSM-5400 SEM microscope operating at 15 kV. The specimens were carried out under liquid nitrogen to avoid plastic deformation of the matrix and to obtain a well-defined fiber-matrix interface. Prior to analysis, all samples were coated with a 10 nm layer of gold to avoid sample charging under the electron beam.

Thermo-Gravimetric Analyzes (ATG)

Thermogravimetric analyzes (ATG) were performed using Perkin Elmer Pyris-1 equipment (USA) under a nitrogen flow up to 800 °C and with a temperature rise rate of 10 °C/min. The objective is to study the thermal stability of our samples (fibers alone, virgin polymer, and composite materials).

Results and Discussion

Morphology and Chemical Composition of Olive Wood Flour

Olive wood flour is similar to a composite material reinforced by longitudinally bonded fibrils of cellulose, lignin and other non-cellulosic materials. Waxes and impurities regularly cover the whole [START_REF] Fiore | Characterization of a new natural fiber from Arundo donax L. as potential reinforcement of polymer composites[END_REF][START_REF] Mwaikambo | Chemical modification of hemp, sisal, jute, and kapok fibers by alkalization[END_REF]. Figure 1 shows the micrographs obtained by the SEM show morphological transformations of the fibers following chemical treatment with silane The surface of raw olive wood flour is smooth, while that of the treated FBO shows exposed fibrils. The reduction in the level of non-cellulosic material such as lignin, hemicellulose and wax, thanks to the treatment performed, has a positive effect on the roughness of the fiber surface, which increases visibly according to the micrograph [START_REF] Narendar | Chemical treatments of coir pith: morphology, chemical composition, thermal and water retention behavior[END_REF].

The chemical composition of the OWF. It confirmed the lignocellulosic character of the filler being composed of lignin, hemicelluloses and cellulose, with content about 21.9%, 29.2% and 45.9%, respectively.

Morphological and Thermal Characterization of Elaborate Biocomposites

By switching from crude fiber-based biocomposite (Fig. 2 (a)) to processed fiber-based fiber (Fig. 2 (b)), the number of these cavities decreases. As a result, the improvement in adhesion between the fibers and the matrix can be distinguished by addition of the coupling agent (the treatment with the silane). From the PP matrix. As shown in photo (a) of Fig. 2, the fiber surfaces are smooth and clean, showing poor inter-facial adhesion between the fibers and the matrix. On the other hand, the photo (b) of Fig. 2 shows that for the bio-composite loaded by the treated FBO the fibers are closely linked to the matrix with a surface characterized by a coating of a layer of the polymer matrix (PP). Indeed, it can be seen that the fracture occurred at the cross-section of the fibers themselves. This confirms that there is a transfer of the stresses from the matrix to the load, because of the adhesion to the contact between fibers and matrix. The fibers therefore participate in this case in the energy required for breaking and thus ensure good cohesion of the entire composite material.

These results are in agreement with previous studies which confirmed that the chemical treatment with silane of natural fibers used as reinforcements improves the cohesion between the fibers and the matrix of the corresponding bio-composite materials.

Biocomposites loaded with the raw and treated OWF were analyzed by thermogravimetry (ATG) in order to know the effect of chemical treatment of the load and the rate of reinforcement on the thermal stability of the material.

The results obtained by thermogravimetric analysis (ATG) are presented in Table 1 as a function of degradation temperature (Tmax).

The addition of 10%, 20%, 30% of untreated fibers significantly increases the degradation temperatures respectively by 5, 3, and 9 °C relative to the matrix alone. This effect is comparable with the treated fiber content, for which the temperature varies from 0, 3, and 6 °C with the addition of 10%, 20%, 30% of treated fibers, respectively. In our case, neither the addition of fiber content nor the chemical treatment significantly affected the thermal stability of the PP/OWF composite.

The literature [START_REF] Shukor | Effect of ammonium polyphosphate on flame retardancy, thermal stability and mechanical properties of alkali treated kenaf fiber filled PLA biocomposites[END_REF][START_REF] Manshor | Mechanical, thermal and morphological properties of durian skin fibre reinforced PLA biocomposites[END_REF] shows that generally the concentration of chemical treatment is the factor that significantly influences the thermal stability of composites because of its ability to remove hemicellulose, lignin and impurities from the fiber. This allows an improvement in the thermal stability of the fibers. However, the modification made to the fibers in our study is not sufficient to affect the thermal stability of the composite in its entirety

Conclusion

In this study we characterized the raw olive wood flour and chemically treated with silane. Indeed, the development of new bio composites requires a prior knowledge of the different properties and the composition of fibers.

The SEM micrographs of both types of flours revealed an increase in the roughness of the treated versus untreated fibers due to the removal of non-cellulosic constituents. The chemical analysis reveals that after chemical treatment, the hemicellulose and lignin contents decrease, which causes the increase of the cellulose content in the treated fibers. Secondly, we have characterized morphologically and thermally elaborate bio-composites. SEM images showed good adhesion between the matrix and treated or untreated fibers. Thus, the cohesion at the fiber/matrix interface seems to improve with the chemical treatment.

Thermogravimetric analysis has shown that the thermal stability of bio-composites does not change much for both treated and untreated fibers.

Fig. 1 .

 1 Fig. 1. SEM images of olive wood flour fibers: (a) untreated OWF, (b) OWF treated

Fig. 2 .

 2 Fig. 2. SEM images of composite: (a) PP/OWF, (b) PP/OWFT

Table 1 .

 1 Degradation temperatures of raw and processed PP/OWF composites.

	Composition Degradation temperature of PP/OWF Degradation temperature of PP/OWF
	PP/	466.61	466.61
	PP/OWF 10% 471.96	466.54
	PP/OWF 20% 469.41	469.82
	PP/OWF 30% 475.54	472.58