Numerical simulation of cylindrical solitary waves in periodic media
Abstract
We study the behavior of nonlinear waves in a twodimensional medium with density and stress relation that vary periodically in space. Efficient approximate Riemann solvers are developed for the corresponding variablecoefficient firstorder hyperbolic system. We present direct numerical simulations of this multiscale problem, focused on the propagation of a single localized perturbation in media with strongly varying impedance. For the conditions studied, we find little evidence of shock formation. Instead, solutions consist primarily of solitary waves. These solitary waves are observed to be stable over long times and to interact in a manner approximately like solitons. The system considered has no dispersive terms; these solitary waves arise due to the material heterogeneity, which leads to strong reflections and effective dispersion.
1 Introduction
Solutions of first order hyperbolic partial differential equations (PDEs) generically lead to formation of shock singularities and subsequent entropy decay. In contrast, nonlinear wave equations with dispersive terms, such as the KortewegdeVries equation (KdV), may exhibit solitary wave solutions [14] . These solutions are remarkable in the context of nonlinear PDEs since they propagate without changing shape and, in many cases, interact only through a phase shift. Thus the longtime solution behavior of nonlinear waves depends critically on the presence of dispersive regularizing terms.
Santosa and Symes [13] showed that solutions of the linear wave equation in a medium with periodically varying coefficients exhibit dispersive behavior, even though the equation contains no dispersive terms. LeVeque and Yong [11] found that solutions of the 1D system in a periodic layered medium form solitary waves (referred to as stegotons due to their discontinuous shape) similar to those arising in nonlinear dispersive wave equations like KdV. This is remarkable since the equations they considered are first order hyperbolic PDEs with no dispersive terms. Recently it has been observed that solutions of first order hyperbolic systems with periodic coefficients may be characterized by shocks or by solitary waves, depending on the medium and the initial conditions [5]. Further experiments in [4, 5] indicate that solitary waves may arise quite generally in the solution of nonlinear hyperbolic PDEs with periodically varying coefficients. Those results have revealed new behaviors of nonlinear waves in 1D heterogeneous materials, and it is natural to ask if such behaviors persist in higher dimensions. The present work, which builds on [12], extends the aforementioned studies to a fully twodimensional setting.
We consider nonlinear wave propagation in a 2D periodic medium, modeled by variablecoefficient first order hyperbolic PDEs. Since 1D solitary waves are often found to be unstable when extended to higher dimensions, a principal question is whether waves like stegotons are stable in 2D, and whether they form in materials that vary in multiple spatial directions. In order to investigate these questions in a general setting, we have chosen to study the nonlinear wave equation
This is perhaps the simplest multidimensional nonlinear wave model that is general in the sense of allowing for wave propagation in all directions. In 1D, it is commonly referred to as the system due to its connection with Lagrangian gas dynamics.
In order to introduce the object of our study, four representative examples of the behavior of cylindrical wavefronts are depicted in Figure c. Each solution shown results from the same initial condition: a small, cylindrically symmetric Gaussian pulse centered at the origin, shown (closeup) in Figure a and given by:
(1) 
Since the solution is symmetric under reflection about the  and axes, a single quadrant is sufficient to characterize the full solution. The top left quadrant of Figure c shows the solution obtained in a linear, homogeneous medium: a smooth pulse expanding at the sound speed of the medium. The top right quadrant corresponds to a nonlinear, homogeneous medium; the leading edge of the pulse steepens into a shock wave. The bottom right quadrant corresponds to a linear, heterogeneous medium shown in Figure b that is composed of alternating square homogeneous regions in a checkerboard pattern. The front travels more slowly than in the homogeneous case, due to the effect of reflections. Finally, the bottom left quadrant shows the subject of the present study. The medium is both nonlinear and heterogeneous (with the same checkerboard structure). In this case, a train of cylindrical pulses forms. As we will see, these pulses appear to behave as solitary waves and we refer to them as cylindrical stegotons. Note that they do not exhibit exact cylindrical symmetry, due to the lack of such symmetry in the medium.
In the remainder of the paper, we investigate the conditions that lead to these cylindrical stegotons, as well as their properties and behavior. The model equations and materials are introduced in Section 2.
Accurate numerical modeling of nonlinear waves in multidimensional heterogeneous media is challenging. Care must be taken to properly handle material inhomogeneities (especially discontinuities) and shocks. We use finite volume methods based on approximate Riemann solvers implemented in Clawpack [8] and SharpClaw [7], New approximate Riemann solvers for the heterogeneous 2D system are discussed in Section 3. Because very large grids (with approximately 7 billion unknowns) are required to resolve this multiscale problem over times and distances sufficient to observe solitary wave formation, we have used the parallel PyClaw framework [6] in order to run on 16,384 cores of the Shaheen system at KAUST.
Motivated by known behaviors of more traditional solitary waves, and also by studies of 1D stegotons, we investigate qualitatively the following questions:

