skip to main content
OSTI.GOV title logo U.S. Department of Energy
Office of Scientific and Technical Information

Title: Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling

Abstract

The mass-lumped method avoids the cost of inverting the mass matrix and simultaneously maintains spatial accuracy by adopting additional interior integration points, known as cubature points. To date, such points are only known analytically in tensor domains, such as quadrilateral or hexahedral elements. Thus, the diagonal-mass-matrix spectral element method (SEM) in non-tensor domains always relies on numerically computed interpolation points or quadrature points. However, only the cubature points for degrees 1 to 6 are known, which is the reason that we have developed a p-norm-based optimization algorithm to obtain higher-order cubature points. In this way, we obtain and tabulate new cubature points with all positive integration weights for degrees 7 to 9. The dispersion analysis illustrates that the dispersion relation determined from the new optimized cubature points is comparable to that of the mass and stiffness matrices obtained by exact integration. Simultaneously, the Lebesgue constant for the new optimized cubature points indicates its surprisingly good interpolation properties. As a result, such points provide both good interpolation properties and integration accuracy. The Courant–Friedrichs–Lewy (CFL) numbers are tabulated for the conventional Fekete-based triangular spectral element (TSEM), the TSEM with exact integration, and the optimized cubature-based TSEM (OTSEM). A complementary study demonstrates themore » spectral convergence of the OTSEM. A numerical example conducted on a half-space model demonstrates that the OTSEM improves the accuracy by approximately one order of magnitude compared to the conventional Fekete-based TSEM. In particular, the accuracy of the 7th-order OTSEM is even higher than that of the 14th-order Fekete-based TSEM. Furthermore, the OTSEM produces a result that can compete in accuracy with the quadrilateral SEM (QSEM). The high accuracy of the OTSEM is also tested with a non-flat topography model. In terms of computational efficiency, the OTSEM is more efficient than the Fekete-based TSEM, although it is slightly costlier than the QSEM when a comparable numerical accuracy is required. - Highlights: • Higher-order cubature points for degrees 7 to 9 are developed. • The effects of quadrature rule on the mass and stiffness matrices has been conducted. • The cubature points have always positive integration weights. • Freeing from the inversion of a wide bandwidth mass matrix. • The accuracy of the TSEM has been improved in about one order of magnitude.« less

Authors:
 [1];  [1];  [1];  [2];  [3]
  1. State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, 100029 (China)
  2. (China)
  3. Physics of the Earth, Sciences B, University of Zaragoza, Pedro Cerbuna 12, 50009 Zaragoza (Spain)
Publication Date:
OSTI Identifier:
22622290
Resource Type:
Journal Article
Resource Relation:
Journal Name: Journal of Computational Physics; Journal Volume: 336; Other Information: Copyright (c) 2017 Elsevier Science B.V., Amsterdam, The Netherlands, All rights reserved.; Country of input: International Atomic Energy Agency (IAEA)
Country of Publication:
United States
Language:
English
Subject:
71 CLASSICAL AND QUANTUM MECHANICS, GENERAL PHYSICS; 97 MATHEMATICAL METHODS AND COMPUTING; ACCURACY; ALGORITHMS; APPROXIMATIONS; COMPARATIVE EVALUATIONS; COMPUTERIZED SIMULATION; CONVERGENCE; DISPERSION RELATIONS; DISPERSIONS; EFFICIENCY; INTERPOLATION; MASS; OPTIMIZATION; QUADRATURES; SCANNING ELECTRON MICROSCOPY; TOPOGRAPHY

Citation Formats

