Rogue waves: New forms enabled by GPU computing

Physics Letters A 378 (2014) 2377–2381

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Rogue waves: New forms enabled by GPU computing Christopher Chabalko a,∗ , Ayan Moitra a , Balakumar Balachandran b a b

M0102 Glenn Martin Mezzanine (089), University of Maryland, College Park, MD 20742-5031, United States 2181 Glenn Martin Hall (088), University of Maryland, College Park, MD 20742-5031, United States

a r t i c l e

i n f o

Article history: Received 24 January 2014 Received in revised form 6 June 2014 Accepted 7 June 2014 Available online 12 June 2014 Communicated by C.R. Doering Keywords: Nonlinear Schrödinger equation Rogue waves GPGPU

a b s t r a c t A method to determine parameters governing periodic Riemann theta function rogue-wave solutions to the nonlinear Schrödinger equation is presented. A map of parameter values leading to candidate solutions is developed. In addition to candidate solutions, an overview of qualitative aspects of the solution space can be gained from this map. Based on these findings, several new extreme wave solutions are presented. Although the computations required to determine the map are quite demanding, it is shown that these computations can be efficiently accelerated with a parallel computing architecture. A general purpose computing on a graphics processor unit (GPGPU) implementation yielded a 400× acceleration over a single threaded high level implementation. This acceleration enabled exploration and examination of the solution space, which otherwise would not have been possible. In addition, the solution methodology presented here can be extended to explore other classes of solutions. © 2014 Elsevier B.V. All rights reserved.

1. Introduction The nonlinear Schrödinger equation (NLSE) has been used to model extreme waves in many domains. The NLSE admits the Benjamin–Feir (modulational) instability [1], which has been proposed as one of the mechanisms for rogue-wave formation [2]. For this reason, it is often used as a model for extreme wave behavior. Several families of analytical solutions have been determined for the NLSE. Much of the motivation in developing these solutions is to model physical rogue waves. These waves can appear on the surface of the ocean, in fiber optical systems, or in other domains as well. Some of the first analytical solutions to the NLSE were presented by Tracy [3]. Other contributions to determining analytical solutions of the NLSE include those due to Akhmediev and Korneev [4], who determined a family of single parameter solutions. Based on finite gap integration, Smirnov [5] constructed a family of two-gap solutions and derived conditions under which they behave as rogue waves. He extended this work to periodic two phase and three phase solutions in another study [6]. He was able to show that these solutions lead to rogue-wave type behavior when the eigenvalues are close and have large imaginary components. In addition, Smirnov [7] presented rogue wave type elliptic solutions to the NLSE. These solutions, based on three parameters, are distinct from Akhmediev’s previously discovered elliptic solutions,

*

Corresponding author. Tel.: +1 301 405 2377.

http://dx.doi.org/10.1016/j.physleta.2014.06.013 0375-9601/© 2014 Elsevier B.V. All rights reserved.

which are not considered to be rogue-wave type solutions. Furthermore, the three parameter elliptic rogue-wave solutions degenerate to Akhmediev breathers and Peregrine solitons for certain choices of parameters. A review of nonlinear optical waves, including exact solutions to the nonlinear Schrödinger equation, nonlinear interference, and soliton behavior in dispersive media is available in the book by Akhmediev and Ankiewicz [8]. Different groups have determined other families of rogue-wave type solutions to the standard NLSE. Notably, Akhmediev, Soto-Crespo, and Ankiewicz [9] identify the interference of Akhmediev breathers (ABs) as leading to a type of rogue-wave solution. They show that properly phased AB collisions can result in rogue waves and suggest it as a method to explain and possibly provoke rogue waves in optical fibers. Akhmediev, Ankiewicz, and Soto-Crespo [10] have also determined a family of rational solutions to the NLSE. The rational solutions are determined by taking a modified Darboux transform of a specially chosen seed solution. They have successfully tested for the presence of rational solutions in a randomly perturbed wave field. A system governed by the NLSE was excited with a plane wave with random perturbations. Regions of large amplitudes were identified and found to match almost identically to the envelope predicted by the rational solutions. Several groups have verified the analytically predicted solutions with experimental results. The Peregrine soliton, which is a limiting form of several families of analytical rogue-wave solutions, has been studied in a fiber optic cable [11]. In turn, several analytically predicted extreme waves have been demonstrated experimentally in optical fibers [12] and water wave tanks [13,14]. In each case,

