A fast-marching like algorithm for geometrical shock dynamics

Journal of Computational Physics 284 (2015) 206–229

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A fast-marching like algorithm for geometrical shock dynamics Y. Noumir a , A. Le Guilcher b , N. Lardjane c,∗ , R. Monneau b , A. Sarrazin b a b c

LRC-MESO, CMLA, ENS Cachan, 61 avenue du Président Wilson, 94235 Cachan, France CERMICS-ENPC, Cité Descartes, 6-8 Avenue Blaise Pascal, Champs-sur-Marne, 77455 Marne-la-Vallée Cedex 2, France CEA, DAM, DIF, F-91297 Arpajon, France

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 4 November 2013 Received in revised form 8 December 2014 Accepted 10 December 2014 Available online 19 December 2014

We develop a new algorithm for the computation of the Geometrical Shock Dynamics (GSD) model. The method relies on the fast-marching paradigm and enables the discrete evaluation of the first arrival time of a shock wave and its local velocity on a Cartesian grid. The proposed algorithm is based on a first order upwind finite difference scheme and reduces to a local nonlinear system of two equations solved by an iterative procedure. Reference solutions are built for a smooth radial configuration and for the 2D Riemann problem. The link between the GSD model and p-systems is given. Numerical experiments demonstrate the efficiency of the scheme and its ability to handle singularities. © 2014 Elsevier Inc. All rights reserved.

Keywords: Geometrical shock dynamics Fast-marching method Level set method Riemann problem Finite difference method Shock wave p-system

1. Introduction In 1957, when G.B. Whitham published the Geometrical Shock Dynamics (GSD) model [29], he qualified it as “a relatively simple approximate method developed for treating problems of the diffraction and stability of shock waves”. The simplicity comes from the fact that the shock front is seen as a surface, Γ (t ) at time t, evolving under its own normal local speed, D n , and curvature, independently of the post-shock flow. The shock adjusts itself to changes in the geometry only [30]. As explained by Best [8], Whitham considered the motion of a shock into a uniform gas at rest, down a tube of slowly varying cross sectional area, A, and under some physically grounded hypothesis, he obtained an expression relating the local shock Mach number, M, to A, now known as the A–M relation [32]. Assuming that the shock front is characterized by a level set function φ in space Ω , i.e. Γ (t ) = {x ∈ Ω; φ(t , x) = 0}, one gets the Hamilton–Jacobi equation

∂t φ + D n |∇φ| = 0.

(1)

For a single pass front, i.e. when each location is reached only once, the shock position can be written as φ(t , x) = α (x) − c 0 t = 0 where c 0 is the constant sound speed of undisturbed air, and α is the shock position, also called the travel time by extension. Following Whitham [30], the GSD model reads





M (x)∇ α (x) = 1,

*



div

n(x)



A ( M (x))

Corresponding author. Tel.: +33 1 69 26 40 00. E-mail address: [email protected] (N. Lardjane).

http://dx.doi.org/10.1016/j.jcp.2014.12.019 0021-9991/© 2014 Elsevier Inc. All rights reserved.

=0

(2)

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

207

∇α where n = |∇ α | is the local normal to the front, and M = D n /c 0 > 1 is the local Mach number. This model is hyperbolic provided that A  ( M ) < 0, and can thus develop disturbances on the front which are the traces of waves, not modeled, behind the shock. For GSD, the relation between A and M is

A (M ) A(M ) with

=−

M λ( M ) M2 − 1

⎧    2 1 − μ2 1 ⎪ ⎪ 0 < λ( M ) = 1 + 1 + 2μ + 2 , ⎪ ⎪ ⎨ γ +1 μ M ⎪ ⎪ (γ − 1) M 2 + 2 ⎪ ⎪1≥μ= , ⎩ 2γ M 2 − (γ − 1)

(3)

(4)

γ > 1 being the gas polytropic coefficient, μ the local post-shock Mach number and by definition M > 1. Thus, the computations only make sense when M ≥ 1, which is compatible with the eikonal equation. As stated in [30, pp. 296–297], this reservation has to come into consideration when some initial/boundary conditions lead to the apparition of Mach numbers lower than 1 and have to be excluded. The main drawbacks of the model are in the weak shock limit, where the speed of perturbations along the front are underestimated by a factor of two [28] and the absence of regular reflection for shock/solid wall interaction. Nevertheless, in practice, GSD has proven to be fairly accurate for diffraction around a corner, non-regular Mach reflection [30], or accelerating shocks and shown only little deviation for expanding decelerating flows [5]. In the last three decades, Whitham’s model has been extended to take into account unsteady flow behind the shock [8–10], non-uniform gases properties [34], and has been applied, among others, to imploding shock waves [11,1], atmospheric propagation [7], detonation in explosives [2,3,6], supersonic engine unstart [28] and astrophysics [14]. This success, linked to the compact model formulation and the dimensional reduction, was supported by the development of three kinds of algorithms. (i) Lagrangian, or front-tracking, methods have first been experimented [15,18]. In such an approach, the shock front is explicitly discretized by markers evolved in time and regularly resampled. This method is natural and quite accurate but difficult to implement in three dimensions, mainly when surface merging or breaking is expected. (ii) Eulerian conservative algorithms [19,20] reduce this difficulty but do not take any advantage of the front locality. Furthermore, they rely on the definition of an a priori propagation direction, not always easy to determine. (iii) Localized level set methods are a good compromise since they handle any kind of surface deformation but in an implicit way. The front shock is obtained from a table of arrival time, also called burn table. A 3D unsteady algorithm, based on the Hamilton– Jacobi form of the GSD system (2) [17], is described in [2–4] for Detonation Shock Dynamics. It compares well with reactive Eulerian model results at a much lower CPU time. Nevertheless, due to the nonlinear nature of GSD equations, nonphysical shocks can form away from the front position and a frequent resampling of the signed distance is mandatory [23]. In this article we propose an alternative approach based on the level set fast-marching paradigm [22], which combines the flexibility of (iii) and the locality of (i) while remaining easy to implement. The algorithm is described in Section 2. Reference solutions for a smooth radial problem and the GSD Riemann problem are derived in Section 3. In Section 4, the comparison to numerical results and the performance of the algorithm are discussed. At last, conclusions are summarized in Section 5. 2. A fast-marching like GSD scheme In 1988 Osher and Sethian [17] introduced the Eulerian level set method to solve Hamilton–Jacobi equations and the eikonal equation in particular. Unlike the Lagrangian approach, the level set method is simple to implement in 3D, highorder extensions are readily derived and topological properties of the front, as the curvature, are easily calculated. However, assuming a cubic domain in a space of dimension d discretized with N grid points per direction, the unsteady level set method has a complexity O ( N d+1 ), which can be quite time consuming when N is large. In the past two decades, several improvements have been proposed in order to reduce this complexity and at same time to enhance the accuracy. As explained in [16], one of the most popular approach is the Fast-Marching Method (FMM) first proposed by Tsitsiklis in [27], rediscovered later by the level set community, and popularized by Sethian [22]. It has since been used with success in a variety of applications. For a single pass front, the complexity of the FMM reduces to O ( N d log N ) and even to O ( N d ) under some further assumptions [31]. In this section, we introduce a method to solve the GSD model (2), on a Cartesian grid, based on the fast-marching paradigm. We first reformulate the model as a coupled eikonal-transport system to facilitate its discretization. The numerical method, boundary treatment and implementation details are then given. 2.1. A modified transport equation As in the work of Besset and Blanc [7] or Aslam [2], the GSD model (2) is rewritten under the local form

208

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

⎧ (a) ⎨ M |∇ α | = 1 ∇α ˙ ( M , κ ) (b) ⎩M · ∇M = M |∇ α |

(5)

with the initial boundary conditions



α|Γ0 = α0 ,

M |Γ = M 0 , 0

(a)

(6)

(b)

where α0 and M 0 are given functions on the hypersurface Γ0 , the initial shock position. Using Eqs. (2)–(4) linking the local Mach number, M, and the mean curvature of the front, κ = div(n), the GSD closure reads

˙ (M , κ ) = − M

M2 − 1

λ( M )

κ.

(7)

Since λ( M ) > 0 for M > 1, Eqs. (7) and (5b) show that M decreases for expanding fronts. In the weak shock limit, assuming 1−d

a radial expansion, one easily show that M (r ) ∼ 1 + ( M a − 1)( rr ) 2 with r ≥ ra , M a > 1 and d = 1, 2, 3 for plane, cylindrical, a or spherical case respectively. This equation indicates that the Mach number remains greater than unity, at least for smooth problems, which is a desirable feature for fast-marching methods where a single pass front is assumed. In practice, the discretization of the mean curvature is difficult due to the mixed derivatives of α , especially in the context of the fast-marching method when neighboring points are not yet assigned a value. For this reason we choose to reformulate the transport equation on the Mach number as a convection–diffusion one. By combining the eikonal equation, M |∇ α | = 1, and the normal definition, n = ∇ α /|∇ α |, one checks that the mean curvature of the front reads κ = ∇ M · ∇ α + M α . Transport equation (5b) is then reformulated as

∇ M · ∇ α = S ( M ) α , ( M −1) where S ( M ) = − M 2M is non-positive and smooth for M > 1. The reformulated boundary value problem to be solved (λ( M )+1)−1 2

is then



M |∇ α | = 1 (a) ∇ M · ∇ α = S ( M ) α (b)

(8)