What kinds of media and initial conditions give rise to solitary waves?

How do cylindrical stegotons interact with each other?

What is the role of shock formation in 2D periodic media with strongly varying impedance?
These questions are addressed through further numerical experiments in Section 4.
2 The 2D spatiallyvarying system
2.1 Governing equations
As mentioned in the introduction, we have chosen to study the system due to its relative simplicity and generality. This system can be viewed as a simplification of models governing more complex wave behavior, such as that of elastic waves. In a solid, an elastic wave is composed of longitudinal or Pwaves and transversal or Swaves [9]. If we assume the stress is hydrostatic (i.e., there is no shear stress and the extensional stress components are equal), the propagation of elastic waves can be modeled by:
(2) 
where represents the strain, is the spatiallyvarying material density, is the stress and is the position vector. Similar to [11], we consider the nonlinear constitutive relation
(3) 
where plays the role of bulk modulus. Equation (2) with the stress relation (3) admits shock formation. In order to determine entropysatisfying weak solutions, let us write (2) as a firstorder hyperbolic system of conservation laws:
(4a)  
where  
(4b) 
Here and are the  and components of velocity, is the vector of conserved quantities, and are the components of the flux in the  and directions, respectively. This form arises naturally through consideration of kinematics and Newton’s second law, and leads to the correct jump conditions across shocks.
2.2 Periodic Media
We consider two types of periodic variation in the density and the bulk modulus . The period of the medium is always taken to be unity. The first medium, shown in figure a, consists of a checkerboard pattern:
(5) 
The second medium, shown in figure b, is similar but smoothly (sinusoidally) varying:
(6a)  
(6b) 
We always take .
3 Numerical Discretization
In this section we describe the numerical methods used to solve the 2D spatially varying system (4). For all computations in this work we use the algorithms implemented in Clawpack and SharpClaw. Both are finite volume methods for solving hyperbolic PDEs based on solving Riemann problems. For the simulations presented here, both Clawpack and SharpClaw are driven through PyClaw [6] which is a lightweight Python framework that calls the lowlevel Fortran routines of Clawpack and SharpClaw, and also interfaces with PETSc [1] to provide parallelism.
After briefly introducing the discretization schemes, we focus on their key ingredient, the approximate solution of the Riemann problem for the spatiallyvarying system (4). The Riemann solvers developed here follow the ideas of [10]. In particular, they are based on use of an allshock solution and wave decomposition [2].
3.1 Secondorder Clawpack discretization
Clawpack is based on LaxWendroff discretization combined with TVD limiters, and is second order accurate in space and time [8]. The multidimensional Clawpack algorithms used here require the propagation of waves in both the normal and transverse directions at each cell edge. The Clawpack discretization takes the form
where and are firstorder fluctuations computed by the normal Riemann solvers described in Section 3.3. The quantities and are secondorder corrections that include the fluctuations computed by the transverse Riemann solvers described in Section 3.4, as well as highresolution approximations to . For details, see [9].
3.2 Highorder WENO discretization using SharpClaw
SharpClaw is based on the method of lines approach and uses fifthorder weighted essentially nonoscillatory (WENO) reconstruction in space with fourthorder strong stability preserving (SSP) RungeKutta integration in time [7]. SharpClaw requires propagation of waves only in the direction normal to each edge.
First, a WENO reconstruction of the solution is computed from the cell averages to give high order accurate point values and just to the left and right (respectively) of each cell interfaces . A Riemann solution is computed at each interface based on the reconstructed values there. The resulting fluctuations are used to update the adjacent cell averages. An additional term appears that is proportional to . For conservative systems like (4), this term can be conveniently computed in terms of a fictitious internal Riemann problem in each cell [7]. The semidiscrete scheme takes the form
SharpClaw employs the same wave propagation Riemann solvers and user interface as Clawpack. In multidimensions SharpClaw requires propagation of waves only in the normal direction to each edge.
3.3 Normal Riemann solvers
We assume that the density and bulk modulus are constant within each grid cell: . In the case of smoothly varying media, this is an approximation. Then the system of conservation laws (4) can be written in the following quasilinear form within cell :
(7) 
where
(8) 
Note that denotes the derivative of with respect to . For concreteness, we consider a Riemann problem in the xdirection, at interface which consists of the hyperbolic system (4) with coefficients
(9) 
and initial condition
(10) 
In practice, may represent a cellaverage (in Clawpack) or a reconstructed value at the cell interface (in SharpClaw). The eigenvectors of are
with corresponding eigenvalues . Here is the impedance and is the sound speed.
In the linear case, each wave in the Riemann solution is a discontinuity proportional to the corresponding eigenvector in the material carrying the wave. Thus the solution can be found by decomposing the difference in terms of the following three eigenvectors:
(11) 
In the nonlinear case, the Riemann solution also consists of three waves, one of which is a stationary shear wave with zero velocity. The other two waves may be rarefaction waves or shock waves, but the solution cannot include transonic rarefactions, since the system is derived in a Lagrangian frame. For reasons of computational efficiency, we will use an approximate allshock solver. Shock waves in the Riemann solution correspond to traveling discontinuities that are proportional to (leftgoing) or (rightgoing). Equations for an exact allshock Riemann solution can be derived in a straightforward manner by considering the RankineHugoniot conditions. However, these lead to a coupled nonlinear system that is expensive to solve numerically. Instead, we again follow [10] and approximate the solution further by replacing the exact wave speeds with the sound speeds in each cell. The waves themselves are found following the wave approach [2], by decomposing the jump in the normal flux in terms of the eigenvectors (11):
where . The second wave is not needed in the numerical solution, since it has velocity zero and thus does not affect the solution value in either cell. The other two waves are accumulated into left and rightgoing fluctuations, which are quantities that consider the net effect of all left and rightgoing waves respectively:
(12)  
(13) 
The normal Riemann solver for the ydirection uses the same approach. To solve the Riemann problem at , in place of (11) we use the eigenvectors
(14) 
The normal flux difference is decomposed in terms of these eigenvectors:
where . Finally, the waves are accumulated into up and downgoing fluctuations:
(15)  
(16) 
As observed in [10] for the 1D system, this approach leads to efficient approximate Riemann solvers that are conservative and provide good accuracy at least for weakly nonlinear problems.
3.4 Transverse Riemann solvers
The multidimensional Clawpack algorithm makes use of a transverse Riemann solver that decomposes the horizontally traveling waves into up and downgoing corrections and the vertical traveling waves into right and leftgoing corrections. These secondorder corrections capture the effect of corner transport.
Transverse corrections to the vertical fluctuations are computed by decomposing the horizontal fluctuations into up and downgoing parts using the eigenvectors (14):
(17) 
The speed of the waves in the vertical direction determines the amount of horizontal fluctuation that must be added to the vertical one. These speeds are given by the eigenvalues of , which are , and . Finally, the corrections to the vertical fluctuations are:
The transverse corrections to the horizontal fluctuations are obtained in an analogous way. Those corrections are:
3.5 Accuracy tests
For the nonlinear, variablecoefficient system studied in the present work, exact solutions are not available and error and convergence estimates are difficult to obtain because the degree of regularity of solutions is not known. In this section some measure of the accuracy of the numerical solutions is obtained by conducting selfconvergence tests on small problems with the same qualitative features as the problems of interest. The principal purpose of these tests is to determine the grid resolution necessary to ensure small relative errors and qualitatively accurate results. Detailed numerical analysis of the schemes’ convergence behavior for these problems is beyond the scope of this work. In all computations here and in the rest of this work, a uniform cartesian grid is used with .
Sinusoidal medium
Considering the nonlinear constitutive relation (4) and the material properties , we perform a selfconvergence study in the sinusoidal medium (6) using the Clawpack discretization. The initial condition is given by (1). The computational domain is restricted to the positive quadrant by imposing reflecting boundary conditions at the left and bottom. The stress at on a grid with is used as reference solution. Table 1 shows selfconvergence rates as well as relative errors, computed by:
(18) 
where is the computed solution and is the solution on the fine grid. The observed convergence rate is roughly second order.
Error  Rate  