2378

C. Chabalko et al. / Physics Letters A 378 (2014) 2377–2381

the observed rogue waves have been modeled after solutions to the NLSE. Finally, Dysthe [15] has introduced a higher order approximation to the wave equation, and this equation is called the Dysthe equation. This equation is considered to provide a more accurate model of extreme wave behavior under certain conditions. To further contribute to the advancement of the physical realization of rogue-wave solutions to the NLSE, in the present work, the authors offer a method to computationally determine parameters of theta functions which result in as of yet undiscovered rogue waves. While the theta function form of rogue waves is well known [3], the realization of particular theta function solutions enabled by this work enables new physical forms of extreme waves not yet available in the literature. One key aspect of the method, the GPU map, provides an overview of the possible solution parameters. This presentation format provides the analyst with a broad view of the parameter space. Such a broad view can stimulate intuition about governing features of the space, and allow the investigation of sparse regions of the space which would be otherwise difficult to discover. The scaled NLSE takes the form [16]

iut − u xx + 2σ |u |2 u = 0

(1)

where u (x, t ) is the complex wave envelope field, t is time, x is the √ spatial variable, i = −1, the subscripts indicate the corresponding partial derivatives, σ = −1 yields the focusing case, and σ = 1 yields the defocusing case. The inverse scattering transform was used by Shabat and Zakharov [17] to develop analytic solutions of the NLSE. The solution space can be viewed as having a nonlinear Fourier structure, which is comprised of stable and unstable modes. Nonlinear interactions can occur between these modes based on associated eigenvalues [18]. The unstable modes are potential “rogue-wave” solutions. These modes can be expressed in terms of Riemann-theta functions. In this work, the authors provide a procedure to solve for certain unstable “rogue-wave” modes. Nayfeh and Balachandran [19] showed how externally forced nonlinear systems may exhibit chaotic dynamics and wave behaviors in a variety of physical systems. Given the current state of understanding of solutions of the NLSE, as of yet unknown roguewave solutions may be critical to further the understanding of instabilities and extreme behaviors of many systems. Several solutions to the NLSE are already known, and motivate the search for more solutions. The Peregrine soliton is a well-known extreme wave solution [20]. Akhmediev, Ankiewicz, and Taki [21] have detailed exact expressions for rational solutions, which can be used to describe rogue waves. Ma and Ablowitz [22] have provided a solution methodology for obtaining spectral solutions for periodic boundary conditions for both the focusing and defocusing cases. Here, the authors consider solutions to the focusing NLSE for periodic boundary conditions. A new approach is offered to determine parameters that correspond to rogue-wave solutions. Furthermore, the fundamental steps outlined in this parameter selection approach are not restricted to theta function based solutions. Major contributions of this effort are the following: (a) a description of a general procedure to choose likely Riemann theta function solutions to the NLSE and (b) exploration of the parameter space governing possible wave solutions through GPU computing. The remainder of the article is organized as follows. In the next section, a map of candidate solutions is developed. Next, the predictor–corrector solution procedure is described. Based on the predictor-correct procedure, several rogue-wave solutions are determined and presented. Following that, concluding remarks are collected and presented together at the end.

Fig. 1. Flow chart illustrating the procedure to discover Riemann theta function described rogue waves.

2. Rogue-wave formations: computational formulation, solutions, and features The goal of this work is to present a procedure, based on a predictor–corrector style framework, which allows the discovery of a certain form of rogue-wave solution to the NLSE. The steps involved are summarized in the flowchart given in Fig. 1 and explained with the aid of Eqs. (2) to (18) included in this section. Through GPU computing, this procedure allows an investigator to explore the parameter space, for possible solutions of Eq. (2). In the predictor step, a map of parameters which result in periodic functions, ψ(x, 0), is generated. A particular combination of parameters that has a high likelihood of resulting in a rogue wave can then be determined. In the corrector step, the parameters are conclusively refined by solving the spectral eigenvalue problem, shown in Eq. (12), based on the initial guess. In the verification step, a candidate solution is formed by substituting the corrected parameters into Eq. (2). The candidate solution is verified by numerically evaluating the NLSE to determine a residual of zero. The candidate solution is put through two numerical tests, a successful completion of which leads to the acceptance of the candidate solution as a solution. These steps are described in more detail below. 2.1. Predictor Space periodic spectral solutions to the NLSE can be described by