where the source term is now easier to discretize. This recasting can be seen as an expansion of Eq. (8.51) given in Whitham’s book [30, p. 280]. Note that this system is not in a conservative form, which could raise difficulties to handle front discontinuities, but as we shall see in Section 4.2, the numerical scheme performs reasonably well in practice. We emphasize that the system (8a)–(8b) is fully coupled and in the discrete setting, both equations will be solved simultaneously. Eq. (8b) shows that the order of the discrete Laplacian operator must be consistent with the order of the discrete gradients on α and M. To limit the size of the Narrow-Band (defined in the next section), we restrict to first order gradients, second order Laplacian, in the sequel of this paper. Before giving the details of our algorithm, we recall the standard fast-marching method for the eikonal equation. 2.2. Standard fast-marching method for the eikonal equation The level set function α is a solution of the eikonal equation (8a) on Ω together with the initial condition (6a) on Γ0 . We seek for a solution on a uniform Cartesian mesh of the domain Ω ≡ [0, Lx] × [0, L y ] × [0, Lz] ⊂ R3 with grid spacings x, y and z. Let αi , j ,k and M i , j ,k be the approximate solution at a grid point, i.e. αi , j ,k = α (xi , y j , zk ) and M i , j ,k = M (xi , y j , zk ) where xi = i x, y j = j y and zk = k z. We denote by ul , v l and w l (resp. u r , v r and w r ) the backward (resp. forward) approximation of the derivatives of α at (xi , y j , zk ) along x, y and z respectively. The discrete form of (8a) reads

ˆ ( u l , u r , v l , v r , w l , w r ) = 1, M i , j ,k H

(9)

ˆ is chosen as the numerical Hamiltonian of Godunov [25]: where H

ˆ (u l , u r , v l , v r , w l , w r ) = H









max2 ul+ , u r− + max2 v l+ , v r− + max2 w l+ , w r− ,

with the notations: x+ = max(x, 0) and x− = max(−x, 0) for x ∈ R. This upwind scheme has been successfully used in several applications (cf. [16,22]). It has the ability to capture viscosity solutions of the eikonal equation [12] and do not smear out sharp discontinuities excessively [25], which are desirable features in the case of GSD where front singularities may appear.

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

209

Fig. 1. Sketch of points used in the classical fast-marching method. Known points: dark circles, NarrowBand points: grey squares, Far points: light circles. The newly accepted point and its local updated neighbors are in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

α are written: 

αi+1, j,k − αi, j,k := ur = D + x α i , j ,k := x x

+  αi, j,k − αi, j−1,k αi, j+1,k − αi, j,k := v r = D y α i , j ,k := y y

+  αi, j,k − αi, j,k−1 αi, j,k+1 − αi, j,k := . w r = D z α i , j ,k := z z

At first order, the discrete derivatives of







ul = D − x α i , j ,k vl = D − y α i , j ,k wl = D − z α i , j ,k

αi, j,k − αi−1, j,k

Following Sethian [22], the classical FMM is now outlined. In this method, the velocity of propagation is assumed to be known and of constant sign. The front is thus “single pass” and one needs to calculate the value of α only in the vicinity of it. The CPU time is then dramatically reduced and, interrelated, the monotonicity of the scheme is guaranteed by following the direction in which the information flows. This is done by propagating the solution from lower to higher values of the level set function α . To this end, the grid points are partitioned into three groups, see Fig. 1: ❶ Known is the set of vertexes where the values of α are known, i.e. the vertexes already intercepted by the front; ❷ NarrowBand is the set of all neighboring vertexes of Known, i.e. the vertexes that are about to be intercepted by the front; ❸ Far is the set of vertexes that are neither in Known nor in NarrowBand. In the Far set, α is assigned a huge initial value INF, greater than the maximal final expected value of α on the domain. In practice we choose INF = 1010 . Think of the NarrowBand set as a buffer zone used to start the calculation from the vertexes of Known and to spread the information to the vertexes of the Far set. For every item in NarrowBand, test values of α are calculated from the points of the Known set only using upsided discrete gradients. The point of NarrowBand corresponding to the minimal test value cannot evolve anymore and is then validated. The second step consists in including this point in the Known set and adding its close Far neighbors as new points of the NarrowBand. By repeating the previous steps, as long as the NarrowBand is non-empty, the level set function is calculated in the whole computational domain. The starting values of this process are set by the initial conditions on Γ0 . It is worth mentioning that performance of the FMM depends on the search for the point with the minimal test value. In practice, a min-heap data structure is used for this purpose. At a point (i , j , k) of NarrowBand, one can deduce from (9) that the test value ϑ of α is calculated through the following quadratic equation: 3 

max

2

 ϑ − tν−

ν =1

i , j ,k

−

ν

+

,

ϑ − t iν, j ,k +

ν

 ,0 =

1 M i2, j ,k

,

(10)

with the notations

± 1 = x, 1±

± 2 = y,

t i , j ,k = αi ±1, j ,k ,

2±

± 3 = z,

t i , j ,k = αi , j ±1,k ,

±

t i3, j ,k = αi , j ,k±1 .

When solving Eq. (10), one must consider only the largest root and deal special cases with care as explained in [16] or [33]. 2.3. Fast marching method applied to the GSD system As outlined previously, the GSD system models a one pass front with positive local velocity. These features motivated us to extend the fast-marching paradigm to the system of Eqs. (8).

210

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

Since M is an unknown of the problem, one has to calculate test values of α (i.e. ϑ ) and of M (i.e. m), at the same time, for each point of the NarrowBand, leading to a local nonlinear system in ϑ and m. The discrete normalized velocity M i , j ,k is now replaced by m in Eq. (10) and the discretization of the transport equation (8b) is done as follows. The advection part of Eq. (8b) reads ∇ M · ∇ α . Following the upwind direction of discretization for ∇ α ensures that only receivable points are used in the computation of ∇ M. More precisely, we have: 3   ϑ − t iν, j ,k + m − νi , j ,k −

(∇ M · ∇ α )|i , j ,k =

−

− ν

ν =1

− ν

−

 t ν + − ϑ − ν + − m  i , j ,k i , j ,k + ν

+ ν

±

where the coefficients ( νi , j ,k )1≤ν ≤3 are defined as: ±

±

1i , j ,k = M i ±1, j ,k ,

2i , j ,k = M i , j ±1,k ,

±

3i , j ,k = M i , j ,k±1 .

Concerning the source term of Eq. (8b), S ( M ) α , we evaluate S implicitly and discretize the Laplacian operator as follows. From our experience, it is crucial to use a centered scheme as much as possible, otherwise no disturbance will develop along the shock curve. To illustrate this point, think about a planar shock reflecting, at some angle, onto a rigid wall aligned with the mesh. This Riemann problem, detailed in Section 3.2, leads to a new intermediate state. If one uses only one sided Laplacians relying on the causality condition, the curvature would remain zero and no new wave can appear. For this reason, we take into account all finite value neighbors of a point of interest, including those of the NarrowBand, in the calculation of the test value. In one space direction, if a left and right values are available, the second derivative is chosen centered, no matter the causality condition in the NarrowBand (i.e. we also use points which are not in the Known set). For this reason, the algorithm is not a fast-marching method 3in the strictest sense, and we say it has fast-marching like properties. ˆ ν , where we express only the first term for the sake of simplicity More precisely, (S ( M ) α )|i , j ,k = S (m) ν =1

if (αi −1, j ,k < INF and

ˆ1=

ˆ 1c

αi+1, j,k < INF) then

else − + if (ϑ ≥ max(t i1, j ,k , t i1, j ,k )) then +

if (max(

ϑ−t i1, j ,k +

1

, 0) ≥ max(

ˆ1=

ˆ 1r

−

ϑ−t i1, j ,k −

1

, 0)) then

else

ˆ1=

ˆ 1

end if else ˆ1=0

end if end if with αi−1, j,k ˆ 1c = αi+1, j,k −2ϑ+ •

x2 ϑ−2αi +1, j ,k +αi +2, j ,k ˆ 1r = s+x s+x •

i , j ,k i +1, j ,k x2 ϑ−2αi −1, j ,k +αi −2, j ,k ˆ 1 = s − x s − x •

i , j ,k i −1, j ,k x2

For the upwind cases, switches are introduced to limit the influence of newly calculated NarrowBand points. Assuming that plane shocks are preserved lead us to set, in the x direction for example, x s± i , j ,k

 =

1 if αi ±1, j ,k < αi , j ,k 0 otherwise.

Note also that upwinding is at first order to enforce stability, which could reduce the accuracy on M. This is not a serious drawback in practice since upwinding principally affects new points entering the NarrowBand. To increase the use of the central formula, we allow the calculation of diagonal points in the update procedure. When a node is accepted, we thus add all its neighbors, not already in the Known set, into the NarrowBand and flag them as active. The small set of active points is then sorted by increasing values of α and updated. The extra diagonal points, called extended NarrowBand points in Fig. 2, are evaluated with an upwind Laplacian until the front moves close to them where they become “regular” NarrowBand points.

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

211

Fig. 2. Sketch of points used in the GSD fast-marching method. Known points: dark circles, NarrowBand points: grey squares, NarrowBand extended points: dark grey triangles, Far points: light circles. The newly accepted point and its local updated neighbors are in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The final form of the local nonlinear system on ϑ and m is then ⎧ 3  ϑ − tν− ϑ − tν+  

⎪ ⎪ ⎪ ⎪ ⎪ ⎨

max2

ν =1

i , j ,k

− ν

,

i , j ,k

+ ν

,0 =

1

m2

(a)

