Numerical simulation of cylindrical solitary waves in periodic media

Numerical simulation of cylindrical solitary waves in periodic media


We study the behavior of nonlinear waves in a two-dimensional medium with density and stress relation that vary periodically in space. Efficient approximate Riemann solvers are developed for the corresponding variable-coefficient first-order 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 Korteweg-deVries 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 long-time 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 two-dimensional setting.

We consider nonlinear wave propagation in a 2D periodic medium, modeled by variable-coefficient 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 (close-up) in Figure a and given by:


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.

(a) Initial condition (close-up)
(b) Checkerboard medium
(c) Solution in four different media
Figure 1: Wave behavior in different media. Top left: linear, homogeneous. Top right: nonlinear, homogeneous. Bottom right: linear, periodic. Bottom left: nonlinear, periodic. The quantity plotted is stress ().

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 spatially-varying -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 P-waves and transversal or S-waves [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:


where represents the strain, is the spatially-varying material density, is the stress and is the position vector. Similar to [11], we consider the nonlinear constitutive relation


where plays the role of bulk modulus. Equation (2) with the stress relation (3) admits shock formation. In order to determine entropy-satisfying weak solutions, let us write (2) as a first-order hyperbolic system of conservation laws:


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:


The second medium, shown in figure b, is similar but smoothly (sinusoidally) varying:


We always take .

(a) Checkerboard medium (equation (5))
(b) Sinusoidal type medium (equation (6))
Figure 2: The domains discussed in Section 2.2. Sections of size are shown here, but much larger domains are used in the numerical simulations.

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 low-level 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 spatially-varying -system (4). The Riemann solvers developed here follow the ideas of [10]. In particular, they are based on use of an all-shock solution and -wave decomposition [2].

3.1 Second-order Clawpack discretization

Clawpack is based on Lax-Wendroff 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 first-order fluctuations computed by the normal Riemann solvers described in Section 3.3. The quantities and are second-order corrections that include the fluctuations computed by the transverse Riemann solvers described in Section 3.4, as well as high-resolution approximations to . For details, see [9].

3.2 High-order WENO discretization using SharpClaw

SharpClaw is based on the method of lines approach and uses fifth-order weighted essentially non-oscillatory (WENO) reconstruction in space with fourth-order strong stability preserving (SSP) Runge-Kutta 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 semi-discrete scheme takes the form

SharpClaw employs the same wave propagation Riemann solvers and user interface as Clawpack. In multi-dimensions 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 :




Note that denotes the derivative of with respect to . For concreteness, we consider a Riemann problem in the x-direction, at interface which consists of the hyperbolic system (4) with coefficients


and initial condition


In practice, may represent a cell-average (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:


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 all-shock solver. Shock waves in the Riemann solution correspond to traveling discontinuities that are proportional to (left-going) or (right-going). Equations for an exact all-shock Riemann solution can be derived in a straightforward manner by considering the Rankine-Hugoniot 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 right-going fluctuations, which are quantities that consider the net effect of all left- and right-going waves respectively:


The normal Riemann solver for the y-direction uses the same approach. To solve the Riemann problem at , in place of (11) we use the eigenvectors


The normal flux difference is decomposed in terms of these eigenvectors:

where . Finally, the waves are accumulated into up- and down-going fluctuations:


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 down-going corrections and the vertical traveling waves into right- and left-going corrections. These second-order corrections capture the effect of corner transport.

Transverse corrections to the vertical fluctuations are computed by decomposing the horizontal fluctuations into up- and down-going parts using the eigenvectors (14):


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, variable-coefficient 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 self-convergence 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 self-convergence 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 self-convergence rates as well as relative errors, computed by:


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
Table 1: Self-convergence test for the sinusoidal medium using the Clawpack discretization.

Checkerboard medium

Next we consider self-convergence 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 self-convergence 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.

80 2.935x10
120 1.921x10 1.045
160 1.358x10 1.205
240 7.293x10 1.533
(a) Clawpack
80 1.223x10
120 7.463x10 1.218
160 5.035x10 1.367
240 2.561x10 1.667
(b) SharpClaw
Table 2: Self-convergence test for the checkerboard medium using (a) Clawpack and (b) SharpClaw.

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.

(a) Checkerboard medium, ,
(b) Checkerboard medium, ,
(c) Sinusoidal medium, ,
(d) Sinusoidal medium, ,
Figure 3: Stress at and for solitary wave trains in the sinusoidal and checkerboard media. The slice plots (on the right) show results along the lines (solid black line) and (dashed red line).

Note that the solitary wave front is noticeably non-circular and that the cross-sectional 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.

Figure 4: (a) Stress at in a checkerboard medium with small impedance ratio: . (b) Slice plots showing results along the lines (solid black line) and (dashed red line).

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 p-system (2) is the total energy:


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 .

Figure 5: (a) Stress at in a homogeneous nonlinear medium. (b) Slice plots showing results along the lines (solid black line) and (dashed red line).
Figure 6: Entropy evolution considering a homogeneous nonlinear medium (solid line) and a heterogeneous nonlinear medium (dashed line). The entropy has been normalized relative to the initial entropy.

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.

(a) Entropy evolution up to
(b) Entropy evolution up to
Figure 7: Entropy evolution up to (a) and (b) considering (a) and (b) . The entropy has been normalized relative to the initial entropy.

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:


Here is chosen so as to separate the leading stegoton from everything to the left as well as possible. Specifically, for each , is the right-most local minimum of :


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 outward-going 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 head-on collision. The same effect can be seen in one-dimensional 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 head-on 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 outward-going 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.

(a) Spatial distribution of stress.
(b) Slices at .
Figure 8: 2D collision of cylindrical stegotons for . The stegotons are initially moving toward each other.
(a) 1D collision with same direction.
(b) 1D collision with opposite direction.
Figure 9: Interaction of 1D stegotons. The dashed line shows the stegoton originally at the left propagating on its own.

5 Conclusion

We have performed direct numerical simulations of the multiscale behavior of nonlinear waves in non-dispersive 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 two-dimensional first-order nonlinear hyperbolic PDEs with spatially periodic coefficients, including both smooth and piecewise-constant 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 variable-impedance media. In media with strongly varying impedance, shock formation is strongly suppressed relative to that occurring in homogeneous media, as demonstrated by near-conservation 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.


  1. 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.
  2. 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.
  3. 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(3-4):324–346, 2004.
  4. David I. Ketcheson. High Order Strong Stability Preserving Time Integrators and Numerical Wave Propagation Methods for Hyperbolic PDEs. PhD thesis, Citeseer, 2009.
  5. David I. Ketcheson and Randall J. LeVeque. Shock dynamics in layered periodic media. Communications in Mathematical Sciences, 10(3):859–874, 2012.
  6. 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.
  7. David I. Ketcheson, Matteo Parsani, and Randall J. LeVeque. High-order wave propagation algorithms for hyperbolic systems. SIAM Journal on Scientific Computing. In press.
  8. Randall J. LeVeque and Marsha J Berger. Clawpack software version 4.5. 2011. Url:
  9. R.J. LeVeque. Finite volume methods for hyperbolic problems. Cambridge Univ Press, 2002.
  10. R.J. LeVeque. Finite-volume methods for non-linear elasticity in heterogeneous media. International Journal for Numerical Methods in Fluids, 40(1-2):93–104, 2002.
  11. R.J. Leveque and D.H. Yong. Solitary waves in layered nonlinear media. SIAM Journal on Applied Mathematics, 63(5):1539–1560, 2003.
  12. Manuel Quezada de Luna. Nonlinear wave propagation and solitary wave formation in two-dimensional heterogeneous media. Master’s thesis, King Abdullah University of Science and Technology (KAUST), 2011.
  13. 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.
  14. 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.
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.