ψ(x, t ) = A

Θ(x, t |τ , δ − ) 2i A 2 T e Θ(x, t |τ , δ + )

(2)

where Θ(x, t |τ , δ ± ) is a Riemann theta function [3,23,24]. A single unstable mode can be considered by taking Θ(x, t |τ , δ ± ) as a two-dimensional theta function defined as



Θ x, t |τ , δ

±



=

∞ 

∞ 



exp i

m1 =−∞ m2 =−∞

+

2 

±

mn δn +

n =1

2 

mn K n x +

n =1 2 2  



m j mk τ jk

2 

mn Ωn T

n =1

(3)

j =1 k =1

The parameters governing the theta function (K n , Ωn , and δ ± ) are defined in terms of five spectral parameters A, λ R , λ I , 0 , and θ . Following the notation used in earlier work [25], the spectral parameters are defined as

1 = 0 e iθ , λ1 = λ R + i λ I ,

2 = 1∗ ,

σ1 = 1 , ∗

λ2 = λ1

σ2 = −1

(4) (5)

C. Chabalko et al. / Physics Letters A 378 (2014) 2377–2381

K n = −2



Ωn = 2λn K n  1 1 δn± = π + i ln λn ∓ K n + i ln σn λn − (−1)n K n

τ11 =

1 2

+

A 2 + λn2 ,





i

π

τ12 = τ21 =

ln i

2π

2

K 12

1 

ln

,

τ22 =

1 + λ 1 λ2 + 1 + λ 1 λ2 −

1 2

+

1 K K 4 1 2 1 K K 4 1 2

i

π



 ln

2

K 22

2379

(6) (7)



2 (8)

In this step, the solution space has been reduced by choosing A = 1, θ = 0, and various reasonable values for 0 < 0.05. The parameter represents the modulation of the carrier wave. A small value is required since the expansion of the period matrix and other quantities related to the computation are carried out to O( ) [25]. Only single unstable modes ( N = 2) were considered for this work. In a future effort, the authors may examine the interactions between stable and unstable modes. The choice of θ and 0 is refined in the corrector step. These choices leave λ R and λ I , as the only free parameters governing the initial function selection. While the mathematical relationships presented in this manuscript refer exclusively to (λ R , λ I ), these values are chosen in the predictor step, and later explicitly computed in the corrector step. Based on these two distinct classifications, in the text, the authors refer to candidate values which are chosen by an investigator as (λRC , λIC ) and those that are determined as solutions to the eigenvalue problem as (λ R , λ I ). The two-dimensional space can be evaluated for periodic functions by direct numerical evaluation of ψ(x, 0) over an interval L. Each function evaluation of ψ(x, 0) is computationally expensive; however, each evaluation is independent of all others. This allows a large domain to be evaluated rapidly through a general purpose computing on a graphics processor unit (GPGPU) implementation. 2.2. GPGPU implementation Eq. (2) with (t = 0) was implemented in a CUDA kernel. The kernel computations were directly modeled after a single threaded Matlab validation code. All matrix operations in Eq. (2) were expanded by hand, and cuDoubleComplex variables were used to handle complex numbers. To utilize the processing capabilities of the GPU, a single thread was assigned to compute the function at a given coordinate (λRC , λIC ). Based on the hardware utilized, a maximum of 832 double precision cores can be harnessed concurrently to compute Eq. (2), each with a distinct value for (λRC , λIC ). The .cu file was compiled in Visual Studio 2010 targeting a 64-bit machine, with code generation options compute_20, sm_20. The use-fast-math option was deliberately excluded to provide the highest possible accuracy. All codes were executed on a Dell E5607 2.27 GHz machine with sufficient RAM running Windows 7 64-bit, Matlab R2013a 64-bit, and CUDA 4.2, on a Tesla k20c CUDA card. An independent version, used for validation, was coded directly in Matlab. Furthermore, the Matlab implementation was translated to a lower level c implementation through Matlab’s built in translator: codegen. The accelerations experienced by the GPU implementation increase with increasing domain size. For typical (λRC , λIC ) domain sizes of 256 × 256, the basic Matlab implementation required 133.0 seconds (baseline), the codegen implementation required 25.8 seconds (5.15×), and the CUDA kernel implementation required 0.0573 seconds (2, 321×). To provide a fair comparison and mitigate Matlab’s asynchronous execution of GPU kernels, the CUDA timing results were measured after memory was transferred back to the host using the gather() command. The periodicity of ψ(x, 0) was also determined as part of the computation, according to the metric presented in Eq. (11).