− − + + ⎪ 3  3 ⎪   ϑ − t iν, j ,k + m − νi , j ,k  t iν, j ,k − ϑ − νi , j ,k − m  ⎪ ⎪ ˆ ν . (b) ⎪ = S (m)

− ⎩ − − + + ν ν ν ν ν =1 ν =1

(11)

Two strategies have been tested for the numerical resolution of (11). Given a current approximation (ϑ ( p ) , m( p ) ), the first approach is a fixed point method which consists in solving the quadratic equation (11a) on ϑ with m = m( p ) to get ϑ ( p +1) which is injected in Eq. (11b). This gives us a new value m( p +1) by resolution with either a new fixed-point iteration or Newton’s method. This procedure is repeated until a desired tolerance is met (10−6 in practice). The second approach is to solve the full coupled system of Eqs. (11) by Newton’s procedure. Indeed, the system (11) can be written in the generic form G (W ) = 0, where G is a nonlinear function in W = (m, ϑ)t which depends on the values of α and M around the point (i , j , k) of the NarrowBand. The resolution of (11) is then done by the following iterative algorithm, for p ≥ 0,

−1 ( p ) 

, W ( p +1 ) = W ( p ) − D W G W ( p ) G W

when the computation of the gradient DW G is allowed, and where the starting point W (0) is determined from the values of M and α on neighboring points. Like for the classical fast-marching method, only the largest root of the discrete eikonal equation is retained. In practice we did not observe any robustness problem and Newton’s procedure is used. One can remark that this numerical scheme is not monotone. For this reason, it is difficult to obtain theoretical results regarding its behavior. In particular, there is no guarantee that the arrival time of an accepted point couldn’t be modified as the computation progresses. Nevertheless, for smooth expanding waves at least, the curvature should not vary significantly from point to point and the NarrowBand is then a good approximation of the front. If M is known at first order, variations on α are then of second order only. Numerical results of Section 4 confirm this and the lack of theoretical monotonicity appears not to be a problem. 2.4. Boundary conditions Two kinds of boundary conditions are commonly required: outgoing and rigid walls. In practice, fictitious points are added outside of the computational domain, on which α and the Mach number are set. The number of fictitious cells depends on the interior numerical scheme and is set to two here to allow for potential upwind Laplacian. The outgoing conditions are used on the artificial limits of a free boundary. We impose the huge INF value on α at fictitious points to make the scheme upwind. When a boundary of the computational domain is bounded by a rigid body Γ R , a wall condition is used. In this case, a Neumann condition applies: (n · ∇ α )|Γ R = 0. The numerical implementation is done by reflecting the values of α and M, by symmetry, on the fictitious points. Let us recall here that the GSD model is able to predict only the irregular shock reflection, namely the Mach reflection or shock–shock singularity in the terminology of Whitham. 2.5. Summary of the algorithm The numerical scheme we designed to solve the GSD model (2) is now summarized. Whitham’s model is first rewritten to avoid the direct discretization of the curvature, which leads to the local nonlinear system of Eqs. (8). The proposed

212

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

algorithm is Eulerian, similar to the fast-marching method performed on a Cartesian grid. It relies on a causal first order finite difference upwind scheme for gradients and a second order discretization of the Laplacian operator. The arrival time and Mach number are obtained by solving the local nonlinear system (11) via Newton’s method. The overall structure of our algorithm is close to the FMM one but has three main differences: update of diagonal points, utilization of values not yet accepted and resolution of the eikonal equation together with the transport equation on the propagation velocity. The fast-marching like algorithm contains two main steps, namely initialization and iteration, organized as follows: ➥ Initialization step ➠ ➠ ➠ ➠

Set all points of the computational domain in the Far set by assigning a huge INF value to α . Add each point of the user initial condition in the Known set. Create the NarrowBand from first neighbors of the Known set. For each point of NarrowBand, compute the test values ϑ and m by solving the local nonlinear system (11) by Newton’s method for instance and apply boundary conditions.

➥ Iterative step until the NarrowBand is empty ➠ Pick P min ≡ (i min , j min , kmin ) ∈ NarrowBand such that

α (i min , j min , kmin ) =

min

(i , j ,k)∈NarrowBand

α (i , j , k).

Add P min in the Known set/Delete it in NarrowBand.

➠ Add the close neighbors of P min , i.e. points (i min ± 1, j min ± 1, kmin ± 1), in NarrowBand if they were in the Far set and flag them as Active.

➠ Sort the Active points by increasing value of α . ➠ For each point of Active, compute the test values ϑ and m and apply boundary conditions. ➠ Unflag the Active points, delete the Active list.

3. Reference solutions to GSD In this section, we provide reference solutions to the GSD equations, that is geometrical configurations where the exact solution can be either calculated analytically or approximated with a quick by-calculation (typically the resolution of a one-dimensional ODE). The numerical solutions of our fast-marching algorithm can then be compared to these reference solutions to evaluate its precision. Ideally, they should cover a wide array of cases, to test the algorithm in the most diverse configurations. Conversely, such reference solutions are often hard to provide, and we only expect to determine them in very simple geometrical configurations. In the case of the GSD equations, it can be noted that infinitely large plane fronts of any normalized velocity M > 1 are exact solutions of the system in the whole space. Reference solutions can also be found in the radial case, when the configuration is invariant under rotations, and more interestingly, for the Riemann problem in R2 . Such solutions enable to test the behavior of the algorithm on smooth solutions as well as in cases where the velocity has discontinuities. The simplest test case for discontinuous solutions is the compression wedge (see Subsections 3.2.5 and 4.3). In the remainder of this section we will use the conservative form (2) of the GSD system, with the integral A–M relation

 M  mλ(m) dm , A ( M ) = A 0 exp − m2 − 1

(12)

M0

where M 0 > 1 is the initial condition and A 0 the corresponding section area. We also introduce the notation

β( M ) =

M2 − 1

λ( M )

.

The quantity β allows to write the results of this section in a simpler manner. 3.1. Radial solutions The simplest configuration in which a solution of the GSD equations can be computed is the radial case where the solution depends only on the radial coordinate r. More precisely, in dimension d = 2, 3, the normalized velocity M is solution of a one-dimensional ODE in the r variable:

⎧ ⎨ M (r )∂r α (r ) = 1

(a)

d − 1 ⎩ M (r )∂r M (r ) = −β M (r ) . (b)

(13)

r

For a given initial condition, this equation can be solved by a high order algorithm to provide a solution of reference, helpful for convergence studies. Analytical solutions are easily obtained in the strong shock limit, M  1.

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

213

Fig. 3. Sketch of the initial condition of the GSD Riemann problem.

3.2. Solutions of the Riemann problem in dimension 2 The solutions of the Riemann problem for the GSD equations are of prime importance since they can exhibit discontinuous velocities and develop new intermediate states. The behavior of the algorithm can then be evaluated on this difficult problem. The simplest test case for discontinuities is the compression wedge, which will be developed in more details later. For general hyperbolic systems, the theory for the resolution of the Riemann problem only covers the case where the two initial states are not far from each other. Solutions of the Riemann problem for far away initial states are known only for very particular cases, such as the p-system. While the GSD equations are in a setting that differs from the usual p-system, they can be linked to the p-system and indeed the full resolution of the Riemann problem is available. In his works [29] and [30], Whitham described the form of simple shocks and simple rarefactions. They are studied respectively as the results of the compression of a plane front by a wedge and the diffraction of a plane front by a corner. While these works fully describe the structure of shocks and rarefactions, the rarefaction solution is only given with implicit coordinates. In this part, we describe more explicitly the geometry of the solution. From our knowledge, this is the first time that a comprehensive solution of the GSD Riemann problem is given, together with a link to the p-system. Nevertheless, one can mention the works of Henshaw, Smyth and Schwendeman [15] and Schwendeman [20] where the Whitham GSD equations are rewritten in conservative form. Schwendeman [20] also mentions the fact that the Riemann problem admits simple solutions, but does not give full details. In this section the GSD equations are recast as a system of M and θ where θ parametrizes the front normal as n = (cos θ)ex + (sin θ)e y , with ex and e y the unit vectors along the x and y axes respectively. For the Riemann problem, the initial condition is made of two planar fronts as sketched in Fig. 3: a shock of velocity M and angle θ interacts with another one of velocity M r and angle θr . We work in polar coordinates (ρ , χ ), taking χ = 0 for the horizontal axis ( O x). Note that the initial position of the first front is then the half-line {χ = χ = θ − π2 }. Similarly, the initial position of the second front is the half-line {χ = χr = θr + π2 }. With these notations, the left front comes before the right one when the angle χ increases. This is coherent with the standard notations for the theory of the Riemann problem and explains why the left front is on the bottom of the picture. 3.2.1. The p-system The full resolution of the Riemann problem is known for the p-system (see [24,21]), which can be written in the form



∂t v − ∂x u = 0 ∂t u + ∂x p ( v ) = 0,

(14)

with t > 0, x ∈ R, along with the properties p  < 0 (and p  > 0, but the sign of p  plays no role in the resolution of the Riemann problem and is only a result of the usual physical interpretation in terms of isentropic gas dynamics). As one can note, the p-system is one-dimensional in space, and we want to solve the stationary GSD equations in a two-dimensional space. Nevertheless, it can be argued that the level set function α , is a disguised time variable. We then ∇α introduce the (σ , τ ) system of coordinates where eσ = n = |∇ α | is the unit vector normal to the front and eτ the unit vector tangent to the front such that (eσ , eτ ) is direct. One can note that Whitham [30] already introduced a similar system of coordinates, α and β , proportional to our σ and τ , but with a non-normalized gradient, leading to a slightly different system. Using the fact that (2) implies

  ⎧ n ⎪ ⎪ = 0, curl ⎨ M   ⎪ n ⎪ ⎩ div = 0, A(M )

