Bharat B Tripathi

Adrian Luca

Sambandam Baskar

François Coulouvrat

Régis Marchiano

Element centered

Keywords: Discontinous Galerkin, Shock capturing, Artificial viscosity, Nonlinear acoustics

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

One of the most spectacular features of nonlinear acoustics is the generation of shock waves along the propagation.

In this case, the speed c of finite amplitude sound waves is not strictly constant, even in homogeneous fluids. It is dependent on the wave instantaneous pressure amplitude p a . At first order, one has c ≈ c 0 + βp a /ρ 0 c 0 where c 0 is the speed of sound of waves of infinitely small amplitude, and ρ 0 is the medium density. The parameter β characterizes the medium nonlinearity and will be defined later on. Accordingly, the parts of the waves with the highest amplitude travel faster than those with smaller values. This results into a distortion of the waveform that cannot keep its profile unchanged. Starting, for instance, from a smooth sine wave, this distortion first leads to the steepening of parts of the wave profile, and ultimately to the formation of acoustic shock waves if the amplitude is high enough, or if the propagation distance is sufficiently long [START_REF] Lighthill | Waves in Fluids[END_REF]. The acoustic shock waves can also be produced directly by the source itself and persist along propagation. Atmospheric acoustic shock waves arise in many situations, like for instance the sonic boom from a supersonic aircraft [START_REF] Maglieri | Sonic boom[END_REF], bolide hypersonic atmospheric entries [START_REF] Henneton | Numerical Simulation of Sonic Boom from Hypersonic Meteoroids[END_REF], or the so-called "buzz saw noise" due to the supersonic rotation of the tip of fan blades [START_REF] Mcalpine | On the prediction of "buzz-saw" noise in acoustically lined aero-engine inlet ducts[END_REF][START_REF] Fernando | Nonlinear waves and shocks in a rigid acoustical guide[END_REF]. Lithotripsy [START_REF] Averkiou | Modeling of an electrohydraulic lithotripter with the KZK equation[END_REF] uses ultrasonic shock waves propagating in biological tissues to fragment kidney stones, while new medical applications are based on HIFU (High Intensity Focused Ultrasound) [START_REF] Canney | Shock-induced heating and millisecond boiling in gels and tissue due to high intensity focused ultrasound[END_REF]. Many of these examples involve shock formation, reflection, diffraction or focusing in propagation media with complex geometries. Most of the numerical methods developed and well studied in the field of nonlinear acoustics are for models that govern "one-way" wave propagation. The most popular model is the KZK equation [START_REF] Kuznetsov | Equation of nonlinear acoustics[END_REF] which is a parabolic approximation of the nonlinear wave equation. Even if this model is valid for many situations, it suffers from an angular limitation (±18 o ). Several improvements have been proposed to go beyond the parabolic approximation (see [START_REF] Christopher | Modeling the Dornier HM3 lithotripter[END_REF] or [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF]). These methods are highly efficient for modeling long distance propagation, for instance, sonic boom in atmosphere. But the intrinsic limitation of these methods, is the inability to model back-scattering/reflections. This problem can only be addressed by the "full-wave" methods which model the nonlinear wave propagation in all the directions. Some of the popular full-wave methods are: "Fullwave" by Pinton et al. [START_REF] Pinton | Nonlinear and von neumann reflection of elastic shock waves in soft solids[END_REF] using a finite difference approach, "K-wave" by Treeby et al. [START_REF] Treeby | Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method[END_REF] using the k-space method, Fourier Continuation method [START_REF] Albin | Fourier continuation methods for high-fidelity simulation of nonlinear acoustic beams[END_REF], these methods are based on Cartesian meshes and can manage weak nonlinear propagation (only the first few harmonics of the wave spectrum). A numerical method handling complex geometries by using unstructured mesh and acoustic shock waves together is still a challenge, which is the aim of this work.

To solve these two issues, we propose to use the discontinuous Galerkin method (DGM) for the propagation of acoustic shock waves on an unstructured mesh. DGM is a kind of hybrid between the finite element and the finite volume methods. Hence, it gives the advantage of local high order accuracy. Moreover, it is suitable to handle complex geometries thanks to the use of unstructured meshes. Another advantage is its parallelization potential similar to the finite volume methods, which further increases the efficiency of the method. The DGM was first proposed by Reed and Hill [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] for solving a steady-state neutron transport equation, with its analysis provided by Lesaint and Raviart [START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF]. At present the DGM is widely applied to many areas [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF]. In acoustics, it has been mainly used for linear acoustics [START_REF] Kopriva | Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers[END_REF], aeroacoustics [START_REF] Toulopoulos | High-Order Discontinuous Galerkin Discretizations for Computational Aeroacoustics in Complex Domains[END_REF][START_REF] Gabard | Discontinuous Galerkin methods with plane waves for the displacement-based acoustic equation[END_REF][START_REF] Gabard | Discontinuous Galerkin methods with plane waves for time-harmonic problems[END_REF], propagation at the interface between moving media and isotropic solids [START_REF] Luca | Numerical simulation of transit-time ultrasonic flowmeters by a direct approach[END_REF], and weakly nonlinear propagation in solids [START_REF] Bou Matar | A nodal discontinuous Galerkin finite element method for nonlinear elastic wave propagation[END_REF]. To our knowledge, DGM has not been used for propagation of acoustic shock waves in complex geometries.

Nevertheless, using DGM for acoustic shock wave propagation requires to pay attention to the shock handling.

One of the main difficulties in shock handling is to avoid Gibbs oscillations. In order to tackle this problem for high-order schemes, many tools are available in the literature based on slope limiters, filters, and artificial viscosity.

