{'fraction_non_alphanumeric': 0.03392448070834728, 'fraction_numerical': 0.011299580448379707, 'mean_word_length': 4.6512, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Alternating direction implicit time integrations for finite difference acoustic wave propagation: Parallelization and convergence

This work studies the parallelization and empirical convergence of two finite difference acoustic wave propagation methods on 2-D rectangular grids, that use the same alternating direction implicit (ADI) time integration. This ADI integration is based on a second-order implicit Crank-Nicolson temporal discretization that is factored out by a Peaceman-Rachford decomposition of the time and space equation terms. In space, these methods highly diverge and apply different fourth-order accurate differentiation techniques. The first method uses compact finite differences (CFD) on nodal meshes that requires solving tridiagonal linear systems along each grid line, while the second one employs staggered-grid mimetic finite differences (MFD). For each method, we implement three parallel versions: (i) a multithreaded code in Octave, (ii) a C++ code that exploits OpenMP loop parallelization, and (iii) a CUDA kernel for a NVIDIA GTX 960 Maxwell card. In these implementations, the main source of parallelism is the simultaneous ADI updating of each wave field matrix, either column-wise or row-wise, according to the differentiation direction. In our numerical applications, the highest performances are displayed by the CFD and MFD CUDA codes that achieve speedups of 7.21x and 15.81x, respectively, relative to their C++ sequential counterparts with optimal compilation flags. Our test cases also allow to assess the numerical convergence and accuracy of both methods. In a problem with exact harmonic solution, both methods exhibit convergence rates close to 4 and the MDF accuracy is practically higher. Alternatively, both convergences decay to second order on smooth problems with severe gradients at boundaries, and the MDF rates degrade in highly-resolved grids leading to larger inaccuracies. This transition of empirical convergences agrees with the nominal truncation errors in space and time.


Introduction
Finite difference (FD) modeling methods have been actively applied to seismic wave propagation problems on the last fifty years. In this area, advanced methods use staggered grids (SG) in combination to fourth-and higherorder FD formulas to resolve the wave equation in first order formulation, and then solving for all dependent variables [2], [17], [19], [26], [35], [38] [41], [42]. In particular, the velocity-pressure formulation of the acoustic wave equation allows stating free surface and rigid wall boundary conditions as simple Dirichlet conditions, and is also well suited for the implementation of absorbing boundaries. In a SG, discrete wave fields may be defined on separate nodal meshes, which are displaced by the half of grid spacing h 2 in one or more coordinate directions. On these grids, an unknown wave field is placed at the center of those it depends, and thus SG differentiation yields more accuracy in comparison to a traditional nodal FD stencil with an error dependent on full h.
improve computational performance in available CPU and GPU scenarios. Section 4 states the numerical tests used to assess accuracy and convergence of these schemes, and to compare the speedup of execution times achieved by all their parallel implementations in Section 5. Our main conclusions are finally given in Section 6.

Numerical methods
In a medium Ω with density ρ and adiabatic compression κ, the propagation of acoustic waves obyes the following system of non-homogeneous differential equations Above, the dependent variables are the scalar pressure u and the particle velocity vector v = (v, w) , while f may correspond to an additional forcing term. Appropriate initial and boundary conditions must be considered in order to find a particular solution to (1), and this velocity-pressure formulation is well suited to state boundary conditions of relevance in seismic applications. For instance, a free surface can be modeled by a homogeneous Dirichlet condition on u, and most absorbing conditions require the computation of u along the domain boundaries. However, we present the following formulation of our numerical methods on the square domain Ω = [0, 1] × [0, 1], where a Dirichlet condition u(x, y, t) = u 0 (x, y, t) is imposed at (x, y) ∈ ∂Ω, and velocities at those edges must be computed during the simulation. Brief comments on the easy reformulation that also computes boundary values of u are appropriately given nonetheless. Initial distributions of all three dependent variables are given at t = 0.

The nodal CFD method
To simplify notation, we use a square N × N lattice to discretize Ω for a given integer N, and then cell sides have a common length h = 1 N . Grid values of the continuous fields u,v, and w are hold by matrices U, V, W, respectively. The result of FD differentiation of discrete pressures are collocated on matrices U x and U y , according to the direction of differentiation, and similar structures are used in the case of particle velocities. Using the CFD operators P and Q presented in Appendix A, pressure differentiation at all grid nodes is computed by the matrix operations U x P T = UQ T , PU y = QU.
However, pressure is known along boundary edges because of Dirichlet conditions, thus this formulation also makes use of the reduced matrixŪ that only holds interior grid values. Time updating of these inner discrete pressures is based on the discretization of momentum conservation in (1), and requires approximations to v x and w y . In space, this discretization in carried out by reduced operatorsP andQ, also introduced in Appendix A, applied to discrete velocitiesV xP T =VQ T ,PW y =QW.
Reduced matricesV andW are obtained from whole matrices after removing first and last rows of V, and first and last columns of W, respectively. By doing so, we avoid unnecessary velocity differentiation along boundaries. The application of the Alternate Direction Implicit (ADI) formulation to the acoustic model (1) starts by defining the discrete vector operators Next, the Crank-Nicolson discretization of time derivatives leads to the fully implicit method where ∆t is the time step, U CN represents the discrete solution [Ū, V, W] T at a given discrete time, and A F = [F, 0, 0] T collects the source contribution f at grid interior. Coupled equations of this sort are efficiently solved after the application of an ADI splitting, that only requires the solution of banded SEL for CFD differentiation along a single coordinate direction. The Peaceman-Rachford ADI decomposition of (4), as used in [10], states a calculation in two stages each of which alternates the application of operators A 1h and A 2h . The intermediate approximation U * CN is common on these ADI formulations, but only approximations at time levels t = m∆t and t = (m + 1)∆t are consistent solutions to (1). According to the definition of A 1h , the first stage allows an explicit calculation of vertical velocity W * in terms of its previous values at time t = m∆t, and the differentiation of pressure U m along the y direction. However, this stage also couples the calculation (row by row) of intermediateŪ * and V * in the following two sets of linear systems We denote the matrices on the right hand side of above equations by A and B, respectively, which can be computed This fixed point iteration stops when (||Ū * k+1 −Ū * k || < ε) and (||V * k+1 − V * k || < ε), for a given tolerance ε > 0. The second stage of the ADI decomposition (6) allows computing all wavefields at time t = (m + 1)∆t through a similar procedure. The vertical velocity V m+1 is computed explicitly using its intermediate value V * and the approximation U * x , according to the definition of A 2h . The columns of pressure U m+1 and velocity W m+1 are successively computed by the fixed point ADI iteration where matrices C =PŪ * − ∆t 2P κV * QT (P T ) −1 − F m+1 and D = PW * are obtained from intermediate variables, as one can see. In above formulation, matrix inversion is only used for notation convenience in auxiliary matrices A and C, but LU is actually employed to resolve the 2N linear systems involved in their calculations. As mentioned earlier, the additional 2N linear systems embedded into each fixed-point iteration in (8) and (9), are also solved by forward/backward substitutions, using the previous LU factorization of P andP. The above CDF formulation is amenable for including absorbing damping layers that require the computation of discrete wave fields along boundaries. By simply replacing reduced operatorsP andQ, and reduced matricesŪ,V, andW, by their full matrix versions in previous formulation steps as convenient, the boundary wave fields can be also obtained. These values can be further damped by the absorbing technique proposed in [7]. Even more, the versatile strategy developed in [44] for the damped wave equation can be also adopted after a straighforward reformulation of our CFD method.
A high-level implementation of above nodal CFD method is sketched in algorithm 1. This algorithm presents two sequential inner loops in lines 11 and 15, corresponding to the two ADI stages that finally calculate (Ū, V, W) m+1 from discrete solutions at t = m∆t. For clarity, each of these stages is broken down into a separate coding procedure given in Appendix C. At each of these procedures, there is no data dependency among the N linear systems solved, either row wise or column wise. Thus, each ADI set of N tridiagonal systems can be solved simultaneously by our parallel methods.
• K, R: Matrices of grid values of material properties κ and ρ −1 , respectively.
• c f l, c max : CFL stability parameter and the maximum value of wave speed c = κ/ρ.
• T sim : Physical time of the global simulation.
Moreover, we use some special index notations in this pseudo-code. MatrixX represents the reduced form of matrix X according to the formulations given in sections 2.1 and 2.2. Matrix indexing by : refers to every element along the corresponding dimension, while: only denotes the interior elements 2, 3, . . . , N − 1. Finally, . * represents the element wise vector/matrix multiplication, as used in Octave and MATLAB.

Complexity analysis of the nodal CFD method
As expected, the computational cost of this method is driven by the solution of the ADI embedded SEL. LU has been traditionally employed for solving tridiagonal linear systems due to its computational speed, in addition to a high numerical stability for diagonally dominant matrices (i.e., the Gaussian elimination can proceed with no pivoting). The algorithm only performs 2N steps to solve a given system, that results in an algorithmic complexity of O(N). In this work, we also omit any pivoting strategy, even though CFD matrices P andP, are not strictly diagonally dominant. Numerical results do not reflect any instability signs arising from this restraint.
The operation complexity depends on the grid density N, that sets the dimension of vectors and matrices used by the nodal CFD method. In case of reduced matricesP,Q,Ū, and others, N − 1 ≈ N, for N >> 1. The total number of time iterations simply results I ∆t = N c max T sim c f l , where T sim is the physical simulation time. In this analysis, let us consider the matrix operations in Table 1.  Table 2 lists the amount of matrix operations performed by the nodal CFD method during I ∆t time iterations. From there, one can easily see that the time complexity of this scheme is of polynomial order, and highly dominated by the row wise and column wise ADI cycles (see statements 11 and 15, in the algorithm 1). In particular, the complexity dominant term is where k max is the maximum number of iterations spent by each ADI stage to converge. In all numerical experiments carried out on this work, k max varies from 6 to 8, and thus we expect k max being of O(1) in similar wave propagation applications.

Section
Maximum Additional comments on the CFD implementation Along our optimization and parallelization process, we implement and test slightly different versions of the ADI embedded iterations (8) and (9), in search for a higher level of parallelism. In particular, matrix V * k+1 might be computed from pressure U * k in the second step of (8), so this calculation can proceeds simultaneously with the first step that updates U * k+1 . A similar strategy can be applied in (9) for a simultaneous computation of U m+1 k+1 and W m+1 k+1 , making this ADI formulation more suitable for parallel platforms. This simple idea is similar to the one behind the iterative Jacobi linear solver, fully parallelizable in comparison to its Gauss Seidel counterpart. However, this Jacobi-like ADI actually requires few more iterations for convergence, than the above implementation that delays V * k+1 computation until U * k+1 is available. In addition, this ADI choice also leads to a reduction of the CFL stability limit of this nodal method around 90%, and this deficiency is even higher in the case of the mimetic SG method. Thus, we discard this implementation choice in current developments.

The staggered MFD method
In this formulation, the grid distribution of discrete wavefields takes place on a rectangular staggered grid. In analogy with the nodal mesh used by the previous CFD method, let us consider the 1-D grid of N cells of length h along the x direction, with nodes x i = i * h, and cell centers x i+1/2 = x i +x i+1 2 for i = 0, · · · , N − 1. These grid nodes are collected in vector X n , while X cb is another vector comprising cell centers in addition to grid edges, that is X cb = (0, x 1 2 , ..., x N− 1 2 , 1). After defining similar grid vectors along the y direction, the staggered distribution of discrete pressures correspond to the mesh locations given by X cb ⊗ Y cb , where ⊗ represents the standard cartesian product operator used in set theory. Alternatively, grid sites of discrete horizontal and vertical velocities are given by X n ⊗ Y cb and X cb ⊗ Y n , respectively. According to these staggered distributions, matrix dimensions of discrete wavefields are U ∈ R (N+2)×(N+2) , V ∈ R (N+2)×(N+1) and W ∈ R (N+1)×N , in contrast to the previous nodal formulation where all these matrices are (N + 1) × (N + 1) square. Due to the boundary conditions, we again make use of reduced matricesV andW, obtained from original V and W after removing rows from the former and columns from the latter, as before.
The staggered FD differentiation employs the mimetic fourth-order operators G 4 and D 4 given in Appendix B, for an approximation of partial derivatives along both directions of all three wavefields Next, we replace the nodal CFD definition of discrete vector operators A 1h and A 2h by the following staggered approach With above definitions at hand, the formulation of this staggered MFD method proceeds by replicating the time discretization and Peaceman-Rachford two-stage ADI decomposition, applied in the previous section. A pseudo-code of this method is given by the algorithm 2. For completeness, the two ADI stages are broken down into separate procedures, and given by the algorithm 4 in Appendix C. The algorithm 2 employs the same parameters used by the nodal algorithm 1, as well as similar matrix and indexing notation.

Complexity analysis of the staggered MFD method
The complexity analysis of this method considers the same matrix operations defined in Table 1. The new Table  3 lists the total amount of these operations performed by the staggered method during I ∆t time iterations. Similar to the nodal scheme, the operational cost is a N-dependent polynomial, with a leading term given by the solution of the ADI stages (referred in statements 11 and 15 of the algorithm 2), that we denote as T MFD . Specifically, T MFD ≈ 28 * I ∆t * k max * N 2 = O(N 3 ). Thus, the total amount of operations is less than those employed by the nodal method. However, the tridiagonal structure of embedded linear systems on the former method, and our LU based solution strategy, make equivalent the operational complexity of both methods.

Computational optimization
The implementation of the methods presented in Section 2 takes place on an Intel(R) Core(TM) i7-4770 CPU with 4 cores running at 3.4 GHz clock speed, and 8 MB of cache memory. This architecture counts with a NVIDIA GTX 960 Maxwell GPU card and the Ubuntu Linux 16.04 LTS operating system.

Sequential versions
Original codes of the proposed methods were developed in Octave and used the LU decomposition of CFD matrices P andP to solve ADI embedded SEL [10]. To better use the capabilities of our CPU architecture and apply

Algorithm 2 The staggered MFD method
∆t ← c f l Ncmax , I ∆t = T sim ∆t 3: Ū m+1 ,W m+1 ←ADI-columns(U * , W * , C, D) 16: t ← t + ∆t 17: end for 18: return U, V, W 19: end procedure some code optimization techniques, we migrate these original codes to C++, and use the 4.8.4 − 2 version of the gcc compiler in this process. The first optimization step is based on using different compilation flags, among which the −O3 flag significantly improved our C++ codes, leading to an immediate performance improvement. This −O3 flag automatically introduces the following software optimizations: loop unrolling, function in-lining and automatic vectorization. A second set of optimizing strategies are related to matrix orientation, and its application highly affects to the new CFD code. The original code in [10] solved the inherent SEL's in a vector oriented fashion, i.e. matrix updates and system solutions are performed a row/column at a time. To take better performance advantage of the algebra libraries, a matrix oriented approach is taken for the C++ implementation. That is, updates and system solutions proceed simultaneously on the full matrix. It can be seen that most of vector operations in both algorithms CFD and MDF decouple, so they can be performed as a single matrix operation. In particular, the CFD SEL's with multiple right hand sides and same matrix, can be solved simultaneously.

Parallel versions
For multiprocessing, we use the Octave parallel package to distribute computation along rows/columns among several processes. At first, one process is used for each row/column, but this solution performs extremely bad. Little work is assigned to each process, resulting in most of computing power being wasted in managing these processes and their communication. A second attempt consists of using a small number of processes (up to 8 in an eight core CPU), and then lists of columns/rows are assigned to each of these processes. Performance results are slightly better, however the single-threaded implementation is still faster for all tested numerical meshes, with the grid of biggest size (N = 1024, as detailed in Section 4), as the only exception. In general, results of using this Octave parallel package are poor, because of the need of data replication as required by the multiple launched processes. Up to date, no mechanism is available for fast memory sharing/synchronization among processes in this parallel package.
Due to results discussed above, we return to the original vector orientation of this problem, as it parallelizes naturally thanks to the decoupled operations at the vector level. The partition proceeds at the row/column level, distributing each vector in a thread pool that performs the required operations. Our C++ implementations follow and exploit this data orientation. The next optimization step applied to these vector-oriented versions is loop parallelization using OpenMP. OpenMP allows us to take full advantage of Intel multi-core CPU architectures. In our OpenMP loop parallelization, we have taken into account two intrinsic aspects of our numerical discretization: (i) data dependencies, and (ii) stencil dependencies. Each parallel loop collects a common set of updating instructions applied on a particular grid zone (boundary, near boundary, or interior nodes). Scheduling of OpenMP threads is established by the flag auto, which invokes as many threads as the number of available CPU cores. We observe in our experiments that this option provides the best performance for a wide range of scenarios. Increasing the number of threads over the number of cores usually entails an overhead due to communication conflicts, whereas using less threads than available cores might lead to an underuse of a multi-core platform.

GPU implementation
Once CPU parallelization has been appropriately exploited, we address a new GPU implementation of our numerical schemes. On the available architecture, the programming language is CUDA C based on the nvcc compiler. In this migration process of our original sequential codes, we also use the Toolkit 7.5.17 and functions of the cuBLAS library, and only develop specific kernels for particular functionalities not supported by cuBLAS.
We next comment on crucial implementation details. Due to the relatively small sizes of the matrices involved in our test cases (up to 1024 × 1024), the vector oriented workload does not yield enough parallelism to fill the GPU. Thus, we return to the matrix approach in our GPU implementations, and instead of partitioning at the row/column level, we focus on optimizing the matrix operations themselves at the element level. This allows achieving a much higher level of parallelism, given that parallelization is performed throughout the full extent of the algorithm, and not only at the row/column solvers. 1 .
As already mentioned, we implement some missing banded matrix functionalities in cuBLAS. For instance, the Level 3 (L3) product of a banded matrix times a dense matrix is not provided, while the Level 2 (L2) multiplication of a banded matrix times a vector is indeed available. We build the L3 operation over the L2 one by decomposing a matrixmatrix multiplication into independent matrix-vector products and accumulation operations. Our first approach, using streams to overlap the matrix-vector products, achieves a very low occupancy due to the cost of kernel launching, relative to the workload of each thread. We found that best results are obtained when the whole matrix is computed in a single launch, as the work per thread increases. Here, a higher number of threads are effectively launched, increasing the opportunities to obtain better occupancy. Same results are achieved by using stream-based distances and matrix norm computations.
On the other hand, our GPU implementations need to apply reduction operations, like for instance, those involved in matrix norm computations. In this case, parallelization of the reduction process is performed in a hierarchical fashion. Each thread adds a pair of elements, and then half of the threads add the partial results computed. This process is repeated for log 2 (N) iterations, until the whole list is reduced to a single value. Several practical considerations must be considered to achieve optimal performance: minimal divergence, array indexing to optimize data locality, unrolling of some operations, etc. Bigger reductions are computed in two steps. First, one reduces in parallel all the columns of the matrix to a vector in shared memory. Then, this vector is reduced again to a single value. Here, we have used this reduction technique to calculate the matrix Frobenius norm in our new kernels, as required to test ADI convergence according to algorithms in Appendix C.
Finally, we want to notice that the stopping criterion of ADI iterations in algorithms 3 and 4 of Appendix C, imposes an important constraint for GPU performance. The number of steps required until convergence is actually limited by the constant parameter k max , that varies from 6 to 8, in our experimental set. However, convergence on the Frobenius norm of matrices U, V, and W, updated by GPU processing, require that a floating point value must be transferred form device to host, for testing against ε tolerance. This is wasteful as the PCIe bus is slow and has a high latency. To circumvent this problem, we study the average number of steps until convergence. In our cases, this number is roughly 6, so we thus perform this minimum number of iterations until convergence start being tested, and saving most of the unnecessary device to host transfers.
Concerning our new CUDA schemes, we want to emphasize that checking for inner ADI convergence forces the synchronization of kernel threads, slowing down the global computation, and leading to a bottleneck that hampers performance. As it was said, convergence is almost always achieved in a constant number of iterations for each test case, so we opt for executing a minimum number of iterations before any checking take place.

Numerical tests
We next present the family of numerical tests that allow assessing accuracy and convergence properties of both the CFD and MDF schemes. Same tests are used to evaluate their computation speedup gained through the optimization stages described previously, but this discussion is given in next section. Under appropriate initial and Dirichlet boundary conditions (BCs), the propagation model (1) has the following exact solution This solution u is time harmonic with period T and presents a spatial behaviour controlled by the polynomial coefficient Γ, being k a possitive integer, and the spatial period λ of the sinusoidal component. These four parameters allow adjusting the difficulty of the particular problem numerically solved. In particular, the test solution can switch from a purely sinusoidal function for Γ = 0, to one with significant gradients at the boundary neighbourhood in the case of Γ and k with large magnitudes. We first solve the case of Γ = 0 where BCs reduce to homogeneous, consider a unitary wave speed c = 1 (i.e., κ = ρ), and take λ = 1 4 and T = 1 √ 2 . We perform ten numerical simulations of each scheme for an increasing number of grid cells N = 16, 24, 32, 48, 64, 96, 128, 256, 512, 1024, so grid resolution improves with the reduction of h = N −1 . In each case, the global simulation time corresponds five time periods, namely T S I M = 5T , and the time step ∆t is taken as the maximum permitting stable calculations, where ∆t = hc −1 c f l max . The term c f l max is the well known CFL bound that varies between our two FD ADI schemes, and were found in [10]. Especifically, we employ the c f l max values of 0.91 and 0.81 in our CFD and MDF simulations, respectively. Next, a second set of harder tests with non homogeneous BCs are also numerically solved, where Γ 0. In this case, we only inspect the representative cases of k = 2 and k = 9, and also assume Γ = k to simplify this current parametrization. Figure 1 shows absolute errors in Frobenius norm found on numerical tests performed under homogeneous and heterogeneous BCs. In general, these errors steadily decay with h validating our new parallel implementations. On both tests with mild difficulty, i.e., Γ = 0 and Γ = k = 2, the MFD method is slightly more accurate on most of the explored grids with N ≤ 512, because staggering halves the grid spacing during numerical differentiation, as compared to the nodal CFD scheme. However, when solving the hardest test with Γ = k = 9, MDF errors on simulation with N = 1024 suffer a clear stagnation, making useless any further grid refining. Thus, we limit this experimental study, and subsequent speedup analysis, to meshes where N ≤ 1024.
In addition, we estimate the convergence rates of the CFD and MDF solutions as h reduces, and list them in tables 4 and 5, respectively. These estimations are based on the standard error weighting from two consecutive simulations. For each test, we also average the resulting rates after omitting the highest and lowest values, to quantify  the overall convergence behaviour avoiding main outliers. For both schemes, convergence degrades as the problem becomes harder, being almost quartic on the homogeneous test (limited by the spatial discretization), and reducing to quadratic on the hardest heterogeneous test (more constrained by the order of CN time integration). We would like to remark that accuracy and convergence results reported in [10], clearly show a MDF superiority on similar tests under homogeneous BCs. This disadvantage of the former CFD method arises from a reduction of nominal accuracy from fourth to third order at grid boundaries. Here, we opt for a fully fourth order set of CFD stencils at all grid points when implementing our current method.

Figures 2a and 2b
show speedups in log2 scale, relatives to the original sequential Octave codes, of each implementation delivered by the optimization process discussed in Section 3. These results were obtained after solving the numerical test under heterogeneous BCs for k = 9. Speedup values observed on the alternative heterogeneous or homogeneous tests, are very similar to the ones shown in these figures, confirming that evaluation of the forcing function f , and related operations at boundaries, minorly affect computation times.
The new CFD and MDF implementations using parallel Octave report a poor performance in most of our tests, being the former better that finally improves to a 1.19x ratio at N=1024. Conversely, sequential C++ implementations yield a global improvement that reaches 5.39x (CFD) and 2.6x (MFD) at the finest grid. This improvement is mainly (a) (b) Figure 2: Speedup values in log2 scale achieved by the sequential C++ and parallel codes, relative to the sequential Octave implementations, of both methods the nodal CFD (a) and the staggered MFD (b), when solving the heterogeneous test for k = 9. due to software optimizations performed by the the compiler. Now, C++ parallel codes that use OpenMP achieve a higher improvement at all grids, and limiting speedups reach 11.1x (CFD) and 6.74x (MFD), at N = 1024. In this case, the main reason that favors speedup is using optimized LAPACK solvers, combined to the parallelization of some operations, and to the reduction of L1-L2 cache access failures. The OpenMP implementation was almost trivial to develop from the vector oriented approach, as vectors decouple and no synchronization is needed among threads. Finally, we observe the best performance results on the GPU implementations, where the CFD GPU speedup grows with N until 38.85x at the best resolved grid. This is a clear indication of the high level of parallelism achieved by our GPU implementation, on the range of grid sizes tackle in this work. On the other hand, speedup of the CUDA MFD method seems to converge steadily on grids with N ≥ 96 to a limit of 40x, with a minor improvement at the limiting mesh, and this behavior reveals a soon high occupancy of GPUs. As a complementary set of speedup results, we list in Table 6 only the ratios achieved by the parallel CFD and MFD methods, developed either in C++ using 8 threads or in CUDA, relative to their sequential C++ implementations. In most of our grids, these speedups improve with the mesh size. In particular, the improvements accomplished by the GPU implementations grow steadily, given the high device occupancy by our bigger matrices. It is worth reminding that grid size is limited in our tests by the numerical precision achieved by our methods, and higher refinements would not lead to more accurate results, at least for the MDF scheme. Thus, we do not test our parallel implementations with bigger matrices. Again, using the sequential C++ computations as reference, the OpenMP implementations only reach limiting improvements of 2.06x (CFD) and 2.59x (MFD) at N = 1024, due to matrix sizes are not big enough to balance expenses on either, the parallel LAPACK LU solutions of the nodal scheme, or the several matrixvector products of the staggered method. On the contrary, performances gained by the GPU implementations achieve maximum speedups of 7.21x (CFD) and 15.80x (MFD). Keeping matrices in GPU memory and minimizing hostdevice transfers were essential to obtain such performances. In these cases, the usage of streams was not effective giving that launching times were too high with respect to the kernel working times.
The testing process of ADI convergence (inner loops of algorithms 3 and 4 in Appendix C) may result a critical restricting factor to performance, as it demands several small device-host DMA transfers. As the number of loop iterations until convergence is roughly constant, we minimize this bottleneck by delaying checking the stopping criterion until a fixed number of iterations is performed.

Conclusions
This paper first focuses on the computer parallelization of two fourth-order FD methods to model 2-D acoustic waves on available multi-core (CPU) and many-core (GPU) architectures. The first scheme uses CFD operators on nodal grids and requires solving a tridiagonal system of linear equations (SEL) to compute coupled derivative values along each coordinate grid line. Alternatively, the second scheme employs MFD stencils on SG for local differentiation free of SEL solution. Both methods share a Crank-Nicolson (CN) time integration combined to a Peaceman-Rachford splitting that lead to an efficient ADI updating scheme. These FD ADI methods have shown better stability and convergence properties than an explicit Leapfrog integration in a previous work [10], but the computing expenses of ADI-embedded iterations demand parallelization to be competitive on practical propagation problems. Along the implementation process of each method, we present a sequential C++ code with optimal compilation flags, and develop three different parallel versions. Precisely, an Octave and a C+ multithreading code, both fully exploiting array-oriented storing for matrix and vector operations, and a CUDA implementation on a NVIDIA GTX 960 Maxwell. Performance speedups of parallel implementations were assessed relative to the computing times of the sequential C++ code, but we additionally report results with respect to the original Octave version. These results can be insightful for someone replicating our optimization procedures for similar numerical methods. Results of GPU methods show a significant speedup gain over their parallel counterparts. Main reasons are a high occupancy of GPUs on explored grid densities, better realized by the CFD method, and a minimization of CPU-GPU transfers by delaying the convergence testing of ADI-embedded iterations, until an expected number of cycles is actually performed. Relative to Octave, the GPU CFD method achieves a speedup of 38.85x at the finest grid, while the GPU MDF reaches a higher 41x. Same performances, but measured with respect to optimal sequential C++ versions, are 7.21x for the former and 15.81x for the latter. Keeping matrices in GPU memory and minimizing host-to-device transfers result essential in obtaining such performances.
Secondly, numerical experiments carried out in this study also enlighten convergence properties of these methods on a test suite of increasing difficulty. On a problem with an exact harmonic solution, the MDF accuracy was slightly better, and both methods exhibit convergence rates close to 4, as the differentiation order used on the spatial discretization. However, both convergences decay to second order on smooth problems with significant gradients at boundaries, that may reveal a stronger dependence on the time integration errors. In any case, these second order experimental rates are consistent to the theoretical global convergence of 2, and limited by the time discretization.

Appendices Appendix A: Compact finite difference matrices
Lele in [25] presents a Taylor-based construction of compact finite difference (CFD) formulas of second, fourth and higher order accuracy on either nodal or staggered grids. In this work, we use the fourth order accurate CFD matrix operators P and Q that allow differentiation of a smooth function v on [0, 1], whose evaluations at nodes x i = i * h, for h = 1 N and i = 0, · · · , N, are the components of vector V. Nodal approximations to dv dx result from solving the banded linear system PV x = QV, for the following (N + 1) × (N + 1) matrices P and Q . . . . . .
Above, matrix sub-indexes reflect the accuracy choice among the alternative CFD operators given in [25]. This formulation is quite general yielding a three-parametric family of CFD operators with fourth order accuracy, at least, and those given above result from reducing matrix bandwidth to assure P tridiagonally. An alternative pair of CFD operators can be derived from each P and Q set, by omitting the approximations to dv dx at grid edges. In this case, the reduced vectorV x ∈ N−1 only collects values of v at interior grid nodes, and the new CFD approximation to dv dx at point x 1 result from the linear combination of original formulas at nodes x 0 and x 1 . The same formula is now used at x N−1 , after switching signs of CFD coefficients to account for backward differentiation. This simple process leads to the reduced matrix operatorsP andQ of dimensions (N − 2) × (N − 2) and (N − 2) × N, respectively, that take the following form in the case of 14 given abovē 0 · · · 0 1 4 1 0 · · · 0 0 1 4 1 0 · · · 0 . . . . . . . . .
Note thatP 4 is also tridiagonal which permits using same algorithm on the calculation of both V andV x .

Appendix B: Staggered finite differences
For fourth order finite differentiation on staggered grids (SG), we employ the 1-D mimetic operators developed in [5] and [6], that have been applied to wave propagation by [10], [36], [37], [38] and [40]. On the interval [0, 1], we retake the SG used in section 2.2 with N cells of length h = 1 N , where nodes x i = i * h and cell centers x i+1/2 = x i +x i+1 2 correspond to the evaluation points of two given scalar functions v and u. Specifically, nodal v evaluations are collected in vector V ∈ N+1 , while evaluations of u at cell centers and grid edges conform vector U ∈ N+2 . Staggered differentiation is carried out by two separate FD matrix operators D and G, where components of DV approximate dv dx at every cell center, while GU yields approximations to du dx at all nodes. For clarity, the following figure illustrates the spatial distributions of components of V, U, DV and GU on a 1-D SG By construction, each D or G operator with fourth order accuracy, comprises an individual set of SG FD stencils dependent on three free parameters α, β and γ. An optimal parametrization has been recently proposed in [9] to preserve the self adjointness between mimetic D or G, an important property satisfied by their continuous differential counterparts divergence and gradient operators. However, we here adopt same parameter values α = 0, β = 0, and γ − 1 24 used in [10], [37], and [38] that lead to D and G with minimun bandwidth, to favor computational efficiency.