(15)

214

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

one can rewrite the GSD system in our system of coordinates as



∂σ θ + ∂τ log ( M ) = 0 ∂τ θ − ∂σ log A ( M ) = 0,

which bears a strong similarity with the p-system. Keeping in mind that the setting is different, to write it as a p-system we can take



u=θ  v = log A ( M ) .

Given that A is a smooth function in M with A  ( M ) < 0, it is invertible. Noting its inverse by A −1 , the expression for the function p is



 .

p ( v ) = log A −1 e v We can then check that

p (v ) =

ev A  ( A −1 (e v )) A −1 (e v )

< 0,

because A  < 0, which guarantees that the results obtained in the case of the p-system still hold. 3.2.2. Geometrical structure of simple rarefactions A rarefaction is a smooth transition between two constant states ( M , θ ) and ( M r , θr ) (see for instance Fig. 4). As a consequence of our analysis, we will show that admissible solutions have to satisfy ∂χ θ ≥ 0 and θr ≥ θ , even if we do not assume it at the beginning of our analysis. Writing (15) in polar coordinates (ρ , χ ), we get

    ⎧ ∂ cos(θ − χ ) ∂ ρ sin(θ − χ ) ⎪ ⎪ − =0 ⎨ ∂ρ M ∂χ M     ⎪ ⎪ ⎩ ∂ ρ cos(θ − χ ) + ∂ sin(θ − χ ) = 0. ∂ρ A(M ) ∂χ A(M )

Keeping in mind that we are looking for solutions that only depend on the angular variable be rewritten as

(16)

χ , we see that Eqs. (16) can

⎧ (∂χ M ) ⎪ ⎪ cos(θ − χ ) ⎨ (∂χ θ) sin(θ − χ ) = − M

A (M ) M ⎪ ⎪ ⎩ (∂χ θ) cos(θ − χ ) = (∂χ M ) sin(θ − χ ) = −(∂χ M ) sin(θ − χ ). A(M ) β( M )

(17)

Combining the equations of system (17), we get the following fundamental relation

  β( M ) (∂χ θ) M sin2 (θ − χ ) − cos2 (θ − χ ) = 0. M

From this equation we deduce that either θ is constant (and then M is also constant), which corresponds to the case where the front is planar, or M and θ satisfy the relation

√

tan(θ − χ ) = σ

σ

β( M ) M

with σ = ±1.

(18)

This relation is satisfied for any rarefaction, as soon as θ is not constant. In particular, we have π2 > σ (θ − χ ) > 0 with = 1 for 1-rarefactions, and with σ = −1 for 2-rarefactions. The combination of the equations of system (17) also leads to the relation

(∂χ θ)2 =

(∂χ M )2 , β( M )

(19)

which gives

(∂χ M ) ∂χ θ = ε √ β( M )

with ε = ±1.

Plugging Eqs. (20) and (18) in system (17) leads to

σ ω( M ) + θ = const

(20)

σ ε = −1. We then have (21)

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

215

Fig. 4. Geometry of simple rarefaction waves.

with

M

ω( M ) =

M

− A  (m) dm = A (m)m

M0

√

1

β(m)

dm.

(22)

M0

It can then be shown with other lengthy calculations involving derivations of (18) and (21) with respect to satisfies the ODE

dM (χ ) dχ

χ , that M

 = −σ F M (χ )

(23)

along the rarefaction, with F ( M ) =

√

2 β( M )( M 2 +β( M )) 2M 2 + M β  ( M )

> 0 (when it is well defined). It can be checked that the denominator

of F ( M ) is positive when A  ( M ) > 0 and this last condition can be checked at least numerically on some values of general for large M. In particular for M (χ ) and θ(χ ) solutions of (23) and (18), the quantity

γ , or in



σ ω M (χ ) + θ(χ ) = const

(24)

is independent of χ (it is indeed a Riemann invariant here). This shows in particular that in any case we have ∂χ θ ≥ 0 and then θr ≥ θl . We also deduce from (24) that

ω( M ) − ω( M r ) = σ (θr − θl ).

(25)

We can then give the geometrical structure of 1-rarefactions and 2-rarefactions in detail. Proposition 3.1 (Structure of simple rarefactions). Assume that 2M + β  ( M ) > 0. Let χr = θr + π2 .

π > θr − θl > 0. We set χ = θ −

(i) 1-rarefaction Let M > M r > 1 be satisfying (25) for σ = 1. Then there exist two angles χ1r − and χ1r + < χr − π2 and such that: ❶ In the sector {χ ≤ χ ≤ χ1r − }, we have M = M and θ = θ . ❷ In the sector {χ1r − ≤ χ ≤ χ1r + }, M satisfies (23) and θ is given by (18) for σ = 1. ❸ In the sector {χ1r + ≤ χ ≤ χr }, we have M = M r and θ = θr . The angles χ1r − and χ1r + are given by (18) for σ = 1, so

tan θ − χ

r− 1



√

=

β( M ) M

and

tan θr − χ

Moreover (24) holds true for σ = 1 for all χ ∈ [χ , χr ].

r+ 1



√

=

β( M r ) Mr

.

π and 2

χ1r + , cf. Fig. 4(a), satisfying χ < χ1r − <

216

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

Fig. 5. Geometry of simple shock waves.

(ii) 2-rarefaction Let M r > M > 1 be satisfying (25) for σ = −1. Then there exist two angles χ2r − and χ2r − < χ2r + < χr and such that: ❶ In the sector {χ ≤ χ ≤ χ2r − }, we have M = M and θ = θ . ❷ In the sector {χ2r − ≤ χ ≤ χ2r + }, M satisfies (23) and θ is given by (18) for σ = −1. ❸ In the sector {χ2r + ≤ χ ≤ χr }, we have M = M r and θ = θr . The angles χ2r − and χ2r + are given by (18) for σ = −1, so



tan θ − χ2r − = −

√

β( M ) M

and



tan θr − χ2r + = −

√

β( M r ) Mr

χ2r + , cf. Fig. 4(b), satisfying χ +

π < 2

.

Moreover (24) holds true for σ = −1 for all χ ∈ [χ , χr ]. Because M satisfies the ODE (23) for σ = 1 along a 1-rarefaction and for σ = −1 along a 2-rarefaction, independently of θ , the shape of the front along a rarefaction is universal. More precisely if we fix a level set value α = t 0 , associated with a time t 0 > 0, there exists a spiral St0 of polar equation ρ = t 0 ρ1 (χ ) such that the level set {α = t 0 } on a rarefaction fan is the image of a portion of St0 by an isometry (a direct isometry for 1-rarefaction, and an indirect isometry for a 2-rarefaction). Moreover, when the time goes on, the shape of the spiral is simply changed by dilation. Finally this spiral is locally a convex curve, because we have shown that ∂χ θ ≥ 0. 3.2.3. Geometrical structure of simple shocks A shock is a discontinuity between two constant states ( M , θ ) and ( M r , θr ). Physically, we expect to have θr < θ for a simple shock and indeed, while it is possible to construct functions with a simple shock that are formal solutions of (15) with θr < θ , such solutions do not satisfy the entropy condition and are thus not physical (see Remark 3.2). Both quantities M and θ are discontinuous through the shock, but the level set α is continuous, albeit not smooth (see Fig. 5). The trajectory of the (punctual) shock is the half-line χ = χ s with χ < χ s < χr . In order to satisfy the divergence equation in (2), we need to have χ s − χ < π and χr − χ s < π and to have either χ < χ s < χr − π2 or χ + π2 < χ s < χr . Otherwise, the flux of nA would be entering the shock from both sides, which would contradict the fact that nA is divergence free. When χ < χ s < χr − π2 , we have a 1-shock and the continuity of α forces M r > M . Using system (15), we can get the relations between M , M r , θ , θr and χ s = χ1s with θ , θr > χ1s . The first equation and the continuity of α lead to

M cos(θ − χ s )

=

Mr cos(θr − χ s )

(26)

.

The second equation gives that the flux of A (nM ) through the boundary of a closed domain is zero. If we consider the flux for a tube going through the interface between the two states, we get

A(M ) sin(θ − χ s )

=

A(Mr ) sin(θr − χ s )

(27)

.

From (26) and (27) we deduce



tan θr − χ s =

A ( M r ) sin(θ − θr ) A ( M ) − A ( M r ) cos(θ − θr )

=−

M − M r cos(θ − θr ) M r sin(θ − θr )

,

(28)

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

217

Fig. 6. Example of a phase diagram for ( M , θ ) = (10, 0) (curves R 1( M ,θ ) , R 2( M ,θ ) , S (1M ,θ ) and S (2M ,θ ) ), for ( M , θ ) = (1.5, 0) (curve R 1(1.5,0) ) and for







( M , θ ) = (1, θ0 ) (curve R 2(1,θ ) ). 0

and then by symmetry r / , we get



tan θ − χ s = −

A ( M ) sin(θ − θr ) A ( M r ) − A ( M ) cos(θ − θr )

=

M r − M cos(θ − θr ) M sin(θ − θr )

.

(29)

Notice that this implies

cos(θ − θr ) =

A(M )M + A(Mr )Mr A(Mr )M + A(M )Mr

.

(30)

