High-performance tsunami modelling with modern GPU technology

PhD ThesisEarthquake-induced tsunamis commonly propagate in the deep ocean as long waves and develop into sharp-fronted surges moving rapidly coastward, which may be effectively simulated by hydrodynamic models solving the nonlinear shallow water equations (SWEs). Tsunamis can cause substantial economic and human losses, which could be mitigated through early warning systems given efficient and accurate modelling. Most existing tsunami models require long simulation times for real-world applications. This thesis presents a graphics processing unit (GPU) accelerated finite volume hydrodynamic model using the compute unified device architecture (CUDA) for computationally efficient tsunami simulations. Compared with a standard PC, the model is able to reduce run-time by a factor of > 40. The validated model is used to reproduce the 2011 Japan tsunami. Two source models were tested, one based on tsunami waveform inversion and another using deep-ocean tsunameters. Vertical sea surface displacement is computed by the Okada model, assuming instantaneous sea-floor deformation. Both source models can reproduce the wave propagation at offshore and nearshore gauges, but the tsunameter-based model better simulates the first wave amplitude. Effects of grid resolutions between 450-3600 m, slope limiters, and numerical accuracy are also investigated for the simulation of the 2011 Japan tsunami. Grid resolutions of 1-2 km perform well with a proper limiter; the Sweby limiter is optimal for coarser resolutions, recovers wave peaks better than minmod, and is more numerically stable than Superbee. One hour of tsunami propagation can be predicted in 50 times on a regular low-cost PC-hosted GPU, compared to a single CPU. For 450 m resolution on a larger-memory server-hosted GPU, performance increased by ~70 times. Finally, two adaptive mesh refinement (AMR) techniques including simplified dynamic adaptive grids on CPU and a static adaptive grid on GPU are introduced to provide multi-scale simulations. Both can reduce run-time by ~3 times while maintaining acceptable accuracy. The proposed computationally-efficient tsunami model is expected to provide a new practical tool for tsunami modelling for different purposes, including real-time warning, evacuation planning, risk management and city planning

PDE-Based Multidimensional Extrapolation of Scalar Fields over Interfaces with Kinks and High Curvatures

We present a PDE-based approach for the multidimensional extrapolation of smooth scalar quantities across interfaces with kinks and regions of high curvature. Unlike the commonly used method of [2] in which normal derivatives are extrapolated, the proposed approach is based on the extrapolation and weighting of Cartesian derivatives. As a result, second- and third-order accurate extensions in the $L^\infty$ norm are obtained with linear and quadratic extrapolations, respectively, even in the presence of sharp geometric features. The accuracy of the method is demonstrated on a number of examples in two and three spatial dimensions and compared to the approach of [2]. The importance of accurate extrapolation near sharp geometric features is highlighted on an example of solving the diffusion equation on evolving domains.Comment: 17 pages, 13 figures, submitted to SIAM Journal of Scientific Computin

S-wave splitting in the offshore South Island, New Zealand : insights into plate-boundary deformation

Author Posting. © American Geophysical Union, 2015. This article is posted here by permission of American Geophysical Union for personal use, not for redistribution. The definitive version was published in Geochemistry, Geophysics, Geosystems 16 (2015): 2829–2847, doi:10.1002/2015GC005882.Local and regional S-wave splitting in the offshore South Island of the New Zealand plate-boundary zone provides constraints on the spatial and depth extent of the anisotropic structure with an enhanced resolution relative to land-based and SKS studies. The combined analysis of offshore and land measurements using splitting tomography suggests plate-boundary shear dominates in the central and northern South Island. The width of this shear zone in the central South Island is about 200 km, but is complicated by stress-controlled anisotropy at shallow levels. In northern South Island, a broader (>200 km) zone of plate-boundary parallel anisotropy is associated with the transitional faulting between the Alpine fault and Hikurangi subduction and the Hikurangi subduction zone itself. These results suggest S-phases of deep events (∼90 km) in the central South Island are sensitive to plate-boundary derived NE-SW aligned anisotropic media in the upper-lithosphere, supporting a “thin viscous sheet” deformation model.United States National Ocean Bottom Seismograph Instrumentation Pool2016-02-2

Adaptive mesh refinement method. Part 2: Application to tsunamis propagation

Numerical simulations of multi dimensional large scale fluid-flows such as tsunamis, are still nowadays a challenging and a difficult problem. To this purpose, a parallel finite volume scheme on adaptive unstructured meshes for multi dimensional Saint-Venant system is presented. The adaptive mesh refinement method is based on a block-based decomposition (called BB-AMR) which allows quick meshing and easy parallelization. The main difficulty addressed here concerns the selection of the mesh refinement threshold which is certainly the most important parameter in the AMR method. Usually, the threshold is calibrated according to the test problem to balance the accuracy of the solution and the computational cost. To avoid " hand calibration " , we apply an automatic threshold method based on the decreasing rearrangement function of the mesh refinement criterion. This method is applied and validated successfully to the one and two dimensional non homogeneous Saint-Venant system through several tsunamis propagation test cases

Efficient surface water flow simulation on static Cartesian grid with local refinement according to key topographic features

Aiming at improving the computational efficiency without accuracy losses for surface water flow simulation, this paper presents a structured but non-uniform grid system incorporated into a Godunov-type finite volume scheme. The proposed grid system can detect the key topographic features in the computational domain where high-resolution mesh is in need for reliably solving the shallow water equations. The mesh refinement is automatically carried out in these areas while the mesh in the rest of the domain remains coarse. The criterion determining the refinement is suggested by a dimensionless number with a fixed value of 0.2 after sensitivity analysis. Three laboratory and field-scale test cases are employed to demonstrate the performance of the model for flow simulations on the new non-uniform grids. In all of the tests, the grid system is shown to successfully generate high-resolution mesh only in those areas with abruptly changing topographic features that dominate the flooding processes. To produce numerical solutions of similar accuracy, the non-uniform grid based model is able to accelerate by about two times comparing with the fine uniform grid based counterpart