Fig. 2. Map of periodic Riemann theta functions as defined in Eq. (2) for ( A = 1, θ = 0, = 0.01, t = 0) generated by GPGPU computations. Light colored locations indicate periodic functions, while dark colored locations indicate aperiodic functions over the interval L. A point of interest is identified by an asterisk *.



ψ( L , 0) − ψ(0, 0)



ψ(0, 0)

C 0 =

(9)

C 1 = ψx (x, 0)|x=0 − ψx (x, 0)|x= L  C1 C 0 ≤ 0.01 Cf = N / A otherwise

(10) (11)

The resulting map of (λRC , λIC ) pairs, which form periodic ψ(x, 0) functions, is shown in Fig. 2. Parameter combinations resulting in periodic functions are displayed in white. Combinations that do not meet the periodicity criteria, are displayed in progressively darker colors. This unique method of exploring the thetafunction space has not been considered previously in the literature to the best of the authors’ knowledge. The authors acknowledge that the thresholds are ad hoc; however, valuable information related to the periodicity of the function is conveyed. The solutions are constructed and verified in the corrector procedure. The GPU map suggests that, given a rogue wave, similar rogue waves may exist for variations in (λ R , λ I ) along the string. Varying (λ R , λ I ) will necessitate revising the other parameters which govern the rogue wave as well. Studying the effect of these variations will be the topic of a future effort. In Fig. 2, the particular value (λRC = 1.2495, λIC = 1.6125) is identified as a solution of interest. This value was chosen only to serve as an example to illustrate the steps of the predictor corrector method. The procedure will be similar for any other (λRC , λIC ) pair along the string. 2.3. Corrector The parameters governing the periodic ψ(x, 0) function chosen above must be refined to yield a solution to the NLSE. The spectral eigenvalues are determined by solving the following spectral eigenvalue problem:



Ψx = Q (λ)Ψ,

Q =

−i λ

σ u∗

u iλ



(12)

The eigenvalue solution procedure follows closely the work of Osborne [25]. To summarize this procedure, the eigenvalue problem of Eq. (12) is recast as a Floquet problem by appropriately discretizing the complex wave train. The choice of λRC and λIC from

2380

C. Chabalko et al. / Physics Letters A 378 (2014) 2377–2381

Fig. 4. Solution to the NLSE for (λ R = 1.2415, λ I = 1.61108), = 0.006834, θ = 1.11439, and L = 4.44. This solution envelope has periodic temporal peaks, which reach a maximum amplitude of ≈ 4.2x the background. Fig. 3. Main spectrum eigenvalues with a detailed view of choice of (λ R , λ I ).

and θ for a particular

the GPU map implies periodic boundary conditions. A constant potential is then used to integrate the above eigenvalue problem to determine the spectral eigenfunction over a given interval.

Ψ (xn + x) = U (un , x)Ψ (xn )

ρ=

Then, along the lines of [25], the exponential of the trace of the vanishing matrix Q (λ) is given by U (un , x) as

U j j = cosh(k x) + (−1) j

iλ k

μ ( j = 1, 2)

U 21 = σ u ∗ μ

U 12 = u μ,

(13)

Here, k = σ |u | − λ and μ = sinh(k x)/k is constant inside an interval x. Floquet analysis [19] then dictates that the monodromy matrix is given by 2

2

T (x0 , λ) =

0

2

U (u j , λ)

(14)

Finally, the main spectrum eigenvalues, λk for k = 1, 2, . . . , 2N, are determined from the trace of the monodromy matrix as

2

T rT =

1 2

r(x, t )2 uˆ (x, t )2

( T 11 + T 22 ) = a R (x0 , λ) = ±1