We notice in particular that the general reordering inequality holds: a+ b+ + a− b− ≥ a+ b− + a− b+ for all a+ ≥ a− ≥ 0 and b+ ≥ b− ≥ 0. Therefore, the fact that A ( M ) is decreasing in M implies in particular that

0<

A(M )M + A(Mr )Mr A(Mr )M + A(M )Mr

≤1

which makes sense for equality (30). When χ + π2 < χ s < χr , we have a 2-shock and the continuity of α forces M r < M . Moreover, we have θ , θr < χ s = χ2s . The same reasoning leads to (26) and (27) which gives the same relations as the ones above. Remark 3.2. Recall that admissible shocks have to satisfy an entropy condition in order to be stable: this is the Liu Econdition (which implies the Lax E-condition), see Paragraph 8.4 in [13]. It is also known that Lax E-condition selects only half of the shock curve (see p. 244, Paragraph 8.3 in [13]). We can see it on Fig. 6. This shows that only the case θr < θ is selected for shocks. With these information in mind, we can give the structure of simple shocks: Proposition 3.3 (Structure of simple shocks). Let 0 < θ − θr < π2 and M r , M > 1 be satisfying (30). (i) 1-shock If M r > M > 1, then there exists χ s = χ1s < θ , θr given by (28) (or equivalently by (29)), such that the situation is pictured on Fig. 5(a). (ii) 2-shock If M > M r > 1, then there exists χ s = χ2s > θ , θr given by (28) (or equivalently by (29)), such that the situation is pictured on Fig. 5(b). Remark 3.4. We note that by replacing cos(θ − θr ) by its expression (30) in (29), one checks that (29) is equivalent to the relation given in [30] for a compression wedge.





tanθ − χ s  =

A M



M r2 − M 2 A 2 − A r2

1/2

with A b = A ( M b ) for b = , r .

(31)

218

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

Fig. 7. Geometrical determination of the intermediate state for a double rarefaction.

Remark 3.5. Here, we have privileged the geometric interpretation of the Riemann problem. However, one can use the conservative form (16) of the GSD model in polar coordinates to obtain in similar manner as in [30] Eqs. (28)–(30). 3.2.4. Complete solutions of the Riemann problem The relations satisfied by the velocities M and M r and the angles θ and θr for rarefactions and for shocks lead us to conclude that, as in the case of the p-system, there always exists a unique solution of the Riemann problem constituted of at most two transitions (a transition being a rarefaction or a shock). This happens to be true in a whole domain except in a region where the model ceases to be physical (see Fig. 6). Supposing the left state ( M , θ ) is fixed, we can define the 1-shock and 2-shock curves functions, S (1M ,θ ) ( M ) and



S (2M ,θ ) ( M ), by stipulating that for M > M , θ = S (1M ,θ ) ( M ) is the only angle such that ( M , θ ) and ( M , θ) satisfy the



shock relation (30) (which implies that there is an acceptable 1-shock between the states ( M , θ ) and ( M , θ)) and for M < M , θ = S (2M ,θ ) ( M ) is the only angle such that ( M , θ ) and ( M , θ) satisfy the shock relation (30) (which implies that



there is an acceptable 2-shock between the states ( M , θ ) and ( M , θ)). We define R 1( M ,θ ) ( M ) and R 2( M ,θ ) ( M ) in a similar



manner. For M < M , θ = R 1( M ,θ ) ( M ) is the only angle such that ( M , θ ) and ( M , θ) satisfy the 1-rarefaction relation (25)



for σ = 1, and for M > M , θ = R 2( M ,θ ) ( M ) is the only angle such that ( M , θ ) and ( M , θ) satisfy the 2-rarefaction relation

(25) for σ = −1. These functions are drawn in a phase diagram on Fig. 6. Nine cases, depending on the position of ( M r , θr ) relatively to ( M , θ ) in this diagram, are then defined. The first four ones are the general cases, containing two transitions and an intermediate state, as in the study of the p-system (see [24, pp. 306–320]). The other five cases are limit cases with at most one transition. Case 1: ( M r , θr ) is in zone I The solution contains a 1-rarefaction between states ( M , θ ) and ( M ∗ , θ ∗ ), and a 2-shock between states ( M ∗ , θ ∗ ) and ( M r , θr ). Case 2: ( M r , θr ) is in zone II The solution contains a 1-shock between states ( M , θ ) and ( M ∗ , θ ∗ ), and a 2-shock between states ( M ∗ , θ ∗ ) and ( M r , θr ). Case 3: ( M r , θr ) is in zone III The solution contains a 1-shock between states ( M , θ ) and ( M ∗ , θ ∗ ), and a 2-rarefaction between states ( M ∗ , θ ∗ ) and ( M r , θr ). Relations (25) and (30) ensure that solutions of this form always exist when ( M r , θr ) is in zone I, II or III. Case 4: ( M r , θr ) is in zone IV Except when ( M r , θr ) is in a particular, excluded region, there exists a solution containing a 1-rarefaction between states ( M , θ ) and ( M ∗ , θ ∗ ), and a 2-rarefaction between states ( M ∗ , θ ∗ ) and ( M r , θr ). The curve R 1(M ,θ ) intersects the vertical



axis at the point (1, θ0 ). The boundary of the excluded region is then given by the curve θ = R 2(1,θ ) ( M ). 0 In zone IV, as the shape of the front along a rarefaction is universal, it is possible to find the shape of the front (for example the level set {α = t 0 }) by noticing that the intermediate plane front must be both tangent to the (sole) image by a direct isometry of St0 tangent to the front of the state ( M , θ ) and tangent to the sole image by an indirect isometry of St0 tangent to the front of the state ( M r , θr ), where we recall that St0 was defined at the end of Subsection 3.2.2. In particular the front is constructed by a simple geometrical procedure (see Fig. 7). We also remark that as in the analysis of [24], some states in the excluded region (as illustrated in Fig. 6) cannot be joined with ( M , θ ) by the means of two rarefactions. The reason is that lim M →1 ω( M ) < +∞, which means that all curves

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

219

Fig. 8. Half-solutions (with a horizontal wall).

Fig. 9. Complete solutions of the Riemann problem.

R 1 intersect the line M = 1 at a finite point (1, θ0 ) depending on ( M , θ ). When M is high enough, θ0 ≥ θ + π and there is a solution for every physically possible initial data ( M r , θr ) (that is with θr < θ + π ). But when M is small, θ0 < θ + π and in this case, when ( M r , θr ) is above the curve R 2(1,θ ) , the Riemann problem does not admit a solution with two rarefactions. 0 The physical interpretation of this is less clear than in the usual interpretation of the p-system in terms of isentropic gas dynamics, because the GSD approximation does not work as well when M is near 1. Case 5: M r > M and θr = S (1M ,θ ) ( M r )

The solution only contains a 1-shock between states ( M , θ ) and ( M r , θr ). Case 6: M r < M and θr = S (2M ,θ ) ( M r )

The solution only contains a 2-shock between states ( M , θ ) and ( M r , θr ). Case 7: M r < M and θr = R 1( M ,θ )( M r )

The solution only contains a 1-rarefaction between states ( M , θ ) and ( M r , θr ). Case 8: M r > M and θr = R 2( M ,θ ) ( M r )

The solution only contains a 2-rarefaction between states ( M , θ ) and ( M r , θr ). Case 9: M r = M and θr = θ The solution is a uniform planar front of characteristics ( M , θ ). 3.2.5. Geometrical interpretation In the preceding subsection we have constructed the complete solutions of the Riemann problem in the phase space before describing their geometry. There is a more geometrical way to view this reconstruction. Indeed, as it is described by Whitham in [30], a simple shock is the result of the compression of a plane front by a wedge, and a rarefaction the result of the diffraction of a plane front by a corner. Such solutions involve a wall boundary, and in both cases the solution is a plane front near the wall, with a normal vector n parallel to the wall. As a consequence, two such solutions can be glued together along the wall boundary to form a complete solution of (2), provided they have the same velocity M ∗ at the wall. The solution thus constructed is naturally the solution of a Riemann problem. Conversely, the solution of any Riemann problem can be seen as the reunion of two such half-solutions involving a wedge or a corner. For a typical solution constituted of two transitions with an intermediary state, we can split the solution by a virtual boundary at the angle χ w for which θ = χ . Then the restrictions of the solution on both sectors {χ < χ w } and {χ > χ w } are solutions which arise from the interaction of a plane front with a wedge or with a corner. Figs. 8 and 9 illustrate this construction. 4. Algorithm validation The algorithm summarized in Section 2.5 is now evaluated on a set of test cases of increasing complexity. For each of them a reference solution is known or has been introduced previously in Section 3. We begin by checking the numerical

220

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

Table 1 Mesh convergence in the L 1 norm for the cylindrical case and measured CPU time. Mesh

Error M

100 × 100 200 × 200 400 × 400 800 × 800 1600 × 1600 3200 × 3200 6400 × 6400

1.93e−02 9.73e−03 4.90e−03 2.54e−03 1.28e−03 6.41e−04 3.20e−04

α

Order M

Error

0.99 0.99 0.95 0.99 1.00 1.00

6.65e−02 3.04e−02 1.46e−02 7.20e−03 3.52e−03 1.74e−03 8.53e−04

Order 1.13 1.06 1.02 1.03 1.02 1.03

α

CPU time (s) 2.80e−02 7.80e−02 2.91e−01 1.12 4.55 19.12 84.01

convergence order of the scheme by comparison with the smooth radial solution, and we show that the algorithm has optimal computational complexity. Discontinuous solutions are first examined through the Riemann problem where reference explicit expressions are developed in the strong shock limit (i.e. M  1) using the simplified A–M relation

