GASE: a high performance solver for the Generalized Nonlinear Schrödinger equation based on heterogeneous computing

GASE: a high performance solver for the Generalized Nonlinear Schrödinger equation based on heterogeneous computing Nuno Miguel Azevedo Silva Mestrado em Física Departamento de Física e Astronomia 2013 Orientador Professor Ariel Guerreiro Todas as correções determinadas pelo júri, e só essas, foram efetuadas. O Presidente do Júri, Porto, ______/______/_________ Agradecimentos Chegando ao m de um ciclo de um ano de trabalho, não posso deixar de aproveitar a ocasião para agradecer a todas as pessoas que contribuíram para esta dissertação. Primeiro, agradeço ao meu orientador, Professor Ariel Guerreiro, por toda a paciência e ajuda prestada tanto durante os trabalhos da dissertação como na sua elaboração. Diria mesmo que, sem o seu auxílio, esta dissertação só estaria completa daqui a dois meses. Todo o trabalho não teria sido possível também sem o apoio inestimável do Paulo, que me ensinou os primeiros passos da computação em placas grácas. A ele, a todos os outros colegas do INESC e pessoas que contribuíram directamente para este trabalho, agradeço por todo o apoio prestado. Finalmente, agradeço à minha família, aos meus pais e irmão, à Rute e a todos os meus amigos, especialmente os que me acompanharam nos anos em que passei na FCUP. Sem todos eles, teria talvez acabado esta dissertação há dois meses atrás. Mas tenho a certeza que sem eles, não teria sido a experiência tão agradável e reconfortante que foi. i Abstract This dissertation presents a research work in the eld of computational physics, namely the development, testing and benchmark of a high performance solver of the Generalized Nonlinear Schrödinger equation that can address problems with high dimensionality and complex geometries, based on massive parallel computing using graphical processing units. The Generalized Nonlinear Schrödinger equation is an active topic of research and has attracted the attention of many researchers during the last decades. A major diculty in this eld is that the model is usually non-integrable and even perturbation methods are not valid in multidimensional and complex geometries. Instead, most of the research is done using numerical simulations to address this class of problems. However, in general, these have a high computational cost and can only be performed eciently in costly computer clusters or supercomputers. In recent years, the eld of computer sciences came up with a new computation concept called heterogeneous computing, that allows to use all the computing resources of a machine in an integrated way to do massive computing. This new computing paradigm is in the core of this dissertation as the enabling technology that supports the development of our solver. This dissertation begins with a general overview of the state-of-the-art in both Generalized Nonlinear Schrödinger equation and graphical processing units computing. The fundamental aspects of this equation are analyzed, presenting the most relevant analytical and numerical methods, which includes the overview of the Split-step Fourier method, the algorithm chosen for the development of the solver. The algorithm was implemented in the CUDA language, which runs on NVIDIA graphical processing units, and the extensive tests performed on the solver revealed that it outperforms the serial version of the solver. It is shown that the solver developed turns the graphical processing units into low-budget solutions for high performance computation, many times faster than state-of-the-art central processing units, and with performances that can compete with expensive supercomputers. Also, two physical problems described by the Generalized Nonlinear Schrödinger equation are considered in iii the last part of the dissertation and even though they have not been fully explored, given the limited duration of this dissertation project. The preliminary studies show some interesting results and illustrate the potential and versatility of the solver developed. Finally, some conclusions and future directions of research are discussed. This dissertation hopes to contribute to the eld of computational physics by showing how the new computing paradigms (such as heterogeneous computing) can be used to improve the performance of existing algorithms and methods and, through that, provide a tool of research to study more complex and demanding problems. iv Resumo Esta dissertação apresenta o desenvolvimento de um trabalho na área de física computacional, mais precisamente, o desenvolvimento, teste e análise de performance de um solver para a equação Não Linear de Schrödinger Generalizada, que seja capaz de resolver problemas de alta dimensionalidade e em geometrias complexas, baseado na computação paralela usando placas grácas. A equação Não Linear de Schrödinger Generalizada é um tópico bastante activo e tem atraído a atenção de muitos investigadores nas últimas décadas. Uma das maiores diculdades na investigação destes sistemas é que geralmente o modelo é não-integrável e mesmo os métodos perturbativos são incapazes de oferecer soluções para estes problemas. Assim, a grande parte da investigação passa pela simulação numérica desta classe de problemas. No entanto, estes têm em geral um elevado custo computacional e só em clusters de computadores ou supercomputadores a sua simulação é eciente. Nos últimos anos, surgiu um novo conceito na ciência de computadores denominado computação heterogénea, que permite a utilização integrada de todos os recursos de um computador para o cálculo de grandes tarefas numéricas. Este paradigma computacional constituí o núcleo desta dissertação ao ser a tecnologia que permite o desenvolvimento do solver . Esta dissertação começa com uma síntese geral do estado da arte da investigação tanto do caso da equação Não Linear de Schrödinger Generalizada como da computação em placas grácas. Os aspectos fundamentais da equação são analisados e os métodos analíticos e numéricos mais importantes são apresentados, com especial atenção ao Split-Step Fourier Method , o algoritmo que escolhemos para utilizar no solver . O algoritmo foi implementado na linguagem CUDA , que corre em placas grácas da NVIDIA , e uma grande variedade de testes foram executados ao solver, revelando performances muito acima das obtidas para versões tradicionais com base no processamento em série. É mostrado também que o solver desenvolvido torna um computador com uma simples placa gráca numa solução de baixo custo para obtenção de performances elevadas, muito mais rápidas que os habituais processadores em série e com resultados capazes de competir com supercomputadores muito mais caros. Além desta análise, dois v problemas descritos pela equação Não Linear de Schrödinger Generalizada são considerados na última parte da dissertação, mesmo não tendo sido completamente explorados devido à duração limitada deste projeto de dissertação. Os estudos preliminares mostram alguns resultados interessantes e acima de tudo ilustram o potencial e versatilidade do solver desenvolvido. Deste trabalho espera-se sair uma contribuição para o ramo da física computacional ao demostrar como novos paradigmas computacionais (como a computação heterogénea) podem ser utilizados para melhorar a performance dos algoritmos e métodos existentes e, assim, oferecer uma ferramenta computacional capaz de resolver problemas computacionalmente custosos. vi Contents 1 Introduction 1 1.1 Asolitonstory ................................ 2 1.2 1+1=3 and the Nonlinear Schrödinger Equation . . . . . . . . . . . . . . 4 1.3 Opticalsolitons................................ 6 1.4 Analytical, variational and numerical methods . . . . . . . . . . . . . . . 10 1.5 GPU computing and prospects for the NLSE . . . . . . . . . . . . . . . . 11 1.6 GASE - GPU Accelerated Soliton Explorer . . . . . . . . . . . . . . . . . 15 1.7 Outline and structure of the dissertation . . . . . . . . . . . . . . . . . . 15 2 Nonlinear Schrödinger equation in a nutshell 17 2.1 NLSE and analytical solutions . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 GNLSE .................................... 20 2.3 Noether's theorem and conservation laws in the GNLSE . . . . . . . . . . 21 2.4 Variationalmethods ............................. 23 2.4.1 Trial functions and solutions . . . . . . . . . . . . . . . . . . . . . 23 2.4.2 Perturbation methods . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.3 Eective particle approach . . . . . . . . . . . . . . . . . . . . . . 26 2.5 Numericalmethods.............................. 27 2.5.1 Explicit and implicit nite dierences methods . . . . . . . . . . . 29 2.5.2 Pseudo-spectral methods and the SSFM . . . . . . . . . . . . . . 31 2.5.3 Boundary conditions for the SSFM . . . . . . . . . . . . . . . . . 34 2.6 Concludingremarks.............................. 36 3 Implementation of the GPU-based GNLSE solver 37 3.1 Simple problem, high computational time . . . . . . . . . . . . . . . . . . 37 3.2 Howtoplowaeld? ............................. 40 3.3 Gaming vs Scientic precision . . . . . . . . . . . . . . . . . . . . . . . . 42 vii List of Tables 3.1 Comparison between dierent top-of-the-line GPUs of the NVIDIA consumer line Geforce. A model from the Fermi line of a professional computing dedicated GPUs is also presented. An increasing computational power is notorious over the years, as well the increase of chip memory and memory bandwidth. It is well patented that evolution of GPUs will reach another level in the next few years, with the new high performance Geforce Titan setting the pace. . . . . . . . . . . . . . . . . . . . . . . . 43 4.1 Error analysis for the simulations with xed integration step h= 0.01 and variable number of points, which introduces a variable discretization ∆x . ...................................... 57 4.2 Error analysis for the simulations with xed number of points N= 28 and variable integration step h . ....................... 57 4.3 Accuracy and performance comparison between GASE solver, based on SSFM, and the CN solver. It can be easily seen that GASE outperforms in every aspect the CN method. . . . . . . . . . . . . . . . . . . . . . . . 58 4.4 Specications of both the two GPUs and the CPU used during the benchmarks...................................... 59 4.5 A collection of results for simulation times and speedup of the solver for the (1+1)-d NLSE using double precision, for both GASE (running in the desktop GPU) and the CPU-based version of the solver. . . . . . . . . . 61 4.6 A collection of results for simulation times and speedup of the solver for the (2+1)-d GNLSE for a cubic-quintic media, using double precision, for both GASE (running in the desktop GPU) and the CPU-based version of thesolver.................................... 63 xv 4.7 A collection of results for simulation times and speedup of the solver for the (2+1+1)-d GNLSE for a cubic-quintic media, using double precision, for both GASE (running in the desktop GPU) and the CPU-based version ofthesolver................................... 65 xvi List of abbreviations GNLSE - Generalized Nonlinear Schrödinger Equation GPU - Graphical Processing Unit GPGPU - General Purpose GPU applications CPU - Central Processing Unit KdV - Korteweg-de-Vries FPU - Fermi-Pasta-Ulam CW - Continuous Wave BEC - Bose-Einstein Condensate IST - Inverse Scattering Transform FD - Finite Dierences FFT - Finite Fourier Transform PS - Pseudo-Spectral SSFM - Split-Step Fourier Method xvii 1 Introduction The central problem of this dissertation is the development of a high performance solver of the Generalized Nonlinear Schrödinger Equation (GNLSE) based on heterogeneous programming using graphical processing units (GPU), that is able to address physical systems with high dimensionality (more than one spatial dimension) and complex geometries or structures with reasonable simulation times and that can run using low cost desktop computers. This constitutes a paradigm in computational physics since until few years ago this type of problems could only be addressed using costly state-of-the-art supercomputers and computer clusters, generally inaccessible to most scientists. GPU computing is bringing a revolution into scientic computing by allowing to do supercomputing by using several hundreds of processing units inside a desktop computers as a massive cluster. However many of the computational models and simulation codes previously developed cannot be straightforwardly adapted to heterogeneous computing, given its distinct computing architecture and the specic programming tools required. Although the focus of this dissertation is on the development of the GNLSE solver, both as a proof of concept and as a simulation tool for future research, in the following chapters we will also present several case studies where we do some earlier exploration of soliton dynamics, mainly as a test and illustration of the potential of this code. These examples were not fully explored in terms of a complete scientic analysis since that goes beyond the scope of this dissertation and would require a longer research time. The existence and propagation of solitons in nonlinear media attracted a substantial research interest for the past 50 years. While thousands of papers were published in this eld and particularly in the study of the Nonlinear Schrödinger Equation (NLSE), new and more complex systems provide substance for prospective explorations. As current investigations focus on multidimensional solitons and spatial distribution of nonlinearity, described by non-integrable models, numerical simulations become mandatory. In a normal computer, simulation times for such systems are usually prohibitive and researchers that do not have access to a supercomputer are either limited in the research to smaller and simple systems. In this context, the development of new and high performance 1 CHAPTER 1. INTRODUCTION computational tools for the study of the propagation of solitons is a subject of great importance for state-of-the-art problems. Heterogeneous computing is one of the newest and fascinating trends in modern physics. It consists on the use of the central processing unit (CPU) in addition to other specialized hardware, usually the graphical processing units to get faster numerical simulations. Particularly, the use of the GPU for general purpose applications (GPGPU) created a buzz in recent years, as researchers from many areas attained overwhelming performances - up to 100 times faster than state-of-the-art CPUs. This potential, that is at the same level of the 1999-2000 timeframe best supercomputers, suggests that GPU based solvers of the NLSE could be the solution for simulating the cutting-edge demanding problems using only inexpensive personal machines, a hypothesis that is the motivation behind this dissertation. 1.1 A soliton story The story of solitons is prolic in coincidences and fortuitous events. It all started in 1834 with a curious scottish engineer, John Scott Russell, hired to investigate how to improve the eciency of boat designs at the Union Canal, near Edinburgh. In a fortuitous accident, a rope pulling a boat broke and Russell observed the formation of a wave that he described accurately in his report [69] as a large solitary elevation, a rounded, smooth and well-dened heap of water, which continued its course along the channel apparently without change of form or diminution of speed . Most probably, this was not the rst time that a solitary wave was observed, but Russell was the rst to report it. Believing that the discovery was important he did extensive experiments in a scale model constructed at his backyard. In 1895, Dutch physicists Diederick Korteweg and Gustav de Vries derived an equation [48] to match the observations reported by Russell. This partial dierential equation, now called Korteweg-de Vries equation (KdV), contained both linear and nonlinear terms and although they were unable at that time to produce general solutions, they found a solitary-wave solution that resembled Russell's wave. Strangely, as it happened to Russell, their work fell into obscurity and was overlooked by mathematicians, physicists and engineers for more than 50 years. The story continues in the early 1950s, at Los Alamos Scientic Laboratory when Enrico Fermi, John Pasta and Stanislaw Ulam used one of the earliest digital computers - the MANIAC (MAthematical Numerical Integrator And Computer) - to investigate the simple nonlinear system of a one dimensional chain of masses connected by springs with 2 1.1. A SOLITON STORY Figure 1.1: a) On 1995 during a conference on nonlinear waves at Heriot-Watt University, the attending scientists recreated the rst reported observation of a solitary wave, as part of a ceremony to honor Russell and name the new aqueduct with his name. Image from [1] b) Figure taken from the original report of the FPU problem [22], showing a simulation performed in the MANIAC. Fermi, Pasta and Ulam initialized the system with all energy in the lowest normal mode and observed the evolution. The energy is transferred into several modes and, after some time, the system returns to a state similar to the initial condition, with energy back to the rst mode, unlike the expected thermalization. linear and small nonlinear restoring forces [22]. Exciting one normal mode of the linear system, they believed that the nonlinearity term would excite dierent modes and then at some point in time the system would thermalize, i. e., the energy would be equally distributed among all the possible normal modes. Nevertheless, the results were unexpected: while it is true that after some periods the energy was shared between several modes, the prolongation of the simulation revealed a near return to the initial mode, as 97% of the energy focused again in the initial mode. This unexplained recurrence, later known as Fermi-Pasta-Ulam problem (FPU), did not convinced everyone, as some thought the system did not run enough time [8]. One of the rst attempts to solve this puzzle was a phenomenological explanation suggested in 1965 by Zabusky and Kruskal [86]. Taking the FPU in the continuum limit, they found that the system was governed by a KdV equation, an equation analytically intractable at that time. Solving the equation using numerical simulations, they observed a breakdown of the initial periodic condition into a train of solitary waves. 3 CHAPTER 1. INTRODUCTION Regardless of the initial condition tested, the recurrence of the initial wave suggested by the results of FPU was not observed during the simulations. Instead, they detected that the solitary waves started to move and collide. They also observed a quasi-recurrence to an intermediate state after a near chaotic behavior. It is a common mistake to attribute the solution of the FPU problem to this study of Zabusky and Kruskal [46]. Despite obtaining a recurrence, their belief that the discovery provided a phenomenological explanation was poorly grounded, as they did not obtain any return to initial mode. A more concise explanation was only reached recently (1997) by Casetti et al. [13] with the discovery of two regimes for the dynamics of the FPU that depends on the energy per oscillator and on the total number of masses. For low oscillator energy and number of masses, dynamics are regular to weakly chaotic, like the FPU results. For higher oscillator energies and number of masses, dynamics are completely chaotic. Indeed, we should be thankful that Fermi and his co-workers did not simulate bigger systems nor used stronger nonlinearities, because if they did, they would have found an equipartition of the energy, and then maybe much of the understanding of nonlinearities that arose from their studies might not have existed. In spite of not succeeding in explaining the FPU problem properly, Zabusky and Kruskal's publication still became very famous. They observed the survival of solitary waves after collisions - a behavior typical of a particle - and led them to coin one of the most successful terms in nonlinear science: the soliton. 1.2 1+1=3 and the Nonlinear Schrödinger Equation When one starts to study physics it is common to think that nonlinear is synonym of anomalous, pertaining to something that diverges from well behaved physics. However, nonlinear systems are far more common in the real world than the linear ones. The study of nonlinear systems is the subject of nonlinear science, and the main idea is that the whole is more than a sum of its parts , or in physicist language, the superposition principle is not valid. Another feature that boosted the development of the nonlinear science was its universality, not only in terms of being present in almost every phenomenon of the Universe, but also because it can be found in many elds of science. Indeed, models explaining phenomena like chaos or coherent structures can be used to describe a wide panoply of problems, not only in theoretical physics but also in mathematics, biology, neurosciences, sociology and more [71]. 4 1.2. 1+1=3 AND THE NONLINEAR SCHRÖDINGER EQUATION The Nonlinear Schrödinger Equation is one of this models. The NLSE was introduced rst in 1964 by Chiao et al. [14] to describe the propagation and self-trapping of continuous wave light beams (CW) incident in a nonlinear Kerr media, both in one and two spatial dimensions. Soon, the rst hints of the universality started to appear, when Hasewaga and Tappert [30] suggested a NLSE equation to describe the propagation of pulsed beams in optical bers. Nowadays, it is normally accepted that NLSE i∂ψ ∂z +1 2∇2 ⊥ψ+|ψ|2ψ= 0 (1.1) and the more global GNLSE i∂ψ ∂z +1 2∇2 ⊥ψ+F(|ψ|2)ψ= 0 (1.2) are universal equations that describe the evolution of wave envelope in a weakly nonlinear medium. This equation arises in many and distinct areas, the most common being: • Nonlinear optics: to describe the propagation of CW and pulsed beams in nonlinear medium [14, 46, 30]; • Bose Einstein Condensates (BEC): to describe the mean-eld dynamics of the BEC [43, 27] where it is commonly known as Gross Pitaevski equation; • Fluid dynamics:to describe for example the instability of Poiseuille ow [76], deep water waves [85] and Couette-Taylor ow [19]. It is commonly known as Complex Ginzburg-Landau equation; • Plasma physics: to describe Langmuir waves [6]; • Protein chemistry: to model the vibrations of molecular chains [16]. This variety of applications reinforced the interest and the search for solutions for the NLSE. In the seminal paper, Chiao and his co-workers said at some point that NLSE appears to have no simple analytical solution [14] leading them to search for numerical solutions. Although the solution had a bell shape similar to a solitary wave, the relationship was not explicitly noticed. Meanwhile, the paper of Zabusky and Kruskal on solitons [86] led to the development of the Inverse Scattering Transform (IST) method in 1967 [26]. This elegant mathematical tool is an analog of the Fourier transform for nonlinear systems, such as the initial value 5 CHAPTER 1. INTRODUCTION Figure 1.4: a) GPUs have long surpassed peak performances of CPUs, which triggered the development of GPU computing. b) Also, they attain this peak performances without prohibitive power consumptions, characteristic of frequency scaling of the CPUs, using the parallel computing paradigm. 12 1.5. GPU COMPUTING AND PROSPECTS FOR THE NLSE ical scientic computing is no longer a nuclear objective of computer developers, it can still much benet from the continuous waves of innovation and improvements constantly occurring in this technology. During the past decade, the frantic demand for faster processors by software industries made the computer engineers to come up with the parallel programming paradigm for common purposes. This included the development of multi-core CPUs with ever increasing number of cores and the development of the necessary software to allow them to work in parallel. This helped to achieve computer performances as never before. The quest for parallel computation is not limited to CPUs. GPUs also experienced a similar evolution. Modern GPUs contain several hundreds of cores and they exploit the highly parallelizable task of calculating the value of a pixel, thus speeding up video games and image processing. The former is so important for the gaming industry that GPUs have long surpassed multi-core CPUs in both number of cores and computer performance, reaching incredible speeds of the order of a Teraop/s, while the best CPUs are limited up to 50 GFLOP/s. GPUs have yet another advantage, as they cost the same or less than a CPU. This fact can get even striking, as the performance of a mid-range GPU is the same of the best supercomputers of year 2000, that cost 110$ million, for example ASCI White. Excited by this technology, science world started to think about the use of GPU for general purposes (GPGPU) such as scientic computations, but the rst approaches were challenging, as no easy and versatile programming framework was available. Recognizing the problem, NVIDIA made an eort to make this potential available to the industry and scientic community by developing a new programming framework for NVIDIA GPUs called CUDA. With CUDA and more recently with OPENCL - a platform-independent framework - the modern researcher can nowadays move his computational codes to the GPU of his personal computer and obtain speedups 2 worthy of a modern supercomputer. However life is not perfect yet! There are some obstacles for the average physicist to become a GPGPU user, as learning the architecture of the GPU and new programming paradigms can be time consuming and a truly jigsaw puzzle. The rst use of GPU computing for solving dierential equations was probably by Mark Harris [29] and ever since GPU computations proliferated in many areas of physics. A few examples of scientic computation using GPU include: • Fluid dynamics: the power of GPU is exploited to simulate large systems showing 2 Speedup is a measure of the relative performance of a code developed in parallel when comparing to serial computations, dened usually as Speedup = serial computing execution time parallel computing execution time. 13 CHAPTER 1. INTRODUCTION turbulence obtaining simulations 22 times faster than CPUs [44]; • Statistical physics: Multidimensional Ising model simulations were done with speedups of 8 [47, 66], while Brownian dynamics and reaction-diusion systems were simulated with speedups of 8 and 55 [65, 82]; • Electromagnetic waves: Maxwell equation was simulated with speedup of 50 and 60 in two [7] and three [57] dimensional systems, respectively. NVIDIA also provides an extensive report [15] that shows the utilization of graphical cards for GPU-ready commercial software. From bioinformatics - sequence mapping software up to 100 times faster - to computational nances - nancial analytic software 500 times faster - the applications are numberless. It was in 2009-2012 timeframe that Ron Caplan developed the rst approach of a NLSE solver using GPU computations during his PhD at San Diego University [11]. The code developed, that later became a package for MATLAB called NLSEmagic, served his purpose obtaining simulations up to 20 times faster than a CPU based script [12]. This code was a conceptual breakthrough by demonstrating the potential power of GPUs in solving the NLSE. However, it had many important limitations: • First, is based on a Runge-Kutta scheme, which is an explicit FD method and then conditionally stable and slower than PS methods for most of the computations; • NLSEmagic does not admit spatial distribution of the nonlinearities, which is very important for current research; • Finally, NLSEmagic package is developed in MATLAB and although he used CMEX - a MATLAB interface for developing part of the code in C - MATLAB is still a scripting language, hence with low performance when compared with C or C++. Based on the experience acquired during the preparation of this master's dissertation, it is my opinion that even if Caplan has achieved a large speedup, he was strongly limited by using a MATLAB script, which can be easily outperformed by a CPU C++ based code. In conclusion, GPU computing has an enormous potential for researchers and engineers. In particular, there is still space for improvements regarding high performance GNLSE solvers and integrators. Regardless of the NLSE simulation in GPU being already developed recently, it is our belief that there is still room for progress and higher performance, because it neither allows to solve the GNLSE nor uses high-performance C++ language, 14 1.6. GASE - GPU ACCELERATED SOLITON EXPLORER nor it is based on spectral methods, which are usually more stable and have higher performance. 1.6 GASE - GPU Accelerated Soliton Explorer The solver of the GNLSE developed during this dissertation project was named GASE. It consists in an executable le compiled from a CUDA C++ code compiled using Microsoft Visual Studio, which computes the numerical solutions of the GNLSE, as well as a series of complementary scripts in Python which do the analysis of the data and produce the graphical outputs. The solver can run in a normal computer having a NVIDIA GPU installed and enabled to use CUDA. This hardware component usually costs about few hundreds euros and is the element responsible for doing most of the massive computation in GASE. The code GASE is capable of simulating physical problems with 1, 2 and 3 spatial dimensions in simulation boxes with a number of sampling points up to 223 , although this value is only limited by the hardware and not by the code itself. GASE can simulate systems with any type of nonlinearity, including cubic, quintic and logarithmic, as well as nonlinearities that have a spatial dependence. Also, a recent upgrade of GASE allows to simulate a system of two coupled GNLSE (this can also be extended for more than two GNLSE). The code is also prepared to simulate problems with periodic, reective and absorbing boundary conditions. In short, this code has a high performance when compared with other sequential algorithms and is designed to be able to address a wide class of problems. 1.7 Outline and structure of the dissertation This dissertation addresses aspects of two immense topics: GPU computing and the NLSE. Its main output is the development of a solver of the GNLSE based on CUDA in C++ framework, that uses GPU computing and is capable of addressing problems with high dimensionality and spatial distribution of nonlinearities, such as optical lattices. We hope that this output can give a contribution to other researchers by providing them with a tool to investigate computationally modern problems in nonlinear science. This is the result of one year of work whereas becoming an expert in both GPU computing and the NLSE requires years of dedication. Not surprisingly, there is still space for further improvement of the simulation code, as well as to explore its full scientic potential. 15 CHAPTER 1. INTRODUCTION It has been a long way to reach this point, a road covered with many hours at the computer, testing dierent algorithms, making few mistakes and learning from them. In the words of Edison, I have many results, I know many things that do not work. Indeed, the hardest tasks during the last year were to learn C++ and CUDA, to become familiar with GPU architecture and with the concepts needed in the development of the code, and to overcome all the diculties that come with unexplored territory, without having any work as reference as it was one of the rst GPU based codes developed at the department. These tasks consumed more than eight to nine full months but however they have no place in this dissertation, as it only reports what worked well. In an analogy, this dissertation is like a building: the outcome can be analyzed and we can tell how we built it, but the hard work and the needed strength can be wrongly overlooked. The dissertation is structured as follows. In this rst chapter, a general overview of the subject was given, in theory of solitons, NLSE and GPU computing. A small motivation and the framework was also discussed. In Chapter 2 a brief synopsis of NLSE and solitons is presented, discussing succinctly the mathematical formulation and examples of variational methods and eective particle approach. The numerical methods, both nite-dierences and a pseudo spectral method called Split Step Fourier Method (SSFM) are introduced, followed by a discussion of the boundary conditions. The code implementation is discussed in chapter 3, and a comparison of the performance of GPU-based versus CPU-based simulations is described in chapter 4. In chapter 5 and 6 we present two case studies as proof of concept of the developed tools. In chapter 5 we analyze a one dimensional chain of spatial solitons, predicting numerically and showing computationally the possibility of having phononlike oscillations. Chapter 6 is devoted to the problem of soliton collision in (2+1)-d system, investigating both the in-phase and out-of-phase soliton collision. Finally, future perspectives and an outline of the main conclusions are provided in chapter 7. 16 2 Nonlinear Schrödinger equation in a nutshell This chapter is devoted to briey review some of the principal aspects of the NLSE and its solution, focusing on the main analytical results and introducing the most relevant numerical methods. Given the tremendous work done over the years and extensive literature in this topic, we restrain this review to the aspects which are most relevant to future chapters and more specically to the development of the code GASE. In particular, we mainly discuss the solutions of the (1+1)-d NLSE, which is (to our knowledge and so far) the only case with a generic method of obtaining exact analytical solutions via IST method. However, we also describe ways of obtaining approximate soliton-like solutions for the GNLSE in cases with higher dimensions. Finally, we discuss the most notorious successful numerical methods used to solve computationally the GNLSE, namely the Finite Dierences (FD) methods and Pseudo-spectral (PS) methods. In particular, the review of the Split-step Fourier method (SSFM) establishes the ground base for the two following chapters, since it corresponds to the numerical method used in the development of our solver of the GNLSE. In the following chapter we describe how this numerical method was adapted and implemented to work on GPUs and make use of its tremendous computing power. 2.1 NLSE and analytical solutions Traditionally, the NLSE refers to the Schrodinger equation where an extra term was added, corresponding to a cubic nonlinearity: i∂ψ ∂z +1 2∇2 ⊥ψ+s|ψ|2ψ= 0. (2.1) This model is widely spread in nonlinear science and describes the evolution of a dimensionless amplitude eld ψ in a dispersive and weakly nonlinear medium. In this formulation, the coordinate z is usually associated to the longitudinal direction, along 17 CHAPTER 2. NONLINEAR SCHRÖDINGER EQUATION IN A NUTSHELL which the eld propagates, while the ∇2 ⊥ is the Laplacian in the transverse directions. Also the number s=±1 refers to the sign of the nonlinearity. The (1+1)-d NLSE corresponds to a simplied version of equation (2.1) and is given by i∂ψ ∂z +1 2 ∂2ψ ∂x2+s|ψ|2ψ= 0. (2.2) It has soliton solutions for both values of s , which can be calculated using IST ([14]). In nonlinear optics equation 2.2 has two major applications. On one hand, when x refers to a spatial coordinate then, the NLSE describes the connement of a CW light beam in a Kerr media [14]. On the other hand, when x refers to a temporal coordinate (sometimes x is replaced by τ in this equation to make the temporal character more evident), then, NLSE describes the propagation of a pulsed beam in an optical ber [30]. From now on, and unless noted otherwise, only spatial solitons are considered. For s= 1 , the media is also called self-focusing and the supported solutions are called bright solitons. The one dimensional bright soliton solution of equation (2.2), centered at constant x= ¯x0 , is given by [46] ψ(x, z)=2ν sech [2ν(x−¯x0)] exp 2iν2z, (2.3) where ν is the amplitude of the soliton. It can be proven that the NLSE is invariant under the Galilean transformation [77] x7→ x0=x−µz z7→ z0=z (2.4) ψ(x, z)7→ ψ0(x0, z0) = ψ0(x−µz, z) exp iµx −iµ2z/2 allowing us to consider a more general solution of a moving soliton, ψ(x, z) = 2ν sech [2ν(x−¯x0−µz)] exp {iµ (x−¯x0) + iδ(z)} δ(z) = (2ν2−µ2/2)z+δ0 (2.5) where ν is the amplitude, ¯x0 is the initial position of the centroid of the eld distribution ψ , µ is the transverse velocity and δ0 is the initial phase of the soliton. Bright solitons can exist for all values of ν and µ , constituting a two-parameter family of solutions. For s=−1 , the media is called self-defocusing and the supported solitons are called 18 2.1. NLSE AND ANALYTICAL SOLUTIONS Figure 2.1: Depending on the positive or negative sign of the nonlinearity, solitons can be either bright - gure a) - or dark - gure b) - respectively. Both were represented in arbitrary units. dark solitons. Using the IST and the boundary condition |ψ|=ψ0 as x→ ∞ , these dark solitons solutions correspond to a localized intensity reduction in an otherwise constant CW background. They can be expressed analytically [46] as ψ(x, z) = ψ0{Btanh [ψ0B(x−Aψ0z)] + iA}exp −iψ2 0z (2.6) where the two parameters A and B obey the relation A2+B2= 1 . Even though it is possible to interpret the parameter A as being related with the velocity of the dark soliton, the similarities end here. Dark and bright solitons have very distinct properties since they are not the dual of each other. Such properties will not be discussed here since they fall out of the scope of this dissertation. From now on, the discussion is focused on the case of bright solitons, that shall be designated simply as solitons. The case of (1+1)-d solitons is well studied in the literature, much due to the development of IST. However, the dimensionality of the physical system has a key role on the nature of the solutions of NLSE. Indeed, it is proven that for the (D+1)-d NLSE, with a Kerr-type nonlinearity, a generic localized envelope-like solution collapses into a singularity for D= 2 [46, 87]. This behavior is indicative of the diculty of nding envelope and soliton-like solutions of the NLSE in systems with high dimensionality. In fact, stable (2+1)-d solitons are only possible considering higher and saturable nonlinearities, described by the GNLSE. 19 CHAPTER 2. NONLINEAR SCHRÖDINGER EQUATION IN A NUTSHELL 2.2 GNLSE The GNLSE is generalization of the NLSE obtained by replacing the cubic nonlinearity with a generic nonlinear term, say i∂ψ ∂z +1 2∇2 ⊥ψ+F(|ψ|2)ψ= 0, (2.7) where F(|ψ|2) is a real valued function describing the nonlinearity. The interest of the GNLSE is that for certain nonlinearities it allows for (D+1)-d soliton-like solutions with D > 1 , and specially for (2+1)-d. Various types of non-Kerr law nonlinearities had been studied, these include the following: • Parabolic law: F(|ψ|2) = |ψ|2+s|ψ|4 For many optical materials and media, the refractive index begins to deviate from the Kerr type for large intensities of ψ . For example, a polydiacetene para-toluene sulfonate (PTS) crystal has a parabolic law dependence for the refraction index with s < 0 [18]. This situation, usually called cubic-quintic media, relies on the competition between the two nonlinearities to stabilize (2+1)-d solitons. In fact, at low intensities the self-focusing dominates the system, but for high intensities the beam collapse is avoided by the self-defocusing eect. • Saturating law: F(|ψ|2) = λ1−1 (1+|ψ|2/|Ψsat|2) Simple two-level atomic systems [9] or photorefractive materials [37] displays a type of nonlinearity which saturates for eld intensities above |ψsat|2 . On the other hand, when the nonlinear term in the (2+1)-d GNLSE depends explicitly on the spatial coordinates (described by a formal dependence of F on the transverse and longitudinal coordinates, say F(|ψ|2;r⊥, z) ) it describes the so called optical lattices, which can also support the formation and propagation of solitons. Typically optical lattices can be classied into two main classes: • Linear lattices: Considering that F(|ψ|2;r, z) = F1(|ψ|2) + V(r, z) we obtain an equation that governs the solitons in a linear lattice i∂ψ ∂z +1 2∇2 ⊥ψ+F(|ψ|2)ψ+V(r, z)ψ= 0 (2.8) described by the potential V(r, z) . In BECs, V(r, z) is the trapping potential and this equation is usually called Gross Pitaevski equation. 20 2.3. NOETHER'S THEOREM AND CONSERVATION LAWS IN THE GNLSE • Nonlinear lattices: Nonlinear lattices are spacial distributions of the nonlinearities. It is not possible to present a general formula for this case but for example, one can consider a cubic-quintic nonlinearity modulated in space by the function R(r, z) , being the problem governed by i∂ψ ∂z +1 2∇2 ⊥ψ+R(r, z)|ψ|2− |ψ|4ψ= 0. (2.9) Unfortunately, it is only for very specic cases that GNLSE constitutes an integrable model with analytical soliton solutions. In a rigorous sense, solutions of nonintegrable systems are not solitons, but it is common to use the term because most of the solutions normally tested are soliton-shaped waves. As nonintegrable GNLSEs are important models in nonlinear physics, it is necessary to develop methods capable of analyzing the properties of solitons in such systems. In the following sections some of these methods are briey reviewed. 2.3 Noether's theorem and conservation laws in the GNLSE Unlike the (1+1)-d NLSE, the GNLSE does not possess an innite number of conserved quantities [14]. However, the existence of some conserved quantities is of great importance for the analysis of the GNLSE, as it can provide clues about soliton behavior. An investigation of the conservation laws is based on the structure of the Lagrangian associated with GNLSE described by equation (2.7). This Lagrangian density is dened as L=i 2(ψ∗ψz−ψψ∗ z)−1 2∇⊥ψ∇⊥ψ∗+G|ψ|2, (2.10) with G(λ) = ˆλ 0 F(λ)dλ, (2.11) where we introduced the notation ∂ψ ∂z =ψz and x refers to the transverse coordinates [77]. 21 CHAPTER 2. NONLINEAR SCHRÖDINGER EQUATION IN A NUTSHELL mon are the Crank-Nicholson implicit scheme and pseudo-spectral methods which are both unconditionally stable [25]. The presence of the nonlinear term changes drastically the situation and the otherwise known methods are no longer unconditionally stable, becoming important to choose carefully the integration step. Usually, integration schemes for nonlinear problems rely on the split of the integration step into a linear and a nonlinear sub-steps. The splitting procedure is a well known mathematical method and the idea is to decompose a complex model into a sequence of simple sub-problems. The method has an associated error that can be theoretically estimated. One of the most popular splitting methods is the Strang-Splitting algorithm which constitutes a second order splitting [2]. To explain how the splitting procedure can be applied to the GNLSE we consider that it can be written as: ∂ψ ∂z =ˆ D+ˆ Nψ, (2.42) where ˆ D=i 2∇2 ⊥ is a linear operator containing the linear terms of the GNLSE and ˆ N=iF(|ψ|2) is its nonlinear counterpart. When applied to this equation, StrangSplitting algorithm is as follows: ∂ψ ∂z =ˆ Dψ, with z∈[z, z +h/2], ψ(z) = ψ(z) (2.43) ∂ψNL ∂z =ˆ NψNL, with z∈[z, z +h], ψNL(z) = ψ(z+h/2) (2.44) ∂ψL ∂z =ˆ DψL, with z∈[z+h/2, z +h], ψL(z+h/2) = ψNL(z+h) (2.45) where h is the integration step. It can be proven that this algorithm is second order accurate. The nonlinear problem is then divided into two simple problems. The linear sub-step requires the integration of the diusion equation which can be done as previously discussed. The nonlinear sub-step can be integrated numerically using the Euler method, the Runge-Kutta or other methods. The splitting method is commonly referred in beam propagation studies as the splitstep. The idea is that the linear and nonlinear parts of the dynamics can be treated separately considering small integration steps. Depending on the strategy to solve the linear sub-problem of the algorithm, beam propagation methods can be either nite dierence or pseudo-spectral. 28 2.5. NUMERICAL METHODS Figure 2.3: Visual scheme describing the split-step algorithm for evolving an initial eld ψ(z) , described in equations (2.43-2.45). 2.5.1 Explicit and implicit nite dierences methods Finite dierence methods can be grouped into two broad categories: explicit or implicit schemes. During this section we briey present an example of each. The fourth-order Runge-Kutta (RK4) scheme was used by Ron Caplan in the rst CUDA solver of the NLSE [11]. The RK4 is an explicit method that can be used also for the GNLSE. Writing the GNLSE as ∂ψ ∂z =f(ψ) = i1 2∇2 ⊥ψ+F(|ψ|2)ψ (2.46) the RK4 scheme is dened by [12] k1=f(ψ(z)) (2.47) k2=f(ψ(z) + h 2k) (2.48) k3=f(ψ(z) + h 2k2) (2.49) k4=f(ψ(z) + h 2k3) (2.50) ψ(z+h) = ψ(z) + h 6(k1+ 2k2+ 2k3+k4). (2.51) To convert the GNLSE into an algebraic equation, the domain of ψ is replaced by a discrete set of points that lay on a regular grid. For example, in the case of (1+1)-d GNLSE, the transverse spatial variable is reduced to a set of spatial points separated by a step ∆x , while the longitudinal variable becomes a discrete set of values separated by the integration step h . Then, after discretization, the solution ψ is replaced by a discrete set of values ψ(j∆x, nh) = ψn j with n and j being integer numbers. This case can be generalized to other spatial dimensions, but we restrict ourselves to the study of the easiest case. In nite dierences the Laplacian is computed by means of a stencil 29 CHAPTER 2. NONLINEAR SCHRÖDINGER EQUATION IN A NUTSHELL (considering a 3-point stencil) ∂2ψn j ∂x2=ψn j−1−2ψn j+ψn j+1 ∆x2. (2.52) Once the Laplacian is computed, nonlinear term corresponds only to a point-to-point vector multiplication and the implementation of the method is complete. It turns out that this method is not only considerably unstable but also non conservative in the sense that it does not conserve the wave energy dened in equation (2.17) during the system evolution[25]. However it can be shown that an implicit method ψn+1 j−ψn j h=θfjψn+1+ (1 −θ)fj(ψn) (2.53) is conservative under the condition θ= 1/2 , which corresponds to a scheme commonly known as the Crank-Nicholson (CN) scheme. The split-step CN scheme can be obtained from the algorithm (2.43-2.45). The most common strategy is to use a simple Euler or RK4 method for the nonlinear sub-step and a CN for solving the linear sub-steps. In (1+1)-d, using the discretization previously discussed and the 3-point stencil, the linear step can be written as ψn+1/2 j−ψn j h=i ψn+1/2 j−1−2ψn+1/2 j+ψn+1/2 j+1 2∆x2!+iψn j−1−2ψn j+ψn j+1 2∆x2. (2.54) Then, grouping the n+ 1/2 and n terms, ψn+1/2 j h−iψn+1/2 j−1−2ψn+1/2 j+ψn+1/2 j+1 2∆x2=ψn j h+iψn j−1−2ψn j+ψn j+1 2∆x2 (2.55) the problem is reduced to the solution of the following linear system for ψn+1/2 Aψn+1/2=A∗ψn (2.56) with A=i 2∆x2              ∆x2/h + 2 −1 0 0 ... 0 −1 ∆x2/h + 2 −1 0 ... 0 0 ... . . ..... . . 0· · · −1 ∆x2/h + 2 −1 0· · · −1 ∆x2/h + 2              (2.57) This system is usually solved using iterative methods, that, depending on the problem, 30 2.5. NUMERICAL METHODS can either be fast or slow. Normally, when the matrix is sparse (for example it is a 3diagonal for (1+1)-d and 5-diagonal for (2+1)-d) there can be considerable speedups if the solver uses sparse matrix algorithms. A solver based on the CN method using GPU computing was implemented by Paulo Alcino in 2012 at INESC Porto. Even tough the CN scheme relies on the use of iterative methods that are slower than matrix multiplications of RK4, the performance of CN is usually better than RK4 [25] because RK4 is not conservative. In general, to improve the solutions obtained by RK4, it is necessary to use a smaller integration step h , which increases the number of integration steps needed for the simulation and drastically reduces the performance of the method. 2.5.2 Pseudo-spectral methods and the SSFM Pseudo-spectral methods rely on the utilization of the decomposition of the eld ψ in an orthogonal basis of functions, where it is easy to compute the linear sub-step [2]. From direct integration, the exact solution of the equation (2.42) is given by ψ(z+h, x) = exp hˆ D+ˆ Nψ(z, x), (2.58) with ˆ D=i 2∇2 ⊥ a linear operator relative to the dispersion and ˆ N=iF(|ψ|2) relative to the nonlinearities of the media. Using the Strang-Splitting algorithm we can reach an approximation for the solution of the GNLSE as ψ(z+h, x)≈exp h 2ˆ Dexp hˆ Nexp h 2ˆ Dψ(z, x). (2.59) This means that computationally the solution ψ(z+h, x) is calculated from ψ(z, x) by applying sequentially the operators exp hˆ D/2 , exp hˆ N and exp hˆ D/2 again. The Baker-Hausdor formula [2] for two operators ˆa and ˆ b is exp (ˆa) exp ˆ b= exp ˆa+ˆ b+1 2hˆa,ˆ bi+1 12 hˆa−ˆ b, hˆa,ˆ bii+... (2.60) and can give us a good insight of the error of the method. Indeed, applying the formula two times using ˆa=hˆ D/2 and ˆ b=hˆ N , it can be obtained exp h 2ˆ Dexp hˆ Nexp h 2ˆ D= exp hˆ D+hˆ N+Oh3hˆ D−ˆ N, hˆ D, ˆ Nii+... (2.61) 31 CHAPTER 2. NONLINEAR SCHRÖDINGER EQUATION IN A NUTSHELL suggesting that the dominant error term is of the order of h3 and that the method is accurate up to the second order. Before advancing, the spatial discretization must be considered. Considering a (3+1)- d system, the discretization on spatial coordinates can be introduced by a grid of integers (j, k, l) where 0≤j < Nx , 0≤k < Ny , 0≤l < Nt . Thus, any point in the continuous spatial space x= (x, y, t) 1 is represented by the corresponding X= (j∆x, k∆y, l∆t) . Is also useful to dene a discretization vector ∆X= (∆x, ∆y, ∆t) and the vector number of points N= (Nx, Ny, Nt) . The Fourier transform of the eld ψ is the decomposition of the eld ψ in an orthogonal basis of plane waves. Usually, the computational Fast Fourier transform (FFT) maps the complex-valued vector ψ into its frequency domain representation by ˆ ψ(z, k) = Nx−1 X j=0 Ny−1 X k=0 Nt−1 X l=0 ψ(z, X) exp (−ik·X), (2.62) where X is the discretized space vector. The discretization of the spatial domain re- ects as a discretization in the k= (kx, ky, kt) frequency domain. For even values of the components of the vector N , the discretization can be done in terms of three integers, (ˆ j, ˆ k, ˆ l) , within the limits dened by N ; however it is not linear like the spatial discretization. For example, kx is discretized under the formula kx=   2πˆ i Nx∆x, for 0≤ˆ i≤Nx 2 2π(ˆ i−Nx) Nx∆x, for Nx 2<ˆ i<Nx . (2.63) This allows to build a complete map between the eld ψ in the discretized direct space and the frequency discretized version ˆ ψ . The advantage of the using the Fourier transforms is that in the frequency space the Laplacian is algebraic, namely ˆ D(k) = −i 2k·k. (2.64) Therefore, it is possible to evaluate the linear sub-step in Fourier space using 1 In this sense we are considering t as a spatial variable. 32 2.5. NUMERICAL METHODS ˆ ψ(z, k) = FT{ψ(z, X)} (2.65) exp h 2ˆ Dψ(z, X) = F−1 T exp −ihk·k 2ˆ ψ(z, k) (2.66) where FT denotes the FFT operation. The nonlinear sub-step can be evaluated in direct space using the formula ψNL(z+h/2,X) = exp ihF(|ψ|2,X)ψ(z+h/2,X) (2.67) which completes the Split-step Fourier method (SSFM). It is important to notice that both the linear and nonlinear sub-steps are computationally solved in a discretized grid. This discussion about the discretization and the vectors X and k concludes that both equations (2.66) and (2.67) can be done by point-to-point calculations. In summary the SSFM algorithm is as follows: ˆ ψ(z, k) = FT{ψ(z, X)} ψ(z+h/2,X) = F−1 T exp −ihk·k 2ˆ ψ(z, k) ψNL(z+h/2,X) = exp ihF(|ψ|2,X)ψ(z+h/2,X) (2.68) ˆ ψNL(z+h/2,k) = FTψNL (z+h/2,X) ψ(z+h, X) = F−1 T exp −ihk·k 2ˆ ψNL (z+h/2,k). Before concluding this section we shall notice three important features of the SSFM. First, in a usual problem the interest is not to do only a single integration step but several of them. In that situation, it can be shown that, except for the rst step, we need only to compute one linear sub-step per integration step instead of two. In fact, the First Same As Last [62] property allows to concatenate the linear sub-steps as follows: ψ(z+ 2h, x)≈exp h 2ˆ Dexp hˆ Nexp h 2ˆ Dexp h 2ˆ D | {z } FSAL exp hˆ Nexp h 2ˆ Dψ(z, x) = exp h 2ˆ Dexp hˆ Nexp hˆ Dexp hˆ Nexp h 2ˆ Dψ(z, x). (2.69) 33 CHAPTER 2. NONLINEAR SCHRÖDINGER EQUATION IN A NUTSHELL Thus, considering several integration steps, the cost of this second order method is basically the same of the rst order method. Secondly, the calculation of the linear sub-step, the most time consuming step of the method, is done by using the FFT. If the dimensions N are all powers of 2, this method has a computational cost of order O(NtotlogNtot) where Ntot =Nx·Ny·Nt . As most FD methods rely in O(N2 tot) matrix operations, SSFM are usually faster, especially for multidimensional or large systems. Last but not least, the SSFM is conservative and normally admits larger integration steps than the FD methods for obtaining the same accuracy. In fact, the only situation where FD are preferable over SSFM is when the system is small and the dynamics of the envelope is fast, introducing a limitation to smaller integration steps [2]. The SSFM is not the faster neither the most accurate method for every situation [2, 25]. However, it constitutes the best performance GNLSE solver for the majority of the problems. Thus, it is our choice for the implementation of a high performance GNLSE solver using GPU computing. 2.5.3 Boundary conditions for the SSFM An important aspect of any solver of the GNLSE are the boundary conditions. In general, a soliton like solution of the GNLSE extends well beyond the limits of the simulation box even though in that portion of space the amplitude of the eld can be very close to zero, and therefore negligible. Also, it is possible that the solitons propagate to close proximity (and scatter through) to the boundaries of the simulation box. The boundary represents a discontinuity in the simulation box and can interact with the soliton-like pulses yielding diverse, and many times unwanted, eects. These may include reection and numerical dispersion, depending on the type of the solver being developed. Therefore, great care must be put in addressing the boundary-eld interaction. This is specially important when considering problems where the medium is supposed to be innite or the relevant interaction is restricted to a small region of space after which the eld evolves into far regions, as occurs during soliton scattering. If no attention is put to boundary conditions then extremely large simulation boxes must be used, which are costly in terms of computational resources, eciency and simulation times. To avoid these problems, several types of numerical solvers for the GNLSE have been developed which control the physics of boundary-eld interaction, namely allowing to produce periodic, reective and absorbing boundaries. Periodic boundary conditions allow that, when the eld reaches one of the boundaries 34 2.5. NUMERICAL METHODS of the simulation box, then it emerges from the opposing boundary. This implies that a two dimensional simulation box corresponds topologically to a torus. This type of boundary condition arises natively from the SSFM algorithm given its calculation of the linear step in the Fourier space. Reecting boundary conditions force the eld that reaches one of the boundaries to bounce back. This type of boundary condition can be implemented by dividing the simulation box into two domains, one at the center where the simulation occurs and another corresponding to a thin layer of points, as shown in gure (2.4). The idea is to replace the original GNLSE with an altered version, say i∂ψ ∂z +1 2∇2 ⊥ψ+F(|ψ|2)ψ+Vrψ= 0 (2.70) where Vr=   vr at the boundary layer 0 elsewhere (2.71) where vr is a real-valued constant larger than any of the other terms in the GNLSE. The absorbing boundary conditions correspond to the case where a eld reaching the boundaries of the box is totally (or almost totally) absorbed and disappears from the simulation box. This boundary condition is specially indicated to simulate the propagation of solitons in innite or very large domains. Again this boundary condition is obtained by replacing the original GNLSE with i∂ψ ∂z +1 2∇2 ⊥ψ+F(|ψ|2)ψ+iVaψ= 0 (2.72) where Va is a positive real valued function, chosen to maximize the absorption of the radiation. The most common absorbing potentials are the Gaussian [62] Va=Ae−(X−X0)2 2w0 (2.73) and the hyperbolic tangent [67] Va=A(1 + tanh (w0(r−x0))) . (2.74) Parameters A , w0 and x0 are characteristic of the absorbing potential whose choice depends on the problem and must be optimized to maximize the absorption of the outgoing radiation. 35 CHAPTER 2. NONLINEAR SCHRÖDINGER EQUATION IN A NUTSHELL Figure 2.4: Visual scheme of the simulation box, describing the concept of boundary layer. The blue interior of the box is where the eld evolves. 2.6 Concluding remarks By the end of this chapter we have presented some of the most important aspects of the structure, solutions and numerical methods of the GNLSE. Particularly, the introduction of the SSFM provides the cornerstone for the next chapter, where we describe the implementation of our GNLSE solver based on GPU computing, GASE. In the following chapter we discuss how the SSFM can be adapted to operate in a GPU architecture, taking to account the allocation of data into the memory available. Also, we address how the resources of a computer can be used in heterogeneous programming to boost code eciency. 36 3 Implementation of the GPU-based GNLSE solver In the previous chapter we reviewed the basic framework of the GNLSE, including the type of equations and situations that have been studied, the analytical methods used, and some of the most relevant numerical methods developed to solve it. In this chapter we focus on the implementation of a GNLSE solver based on a SSFM using GPU computing. We start by analyzing the advantages and disadvantages of this approach, especially when compared with CPU-based computing, which constitutes the most common base approach in the past. The core of this chapter is devoted to the description of the algorithm used, its computational implementation and to its performance analysis. As it shall be shown, we have obtained computational speedup factors of almost 100 when compared to the same CPU-based solver, demonstrating the high potential of GPU computing for numerical analysis of the GNLSE. 3.1 Simple problem, high computational time The calculation of solutions of the GNLSE in systems with dimensionality higher than (1+1)-d and specially, if it involves complex geometries and higher order nonlinearities, is a very large computational problem. This is mostly due to the large number of points of the spatial mesh used to sample the eld and the local optical properties, which need to be determined and computed on each time step. Using serial programming, such as used in single core CPUs, the solution of this type of problems requires vast running times as most calculations must be done sequentially for a very large set of sampling points of the eld. To illustrate the immense challenge at hands, consider a simple problem, consisting of the investigation of the dynamics of a supergaussian (2+1)-d soliton in a cubic-quintic media limited in a small square domain with side length a having a small circle - we will call it a hole - of linear material at the center, with diameter of value a/25 , as shown in gure 3.1. Consider also that the soliton has a characteristic size of a/4 and 37 CHAPTER 3. IMPLEMENTATION OF THE GPU-BASED GNLSE SOLVER based on GPU appears to have a bright future ahead. 3.4 Some memories must be kept closer than others: memory considerations for GPU computing Other key aspect of hardware that is necessary to account when programming is memory allocation. In more detail, dierent types of data needed for the calculations must be stored in distinct types of memory of the GPU architecture, if optimal (or close to optimal) performance is to be achieved. Inside the GPU card there are three types of memory, with dierent access times: • Shared memory: small amount of memory (48 kB) with fast access but not accessible by all cores in a GPU; • Constant memory : small amount of memory (64 kB = 8192 double precision oats) with fast access and accessible by all cores. The speed makes it the ideal to store constant data, such as physical constants used by the solver. • Global memory: big amount of memory (few GB, see table 3.1) with slow access and accessible by all cores. Its size makes it ideal to store the eld distribution and the nonlinearities. As mentioned, global memory, the RAM of the device, is the memory used to store elds. It should be noticed that the amount of memory of the GPU is xed from factory and cannot be augmented afterward, unlike the CPU, where one can always add larger RAM memories. This can limit the capabilities of the GPU in larger simulation systems. Generally, for each GB of RAM memory, 222 double precision numbers can be stored and operated in the GPU. Hopefully, future GPU models will have increasing memory that will make GPUs more capable of addressing larger simulations. Finally let us give a word about data storage, namely about saving the numerical data produced during calculation in the GPU to a more permanent storage, in this case a hard disk, for later analysis. The transference of data between dierent groups of cores in the GPU is fast (the bandwidth concept mentioned in the last section) resulting in a short latency time between two time steps of integration of the GNLSE. However, transferring that information from the GPU to be stored in an external hard disk is time consuming, and must be avoided during the calculations. This problem will be addressed in more detail in the next sections. 44 3.5. COMPATIBILITY ISSUE 3.5 Compatibility issue The use of GPUs for general purpose computing, and more specically for scientic computing, is something new in computer sciences and is a part of a broader concept called heterogeneous computing. To put it simple, heterogeneous computing aims to use all the resources available in a machine (GPUs, CPUs or others) in an integrated way to do massive computing. The early precursors of this concept included the engineers of NVIDIA, one of the most important graphics cards company and the developer of CUDA. Although CUDA is currently perhaps most advanced platform for heterogeneous computing, it operates only on specic GPU cards from NVIDIA. There is some irony in this fact: heterogeneous computing has a main goal in promoting portability between dierent devices but CUDA and programs written in CUDA can only be used in hardware from a specic manufacturer. Therefore, our code is not portable to machines that have other types of GPU cards. To address this limitation, an alternative to CUDA is being developed in recent years by a consortium of GPU card and CPU manufacturers (especially AMD) called OPEN CL. The idea is that OPEN CL can become a standard language for heterogeneous computing capable of operating in any type of device. Unfortunately, OPEN CL is still in its earlier versions and even though it can already compete with CUDA in terms of managing parallel computation in distinct devices, it still lacks numerical packages that support scientic computing. This short-come of OPEN CL is expected to be overcome in future years as it becomes more used and software developers bridge the gap between it and CUDA. It is not possible to consider OPEN CL to support the development of a GNLSE solver yet, however the structure of GASE should (in principle) be easily transposed from CUDA as OPEN CL reaches later stages of development. 3.6 Implementation of the GNLSE solver This section is devoted to explain the general structure of the numerical code developed during this dissertation's research. 3.6.1 Outline of the code The code is basically composed of three parts, each corresponding to a specic stage of the calculation as shown in gure 3.3. 45 CHAPTER 3. IMPLEMENTATION OF THE GPU-BASED GNLSE SOLVER The rst part encompasses the initialization of the data structures that store the data. These include the lists containing the spatial coordinates, the value of the eld and the optical properties and parameters in each point of the mesh of the simulation box. Within these structures are also included the data necessary to implement the boundary conditions. During the second part of the code the actual numerical calculations necessary to integrate the GNLSE are performed using the SSFM. In a very simple way, it consists in a loop which repeats the integration step. This is the most time consuming part of the code and the duration of each step will be discussed in the following sections. During this second stage, the results are stored into the hard disk to be used in the third and last stage of the code where the data is analyzed. In order to reduce both the memory requirements of the hard disk and the latency of the process of transferring data from the memory to the hard disk, data is only stored after doing an user-dened number of integration steps. While the rst and second part are written in C++ with CUDA, this last part of the code is written in Python (instead of CUDA) by reasons of convenience, since not only Python allows to produce graphics with high quality but also is much easier to program, given that it is a scripting language. It should be noticed that CUDA and the associated numerical toolkit allow us to do most of the calculations both in the GPU and CPU without signicative changes to the structure of the code. In order to compare our solver GASE with a CPU version of a GNLSE solver we adapted and developed also a CPU solver using the same structure and the Fftw library for the Fourier transforms. Therefore, and to benchmark the GPU computations relative to the CPU version, this code is prepared to operate in both platforms with specications that we discuss in the following section. 3.6.2 Integration step routine One of the most important features of GASE is that it can address systems with higher dimensionality, i.e., (2+1)-d and (2+1+1)-d domains (although it can be easily extended for higher dimensions if necessary). This implies the discretization of the domain into a regular mesh of points on which the eld and the optical properties of the medium must be evaluated. In principle, these values could be stored in multidimensional lists of data (similar to tensors) but this is not the most ecient way to use data in GPU computing. In fact, the numerical libraries built on CUDA and used to develop GASE operate only with one dimensional lists, including the numerical library which takes care of the FFT. Therefore, the multidimensional list must be spanned into a single 46 3.6. IMPLEMENTATION OF THE GNLSE SOLVER Figure 3.3: Succinct description of the code structure. The code is divided in three parts: the rst initializes the data and the simulation box, the second is the integration routine and the third is the post-simulation analysis of the data. 47 CHAPTER 3. IMPLEMENTATION OF THE GPU-BASED GNLSE SOLVER one dimensional list (which is quite simple to do). More importantly, it is necessary to convert the identifying index of the data in the multidimensional list (j, k, l) into the index in the one dimensional list I=j+k×Nx+l×Nx×Ny . The conversion between the index I to (j, k, l) is easily done with the help of oor and mod functions, using the expressions j=Imod (Nx) k= oor I Nxmod (Ny) l= oor I NxNy In the simulation mesh, this allows to compute the spatial coordinates of the point to which the data pertains to. For example in three dimensional mesh, coordinates are obtained as X= (j∆x, k∆y, l∆t) . Also, having this methodology in mind, the nonlinear step can be easily calculated from expression (2.68). The FFT routine transforms the eld 1D lists with the eld data Ψ(I) into another, say ˆ Ψ(I) , where each element corresponds to a specic spatial frequency. According to the documentation of Cut , the relation between true index of the list I and the corresponding spatial frequency is given computing (j, k, l) and then using the correspondence formulas (here shown for kx but the same is extended to the other dimensions) kx=   2πj Nx∆x, for 0≤j≤Nx 2 2π(j−Nx) Nx∆x, for Nx 2< j < Nx allows us to compute the Laplacian of the eld and then computing the linear step using expression (2.68). The second stage of the code is basically a loop which operates consecutive linear and nonlinear steps of the SSFM on the eld data and is depicted in gure . However, it should be noticed that these steps are grouped (blocks) into sequences of about a few hundred steps (steps per block) during which no data is registered in the hard disk. As explained before, the transference of data from the GPU to the disk is time consuming and if done after each step it would eat away the performance of doing parallel computing. Also, the data on GPU cannot be transferred directly to the hard disk. Rather, it is rst copied to the RAM of the host and then it undergoes some simple processing in the CPU. Namely, the eld intensity and phase are computed from the complex eld 48 3.7. CONCLUDING REMARKS amplitude prole. Only then, the data is transferred from the RAM to the hard disk. Together, the computation of eld intensity and phase and the transfer of the data from the RAM to the hard disk, constitute a very slow process. This would kill some of the performance of the GPU code so a new feature was added to GASE, allowing to perform this operations in the CPU while GPU is already and simultaneously running the next step of the SSFM. It is important, but not mandatory, to choose sequences of SSFM steps suciently long to give time for the CPU to nish its task, but short enough to allow a good insight of the evolution of the eld. As a result of using the CPU for part of the workload, the GPU is relieved from part of the calculations and from latency times during data storage, resulting in faster computing processes as a whole. In fact, this is a good example of heterogeneous computing where all the resources of the machines are used to promote overall eciency. 3.6.3 Code features Before ending this chapter it is important to make a synthesis of the capabilities of GASE. Therefore, summarizing all the work developed, GASE is currently able to perform integrations of any given GNLSE in (1+1)-d, (2+1)-d and (2+1+1)-d geometry, with any user given initial condition. The nonlinearities are chosen by the user and might be of arbitrary power and can be either a constant number or a spatial distribution. The nonlinearities can also be given by a function to have a dependance on the propagation distance. All of these systems can be simulated in a simulation box either using periodic, re- ective or absorbing boundary conditions. Moreover, there is also a recent and still under development, feature that allows the simulation of two coupled GNLSE. Therefore, GASE is a powerful tool and we believe it to be ready to investigate the majority of the most state-of-the-art problems in solitons and in GNLSE subject. 3.7 Concluding remarks In the beginning of this chapter we have introduced a problem and conclude that simulating such system would imply very high computational running times, which make the problem very hard to investigate. We proceed trying to develop a new tool to address the problem and presented the GPUs as a recent technology that could be useful for such systems. 49 CHAPTER 3. IMPLEMENTATION OF THE GPU-BASED GNLSE SOLVER Figure 3.4: Integration routines for both versions of the solver. Figure a) shows the structure of the evolution routine for GASE, where it is possible to see the additional memory transfer needed but also the parallel structure of computations running both in CPU and GPU. Figure b) describes the integration procedure for the CPU version of the solver. 50 3.7. CONCLUDING REMARKS After a small discussion on GPU computing framework, we described succinctly the structure of GASE, the GNLSE solver based on GPU computing developed during this dissertation. Full details were not given as this description is intended to be simple to the reader, avoiding to enter in the complicated world of CUDA programming. At the end of this chapter an obvious question arises: what is the speedup, if any at all, that GASE can achieve? To answer the question we recover the initial problem introduced at the begin of the chapter. The propagation of a supergaussian of order m= 1.9 with parameters described in section 2.4.1 was done using GASE. In conformity with the initial description of the problem, shown in gure (3.1), we choose a simulation box with limits [0,120] ×[0,120] . Also, the center hole has a radius of 2.4 and refractive index n= 0 . After some initial numerical tests we found that the system could be solved in single precision using an integration step of h= 0.02 . The total number of integration steps used is 100 000. Figure 3.5: Sequence a)-d) shows a collision of a supergaussian state with velocity µ= 0.3 with a hole of radius 2.4 . The light state emerges as two smaller and low intensity beams after the scattering. ( Pdf version only - click twice on sub-gure a) for a small clip of the simulation) A series of simulations were performed and results are presented in gures (3.5 - 3.9). The physics of the results are interesting as they resemble a collision of a drop of liquid with a circular object, which is not unexpected given the liquid behavior of high power supergaussian states [58]. Also, the results change depending on the initial velocity and the hole radius, and a plethora of dierent behaviors is achieved. For the initial hole radius, a soliton with velocity µ= 0.15 rebound on the hole, but with higher velocity, such as µ= 0.3 and µ= 0.5 , is decomposed in two smaller and lower intensity light beams. Diminishing the hole radius to 2.0 , a light beam with a velocity µ= 0.3 is momentarily divided into two low power light beams but regroup as one after the collision with the hole. A smaller hole of 1.0 seems even to not aect the propagation. These results are interesting but are not our main goal in dissertation. Instead, they 51 CHAPTER 3. IMPLEMENTATION OF THE GPU-BASED GNLSE SOLVER illustrate the power of the developed GNLSE solver. As a matter of fact, we have not performed only the 5 presented simulations but a series of simulations that, including the earlier investigations and those which were wrongly set up, took more than one day to run in our computer. An initial comparison between GPU and CPU versions of the solver show that GASE runs the problem almost 100 times faster than the CPU version. The conclusion is obvious: if we had chosen to investigate this problem in a normal CPU then it would have taken almost 100 days, a third of the duration of this dissertation. But how fast can we go and for which problems? The answer is addressed in the next chapter, where we perform the benchmarking of GASE. Figure 3.6: Sequence a)-d) shows a collision of a supergaussian state with velocity µ= 0.5 with a hole of radius 2.4 . The light state emerges as two smaller and low intensity beams after the scattering, with a dierent angle than the situation with velocity µ= 0.3 . ( Pdf version only - click twice on sub-gure a) for a small clip of the simulation) Figure 3.7: Sequence a)-d) shows a collision of a supergaussian state with velocity µ= 0.15 with a hole of radius 2.4 . The light state collides and is reected by the hole. ( Pdf version only - click twice on sub-gure a) for a small clip of the simulation) 52 3.7. CONCLUDING REMARKS Figure 3.8: Sequence a)-d) shows a collision of a supergaussian state with velocity µ= 0.3 with a hole of radius 2.0 . The light state after an intermediate division in two light states collapses again in one high power state. ( Pdf version only - click twice on sub-gure a) for a clip of the simulation, where can also be seen that the state became trapped between two consecutive holes, that we simulate in the same box using periodic conditions) Figure 3.9: Sequence a)-d) shows a collision of a supergaussian state with velocity µ= 0.3 with a hole of radius 1.0 . The light state emerges almost as if not been scattered. ( Pdf version only - click twice on sub-gure a) for a small clip of the simulation) 53 CHAPTER 4. BENCHMARK OF GASE Figure 4.2: Single (top) and Double (bottom) precision benchmarks for simulations of the (1+1)-d NLSE, with the results for computational time per step (left) and the corresponding speedup in comparison with the CPU version of the code (right). 60 4.3. BENCHMARK OF THE CODE: GPU VERSUS CPU Number of Steps h N GPU (s) CPU (s) Speedup 100 000 0.01 28 3.4 6.2 1.8 10 000 0.1 212 0.8 10.8 13.5 1 000 1.0 218 2.8 82.1 29.3 1 000 1.0 222 48.4 1575.6 32.5 Table 4.5: A collection of results for simulation times and speedup of the solver for the (1+1)-d NLSE using double precision, for both GASE (running in the desktop GPU) and the CPU-based version of the solver. understand the performance of this code in comparison with another code written in a scripting language, we developed a similar version of the code in Python. For the same simulation, Python took more than 11 days, more than 20 000 times slower than our GASE code, which provides a demonstration that Ron Caplan's NLSEmagic [12] is strongly limited by using a scripting language. It is noticeable that the performance when using double precision is reduced to less than a half when compared with single precision performance. This was expected but still good speedup results are obtained. Single performance can be used for some problems but should be used carefully, and only after being sure that results obtained are the same to those obtained using double precision. Simulations involving many integrations steps should be done using double precision, in order to keep numerical errors under control. 4.3.2 (2+1)-d speedup results For the (2+1)-d speedup tests we choose to run a simulation of a cubic-quintic medium described by the GNLSE of equation (2.2) with F(|ψ|2) = |ψ|2− |ψ|4 . For the initial condition, it is used a supergaussian of order m= 1 , with parameters and form given by the equations introduced in section (2.4.1). The nal propagation distance is considered ztotal = 10 and the integration step is now xed at h= 0.01 . The simulation box has a constant spatial discretization ∆x= ∆y= 0.2 and is dened over a grid of Nx×Ny points, corresponding to the number of the mesh points in the x and y dimension, respectively. Again, speedup results are analyzed with the help of a time per step variable, as done for the (1+1)-d case. Running times and speedup results are shown in gure (4.3) and table (4.6) and again signicant speedups were obtained, reaching in double precision a top speedup factor of just over 40. For example, for a computational grid with N= 211 ×211 points, GPU took just one minute to solve the system while CPU took more than 40 minutes. If we consider a simulation that took a day to solve in GPU 61 CHAPTER 4. BENCHMARK OF GASE Figure 4.3: Double (top) and single (bottom) precision benchmarks for simulations of the (2+1)-d GNLSE for a cubic-quintic media, with the results for computational time per step (left) and the corresponding speedup in comparison with the CPU version of the code (right). - just consider for example the same mesh and nal propagation time of ztotal = 14400 - solving in the CPU will took more than one month, which is a prohibitive duration for the majority of research. It is important to note that these results are even better than the (1+1)-d case. This is related to the fact that now we consider two nonlinear terms and the system is then more complex. This becomes more clear with the results of gure (4.4), corresponding to the same initial value problem but now solved in a simulation box with boundary conditions as introduced in section 2.5.3. For this it is used an absorbing potential Va= (1 + tanh (r−L)) , with L the limits of the simulation box . Speedup factors for this problem show a top value of over 80 which is really signicative. 62 4.3. BENCHMARK OF THE CODE: GPU VERSUS CPU Figure 4.4: Single precision benchmarks for simulations of the (2+1)-d GNLSE for a cubic-quintic media with absorbing boundaries, with the results for computational time per step (left) and the corresponding speedup in comparison with the CPU version of the code (right). Number of Steps h N GPU (s) CPU (s) Speedup 1 000 0.01 28×28 1.1 34.1 31.0 1 000 0.01 29×29 3.8 142.9 37.6 1 000 0.01 210 ×210 14.7 587.5 39.9 1 000 0.01 211 ×211 61.0 2506.5 41.1 Table 4.6: A collection of results for simulation times and speedup of the solver for the (2+1)-d GNLSE for a cubic-quintic media, using double precision, for both GASE (running in the desktop GPU) and the CPU-based version of the solver. As obtained for (1+1)-d, gure (4.3) shows that single precision simulations are twice faster than double precision, with overwhelming speedup factor of almost 80. In (2+1)-d cases is common to be interested in evolving the system much less steps than the (1+1)-d problems. Therefore, single precision can be used for some cases but again, this must be done with caution. 4.3.3 (2+1+1)-d speedup results For the three dimensional simulations, or as called before, the (2+1+1)-d GNLSE, we use a supergaussian state of order m= 1 as one in the last section. The nal propagation distance is considered zfinal = 10 and the integration step is xed as h= 0.01 . We use the discretization of the space as before, ∆x= ∆y= ∆t= 0.2 and dene a simulation box represented by the mesh Nx×Ny×Nt . 63 CHAPTER 4. BENCHMARK OF GASE Figure 4.5: Double (top) and single (bottom) precision benchmarks for simulations of the (2+1+1)-d GNLSE for a cubic-quintic media, with the results for computational time per step (left) and the corresponding speedup in comparison with the CPU version of the code (right). 64 4.4. CONCLUDING REMARKS Number of Steps h N GPU (s) CPU (s) Speedup 1 000 0.01 25×25×25 0.7 16.4 23.4 1 000 0.01 26×26×26 3.7 133.5 36.1 1 000 0.01 27×27×27 29.8 1112.7 37.3 Table 4.7: A collection of results for simulation times and speedup of the solver for the (2+1+1)-d GNLSE for a cubic-quintic media, using double precision, for both GASE (running in the desktop GPU) and the CPU-based version of the solver. Results for speedups compared with serial CPU performances are shown in gure (4.5) and table (4.7). It can be observed that speedups are slight smaller than those obtained for (2+1)-d case, which is possibly related with the lower performances of the GPU for computing Fourier transforms in three dimensions. Still, a considerable speedup factor of over 35 is obtained for the larger simulation boxes. For single precision, speedup factor is almost twice of that obtained for the double precision, which as said in sections before, was already expected. 4.4 Concluding remarks Along this chapter we presented the tests done to validate the GNLSE solver developed during this dissertation - GASE. We started by comparing the results obtained with GASE with the exact solution, in situations where an analytical solution exists. Afterward, we proceed to a series of tests intended to benchmark the performance of GASE in comparison with the CPU version developed. We have obtained maximum speedups of around 80 for single precision and around 40 for double precision. Also we have observed that increasing the dimensionality of the system or the complexity leads to a larger speedup factor. Therefore, the results of this chapter are important as they demonstrate the realization of the goal of this dissertation: the development of a high performance GNLSE solver, capable of addressing eciently problems in multidimensional or complex systems, using only low-cost computers and based on GPU computing. 65 5 Case study 1: Lightons: phonons with Light In the present and the following chapter we present two case studies used to evaluate the potential of our implementation of a GNLSE solver based on GPU computing. These case studies have not been fully explored in terms of scientic analysis for two reasons. First, because the main objective of this dissertation is the development of a GPU based GNLSE solver, which is by itself very challenging and required much eort. Secondly, each of these topics can constitute a subject of a master thesis on their own. However, we already present many and interesting new results. In this chapter we analyze a chain of interacting solitons which can support collective excitations similar to phonons in chains of masses. This problem is exceptionally dicult to address using conventional SSFM implemented on a single CPU due mainly to the large number of discretization points and integration steps required to simulate several dozens of solitons and investigate the continuous limit of this excitations. Moreover, we also present results for (2+1)-d chains of solitons, not just (1+1)-d, which truly demonstrate the potential of GPU computing. 5.1 Motivation At a fundamental point of view, soliton appear as stable solution of nonlinear wave equations, such as the Nonlinear Schrodinger equation (NLSE), balancing wave dispersion with nonlinear wave interaction, and resulting in a coherent wave package that can behave both like a particle and a wave. In some sense, solitons bridge the gap between waves and particles. Being a wave phenomenon in its fundamental aspects, solitons can behave like a particle by avoiding the eects of dispersion and keeping most of its intensity conned to a small region of space. Also, like particles, in some situations they can scatter each other as if they where solid blocks of light. 67 CHAPTER 5. CASE STUDY 1: LIGHTONS: PHONONS WITH LIGHT Just as waves can exhibit particle like behavior, particles can also support waves, such as mechanical waves. For example, a chain of interacting particles can display collective oscillatory motions which depend on the nature of the interaction forces. Elastic linear interaction allows small displacements of the particles from their equilibrium positions to form mechanical waves, known as phonons, that propagate throughout the chain. As waves, phonons can not only propagate but also superpose and interfere. From this duality, where nonlinear wave packages display particle behavior and particles can support waves, arises a question: can solitons aligned in a chain support mechanical waves similar to phonons? In this chapter, we focus our attention in the collective oscillations of a 1-dimensional chain of N optical solitons, supported by a medium with a cubic nonlinearity as described by the NLSE. We investigate the nature of these oscillations using GASE and extend the models existing in the literature not only by identifying its limitations but also by presenting illustrations of their extension to (2+1)-d systems. 5.2 Physical model The evolution of optical solitons in a media with cubic non-linearity is governed by the (1+1)-dimensional NLSE i∂ψ ∂z +1 2 ∂2ψ ∂x2+|ψ|2ψ= 0 (5.1) where z is the propagation distance and x is coordinate along the transverse direction. The single soliton solution for this equation is ψj(x, z) = 2νjsech {2νj(x−¯xj)}exp {i2µj(x−¯xj) + iδj} (5.2) where the parameters νj , ¯xj , µj , and δj refer to the amplitude, the position, the frequency and the phase of the soliton, respectively. We consider the limit of large distance and small overlap between consecutive solitons, where the N-soliton solution for the equation (5.1) can be given as a sum of N single soliton solutions ψ(x, z) = N X j=1 ψj(x, z) (5.3) When replacing (5.3) into equation (5.1) it becomes clear that the nonlinear term is responsible for the interaction between solitons. To solve the resulting equation, we use an ansatz and separate equation (5.1) into the following N coupled equations, one for 68 5.2. PHYSICAL MODEL each single soliton: i∂ψj ∂z +1 2 ∂2ψj ∂x2+|ψj|2ψj=−Rj[ψ] (5.4) It is simple to see that the sum of the solutions of each of these equations is a solution of equation (5.1). Each of these equations is just a NLSE with a small perturbation Rj associated to solitonic interaction. Considering that only rst neighbors can interact and expanding the rst order in the overlap Oψjψ∗ j+1, ψjψ∗ j−1 we get Rj[ψ] = 2 |ψj|2(ψj+1 +ψj−1) + ψ2 jψ∗ j+1 +ψ∗ j−1 (5.5) The initial problem is now reduced to a form which is suitable to be analyzed using a quasiparticle approach which results in the following evolution equations for the parameters of each j-th soliton[81] dµj dz = 16ν3[cos (φj+1,j) exp (−∆j,j+1)−cos (φj−1,j) exp (−∆j−1,j)] (5.6) d¯xj dz = 2µj−4ν[sin (φj+1,j) exp (−∆j,j+1)−sin (φj−1,j) exp (−∆j−1,j)] (5.7) dδj dz = 2 ν2+µ2 j+ 2µd¯xj dz −2µj+ (5.8) +24ν2[cos (φj+1,j) exp (−∆j,j+1) + cos (φj−1,j) exp (−∆j−1,j)] Here φj,l = 2µ( ¯xl−¯xj) + Ψj,l is the complex phase between solitons, with Ψj,l =δj−δl , the quantity ∆j,l = 2ν|¯xl−¯xj| is the spacing between two solitons and ν and µ are the mean values for the amplitude and frequency, respectively. The derivation of these equations takes into account the assumptions that the overlap is small, i.e. ν|¯xl−¯xj|  1 , and that the solitons have similar frequencies and amplitudes, i.e. |µj−µl|  µ , |νj−νl|  ν . Also it is assumed that the amplitude νj of each soliton is approximately constant, a fact which occurs in the adiabatic limit and is supported by small amplitude variation observed in numerical simulations. We also adopt a conjecture commonly made [81], where it is assumed that if the system is initialized with consecutive solitons having a phase dierence of 0 or π , this dierence remains constant during the motion of the solitons. This assumption leads us 69 CHAPTER 5. CASE STUDY 1: LIGHTONS: PHONONS WITH LIGHT Figure 5.6: Evolution of the collision of a soliton with velocity µ= 0.2 with a soliton chain with 4 solitons initially separated by ∆ = 10 . The trapping potential allows to obtain results that resemble the Newton's cradle. V(x) = 3 ×10−5(x−100)2 is added to the system. The results are shown in gure (5.6). The chain now responds similar to a Newton's cradle, where each soliton transfers its momentum to the following and replaces it in the chain. This behavior is closer to a particle than to the mechanical-like waves, previously observed in the simulation. This appears to suggest that the system can exhibit many types of behavior, from wave-like to particle-like. This is a result of the nonlinear nature of the system which allows for distinct types of phenomena. 5.3.2 Soliton chains in (2+1)-d Finally we show that these eects are not limited to soliton chains in (1+1)-d but can be generalized to higher spatial dimensions. In gure (5.7) we show a chain in a (2+1)-dimensional space. Here the spatial solitons are Gaussian shapes supported by a cubic-quintic media and the early results show that it is possible to propagate energy and momentum with wave-like phenomena in a chain of solitons. 76 5.4. CONCLUDING REMARKS 5.4 Concluding remarks We have explored the dynamics of an 1-dimensional chain of optical solitons in a cubic nonlinear media described by the NLSE. Using perturbative methods we show that the interaction between consecutive solitons could support wave-like oscillations of the positions of the optical solitons in the chain. We name these waves lightons: phonons of light. Numerical simulations show it is possible to create standing wave oscillations in the considered system that follows a predicted relation of dispersion. For small propagation distances and displacement amplitudes, the system revealed a controllable sinusoidal motion for the soliton position which is obtained only with light-light interaction. For higher amplitudes, the solitons enter a collisional regime and their behavior changes from a phonon type excitation to something that resembles a Newton cradle. 77 CHAPTER 5. CASE STUDY 1: LIGHTONS: PHONONS WITH LIGHT Figure 5.7: Evolution of an 1-dimensional chain of 2-dimensional spatial solitons with ∆ = 20 . 78 6 Case study 2: Soliton-soliton scattering in (2+1)-d In the previous chapter we analyzed an example of how the developed GNLSE solver can be used to study problems with high complexity, specically, with a large number of solitons. In this chapter we demonstrate the potential of this computing approach in addressing systems with higher dimensionality. In particular, we study the collisions of (2+1)-d solitons in the physics of the dynamics in two dimensional plane. Unlike the case with (1+1)-d collisions, it is possible to study a wide range of collision parameters, including the impact parameter, the phase dierence, the initial energies, etc.. We have chosen to analyze the transition between solid and liquid-light behavior and their inuence on the collision dynamics. 6.1 Motivation The example of the previous chapter illustrated how solitons can have a versatile nature, exhibiting both wave and particle-like behavior. However, the diversity of soliton behavior extends well beyond that. In systems with higher dimensionality, (2+1)-d or more, soliton-like solutions of the GNLSE must be supported by a mix of higher nonlinearities, as discussed in chapter 2. Depending on the relative global phase of the soliton and type of nonlinear media, two solitons can scatter each other like rigid bodies, go right through each other like a wave, or coalesce in a wider soliton, like the coalescence of two droplets of light. Exploring the dynamics of interacting solitons is a process with major interest in nonlinear optics. However, most studies have been restricted to the one dimensional simulations, where the simulations can be performed in a state-of-the-art computer. Exploring systems with higher dimensions was limited to the use of clusters and supercomputers. The use of GPU computing in GASE allows to explore these situations and take into account a wider and richer diversity of situations. In this chapter we explore 79 CHAPTER 6. CASE STUDY 2: SOLITON-SOLITON SCATTERING IN (2+1)-D the scattering of solitons in cubic-quintic media in (2+1)-d scenarios. As a result, we are able to explore the interplay between solitons phase, their original velocities and impact parameters, thus providing a better insight of the physical processes. 6.2 Physical model The starting point of analysis is again the GNLSE i∂ψ ∂z +1 2∇2 ⊥ψ+F(|ψ|2)ψ= 0, (6.1) which describes the evolution of the dimensionless amplitude of a light eld ψ in a nonlinear media with properties given by the function F(|ψ|2) . Here z is the longitudinal coordinate parallel to the propagation and ∇2 ⊥ is the Laplacian in the transverse directions. As discussed in chapter 2, for systems with (2+1) or more dimensions it is necessary to consider special types of nonlinearities for soliton-like solutions is assumed that the solutions of equation (6.1) can be described by a general soliton solution q(x, z) = A(z)g[B{x−¯x(z)}] exp(−ik(z){x−¯x(z)}+iθ(z)), (6.2) where A , B , g , k , θ , ¯x(t) are the amplitude, the width, the shape, the frequency, the phase and the center of the soliton, respectively. The existence and the dynamics of the soliton can be then be investigated using variational methods in terms of the variation of these parameters for specic forms of F , corresponding to dierent nonlinearities. Media with cubic-quintic nonlinearities are known to support solitons in more than one dimension and have been extensively studied because they are described by a simple nonlinear potential of the form F|ψ|2=|ψ|2− |ψ|4 . In this case, it is possible to obtain approximate solutions in the form of supergaussians pulses [18], dened by ψ(r, z) = Aexp −B2(r−¯r)2mexp (iδz), (6.3) where r is the vector with transverse coordinates and m a parameter related to pulse energy. For the soliton to be stable, these parameters must be mutually related by the following conditions, 80 6.3. SCATTERING OF COLLIDING (2+1)-DIMENSIONAL SPATIAL SOLITONS A=s3 21/m 3 2 m−ln 2 2m−ln 3 (6.4) B= 1 As21/m Γ(1 + 1/m) 2m−ln 3 ln(4/3) !−1 (6.5) δ=2B2m Γ(1 + 1/m) m+ ln(2/3) ln(4/3) , (6.6) expressions derived from those presented in section 2.4.1 and from [18]. As the system is non integrable, studying the dynamics of the solitons is only possible using either perturbative methods or numerical simulations. Here, we focus on the numerical simulations because they provide a more direct way of studying the interactions of solitons in wider range of situations. 6.3 Scattering of colliding (2+1)-dimensional spatial solitons We are interested in the computational analysis of two (2+1)-dimensional soliton collisions in the xy plane. The spatial solitons dynamics is described by the cubic-quintic GNLSE . The solitons have parameter m= 1 , and the other parameters set according to equations (6.4-6.6). In this study, we initialize the solitons with a phase dierence δ1−δ2= 0 or π and with opposing but equal velocities |k1|,|k2|=k (see gure (6.1)), by multiplying the supergaussian shape () by exp (−ik(z){x−¯x(z)}+iθ(z)) . Figure 6.1: Graphical description of the problem analyzed. Here, impact parameter b was exaggerated for better comprehension. Simulation box has limits [0,80]× [0,80] and a mesh of N= 210 ×210 points was used. 81 CHAPTER 6. CASE STUDY 2: SOLITON-SOLITON SCATTERING IN (2+1)-D The global phase dierence δ1−δ2 determines the nature of the interaction between the two solitons and ultimately their behavior, switching between particle, wave and liquid-like. As shown in the literature [46], the phase dierence between solitons determines whether the interaction is attractive or repulsive. In the situations considered here, the two solitons and their trajectories are completely symmetric (point reective symmetry relative to the origin), therefore their phase dierence and the character of their interaction remains constant throughout the simulations. Another relevant parameter in the soliton-soliton scattering is the angular momentum, which introduces an eective repulsion between the solitons, and is determined by the initial velocity k and impact parameter b of the solitons. For large impact parameters and large velocities, the interaction is very weak and both solitons almost are not deected from a straight trajectory, being nearly impossible to classify their behavior. Instead, for small impact parameters, the interaction is stronger and a wide variety of soliton behaviors can be observed, as shown in gure (6.2). To help understand the dierent types of phenomena found, we introduce the following classication of soliton scattering: • hard soliton scattering occurs for δ1−δ2=π , where their interaction is mainly repulsive and they interact as if they were rigid spheres for high velocities; • soft soliton scattering occurs for δ1−δ2= 0 and the solitons have distinct behaviors, dominated by mutual attraction. Soft solitons can exhibit a wide variety of behaviors, from liquid light (when they coalesce into a single droplet of light) and planetary-like, when they orbit for a few moments around each other. When the collision velocity is very large, solitons can also destroy each other giving rise to radiation. The designation of hard and soft light introduced here follows a trend found in the literature and originated by several authors, including Michinel et. al [64] who introduced the concept of liquid light. In the following sections we describe in more detail the behavior identied in the simulations. 6.4 Hard soliton scattering For δ1−δ2=π , the two solitons are set to be out-of-phase and their interaction is strongly repulsive, as seen in gure (6.2) in sequences (a-j). This is conrmed in gure 82 6.4. HARD SOLITON SCATTERING Figure 6.2: Typical numerical results for the evolution of two colliding solitons. Sequence a)-e) shows a collision between two out-of-phase solitons with k1= 0.2 and b= 0 , sequence f)-j) displays the results for b= 4 and k)-o) for b= 9 . Sequence p)-t) displays the coalescence of two colliding in-phase solitons with b= 0 and k= 0.3 . Sequence u)-y) shows the results for b= 5 and k= 0.3 . Sequence z)-dd) displays the destruction of two colliding in-phase solitons with b= 0 and k= 0.8 . 83 CHAPTER 6. CASE STUDY 2: SOLITON-SOLITON SCATTERING IN (2+1)-D (6.3), which shows the dependence of the scattering angle θ (measured as a deection of the original straight trajectory) on the impact parameter. Figure 6.3: Computational results for the relation between the scattering angle and the impact parameter, θ(b) . This gure displays the results for colliding outof-phase solitons for k= 0.05 (full line with markers), k= 0.2 (dashed line with markers) and k= 0.3 (pointed line with markers). The full line without markers shows the hard-sphere limit for a sphere with radius of 6 . The increase of collision velocity further strengths the repulsive nature of the interaction but also forces the solitons to come closer to each other. As a result, in this limit the scattering angle approaches the results predicted by the scattering model for hard-spheres [45], revealing the particle-like behavior of solitons. 6.5 Soft soliton scattering For δ1−δ2= 0 , the two solitons are set to be in-phase and their interaction is attractive, as seen in gure (6.2) in sequences (k-y). However the anticipated particle-like dynamics does not hold for small impact parameter as the light pulses tend to coalesce, revealing a liquid-like behavior of solitons. In gure (6.4) it is seen that with the increase of the impact parameter the system undergoes a transition to a particle-like dynamics. For large impact parameters it is seen that the solitons almost do not interact, and thus the scattering angle is approximately null. With the decreasing of the impact parameter we see the expected increase of θ . The peak in gure (6.4) seems to be related with the formation of a quasi-stable state with angular momentum. Unfortunately, the complete analysis of these processes is dicult since using the 84 6.5. SOFT SOLITON SCATTERING Figure 6.4: Computational results for the relation between the scattering angle and the impact parameter, θ(b) . This gure shows the results for colliding in-phase solitons with k= 0.3 (circles), k= 0.25 (crosses), k= 0.23 (triangles) and k= 0.2 (squares). Dierent behaviors are obtained depending on both the collision velocity and the impact parameter. Shaded region is the zone of coalescence for k= 0.3 , where solitons reveal liquid-like behavior. basic algorithm of GASE it is impossible to separate the light intensity pertaining to each soliton. Only the total intensity in each point can be calculated. To overcome this limitation, we developed a version of GASE algorithm where the original GNLSE is decomposed into two equations, say i∂ψ1 ∂z +1 2∇2 ⊥ψ1+F(|ψ|2)ψ1= 0, (6.7) i∂ψ2 ∂z +1 2∇2 ⊥ψ2+F(|ψ|2)ψ2= 0, (6.8) each describing the evolution of the eld associated with each soliton. Notice that the value of F is computed from ψ=ψ1+ψ2 corresponding to the sum of the elds of both solitons. It is trivial to show that if ψ1 and ψ2 satisfy the equations (6.7) and (6.8) respectively then, ψ=ψ1+ψ2 satises the original GNLSE (6.1). This upgrade to the original GASE code allows to identify the evolution of each eld. Figure (6.5-6.8) shows the results obtained for the scattering of two solitons with planetary-like trajectories. The results show a transfer of intensity between these two solitons. This suggests that during the interaction, the two solitons exchange energy 85 CHAPTER 7. CONCLUSIONS computer clusters is dierent, being necessary to re-adapt the numerical tools to work in GPU. In particular, it must be taken into account in which type of memory the data is saved in other to avoid possible memory bottlenecks. On the other hand, the GPU computing belongs to a more general computing paradigm: the heterogeneous computing. The heterogeneous computing aims to use every resources of a computer (CPU, GPU and others) to do massive parallel computations. In the GASE, this concept was used to improve even more the performance of the solver by allocating dierent jobs to the CPU and the GPU, having in mind the times associated to the data transfer between them. This way, all of the resources of the machine are used simultaneously during the process and yields speedups in excess of 100 when comparing to simulations run using a top-of-the-line CPU. With the future improvement of these technologies it is expected that this optimal value of GASE could be even larger. Besides the performance results of GASE, we should also underline its versatility in the analysis of physical problems described by the GNLSE. The code is capable of addressing simulation boxes with 1, 2 and 3 dimensions and the number of sampling points of the eld in excess of 223 , which allows reasonably large simulation domains. Also, several types of nonlinearities can be used simultaneously. These include not only quadratic, cubic, quintic, logarithmic, etc. nonlinearities but also spatial distributions (including the dependence on the propagation coordinate) of the nonlinear refraction index - called nonlinear lattices. Basically, any type of nonlinearity could be simulated. We have also included three types of boundary conditions, namely periodic, reective or absorbing boundary conditions. Finally, in the last stage of the development, we included the possibility of simulating several elds satisfying their own GNLSE but coupled to each other. The development of the code GASE included extensive testing to validate the numerical accuracy of the results and to benchmark the code in relation to other numerical approaches: Crank-Nicholson GPU-based solver and CPU SSFM. In all these cases, GASE was superior in performances. We concluded this dissertation by presenting the application of GASE to the study of two physical problems. Although the results obtained in both problems correspond to a preliminary study, they illustrate the versatility and potential of GASE and GPU computing. Had we more time and we could have explored both these problems in more depth, as well as investigate others. In the following section we give a brief overview of future work that could be developed with the help of GASE. 92 7.1. FUTURE WORK 7.1 Future work At the end of this dissertation, there is a multitude of extensions to this work that could be done. From the point of view of the solver, GASE can be extended to being capable of addressing the problems with many coupled GNLSE, a feature that can be useful for problems where the eld has a vector character. As novices in GPU computing, we believe that there are still certain optimizations that could be done and lend to an even faster solver, such as asynchronous memory transfer from the GPU to the CPU or others. The extension of GASE to work in other GPUs from other manufacturers rather than NVIDIA can also be done, passing the developed CUDA code to the OPEN CL language. The extension of GASE to multi-GPU architectures is also another point of interest. From the point of view of the study of the GNLSE, there is still plenty of room for further developments. In fact, as GASE is capable of addressing almost any problem described by the GNLSE, the solver is ready to tackle a multitude of research problems in the subject of solitons, either in multidimensional systems or in the case of optical lattices. As the GNLSE has a certain character of universality, being present in many elds of physics, the same solver could also constitute a tool for studying phenomena in Bose-Einstein condensates, plasma physics and uid dynamics, among many others. The solver can also be proven useful in other areas of physics, such as general relativity, where it can be used to analyze the propagation of gravitational waves, a high nonlinear type of wave. It is our hope that in the short term GASE will become an important tool of research in these problems. 7.2 Publications Conference proceedings • Nuno A. Silva, M. I. Carvalho, A. Guerreiro, Lightons: Phonons with light, Proceedings of RIAO/OPTILAS 2013, Porto, Portugal, 2013 • Nuno A. Silva, M. I. Carvalho, A. Guerreiro, Spatial soliton dynamics in cubicquintic media, Proceedings of RIAO/OPTILAS 2013,Porto, Portugal, 2013 Poster presentations • Nuno A. Silva, M. I. Carvalho, A. Guerreiro, Lightons: Phonons with light, NANOPT 2013, Porto, Portugal, 2013 93 CHAPTER 7. CONCLUSIONS • Nuno A. Silva, M. I. Carvalho, A. Guerreiro, Lightons: Phonons with light, RIAO/OPTILAS 2013, Porto, Portugal, 2013 • Nuno A. Silva, M. I. Carvalho, A. Guerreiro, Spatial soliton dynamics in cubicquintic media, RIAO/OPTILAS 2013, Porto, Portugal, 2013 94 Bibliography [1] Soliton wave receives crowd of admirers. Nature , 376:373, August 1995. [2] G. Agrawal. Nonlinear Fiber Optics . Academic Press, 3 edition, January 2001. [3] Adrian Alexandrescu, Jose R. Salgueiro, Victor M. Perez-Garcia, and Humberto Michinel. Nonperturbative vector solitary waves in four-level coherent media. Phys. Rev. A , 79:063843, Jun 2009. [4] D. Anderson. Variational approach to nonlinear pulse propagation in optical bers. Phys. Rev. A , 27:31353145, Jun 1983. [5] David W. Aossey, Steven R. Skinner, Jamie L. Cooney, James E. Williams, Matthew T. Gavin, David R. Andersen, and Karl E. Lonngren. Properties of soliton-soliton collisions. Phys. Rev. A , 45:26062610, Feb 1992. [6] Naruyoshi Asano, Tosiya Taniuti, and Nobuo Yajima. Perturbation method for a nonlinear wave modulation. ii. Journal of Mathematical Physics , 10(11):20202024, 1969. [7] A. Balevic, L. Rockstroh, A. Tausendfreund, S. Patzelt, G. Goch, and S. Simon. Accelerating simulations of light scattering based on nite-dierence time-domain method with general purpose gpus. In Computational Science and Engineering, 2008. CSE '08. 11th IEEE International Conference on , pages 327334, 2008. [8] G. P. Berman and F. M. Izrailev. The Fermi-Pasta-Ulam problem: Fifty years of progress. Chaos: An Interdisciplinary Journal of Nonlinear Science , 15(1):015104, 2005. [9] Anjan Biswas. Adiabatic dynamics of non-Kerr law solitons. Applied Mathematics and Computation , 151(1):41  52, 2004. [10] Alexander V. Buryak, Yuri S. Kivshar, Ming-feng Shih, and Mordechai Segev. Induced coherence and stable soliton spiraling. Phys. Rev. Lett. , 82:8184, Jan 1999. 95 Bibliography [11] Ron Caplan. Study of vortex ring dynamics in the nonlinear Schrodinger equation utilizing gpu-accelerated high-order compact numerical integrators, 2003. [12] Ronald M. Caplan. Nlsemagic: Nonlinear Schrodinger. CoRR , abs/1203.1263, 2012. [13] Lapo Casetti, Monica Cerruti-Sola, Marco Pettini, and E. G. D. Cohen. The FermiPasta-Ulam problem revisited: Stochasticity thresholds in nonlinear hamiltonian systems. Phys. Rev. E , 55:65666574, Jun 1997. [14] R. Y. Chiao, E. Garmire, and C. H. Townes. Self-trapping of optical beams. Phys. Rev. Lett. , 13:479482, Oct 1964. [15] NVIDIA Corporation. Popular gpu accelerated applications, 2012. [16] A S Davydov. Solitons in molecular systems. Physica Scripta , 20(3-4):387, 1979. [17] A S Desyatnikov and Andrei I Maimistov. Conservation of the angular momentum for multidimensional optical solitons. Quantum Electronics , 30(11):1009, 2000. [18] Kristian Dimitrevski, Erik Reimhult, Erik Svensson, Anders Ohgren, Dan Anderson, Anders Berntson, Mitek Lisak, and Manuel L. Quiroga-Teixeiro. Analysis of stable self-trapping of laser beams in cubic-quintic nonlinear media. Physics Letters A , 248(5-6):369376, 1998. [19] R. C. Diprima, W. Eckhaus, and L. A. Segel. Non-linear wave-number interaction in near-critical two-dimensional ows. Journal of Fluid Mechanics , 49:705744, 10 1971. [20] Sergey V. Dmitriev, Yuri S. Kivshar, and Takeshi Shigenari. Fractal structures in multi-soliton collisions. Physica B: Condensed Matter , 316-317(0):139142, 2002. Proceedings of the 10th International Conference on Phonon Scattering in Condensed Matter. [21] Galen C. Duree, John L. Shultz, Gregory J. Salamo, Mordechai Segev, Amnon Yariv, Bruno Crosignani, Paolo Di Porto, Edward J. Sharp, and Ratnakar R. Neurgaonkar. Observation of self-trapping of an optical beam due to the photorefractive eect. Phys. Rev. Lett. , 71:533536, Jul 1993. [22] J.; Ulam S. Fermi, E.; Pasta. Studies of nonlinear problems. Technical report, 1955. 96 Bibliography [23] Jason W. Fleischer, Mordechai Segev, Nikolaos K. Efremidis, and Demetrios N. Christodoulides. Observation of two-dimensional discrete solitons in optically induced nonlinear photonic lattices. Nature , 422(6928):147150, March 2003. [24] Michael Fleischhauer, Atac Imamoglu, and Jonathan P. Marangos. Electromagnetically induced transparency: Optics in coherent media. Rev. Mod. Phys. , 77:633673, Jul 2005. [25] Victor M. Perez Garcia and Xiao yan Liu. Numerical methods for the simulation of trapped nonlinear Schrodinger systems. Applied Mathematics and Computation , 144(2-3):215235, 2003. [26] Cliord S. Gardner, John M. Greene, Martin D. Kruskal, and Robert M. Miura. Method for solving the Korteweg-deVries equation. Phys. Rev. Lett. , 19:10951097, Nov 1967. [27] E.P. Gross. Structure of a quantized vortex in boson systems. Il Nuovo Cimento Series 10 , 20(3):454477, 1961. [28] Chao Hang and V. V. Konotop. Spatial solitons in a three-level atomic medium supported by a laguerre-gaussian control beam. Phys. Rev. A , 83:053845, May 2011. [29] Mark Jason Harris. Real-time cloud simulation and rendering, 2003. [30] Akira Hasegawa and Frederick Tappert. Transmission of stationary nonlinear optical pulses in dispersive dielectric bers. i. anomalous dispersion. Applied Physics Letters , 23(3):142144, 1973. [31] Rainer W. Hasse. Schrödinger solitons and kinks behave like newtonian particles. Phys. Rev. A , 25:583584, Jan 1982. [32] Hermann A. Haus and William S. Wong. Solitons in optical communications. Rev. Mod. Phys. , 68:423444, Apr 1996. [33] Simon Huang, Peng Zhang, Xiaosheng Wang, and Zhigang Chen. Observation of soliton interaction and planetlike orbiting in Bessel-like photonic lattices. Opt. Lett. , 35(13):22842286, Jul 2010. [34] A.A. Kanashov and A.M. Rubenchik. On diraction and dispersion eect on three wave interaction. Physica D: Nonlinear Phenomena , 4(1):122  134, 1981. 97 Bibliography [35] V.I. Karpman and V.V. Solov'ev. A perturbational approach to the two-soliton systems. Physica D: Nonlinear Phenomena , 3(3):487  502, 1981. [36] Y. V. Kartashov, V. A. Vysloukh, and L. Torner. Solitons in complex optical lattices. The European Physical Journal Special Topics , 173(1):87105, 2009. [37] Yaroslav V. Kartashov, Boris A. Malomed, and Lluis Torner. Solitons in nonlinear lattices. Rev. Mod. Phys. , 83:247305, Apr 2011. [38] Yaroslav V. Kartashov, Victor A. Vysloukh, and Lluis Torner. Rotary solitons in Bessel optical lattices. Phys. Rev. Lett. , 93:093904, Aug 2004. [39] Yaroslav V. Kartashov, Victor A. Vysloukh, and Lluis Torner. Stable ring-prole vortex solitons in Bessel optical lattices. Phys. Rev. Lett. , 94:043902, Feb 2005. [40] Yaroslav V. Kartashov, Victor A. Vysloukh, and Lluis Torner. Soliton control in fading optical lattices. Opt. Lett. , 31(14):21812183, Jul 2006. [41] Yaroslav V. Kartashov, Victor A. Vysloukh, and Lluis Torner. Brownian soliton motion. Phys. Rev. A , 77:051802, May 2008. [42] Yaroslav V. Kartashov, Victor A. Vysloukh, and Lluis Torner. Soliton modes, stability, and drift in optical lattices with spatially modulated nonlinearity. Opt. Lett. , 33(15):17471749, Aug 2008. [43] P.G. Kevrekidis, D.J. Frantzeskakis, and R. Carretero-González. Emergent Nonlinear Phenomena in Bose-Einstein Condensates: Theory and Experiment . Atomic, Optical, and Plasma Physics. Springer, 2007. [44] Ali Khajeh-Saeed and J. Blair Perot. Direct numerical simulation of turbulence using {GPU} accelerated supercomputers. Journal of Computational Physics , 235(0):241  257, 2013. [45] Tom W. B. Kibble and Frank H. Berkshirel. Classical Mechanics . Imperial College Press, London, fth edition, 2004. [46] Yuri S. Kivshar, Govind P. Agrawal, and Govind P. Agrawal. Optical Solitons . Academic Press, Burlington, 2004. [47] Yukihiro Komura and Yutaka Okabe. Gpu-based single-cluster algorithm for the simulation of the ising model. Journal of Computational Physics , 231(4):1209  1215, 2012. 98 Bibliography [48] D. J. Korteweg and G. de Vries. Xli. on the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves. Philosophical Magazine Series 5 , 39(240):422443, 1895. [49] Andrei I Maimistov. Solitons in nonlinear optics. Quantum Electronics , 40(9):756, 2010. [50] Andrey Maimistov, Boris Malomed, and Anton Desyatnikov. A potential of incoherent attraction between multidimensional solitons. Physics Letters A , 254(34):179184, 1999. [51] Boris A. Malomed. Bound solitons in the nonlinear schrödingerginzburg-landau equation. Phys. Rev. A , 44:69546957, Nov 1991. [52] Boris A. Malomed. Potential of interaction between twoand three-dimensional solitons. Phys. Rev. E , 58:79287933, Dec 1998. [53] Boris A Malomed, Dumitru Mihalache, Frank Wise, and Lluis Torner. Spatiotemporal optical solitons. Journal of Optics B: Quantum and Semiclassical Optics , 7(5):R53, 2005. [54] S. Maneuf, R. Desailly, and C. Froehly. Stable self-trapping of laser beams: Observation in a nonlinear planar waveguide. Optics Communications , 65(3):193  198, 1988. [55] Robert McLeod, Kelvin Wagner, and Steve Blair. (3+1)-dimensional optical soliton dragging logic. Phys. Rev. A , 52:32543278, Oct 1995. [56] Robert McLeod, Kelvin Wagner, and Steve Blair. Variational approach to orthogonally-polarized optical soliton interaction with cubic and quintic nonlinearities. Physica Scripta , 59(5):365, 1999. [57] David Michea and Dimitri Komatitsch. Accelerating a three-dimensional nitedierence wave propagation code using gpu graphics cards. Geophysical Journal International , 182(1):389402, 2010. [58] H. Michinel, J. Campo-Taboas, R. Garcia-Fernandez, J. R. Salgueiro, and M. L. Quiroga-Teixeiro. Liquid light condensates. Phys. Rev. E , 65:066604, Jun 2002. [59] Humberto Michinel, Maria J. Paz-Alonso, and Victor M. Perez-Garcia. Turning light into a liquid via atomic coherence. Phys. Rev. Lett. , 96:023903, Jan 2006. 99 Bibliography [60] Humberto Michinel, Jose R. Salgueiro, and Maria J. Paz-Alonso. Square vortex solitons with a large angular momentum. Phys. Rev. E , 70:066605, Dec 2004. [61] L. F. Mollenauer, R. H. Stolen, and J. P. Gordon. Experimental observation of picosecond pulse narrowing and solitons in optical bers. Phys. Rev. Lett. , 45:1095 1098, Sep 1980. [62] Gaspar D. Montesinos and Victor M. Perez-Garcia. Numerical studies of stabilized townes solitons. Math. Comput. Simul. , 69(5):447456, August 2005. [63] David Novoa, Boris A. Malomed, Humberto Michinel, and Victor M. Perez-Garcia. Supersolitons: Solitonic excitations in atomic soliton chains. Phys. Rev. Lett. , 101:144101, Sep 2008. [64] David Novoa, Humberto Michinel, and Daniele Tommasini. Pressure, surface tension, and dripping of self-trapped laser beams. Phys. Rev. Lett. , 103:023903, Jul 2009. [65] Carolyn L. Phillips, Joshua A. Anderson, and Sharon C. Glotzer. Pseudo-random number generation for brownian dynamics and dissipative particle dynamics simulations on {GPU} devices. Journal of Computational Physics , 230(19):7191  7201, 2011. [66] Tobias Preis, Peter Virnau, Wolfgang Paul, and Johannes J. Schneider. {GPU} accelerated monte carlo simulation of the 2d and 3d ising model. Journal of Computational Physics , 228(12):4468  4477, 2009. [67] M. L. Quiroga-Teixeiro, A. Berntson, and H. Michinel. Internal dynamics of nonlinear beams in their ground states: shortand long-lived excitation. J. Opt. Soc. Am. B , 16(10):16971704, Oct 1999. [68] F. Reynaud and A. Barthelemy. Optically controlled interaction between two fundamental soliton beams. EPL (Europhysics Letters) , 12(5):401, 1990. [69] J. Scott Russel. Report on waves. Technical report, 1844. [70] M. Salsi, O. Bertran-Pardo, J. Renaudier, W. Idler, H. Mardoyan, P. Tran, G. Charlet, and S. Bigo. Wdm 200gb/s single-carrier pdm-qpsk transmission over 12,000km. In Optical Communication (ECOC), 2011 37th European Conference and Exhibition on , pages 13, 2011. 100 Bibliography [71] A. Scott. Encyclopedia of Nonlinear Science . Routledge and Sons, 2005. [72] Wen-Rui Shan, Feng-Hua Qi, Rui Guo, Yu-Shan Xue, Pan Wang, and Bo Tian. Conservation laws and solitons for the coupled cubic quintic nonlinear Schrodinger equations in nonlinear optics. Physica Scripta , 85(1):015002, 2012. [73] Ming Feng Shih and Mordechai Segev. Incoherent collisions between twodimensional bright steady-state photorefractive spatial screening solitons. Opt. Lett. , 21(19):15381540, Oct 1996. [74] Ming-feng Shih, Mordechai Segev, and Greg Salamo. Three-dimensional spiraling of interacting spatial solitons. Phys. Rev. Lett. , 78:25512554, Mar 1997. [75] Mário G. Silveirinha. Eective medium response of metallic nanowire arrays with a Kerr-type dielectric host. Phys. Rev. B , 87:165127, Apr 2013. [76] K. Stewartson and J. T. Stuart. A non-linear instability theory for a wave system in plane poiseuille ow. Journal of Fluid Mechanics , 48:529545, 8 1971. [77] C. Sulem and P.L. Sulem. The Nonlinear Schrodinger Equation: Self-Focusing and Wave Collapse . Number v. 139 in Applied Mathematical Sciences. Springer, 1999. [78] Qingqing Sun, Yuri V. Rostovtsev, and M. Suhail Zubairy. Optical beam steering based on electromagnetically induced transparency. Phys. Rev. A , 74:033819, Sep 2006. [79] Zhi-Yuan Sun, Yi-Tian Gao, Xin Yu, Wen-Jun Liu, and Ying Liu. Bound vector solitons and soliton complexes for the coupled nonlinear schrödinger equations. Phys. Rev. E , 80:066608, Dec 2009. [80] Morikazu Toda. Vibration of a chain with nonlinear interaction. Journal of the Physical Society of Japan , 22(2):431436, 1967. [81] I.M. Uzunov, V.S. Gerdjikov, M. Golles, and F. Lederer. On the description of nsoliton interaction in optical bers. Optics Communications , 125(4):237242, 1996. [82] Sibo Wang, Junbo Xu, and Hao Wen. Accelerating dissipative particle dynamics with multiple {GPUs}. Computer Physics Communications , 184(11):2454  2461, 2013. 101