CC BY
6
DOI: 10.17516/1997-1397-2021-14-4-433-444 УДК 519.684.4, 550.394.2
Fast Modelling of Tsunami Wave Propagation using PC with Hardware Computer Code Acceleration
Mikhail M. Lavrentiev*
Institute of Automation and Electrometry SB RAS Novosibirsk, Russian Federation
Andrey G. Marchuk
Institute of Computational Mathematics and Mathematical Geophysics SB RAS
Novosibirsk, Russian Federation
Received 10.03.2021, received in revised form 05.04.2021, accepted 20.05.2021 Abstract. The field programmable gates array (FPGA) microchip is applied to achieve considerable performance gain in simulation of tsunami wave propagation using personal computer. The two-step Mac-Cormack scheme was used for approximation of the shallow water equations. An idea of PC-based tsunami wave propagation simulation is described. Comparison with the available analytic solutions and numerical results obtained with the reference code show that developed approach provides good accuracy in simulations. It takes less then 1 minute to compute 1 hour of the wave propagation in computational domain that contains 3000 x 2500 nodes. Using the nested greed approach, it is possible to decrease the size of space step from about 300 meters to 10 m. Using the proposed approach, the entire computational process (to calculate the wave propagation from the source area to the coast) takes about 2 min. As an example the distribution of maximal heights of tsunami wave along the coast of the Southern part of Japan is simulated. In particular, the interrelation between maximal wave heights and location of tsunami source is studied. Model sources of size 100 x 200 km have realistic parameters for this region. It was found that only selected parts of the entire coast line are exposed to tsunami wave with dangerous height. However, the occurrence of extreme tsunami wave heights at some of those areas can be attributed to the local bathymetry. The proposed hardware acceleration to compute tsunami wave propagation can be used for rapid (say, during few minutes) evaluation of danger from tsunami wave for a particular location of the coast.
Keywords: numerical modelling, tsunami wave propagation, computer code acceleration. Citation: M.M. Lavrentiev, A.G.Marchuk, Fast Modelling of Tsunami Wave Propagation using PC with Hardware Computer Code Acceleration, J. Sib. Fed. Univ. Math. Phys., 2021, 14(4), 433-444. DOI: 10.17516/1997-1397-2021-14-4-433-444.
Early warning of dangerous tsunami waves at a particular coastal location is crucially important to reduce human losses and minimize possible impact on economy. Unfortunately, the problem of tsunami early warning after the major offshore earthquake is still unresolved, despite the rather large number of publications on this issue (see, for example, [1]). In the case of the seismic event offshore Japan, tsunami wave approaches the nearest point at the coast in approximately 20 minutes. It means that just a few minutes are available for the analysis to provide the authorities with evaluation of the expected tsunami wave danger. In the case of the strong earthquake the electric power supply may be disrupted, so it would be better not
* mmlavrentiev@gmail.com tmag@omzg.sscc.ru © Siberian Federal University. All rights reserved
to use supercomputer facilities. Advantages of the modern computer architectures to accelerate numerical simulation of tsunami wave propagation can be used. An approach based on the specialized FPGA (Field Programmable Gates Arrays) have been developed and tested. A number of numerical tests demonstrate the advantages of the new fast method for modelling tsunami wave propagation. In the present paper we briefly describe the proposed approach and show the numerical results obtained in the recent years.
Our point is to use just a personal computer (PC) and to achieve performance gain using Graphic Processing Units (GPUs) and FPGAs as co-processors. So, we propose a specific hardware configuration and the corresponding code.
Robust evaluation of tsunami wave danger should be based on the correct process simulation: wave generation, wave propagation, and inundation of dry land. In the study we deal with the stages of wave formation and propagation only to decrease computation time and keep at the same time sufficient accuracy. So, we suppose that the tsunami wave is caused by a certain disturbance of sea surface. From this initial disturbance follows initial conditions for the governing evolution type equations (shallow water equations). The issue of the inundation mapping was not considered. Therefore, we do not compute waves when depth is small (below 5 m) where reflection type boundary conditions are suggested at such depth to estimate the wave height in the near-shore area and to account for reflected waves. On the parts of the boundary which separate our computational domain from the ocean, conditions for free passage of the wave out of the domain are used. There are several application programs to simulate the wave propagation over the real digital bathymetry [2-6].
Among the most popular programs the MOST (Method of Splitting Tsunamis) package of programs should be mentioned. This package is used by the USA NOAA tsunami warning centres to simulate all tsunami phases - generation, propagation, and inundation of the dry land [2,4]. Simulation of the wave propagation over chosen water area is based on the numerical solution of linear or non-linear shallow water differential equations.
Alternatively, the Mac-Cormack scheme for numerical approximation of the shallow water equations was implemented [7,8]. Comparison with the exact solutions (in special cases of sea bed relief) shows a very good accuracy of the implemented method [9,10]. It shows better tracking of the wave front in comparison with the MOST program.
A number of numerical tests where real digital bathymetry of the offshore Japan and Kamchatka Peninsula was used prove that it takes about 50 sec to solve numerically the shallow water equations for the computational domain with 107 nodes [11,12]. Nested grid approach was also tested [13,14].
1. Formulation of the problem
The referred program MOST (like many other tools) uses the following equivalent form of the shallow water equations which does not take into account such external forces as sea bed friction, Coriolis force and others [4]:
dH duH dvH _ 0 dt dx dy ' du + du + du + dH dD dt dx dy dx dx ' dv dv dv dH dD dt dx dy dy dy
where H(x, y, t) = n(x, y, t) + D(x, y, t) is the entire height of water column, n(x, y, t) is the sea surface disturbance (wave height), D(x,y) — depth (which is supposed to be known at all grid points), u and v are components of flow velocity vector, g — acceleration of gravity.
The system of shallow-water equations can be solved using the difference scheme. In this case values of tsunami wave parameters n(x,y,t), u, and v are defined in nodes of the regular grid linked to geographical coordinates. To begin with, the initial conditions are set in all grid nodes. For example, everywhere besides the source area, the values of the grid variables n0j, U0, vj, (i = 1,... ,N, j = 1,..., M) are equal to zero. Further, according to the difference scheme that approximates system of differential equations (1), the wave parameters nnj, un, vn at the grid nodes on subsequent time layers tn = nr are calculated. Here, the value of the time step t is usually determined from the stability condition. This condition requires that wave can not travel more than one spatial step (Ax or Ay) in one time step.
Shallow water equations (1) are approximated at the grid nodes on the n-th time step with the help of explicit two-step Mac-Cormack finite difference scheme of the second order of approximation [15]:
First half-step:
n+1 TTn un„,n Tjn „,n Un „,n Tjn „,n
ij - Hij , HijUij - Hi-1jUi-1j , Hijvij - Hij-1vij-1
+ j j_i—ij + j j_j— 1 j—1 _o
t Ax Ay '
f,n+1 - un un — un un — un nn — nn
aij uij , n uij ai- 1j , n uij u'ij-1 , 'lij 'k- 1j „ /0\
-+ uij Ax + vij +g_ 0 (2)
j - v^nj , ^n vij- - vn-ij , vn vij- - vj-1 , - j ij
+ vn.-j_—j + vn.-j_j-1 + g-^j_j-1 _ 0
+ Uij Ax + vj Ay + g Ay _0-
Second half-step:
j - {Èn + 1 + Hn )/2 + Èn++1 Un+1 - Èn + j j - H j'y j1 _o
t/2 + Ax + Ay '
n+1 ( ^n+1 I n \ /o ^n+1 ^ n+1 ^n+1 ^n+1 ^n+1 ^n+1
Uij - (Uij + Uij)/2 + Un Ui+1j - Uij + vn Uij+1 - Uij + 3+1j - Vij _0 (3)
t/2 +Uij Ax +vij Ay +g Ax _0, (3)
vIT1 - ( j+vij)/2 + „ *r++j -j + vn j -j + j -j _0
+ Uij \ + vij A + g A _ 0.
"j" " I un/i+1j "j + vn/ij+1 "j + g' lij+1 n i' ij Ax ij Ay Ay
Here Fij+1 are intermediate values of wave parameters after the first time half-step.
Usually, the real tsunami wave simulation is performed in the spherical or geodetic coordinate system (X,4), where A is the longitude and $ is the latitude in arc degrees. Accordingly, the following relations are used to calculate the differences Ax and Ay:
n(Ai+1 - Ai) n($i+1 - u
Axij = -Rr cos(ffli), Ayij =-RE,
ij 180° E KYih yij 180°
where RE is the Earth radius. This scheme looks similar to the splitting method (with respect to space variables) which is used in the referred MOST program package. Indeed, in order to calculate the values of the sought functions at grid-point (i,j,n + 1) the values at 3 points of the previous time step, (i,j,n), (i — 1,j,n), and (i,j — 1,n) are used during the first half-step in (2), and the values at the points (i,j,n), (i + 1,j,n), and (i,j + 1,n) are used in the second half-step in (3). Comparison of the known analytic solutions with the numerical solutions shows that the proposed attempt to realize the three-points calculation stencil (Mac-Cormack scheme) seems to be preferable compared to the one from the MOST software package [9,10]
2. FPGA based calculator
In order to achieve performance gain, the FPGA-based Calculator has been developed.
To employ advantages of the FPGA microchip features, the stream processor architecture was proposed for this algorithm implementation. The proposed Calculator contains several processor elements (PEs). Each elements performs a pipeline with a sequential data stream. "On board" memory contains all the necessary information. The calculation speed-up by FPGA architecture is based on the inner memory (BRAM) access for implementing stencil buffer.
The Calculator architecture allows one 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 can be optimized considering features of the FPGA microchip in use. The Mac-Cormack finite difference approximation fits very well with the Calculator architecture presented in Fig. 1, processing 1 node at one computer clock cycle.
Calculator based of FPGA microchip Xilinx Virtex-7 VC709 was used for numerical tests (see [7,8] for details).
Fig. 1. Calculator architecture. To implement the FPGA algorithm the following architecture of the stream processor was proposed. It consists of processor elements (PE). Such PE executes a version of 2-dimension run, a pipeline with a sequential data stream. In addition to the calculator itself, the processor has memory controllers DDR3, PCIe controller, and DMA module responsible for the interaction between the calculator and the memory of the host computer. Such interaction is arranged as a direct memory access (DMA)
3. Comparison with exact solutions and reference code
In order to verify results obtained by the described method, a number of numerical tests have been carried out. The first test consists of simulation of the tsunami wave propagation from a round source in the area with the sloping bottom topography. The water area of 1000 x 1000 km was considered. Computational grid has equal spatial steps in both directions, namely, Ax = Ay = 1000 m. The centre of the circular tsunami source with the radius of 50 km was located in the middle of the region (horizontally) 300 km from the lower boundary where the depth was vanishing. The depth is linearly increased according to the formula D(x,y) = 0.01y, where y is the distance from the lower boundary of the region. The initial vertical displacement h of the tsunami source is determined by the formula
f nr\ Vro/ '
h(r) = 1 + cos ( — ) ' 0 < r < r0. (4)
Here, r is the distance to the centre of the source of radius r0. Thus, in the centre of the source area the initial displacement of the water surface was +2 m. This source generates a circular wave with the height of 0.95 m at 50 km from the centre. It is the wave height that was used as initial circular wave front with the radius of 50 km to estimate the wave amplitude at all points of the region according to the ray approximation [17,18]. This distribution of tsunami wave amplitude (given in the analytic form) was compared with the distribution obtained with the use of the MOST program and the proposed Mac-Cormack algorithm (Fig. 2).
Fig. 2. Left: Isolines of maximal wave heights distribution above bottom slope: exact solution of shallow water equations [17,18] (brown lines), numerical solution with the FPGA based Calculator (red dashed lines), the MOST program (blue lines). The offshore distance is measured along vertical axis relative to the figure bottom boundary. Right: Isolines of maximal wave heights distribution above the parabolic bottom topography [16]: the wave-ray solution to shallow water equations (brown lines), numerical solution with the FPGA based Calculator (dashed lines), the MOST program (blue lines). The offshore distance is measured along vertical axis relative to the figure bottom boundary
Fig. 2 (left) shows that at sufficiently large depths (exceeding 500 m) the contours of all three distributions of tsunami height maxima are quite close to each other. The proximity of results of numerical calculations for the two algorithms under consideration is also preserved near the coast. The discrepancy between numerical and wave-ray approach results in the coastal zone is caused by the neglect of the effect of increasing wave height due to reflection from the coast.
Another numerical test is similar to the first one and considers the case of the parabolic bottom relief. Let us consider the same computational area of 1000 x 1000 km with the computational grid having spatial steps Ax = Ay = 1000 m. The centre of the circular source with 50 km radius is also located in the middle of the area at 300 km from the lower boundary where the depth is equal to zero. The depth increased according to the formula D(x,y) = 10~8y2, where y is the distance from the lower boundary of the region. The initial vertical displacement inside the circular source is determined by (4). Fig. 2 (right) presents the isolines of distributions of tsunami height maxima calculated by the MOST program and the Mac-Cormack algorithm together with the estimates of these maxima obtained in the framework of the ray model [16].
Fig. 2 (right) shows that at sufficiently large depths (more than 200 m) the contours of all three distributions of tsunami height maxima are quite close to each other. Here, the similarity of the results of numerical simulations with two algorithms under consideration is observed up
to the coastline (the lower boundary of the region). The discrepancy between numerical and wave-ray approximation results is increased in the coastal zone.
Based on the results of the wave propagation over a bottom slope the correctness of the wave front kinematics modelling is also estimated. Fig. 3 shows the comparison of the wave front position with an interval of 5 minutes obtained with the McCormack scheme with the exact solution of the kinematic problem at the same time points. In order to prevent the points from merging (as it happens at the initial moment) the moments of output of the points of the calculated wave front are taken 3 seconds later than the moment of the corresponding exact solution of the kinematic problem [17,18]. The wave front positions obtained with the MOST program are not presented since these positions of the front points exactly coincide with the positions obtained with the McCormack scheme.
Fig. 3. Comparison of the tsunami wave isochrones over the sloping bottom: numerical experiments with the MOST and Mac-Cormack algorithms (grey lines) and exact solution (blue points)
Additional numerical test was carried out in order to verify the correctness of modelling the reflection of wave from a completely reflecting boundary positioned at 45 degrees to the direction of motion of the flat wave front. Let us consider rectangular 1000 x 2000 nodes computational domain. A long wave about 1 m height generated by a one-dimensional source parallel to the lower boundary of the region propagates over the region and it is reflected from the inner boundary positioned at 45 degrees to ordinate axis (Fig. 4).
Fig. 4 shows the distribution of the vertical displacement of the water surface in the entire region calculated with the Mac-Cormack scheme (left figure) and the MOST algorithm (right figure). The dark line shows the tsunami height isoline corresponding to the value of 0.4 m. The grey line outlines the water area with the surface displacement less then —0,4 m. Both
Fig. 4. Water surface after 3,000 sec of wave propagation using the Mac-Cormack scheme with the FPGA (left figure) and the MOST program (right figure). Thin black line indicates the 0.4 m isoline, and the gray line shows —0.4 m isoline
figures confirm the correctness of numerical modelling of the process of wave reflection from a completely reflecting boundary. One can see that the direction of motion of the reflected wave in both cases is orthogonal to the direction of the incident wave.
4. Distribution of maximal wave heights along the shoreline
The acceleration of numerical calculations of the tsunami propagation is required, first of all, by tsunami warning services to fast estimate the expected wave height at various points on the coastline. Therefore, this estimation is required before tsunami wave reaches the shore. The ability of the proposed approach to solve this problem in the area with real bathymetry within a few tens of seconds is demonstrated in this section.
The series of numerical experiments were performed for the areas around Kii Peninsula and Shikoku Island (southern part of Japan). Japanese bathymetric data produced by the Japan Oceanographic Data Center (JODC) (see [19]) were used, and they are presented in Fig. 5.
The above bathymetry and the computational grid have the following characteristics:
(1) Computational domain contains 3000 x 2496 nodes; (2) Grid steps are 0.003 and 0.002 degrees (which means 280.6 and 223 meters, respectively); (3) The grid covers the area between 131° and 140° E, 30.01° and 35° N; (4) Time step used in computations is equal to 0.5 sec.
The shape of model tsunami sources used in numerical experiments are based on the available geological and geophysical information. The typical for subduction zone seabed displacement area
Fig. 5. Digital bathymetry around Kii Peninsula and Shikoku Island (Japan). Positions of model tsunami sources are indicated
for 8.0 M earthquake was approximated by 100 x 200 km rectangle having maximum height 300 cm. The initial seabed displacement for the model source is shown in Fig. 6.
Fig. 6. 3D image of the model tsunami source used in numerical simulation
Positions and shapes of the sources used for tsunami modelling are shown in Fig. 5, where geography of the computational domain is also presented. In this figure only the closest to the coast sources Si — a (i = I,... 4) and most distant to the coast sources Si — c (i = !,... 4) are indicated. Their positive parts are shaped by pink colour, and negative parts (water surface depression) are outlined by yellow colour. Intermediate sources Si — b (i = I,... 4) are located between sources Si — a and Si — c.
The distributions of the wave height maxima in the entire area generated by some model sources are presented in Fig. 7. In the right part of each drawing the legend for colour-height relation is presented.
As is observed from Fig. 7 (left), the same shape of the initial sea surface displacement causes the tsunami wave heights of up to 6 m at certain areas of the Shikoku Island and the Kyushu Island coasts. At the same time, the Kii Peninsula coast is practically safe with wave heights limited by 0.5 m. Otherwise, the source located opposite Kii peninsula seriously affects only its coast and it has no effect on Shikoku Island (see Fig. 7 (right)).
Distribution of tsunami wave maxima along the shoreline is also important information for
Fig. 7. Calculated height maxima in entire computational domain from tsunami wave generated by the source S0 — c (left) and S3 — b (right) (see Fig. 3)
tsunami warning. Figs. 8(A) and 8(D) show such distribution for 6 model sources S0 — a,b,c and S3 — a, b, c. Numerical experiments were performed for the same shape of the initial see surface displacement (given in Fig. 6). Positions of the model source along the shore and the distance of the model source from the shore were varied. Numbers along the horizontal axis indicate the horizontal indexes of coastal computational grid points.
Fig. 8. Distributions of tsunami wave maxima along the shore generated by the sources S0 — a, b, c (A) and S3 — a,b,c (D). Wave height maxima: yellow lines — Si — a sources; blue lines — Si — b sources; pink lines — Si — c sources
Let us say few words about digital bathymetry. The grid step 250 m used in our numerical
experiments is considered too large these days. However, the size of grid step depends on the goal of simulation. If a detailed evaluation of the expected wave heights along the entire shore line is needed then it is necessary to carry out numerical simulation with the corresponding fine mesh size in the near-coastal regions. It can be done by choosing the sufficiently small grid step in the whole area. However, the number of computational nodes is increased by 2-3 orders. It results in the necessity to extend computational facilities or, alternatively, in a dramatic increase of the CPU time required for simulation. As for the proposed FPGA-based Tsunami Wave Calculator, the available memory resources permits the use of approximately 50 millions computational nodes.
The Calculator with a regular modern PC needs 25 sec to simulate wave propagation from the southern edge of the computational domain shown in Fig. 5 to the shore. The estimated travel time for the wave is 3200 sec. So, realization of Mac-Cormack scheme on FPGA hardware allows one to estimate the expected wave height distribution along the coastline before tsunami arrival.
Conclusion
In order to accelerate the calculation of tsunami wave propagation over the deep water area, a special FPGA based Calculator has been developed. The Mac-Cormack scheme was used for numerical solution of the shallow water equations. Accuracy of the solution obtained with the use of the Calculator was tested by comparison with the known analytic solutions. Similar or even better accuracy is achieved in comparison with the MOST program. These results show the possibility of tsunami danger forecast in the real time mode.
This work was supported by ICMMG SB RAS (state contract 0315-2019-0005) and by IAE SB RAS.
References
[1] H.Tsushima, Y.Yusaku, Review on Near-Field Tsunami Forecasting from Offshore Tsunami Data and Onshore GNSS Data for Tsunami Early Warning, J. Disaster Res., 9(2014), no. 3, 339-357.
[2] E.Gica, M.Spillane, V.Titov, C.Chamberlin, J.Newman, Development of the forecast propagation database for NOAA's short-term inundation forecast for tsunamis (SIFT), NOAA Technical Memorandum, 2008. URL: http://www.ndbc.noaa.gov/dart.shtml (access date: 15.06.2016).
[3] H.Kensaku, A.P.Vazhenin, A.G.Marchuk, Trans-Boundary Realization of the Nested-Grid Method for Tsunami Propagation Modeling, Proceedings of the Twenty-fifth (2015) International Ocean and Polar Engineering Conference Kona, Big Island, Hawaii, USA, 3(2015), 741-747.
[4] V.Titov, F.Gonzalez, Implementation and testing of the method of splitting tsunami (MOST) model, NOAA Technical Memorandum ERL PMEL-112, 1997.
[5] N.Shuto, C.Goto, F.Imamura, Numerical simulation as a means of warning for near field tsunamis, Coastal Engineering in Japan, 33(1990), no. 2, 173-193.
[6] A.C.Yalciner, B.Alpar, Y.Altinok, I.Ozbay, F.Imamura, Tsunamis in the Sea of Marmara: Historical Documents for the Past, Models for Future, Marine Geology, 190(2002), 445-463.
[7] M.M.Lavrentiev, A.A.Romanenko, K.K.Oblaukhov, An.G.Marchuk, K.F.Lysakov, M.Yu.Shadrin, FPGA Based Solution for Fast Tsunami Wave Propagation Modeling, The 27th International Ocean and Polar Engineering Conference, (2017), 25-30 June, San Francisco, California, USA, 924-929.
[8] M.M.Lavrentiev, K.F.Lysakov, An.G.Marchuk, K.K.Oblaukhov, M.Yu.Shadrin, Hardware Acceleration of Tsunami Wave Propagation Modeling in the Southern Part of Japan, Appl. Sci., 10(2020), no. 12, 4159. DOI:10.3390/app10124159
[9] M.M.Lavrentiev, An.G.Marchuk, K.K.Oblaukhov A.A.Romanenko, Comparative testing of MOST and Mac-Cormack numerical schemes to calculate tsunami wave propagation, J. Phys.: Conf. Ser, 1666(2020), 012028. DOI: 10.1088/1742-6596/1666/1/012028
[10] M.M.Lavrentiev, A.A.Romanenko, K.K.Oblaukhov, An.G.Marchuk, K.F.Lysakov, M.Yu.Shadrin, Implementation of Mac-Cormack scheme for the fast calculation of tsunami wave propagation, Oceans'17 MTS/IEEE, Aberdeen, June 19-22, 2017.
[11] M.M.Lavrentiev, K.F.Lysakov, An.G.Marchuk, K.K.Oblaukhov, and M.Yu.Shadrin, Fast evaluation of tsunami waves heights around Kamchanka and Kuril Islands, Science of Tsunami Hazards, 38(2019), no. 1, 1-13.
[12] M.Lavrentiev, K.Lysakov, An.Marchuk, K.Oblaukhov, M.Shadrin, FPGA-based Modelling of the Tsunami Wave Propagation at South Japan Water Area, Proc. OCEANS'18 MTS/IEEE Kobe, Japan, May 28-31, 2018.
[13] M.M.Lavrentiev, K.F.Lysakov, An.G.Marchuk, K.K.Oblaukhov, M.Yu.Shadrin, FPGA Based Tsunami Wave Propagation Calculator, J. Phys.: Conf. Ser., 1789(2021), 012011. DOI: 10.1088/1742-6596/1789/1/012011
[14] M.M.Lavrentiev, K.F.Lysakov, An.G.Marchuk, K.K.Oblaukhov, M.Yu.Shadrin, FPGA based modeling of Tohoku tsunami using nested grids, Proceedings of the Global Oceans 2020: Singapore - U.S. Gulf Coast, October 5-14, 2020.
[15] R.W.MacCormack, A.J.Paullay, Computational Efficiency Achieved by Time Splitting of Finite-Difference Operators, AIAA paper, (1972), 72-154.
[16] An.G.Marchuk, Benchmark solutions for tsunami wave fronts and rays. Part 2: Parabolic bottom topography, Science of Tsunami Hazards, 36(2017), no. 2, 70-85.
[17] An.G.Marchuk, Estimating Tsunami Wave Height over a Sloping Bottom in the Ray Approximation, Numerical Analysis and Applications, 8(2015), no. 4, 304-313.
DOI: 10.1134/S1995423915040047
[18] An.G.Marchuk Benchmark solutions for tsunami wave fronts and rays. Part 1: sloping bottom topography, Science of Tsunami Hazards, 35(2016), no. 2, 34-48.
[19] https://jdoss1.jodc.go.jp/vpage/depth500_file.html
Быстрое моделирование распространения волны цунами на ПК за счет аппаратного ускорения исполнения кода
Михаил М. Лаврентьев
Институт автоматики и электрометрии СО РАН Новосибирск, Российская Федерация
Андрей Г. Марчук
Институт вычислительной математики и математической Новосибирск, Российская Федерация
Аннотация. За счет применения микросхемы вентильной матрицы, программируемой пользователем (Field Programmable Gates Array - FPGA), достигается значительное ускорение расчета распространения волн цунами на современном обычном персональном компьютере. Для численной аппроксимации системы уравнений мелкой воды использовалась двухступенчатая схема Мак-Кормака. На базе проведенных численных тестов авторы описывают идею моделирования распространения волн цунами на базе ПК. Проведенное сравнение с известными аналитическими решениями и с эталонным кодом показывает хорошую точность разработанного программного приложения. Расчет одного часа распространения волны занимает меньше 1 минуты на сетке 3000 х 2500 узлов. Используя технологию вложенных сеток, можно перейти от расчетной сетки с шагом 300 м до сетки с шагом 10 м. При использовании предложенного Калькулятора, весь вычислительный процесс (для расчета распространение волны от очага до берега) занимает около 2 мин. Получено распределение максимальных высот волн цунами вдоль побережья южной части Японии. В частности, исследуется зависимость максимальных высот волн от конкретного местоположения источника цунами. Модельный источник размером 100 х 200 км имеет реалистичные параметры для этого географического региона. Результаты численных экспериментов показывают, что только на отдельных участках всей береговой линии наблюдаются опасные амплитуды волн цунами. Наличие аномально высоких волн цунами в некоторых из этих районов могут быть вызваны особенностями локальной батиметрии. Предлагаемое аппаратное ускорение вычисления распространения волны цунами может быть использовано для быстрой (скажем, за несколько минут) оценки опасности цунами для конкретного населенного пункта или промышленного объекта на побережье.
Ключевые слова: численное моделирование, распространение волны цунами, ускорение исполнения программного кода.