A A0

 =

M0 M

λ∞ ,

with λ∞ = 5.0743 for γ = 1.4.

(32)

The analyses of the compression wedge case (shock–shock perturbation) completes this experimental study of the numerical scheme. We emphasize that all numerical experiments are performed with the full A–M relation (3) but for a large enough Mach number value such that the infinite limit is asymptotically reached. 4.1. Order of the scheme and algorithm complexity We start with a smooth radial configuration. The reference solution is obtained by numerically integrating the radial system (13) with a fourth-order Runge–Kutta scheme in cylindrical and spherical coordinates. The convergence study is made on the discrete L 1 norm.

• Cylindrical expansion. For this two dimensional case, the computational domain is reduced to the square [0, 50] × [0, 50] due to the symmetry of the problem. The initial conditions on α and M are given by the radial solution (13) on the quarter-circle of radius 4.9 centered at (0.0). We impose M = 10 for r < 1 to avoid the singularity at the origin. The boundary conditions are of outgoing type on the edges {x = 50} and { y = 50}, and of rigid wall type on the edges {x = 0} and { y = 0}. Table 1 gives the norm of the error versus the grid spacing and the CPU time measured for a single core of the Intel i3-3220 CPU @ 3.30 GHz on a Linux box with 8 GB of RAM. Note that N is the number of nodes in each direction and that x = y. The discrete order of convergence is defined as: order =

log( E 2 / E 1 ) log( N 1 / N 2 )

,

where E 1 and E 2 are respectively the errors related to the meshes with N 1 × N 1 and N 2 × N 2 elements. Columns three and five of Table 1 show that the expected order of convergence is recovered. Fig. 10 demonstrates that the algorithm has optimal complexity, i.e. it scales as N 2 log( N ), once the number of grid points per direction, N, is large enough.

• Spherical shock. The analysis of the numerical order of the scheme is also carried on in 3D. Using symmetry conditions, the computational domain reduces to the cube [0, 10] × [0, 10] × [0, 10]. The initial front is given by the radial solution on the sphere of radius 2 centered at the origin. As in the 2D case, we impose M = 10 for r < 1 to avoid the singularity r = 0. The boundary conditions are of outgoing type on the edges {x = 10}, { y = 10}, { z = 10}, and of rigid wall type on the edges {x = 0}, { y = 0}, { z = 0}. The results of the numerical scheme are compared with the semi-analytical radial solution considered previously. Table 2 gives the error in the discrete L 1 norm between the two solutions and the measured CPU time. For this 3D case, we do not recover the expected order of convergence. A detailed analysis of these simulations reveals that the maximum error is located in the vicinity of the initial front position, where some early accepted points use a mix of central and upwind terms in the discrete Laplacian operator. As previously explained, and observed here, this can lead to a decrease of the overall accuracy of the numerical scheme. We found that looping a few times over the active list, i.e. locally repeating the update procedure for the test values ϑ and m until the maximal local variation is lower than a threshold (10−8 here) or a maximum iteration count is reached, helps to recover the optimal rate of convergence. Fig. 11 shows that a maximum of three iterations is sufficient. There is nearly no difference with results obtained for five iterations. The local orders of convergence and CPU times are presented in Tables 3 and 4 for the three and five iterations cases respectively. In

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

221

Fig. 10. Algorithm complexity evaluated as the normalized ratio of CPU time to N 2 log( N ) where N is the number of grid points, for a two-dimensional radial smooth case.

Table 2 Mesh convergence in the L 1 norm for the spherical case and measured CPU time. Mesh

Error M

100 × 100 × 100 200 × 200 × 200 300 × 300 × 300 400 × 400 × 400 500 × 500 × 500

8.02e−03 4.53e−03 3.44e−03 2.86e−03 2.53e−03

α

Order M

Error

0.82 0.68 0.64 0.55

1.74e−02 8.92e−03 6.12e−03 4.70e−03 3.85e−03

Order 0.96 0.93 0.92 0.89

α

CPU time (s) 8.9 81.7 297.3 780.9 1852.6

Fig. 11. Mesh convergence in the L 1 norm for the spherical case for α , straight lines, and M dashed lines. Results for one, three and five iterations on points of the active list are represented by stars, squares and circles respectively.

comparison to the initial, one iteration, case reported in Table 2, the CPU time is only increased by a factor about 1.7 and 2.3 respectively. These smooth test cases have proven that our fast-marching like algorithm doesn’t suffer from using points in the NarrowBand and that the expected theoretical order of convergence and optimal algorithm complexity are reached, at least in two dimensions. The three dimensional case is more complex due to an increased importance of Laplacian upwinding. In this case, looping over the active list may be used without penalizing the CPU time. This procedure was found not necessary for the next two dimensional cases reported.

222

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

Table 3 Mesh convergence in the L 1 norm for the spherical case and measured CPU time with three iterations on the active list. Mesh

Error M

100 × 100 × 100 200 × 200 × 200 300 × 300 × 300 400 × 400 × 400 500 × 500 × 500

7.66e−03 3.83e−03 2.62e−03 2.02e−03 1.65e−03

α

Order M

Error

1.00 0.94 0.90 0.91

1.71e−02 8.49e−03 5.67e−03 4.28e−03 3.43e−03

Order

α

CPU time (s) 16.6 142.4 511.8 1256.3 2953.8

1.01 1.00 0.98 0.99

Table 4 Mesh convergence in the L 1 norm for the spherical case and measured CPU time with five iterations on the active list. Mesh

Error M

100 × 100 × 100 200 × 200 × 200 300 × 300 × 300 400 × 400 × 400 500 × 500 × 500

8.09e−03 3.94e−03 2.65e−03 2.01e−03 1.62e−03

α

Order M

Error

1.04 0.98 0.96 0.97

1.73e−02 8.52e−03 5.66e−03 4.25e−03 3.40e−03

Order 1.02 1.01 1.00 1.00

α

CPU time (s) 24.2 201.9 703.0 1721.5 3549.4

Fig. 12. (a) 1-rarefaction case: ( M , θ ) = (10, 0) and ( M r , θr ) = (8.5, 20.97). (b) 2-rarefaction case: ( M , θ ) = (10, 0) and ( M r , θr ) = (12, 23.53). The characteristic curves are displayed for convenience.

4.2. Test cases for the Riemann problem A series of numerical experiments is now conducted to assess in more details the ability of the numerical scheme to solve the two dimensional GSD Riemann problem. In Section 3 we have given the mathematical solutions and we refer the reader to Appendix A for the detail of their explicit construction in the strong shock limit. We compare these reference solutions to the numerical ones, by taking varying values of the initial constant states to cover all configurations of Paragraph 3.2.4. The first tests use the left state ( M , θ ) connected directly to the right state ( M r , θr ) by an elementary wave (either a shock or a rarefaction), i.e. ( M r , θr ) lies on one of the curves R1 , R2 , S1 , S2 issued from ( M , θ ) (see Fig. 6). More general initial data are then considered to cover the different zones of Fig. 6: zone I, zone II, zone III and zone IV, i.e. solutions with an intermediate state. The numerical solution is obtained on the domain [0, 5] × [0, 5] with a 500 × 500 grid. In all cases, the left state ( M , θ ) = (10, 0) is given and we make the right state ( M r , θr ) vary in the plane ( M , θ). We display a comparison between the numerical isolines of α and the exact solution in each configuration. In these figures, the red curve represents the numerical solution, the blue one represents the exact solution and the green lines correspond to the characteristic curves. The following figures gather all the results obtained for the different solutions of the Riemann problems. 4.2.1. Simple waves Simple waves correspond to cases 5 (1-shock), 6 (2-shock), 7 (1-rarefaction) and 8 (2-rarefraction) of Paragraph 3.2.4. They are difficult to compute in practice since numerical rounding on the initial condition may lead to the emergence of an intermediate state. The initial right states for the 1-rarefaction and the 2-rarefaction cases are ( M r , θr ) = (8.5, 20.97) and ( M r , θr ) = (12, 23.53) respectively. Results are drawn in Fig. 12. As expected from the smooth radial case, the numerical and analytical solutions match almost perfectly.

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

223

Fig. 13. (a) 1-shock case with ( M , θ ) = (10, 0) and ( M r , θr ) = (12, −22.45). (b) 2-shock case: ( M , θ ) = (10, 0) and ( M r , θr ) = (8.5, −20.24). Comparison of the numerical solution and the exact solution . Same drawing convention.

Fig. 14. 1-rarefaction + 2-shock: ( M , θ ) = (10, 0) and ( M r , θr ) = (5, 30). Same drawing convention.