Liu, Youshan, E-mail: ysliu@mail.iggcas.ac.cn, Teng, Jiwen, E-mail: jwteng@mail.iggcas.ac.cn, Xu, Tao, E-mail: xutao@mail.iggcas.ac.cn, CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing, 100101, and Badal, José, E-mail: badal@unizar.es. Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling. United States: N. p., 2017. Web. doi:10.1016/J.JCP.2017.01.069.
Liu, Youshan, E-mail: ysliu@mail.iggcas.ac.cn, Teng, Jiwen, E-mail: jwteng@mail.iggcas.ac.cn, Xu, Tao, E-mail: xutao@mail.iggcas.ac.cn, CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing, 100101, & Badal, José, E-mail: badal@unizar.es. Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling. United States. doi:10.1016/J.JCP.2017.01.069.
Liu, Youshan, E-mail: ysliu@mail.iggcas.ac.cn, Teng, Jiwen, E-mail: jwteng@mail.iggcas.ac.cn, Xu, Tao, E-mail: xutao@mail.iggcas.ac.cn, CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing, 100101, and Badal, José, E-mail: badal@unizar.es. Mon . "Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling". United States. doi:10.1016/J.JCP.2017.01.069.
@article{osti_22622290,
title = {Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling},
author = {Liu, Youshan, E-mail: ysliu@mail.iggcas.ac.cn and Teng, Jiwen, E-mail: jwteng@mail.iggcas.ac.cn and Xu, Tao, E-mail: xutao@mail.iggcas.ac.cn and CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing, 100101 and Badal, José, E-mail: badal@unizar.es},
abstractNote = {The mass-lumped method avoids the cost of inverting the mass matrix and simultaneously maintains spatial accuracy by adopting additional interior integration points, known as cubature points. To date, such points are only known analytically in tensor domains, such as quadrilateral or hexahedral elements. Thus, the diagonal-mass-matrix spectral element method (SEM) in non-tensor domains always relies on numerically computed interpolation points or quadrature points. However, only the cubature points for degrees 1 to 6 are known, which is the reason that we have developed a p-norm-based optimization algorithm to obtain higher-order cubature points. In this way, we obtain and tabulate new cubature points with all positive integration weights for degrees 7 to 9. The dispersion analysis illustrates that the dispersion relation determined from the new optimized cubature points is comparable to that of the mass and stiffness matrices obtained by exact integration. Simultaneously, the Lebesgue constant for the new optimized cubature points indicates its surprisingly good interpolation properties. As a result, such points provide both good interpolation properties and integration accuracy. The Courant–Friedrichs–Lewy (CFL) numbers are tabulated for the conventional Fekete-based triangular spectral element (TSEM), the TSEM with exact integration, and the optimized cubature-based TSEM (OTSEM). A complementary study demonstrates the spectral convergence of the OTSEM. A numerical example conducted on a half-space model demonstrates that the OTSEM improves the accuracy by approximately one order of magnitude compared to the conventional Fekete-based TSEM. In particular, the accuracy of the 7th-order OTSEM is even higher than that of the 14th-order Fekete-based TSEM. Furthermore, the OTSEM produces a result that can compete in accuracy with the quadrilateral SEM (QSEM). The high accuracy of the OTSEM is also tested with a non-flat topography model. In terms of computational efficiency, the OTSEM is more efficient than the Fekete-based TSEM, although it is slightly costlier than the QSEM when a comparable numerical accuracy is required. - Highlights: • Higher-order cubature points for degrees 7 to 9 are developed. • The effects of quadrature rule on the mass and stiffness matrices has been conducted. • The cubature points have always positive integration weights. • Freeing from the inversion of a wide bandwidth mass matrix. • The accuracy of the TSEM has been improved in about one order of magnitude.},
doi = {10.1016/J.JCP.2017.01.069},
journal = {Journal of Computational Physics},
number = ,
volume = 336,
place = {United States},
year = {Mon May 01 00:00:00 EDT 2017},
month = {Mon May 01 00:00:00 EDT 2017}
}
  • We implement a high-order finite-element application, which performs the numerical simulation of seismic wave propagation resulting for instance from earthquakes at the scale of a continent or from active seismic acquisition experiments in the oil industry, on a large cluster of NVIDIA Tesla graphics cards using the CUDA programming environment and non-blocking message passing based on MPI. Contrary to many finite-element implementations, ours is implemented successfully in single precision, maximizing the performance of current generation GPUs. We discuss the implementation and optimization of the code and compare it to an existing very optimized implementation in C language and MPI onmore » a classical cluster of CPU nodes. We use mesh coloring to efficiently handle summation operations over degrees of freedom on an unstructured mesh, and non-blocking MPI messages in order to overlap the communications across the network and the data transfer to and from the device via PCIe with calculations on the GPU. We perform a number of numerical tests to validate the single-precision CUDA and MPI implementation and assess its accuracy. We then analyze performance measurements and depending on how the problem is mapped to the reference CPU cluster, we obtain a speedup of 20x or 12x.« less
  • The scalability of computational applications on current and next-generation supercomputers is increasingly limited by the cost of inter-process communication. We implement non-blocking asynchronous communication in the High-Order Methods Modeling Environment for the time integration of the hydrostatic fluid equations using both the spectral-element and discontinuous Galerkin methods. This allows the overlap of computation with communication, effectively hiding some of the costs of communication. A novel detail about our approach is that it provides some data movement to be performed during the asynchronous communication even in the absence of other computations. This method produces significant performance and scalability gains in large-scalemore » simulations.« less
  • This paper presents the development of parallel direct Vlasov solvers using the Spectral Element Method (SEM). Instead of the standard Particle-In-Cell (PIC) approach for kinetic space plasma simulation, i.e. solving the Vlasov-Maxwell equations, the direct method has been used in this paper. There are several benefits to solve the Vlasov equation directly, such as avoiding noise associated with the finite number of particles and the capability to capture the fine structure in the plasma, etc. The most challenging part of direct Vlasov solver comes from high dimension, as the computational cost increases as N{sup 2d}, where d is the dimensionmore » of the physical space. Recently due to fast development of supercomputers, the possibility of high dimensions becomes more realistic. A significant effort has been devoted to solve the Vlasov equation in low dimensions so far, now more interests focus on higher dimensions. Different numerical methods have been tried so far, such as finite difference method, Fourier spectral method, finite volume method, etc. In this paper SEM has been successfully applied to construct these solvers. SEM has shown several advantages, such as easy interpolation due to local element structure and long time integration due to its high order accuracy. Domain decomposition in high dimensions have been used for parallelization, these include scalable parallel 1D and 2D Poisson solvers. Benchmark results have been shown and simulation results have been reported for two different cases: one dimension (1P1V), and two dimensions (2P2V) in both physical and velocity spaces.« less
  • In the variational multiscale (VMS) approach to large eddy simulation (LES), the governing equations are projected onto an a priori scale partitioning of the solution space. This gives an alternative framework for designing and analyzing turbulence models. We describe the implementation of the VMS LES methodology in a high order spectral element method with a nodal basis, and discuss the properties of the proposed scale partitioning. The spectral element code is first validated by doing a direct numerical simulation of fully developed plane channel flow. The performance of the turbulence model is then assessed by several coarse grid simulations ofmore » channel flow at different Reynolds numbers.« less
  • The boundary element-image method is proposed as a method for generating synthetic seismograms in a system of piecewise homogeneous layers. In contrast with many other methods (such as finite difference) in which the boundary conditions are represented simply by a change of material properties, the present approach uses the boundary conditions explicitly in order to develop a system of integral equations. As in the indirect boundary element method, the scattered seismic waves are presumed to result from the action of fictitious sources. However, as in the image method, these fictitious sources are located outside of the domain through which themore » scattered waves propagate. This combination of the image and boundary element methods thereby avoids introducing additional singularities (in the form of fictitious sources) into the wave equation for the scattered waves. Although this method can be formulated in a general way, specific computational advantages occur when the basis and testing functions are chosen as harmonically related complex sinusoids. Physically, this corresponds to a plane wave expansion of the scattered waves. Computationally, when considering piecewise linear interfaces, this results in a (possibly) large system of simultaneous algebraic equations having matrix elements which can be calculated analytically. This method is applied to the reflection of seismic waves from a V-shaped reflector subjected to two different incident waves: a plane wave and a cylindrical wave created by a line source. The synthetic seismograms for both cases are presented (in the form of time sections) and then processed in an attempt to reconstruct the geometry of the reflector (in the form of a depth section).« less