A lot of work has been done in the area of slope limiters by Cockburn and co-workers [START_REF] Cockburn | TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws II : General Framework[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems[END_REF][START_REF] Cockburn | The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case[END_REF][START_REF] Cockburn | The Runge -Kutta Discontinuous Galerkin Method for Conservation Laws V[END_REF] and further extended by Biswas et al. [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF] and Burbeau et al. [START_REF] Burbeau | A Problem-Independent Limiter for High-Order Runge-Kutta Discontinuous Galerkin Methods[END_REF]. Recent work in the direction of slope limiters are for high-order WENO methods like in ADER-WENO approach [START_REF] Dumbser | A posteriori subcell limiting of the discontinuous galerkin finite element method for hyperbolic conservation laws[END_REF] with a subcell shock limiter, and Zhu et al. [START_REF] Zhu | Runge-kutta discontinuous galerkin method with a simple and compact hermite weno limiter on unstructured meshes[END_REF] proposed the use of Hermite WENO limiter for DGM by reconstructing the polynomial using the original DG solution from the neighboring cells. Nevertheless, slope limiters are not the optimal choice for high-order methods as they flatten the smooth extrema so that accuracy is lost. The use of modal filters in spectral methods was initially proposed by Majda et al. [START_REF] Majda | The Fourier Method for Nonsmooth Initial Data[END_REF]. Relevant works on filtering can be found in [START_REF] Hesthaven | Filtering in Legendre spectral methods[END_REF][START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF]. Here also, the problem could be the flattening of smooth extrema, if filters are applied globally. Coupling with shock-capturing tools (based on the spectral modes) is not efficient as filtering destroys the natural evolution of the modes. Recent application of spectral filtering for discontinuous Galerkin methods is achieved by Meister et al. [START_REF] Meister | Application of Spectral Filtering to Discontinuous Galerkin Methods on Triangulations[END_REF], based on the idea of spectral viscosity initially proposed by Maday et al. [START_REF] Maday | Legendre Pseudospectral Viscosity Method for Nonlinear Conservation Laws[END_REF] in coupling with sub-cell shock detection designed by Persson and Peraire [START_REF] Persson | Sub-Cell Shock Capturing for Discontinuous Galerkin Methods[END_REF]; finally the solution is post-processed using what is called as the digital total variation filtering.

We choose here the method of artificial viscosity introduced by von Neumann and Richtmyer [START_REF] Neumann | A Method for the Numerical Calculation of Hydrodynamic Shocks[END_REF] as an efficient method of shock capturing for acoustic shock waves. This approach has been proposed and used with success also to stabilize the Euler equations on arbitrary geometric domains by Jameson et al. [START_REF] Jameson | Numerical solutions of the euler equations by finite volume methods using runge-kutta time-stepping schemes[END_REF]. Relying on this approach, Hughes and co-workers [START_REF] Brooks | Streamline Upwind/Petrov-Galerkin Formulations For Convection Dominated Flows With Particular Emphasis On The Incompressible Navier-Stokes Equations[END_REF][START_REF] Hughes | A New Finite Element Formulation For Computational Fluid Dynamics: I. Symmetric Forms Of The Compressible Euler And Navier-Stokes Equations And The Second Law Of Thermodynamics[END_REF][START_REF] Hughes | A New Finite Element Formulation For Computational Fluid Dynamics: II. Beyond SUPG[END_REF][START_REF] Hughes | A New Finite Element Formulation For Computational Fluid Dynamics: III. The Generalized Streamline Operator For Multidimensional Advective-Diffusive Systems[END_REF][START_REF] Hughes | A New Finite Element Formulation For Computational Fluid Dynamics: IV. A Discontinuity-Capturing Operator For Multidimensional Advective-Diffusive Systems[END_REF] introduced the streamline diffusion method which was successful in damping the oscillations. For DGM in past few years, the local artificial viscosity method has gained significant importance. It is possible to couple it with the sub-cell shock detection algorithm, which is particularly important for unstructured meshes. Persson and Peraire [START_REF] Persson | Sub-Cell Shock Capturing for Discontinuous Galerkin Methods[END_REF] implemented this idea of sub-cell shock detection using the highest-order coefficients in an orthonormal representation of the solution. Once a shock is sensed in a particular element, a piecewise-constant artificial viscosity is introduced depending on the mesh and the solution. This local approach makes it highly adaptable for parallelization, which is of key importance for DG implementation. The problem with this method are the jump discontinuities occurring in the viscosity map associated to the solution, which induce oscillations at the element boundaries. To dampen these oscillations, Barter and Darmofal [START_REF] Barter | Shock capturing with PDE-based artificial viscosity for DGFEM: Part I. Formulation[END_REF] used a smooth artificial viscosity resulting from a diffusion equation. They combined hybrid meshes (structured ones near the shock and unstructured ones otherwise) for solving compressible Navier-Stokes equations. They also used an inter-element jump indicator proposed by Dolejsi et al. [START_REF] Dolejsi | On some aspects of the discontinuous Galerkin finite element method for conservation laws[END_REF]. Klockner et al. [START_REF] Klöckner | Viscous Shock Capturing in a Time-Explicit Discontinuous Galerkin Method[END_REF] extended the work of Persson, actually trying to smoothen the viscosity through its linear interpolation in the neighboring elements, while choosing the parameters more methodically. Alternatively, Guermond et al. [START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF] proposed to use the entropy viscosity, for which the viscosity coefficient depends on the entropy production. The viscosity is then further smoothened, as theoretically justified by Bonito et al. [START_REF] Bonito | Stability Analysis of Explicit Entropy Viscosity Methods for Non-Linear Scalar Conservation Equations[END_REF]. Recent approaches of shock capturing by means of residual-based artificial viscosity are implemented by Kurganov et al. [START_REF] Kurganov | New adaptive artificial viscosity method for hyperbolic systems of conservation laws[END_REF] based on the concept of weak local residual error [START_REF] Karni | A Smoothness Indicator for Adaptive Algorithms for Hyperbolic Systems[END_REF][START_REF] Karni | Local error analysis for approximate solutions of hyperbolic conservation laws[END_REF]. Hartmann and Houston [START_REF] Hartmann | Adaptive Discontinuous Galerkin Finite Element Methods for the Compressible Euler Equations[END_REF][START_REF] Hartmann | Adaptive discontinuous Galerkin methods with shock-capturing for the compressible Navier-Stokes equations[END_REF] relied on the interior penalty DGM for compressible Navier-Stokes while capturing the shocks using local residuals. Reisner et al. [START_REF] Reisner | A space-time smooth artificial viscosity method for nonlinear conservation laws[END_REF] worked in a direction similar to Barter et al. [START_REF] Barter | Shock capturing with PDE-based artificial viscosity for DGFEM: Part I. Formulation[END_REF] by modeling the viscosity coefficient through a linear scalar-diffusion equation where a gradient-based source term is introduced to trigger the viscosity. Reisner et al. [START_REF] Reisner | A space-time smooth artificial viscosity method for nonlinear conservation laws[END_REF] clearly outlined as a perspective that: "In the future, the gradient-based source term used in the current implementation of the C-method may be combined with a noise-indicator that turns off the current gradient-based source term when it is not needed." Therefore, a detailed study of the combined effect of gradient and noice indicators is essential step forward. Moreover, none of these works are directly related to nonlinear acoustics.

In the present work, we derive nonlinear equations of acoustics under a conservative form. Then, we propose a new sub-cell shock capturing for the DG method tool based on the gradient and the noise sensors, which we call as the shock sensor (SS). A detailed study of the combined effect of the gradient and noise indicators is presented, highlighting the motivation behind the development of this very sensitive and robust shock sensor. Though, it still depends on global parameters which is the topic of further research. Once an acoustic shock is sensed in a particular element, an element centered smooth artificial viscosity (ECSAV) is introduced there locally into the system. Note that, as we do not model the viscosity coefficient using a PDE, this makes it comparatively simpler to implement while still having only three empirical parameters. The validity of the method is demonstrated throughout different problems of nonlinear acoustics. Planar wave propagation is used to study the formation of N-wave, along with the illustration on the motivation behind the shock sensor. Then, propagation of cylindrical shock waves is studied to demonstrate the ability of the method to handle geometrical spreading and nonlinear effects. The last test deals with the nonlinear radiation of a piston source. It demonstrates the ability of the code to handle diffraction combined with nonlinearity.

Theoretical Model

Using the fundamental equations of fluid dynamics [START_REF] Landau | Fluid mechanics[END_REF], a first-order system of equations is derived, written in conservative form relevant for nonlinear propagation of acoustic waves. It is equivalent to the Kuznetsov equation [START_REF] Kuznetsov | Equation of nonlinear acoustics[END_REF] which is a second-order, nonlinear scalar wave equation for potential. The fluid is assumed to be homogeneous in composition and quiescent i.e., therefore its ambient density and pressure are uniform and there is no ambient flow. Moreover, it is assumed to be lossless i.e., the viscous and thermal effects are negligible. The state variables, density, velocity, pressure are denoted as ρ, v = (u, v), p, respectively.

State variables can be written as the sum of the ambient state and the acoustic perturbation [START_REF] Pierce | Nonlinear Acoustics[END_REF][START_REF] Coulouvrat | On the equations of nonlinear acoustics[END_REF][START_REF] Hamilton | Nonlinear Acoustics[END_REF]: one has p(x, t) = p 0 + p a (x, t) and ρ(x, t) = ρ 0 + ρ a (x, t), where the subscripts 0 and a indicate the ambient quantities and acoustic perturbations, respectively. We identify v = v a as the medium is quiescent. Substituting these expressions in the conservation laws of fluid dynamics, and retaining only terms up to second order whereas neglecting the cubic O(ρ 3 a ) and higher order terms, one gets a system of three, first-order equations exact up to the quadratic nonlinear terms included. Its non-dimensionalized formulation is

∂ ρa ∂ t + ∂ ∂ x (1 + ρa ) ūa + ∂ ∂ȳ (1 + ρa ) va = O( 2 ) (1)
∂ ∂ t (1 + ρa ) ūa + ∂ ∂ x ū2 a + ρa + B 2A ρ2 a + ∂ ∂ȳ (ū a va ) = O( 2 ) (2)
∂ ∂ t (1 + ρa ) va + ∂ ∂ x (v a ūa ) + ∂ ∂ȳ v2 a + ρa + B 2A ρ2 a = O( 2 ). (3)
The acoustic pressure p a has been eliminated by means of an expansion of the state equation at the same order

pa = ρa + B 2A ρ2 a + O( 2 ). ( 4
)
The non-dimensionalized variables are: pa = for the state variables. The spatial and temporal variables are transformed using

x = x L , ȳ = y L
, with L = c 0 ω 0 and t = ω 0 t, respectively. Here, ω 0 is the angular frequency of the initial wave profile and 2πL is the corresponding wavelength. The acoustic Mach number , in front of all nonlinear quadratic terms, is defined as:

= max x {u a } c 0 = p m a ρ 0 c 2 0 . ( 5
)
For most of the applications cited in introduction, the acoustic Mach number rarely exceeds 10 -2 . The assumption of weak nonlinearity and the neglecting of cubic and higher order terms is therefore fully justified.

The ratio B/A is the fluid nonlinear parameter, measuring the quadratic nonlinearity of the state equation. For a perfect gas, it is equal to (γ -1)/2 where γ is the ratio of specific heats. One also has β = 1 + B/2A. This ratio plays a key role in nonlinear acoustics. Values for various non-gaseous media are collected in [START_REF] Beyer | Nonlinear Acoustics[END_REF][START_REF] Hamilton | Nonlinear Acoustics[END_REF]. Note that Sparrow

and Raspet [START_REF] Sparrow | A numerical method for general finite amplitude wave propagation in two dimensions and its application to spark pulses[END_REF] used the system (1)-( 4) for the first time to simulate directly the propagation of nonlinear acoustic waves. This set of equations can be used to derive the inviscid Kuznetsov equation with just algebraic manipulations and without any additional assumption. In 1D frame, assuming a one-way propagation, the system (1)-( 4) can also be reduced to the inviscid Burgers equation which is a nonlinear, scalar and 1D model.

The method of artificial viscosity involves a parabolic regularization of hyperbolic conservation laws, by adding one or several dissipative terms on the right hand side of each conservation law. The choice of the dissipative term(s) is far from unique. Introducing the viscosity coefficient η(x), the most intuitive choice would be η ∂ 2 ∂x 2 . However, the choice of the parabolic term ∂ ∂x η ∂ ∂x is more consistent with the form of conservation equations with a variable viscosity, and will help to develop a first order system of equations as a prerequisite for the DGM implementation.

Therefore, the parabolic-regularization of the above dimensionless system of equations ( 1), ( 2), (3) becomes

          ρa (1 + ρa )ū a (1 + ρa )v a           t + ∇ •            (1 + ρa )ū a ū2 a + ρa + B 2A ρ2 a ūa va (1 + ρa )v a ūa va v2 a + ρa + B 2A ρ2 a            = ∇ •              η 1 ∂ ∂x ρa η 2 ∂ ∂x ((1 + ρa ) ūa ) η 3 ∂ ∂x ((1 + ρa ) va ) η 1 ∂ ∂y ρa η 2 ∂ ∂y ((1 + ρa ) ūa ) η 3 ∂ ∂y ((1 + ρa ) va )              . (6)
Here, the viscosity coefficients η i = η i (x, y, t), i = 1, 2, 3 are functions of space variables and time, and are non-zero only over a small neighborhood of each shock. Details about the viscosity coefficients are given in section 3. The above system can be written in a generic conservative form for variables q 1 = ρa , q 2 = (1 + ρa )ū a and

q 3 = (1 + ρa )v a as ∂q m ∂t + ∂ f m ∂x + ∂g m ∂y = ∂ ∂x η m ∂ ∂x q m + ∂ ∂y η m ∂ ∂y q m ; for m = 1, 2, 3. (7)
Flux terms f m and g m are obviously deduced from [START_REF] Bonito | Stability Analysis of Explicit Entropy Viscosity Methods for Non-Linear Scalar Conservation Equations[END_REF].

The discontinuous Galerkin formulation of such a convective-diffusive system is achieved using the so-called local discontinuous Galerkin method, initially proposed by Bassi and Rebay [START_REF] Bassi | A High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations[END_REF] for compressible Navier-Stokes equations.

It was further studied by Cockburn and Shu [START_REF] Cockburn | The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems[END_REF]. It involves the splitting of the convective-diffusive equation into a system of first-order equations. Its weak formulation is written using Dubiner basis [START_REF] Dubiner | Spectral Methods on Triangles and other Domains[END_REF] (constructed using 2D Jacobi polynomials) as test functions. The two-dimensional complex domains are discretized using an unstructured mesh made up of triangular elements. The key feature of the discontinuous Galerkin method is the element-centric approach implying that almost all the computations are done independently within each element. The local Lax-Friedrichs flux is used for the connectivity within different elements. This gives the motivation to parallelize the computation within each element. This is done using the Nvidia graphic cards with the pycuda environment [START_REF] Klöckner | Viscous Shock Capturing in a Time-Explicit Discontinuous Galerkin Method[END_REF].

The temporal advancement is achieved using the low storage explicit fourth-order Runge-Kutta method. Additional details about the method and its implementation can be found in the textbook [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF].

Shock Management

As mentioned above, the problem of spurious oscillations is tackled by localizing the regions with non-physical oscillations and introducing there the appropriate amount of viscosity (η i (x), for i = 1, 2, 3) into the system (6).

Shock Sensor

In order to define the shock sensor, the primitive variables ρ a , u a , v a are easily computed from the conserved quantities q 1 , q 2 , q 3 . In discontinuous Galerkin method the approximate solution can be represented in two different forms, namely, the nodal one and the modal one. For the management of shock, the modal solution is of most importance. Indeed, with motivation from the work of Persson and Peraire [START_REF] Persson | Sub-Cell Shock Capturing for Discontinuous Galerkin Methods[END_REF], our new sub-cell shock detection tool is developed along similar lines using the coefficients of the spectral solution (modal solution). The interpolating polynomials ψ i (ξ) being two dimensional, here, ξ is the coordinate system in the reference element. The modal solution can be written using two indices as

ρ a (ξ, t) = N i=0 N-i j=0 ( ρa ) i j (t)ψ i j (ξ). ( 8
)
where i and j denote the order of the interpolating polynomial ψ i j (ξ) with respect to ξ and η, respectively. ψ i j (ξ) is the 2D orthonormal Dubiner basis [START_REF] Dubiner | Spectral Methods on Triangles and other Domains[END_REF][START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF] with the scope in a particular element only. With N as the order of interpolating polynomial, the total number of points in one element turns out to be N p = 1 . . . (N + 1)(N + 2)/2, see [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF] for further details. The coefficients ( ρa ) i j (t) of the basis functions are often referred as the modes of the DG solution. They play a key role in our method for sensing the shock. We intend to exploit the modes of k th element which is ( ρa ) k i j (t), this helps us estimate the shape of the solution vector in each and every element of the mesh.

Our shock sensor is based on the linear components i.e., the coefficients of ψ 01 (ξ) and ψ 10 (ξ), and the highest order components of the modal solution of order N, i.e.the coefficients of ψ 0N (ξ) and ψ N0 (ξ). The reason for this choice of coefficients is that the linear component measures the gradient of the waveform, whereas the highest order coefficients indicate the presence of spurious oscillations associated to the nonlinear generation of higher harmonics. In various tests not reproduced here, the choice of only linear coefficients was insufficient as they were not sensitive to spurious oscillations. Different combinations of ( ρa ) i j (t) were considered but no improvement was evident, implying that the first and last modes are the first to respond to any change in the waveform in nonlinear acoustics. Moreover, this idea was also expressed by Reisner et al. [START_REF] Reisner | A space-time smooth artificial viscosity method for nonlinear conservation laws[END_REF] as an important perspective for future research. Further details of this choice are explained in section 5.1.1. Now, we define the Shock Sensor (SS) in the k th element associated to any physical variable (ρ a , u a , v a ) of the system (6), as

(S S ) k ρ a (t) = (S S 1) k ρ a (t) max k {(S S 1) k ρ a (t)} + (S S N) k ρ a (t) max k {(S S N) k ρ a (t)} , (9)
(S S ) k u a ,v a (t) = (S S 1) k u a ,v a (t) max{max k {(S S 1) k u a (t)}, max k {(S S 1) k v a (t)}} + (S S N) k u a ,v a (t) max{max k {(S S N) k u a (t)}, max k {(S S N) k v a (t)}} , (10)
The notation (S S ) k u a ,v a (t) implies that it is either (S S ) k u a (t) or (S S ) k v a (t). It is important to mention here that, the maximum is calculated over all the elements. Further, (S S 1) f a is the first-order sensor of variable f a equal to either ρ a , u a or v a and defined as

(S S 1) k f a (t) =| ( fa ) k 01 (t) | + | ( fa ) k 10 (t) |, (11)
and (S S N) f a is similarly the Nth-order sensor of variable f a

(S S N) k f a (t) =| ( fa ) k 0N (t) | + | ( fa ) k N0 (t) | . (12)
Here ( fa ) k i j (t) are the coefficients of the modal solution (8) in the k th element. For the sake of brevity we denote the shock sensors (S S ) k m (t) with m = 1 for (S S ) k ρ a (t), m = 2 for (S S ) k u a (t) and m = 3 for (S S ) k v a (t). Note that the L ∞ -norm is used to calculate the shock sensor as the objective is to give the maximum importance to the element with the maximum gradient or maximum oscillations. Also, it is important to note that ( 11)-( 12) gives a relative value, and this choice works efficiently for sensing acoustic shock waves when represented using the Dubiner orthonormal basis.

Once shock sensors are calculated for each element, the need of viscosity η m in the k th element is checked provided the condition

(S S ) k m (t) ≥ max k {(S S ) k m (t)} α 1 , k = 1, • • • , K, (13)
is satisfied. Here α 1 is a user-given parameter quantifying the minimum value of (S S ) k m (t) above which a region of high gradient is sensed for m component of the solution at time t in element k. When satisfied, the corresponding element is tagged as an infected element. (S S ) k m (t) ∈ [0, 1] with the highest value representing the element with highest contribution from linear and highest order modes. This shock check condition helps identifying the elements needing artificial viscosity. All the numerical results in this paper are computed using α 1 = 10. According to our numerous tests α 1 = 10 gives a sufficiently broad range of elements in and around the region of shock, and therefore remains fixed. In order to reduce the spread of the viscosity around the shock, the value of the parameter α 1 can be reduced or vice versa. A physical explanation for our choice of shock sensor will be explained by means of various numerical experiments in section 5.1.1.

It is important to sense regions of high gradient at all time steps irrespective of the presence or not of a shock.

Therefore it is important to calibrate the amount of viscosity to be introduced in the domain at each time step. This is done using the Gradient Factor (GF), defined as

GF(t) = exp           max k {(S S 1) k ρ a (t)} max k {(S S 1) k ρ a (0)} -1           . (14)
The function GF(t) measures the steepening of the wave profile with respect to the initial condition. Note that the GF(t) is chosen as an exponential function, because the evolution of the modes in time tends to behave exponentially.

This keeps the artificial viscosity sufficiently low before the acoustic shock is formed so that the smooth profile is dissipated as less as possible. In order to keep the gradient factor GF(t) under control, it is tapped by an upper limit as

GF(t) ≤ α 2 , (15)
where α 2 is a user given parameter. All the numerical results in this paper are computed using α 2 = 20.

Element Centered Smooth Artificial Viscosity

Once shocks are localized (i.e. the infected elements are flagged), an appropriate amount of smooth artificial viscosity is introduced in the k th infected element as a two-dimensional Gaussian distribution

η k m (x k i ) = (η 0 ) k m exp         -       x k i -(x 0 ) k (σ 0 ) k       2 -       y k i -(y 0 ) k (σ 0 ) k       2         , i = 1, ..., N p . (16)
Here

x k i = (x k i , y k i ), i = 1 . . . N p ,
are the points in the k th element. We label [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems[END_REF] as the element centered smooth artificial viscosity (ECSAV). The ECSAV reaches its maximum in the element at position (x k 0 , y k 0 ), naturally chosen as the centroid of the considered k th element. The width σ k 0 of the ECSAV beyond which artificial viscosity exponentially decays is defined as the circumradius of the respective infected element. Finally, the choice of the amplitude of ECSAV (η 0 ) k m in the k th infected element is very crucial because it has to be nonzero only around the shock and zero elsewhere.

It is defined as

(η 0 ) k m (t) =            α 3 .GF(t).(S S ) k m (t) if (S S ) k m (t) ≥ max k {(S S ) k m (t)} α 1 0 otherwise. ( 17
) Numerical Paramters Values α 1 10 α 2 20 α 3 O(2 × 10 -2 ) Table 1: Numerical Parameters
We recall here that the inequality ( 13) is used as an indicator function to turn on/off the viscosity. Parameter α 3 is empirically chosen. At present, we do not have an expression/bound for this parameter. Nevertheless, based on our experience relying on numerous numerical tests, we propose

α 3 ≈ O(2 × 10 -2 ), (18)
where is the acoustic Mach number. Note that proportionality of numerical viscosity to wave amplitude is nevertheless coherent with the expansion order of the system (1)-(4). All the numerical parameters in this implementation are summarized in the Table 1. It is important to mention that the use of smooth viscosity instead of a piecewise constant viscosity gives a significantly better solution. Sudden inter-element jumps in the viscosity induce oscillations that will propagate and may even lead to new unphysical tiny shocks due to cumulative nonlinearity. Further details are given in the forthcoming sections, and especially in subsection 5.1.4, different ways of implementing ECSAV are discussed so as to make viscosity as smooth as possible.

One-dimensional validation tests

First, the method is implemented in 1D. This preliminary stage allows to focus on the design of the shock sensor and to validate the method carefully by comparing the numerical results to a quasi-analytical solution [START_REF] Hayes | Sonic Boom Propagation In A Stratified Atmosphere, With Computer Program[END_REF][START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF]. In 1D, the system of equations ( 6) can be simplified into the Burgers equation in retarded time :

∂p ∂σ - ∂ ∂τ p 2 2 = ∂ ∂τ η(τ, σ) ∂p ∂τ , (19)
with the following dimensionless variables:

σ = x L sh , τ = ω 0 t - x c 0 . (20)
The characteristic length L sh is the shock length i.e., the distance required for an initially sine wave to become an acoustic shock wave. It is given by

L sh = 1 β k , (21)
where k = ω 0 /c 0 is the wavenumber. We choose as initial condition a single sine wave period

p(0, τ) = sin(π(τ -0.05)) if -2π ≤ τ -0.05 ≤ 2π 0 otherwise. . (22)
Quasi-analytical solution of Burgers equation is given by the so-called Burgers-Hayes method [START_REF] Hayes | Sonic Boom Propagation In A Stratified Atmosphere, With Computer Program[END_REF][START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF], which consists in Poisson's implicit solution expressed for potential rather than pressure. In case of shock formation, Poisson's solution gets multivalued. Then physically admissible solution for potential is the maximum value, in order to satisfy the entropy condition. Analytical Burgers-Hayes solution will be referred as 'Quasi-Analytical' in the following. The numerical parameters used for DGM are chosen as τ ∈ [-2π, 2π], σ ∈ [0, 1.5].

Shock Sensor Vs Smoothness Indicator

The efficiency of our shock sensor (SS) is highlighted over the Persson's smoothness indicator (SI) in this section.

Figure 1 shows the comparison of the two using the same unstabilized DG solution computed (left subplot) slightly beyond shock formation ( σ=1.07) with no viscosity. In Persson's approach (center subplot), the value of S I must be greater than (note S I is a negative axis) the value of the black dotted line in order to turn on the viscosity. The method responds well but at the central element. However, the value of S S (right subplot) is obviously well above the threshold for all three central elements located around the shock, and much below elsewhere. Therefore, SS is more sensitive to such numerical oscillations, which makes it more susceptible to shocks than to smoothness.

Smooth Artificial Viscosity Vs Piecewise constant viscosity

Figure 2 shows the comparison of the solutions obtained using the different methods of artificial viscosity allocation. Persson's algorithm combining smoothness indicator and constant viscosity per element is denoted by 'SI+CV'.

The second approach is a kind of hybrid of Persson's approach and ours, using SS to detect the shock but introducing a constant viscosity (equal to the maximum of ECSAV) when ( 13) is satisfied. This is denoted by 'SS + max(ECSAV)'.

Lastly, Figure 2 displays the full method combining SS shock sensor and ECSAV smooth viscosity. These three ways of introducing artificial viscosity are compared with quasi-analytical solution. Computation is shown at position σ = 1.455 close to maximum shock amplitude, and is performed with K = 50 elements and fourth polynomial order approximation (N = 4). The overall wave profile (left upper subplot) shows that all methods well localize the shock and simulate the waveform away from it. The artificial viscosity map (left lower subplot) nevertheless shows for ECSAV a smooth viscosity distribution well localized around the shock, while the two other methods with piecewise constant viscosity necessarily introduce some discontinuities and a more widespread distribution. From the zoom-in of the lower part of the discontinuity (right subplot), the full 'SS+ECSAV' method shows a very smooth solution with no oscillations at all, thanks to the smooth viscosity distribution, while other methods are not as efficient and lead to oscillations just before shock. Same conclusions could be drawn by examining the solution just after shock.

Full propagation of this case is shown in Figure 3, i.e., before and after shock formation till two shock lengths. It is important to observe that even before the shock is formed, ECSAV is present with a very small amplitude.

Handling of Multiple Shocks and Shock Merging

We now consider an initial condition leading to multiple moving shocks, chosen as Figure 4 shows the initial condition (left column) and the waveform after propagating over two shock lengths (right column). The three initial shocks are now moving, the two first ones leftwards, the last one rightwards, which leads to a lengthening of the waveform. All these features are well simulated, the only approximation being a slight spreading of the numerical solution resulting from the introduction of the stabilizing artificial viscosity. One can observe that there is no viscosity initially (center-left plot), however the viscosity allocation during the propagation (center-right plot) is proportional to the shock strength. Shock sensors (two lower plots) and viscosity maps well follow the shock displacement. The two first shocks keep perfectly separated from one another, despite the fact they are quite close from one another. Figure 5 shows the same figures (waveforms, viscosity maps, shock sensors) after traveling over three shock lengths (left) and five shock lengths (right). At three shocks lengths, the two first shocks are about to merge, the shock sensor is spread over 3 elements only instead of 2, and the two viscosity maps now slightly interfere locally. At five shock lengths, the two head shocks have merged, wave profile has stabilized into an N-wave, with two shocks moving in opposite directions. These ones are perfectly well-tracked over one single element by Shock Sensor.

p(0, τ) =                    - (τ + 0.5π) π , if -π ≤ τ ≤ -0.7π - τ π , if -0.7π < τ ≤ π 0, otherwise (23)

Two-dimensional validation tests

In this section, three numerical tests are presented for the validation of the 2D method. First, the numerical solution for the propagation of a plane wave is compared to the quasi-analytical solution of the inviscid Burgers equation for plane waves [START_REF] Hayes | Sonic Boom Propagation In A Stratified Atmosphere, With Computer Program[END_REF][START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF]. This test is also used to underline the motivations behind the shock sensor, along with the different approaches of implementation of ECSAV on an unstructured mesh. The second test deals with the nonlinear propagation of cylindrical waves. The numerical results are compared to a quasi-analytical solution of the Burgers equation for cylindrical waves [START_REF] Enflo | Saturation of nonlinear spherical and cylindrical sound waves[END_REF]. The last test is a comparison between the numerical results and results computed by another method: a one-way method called HOWARD (Heterogeneous One-Way Approximation for the Resolution of Diffraction) presented in Dagrau et al. [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF] for simulating nonlinear radiation from a piston-transducer.

Propagation of Nonlinear Plane Waves

In this part, a simulation of a plane-wave along x-axis is done on a 2D unstructured mesh (boundary conditions are rigid walls). For this configuration it is possible to compare the numerical results to the quasi-analytical solution of the Burgers equation, as the initial and the boundary conditions are consistent with the assumptions underlying the Burgers equation.

Numerical motivations behind the shock sensor

To illustrate the interest of the composite shock sensor (between the first and the highest orders), its behavior is analyzed for the case of a plane wave propagating in a 2D unstructured mesh. The importance of SS1 is highlighted in subsection 5.1.2, and that of SSN is discussed in subsection 5.1.3. The computational domain is rectangular ( x, ȳ) = [-21, 21] × [-4.2, 4.2]. The unstructured mesh consists in 1308 triangular elements and the order of polynomial approximation is 8. The initial conditions are set to correspond to a plane wave propagating toward + x and the waveform is an inverted sine-period:

pa ( x, ȳ, t = 0) = -sin ( x) ; if | x| ≤ π, ūa ( x, ȳ, t = 0) = pa ( x, ȳ, t = 0), va ( x, ȳ, t = 0) = 0. ( 25
)
The acoustic Mach number is set equal to = 2.2 × 10 -4 . The numerical parameters related to ECSAV taken for this configuration are: α 1 = 10, α 2 = 20, and α 3 = 6 × 10 -6 .

First-Order Contribution to the Shock Sensor

Figure 6 presents three different shock sensors corresponding to ρ a , along with a zoom-in of the plot over x-axis of p a . All these plots are made after the propagation over around three shock lengths.

The subplot-(a) shows (S S ) k ρ a (t), k = 1, • • • , K, when only the linear contribution of the modal solution ( 8) is considered to construct the shock sensor. Recall that the S S k ρ a is constant in each element as it is computed from the modal solution of that element, and this constant is assigned to each node in that particular element. In this case, the shock sensor takes the form:

(S S ) k ρ a (t) = 2(S S 1) k ρ a (t) max k {(S S 1) k ρ a (t)} (S S ) k u a ,v a (t) = 2(S S 1) k u a ,v a (t) max{max k {(S S 1) k u a (t)}, max k {(S S 1) k v a (t)}} . ( 26
)
The subplot-(b) shows (S S ) k ρ a (t), k = 1, • • • , K, when only the highest-order contribution of the modal solution ( 8) is considered to construct the shock sensor. In this case, the shock sensor takes the form: Lastly, the subplot-(c) shows (S S ) k ρ a (t), k = 1, • • • , K, with the actual definition of S S as given by ( 9)- [START_REF] Burbeau | A Problem-Independent Limiter for High-Order Runge-Kutta Discontinuous Galerkin Methods[END_REF]. Note, the factor 2 in ( 26) and ( 27) is due to the purpose of normalization.

(S S ) k ρ a (t) = 2(S S N) k ρ a (t) max k {(S S N) k ρ a (t)} (S S ) k u a ,v a (t) = 2(S S N) k u a ,v a (t) max{max k {(S S N) k u a (t)}, max k {(S S N) k v a (t)}} . ( 27
)
After propagation over around three times the shock length, the sine-period is transformed into a sawtooth waveform due to nonlinear effects. Shock is clearly located by the first order shock sensor (26) as evident from Figure 6-(a). Also, the value of the shock sensor is very high around the shock, thus inducing as expected a high viscosity introduced into the system to suppress the oscillations.

In the case of a highest-order sensor [START_REF] Gabard | Discontinuous Galerkin methods with plane waves for time-harmonic problems[END_REF], the shock front is not as clearly visible on Figure 6-(b), and the value of the shock sensor is not as high as previously in the elements around the shock. Consequently, the viscosity introduced is much smaller and nonuniform around the shock. Therefore, there remains a possibility of spurious oscillations due to insufficient dissipation.

In the case of full shock sensor S S , a clear shock front is once again captured by the shock sensor as evident from Figure 6-(c). This case is similar to the first order sensor, with high values of shock sensor implying the sufficiently high viscosity required for dissipating spurious oscillations.

All the above remarks are supported by the resulting axial pressure pa zoomed around the shock and shown in Figure 6-(d). As expected pa corresponding to SSN (green) keeps some unphysical spurious oscillations, whereas the other two curves corresponding to the first-order sensor and full SS display no oscillations at all and are close to one another. It can therefore be concluded that highest-order sensor SSN is not sufficient alone to manage shocks, at least in the framework of weak acoustic shock waves on unstructured mesh. On the other hand, first-order sensor SS1 seems to be the key ingredient for tackling such situations. However, there are some numerical artifacts which the sensor SSN is capable of controlling, as will be discussed in the next subsection.

Highest-Order Contribution to the Shock Sensor

Since we are simulating the propagation of a plane wave along the x-axis in a 2D domain, va should remain zero.

But due to non-smooth ICs or discontinuities, it may become significant in our method. Therefore, it is required to locate resulting mild oscillations right at the beginning of the simulation so as to damp them as soon as possible.

Otherwise, they could get dominant in long propagation due to the cumulative nonlinear effect. In this subsection, the importance of highest-order sensor SSN is highlighted for this purpose. Figure 7 presents the same three different shock sensors as in previous subsection, but now corresponding to va , along with a zoom-in of the plot over x-axis of va around shock. Recall that the S S k v a is constant in each element as it is computed from the modal solution of that element, and this constant is assigned to each node in that particular element.

All these plots are made after the propagation over around a half shock length. As before, the subplot-(a) shows the value of (S S ) k v a (t), k = 1, • • • , K, when only the first-order sensor ( 26) is considered, whereas the subplot-(b)

shows the map (S S ) k v a (t), k = 1, • • • , K, when considering only the highest-order sensor [START_REF] Gabard | Discontinuous Galerkin methods with plane waves for time-harmonic problems[END_REF]. Lastly, the subplot-(c)

shows the map of (S S ) k v a (t), k = 1, • • • , K when the full shock sensor is used.

It is important to observe that there is no clear pattern in Figure 7-(a) for a first-order sensor. Moreover, the value of the shock sensor keeps very small compared to the previous case, and consequently, the viscosity imposed is feeble. It can therefore be concluded that mild, noisy oscillations of the velocity va are not detected by the shock sensor [START_REF] Gabard | Discontinuous Galerkin methods with plane waves for the displacement-based acoustic equation[END_REF]. On the other hand, a much better pattern is visible on Figure 7-(b), with also a large value of the Nth-order shock sensor.

As a result the imposed viscosity is stronger. However, the shock sensor is distributed over all the numerical domain.

Because of this non-distinguishable pattern of S S map, the viscosity is spread almost everywhere. This could lead to unwanted dissipation. Nevertheless, SSN is definitely required as it senses the oscillations. In the case of Figure 7-(c), a clearly distinguishable pattern is obviously visible where the oscillations are important. Also, SS value is significant enough to impose the required viscosity. As the contrast in the SS pattern is significant, introduction of viscosity is more localized near the oscillations compared to the previous case.

All the above observations are supported by the plot over x-axis of va shown in Figure 7-(d). The value of va corresponding to SS1 (blue) is, as expected, greater than other two, because there is almost no viscosity damping it. On the contrary, the dissipation is maximum in the case corresponding to SSN (green) as the shock sensor is the largest. However, in the case corresponding to the full shock sensor (red), the dissipation keeps nevertheless evident compared to the first-order sensor. Therefore, the first-order sensor is not able to sense the mild oscillations caused by non-smooth ICs or discontinuities, contrarily to the highest-order sensor and the full shock sensor. We can conclude this section with the inference that neither SS1 nor SSN are by themselves sufficient to capture weak acoustic shock waves in the numerical method based on DGM using fully unstructured meshes. However, the full shock sensor (9)-( 10) which is the amalgamation of SS1 and SSN, is able to simultaneously capture shocks as well as detect mild oscillations caused by the non-smooth part of the waveform.

Implementation of ECSAV

In this part, the influence of the implementation of the maps of artificial viscosity is analyzed for different strategies. At this point, it is important to highlight that the method is implemented to use GPU (Graphical Processing Units) and especially Nvidia graphic cards with the CUDA framework. Parallelization is done in an element-centric way where all the computations are made in parallel in each element. If the artificial viscosity is not taken into account, communication between elements arise only for one stage (exchange of the numerical flux across the boundaries) and requires only for each element to communicate with their closest neighbors. As mentioned in section 3.2, the key idea in ECSAV is to have smooth maps of viscosity. It is important to recall that DGM involves elements (triangle in this work) with quadrature points inside, which depend on the order of polynomials (see Figure 8-a for an illustration).

Formulation given in section 3.2 provides a smooth viscosity inside each element (on the quadrature points) but do not necessarily between elements.

Four different ways of implementing artificial viscosity are proposed here. Figure 8 provides an illustration of the allocation through the unstructured mesh. Let us assume the r th element is an infected element i.e., (η 0 ) r m 0. It shares its edges with three different elements, say, r e 1 , r e 2 , r e 3 , we call these neighbors edge neighbors (ENs) (in case of a boundary element, it will have two edge neighbors only). Moreover, its vertices are shared by elements other than just the edge neighbors, let us label them , r v 1 , • • • , r v R , (say), and we call them vertex neighbors (VNs).

The first approach is when a piecewise constant viscosity denoted as 'CV' (similar to [START_REF] Persson | Sub-Cell Shock Capturing for Discontinuous Galerkin Methods[END_REF]) is introduced instead of a smooth artificial viscosity, as shown in Figure 8-(b):

η r m (x) = (η 0 ) r m . (28)
The second approach is when ECSAV is introduced only in the respective element without any interaction with the neighbors, as shown in Figure 8-(c):

η r m (x r i ) = (η 0 ) r m exp        - x r i -(x 0 ) r (σ 0 ) r 2 - y r i -(y 0 ) r (σ 0 ) r 2        (29)
The third approach is when the ECSAV in the r th element interacts with the ECSAVs of its three neighbors elements. In the following it is referenced as ECSAV + EN, as shown in Figure 8-(d): ( 29) is appended by the edge contributions, given by

η r m (x r i ) = (η 0 ) r m exp - x r i -(x 0 ) r (σ 0 ) r 2 - y r i -(y 0 ) r (σ 0 ) r 2 +(η 0 ) re 1 m exp       - x r i -(x 0 ) re 1 (σ 0 ) re 1 2 - y r i -(y 0 ) re 1 (σ 0 ) re 1 2       +(η 0 ) re 2 m exp       - x r i -(x 0 ) re 2 (σ 0 ) re 2 2 - y r i -(y 0 ) re 2 (σ 0 ) re 2 2       +(η 0 ) re 3 m exp       - x r i -(x 0 ) re 3 (σ 0 ) re 3 2 - y r i -(y 0 ) re 3 (σ 0 ) re 3 2       (30)
This makes the viscosity map smoother than the two previous approaches. This is important because, as shown in 1D, discontinuities in the viscosity function could induce oscillations at the element boundaries. To further smoothen the viscosity function, it is convenient to take also into account the viscosity contributions of the vertex neighbors, as shown in Figure 8-(e) and referenced by ECSAV + EN + VN. This is achieved by appending the viscosity function in [START_REF] Hartmann | Adaptive discontinuous Galerkin methods with shock-capturing for the compressible Navier-Stokes equations[END_REF] by viscosity of vertex neighbors, given by

η r m (x r i ) = (η 0 ) r m exp - x r i -(x 0 ) r (σ 0 ) r 2 - y r i -(y 0 ) r (σ 0 ) r 2 +(η 0 ) re 1 m exp       - x r i -(x 0 ) re 1 (σ 0 ) re 1 2 - y r i -(y 0 ) re 1 (σ 0 ) re 1 2       + • • • • • • • • • • • • +(η 0 ) rv 1 m exp       - x r i -(x 0 ) rv 1 (σ 0 ) rv 1 2 - y r i -(y 0 ) rv 1 (σ 0 ) rv 1 2       + • • • • • • • • • • • • +(η 0 ) rv R m exp       - x r i -(x 0 ) rv R (σ 0 ) rv R 2 - y r i -(y 0 ) rv R (σ 0 ) rv R 2       (31)
These different strategies are tested on the propagation of a plane wave on a 2D domain with unstructured mesh.

Parameters are the same as in section 5.1.1 except the initial waveform which is chosen to be the opposite : pa ( x, ȳ, t =

0) = sin ( x) ; if | x| ≤ π.
Note that this change of sign in the initial waveform has an important consequence in nonlinear acoustic, since theoretically the solution after a shock distance is no longer a sawtooth wave, as in previous section, but an N-wave. see that the viscosity maps are smoother and smoother as expected. The locations of the high values of viscosity are also remarkable, on each subplot, they correspond to the locations of the two shocks of the N-wave. Subplot (e)

shows the waveform extracted at the center of the computational domain (ȳ = 0) for the different implementations.

Each of them allows to recover the classical N shape. Nevertheless, the zooms provided in subplots (f) and (g) show that spurious oscillations are reduced if the viscosity is smoother. It demonstrates that the strategies with the edge neighbors and with the vertex neighbors are better than the ones without interaction between elements. The additional cost of communication between elements is not very important and this strategy can easily be implemented in practice.

In the following parts, the viscosity maps are computed with the fourth option (ECSAV + EN + VN). Figure 10 shows the comparison between the numerical solution with this implementation and the 1D reference solution for the dimensionless pressure. There is an excellent agreement concerning the position of the shocks and their amplitude.

Moreover the zooms (second and third subplots) show that spurious oscillations are very weak.

To further demonstrate the robustness of the shock sensor, we consider the following initial condition made of two signals with different amplitudes:

pa ( x, ȳ, t = 0) =          R sin( x + π) if | x + π| ≤ π sin( x -π) if | x -π| ≤ π 0 otherwise , ūa ( x, ȳ, t = 0) = pa ( x, ȳ, t = 0), va ( x, ȳ, t = 0) = 0,
where R < 1 is the amplitude of the first sine-pulse, the amplitude of the second sine-pulse is assumed to be 1 (ie the reference for the computation). This initial condition is an extension of the initial condition taken in the previous test.

Three different cases are considered with R = {0.05, 0.10, 0.20}. All the other numerical parameters and the mesh are the same as in the previous test. This maintains the same ECSAV parameters as in the above case, and thus can be used to study the effect of ECSAV (triggered by a larger amplitude) on smaller amplitudes. Also, unlike the previous test, the distance of propagation is taken to be 20 shock lengths (Eq. 21) for the bigger sine-pulse (previously it was 5 shock lengths). Figure 11 presents the comparison of the three cases with their quasi-analytical solutions. The shock length is inversely proportional to the amplitude of the wave: the higher the amplitude, the shorter the shock length.

Therefore, the values of R correspond to 1, 2 and 4 shock lengths based on the smallest signal respectively. So the nonlinear effects are up to 20 times smaller for the small sine-pulse than for the high sine-pulse. Nevertheless, in all the cases, none of the smaller amplitudes are excessively dissipated. This test demonstrates the ability of the method to handle shocks with very different amplitudes. Indeed, the method does not impose extra dissipation in smooth regions of the waveform, and is capable to model propagation of signals of different amplitudes for long distances.

Results of this section illustrate the ability of the proposed method to simulate the propagation of a plane acoustic shock waves on an unstructured mesh. In the next section, more complex 2D configurations are investigated.

Propagation of Nonlinear Cylindrical Waves

The goal of this section is to compare the numerical results to a quasi-analytical solution in the case of nonlinear propagation of cylindrical waves. This test is meant to demonstrate the shock formation along with the geometrical spreading away from the center. We consider a circular domain and choose the physical parameters to be those of water: ρ 0 = 1000 kg.m -3 , c 0 = 1500 m.s -1 , β=3.5. The cylindrical waves are generated by imposing a Gaussian pulse as initial condition:

p a (x, y, t = 0) = p m a exp -

x 2 + y 2 λ 2 , ( 32
)
where the amplitude p m a = 1.5 × 10 8 Pa and λ = 1.5 × 10 -3 m. Note that the width of the pulse is approximatively λ. The computational domain has a radius of 40λ as shown in Figure 12-(a). It is meshed anisotropically, the left part of the mesh is discretized using 1 element per wavelength, whereas 4 elements per wavelength are used on the right side. The polynomial order is N=6. This mesh allows to have a very good accuracy along the positive x-axis and to limit the global number of nodes. In Figure 12-(b) the x-component of the velocity u a is shown, the global shape of the waveform remains the same in all directions, however some numerical dispersion is clearly visible on the left side of the domain, due to a poor discretization. Consequently, the distribution of ECSAV, η 2 , is significantly different in the left and right semi-circles, as seen in Figure 12-(c). This difference highlights the sensitivity of the shock sensor to small amplitude oscillations (like noise) due to dispersion in the coarser mesh, which increases the magnitude and the spread of the ECSAV.

The numerical results can be compared to a quasi-analytical solution. The nonlinear propagation of cylindrical waves can be described by the generalized inviscid Burgers equation for cylindrical waves [START_REF] Rudenko | Theoretical Foundations of Nonlinear Acoustics/English translation by RT Beyer[END_REF][START_REF] Enflo | Saturation of nonlinear spherical and cylindrical sound waves[END_REF]:

∂v ∂r + v 2r - β c 2 0 v ∂v ∂t = 0, ( 33
)
where v is the particle velocity, r is the propagation distance.

This equation is valid outside the source region for small acoustic Mach number, for more details see [START_REF] Rudenko | Theoretical Foundations of Nonlinear Acoustics/English translation by RT Beyer[END_REF][START_REF] Enflo | Saturation of nonlinear spherical and cylindrical sound waves[END_REF]. This equation can be rewritten under a dimensionless formulation:

∂v ∂ σ -µv ∂v ∂τ = 0, (34)
where

v = v v 0 r r 0 1/2
, and σ = r r 0 1/2 , with r 0 and v 0 are the radius of the source and the maximum particle velocity at r = r 0 . The dimensionless retarded time is τ = ω(tr/c 0 ), and µ = 2βk r 0 measures the importance of the nonlinear effects, smaller values imply weaker nonlinearity. This dimensionless formulation is analogous to the Burgers equation for plane wave. Therefore, the Burgers-Hayes quasi-analytical solution [START_REF] Hayes | Sonic Boom Propagation In A Stratified Atmosphere, With Computer Program[END_REF][START_REF] Coulouvrat | A quasi-analytical shock solution for general nonlinear progressive waves[END_REF] of the previous section can also be used here by changing the variables.

As mentioned above, the cylindrical Burgers equation is only valid away from the source. Therefore, the velocity field is extracted at a distance greater than λ, the approximate width of the source. We choose to extract the velocity field at r 0 = 2λ. In practice, a probe located at (x = r 0 , y = 0) is used during the numerical simulation to store the velocity field every time step. This signal is used as an initial condition for the quasi-analytical method, that is why it is referred to "initial condition" in the following paragraphs. Then, the quasi-analytical solution is compared to the numerical signals extracted with other probes located every λ on the positive x-axis. Note that the value of the initial pressure of the Gaussian pulse has been chosen in order to have µ ≈ 1 with r 0 = 2λ. It is important that the initial conditions have a high amplitude to observe nonlinear effects because the geometrical spreading for cylindrical waves scales as r 1/2 .

Figure 13 shows the comparison between the quasi-analytical solution and the numerical solution at three different distances of propagation 1λ, 17λ, 35λ, respectively. Here, the distances of propagation are along the positive x-axis, defined from the the point where the initial condition has been extracted. The waveforms are presented in the left column and their respective spectra are presented in the right one. Each subplot shows 1) the initial condition (black),

2) the linear analytical solution (green with µ = 0, models geometrical spreading only), 3) the nonlinear quasianalytical solution (red), and 4) the numerical DG solution (blue). Concerning the waveforms, it is important to steepening is very well-matched with the quasi-analytical solution in all cases. Concerning the spectrum, the three subplots show the transfer of energy to the higher harmonics with the propagation distance. The steeper the shock, the greater the energy in the higher harmonics. A close match of both nonlinear and quasi-analytical spectrum is observed, while the linear spectrum remains similar in all three cases without any transfer of energy to higher harmonics. It is important to note here that, the propagation of the acoustic shock is from right to left, this is counter-intuitive, as the wave is propagating outwards i.e., from left to right. This is due to the fact that we are working in a retarded time frame (i.e., when we move with the wave in a time window). Further, the combined effects of nonlinearity and geometrical spreading is illustrated using Figure 14. It shows the maximum amplitude of the first shock, in retarded time, for both numerical and quasi-analytical solutions. Since the comparison is made in retarded time, the shock moves from right to left. Since, the figure is presented in retarded time, if the same curve was plotted for a linear propagation, then the position would remain unchanged (no shift just a decrease of the amplitude due to the geometrical spreading). There is a good agreement between the two curves, which is a proof of the ability of the numerical method to accurately reproduce the shock speed, even after the introduction of artificial viscosity, though some noise is visible around 1.4µs.