The initial right states for the 1-shock and the 2-shock cases are ( M r , θr ) = (12, −22.45) and ( M r , θr ) = (8.5, −20.24) respectively. Here again numerical and analytical results, drawn in Fig. 13, are in excellent agreement. Only a small deviation of the shock–shock position is noticed which can be due to the non-conservative form of the mathematical system (8). 4.2.2. 1-rarefaction, 2-shock This problem corresponds to case 1 of Paragraph 3.2.4. The solution consists of a 1-rarefaction connecting the left state ( M , θ ) = (10, 0) to a constant intermediate state ( M ∗ , θ ∗ ) followed by a 2-shock connecting ( M ∗ , θ ∗ ) to the right state ( M r , θr ) = (5, 30), i.e. ( M r , θr ) is within zone I on Fig. 6. Numerical results are shown on Fig. 14. One can note that the numerical scheme allows a good capture of the intermediate state, retrieving a nearly plane front. A small deviation can be observed at the shock–shock locus. 4.2.3. 1-shock, 2-shock configuration This problem corresponds to case 2 of Paragraph 3.2.4. The solution consists of a 1-shock connecting the left state ( M , θ ) = (10, 0) to a constant intermediate state ( M ∗ , θ ∗ ) followed by a 2-shock connecting ( M ∗ , θ ∗ ) to the right state ( M r , θr ) = (10, −50), i.e. ( M r , θr ) is within zone II on Fig. 6. Numerical results are shown on Fig. 15. As previously noticed, a small deviation is observed at shock–shock locus. 4.2.4. 1-shock, 2-rarefaction configuration This problem corresponds to case 3 of Paragraph 3.2.4. The solution consists of a 1-shock connecting the left state ( M , θ ) = (10, 0) to a constant intermediate state ( M ∗ , θ ∗ ) followed by a 2-rarefaction connecting ( M ∗ , θ ∗ ) to the right state ( M r , θr ) = (20, 10), i.e. ( M r , θr ) is within zone III on Fig. 6. Numerical results are shown on Fig. 16. As expected from the simple wave results, the 2-rarefaction is very well simulated while the 1-shock shows little deviation.

224

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

Fig. 15. 1-shock + 2-shock: ( M , θ ) = (10, 0) and ( M r , θr ) = (10, −50). Same drawing convention.

Fig. 16. 1-shock + 2-rarefaction: ( M , θ ) = (10, 0) and ( M r , θr ) = (20, 10). Same drawing convention.

4.2.5. 1-rarefaction, 2-rarefaction configuration This problem corresponds to case 4 of Paragraph 3.2.4. The solution consists of a 1-rarefaction connecting the left state ( M , θ ) = (10, 0) to a constant intermediate state ( M ∗ , θ ∗ ) followed by a 2-rarefaction connecting ( M ∗ , θ ∗ ) to the right state ( M r , θr ) = (10, 50), i.e. ( M r , θr ) is within zone IV on Fig. 6. Numerical results are shown on Fig. 17. The good behavior of the GSD fast-marching scheme is clearly stated on this smooth case.

4.3. Compression wedge test case

Attention is now paid to discontinuous solutions only. In this section we quantify the previously noticed deviation of the shock–shock locus. We deal here with a well-know test case, studied by many authors [30,2], consisting in the diffraction of an oblique shock by a wedge. It can be seen as an application of the study of simple shocks. For an incoming front of velocity M 0 and angle θ0 = β interacting with a wedge of angle θ w = 0, we are looking for a 2-shock joining the states ( M , θ ) = ( M w , 0) and ( M r , θr ) = ( M 0 , β), the velocity M w of the Mach front and the angle χ2s of the trajectory of the triple point being unknown. The analysis of the structure of 2-shocks then gives us that M w is governed by the relation (30), and χ2s by (28).

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

225

Fig. 17. 1-rarefaction + 2-rarefaction: ( M , θ ) = (10, 0) and ( M r , θr ) = (10, 50). Same drawing convention.

Fig. 18. Shock–shock angle versus angle of incidence for the GSD fast-marching method applied to the compression wedge, in comparison with Whitham’s analytical solution. Angles are in degrees.

In the strong shock limit, Eqs. (31) and (30) can be simplified in the following implicit relationship between β and

⎧  1/2 1 − m2 ⎪ λ∞ ⎪ s ⎪ tan χ = m ⎨ 2 1 − m 2λ∞ ⎪ λ∞ ⎪ ⎪ ⎩ cos β = m + m , 1 + m1+λ∞

M

χ2s

(33)

where m = M 0 < 1 denotes the Mach number ratio on either side of the shock–shock and λ∞ = 5.0743. The subscripts w w and 0 characterize the quantities on the wall and those of the incident (or initial) shock respectively, χ2s is the angle of the shock–shock line with the horizontal axis ( O x). This analytical solution is displayed as a solid curve in Fig. 18. Here, the wedge is chosen aligned with the ( O x) axis to avoid a difficult boundary treatment. The rigid wall condition is then imposed on the edge { y = 0}. The computation was performed on the domain [0, 3] × [0, 3] with a 300 × 300 mesh. The initial shock position is given by

226

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

α0 (x, y ) =

(x − 0.5) cos β − y sin β M0

,

where β ∈ [0, π /2] is the angle between the front normal and the wedge, and M 0 = 10. The values of M and α for all points inside the half space on the left of the initial shock position (i.e. α0 (x, y ) ≤ 0) are already intercepted by the front and easily known since the shock is plane. An outflow condition is imposed on the {x = 3} edge. Special care must be taken for the { y = 3} edge since it corresponds to an inflow condition, so we set the exact solution on it. For a range of incident shock angle, a comparison between this theoretical shock–shock locus and the numerical one obtained by the GSD fast-marching method is shown in Fig. 18. One can notice that the numerical scheme underestimates the value of χ2s at large angles of incidence. The shock–shock locus follows a linear behavior, lying on the tangent (dashed line) to the theoretical curve at M = 1. A mesh refinement does not improve these results. This may be due to the non-conservative formulation we used and could be improved by solving system (2) directly. Above nearly 45 degrees, the algorithm doesn’t develop any shock–shock angle. A better boundary treatment could help improve this behavior. Nevertheless, this are not severe limitations in practice since the GSD model is not expected to give accurate results in such an extreme case [30]. 4.3.1. Synthesis In view of the results obtained in this series of test cases, we note that the GSD fast-marching coincides quite strictly with the exact solution for all cases. A small deviation of the shock–shock locus is nevertheless observed, due to the non-conservative form of the system. On regular smooth cases, the algorithm demonstrated good CPU performances with a complexity of N 2 log( N ) for a plane expanding wave. 5. Conclusion In this paper, we described a new algorithm for solving the Geometrical Shock Dynamics model. In this hyperbolic two equations model, the arrival time, α (x), of a shock wave is given by a nonlinear eikonal equation whose local normalized speed, or Mach number M (x), is governed by a differential equation depending on flow parameters and the front curvature. By reformulating the transport equation on M as a convection–diffusion one, its approximation is made easier and the lengthy discretization of the mean curvature is avoided. The new algorithm follows the fast-marching method on a Cartesian grid but uses trial values when updating the current node. First derivatives, appearing in the eikonal and transport equations, are upwind sided at first order. A key feature of the method is to use a centered discrete Laplacian for the diffusion part, as much as possible, without taking into account the causality condition on α (i.e. using points which are not in the Known set). For this reason, the algorithm is not a fast-marching method in the strictest sense, and we say it has fast-marching like properties. The resulting algebraic nonlinear system is solved by an iterative procedure, fixed point or Newton’s method, at each node. For validation purpose, reference solutions are built for a smooth radial expansion wave and the interaction of two shocks. This latter case is a Riemann problem for which a link is made with the p-system in two dimensions. A detailed analysis of this problem is provided and simple front characterization is given in the strong shock limit for any left and right states. In general, those states are connected by shock or rarefaction waves developing on the front. Numerical experiments have shown that, as expected, the algorithm is first order for smooth fronts and performs remarkably well on discontinuous solutions, although the non-conservative model equations are solved. We demonstrated that the GSD fast-marching method has optimal complexity for a smooth expanding wave, which makes it very fast. The three dimensional smooth case showed that a first order upwind Laplacian reduces the overall accuracy of the scheme unless a few additional iterations over the active list are done. In the future we plan to extend our algorithm to higher order, to deal with immersed rigid bodies and to work on a conservative method. A thorough comparison with the efficient Fast Sweeping method [33] would also be interesting. At last, we infer that this new algorithm could also be valuable to the detonic community, mainly for Detonation Shock Dynamics [26]. Acknowledgements The authors thank the referees for their comments that improved the quality of this paper. N. Lardjane acknowledges the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-12-ASTR-0026. Appendix A. Reference solutions of the Riemann problem in the strong shock limit Based on the geometric description of Section 3, we outline the construction of explicit solutions, in the strong shock limit (M  1), for the 2D Riemann problem. For the sake of simplicity, we give the explicit solution for specific configurations only: simple waves and a particular case with an intermediate state. The key tool is to use the simplified relation (32) between the area section and the Mach number. We consider initial constant left and right states joined at the origin O = (0, 0) and use the polar coordinates x = ρ cos χ and y = ρ sin χ .

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

227

Fig. 19. Example of geometric construction of 1-shock.

A.1. Simple shock wave From the limit form (32) of the A–M relation in the strong shock limit, Eqs. (29)–(30) and (28)–(30) become (see also (33))

⎧ m λ∞ + m ⎪ ⎪ ⎪ ⎨ cos(θ − θr ) = 1 + m λ ∞ +1  1/2 ⎪   λ 1 − m2 ⎪ s ∞ ⎪ , ⎩ tan(χν − θ ) = −1 + ν (ν − 1) m 1 − m 2λ∞

Mr where m = M and ν = 1, 2 for the 1-shock or 2-shock cases. Given 0 < θ − θr < π2 , from the first equation of this algebraic

system, one deduces, by an iterative procedure, the Mach number ratio m on each side of the shock–shock singularity. This singularity is characterized by the angle χνs which is obtained by injecting m in the second equation. For this case, the geometric construction is straightforward and proceeds as follows. As shown in Fig. 19, we consider the two angular sectors separated by the shock line extending from the origin O and having the slope tan χνs . At a given time, one starts from a current position P 1 ≡ (x1 , y 1 ) on the left state and draws the line segment whose normal makes the angle θ with the ( O x) axis. The other endpoint of this segment, P 2 ≡ (x2 , y 2 ), lies on the shock line and so has the polar angle χνs . This process is repeated to draw the front shape at different times to cover the computational domain. The scalar function α , that characterizes the successive front positions, is then (see also Fig. 5)