80  1.737x10  — 
120  8.943x10  1.638 
160  5.306x10  1.814 
240  2.235x10  2.132 
Checkerboard medium
Next we consider selfconvergence of the solution of (4) in the checkerboard medium (5) using the Clawpack discretization and SharpClaw. We use the nonlinear constitutive relation given by (3) and material parameters . The initial condition is given by (1). The stress is computed at and the solution on a grid with is taken as the reference solution for the selfconvergence study. Tables a and b show the convergence rates and relative errors following (18) using Clawpack and SharpClaw respectively. The achieved order of convergence is between 1 and 2 for both discretizations, which seems reasonable given the discontinuous coefficients. The SharpClaw solution is noticeably more accurate, even though the two methods exhibit similar convergence rates.


4 Computational results
4.1 Formation of solitary waves
We now investigate in detail the solitary wave trains mentioned in Section 1. These waves are found to arise in both the checkerboard and the sinusoidal domain. We take a symmetric Gaussian hump (1) as initial stress and solve the system (2) with the nonlinear constitutive relation (3). Reflecting boundary conditions are used at the left and bottom boundaries in order to take advantage of symmetry and compute only in the first quadrant. The mesh is uniform with . Based on the accuracy studies above and additional tests, we expect this to be sufficient computational resolution to obtain qualitatively accurate solutions.
For the checkerboard domain, we take and compute the solution using SharpClaw with . For the sinusoidal domain, we take and compute the solution using Clawpack with . Figure d shows the stress at and for both materials, including slices along the lines and . The persistence of the solitary waves after long times strongly suggests that they are stable, attracting solutions.