A Numerical Method for Sharp-Interface Simulations of Multicomponent Alloy Solidification

We present a computational method for the simulation of the solidification of multicomponent alloys in the sharp-interface limit. Contrary to the case of binary alloys where a fixed point iteration is adequate, we hereby propose a Newton-type approach to solve the non-linear system of coupled PDEs arising from the time discretization of the governing equations, allowing for the first time sharp-interface simulations of the multialloy solidification. A combination of spatially adaptive quadtree grids, Level-Set Method, and sharp-interface numerical methods for imposing boundary conditions is used to accurately and efficiently resolve the complex behavior of the solidification front. The convergence behavior of the Newton-type iteration is theoretically analyzed in a one-dimensional setting and further investigated numerically in multiple spatial dimensions. We validate the overall computational method on the case of axisymmetric radial solidification admitting an analytical solution and show that the overall method's accuracy is close to second order. Finally, we perform numerical experiments for the directional solidification of a Co-Al-W ternary alloy with a phase diagram obtained from the PANDAT database and analyze the solutal segregation dependence on the processing conditions and alloy properties

A comparison of spectral element and finite difference methods using statically refined nonconforming grids for the MHD island coalescence instability problem

A recently developed spectral-element adaptive refinement incompressible magnetohydrodynamic (MHD) code [Rosenberg, Fournier, Fischer, Pouquet, J. Comp. Phys. 215, 59-80 (2006)] is applied to simulate the problem of MHD island coalescence instability (MICI) in two dimensions. MICI is a fundamental MHD process that can produce sharp current layers and subsequent reconnection and heating in a high-Lundquist number plasma such as the solar corona [Ng and Bhattacharjee, Phys. Plasmas, 5, 4028 (1998)]. Due to the formation of thin current layers, it is highly desirable to use adaptively or statically refined grids to resolve them, and to maintain accuracy at the same time. The output of the spectral-element static adaptive refinement simulations are compared with simulations using a finite difference method on the same refinement grids, and both methods are compared to pseudo-spectral simulations with uniform grids as baselines. It is shown that with the statically refined grids roughly scaling linearly with effective resolution, spectral element runs can maintain accuracy significantly higher than that of the finite difference runs, in some cases achieving close to full spectral accuracy.Comment: 19 pages, 17 figures, submitted to Astrophys. J. Supp

Adaptive modelling of long-distance wave propagation and fine-scale flooding during the Tohoku tsunami

The 11 March 2011 Tohoku tsunami is simulated using the quadtree-adaptive Saint-Venant solver implemented within the Gerris Flow Solver. The spatial resolution is adapted dynamically from 250 m in flooded areas up to 250 km for the areas at rest. Wave fronts are tracked at a resolution of 1.8 km in deep water. The simulation domain extends over 73° of both latitude and longitude and covers a significant part of the north-west Pacific. The initial wave elevation is obtained from a source model derived using seismic data only. Accurate long-distance wave prediction is demonstrated through comparison with DART buoys timeseries and GLOSS tide gauges records. The model also accurately predicts fine-scale flooding compared to both satellite and survey data. Adaptive mesh refinement leads to orders-of-magnitude gains in computational efficiency compared to non-adaptive methods. The study confirms that consistent source models for tsunami initiation can be obtained from seismic data only. However, while the observed extreme wave elevations are reproduced by the model, they are located further south than in the surveyed data. Comparisons with inshore wave buoys data indicate that this may be due to an incomplete understanding of the local wave generation mechanisms

Automatic parallel octree grid generation software with an extensible solver framework and a focus on urban simulation

The development of an automatic, dynamic, parallel, Cartesian, linear forest-of-octree grid generator and partial differential equation (PDE) solver framework is presented. This research is bundled into an application programmed with C++ which uses MPI for distributed parallelism. The application is named paros which stands for PARallel Octree Solver. In its current implementation, the application provides a \u27zeroth\u27 order representation of the target geometry, and as such, no cut-cell algorithm, projection method, or immersed boundary condition are implemented. In this case, \u27zeroth\u27 order means that no geometry is ever exactly represented in the final computational mesh: an octree element is either completely in the domain or entirely outside of it. Any element that contains or is intersected by a geometry facet is removed from the final mesh which results in a \u27blocky\u27 or \u27stepped\u27 geometry representation and simplifies boundary computations. The computational octree mesh creation is completely parallel and automated. The algorithm is dynamic in the sense that it is repartitioned dynamically throughout the grid generation process to maintain optimal load balancing during all phases of the mesh genesis. A linear octree data structure is used to store the octree mesh elements and is leveraged for optimal load balancing. An additional hierarchical octree is used to significantly improve algorithms that suffer from this linear storage paradigm. This work focuses on, but is not limited to, applications related to urban simulations and may be applied to plume/contaminant propagation. Within the PDE solution framework a cell-centered, incompressible, unsteady, Navier Stokes solver with an energy term to account for thermal buoyancy is implemented and validated using canonical test cases. Turbulence closure is implemented in the form of the Smagorinski large eddy simulation (LES) model. The parallel grid generation and solution process is tested on a large scale cityscape geometry and shown to be robust and efficient. Additionally, an implementation of the compressible Navier-Stokes equations is coded within the framework. The framework is extensible such that adding other types of numerical PDE solvers should not be difficult. Other features including adaptive mesh refinement (AMR) and contaminant transport functionality are included