GPU-accelerated time integration of Gross-Pitaevskii equation with discrete exterior calculus
Full text
This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ GPU-accelerated time integration of Gross-Pitaevskii equation with discrete exterior calculus © 2022 The Author(s). Published version Kivioja, Markus; Mönkölä, Sanna; Rossi, Tuomo Kivioja, M., Mönkölä, S., & Rossi, T. (2022). GPU-accelerated time integration of Gross-Pitaevskii equation with discrete exterior calculus. Computer Physics Communications, 278, Article 108427. https://doi.org/10.1016/j.cpc.2022.108427 2022
Computer Physics Communications 278 (2022) 108427 Contents lists available at ScienceDirect Computer Physics Communications www.elsevier.com/locate/cpc GPU-accelerated time integration of Gross-Pitaevskii equation with discrete exterior calculus✩ Markus Kivioja∗, Sanna Mönkölä, Tuomo Rossi University of Jyvaskyla, Faculty of Information Technology, P.O. Box 35, FI-40014 University of Jyvaskyla, Finland a r t i c l e i n f o a b s t r a c t Article history: Received 18 November 2021 Received in revised form 5 May 2022 Accepted 19 May 2022 Available online 24 May 2022 Keywords: Gross-Pitaevskii Discrete exterior calculus GPGPU The quantized vortices in superfluids are modeled by the Gross-Pitaevskii equation whose numerical time integration is instrumental in the physics studies of such systems. In this paper, we present a reliable numerical method and its efficient GPU-accelerated implementation for the time integration of the three-dimensional Gross-Pitaevskii equation. The method is based on discrete exterior calculus which allows us the usage of more versatile spatial discretization than traditional finite difference and spectral methods are applicable to. We discretize the problem using six different natural crystal structures and observe the correct choices of spatial tiling to decrease the truncation error and increase the reliability compared to Cartesian grids. We pay attention to the computational performance optimizations of the GPU implementation and measure speedups of up to 152-fold when compared to a reference CPU implementation. We parallelize the implementation further to multiple GPUs and show that 92% of the computation time can fully utilize the additional resources. ©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction Gross [1] and Pitaevskii [2]introduced a mathematical model for the physics of quantized vortices in superfluids. The model is expressed as the Gross-Pitaevskii equation (GPE), also known as the non-linear Schrödinger equation, which in the case of dilute gases of bosonic atoms [3]reads i¯ h∂t(r,t)=−¯ h2 2m∇2+V(r)+g|(r,t)|2(r,t), (1) where is the complex-valued wave function, iis the imaginary unit, ¯ his the reduced Planck constant, mis the atom mass, Vis the external potential, and gis the effective interaction strength. The wave function is normalized so that |(r, t)|2d3r =N, where Nis the number of atoms. The numerical time integration of the GPE is of great interest in physics e.g. for simulating and studying the behavior of Bose-Einstein condensates as shown by recent studies [4–6]. The numerical solutions require the spatial domain to be discretized and a comprehensive survey of GPE time integration methods [7] shows that it is widely done by using evenly spaced square and cubic grids, in 2D and 3D domains, respectively. In addition to the ✩The review of this paper was arranged by Prof. Hazel Andrew. *Corresponding author. E-mail addresses: markus.i.ki[email protected].fi (M. Kivioja), sanna.monkola@jyu.fi (S. Mönkölä), tuomo.j.rossi@jyu.fi (T. Rossi). finite difference method, recently, the two mostly used discretization methods for the GPE have been the Fourier spectral method [8,9] and the finite element method which has enabled the usage of axial symmetric [10] and triangle [11]grids as well. It has been observed however that in the case of Maxwell’s equations the choice of more diverse discretization grids can improve the accuracy of the solution [12,13]. In this paper, we will show that this is the case with the three-dimensional GPE as well, and that the accuracy of the time integration and the reliability of numerical simulations can be improved by increasing the isotropy of the spatial discretization. We measure the error of the time integration by using six different natural crystal structures for the discretization and discover the cubic grid to produce the largest error of them all. We enable the usage of more sophisticated grids than what has previously been used with the GPE, by applying the concepts of differential geometry and exterior calculus. Particularly the employment of differential forms plays an important role by removing the metric from the differential operators and allowing their exact presentation also at the discrete level. The discrete extension to differential forms is given by the discrete exterior calculus (DEC) [14,15] which we will concentrate on in this paper. In DEC the only source of discretization error is the discrete Hodge star operator, which defines the relations between primal and dual discrete differential forms, and thus is determined by the relation of the primal and dual spatial grids. DEC has previously been successfully utilized in elastodynamics [16], fluid dynamics [17], electromaghttps://doi.org/10.1016/j.cpc.2022.108427 0010-4655/©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
M. Kivioja, S. Mönkölä and T. Rossi Computer Physics Communications 278 (2022) 108427 netism [18], and in quantum mechanics in the form of the standard linear Schrödinger equation [19]. Our work is the first one to apply DEC to the non-linear Schrödinger equation. The computational performance of the time integration is of high importance, as the computation time determines how fine the discretization, and how large the spatial and time domains can be in practice. Those factors may determine the set of physics research problems the time integrator is applicable to, as the physical properties of the simulated systems can compose requirements for the accuracy and domain size [20]. The independent nature of a time integration step of an individual grid vertex, makes the method ideal for the massively parallel computing model [21]. In [22], [23], [24], and [25], general-purpose GPU computing (GPGPU) [26] has been utilized for solving the non-linear Schrödinger equation with promising results. Though this has only been done using Cartesian spatial discretization and in many cases by utilizing a ready-made linear algebra library. We implement our DEC-based method on multiple GPUs by giving attention to the lower-level details and common performance bottlenecks [27]in the modern GPU hardware architectures. We note that the memory utilization and access patterns in particular play an important role in the performance optimizations of the GPU implementation, and we concentrate on them especially. The implementation shows up to 152-fold performance improvement when benchmarked against a reference CPU version, and good scalability over the number of GPUs. There is no previous research on the topic of using GPUs in DEC. Esqueda et al. [28]mention the usage of a single GPU in their numerical experiments, but give no description of the implementation nor its computing time performance, therefore this paper aims to fill those gaps. The rest of the paper is organized as follows. In section 2we describe the method for discretizing the spatial domain of the time integration and the six different grid types used in our experiments. Section 3presents the application of the concepts of discrete exterior calculus to the GPE. We give a detailed description of our implementation and performance optimizations of the method on multiple GPUs in section 4. The numerical behavior and accuracy, as well as the execution time performance and its scalability, are examined in section 5. Last, we summarize the paper in section 6. 2. Spatial discretization The three-dimensional spatial domain is discretized by using cell complexes consisting of linear and convex oriented k-cells, where k ∈{0, 1, 2, 3}. Hence in our case, 0-cells are defined as points x ∈R3, oriented 1-cells as line segments between two 0cells, oriented 2-cells as convex two-dimensional linear surfaces enclosed by finite sets of 1-cells, and oriented 3-cells as convex three-dimensional volumes enclosed by finite sets of 2-cells. The (k −1)-cells enclosing a k-cell are called boundary cells and the orientations of k-cells are determined by the order of their boundary cells. We use sets of two cell complexes which consist of a primary complex and a dual complex. The dual complex is constructed by assigning a dual (3 −k)-cell for each primary k-cell and by setting the positions of the dual 0-cells to be the circumcenters of their corresponding primary 3-cells, which ensures that the dual 1-cells are always orthogonal to the primary 2-cells and vice versa. We choose the circumcentric positioning of the dual cells, since it doesn’t yield significant losses in the accuracy [29], but is computing timewise an optimal option, as we will show in section 3. Since the quantized vortices modeled by the GP equation can advance in all directions, we make a hypothesis that the best accuracy for the time integration is achieved by maximizing the spatial isotropy of the cell complex. For this purpose, we construct the Fig. 1. The grid types used for the spatial discretization [13,32]. cell complexes using the six different natural crystal structures that were applied to the time integration of Maxwell’s equations by Räbinä [30]. These structured grids are the cubic crystal systems: cubic, face-centered cubic (FCC), and body-centered cubic (BCC), and the tetrahedrally close-packed (TCP) structures: A15, C15, and Z. The standard cubic tiling acts as the basis for the other cubic crystal systems. The FCC and BCC grids are constructed by adding vertices at the center of the faces and bodies of the cubic grid, respectively. The FCC grid-based tiling of the spatial domain consists of alternating bodies of regular octrahedra and tetrahedra, and the dual bodies are Kepler’s rhombic dodecahedra. In the BCC grid tiling, the primary bodies are congruent tetrahedra, whose faces are isosceles triangles with the relation of 2 :√3 between the bottom and side edges. The dual bodies of the BCC tiling are truncated octahedra with equal length edges. The motivation for the TCP structures is to achieve small dihedral angles, which is shown to be a preferred quality of a grid in the time integration of Maxwell’s equations [31]. The A15 and C15 grids are constructed by adding vertices to the BCC and FCC grids, respectively. The dual bodies of the A15 grid are irregular dodecahedra, centered at the BCC grid vertices, and tetrakaidecahedra surrounding the other primal vertices. The dual C15 grid consists of 12-hedra and 16-hedra. Unlike the other structures, the Z grid is asymmetric in that its z-direction is divergent from the x-, and y-directions. Though, the xy-plane is symmetric in 60-degree increments. Its dual grid is composed of 12-hedra, 14-hedra, and 15-hedra. The BCC and C15 grids are expected to perform well since C15 has the smallest dihedral angle and BCC least variance in the edge lengths. Additionally, both grid types are previously shown to be the most numerically isotropic [32]. All of the aforementioned cubic crystal systems and TCP structures are visualized in Fig. 1. With each of the grid types, we construct the full spatial tiling by repeating a basis replicable structure, and during the rest of this paper we call their separate occurrences replicable structure instances. 3. Computational model We consider the more general dimensionless form of the timedependent GPE (see [33]) i∂¯ t¯ (¯ r,¯ t)=−1 2∇2+¯ V(¯ r)+¯ g|¯ (¯ r,¯ t)|2¯ (¯ r,¯ t), (2) where ¯ , ¯ V, ¯ g, ¯ r, and ¯ tare suitably scaled versions of their barless counterparts, and |(¯ r, ¯ t)|2d3¯ r=1. In this paper, we concentrate only on the time integration of the GPE and omit the computation of the initial value, in which case the exact scales can be ignored and ¯ Vand ¯ gconsidered as some adjustable real-valued scalar parameters that are constant w.r.t. time. In order to apply discrete exterior calculus, we first transform the equation in smooth setting into a more suitable format, which 2
M. Kivioja, S. Mönkölä and T. Rossi Computer Physics Communications 278 (2022) 108427 assumes familiarity with differential geometry and exterior algebra. We can regard the scalar wave function ¯ as a smooth differential 0-form ˆ , which generalizes the equation to smooth differentiable manifolds of arbitrary dimensions i∂¯ tˆ =1 2δd+¯ V+¯ g|ˆ |2ˆ , where dis the exterior derivative, and δ=− −1dis the codifferential. Here is the Hodge star operator which, in n-dimensional space, maps k-forms to (n −k)-forms and is defined by the relation α∧(β) =α, β e1∧···∧enfor all pairs of k-forms α, β, when e1···enare orthonormal basis forms and ∧the wedge product. We introduce a differential 1-form ˆ uand write the above equation as a pair of equations i∂¯ tˆ =1 2δˆ u+¯ Vˆ +g|ˆ |2ˆ , ˆ u=dˆ . Considering the pair in 3-dimensional space, then by applying the Hodge star operator on both sides of the upper equation, integrating the upper and lower equations over a 3-manifold Vand a 1-manifold C, respectively and applying the generalized Stokes’ theorem dω= ∂ ω to both equations, the system takes an integral form of i∂¯ t V ˆ =1 2 ∂V ˆ u+ VV+g|ˆ |2ˆ , C ˆ u= ∂C ˆ , (3) where ∂denotes the boundary of a manifold . Hirani et al. [14]defined the discrete differential k-form on a kcell ckby the de Rham map ˜ α=α,ck= ck α, where αis a smooth differential k-form, and ˜ αits corresponding discrete counterpart. Furthermore the diagonal discrete Hodge star operator was introduced and defined by the relation 1 |∗ck|α,∗ck= 1 |ck|α,ck,(4) where ∗ckis the circumcentric dual cell of ck, and |ck|denotes the k-volume of the cell. These definitions apply naturally to the integral form eq. (3)of the GPE when its smooth manifolds are replaced by oriented complex cells and the smooth Hodge star by the discrete version. More precisely, if Bis an oriented 3-cell and its boundary ∂B consists of oriented 2-cells Fi, and we denote their dual 1-cells by E∗ iand further their boundaries ∂E∗ iby dual 0-cells Ni∗ j, we see that ∂B ˆ u= i Fi ˆ u= iˆ u,Fi= i |Fi| |E∗ i|ˆ u,E∗ i, ˆ u,E∗ i= E∗ i ˆ u= ∂E∗ i ˆ = j Ni∗ j ˆ = jˆ ,Ni∗ j. Then by noting that Bˆ =ˆ , B =|B| |B∗|ˆ , B∗and that the 0-volume |B∗| =1, we get i∂¯ tˆ ,B∗=1 21 |B| i |Fi| |E∗ i| jˆ ,Ni∗ j+[¯ V+g|ˆ |2]ˆ ,B∗.(5) Note that here ˆ , B∗is the value of ¯ at the circumcenter of B, and ˆ , Ni∗ jare the signed values at the positions of Ni∗ j, where the signs are defined by the orientations of E∗ i. When the eq. (5)is applied to every dual 0-cell of a cell complex consisting of m ∈N3-cells, the resulting linear system of equations can be presented in a matrix form of i∂¯ tψ=1 23d2−1 2dT 2+Uψ, (6) where d2∈Zm×mFis a sparse incidence matrix whose elements are defined as (d2)i,j=±1if the jth 2-cell is a boundary cell of the ith 3-cell and 0otherwise. The sign is defined by the orientation of the boundary 2-cell w.r.t. the orientation of the 3-cell. Here mF=m i#∂Bi, where #∂Bidenotes the cardinality of the set of boundary 2-cells of the ith 3-cell. ψ∈Cmis a column vector consisting of the wave function values at the positions of the dual 0-cells, and the diagonal matrices 2∈RmF×mFand 3∈Rm×m defined as 2=diag|E∗ 1| |F1|,|E∗ 2| |F2|,..., |E∗ #∂Bm| |F#∂Bm|, 3=diag1 |B1|,1 |B2|,..., 1 |Bm|. U∈Rm×mis a diagonal matrix with elements Ui,i=¯ Vi+¯ g|ψi|2, where ¯ Videnotes the value of ¯ V(¯ r)at the position of the ith dual 0-cell. The diagonality of the matrices kis a direct consequence of the definition eq. (4) which only applies to circumcentric dual complexes. The truncation error of eq. (6)is fully contained in the matrices kand is numerically shown to be O(|E∗|2)[34]. 3.1. Time discretization The time domain is discretized by using the standard centraldifference method [35] which, when applied to the eq. (6), yields iψ(n+1)−ψ(n−1) 2¯ t=1 23d2−1 2dT 2+U(n)ψ(n), and furthermore gives us the formula for the time integration steps: ψ(n+1)=ψ(n−1)−2i¯ tM(n)ψ(n),(7) where ¯ tis the length of the time step, ψ(n)and U(n)the wave function value vector and its corresponding matrix U, respectively, at a time instant n¯ tso that n ∈Z, and M(n)=1 23d2−1 2dT 2+U(n)∈Rm×m. In a general case, when the edge lengths of the dual grid vary, as is the case with the FCC, A15, C15, and Z grids, the ith row of the eq. (7)is written out as ψ(n+1) i=ψ(n−1) i−i¯ t ×⎡ ⎣ 1 |Bi| j |Fi,j| |E∗ i,j|(ψ(n) i−ψ(n) i,j)+(¯ V+g|ψ(n) i|2)ψ(n) i⎤ ⎦, (8) 3
M. Kivioja, S. Mönkölä and T. Rossi Computer Physics Communications 278 (2022) 108427 where Fi,jdenotes the jth boundary cell of the ith 3-cell and E∗ i,j its dual cell. ψi,jdenotes the element corresponding to the dual 0-cell that is a boundary cell of E∗ i,jbut is not the ith dual 0-cell. The summation is over all boundary 2-cells of the ith 3-cell. The cubic and BCC grids are special cases in that their dual edge lengths and the numbers of boundary faces of dual bodies don’t vary, which allows us to simplify the eq. (8)to ψ(n+1) i=ψ(n−1) i−i¯ t ×⎡ ⎣|F| |E∗||B|(#∂Bψ(n) i− j ψ(n) i,j)+(¯ V+g|ψ(n) i|2)ψ(n) i⎤ ⎦. The time discretization related truncation error of eq. (7)is the well known O(¯ t2)error of the central-difference method. The time discretization method is conditionally stable, in that there is a |E∗|dependent upper limit for the time step size and when the limit is exceeded the method becomes unstable. The instability is generally caused by a real or imaginary part of an element of ψchanging its sign on every time step which will eventually cause the absolute value to increase without a limit as shown in [32]. We try to avoid this by limiting the maximum absolute difference of changes of ψR:=Re(ψ) (or, equivalently, of ψI:=Im(ψ)) on subsequent time steps with |(ψR(n+2) i−ψR(n) i)−(ψR(n) i−ψR(n−2) i)|≤CψRmax i,(9) for some constant C. Here ψRmax is a vector so that |ψR(n) i| ≤ ψRmax ifor all iand n. We take a linear approximation of the eq. (7) by using a matrix M:= M(0)and with that we get (ψ(n+2) R−ψ(n) R)−(ψ(n) R−ψ(n−2) R)=2¯ tMψ(n+1) I−2¯ tMψ(n−1) I =2¯ tM(ψ(n+1) I−ψ(n−1) I)=−4¯ t2M2ψ(n) R, which when applied to eq. (9), with the definition of ψRmax, gives us ¯ ti≤CψRmax i |4M2i,∗ψmax R|, where M2i,∗denotes the ith row of the matrix M2. We note that Mi,i>> Mi,jwhen j =iand hence approximate Mby considering only its diagonal elements which reduces the above condition to ¯ ti≤√C 2Mi,i.(10) For a global ¯ twhich holds for all i, the greatest matrix diagonal element must be picked. Since eq. (10)is an approximation we used numerical experiments to find a sufficient value for the constant C. With the grid types introduced in section 2, the time integration was found to be stable when C=1 which gives us the final condition ¯ t<(2 max 1≤i≤mMi,i)−1.(11) 4. GPU implementation Since the matrix Min the time integration step eq. (7)is usually large, the vector multiplication may be infeasible slow, when computed in a serial manner. We, therefore, take advantage of parallel computing and utilize GPGPU. This is done by implementing the DEC-based GPE time integrator using Nvidia’s CUDA API [36] and targeting hardware with CUDA compute capability of 5.2 and above. In this and the following sections, we will also use the terminology introduced with CUDA. We assign one GPU thread for each dual 0-cell in the cell complex so that the ith thread will effectively compute the value of the expression ψ(n+1) i=ψ(n−1) i−2i¯ tM(n) i,∗ψ(n),(12) where Mi,∗denotes the ith row of the matrix M. Since Mis sparse we precompute, for each row i, a set Ji={j ∈N:Mi,j= 0}of the column indices of the non-zero elements, and only include the columns j ∈Jito the computation of eq. (12). From the derivation of the matrix Min section 3, it follows that #Ji=#∂Bi+1, where Biis the 3-cell enclosing the ith dual 0-cell. Since the full spatial grid is constructed by duplicating a replicable structure, we also generate the sets of the columns of the non-zero matrix elements only for one structure and reuse them with the other structure instances. This is achieved by replacing the sets Jiby defining new sets ˆJi=ˆj∈Z4:ˆj=cj l,sj x−si x,sj y−si y,sj z−si z, Mi,j=0,i= j, where cj lis the local index of the jth dual 0-cell in the replicable structure, and sj x, sj y, sj zdenote the x-, y-, and z-directional indices of the cell’s structure instance, respectively. In other words, we express the columns of the non-zero matrix elements using four-component vectors consisting of the local index in the replicable structure and the x-, y-, and z-directional offsets of the indices of the structure instance w.r.t. the indices of the structure instance of the ith dual 0-cell. With this approach, the number of required different sets ˆJiis decreased to ms, where msdenotes the number of the dual 0-cells in a single replicable structure instance. We know that every set Jiholds an element with a value of i, so we can drop the cases of i =jfrom the sets ˆJiand thereby reduce their cardinality to #ˆJi=#∂Bi. Likewise the matrix Mis compressed by storing only the nonzero elements corresponding to one structure instance so that the stored element count is ms i=0#ˆJi. The matrix M and the sets ˆJi depend on the grid type and are computed only once in advance for each of the grids. Apart from the optimizations mentioned later in the paper, Mand ˆJiare the only grid-dependent components of the implementation and hence the runtime part of the integrator does not require changes between the different spatial discretizations. The solution vector ψis double buffered in the global GPU memory so that one of the buffers holds the values for the time instants t(n)of odd n, and the other for the time instants of even n. The GPU kernel function treats one of the buffers as the vector ψ(n)in eq. (12) and the other one as both ψ(n+1)and ψ(n−1)by updating its value cumulatively. The same kernel is invoked for every time step and the two value buffers are swapped in between the invocations. The GPU kernel function is described in pseudocode in Algorithm 1, where Sis a three component vector, holding the domain size, measured in replicable structure instance counts, Js an array of all unique sets ˆJ, sa dense matrix holding the absolute values of the non-zero elements of 3d2−1 2dT 2corresponding to the indices in Js, and dt the time step size. 4.1. Memory access pattern The thread group dimensions of the kernel launches are configured in a way that the global yand z-directional indices of the threads match with the indices of the replicable structure instances, whose value elements the threads are assigned to compute. The global thread count in the x-direction is multiplied by 4
M. Kivioja, S. Mönkölä and T. Rossi Computer Physics Communications 278 (2022) 108427 Algorithm 1 GPU kernel for computing one time step. 1: procedure TimeStep(ψ(n+1), ψ(n−1), ψ(n), s, Js, V, g, dt) 2: cl←(blockIdx.x×blockDim.x+threadIdx.x)mod ms 3: sx←(blockIdx.x×blockDim.x+threadIdx.x)/ms 4: sy←(blockIdx.y×blockDim.y+threadIdx.y) 5: sz←(blockIdx.z×blockDim.z+threadIdx.z) 6: i ←ms×(sz×S.x×S.y+sy×S.x+sx) +cl 7: J←Js[cl] 8: ←s[cl]Take one row from the matrix s 9: ←(0 +0i) ∈C 10: for ∀f∈{0, 1, ..., #J−1}do 11: j ←J[f]jis a 4 element array 12: j ←ms×(j[3] ×S.x×S.y+j[2] ×S.x+j[1]) +j[0] 13: ← +[f] ×(ψ(n)[i] −ψ(n)[j]) 14: ψ(n+1)[i] ←ψ(n−1)[i] −dt ×i × +(V[i] +g|ψ(n)[i]|2)ψ(n)[i] msso that the corresponding x-directional structure index is determined by sx=tx ms, where tx∈Ndenotes the global x-directional thread index. The local index of the dual 0-cell in the structure is then given by cl=txmod ms. The memory layout of the vectors ψis such that the elements corresponding to the dual 0-cells in the same replicable structure instance are stored in a contiguous interval, and the structure instances themselves in row-major order, where a row is considered to coincide with the x-direction. Therefore by structuring the thread groups in an aforementioned way, we maximize the coalescing of the global memory load and store access patterns. We use double-precision floating-point numbers for all of the data and arithmetics, so a complex-valued element of ψwill occupy 16 bytes of global memory. Hence on our target hardware with L2-cached global memory accesses with 32-byte cache lines, every pair of two subsequent elements can be accessed with the same transaction. Therefore when all of the threads access the element of their assigned dual 0-cells simultaneously, a warp of 32 threads makes 16 transactions per load/store. A set ˆJicontains duplicates of elements of other sets ˆJl, with l = i, in which case every element of ψ(n)will be read multiple (# Ji) times during the computation of one time step. We avoid increasing the number of global memory accesses by utilizing the local on-chip memory on modern GPU hardware, called shared memory. This is accomplished by making every thread first load its corresponding element of ψ(n)from the global memory into the shared memory so that during the computation of ψ(n+1)most of the elements of ψ(n)can be loaded from the shared memory with low latency. With the current GPU hardware architectures the shared memory can be shared only between the threads inside the same thread group, and because some of the needed elements of ψ(n)might be loaded into the shared memory by threads living in a different thread group than the one computing the element of ψ(n+1), those still need to be read from the global memory. The number of elements needed to be read from the global memory during the computation can be reduced by maximizing the size of the thread groups, so that when the threads are associated with the dual 0-cells, the area of the borders between thread groups is minimized. In the GPU hardware the shared memory is divided into banks so that each bank can be accessed simultaneously. However when multiple threads in the same warp make a request to the same bank the accesses conflict and are serialized. We avoid this by properly permuting the elements of the sets ˆJiand by adding memory padding in between the elements of the vectors ψ. We optimize the method for hardware with 32 banks and a bank width of 32 bits. Therefore each double-precision complex-valued element of 128 bits will occupy 4 consecutive banks in the shared memory. Since one bank can only transfer 32 bits per transaction, a load/store of 128 bits per thread will always need at a minimum of four transactions. For this reason, GPUs with CUDA compute capability of 5.2 and above will access 128 bits per thread for 8 threads at a time so that the whole warp access creates exactly four transactions and causes no bank conflicts. In consequence, we can design our shared memory layout and accesses for eliminating bank conflicts by considering an artificial system with a warp size of 8 threads and a shared memory with 8 banks of 128 bits per bank. From this, it follows that the element with the local dual 0-cell index of clis stored in the same bank as all the elements with an index (cl+8k)mod ms, where k ∈N+. Consequently, the number of elements of different dual 0-cell indices stored in the same bank is lcm(ms,8) 8, where lcm(·, ·)denotes the least common multiple. When msis even, the previous expression can have values of ms 2, ms 4, or ms 8. When msis odd every bank would hold elements of all the different dual 0-cell indices, but in that case we add a 128 bit-sized padding in between the last and first elements of different replicable structure instances, which reduces the case back to the even ms. The elements of the sets ˆJiare strove to be permuted so that for any fixed l ∈{0, 1, ···, #ˆJi}, it holds that ˆJi l,0=(ˆJj l,0+8k)mod ms, ∀i,j∈{0,1,···,ms}:i= j,i 32= j 32, (13) where ˆJi l,0denotes the first component of the lth element in the set ˆJi. The equality of the integer parts of the divisions by 32 relaxes the requirement to apply only inside the same warp. When ms>32 there might be no possible permutations for the elements so that the inequality eq. (13)would hold, without also permuting the order of the sets ˆJithemselves. We avoid doing that since it would decrease the number of coalesced accesses to the global memory, which has a greater performance cost. For that reason with some grid types, with large ms, we settle for only minimizing the number of pairs (i, j)breaking the rule. A heuristic iterative algorithm is used for the generation of the permutations. We are able to allow for ˆJi l,0=ˆJj l,0, because of the shared memory broadcasting capability offered by our target GPU hardware, which permits multiple threads in the same warp to access the same 32-bit word without causing bank conflicts. However it is not always the case that the target elements of the loads made by the ith and jth thread are the same, even though ˆJi l,0=ˆJj l,0. This happens when loads are made across the borders of structure instances, and when msmod 32 =0, so that different threads in the same warp may not be assigned with dual 0-cells within the same structure instance. We eliminate the bank conflicts caused by these cases by first squeezing the thread groups to extend only in the x-direction so that one warp is able to access at a maximum of three structure instances, which are also guaranteed to be consecutive. Then by adding memory padding in between all the elements of ψ, if necessary, to increase the number of structure instances in between the two closest elements which have the same local dual 0-cell index and share the same bank, to be greater than two. Fig. 2presents the layout of the shared memory after 64-bit padding has been added, in the case of the BCC and FCC grids with ms=12. A matrix cell in the figure represents a combination of two consecutive 32-bit banks, and the numbers inside the cells are the local dual 0-cell indices corresponding to the elements stored in the banks. From the figure, it can be seen how there are 4ms elements in between the cases where two elements of the same clare stored in the same bank. Without the padding, the interval would be only 2ms, and the cases of ˆJi l,0=ˆJj l,0able to cause bank conflicts. 5
M. Kivioja, S. Mönkölä and T. Rossi Computer Physics Communications 278 (2022) 108427 0246 8 10 12 14 16 18 20 22 24 26 28 30 0 0 1 2 345 1678 9 10 211 0123 . . . 878 9 10 11 9 0 1 2 345 Bank Address / 1024 B Fig. 2. The layout of ψin the shared memory banks after the padding has been added, in the case of the BCC and FCC grids. ngpu −1 ngpu ngpu +1 Fig. 3. The stored elements of ψby different GPUs, and the memory copies between them, presented in the units of z-slices. The gray cells denote the z-slices the GPUs are assigned to compute. 4.2. Multiple GPUs The implementation is further parallelized to multiple GPUs by using domain decomposition [37]. The buffers for the vectors ψ are laid out in the memory in a way that the elements corresponding to the replicable structure instances on the same xy-plane (zslice) will be in a contiguous memory area. Therefore we distribute the work across multiple GPUs by splitting the spatial domain across the z-axis, to minimize the number of data copy operations between different GPUs, and hence the overhead caused by them. The splitting is done as evenly as possible so that the maximum difference in the number of assigned z-slices for different GPUs is one. The computation domains assigned to different GPUs are separate, though two consecutive GPUs store elements of ψin a way that the subdomains covered by the stored elements overlap over two z-slices. Hence if we denote the total number of z-slices with Nz, and the number of GPUs with Ngpu (and we assume that Nzmod Ngpu =0), the GPU of index ngpu ∈Nwill be assigned to compute the elements of the z-slices with indices from ngpu Nz Ngpu to ngpu +1Nz Ngpu −1, inclusive, but store the elements from ngpu Nz Ngpu −1to ngpu +1Nz Ngpu . After each time step the ngputh GPU will copy the elements of its second z-slice to the memory of the last z-slice of the (ngpu −1)th GPU. Similarly the second to last z-slice will be copied to the first zslice of the (igpu +1)th GPU, as shown in Fig. 3. Each array cell in the figure represents a whole z-slice and the cells colored with gray denote the assigned computational domains of the GPUs. To avoid unnecessary synchronizations between the CPU and GPUs, and between and within individual GPUs, we utilize the concepts of streams, events, and asynchronous memory copies and kernel launches introduced by the CUDA API. For every GPU there are three streams and events created, one for the kernel execution and one for each of the two memory copies. The two memory copies of a single GPU are independent of each other, and therefore can be initiated concurrently, which furthermore is achieved by placing them in separate streams. The kernel invocations are placed in their own stream, with an event after each invocation, to signal the memory copy streams when they can start their work. An event is also added to each of the memory copy streams, to signal the kernel streams of the neighboring GPUs of the next time Table 1 The lengths of dual 1-cells after the scales, for equalizing the operation counts, have been applied. Cubic FCC BCC A15 C15 Z 1 1.09 0.80 0.93 0.68 0.99 step. The CPU is not blocked by any of the GPU work, up until the command buffers of the GPUs are fully occupied, or we need to access the time integration results on the CPU side. 5. Numerical experiments We test the method by time integrating the dimensionless GPE eq. (2) and initializing the wave functions using stationary vortex solutions ¯ λ,¯ g,κ, arising from the dynamics of Bose-Einstein condensates [38]. The initial values take a form of ¯ λ,¯ g,κ(¯ r,¯ t)=f(¯ ρ,¯ z)eiκ¯ φ−i¯ μ¯ t, where fis a real-valued function satisfying the time-independent equation 1 2κ2 ¯ ρ2−∂2 ¯ ρ−1 ¯ ρ∂¯ ρ−∂2 ¯ z+¯ V+¯ gf2f=¯ μf.(14) Here ¯ r=(¯ ρ, ¯ φ, ¯ z)is the dimensionless position vector presented in cylindrical coordinates, κthe so-called winding number, and ¯ μ the dimensionless chemical potential. The eq. (14)can be solved with e.g. a relaxation method [39]. An example of an initial value is visualized in Fig. 4. We initialize the vectors of the discretized wave functions with ψ(n) i=¯ λ,¯ g,κ(¯ ri,n¯ t),n∈{−1,0}. The computational domain is the smallest cuboid so that |¯ λ,¯ g,κ(¯ r,0)|>10−5max ¯ r|¯ λ,¯ g,κ(¯ r,0)|,(15) for all ¯ rinside the cuboid. As the boundary condition, we use ¯ (¯ r) =0, for all ¯ routside of the cuboid. All of the source codes used in the measurements of this section are available in a public repository [40]. 5.1. Accuracy We compare the accuracy of different spatial discretizations by scaling the grids in a way that the lengths of the dual 1-cells are less than 5% of the effective wavelength of ¯ λ,¯ g,κ(¯ r, 0), and that the number of arithmetic operations per integration of one time unit m i=0#Ji1 ¯ t, stays constant between the different grid types, and is 6×109. The resulting distances of two adjacent dual 0-cells w.r.t. the cubic grid are presented in Table 1. The accuracies are estimated by comparing the time integration results to a stationary vortex state ¯ λ,¯ g,κ. We use a stationary vortex state with λ =1, ¯ g=300, κ=10, that can be shown by the Bogoliubov equations (see [41]), to be dynamically unstable. This causes the vortex to split up into multiple smaller parts, by even a small perturbation to the stationary state. Hence, we first measure the amount of time ¯ tthe equation can be integrated before the perturbation caused by the numerical error initiates the splitting. We do this by measuring an error defined as E(¯ t) =1 −| ¯ ∗ λ,¯ g,κ(¯ r, 0)¯ (¯ r, ¯ t)d3¯ r|, and considering the split to happen when E(¯ t) >0.01. From Fig. 5a it can be perceived that the simulated vortex stays intact for the longest when the BCC and C15 grids are used. The cubic grid not only causes the quickest vortex splitting but also 6
M. Kivioja, S. Mönkölä and T. Rossi Computer Physics Communications 278 (2022) 108427 Fig. 4. A visualization of an initial value ¯ λ,¯ g,κ(¯ r, 0), with λ =0.1, ¯ g=5000, and κ=10. Three cross sections are highlighted at z=−10, 0, 10. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) 0 20406080100 10−6 10−5 10−4 10−3 10−2 ¯ t Cubic A15 Z FCC BCC C15 (a) E(¯ t) 0204060 0 1 2 3 4·10−2 ¯ t (b) RMSE(¯ t) Fig. 5. Errors E(¯ t)and RMSE(¯ t)produced by the numerical time integration with different spatial discretizations. results in increasing deviation from the stationary state right from the beginning, whereas almost all of the other grids keep the error relatively constant at first. A15 is the only other grid that shows similar deviation behavior as the cubic grid. Since E(¯ t)only takes the magnitude of the wave function into account, we also measure the root mean square error as RMSE(¯ t)= 1 m m i=0|¯ λ,¯ g,κ(¯ ri,n¯ t)−ψ(n) i|2, which considers the phase as well. The measurement outcomes are presented in Fig. 5b, which shows similar results as were observed earlier, as in the cubic grid being the least accurate and C15 and BCC the most. In this measurement, the FCC grid produces similar accuracy with the cubic grid at first, but closer to the end starts to deviate, by showing a slightly slower increase in the error. Interestingly, when the phase is taken into account, the error increases slower with the BCC grid than with C15, even though C15 maintained the vortex shape the longest. The Bogoliubov equations are also able to predict the splitting symmetries of the dynamically unstable stationary states, as in how many parts they split up into, when perturbed. With the stationary state used in our accuracy measurements, the prediction is that a 3-fold splitting should occur. Fig. 6shows visualizations of a cross-section of the wave function at ¯ z=0, after time integrated with the cubic and BCC grids, until E(¯ t) >0.01. The figure clearly shows how the splitting symmetry produced by the cubic grid doesn’t match with the prediction of the Bogoliubov equations but is 4-fold instead. With the BCC grid, a 3-fold splitting can be observed, which matches the prediction. Fig. 6. Visualizations of a cross-section of |¯ |2at ¯ z=0, after time integrations with the cubic and BCC grids, until E(¯ t) >0.01. 5.2. Computational performance At first, we test the computational performance of the GPU implementation by using the same stationary vortex state as in the accuracy measurements, but by increasing the dimensions of the computational domain, so that it becomes a cube, still satisfying the requirement eq. (15). We vary the problem size by scaling the spatial grid and monitoring its effect on the speedups the GPU provides compared to a serial CPU implementation. This experiment was performed with the BCC grid, and using an NVIDIA® GeForce®RTX 2080 Ti GPU, and an AMD Ryzen™9 3900 CPU. The length of an edge of the computation grid cube was varied from 14 to 270 when measured in replicable structure instance counts. Hence, the number of degrees of freedom varied between 143×12 and 2703×12, inclusive. 7
M. Kivioja, S. Mönkölä and T. Rossi Computer Physics Communications 278 (2022) 108427 50 100 150 200 250 10 20 30 40 50 60 Grid width Speedup Fig. 7. The computational performance improvements provided by the GPU against the CPU, as a function of problem size. The results are presented in Fig. 7, which shows a rapid increase in the speedups at first, when the problem size increases, but then reaches the maximum of 57-fold speedup, at the problem size of 1893×12 degrees of freedom, and stays approximately on that level from thereon. From that, we predict that this is the maximum reachable speedup with the used GPU and CPU. We look more deeply into the reasons behind this behavior by using the Roofline model [42]. For that, we pick four cases from the previous experiment, the grid widths of 14, 27, 54, and 189, and examine their computational performances in relation to their arithmetic intensities. In Fig. 8a the aforementioned four cases are numbered from 1to 4, respectively. The horizontal ceiling is the maximum peak performance of the double-precision pipeline of the used GPU, of 368 GFLOP s. The sloped ceiling is produced by the maximum memory bandwidth of 648 GB s. It can be observed that with the smallest problem sizes the GPU is not fully utilized and the performance falls significantly behind the possible maximum. When the problem size increases the performance approaches the ceiling, but at the same time, the arithmetic intensity decreases, making the performance memory bandwidth bound. There is an anomaly when it comes to the arithmetic intensity of the problem size number three, which might be caused by a lucky alignment of the borders of thread groups, minimizing the number of the required compute time global memory loads, mentioned in section 4.1. We repeat the performance measurements with the other most accurate grid type C15 and examine the results likewise using the Roofline model, shown in Fig. 8b. With the C15 grid, the distribution of the arithmetic intensities between different problem sizes is smaller than with BCC, but still indicating similar behavior of the largest problem size moving the performance from being computebound to bandwidth bound. Also as with BCC, the same anomaly in arithmetic intensities in the case of the problem size number three is observable. With the BCC grid, we were able to eliminate all of the shared memory bank conflicts, but couldn’t achieve that with C15, leaving some conflicts to the load transactions. This is likely one of the reasons why the C15 grid only achieves 89% of the performance of the BCC grid. Next we examine the computational performance on multiple GPUs. The measurements were performed using up to 8 NVIDIA® Tesla®P100 GPUs. As an initial value, we used a stationary vortex solution with λ =0.1, ¯ g=5000, and κ=20. For the discretization, the BCC grid was used, and a scale which produces 87 ×87 ×302 × 12 degrees of freedom. We examine the speedups provided by the additional GPUs compared to the performance on a single GPU, by first retaining the problem size constant and measuring T1 TNgpu , where TNgpu is the computation time on Ngpu GPUs. The results are shown in Fig. 9a, where we also present the function from Amdahl’s law [43] St(Ngpu)=1 (1−pt)+pt Ngpu ,(16) fitted to our observed data, yielding pt=0.92. St(Ngpu)is the speedup, and ptthe proportion of the computation time that the part benefiting from the additional GPUs originally occupied. With this, we can then calculate the maximum achievable speedup to be limNgpu→∞ St(Ngpu) =1 1−pt=12.5. The majority of the 8% of the computation time that doesn’t gain from the added GPUs, is likely caused by the latency of the memory copies between different GPUs. We then measure the speedups by fixing the computation time to T1on all of the considered GPU counts and adjusting the problem sizes accordingly. In this case, the speedups are defined as WNgpu W1, where WNgpu denotes the computation workload on Ngpu GPUs. In turn, to these measurements, we fit Gustafson’s law [44] Sw(Ngpu)=1−pw+Ngpu pw,(17) where pwis the proportion of the workload that the part benefiting from the additional GPUs originally occupied. The fitting produces pw=0.78 and both the measured and fitted speedups are presented in Fig. 9b. The most significant cause for pwbeing smaller than ptis that we scale the problem size equally in all spatial dimensions and if we denote this scale by sw, the amount of transferred data between the GPUs increases in O(s 2 3 w). Last, we compare the multi-GPU performance against a parallelized CPU implementation. The parallelization was done with domain decomposition by utilizing the MPI communication protocol, and the performance was measured using up to 48 12-core Intel®Xeon®E5-2690 v3 CPUs. The same GPUs, initial value, and a number of degrees of freedom were used as in the measurements presented in Fig. 9a. This time, we measure the computation speed in how many time units of ¯ tcan be integrated in one second, and present both the CPU and the GPU speeds in Fig. 10. It is observed that one GPU outperforms 144 CPU cores, and in fact, provides 152-fold speedup against one CPU core. The scalability of the CPU implementation exceeds the scalability provided by the GPUs, which we believe to be caused by the lower computational performance of the CPUs, making the computation occupy a larger part of the total time and workload, i.e. causing greater ptand pw in eq. (16) and eq. (17), respectively, than on GPUs. 6. Conclusions We presented a discrete exterior calculus-based numerical time integration method for the three-dimensional Gross-Pitaevskii equation. The spatial domain was discretized with six different natural crystal structures and the method implemented on multiple GPUs. We measured the accuracy provided by the different discretization grids and the computation time performance of the GPU implementation. All of the discretization grid types showed to provide better accuracy than the cubic grid. Especially the C15 and BCC grids were proved to be the most accurate ones in all of the experiments, which coincides with previous work showing them to be most numerically isotropic. The cubic grid also provided unreliable simulation results that didn’t agree with the predictions yielded by the physics theory, whereas with all of the other grids the results did match the predictions. The GPU implementation of the method was described with the emphasis on the performance optimizations. In which, we paid particular attention to the memory utilization, as we then showed 8