Note that the solitary wave front is noticeably noncircular and that the crosssectional shape of the solitary wave pulse varies with respect to the angle of propagation.
4.2 Entropy evolution and shock formation
As discussed in the introduction, solutions of nonlinear hyperbolic PDEs generically develop shock discontinuities, leading to irreversibility and entropy decay. It has been observed that spatially varying materials can inhibit the formation of shocks [3, 5]. Meanwhile, dispersive nonlinear wave equations commonly exhibit solutions that are regular for all times. Here we are solving a PDE without dispersive terms; nevertheless, the spatially varying coefficients create reflections, which yield effective dispersion.
No visible shocks appear in the solutions shown in Figure d. On the other hand, figure 4 shows that if the impedance ratio is small enough, the solution develops shocks. This is in qualitative agreement with the theory for 1D systems developed in [5]. We now study shock formation in the 2D stegotons by considering the entropy evolution and the reversibility of the solution.
An entropy function for a hyperbolic PDE is a function that is conserved while the solution is smooth, but decreases in time when shocks form [9]. Thus, entropy decay is a signature of shock formation. A suitable entropy function for the 2D psystem (2) is the total energy:
(19) 
We now compare the solution in a homogeneous domain with (shown in Figure 5) to that obtained in the sinusoidal domain with (shown in Figure c). The difference in the behavior is clear: the heterogeneity introduces reflections that effectively yield dispersion, breaking the initial profile into solitary waves. For the homogeneous domain, no reflections are present; consequently, there is no dispersion to prevent shock formation. In Figure 6, we see the entropy evolution for both cases. For the simulation with homogeneous domain, the entropy starts to decrease as soon as the shock starts to happen. On the other hand, the entropy using the sinusoidal type medium remains almost constant. In both cases we use a resolution given by .
Indeed, in both cases, there is loss of entropy due to numerical dissipation, but this converges to zero, as the grid is refined. To better determine if there is shock formation in the sinusoidal case, we present, in figure a, the normalized entropy evolution considering different resolutions given by . We see the entropy lost is indeed small. However, it is evident that the rate of entropy lost changes before making the entropy lost not to converge to zero as the grid is refined. This suggests the existence of a shock. To corroborate this, we focus up to , study the entropy for even finer grids and see its convergence. In figure b, we show the entropy evolution up to for different resolutions given by . We see the entropy at converges to a value different than the initial entropy; i.e., there is a loss of entropy, which corroborates there existence of a shock.
4.3 Interaction of cylindrical stegotons
In order to study the interaction of these waves, we start with the solution depicted in Figure b. We attempt to isolate the leading pulse in that solution by setting the solution to zero outside of a narrow band. Let denote the solution of (4) and denote the corresponding isolated solution. Then is given by:
(20) 
Here is chosen so as to separate the leading stegoton from everything to the left as well as possible. Specifically, for each , is the rightmost local minimum of :
(21) 
where is the discretization of the stress . This isolation is imperfect since it seems clear that the tails of the solitary waves still overlap. Obtaining completely separated solitary waves would require even greater computational resources or a different modeling approach. By extracting just the leading pulse in this manner at two different times ( and ), we obtain a pair of nearly isolated solitary waves. The interaction simulation is initialized using the sum of these two solutions as initial condition, but with the velocity negated in the solution so that it will propagate inward. Figure a shows surface plots of the solution at . Corresponding slices at are shown in Figure b (solid line). For comparison, we also simulate just the outer pulse (without the inner pulse) over the same time (dotted line).
The results are typical of solitary waves; the two pulses retain their identity after the interaction. Furthermore, the position of the outwardgoing pulse after the interaction is essentially the same as in the simulation with no interaction. The apparent lack of a noticeable phase shift seems to be a result of the very short interaction time for this headon collision. The same effect can be seen in onedimensional stegoton interactions, as shown in Figure b. Two stegotons traveling in the same direction exhibit a phase shift after interaction, but no phase shift is discernible after a headon interaction. In order to simulate the interaction of two cylindrical stegotons traveling in the same direction (and the associated phase shift), much greater computational resources would be required.
Some small oscillations are visible in the final solution in Figure b. Since similar oscillations are seen in the solution obtained with just the outwardgoing pulse, it seems that these may be attributed to the fact that the tails of the solitary waves were not accurately captured in the initial condition. However, it might also be that some small oscillations are radiated during the interaction.