In this configuration, geometrical spreading is a much more dominant effect than nonlinearity, this is why a very high amplitude initial pulse is chosen to be able to see the nonlinear effects. In the following test case, a configuration where nonlinear effects are stronger is investigated.

Nonlinear radiation from a piston-transducer

In order to further validate the method, the DG results are compared with the numerical results obtained with the HOWARD method [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF]. It is a one-way method (the propagation is modeled only in one direction and cannot model reflections or back-scattering) which is able to simulate the propagation of nonlinear acoustic waves in weakly heterogeneous media. One-way methods are very popular in nonlinear acoustics because they give good results in a short time. Nevertheless, they suffer from intrinsic limitation: due to the one-way approximation they cannot simulate the backscattering, if any.

The chosen validation test is based on Dagrau's test case [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF]. It deals with the radiation of a pure tone piston in a homogeneous medium. The Rayleigh distance and the shock formation distance are chosen to be of the same order of magnitude in order to provide a test case involving diffraction and nonlinearity with equal importance.

The piston is located on the left side of a two-dimensional computational domain. It radiates a pure tone at the amplitude p m a = 1360 Pa and a frequency f 0 = 1000 Hz for t = 0.03 s (30 periods). The medium is homogeneous, its physical parameters are being taken from Dagrau et al. [START_REF] Dagrau | Acoustic shock wave propagation in a heterogeneous medium: a numerical simulation beyond the parabolic approximation[END_REF]: c 0 = 340 m/s, ρ 0 = 1.2 kg/m 3 and the nonlinear parameter β = 1.2 (these values correspond to propagation in air, the acoustic Mach number is = 0.0098). The piston radius a is equal to 2λ. So, the Rayleigh distance is L R = ka 2 /2 = 4πλ. The choice of these parameters implies comparable values for the Rayleigh distance L R and the shock formation distance L sh (Eq. 21): L sh /L R ≈ 1.077.