(15)

Although many eigenvalues are determined, a pair of eigenvalues will be close to the single complex eigenvalue, which was used to construct the original ψ(x, 0), as shown in Fig. 3. From this pair of complex eigenvalues, the parameter is determined to be half the distance between the pair of eigenvalues and θ is determined as the angle between the pair and the horizontal, beginning from the line adjoining the pair, as shown in the expanded portion of Fig. 3. Based on experience, is recognized to be of the same order as 0 . 2.4. Verify A candidate solution, uˆ (x, t ), is determined by evaluating Eq. (2) with the new set of parameters. The candidate solution is then verified by numerical integration into the NLSE as

i uˆ t − uˆ xx + 2σ |uˆ |2 uˆ = r

(16)

(17)

All solutions presented in this work passed this test with ρ 1.00%. Second, a grid convergence test is performed. In this test, the quantity r(x, t )2 x t is verified to decrease as the grid size is refined in time and space for all solutions presented. Continuing with the procedure, after passing both tests, the candidate solution, uˆ (x, t ), is accepted as a solution u (x, t ) to the nonlinear Schrödinger equation, as shown in Fig. 4. The solution can be analyzed in terms of its maximum amplitude defined in [25] as

γ = A + 2λ I

j=M

1

Eighth order central finite differences are used to evaluate both the spatial and temporal derivatives. Due to discretization and finite precision approximations, the numerical solution may not identically satisfy the equation, resulting in a non-zero residual. Two criteria are used to verify the solution. First, the r(x, t )2 is compared to uˆ (x, t )2 , where  · 2 denotes the 2 norm, as

(18)

where A = 1 for all cases presented in this work. As of yet, unknown solutions with maximum amplitudes even slightly greater than the background can be critically important in systems where extreme waves are of interest. Solutions with large maximum amplitudes are of interest for other practical reasons. The GPU map presented in this effort allows identification of likely regions of solutions based on maximum amplitudes. 2.5. Example solutions The GPU map shown in Fig. 2 contains several branches of candidate periodic solutions. Several horizontal bands can be discerned in this map, the major features of which are symmetric about the x and y axes. Although the authors have computed many solutions by using this procedure, only three solutions are featured here for brevity. First, focusing on the point (λRC = 1.2495, λIC = 1.6125) with L = 4.44, the authors obtained the corrected parameter set of (λ R = 1.2415, λ I = 1.61108), = 0.006834, and θ = 1.11439. Based on this parameter set, the solution, shown in Fig. 4 has been generated. The maximum amplitude is approximately 4.2× the background modulation ( A = 1). The peaks are more compact in time than in space, dropping to the background level between adjacent crests as time increases. The next solution examined is for L = 7. The GPU map for L = 7 is not shown here in the interest of space; however, it is similar to

C. Chabalko et al. / Physics Letters A 378 (2014) 2377–2381

2381

for which the spectral eigenvalue problem yields the corrected parameters as (λ R = 0.46198, λ I = 0.79017), = 0.0025291, and θ = 0.660511086. This solution is illustrated in Fig. 6. The peaks reach a maximum amplitude of 2.58× the background and are less compact in time and space than the other solutions presented. 3. Concluding remarks

Fig. 5. Solution to the NLSE for (λ R = 0.6616, λ I = 1.23660), 0.87385, and L = 7.

= 0.00276, θ =

The authors have used a single predictor–corrector type procedure to determine new solutions to the nonlinear Schrödinger equation, based on periodic theta functions. A map of parameters, which likely lead to solutions, is used to characterize the solution space. Such a map can also assist in making decisions about solution features and determining what parameters are not likely to lead to rogue waves. The procedure presented here is not restricted to Riemann theta function solutions and can be extended for application to solutions based on other function spaces. The requirement of these solutions to be periodic makes parameter space exploration a computationally intensive task. As shown in this work, GPGPU computing is a natural fit for this type of computation, accelerating the computation by 400x over a single threaded high level language on commodity hardware. It is expected that other investigators can use this computing capability to explore this particular solution space and find as of yet unobserved nonlinear wave behaviors. Acknowledgement Support received for this work through NSF Grant No. CMMI1125285 is gratefully acknowledged. References

