Hardware Acceleration of Tsunami Wave Propagation Modeling in the Southern Part of Japan

: In order to speed up the calculation of tsunami wave propagation, the ﬁeld-programmable gate array (FPGA) microchip is used. This makes it possible to achieve valuable performance gain with a modern regular personal computer. The two half-step MacCormack scheme was used herein for numerical approximation of the shallow water system. We studied the distribution of tsunami wave maximal heights along the coast of the southern part of Japan. In particular, the dependence of wave maximal heights on the particular tsunami source location was investigated. Synthetic 100 × 200 km sources have realistic parameters corresponding to this region. As observed numerically, only selected parts of the entire coast line are subject to dangerous tsunami wave amplitudes. The particular locations of such areas strongly depend on the location of the tsunami source. However, the extreme tsunami heights in some of those areas can be attributed to local bathymetry. The proposed hardware acceleration to compute tsunami wave propagation can be used for rapid (say, in a few minutes) tsunami wave danger evaluation for a particular village or industrial unit on the coast.


Introduction
To address the problem of timely reporting of the danger of a tsunami wave at a particular coastal zone location, we apply the facilities of modern computer architectures to accelerate the numerical simulation of tsunami wave propagation. The main goal of this paper is to demonstrate the advantages of a new super rapid method for modeling tsunami wave propagation. We do not wish to obtain the detailed wave height maxima distribution along the specified coast. According to [1], rather strong earthquakes (and, therefore, dangerous tsunami waves) happen with an approximately 150-year period in the part of the Nankai subduction zone near Shikoku Island and Kii Peninsula. As the last high-impact event in the region was in 1854 (166 years ago), the next one can be expected in the near future. We are therefore running numerical experiments via digital bathymetry which approximates the real bottom topography of a corresponding water area.
Our aim is to use just a regular modern personal computer (PC) and to achieve valuable performance gain exploiting graphic processing units (GPUs) and field-programmable gate arrays (FPGAs) as co-processors. We propose a specific hardware configuration and the corresponding modeling code.
The following general scheme for tsunami danger prediction is assumed. The measurements of the propagating tsunami wave (time series) are obtained in real-time mode from DART (Deep-ocean Assessment and Reporting of Tsunami) buoys (https://nctr.pmel.noaa.gov/Dart/) or any other real-time

Shallow Water System
Following most tsunami phenomena researchers, we simulate wave propagation with the help of a version of a shallow water system in the following form: v t + uv x + vv y + gH y = gD y where H(x,y,t) = η(x,y,t) + D(x,y) is the distance from the sea surface to the bottom; η is the sea surface displacement relative to the mean sea level (wave height); D(x,y) is the water depth (which is assumed to be known at all grid points); u and v are the water flow velocity components in the x and y directions, respectively; and g is gravity force acceleration, see [6,7]. Note that lower indices t, x, y in (1) indicate partial derivatives with respect to time and space variables. The shallow water system (1) does not take into account dispersive effects because in the case of local propagation of sufficiently long tsunami waves, these effects have no time to valuably change the wave parameters at the coast. The use of the shallow water approximation is typical among tsunami researchers.
An extended reference list exists concerning the numerical treatment of nonlinear PDE systems of hyperbolic type, to which system (1) belongs. In the aforementioned software MOST, the splitting method with respect to spatial variables is widely used for calculation.

MacCormack Scheme
The two half-step MacCormack scheme (see [11]) was used for numerical approximation of the shallow water system (1). This algorithm is an explicit difference scheme which uses the three-point stencil of a "cross" type. The calculation at each time iteration is arranged in two steps. Each turn, the values of wave height and components of the velocity vector are computed using the necessary data at the same point and two adjacent points of the grid at the previous time value; see [12][13][14] for the calculation stencil description. Below is the version of the MacCormack scheme we used for numerical experiments.
Appl. Sci. 2020, 10, 4159 Here, H n ij , u n ij , and v n ij are the gridded variables which correspond to the H, u, and v functions in differential system (1). The parameters τ, ∆x, and ∆y are the time step and spatial steps of the computational grid. In order to account for the spherical shape of the Earth, we use a decreasing grid step with respect to longitude for larger values of latitude. The notation F n ij represents variables at time layer n, F n+1 ij represents intermediate values, and F n+1 ij corresponds to the variables at time layer n + 1. Our implementation of the MacCormack scheme was previously numerically tested against the available analytical solutions of shallow water ray approximation for the cases of model-type bottom topographies; see [15]. In the littoral areas that are considered in the referenced paper, the size of the computation domain is 1000 × 1000 km in both cases. In the first case, the depth linearly increases according to the formula D(x,y) = 0.01y. In the second case of parabolic bottom topography, the depth increases according to the formula D(x,y) = 10 −8 ·y 2 , where y stands for the offshore distance. In Figure 1 we present a comparison of the shallow water system solutions obtained by our FPGA-based "Calculator" and by MOST software with the exact solution in the case of sloping bottom topography (where depth linearly increases along the ordinate axis). The wave is initiated by rounded water surface elevation with radius equal to 50 km and height limited to 2 m. Its center is situated 300 km off the shore. As can be observed in Figure 1 (cf. [12]), the similarity is rather good, while the proposed scheme is better in terms of wave front simulation when compared to MOST software.
Here, H n ij, u n ij, and v n ij are the gridded variables which correspond to the H, u, and v functions in differential system (1). The parameters τ, Δx, and y are the time step and spatial steps of the computational grid. In order to account for the spherical shape of the Earth, we use a decreasing grid step with respect to longitude for larger values of latitude. Our implementation of the MacCormack scheme was previously numerically tested against the available analytical solutions of shallow water ray approximation for the cases of model-type bottom topographies; see [15]. In the littoral areas that are considered in the referenced paper, the size of the computation domain is 1000 × 1000 km in both cases. In the first case, the depth linearly increases according to the formula D(x,y) = 0.01y. In the second case of parabolic bottom topography, the depth increases according to the formula D(x,y) = 10 −8 ·y 2 , where y stands for the offshore distance. In Figure  1 we present a comparison of the shallow water system solutions obtained by our FPGA-based "Calculator" and by MOST software with the exact solution in the case of sloping bottom topography (where depth linearly increases along the ordinate axis). The wave is initiated by rounded water surface elevation with radius equal to 50 km and height limited to 2 m. Its center is situated 300 km off the shore. As can be observed in Figure 1 (cf. [12]), the similarity is rather good, while the proposed scheme is better in terms of wave front simulation when compared to MOST software.  [15] (brown lines) and the numerical solutions using the field-programmable gate array (FPGA)-based "Calculator" (red dashed lines) and using the Method of Splitting Tsunamis (MOST) package (blue lines). The offshore distance is measured along the vertical axis relative to the figure's bottom boundary.

Hardware Code Acceleration
According to the MacCormack finite difference approximation, one needs to obtain the values of the three functions H, u, and v at each mesh point. To leverage the advantages of the FPGA microchip features, the stream processor architecture was proposed as an implementation of the numerical algorithm. The proposed Calculator contains several processor elements (PEs). Each of these elements forms a pipeline with a sequential data stream. The "onboard" memory contains all  [15] (brown lines) and the numerical solutions using the field-programmable gate array (FPGA)-based "Calculator" (red dashed lines) and using the Method of Splitting Tsunamis (MOST) package (blue lines). The offshore distance is measured along the vertical axis relative to the figure's bottom boundary.

Hardware Code Acceleration
According to the MacCormack finite difference approximation, one needs to obtain the values of the three functions H, u, and v at each mesh point. To leverage the advantages of the FPGA microchip features, the stream processor architecture was proposed as an implementation of the numerical algorithm. The proposed Calculator contains several processor elements (PEs). Each of these elements forms a pipeline with a sequential data stream. The "onboard" memory contains all the necessary Appl. Sci. 2020, 10, 4159 5 of 14 information. The acceleration of calculations by the FPGA architecture is due to the inner Block Random Access Memory (BRAM) access for implementing the stencil buffer.
The Calculator architecture makes it possible to process several nodes in parallel. At the same time, the user can connect a number of PEs to make several iterations. So, the computation pipeline could be optimized depending on the features of the FPGA microchip in use. MacCormack finite difference approximation fits very well with the Calculator architecture processing one node in one computer clock cycle.
The originally designed "Calculator" hardware, based on a Xilinx Virtex-7 VC709 FPGA microchip, was used for numerical tests; see [13,14] for details. A general view of the designed printed board is given in Figure 2.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 15 the necessary information. The acceleration of calculations by the FPGA architecture is due to the inner Block Random Access Memory (BRAM) access for implementing the stencil buffer. The Calculator architecture makes it possible to process several nodes in parallel. At the same time, the user can connect a number of PEs to make several iterations. So, the computation pipeline could be optimized depending on the features of the FPGA microchip in use. MacCormack finite difference approximation fits very well with the Calculator architecture processing one node in one computer clock cycle.
The originally designed "Calculator" hardware, based on a Xilinx Virtex-7 VC709 FPGA microchip, was used for numerical tests; see [13,14] for details. A general view of the designed printed board is given in Figure 2.  Table 1 summarizes the achieved performance results of the proposed Calculator (using MOST software) in comparison with those of sequential and several alternative parallel architectures.

Numerical Results
A number of numerical experiments using the MacCormack scheme were arranged to study the detailed distribution of tsunami wave maximal heights along the coast for several particular locations of the tsunami source. We moved the model tsunami source (sea surface displacement at the source) in the alongshore (four different positions i = 0, 1, 2, 3) and cross-shore (three different locations indicated by j = a, b, c) directions to compare the distributions of tsunami wave maxima along the coast line. The positions of all these 12 model tsunami sources are given in Figure 3. The position of the Si-b (i = 0, ..., 3) source's negative wing (part of the initial water surface displacement) is the same as that of the Si-a positive wing, and the position of the Si-b positive wing is the same as that of the Si-c source's negative wing (i = 0, ..., 3). Next, we describe how these experiments were conducted and discuss the obtained results.

Digital Bathymetry
In Figure 3 we show the digital bathymetry of the water area under consideration (the region of Shikoku Island and Kii Peninsula, south Japan). It was developed in [13] on the basis of data produced by the Japan Oceanographic Data Center (JODC) [16].  Table 1 summarizes the achieved performance results of the proposed Calculator (using MOST software) in comparison with those of sequential and several alternative parallel architectures.

Numerical Results
A number of numerical experiments using the MacCormack scheme were arranged to study the detailed distribution of tsunami wave maximal heights along the coast for several particular locations of the tsunami source. We moved the model tsunami source (sea surface displacement at the source) in the alongshore (four different positions i = 0, 1, 2, 3) and cross-shore (three different locations indicated by j = a, b, c) directions to compare the distributions of tsunami wave maxima along the coast line. The positions of all these 12 model tsunami sources are given in Figure 3. The position of the Si-b (i = 0, ..., 3) source's negative wing (part of the initial water surface displacement) is the same as that of the Si-a positive wing, and the position of the Si-b positive wing is the same as that of the Si-c source's negative wing (i = 0, ..., 3). Next, we describe how these experiments were conducted and discuss the obtained results.
According to the stability condition, the time step used in the computations was 0.5 sec. Mostly, regular rectangular computational grids are used for tsunami numerical modelling. The grid could also be matched with one of a geography projection. We took into consideration the spherical shape of the Earth by decreasing the grid step with respect to longitude for larger values of latitude. Matching the grid to geographic coordinates, we suppose that the length of one arc degree at latitude λ is equal to the length of one equator arc degree multiplied by cos(λ). Using the advantages of the proposed Calculator, it takes only 25 s to compute tsunami wave propagation over the water area in Figure 3. In practice, the travelling time of a real wave through this water area is 3200 s.
During the numerical experiments, we moved an artificial tsunami source along the Tokai-Nankai deep water trench. Note that, according to [1], a strong seismic event is expected in this area. The authors discussed the disastrous tsunami that happened in this region nearly 300 years ago and provided a reconstruction of the tsunami source based on historical records (see Figure 4).

Digital Bathymetry
In Figure 3 we show the digital bathymetry of the water area under consideration (the region of Shikoku Island and Kii Peninsula, south Japan). It was developed in [13] on the basis of data produced by the Japan Oceanographic Data Center (JODC) [16].
We conducted 12 numerical experiments for model tsunami sources (indicated as Si-j in Figure 3). The sources are of the same shape with elevation and subsidence parts (details are given in the next subsection). The selected shape of the initial sea surface displacement has four different positions along the coast (indicated with subscript i = 0, 1, 2, 3) and three values of offshore distance (with the indication j = a, b, c).
In Figure 3, only the sources closest to the coast, Si-a (i = 0, ..., 3), and those more distant from the coast, Si-c (i = 0, ..., 3), are presented. Their positive wings (parts) are outlined in pink, and negative wings (water surface depression) are outlined in yellow. The respective middle-posed sources Si-b (i = 0, . . . , 3) are located between sources Si-a and Si-c.
The digital bathymetry and, correspondingly, the grid for computation are as follows.
The size of the computational domain is 3000 × 2496 points (approximately 842 × 556 km). The grid steps in the west-east and south-north directions are 0.003 and 0.002 arc degrees (which means 280.6 and 223 m, respectively).
The computational domain covers the area between 131 • and 140 • east longitude and between 30.01 • and 35 • north latitude.
According to the stability condition, the time step used in the computations was 0.5 s. Mostly, regular rectangular computational grids are used for tsunami numerical modelling. The grid could also be matched with one of a geography projection. We took into consideration the spherical shape of the Earth by decreasing the grid step with respect to longitude for larger values of latitude. Matching the grid to geographic coordinates, we suppose that the length of one arc degree at latitude λ is equal to the length of one equator arc degree multiplied by cos(λ).
Using the advantages of the proposed Calculator, it takes only 25 s to compute tsunami wave propagation over the water area in Figure 3. In practice, the travelling time of a real wave through this water area is 3200 s.
During the numerical experiments, we moved an artificial tsunami source along the Tokai-Nankai deep water trench. Note that, according to [1], a strong seismic event is expected in this area.
The authors discussed the disastrous tsunami that happened in this region nearly 300 years ago and provided a reconstruction of the tsunami source based on historical records (see Figure 4).

Parameters of the Artificial Tsunami Source
Tsunami sources of area 200 × 100 km presenting sea bed displacement due to an earthquake of magnitude 8.0 were used in our computational experiments. The suggested seismic mechanism is similar to the historical tsunamigenic earthquakes at the Nankai-Tokai subduction zone. The initial water surface displacement (taking the same as in [5]) with a maximum elevation of 400 cm and lowest subsidence of −200 cm caused by such an earthquake is presented in Figure 5. The model source was built according to the elastic-plastic model proposed in [17]. This model is currently used by the NOAA tsunami warning services to calculate the initial displacement in a tsunami source area on the basis of seismic monitoring data. The input data for calculating the ocean floor displacement field are the basic parameters describing the earthquake seismic mechanism. In this paper we used this displacement as the sea surface disturbance at the initial time moment. Figure 5. A 3D image of the model tsunami source used for numerical modeling; see [5].

Numerical Results
In this paper we checked the possibility of very rapid numerical simulation of tsunami wave propagation on a modern PC with FPGA-based acceleration. So, numerical simulation was performed up to the depth of 10 m. Full reflection boundary conditions were used at this depth. At the parts of the boundary that separate our computational domain from the ocean, we used free passage of the wave out of the domain. To obtain detailed wave height maxima along the coast, in future we will use a nested grid approach with a mesh resolution of several meters at coastal areas.

Parameters of the Artificial Tsunami Source
Tsunami sources of area 200 × 100 km presenting sea bed displacement due to an earthquake of magnitude 8.0 were used in our computational experiments. The suggested seismic mechanism is similar to the historical tsunamigenic earthquakes at the Nankai-Tokai subduction zone. The initial water surface displacement (taking the same as in [5]) with a maximum elevation of 400 cm and lowest subsidence of −200 cm caused by such an earthquake is presented in Figure 5. The model source was built according to the elastic-plastic model proposed in [17]. This model is currently used by the NOAA tsunami warning services to calculate the initial displacement in a tsunami source area on the basis of seismic monitoring data. The input data for calculating the ocean floor displacement field are the basic parameters describing the earthquake seismic mechanism. In this paper we used this displacement as the sea surface disturbance at the initial time moment.

Parameters of the Artificial Tsunami Source
Tsunami sources of area 200 × 100 km presenting sea bed displacement due to an earthquake of magnitude 8.0 were used in our computational experiments. The suggested seismic mechanism is similar to the historical tsunamigenic earthquakes at the Nankai-Tokai subduction zone. The initial water surface displacement (taking the same as in [5]) with a maximum elevation of 400 cm and lowest subsidence of −200 cm caused by such an earthquake is presented in Figure 5. The model source was built according to the elastic-plastic model proposed in [17]. This model is currently used by the NOAA tsunami warning services to calculate the initial displacement in a tsunami source area on the basis of seismic monitoring data. The input data for calculating the ocean floor displacement field are the basic parameters describing the earthquake seismic mechanism. In this paper we used this displacement as the sea surface disturbance at the initial time moment. Figure 5. A 3D image of the model tsunami source used for numerical modeling; see [5].

Numerical Results
In this paper we checked the possibility of very rapid numerical simulation of tsunami wave propagation on a modern PC with FPGA-based acceleration. So, numerical simulation was performed up to the depth of 10 m. Full reflection boundary conditions were used at this depth. At the parts of the boundary that separate our computational domain from the ocean, we used free passage of the wave out of the domain. To obtain detailed wave height maxima along the coast, in future we will use a nested grid approach with a mesh resolution of several meters at coastal areas.

Numerical Results
In this paper we checked the possibility of very rapid numerical simulation of tsunami wave propagation on a modern PC with FPGA-based acceleration. So, numerical simulation was performed Appl. Sci. 2020, 10, 4159 8 of 14 up to the depth of 10 m. Full reflection boundary conditions were used at this depth. At the parts of the boundary that separate our computational domain from the ocean, we used free passage of the wave out of the domain. To obtain detailed wave height maxima along the coast, in future we will use a nested grid approach with a mesh resolution of several meters at coastal areas.
We began our study with the given shape of initial sea surface displacement at the tsunami source (see Figure 5). Then, the shape of the sea surface and components of the velocity vector were calculated according to shallow water Equation (1) in the entire water area under consideration.
In Figure 6 is the distribution of the wave height maxima from the computed tsunami, generated by the model source S0-c (cf. Figure 3). Figure 7a,b shows two frames of the time history of the free surfaces near the southwest of Shikoku Island as a result of the tsunami generated by source S0-c.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 15 We began our study with the given shape of initial sea surface displacement at the tsunami source (see Figure 5). Then, the shape of the sea surface and components of the velocity vector were calculated according to shallow water Equation (1) in the entire water area under consideration.
In Figure 6 is the distribution of the wave height maxima from the computed tsunami, generated by the model source S0-c (cf. Figure 3). Figure 7a,b shows two frames of the time history of the free surfaces near the southwest of Shikoku Island as a result of the tsunami generated by source S0-c. In order to study the grid size sensitivity of the proposed method, the numerical experiment with source S0-c was conducted on a new computational grid. In the coarse mesh, we used half the number of points in each direction. The numerical results are almost the same. Figure 8 shows the quantitative difference (in meters) between the distributions of wave maxima for these two cases over the entire computational domain. Almost everywhere, the difference between the wave height maxima does not exceed 0.2 m. Only near the coast, in the area opposite the source, this difference reaches almost 1 m. We began our study with the given shape of initial sea surface displacement at the tsunami source (see Figure 5). Then, the shape of the sea surface and components of the velocity vector were calculated according to shallow water Equation (1) in the entire water area under consideration.
In Figure 6 is the distribution of the wave height maxima from the computed tsunami, generated by the model source S0-c (cf. Figure 3). Figure 7a,b shows two frames of the time history of the free surfaces near the southwest of Shikoku Island as a result of the tsunami generated by source S0-c. In order to study the grid size sensitivity of the proposed method, the numerical experiment with source S0-c was conducted on a new computational grid. In the coarse mesh, we used half the number of points in each direction. The numerical results are almost the same. Figure 8 shows the quantitative difference (in meters) between the distributions of wave maxima for these two cases over the entire computational domain. Almost everywhere, the difference between the wave height maxima does not exceed 0.2 m. Only near the coast, in the area opposite the source, this difference reaches almost 1 m. In order to study the grid size sensitivity of the proposed method, the numerical experiment with source S0-c was conducted on a new computational grid. In the coarse mesh, we used half the number of points in each direction. The numerical results are almost the same. Figure 8 shows the quantitative difference (in meters) between the distributions of wave maxima for these two cases over the entire computational domain. Almost everywhere, the difference between the wave height maxima does not exceed 0.2 m. Only near the coast, in the area opposite the source, this difference reaches almost 1 m.  Figures 6, 9, and 10. The color scale in Figures 6, 9, and 10 shows the height-color relation for wave height values between 0 and 5 m. In cases where the tsunami height exceeds 5 m at certain grid points, we used for visualization the same color as the points with maximum level elevation of 5 m. Figure 6 shows the tsunami height maxima over the entire computational domain in the case of the most distant offshore source position at the south of Shikoku Island (model source S0-c). In Figure 9, the results for the closest source position just opposite the center of Shikoku Island are presented (source S1-a). Figure 10 Figures 6 and 9, and Figure 10. The color scale in Figures 6 and 9, and Figure 10 shows the height-color relation for wave height values between 0 and 5 m. In cases where the tsunami height exceeds 5 m at certain grid points, we used for visualization the same color as the points with maximum level elevation of 5 m. Figure 6 shows the tsunami height maxima over the entire computational domain in the case of the most distant offshore source position at the south of Shikoku Island (model source S0-c). In Figure 9, the results for the closest source position just opposite the center of Shikoku Island are presented (source S1-a). Figure 10 shows the distribution of tsunami height maxima in the case of the intermediate tsunami source location near Kii Peninsula (source S3-b).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 15 Figure 9. Distribution of wave height maxima as the result of a tsunami wave generated by source S1-a (see Figure 3).
The computational grid in use had 3000 × 2496 nodes. According to our numerical results, the wave appears at the nearest dry land in 1000 s (approximately 16 min). The computation time to obtain the maximal wave height distribution in the entire water area is on the order of tens of seconds, provided that the proposed Calculator is applied to enhance the performance of a modern PC. When the Calculator is based on a Virtex-6 microchip, the results are obtained in 39.5 s, while in the case of Figure 9. Distribution of wave height maxima as the result of a tsunami wave generated by source S1-a (see Figure 3).  Figure 3).
The numerical simulation shows that the new (compared to the result given in Figure 9) position of the initial surface displacement of the same shape results in 6 m high tsunamis at Kii Peninsula, while all the Shikoku Island and the Kyushu Island coasts are practically safe with wave heights not exceeding 0.5 m. Figure 11A-D shows the computed distributions of tsunami wave maxima along the shore for all 12 model sources considered. We arranged the numerical experiments according to the same shape of initial sea surface displacement (given in Figure 5), varying the position of such model source along the shore and the distance of the model source from the shore. The colors of the wave maxima correspond to this distance. In Figure 11 we display in one picture the wave maxima distributions from three model sources located at different distances from the shore. Pictures A, B, C, and D of Figure 11 present the tsunami height distributions according to shore variation of the model source position. The computational grid in use had 3000 × 2496 nodes. According to our numerical results, the wave appears at the nearest dry land in 1000 s (approximately 16 min). The computation time to obtain the maximal wave height distribution in the entire water area is on the order of tens of seconds, provided that the proposed Calculator is applied to enhance the performance of a modern PC. When the Calculator is based on a Virtex-6 microchip, the results are obtained in 39.5 s, while in the case of a VC709 microchip, one needs only 7.85 s. So, the achieved acceleration compared to the tsunami travel time is nearly 127 times faster, and the achieved performance gain as compared to a regular PC is by approximately 380 times (see Table 1). Figures 6 and 9 show that for such locations of the initial sea surface displacement (sources S0-c and S1-a), the largest tsunami wave heights (5+ m) are obtained at Shikoku Island, while the coast of Kii Peninsula is practically safe as the expected tsunami wave heights do not exceed 1 m.
As follows from the numerical results, the given position of initial sea surface displacement results in an almost 8 m tsunami wave at a certain part of Shikoku Island (Figure 9). At the same time, the height of the tsunami wave is less than 1 m along the entire coastline of Kii Peninsula (see also points with horizontal index 1340-200 in Figure 11A). This is less than the height of the wave going toward the coast orthogonal to the source's long axis direction. Figure 10 presents the computed maximal heights of a tsunami wave from source S3-b with profile given in Figure 5, which is positioned at the Tonankai zone (the east-most source location among all regarded).
The numerical simulation shows that the new (compared to the result given in Figure 9) position of the initial surface displacement of the same shape results in 6 m high tsunamis at Kii Peninsula, while all the Shikoku Island and the Kyushu Island coasts are practically safe with wave heights not exceeding 0.5 m. Figure 11A-D shows the computed distributions of tsunami wave maxima along the shore for all 12 model sources considered. We arranged the numerical experiments according to the same shape of initial sea surface displacement (given in Figure 5), varying the position of such model source along the shore and the distance of the model source from the shore. The colors of the wave maxima correspond to this distance. In Figure 11 we display in one picture the wave maxima distributions from three model sources located at different distances from the shore. Pictures A, B, C, and D of Figure 11 present the tsunami height distributions according to shore variation of the model source position.
shape of initial sea surface displacement (given in Figure 5), varying the position of such model source along the shore and the distance of the model source from the shore. The colors of the wave maxima correspond to this distance. In Figure 11 we display in one picture the wave maxima distributions from three model sources located at different distances from the shore. Pictures A, B, C, and D of Figure 11 present the tsunami height distributions according to shore variation of the model source position.