The computational domain is rectangular but its size depends on the numerical method. Indeed, the HOWARD method is based on Fourier transforms and requires a large lateral extension while it is not necessary (and numerically costly) for the DG method. The numerical domain considered for the DG method is presented in the Figure 15-left, the piston transducer is shown in y ∈ [-2λ, 2λ]. The mesh is built using the 6 elements per wavelength along the central axis and 1 elements per wavelength elsewhere. The polynomial order of approximation is taken to be 8 throughout the mesh. The ECSAV parameters are taken to be α 1 = 10, α 2 = 20, α 3 = 1e -4. For DG, the boundary conditions are non reflecting conditions (using characteristic method [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF]) on all the boundaries except along the left boundary where the pressure is imposed:

p(x = 0, y, t) = w(y) sin(2π f 0 t),

Conclusions

The present paper introduces three new ingredients for the numerical simulation of weak acoustic shock wave propagation by means of discontinuous Galerkin method. The first one is a novel sub-cell shock capturing tool (Shock Sensor), coupling a gradient sensor and a noise sensor. Though here implemented in the framework of the DG method, we believe this tool is more general and can be made independent from the choice of the numerical method. The second innovation is the Element Centered Smooth Artificial Viscosity -ECSAV introduced to smoothen the solution in case a shock is detected by the shock sensor in one mesh element. The last element is the Gradient Factor measuring the steepening of the waveform and accordingly scaling the artificial viscosity. Special attention is paid to discuss the implementation of the shock sensor and the element centered artificial viscosity on unstructured meshes. The resulting numerical solver is equipped for simulating the propagation of weak acoustic shock waves on unstructured meshes, as modeled by a first-order system of equations in conservative form. Validation tests are performed for 1D and 2D configurations, highlighting that the method can model the nonlinear effects well along with the other physical effects like geometrical spreading and diffraction.