Fig. 6. Solution to the NLSE for (λ R = 0.46198, λ I = 0.79017), 0.660511086 and L = 100.44.

= 0.0025291, θ =

a compressed version of the parameter space map for L = 4.44. This solution exists on the second branch above the real axis in the corresponding GPU map. The starting point of (λRC = 0.6653, λIC = 1.2368) yields the corrected parameters (λ R = 0.6616, λ I = 1.23660), = 0.00276, and θ = 0.87385 through solution of the spectral eigenvalue problem. The corresponding solution is shown in Fig. 5. The maximum amplitude is 3.47× the background amplitude. This solution is characterized by localized peaks in time separated by flat regions of unit amplitude. The fluctuations again have compact support in the temporal domain. In order to provide solutions over a large spectral range, the final example shown has been computed for L = 100.44. The parameter space map for L = 100.44 is significantly compressed as compared to that shown in Fig. 2, and contains eight branches above the x-axis, rather than six, as in the case for L = 4.44. The initial point was chosen near (λRC = 0.4594, λIC = 0.8036)

[1] T.B. Benjamin, J. Feir, J. Fluid Mech. 27 (1967) 417–430. [2] C. Kharif, E. Pelinovsky, Eur. J. Mech. B, Fluids 22 (2003) 603–634. [3] E.R. Tracy, H.H. Chen, Phys. Rev. A 37 (1988) 815–839, http://dx.doi.org/10.1103 /PhysRevA.37.815. [4] N. Akhmediev, V. Korneev, Theor. Math. Phys. 69 (1986) 1089–1093. [5] A.O. Smirnov, Theor. Math. Phys. 173 (2012) 1403–1416. [6] A.O. Smirnov, Math. Notes 94 (2013) 897–907. [7] A. Smirnov, J. Math. Sci. 192 (2013) 117–125. [8] N.N. Akhmediev, A. Ankiewicz, Solitons: Nonlinear Pulses and Beams, vol. 4, Chapman & Hall, London, 1997. [9] N. Akhmediev, J. Soto-Crespo, A. Ankiewicz, et al., Phys. Rev. A 80 (2009). [10] N. Akhmediev, A. Ankiewicz, J. Soto-Crespo, Phys. Rev. E 80 (2009) 026601. [11] B. Kibler, J. Fatome, C. Finot, G. Millot, F. Dias, G. Genty, N. Akhmediev, J.M. Dudley, Nat. Phys. 6 (2010) 790–795. [12] D. Solli, C. Ropers, P. Koonath, B. Jalali, Nature 450 (2007) 1054–1057. [13] A. Chabchoub, N.P. Hoffmann, N. Akhmediev, Phys. Rev. Lett. 106 (2011) 204502, http://dx.doi.org/10.1103/PhysRevLett.106.204502. [14] A. Chabchoub, N. Hoffmann, M. Onorato, N. Akhmediev, Phys. Rev. X 2 (2012) 011015, http://dx.doi.org/10.1103/PhysRevX.2.011015. [15] K.B. Dysthe, Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci. 369 (1979) 105–114. [16] A. Osborne, J. Comput. Phys. 109 (1993) 93–107. [17] A. Shabat, V. Zakharov, Sov. Phys. JETP 34 (1972) 62–69. [18] A. Osborne, Mar. Struct. 14 (2001) 275–293. [19] A.H. Nayfeh, B. Balachandran, Applied Nonlinear Dynamics: Analytical, Computational and Experimental Methods, Wiley, 1995. [20] A. Dyachenko, V.E. Zakharov, JETP Lett. 81 (2005) 255–259. [21] N. Akhmediev, A. Ankiewicz, M. Taki, Phys. Lett. A 373 (2009) 675–678. [22] Y.-C. Ma, M.J. Ablowitz, Stud. Appl. Math. 65 (1981) 113–158. [23] O. Its, V. Kotliarov, Akad. Nauk Ukr. RSR Dopov. Ser. Fiz. Mat. Tekh. Nauki 1 (1976) 965–968. [24] I.M. Krichever, Russ. Math. Surv. 32 (1977) 185–213. [25] A. Osborne, Nonlinear Ocean Waves & the Inverse Scattering Transform, Int. Geophys. Ser., vol. 97, 2010, Access online via Elsevier.