5 Conclusion
We have performed direct numerical simulations of the multiscale behavior of nonlinear waves in nondispersive periodic materials. The largest simulations involve unknowns and were performed on 16K cores of a BlueGene/P system. Our computational results indicate that solitary waves may arise in solutions of twodimensional firstorder nonlinear hyperbolic PDEs with spatially periodic coefficients, including both smooth and piecewiseconstant media. As in 1D, these waves apparently result from the combination of effective (material) dispersion and nonlinear steepening. The effective dispersion is a macroscopic effect caused by reflections in variableimpedance media. In media with strongly varying impedance, shock formation is strongly suppressed relative to that occurring in homogeneous media, as demonstrated by nearconservation of entropy. The cylindrical solitary waves that we have studied appear to behave approximately like solitons in their interactions.
Studying the properties of these waves in detail is difficult due to the multiscale nature of the problem. Future work may include simulations using massively parallel adaptive mesh refinement or more computational resources in order to more fully isolate the solitary waves and study longer interactions. A complementary tool in understanding these waves is the application of homogenization theory to multidimensional nonlinear hyperbolic systems. We have focused on localized perturbations in media with uniform sound speed and strongly varying impedance. Many other geometries, for both the medium and the perturbation, could also be considered. All of these topics are the subject of ongoing work.
References
 Satish Balay, Jed Brown, Kris Buschelman, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Barry F. Smith, and Hong Zhang. PETSc Web page, 2012. http://www.mcs.anl.gov/petsc.
 D.S. Bale, R.J. LeVeque, S. Mitran, and J.A. Rossmanith. A wave propagation method for conservation laws and balance laws with spatially varying flux functions. SIAM Journal on Scientific Computing, 24(3):955–978, 2003.
 J.P. Fouque, Josselin Garnier, and A. Nachbin. Shock structure due to stochastic forcing and the time reversal of nonlinear waves. Physica D: Nonlinear Phenomena, 195(34):324–346, 2004.
 David I. Ketcheson. High Order Strong Stability Preserving Time Integrators and Numerical Wave Propagation Methods for Hyperbolic PDEs. PhD thesis, Citeseer, 2009.
 David I. Ketcheson and Randall J. LeVeque. Shock dynamics in layered periodic media. Communications in Mathematical Sciences, 10(3):859–874, 2012.
 David I. Ketcheson, Kyle T. Mandli, Aron Ahmadia, Amal Alghamdi, Manuel Quezada de Luna, Matteo Parsani, Matthew G. Knepley, and Matthew Emmett. PyClaw: Accessible, Extensible, Scalable Tools for Wave Propagation Problems. SIAM Journal on Scientific Computing, 34(4):C210–C231, 2012.
 David I. Ketcheson, Matteo Parsani, and Randall J. LeVeque. Highorder wave propagation algorithms for hyperbolic systems. SIAM Journal on Scientific Computing. In press.
 Randall J. LeVeque and Marsha J Berger. Clawpack software version 4.5. 2011. Url: www.clawpack.org.
 R.J. LeVeque. Finite volume methods for hyperbolic problems. Cambridge Univ Press, 2002.
 R.J. LeVeque. Finitevolume methods for nonlinear elasticity in heterogeneous media. International Journal for Numerical Methods in Fluids, 40(12):93–104, 2002.
 R.J. Leveque and D.H. Yong. Solitary waves in layered nonlinear media. SIAM Journal on Applied Mathematics, 63(5):1539–1560, 2003.
 Manuel Quezada de Luna. Nonlinear wave propagation and solitary wave formation in twodimensional heterogeneous media. Master’s thesis, King Abdullah University of Science and Technology (KAUST), 2011.
 F. Santosa and W.W. Symes. A dispersive effective medium for wave propagation in periodic composites. SIAM Journal on Applied Mathematics, 51(4):984–1005, 1991.
 N.J. Zabusky and M.D. Kruskal. Interaction of "solitons" in a collisionless plasma and the recurrence of initial states. Physical Review Letters, 15(6):240–243, 1965.