,

  where p m a = max |p a |; ūa =

Figure 1 :

 1 Figure 1: Comparison between smoothness indicator (SI) and shock sensor (SS). Left: the unstabilized DG solution at distance σ=1.07 computed with 40 elements and polynomial order 8 with no viscosity, compared to the quasi-analytical solution. Center: Persson's SI for the unstabilized DG solution. Right: Shock sensor SS for the unstabilized solution. The black dotted line shows the threshold above which a shock is sensed.

Figure 2 :

 2 Figure 2: Comparison of the three viscosity implementations namely, Persson's approach (SI+CV), constant viscosity based on SS (SS+max(ECSAV)), ECSAV based on SS (SS+ECSAV) with the quasi-analytical solution in a domain with K = 50 elements and order of polynomial N = 4 at distance σ = 1.455. Top left-subplot shows the pressure variation with corresponding viscosity distribution in lower subplot. The zoom-in of the lower part of discontinuity is shown in the right plot.

Figure 3 :

 3 Figure 3: Full propagation of a sine-period in a domain with K = 50 elements and order of polynomial N = 4 at different distances, the pressure(top) and its respective ECSAV(bottom) is presented. Top Left: σ = 0.182; Top Right: σ = 0.545; Bottom Left: σ = 1.091; Bottom Right: σ = 2.0.