α (x, y ) =

x cos θ + y sin θ , M x cos θr + y sin θr , Mr

χ ≤ χνs χ ≥ χνs .

A.2. Simple rarefaction wave As mentioned in Section 3.2, the rarefaction locus is a spiral

ρ = t 0 ρ1 (χ ) where t 0 > 0, connecting the two constant states. In the strong shock limit (M  1), this spiral is of logarithmic type, and the angle between its normal and the radial direction is constant (see (18)) and is given by



δ∞ = arctan

1

λ∞

.

In other words, the angles in Proposition 3.1 are written

χ1r,−2 = θ ∓ δ∞ and χ1r,+2 = θr ∓ δ∞ , where the signs ∓ refer to the 1-rarefaction and the 2-rarefaction respectively. Moreover Eq. (25) simplifies to give

   Mr θr − θ = ∓ λ∞ log M

with

π > θr − θ > 0 (see Proposition 3.1).

228

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

α , that characterizes the successive positions of the shock front, is then written (see Fig. 4)

The scalar function

α (x, y ) =

⎧ ⎪ ⎪ ⎨

x cos θ + y sin θ , M √χ

±

K ρ e λ∞ , ⎪ ⎪ ⎩ x cos θr + y sin θr Mr

with K =

χ1r − ,2 cos δ∞ ∓ √λ∞ e M

χ ≤ θ ∓ δ∞ θ ∓ δ∞ ≤ χ ≤ θr ∓ δ∞

, χ ≥ θr ∓ δ∞ .

.

A.3. Complete solutions of the Riemann problem Whatever the left and right initial conditions, the intermediate state (θ ∗ , M ∗ ) can be written in the generic form:

 ∗ ⎧ M ∗−θ =G ⎪ ⎪ θ

R,S ⎨ M

 ∗ ⎪ ⎪ ⎩ θr − θ ∗ = GR,S M ,

(A.1)

Mr

where, in the strong shock limit, the function GR,S is given by

  G (m) = −√λ log m for a rarefaction ∞  R

m+mλ∞  GR,S (m) =  for a shock.  GS (m) = − arccos 1+m1+λ∞

For illustration purpose, we consider only the 1-rarefaction followed by 2-rarefaction case, for which system (A.1) can be solved explicitly. One then checks that

⎧ √ Mr ⎪ ) θr + θ − λ∞ log( M ⎪

∗ ⎪ ⎨θ = ⎪ ⎪ ⎪ ⎩ M∗ =

2



Mr M

θ −θr √ e λ∞

.

Remark A.1. For any other configuration, an iterative procedure is required to determine the intermediate state. For the example under consideration, the front is compound of five angular sectors, according to the values of the polar angle χ (see Fig. 7). The Riemann solution consists in two spiral arcs, characterized respectively by the angles ±δ∞ , connecting the intermediate state to the plane shocks on the left and on the right respectively. Note that in the region between the two spiral arcs, the shock is plane (i.e. M = M ∗ and θ = θ ∗ ) and tangent to both curves. More precisely, the function α is written

⎧ x cos θ + y sin θ

⎪ , ⎪ M ⎪ ⎪ χ ⎪ √ ⎪ ⎪ K ρ e λ∞ , ⎪ ⎪ ⎨ 1 ∗ ∗ α (x, y ) = x cos θ M+∗y sin θ , ⎪ ⎪ ⎪ − √χ ⎪ ⎪ K 2 ρ e λ∞ , ⎪ ⎪ ⎪ ⎪ ⎩ x cos θr + y sin θr , Mr

χ ≤ θ − δ∞ θ − δ∞ ≤ χ ≤ θ ∗ − δ∞ θ ∗ − δ∞ ≤ χ ≤ θ ∗ + δ∞ θ ∗ + δ∞ ≤ χ ≤ θr + δ∞

χ ≥ θr + δ∞ ,

with

K1 =

−δ∞ ) cos δ∞ − (θ√ λ∞ , e M

such that the function

K2 =

cos δ∞ Mr

e

(θ +δ∞ ) √ λ∞

,

α is obviously continuous globally.

References [1] [2] [3] [4] [5] [6]

R.K. Anand, Shock dynamics of strong imploding cylindrical and spherical shock waves with non-ideal gas effects, Wave Motion 50 (2013) 1003–1015. T.D. Aslam, Investigations on detonation shock dynamics, PhD thesis, University of Illinois at Urbana-Champaign, USA, 1996. T.D. Aslam, J.B. Bdzil, D.S. Stewart, Level set methods applied to modeling detonation shock dynamics, J. Comput. Phys. 126 (1996) 390–409. T.D. Aslam, D.S. Stewart, Detonation shock dynamics and comparison with direct numerical simulation, Combust. Theory Model 3 (1) (1999) 77–101. S. Baskar, P. Prasad, Propagation of curved shock fronts using shock ray theory and comparison with other theories, J. Fluid Mech. 5 (2005) 171–198. J.B. Bdzil, T.D. Aslam, Detonation front models: theories and methods, LA-UR 00-942, 1999.

Y. Noumir et al. / Journal of Computational Physics 284 (2015) 206–229

[7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]

229

C. Besset, E. Blanc, Propagation of vertical shock waves in the atmosphere, J. Acoust. Soc. Am. 95 (4) (1994) 1830–1839. J.P. Best, A generalisation of the theory of geometrical shock dynamics, Shock Waves 1 (4) (1991) 251–273. J.P. Best, A generalisation of the theory of geometrical shock dynamics, Shock Waves 2 (1992) 125, Erratum. J.P. Best, Accounting for transverse flow in the theory of geometrical shock dynamics, Proc. R. Soc. A. 442 (1916) (1993) 585–598. J.E. Cates, B. Sturtevant, Shock wave focusing using geometrical shock dynamics, Phys. Fluids 9 (10) (1997) 3058–3068. M.G. Crandall, P.L. Lions, Two approximation solutions of Hamilton–Jacobi equations, Math. Comput. 43 (1984) 1–19. C.M. Dafermos, Hyperbolic Conservation Laws in Continuum Physics, 3rd ed., Grundlehren Math. Wiss., vol. 325, Springer, Heidelberg, New York, 2010. J. Goodman, A. MacFayden, Ultra-relativistic geometrical shock dynamics and vorticity, J. Fluid Mech. 604 (2008) 325–338. W.D. Henshaw, N.F. Smyth, D.W. Schwendeman, Numerical shock propagation using geometrical shock dynamics, J. Fluid Mech. 171 (1986) 519–545. S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Appl. Math. Sci., vol. 153, Springer, 2003. S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. D.W. Schwendeman, A numerical scheme for shock capturing in three dimensions, Proc. R. Soc. Lond. A 416 (1850) (1988) 179–198. D.W. Schwendeman, A new numerical method for shock wave propagation using geometrical shock dynamics, Proc. R. Soc. Lond. A 441 (1912) (1993) 331–341. D.W. Schwendeman, A higher order Godunov method for the hyperbolic equations modeling shock dynamics, Proc. R. Soc. Lond. A 455 (1984) (1999) 1215–1233. D. Serre, Systems of Conservation Laws. 1. Hyperbolicity, Entropies, Shock Waves, Translated from the 1996 French original by I.N. Sneddon, Cambridge University Press, Cambridge, 1999. J.A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science, Cambridge University Press, 1999. J.A. Sethian, P. Smereka, Level set methods for fluid interface, Annu. Rev. Fluid Mech. 35 (2003) 341–372. J. Smoller, Shock Waves and Reaction–Diffusion Equations, Grundlehren Math. Wiss. (Fundamental Principles of Mathematical Sciences), vol. 258, Springer-Verlag, New York, 1994. C.W. Shu, S. Osher, High-order essentially non-oscillatory schemes for Hamilton–Jacobi equations, SIAM J. Numer. Anal. 28 (4) (1991) 907–922. D.S. Stewart, J.B. Bdzil, The shock dynamics of stable multidimensional detonation, Combust. Flame 72 (1988) 311–323. J. Tsitsiklis, Efficient algorithms for globally optimal trajectories, IEEE Trans. Autom. Control 40 (1995) 1528–1538. P.A. Varadarajan, P.L. Roe, Geometrical shock dynamics and engine unstart, in: 41st AIAA Fluid Dynamics Conference and Exhibit, Honolulu, Hawaii, 27–30 June 2011. G.B. Whitham, A new approach to problems of shock dynamics. Part I, Two dimensional problems, J. Fluid Mech. 2 (1957) 146–171. G.B. Whitham, Linear and Nonlinear Waves, Wiley, New York, 1974. L. Yatziv, A. Bartesaghi, G. Sapiro, A fast O ( N ) implementation of the fast marching algorithm, J. Comput. Phys. 212 (2006) 393–399. G.A. Zakeri, Analysis of ray tube area in geometrical shock dynamics, Appl. Math. Comput. 217 (2) (2010) 830–836. H. Zhao, A fast sweeping method for eikonal equations, Math. Comput. 74 (2005) 603–627. Z.-Y. Han, X.-Z. Yin, Geometrical shock dynamics, in: Handbook of Shock Waves, vol. 1, Academic Press, 2001 (Chap. 3.7).