Discussion
A unified computational grid for tsunami simulation is currently not available, so professionals have to apply many alternative solutions. For transoceanic events, the step size is usually comparable to 2 arc minutes (that is, nearly 3.6 km at the latitude of Japan). For inundation mapping, one can use a 5-10 m grid step. In other words, the choice of grid step should be correlated with the goals of the numerical modeling. In our case, we studied the distribution of tsunami maximal heights along the coastline for near-field events. So, we calculated wave propagation in a water area within an approximately 1000 × 2000 km rectangle. The simulation was carried on up to a depth of 10 m. We considered grid step sizes of 0.003 arc degrees in the east-west direction (nearly 280.6 m) and 0.002 degrees in the north-south direction (about 223 m), suitable for our purposes. The mesh points were matched with geographic coordinates. In particular, this means that the length of the east-west direction grid step slightly changes from the southern to northern edges of the considered water area. The Japanese bathymetric data produced by the Japan Oceanographic Data Center (JODC) and distributed by the Data On-line Service System (see [16] http://jdoss1.jodc.go.jp/cgibin/1997/depth500_file) were used for the grid development.
The main goal of our study was to provide fast and accurate evaluation of tsunami wave maximal heights along the coastline. We neglected fine analysis of the processes at depths less than

Discussion
A unified computational grid for tsunami simulation is currently not available, so professionals have to apply many alternative solutions. For transoceanic events, the step size is usually comparable to 2 arc minutes (that is, nearly 3.6 km at the latitude of Japan). For inundation mapping, one can use a 5-10 m grid step. In other words, the choice of grid step should be correlated with the goals of the numerical modeling. In our case, we studied the distribution of tsunami maximal heights along the coastline for near-field events. So, we calculated wave propagation in a water area within an approximately 1000 × 2000 km rectangle. The simulation was carried on up to a depth of 10 m. We considered grid step sizes of 0.003 arc degrees in the east-west direction (nearly 280.6 m) and 0.002 degrees in the north-south direction (about 223 m), suitable for our purposes. The mesh points were matched with geographic coordinates. In particular, this means that the length of the east-west direction grid step slightly changes from the southern to northern edges of the considered water area. The Japanese bathymetric data produced by the Japan Oceanographic Data Center (JODC) and distributed by the Data On-line Service System (see [16] http://jdoss1.jodc.go.jp/cgi-bin/1997/depth500_ file) were used for the grid development.
The main goal of our study was to provide fast and accurate evaluation of tsunami wave maximal heights along the coastline. We neglected fine analysis of the processes at depths less than 10 m, even though such analysis is very important for mapping possible flooding areas. To our understanding, the computed distribution of the wave heights provides the basis for such research, and we will arrange the corresponding studies in the future. Note also that, according to [18,19], it is possible to use formulae for recalculating the obtained tsunami wave height at the vertical wall at a depth of 10-20 m for rough determination of the long wave run-up height. Such a coefficient is predetermined for each coast point and is obtained by tsunami numerical modeling using more accurate (in some cases, three-dimensional) hydrodynamic models.
We note that there is no information for fine-resolution (compared to 10 m) grids in the global bathymetry databases. In case one would like to refine the resolution by interpolation, fine details of the real depths could just be missed. This may lead to the divergence of numerical simulation results from reality. The nested grid approach, where the mesh step depends on the distance from the shore, is able to resolve this problem-provided that the distribution of depths is known with the required precision. Hardware code execution acceleration (with the help of "graphics cards"-GPUs) for a nested grid was previously used by the authors for modeling tsunami in the Tohoku/Fukushima area of Japan [20].
Pictures A, B, C, and D of Figure 11 present the tsunami height distribution according to the alongshore variation of the model source position. The longitude position of the coastal grid points in these distributions is counted along the horizontal axis from 1 at the left boundary up to 3000 on the right (see the scale at the bottom of Figure 3).
Analysis of the distribution of wave height maxima along the coast leads to the following conclusions. Close to the coast sources Si-a and intermediately posed sources Si-b give higher tsunami heights when compared to the sources Si-c, which are most remote from the coast. However, the latter sources give higher average levels of tsunami heights outside of the projection of the source to the coastline. Out of these areas, the wave height generated by the remote sources Si-c (drawn in Figure 11A-D in blue color) exceeds the wave height from sources Si-a and Si-b (i = 0, . . . , 3) closer to the coast. This is expected according to the ray theory of long wave propagation. Also, it can be explained by the fact that the source more distant from the coast is situated at a greater depth, and according to Green's law [15], this source provides bigger tsunami height growth approaching the coast. Sharp peaks in the distributions of height maxima of a tsunami are caused, as a rule, by local features of the coastal topography, including the bottom and the configuration of the coastline. In particular, the underwater ridges coming out to some capes work as waveguides, sharply increasing the amplitude of waves reaching these points of the coast. This is clearly visible in Figures 9 and 10. Note that unexpectedly high tsunami wave heights are observed in the narrowing gulfs and bays. This is clearly visible in Figure 11D at the east part of Kii Peninsula for tsunami sources S3-a, S3-b, and S3-c. The horizontal indices of these points are in the range from 1600 to 2000 (see the horizontal scale in Figure 3).
The results of our numerical experiments show that rather small changes to the particular position of a tsunami source of the same shape may result in a dramatic difference in the tsunami wave height maxima distribution along the coast. For example, the source positions in Figure 11A,B differ by just 200 km (size of the source area). However, the maximal wave heights for the case in Figure 11B are notably smaller and are located at different points along the coast when compared to those for the case in Figure 11A. Thus, fast and reliable simulation results are needed as soon as possible after a near-shore seismic event to suggest adequate evacuation measures.

Conclusions
The specific hardware configuration and corresponding modeling code we previously proposed were used to study the distribution of tsunami wave height maxima. Only the wave propagation over sufficiently large depth (10 m and more) was considered. Therefore, near-coast phenomena like dispersive effects, wave breaks, and others were not accounted for. Certainly, these are the limitations of the presented approach. A number of numerical experiments were conducted to understand the dependence on the source location of the tsunami wave maximal height distribution along the coastline. The so-called MacCormack finite difference version of the shallow water approximation of equations of hydrodynamics was applied. As observed, even a 200 km shift in artificial tsunami source position dramatically changes the distribution of tsunami wave maxima along the coast.
A special tool we previously proposed for hardware acceleration of computer code execution on a PC was applied. The achieved performance gain was on the order of several hundred times in the case of a VC709 FPGA microchip versus the same PC without hardware acceleration. This means that one needs just a few seconds to simulate tsunami wave propagation from a source in the Nankai-Tonankai (or another) subduction zone to the coast. The mesh we used contained 3000 × 2496 nodes. Such high performance shows the availability of "real-time" tsunami danger forecast in the very near future.
The proposed software tools and hardware code acceleration devices should be used together with the wave generation models and special approaches for inundation mapping. We note that the suggested tools should be used to study tsunami only at the regional scale, for computation areas within a few geographical degrees. For these situations, the sphericity of the Earth is not able to have a notable impact on the tsunami wave parameters. In forthcoming studies, we will present an approach to determine the initial sea surface disturbance at the tsunami source and an implementation of the code acceleration on nested grids. These two steps will contribute to overcoming the aforementioned limitations of our approach. However, even now, the proposed tool can be used as an input data generator for more precise inundation modeling.

Conflicts of Interest:
The authors declare no conflict of interest.