Figure 4 :

 4 Figure 4: Top: Comparison of the complete numerical method (DG+SS+ECSAV) with the quasi-analytical solution for initial condition with three shocks. Middle: Viscosity allocation in the domain. Bottom: SS in each element, the dotted line indicates the threshold above which a shock is sensed. Left: Initial condition. Right: Simulation at around 2 shock lengths.

Figure 5 :

 5 Figure 5: Same as Figure4. Left: at around 3 shock lengths. Right: at around 5 shock lengths.

Figure 6 :

 6 Figure 6: Comparison of the different numerical solutions after the propagation of a sine-period (25) over three shock lengths. S S maps corresponding to ρ a are shown in subplots (a), (b), (c) using the definition 'SS1':(26), 'SSN': (27), 'SS1+SSN': (9)-(10), respectively. A zoom of pressure plot pa over x-axis ȳ = 0 of pa of all the three cases is shown in subplot-(d).

Figure 7 :

 7 Figure 7: Comparison of the different numerical solutions after the propagation of a sine-period (25) till half shock length. Shock Sensor maps (S S ) corresponding to va is shown in subplots (a), (b), (c) using the definition 'SS1':(26), 'SSN': (27), 'SS1+SSN': (9)-(10), respectively. Plot over x-axis of va of all the three cases is shown in subplot-(d).

Figure 8 :

 8 Figure 8: Illustration of different approaches for viscosity allocation in an unstructured mesh. Left plot (a) : example of unstructured mesh. Plots (b) to (e): examples of viscosity maps for CV (b), ECSAV (c), ECSAV+EN (d) and ECSAV+EN+VN approaches (e). See text for more details.

Figure 9 :

 9 Figure 9: Comparison of the different numerical solutions after the propagation of a sine-period around five shock lengths. Viscosity corresponding to ρa is shown in subplots (a), (b), (c), (d) using the definition 'CV':(28), 'ECSAV': (29), ECSAV+EN: (30), ECSAV+EN+VN: (31), respectively. Plot over x-axis of pa of all the four cases is shown in subplot-(e) with zoom-in near the discontinuities in subplots (f),(g).

Figure 10 :

 10 Figure 10: Comparison of the plot over x-axis of the numerical solution (ECSAV+EN+VN) after the propagation of a sine-period around five shock lengths with a quasi-analytical solution of the 1D inviscid Burgers equation. The numerical solution is obtained using the approach (31) for the introduction of viscosity. Zoom-in of the left and the right shock are shown in the center and the right subplots, respectively.

Figure 9

 9 presents results for the four different implementations of artificial viscosity. The viscosity is computed on each point of the element. Subplot (a) takes the constant state as each point of the element, subplot (b) computes a Gaussian within an element with centroid of that element as the center of the Gaussian, subplot (c) constructs the Gaussian like in (b) and also adds the overlapping three Gaussian adjacent to the edges, whereas subplot (d) does everything as (c) and also adds the overlaps from the vertices (respectively CV, ECSAV, ECSAV + EN and ECSAV + EN + VN). The results are shown for a distance of propagation of around five shock lengths. One can

Figure 11 :

 11 Figure11: Comparison of the plot over x-axis of the numerical solution (ECSAV+EN+VN) after the propagation of a two sine-periods with amplitudes: 5 MPa and 5%, 10%, 20% of the maximum amplitude with their quasi-analytical solution in subplot (left), (center), (right), respectively. The distance of propagation is taken to be 20 shock lengths of the bigger sine-pulse.

Figure 12 :

 12 Figure 12: (a) Mesh used for the simulation, (b) velocity field, and (c) ECSAV map for u a .

Figure 13 :

 13 Figure13: Temporal waveforms(Left) and their corresponding spectra(Right) at three distances from the point on the +x-axis where the initial condition is extracted : (d = {1λ, 17λ, 35λ}, from top to bottom). The black curve is the initial condition (i.e., the signal extracted at r 0 = 2λ in the DG simulation), the green curve corresponds to the solution after a linear propagation of the initial condition, the red curve corresponds to the solution after a nonlinear propagation of the initial condition, the blue curve is the numerical solution computed by the DG method.

Figure 14 :

 14 Figure 14: Shock position along the x-axis vs the amplitude, note the shock is moving leftward by the virtue of retarded time formulation. The data are extracted along the positive x-axis.

Figure 15 :

 15 Figure 15: left: Geometry and its corresponding mesh; Right: Maximum pressure for each point of the computational domain, the color map is computed by DG and the contours by Howard

Figure 15 -

 15 Figure15-right shows the maximal pressure for each point of the computational domain. The color map is computed by the DG method and the contour lines are computed using the HOWARD method. There is an excellent agreement between the two methods concerning the positions of the maxima and minima. Figure16presents the maximum pressure along the x-axis (y = 0) in both the a) linear and b) nonlinear regimes. Linear regime is simulated by forcing = 0 which deactivates nonlinearities in both methods. Therefore, we can see the comparison without nonlinear effect. The two curves are well superimposed. Nonlinear effects increase the maximum value of the pressure field after the near-field. This is due to the formation of shock waves which are interacting with diffraction.Acoustical shock waves are visible in Figure17-(a) which displays a snapshot of the pressure field. They correspond to very sharp transitions between positive and negative parts near the central region. Figures17-(b) and (c) shows the corresponding ECSAV maps for ρ a and u a . We can see that the artificial viscosity is mainly located at the position of the acoustic shocks. It is also evident that the two maps of artificial viscosity are not exactly the same.

Figure 16 :

 16 Figure 16: Maximum pressure along the axis computed by DG (line) and HOWARD (dash) methods for Dagrau's test case [20]: (Left) in linear regime, (Right) in nonlinear regime

Figure 17 :

 17 Figure 17: Snapshots of the pressure (a), and the ECSAV introduced for ρ (b) and u (c).

Acknowledgement

This work has been accomplished under the framework of an Indo-French project (No. 4601-1) funded by CE-FIPRA (Indo-French Center for the Promotion of Advance Research) and partially aided by EGIDE (Campus France).