Modeling 3D seismic problems using high-performance computing systems

This paper examines some issues in numerical modeling of seismology in three-dimensional space on high-performance computing systems. As a method of modeling, the grid-characteristic method is used. This method allows accurate staging of different contact conditions and is suitable for the most physically correct solutions of problems of seismology and seismic prospecting in complex heterogeneous media. We use the grid-characteristic schemes up to the 4th order accuracy inclusive. The software package is parallelized for work in a distributed clustered medium using the MPI technology. We present the results of the simulation of the Love and Rayleigh surface seismic waves, as well as the passage of seismic waves initiated by an earthquake’s hypocenter to the earth’s surface through a multilayer geological formation.     

An anisotropic unstructured triangular adaptive mesh algorithm based on error and error gradient information

In this paper, an algorithm based on unstructured triangular meshes using standard refinement patterns for anisotropic adaptive meshes is presented. It consists of three main actions: anisotropic refinement, solution-weighted smoothing and patch unrefinement. Moreover, a hierarchical mesh formulation is used. The main idea is to use the error and error gradient on each mesh element to locally control the anisotropy of the mesh. The proposed algorithm is tested on interpolation and boundary-value problems with a discontinuous solution.     

Mesh Quality: A Function of Geometry, Error Estimates or Both?

The issue of mesh quality for unstructured triangular and tetrahedral meshes is considered. The theoretical background to finite element methods is used to understand the basis of present-day geometrical mesh quality indicators. A survey of more recent research in the development of finite element methods describes work on anisotropic meshing algorithms, and on providing good error estimates that reveal the relationship between the error and both the mesh and the solution gradients. The reality of solving complex three-dimensional problems is that such indicators are presently not available for many problems of interest. A simple tetrahedral mesh quality measure using both geometrical and solution information is described. Some of the issues in mesh quality for unstructured tetrahedral meshes are illustrated by means of two simple examples.     

An automatic 3D mesh generation method for domains with multiple materials

This paper describes an automatic and efficient approach to construct unstructured tetrahedral and hexahedral meshes for a composite domain made up of heterogeneous materials. The boundaries of these material regions form non-manifold surfaces. In earlier papers, we developed an octree-based isocontouring method to construct unstructured 3D meshes for a single material (homogeneous) domain with manifold boundary. In this paper, we introduce the notion of a material change edge and use it to identify the interface between two or several different materials. A novel method to calculate the minimizer point for a cell shared by more than two materials is provided, which forms a non-manifold node on the boundary. We then mesh all the material regions simultaneously and automatically while conforming to their boundaries directly from volumetric data. Both material change edges and interior edges are analyzed to construct tetrahedral meshes, and interior grid points are analyzed for proper hexahedral mesh construction. Finally, edge-contraction and smoothing methods are used to improve the quality of tetrahedral meshes, and a combination of pillowing, geometric flow and optimization techniques is used for hexahedral mesh quality improvement. The shrink set of pillowing schemes is defined automatically as the boundary of each material region. Several application results of our multi-material mesh generation method are also provided.     

基于四面体剖分的并行地质块体建模方法

针对传统三维地质建模面临的内存消耗大,运行效率低的问题,提出了一种基于非结构四面体网格的并行地质块体建模方法。该方法采用“分治合并”的思想。首先将地质层位散点融合分割成若干个独立封闭块体;然后对每个块体进行表面三角形网格剖分,限定四面体剖分等步骤,形成块体的四面体网格剖分;最后将所有块体网格合并成最终地质模型。该方法中的块体剖分步骤应用多进程并行进行,提高了方法效率,并分摊计算机内存压力,可满足大尺度地质构造块体建模需求。该方法可为基于非结构网格的数值方法（如有限体积法、有限元方法等）的地震正演、偏移等算法提供合适的模型数据。     

Adaptive and Feature-Preserving Subdivision for High-Quality Tetrahedral Meshes

We present an adaptive subdivision scheme for unstructured tetrahedral meshes inspired by the       -subdivision scheme for triangular meshes. Existing tetrahedral subdivision schemes do not support adaptive refinement and have traditionally been driven by the need to generate smooth three-dimensional deformations of solids. These schemes use edge bisections to subdivide tetrahedra, which generates octahedra in addition to tetrahedra. To split octahedra into tetrahedra one routinely chooses a direction for the diagonals for the subdivision step. We propose a new topology-based refinement operator that generates only tetrahedra and supports adaptive refinement. Our tetrahedral subdivision algorithm is motivated by the need to have one representation for the modeling, the simulation and the visualization and so to bridge the gap between CAD and CAE. Our subdivision algorithm design emphasizes on geometric quality of the tetrahedral meshes, local and adaptive refinement operations, and preservation of sharp geometric features on the boundary and in the interior of the physical domain.     

Numerical study of dynamic processes in a continuous medium with a crack initiated by a near-surface disturbance by means of the grid-characteristic method

The purpose of this work is to study the problem of the near-surface disturbance propagation in a massive rock containing various heterogeneities, i.e., empty or filled cracks. Numerical solutions have been obtained for problems of wave propagation in such highly heterogeneous media, including those taking into account the plastic properties of the rock that can be manifested in the vicinity of a seismic gap or a well bore. All kinds of elastic and elastoplastic waves are analyzed resulting from the propagation of the initial disturbance and the waves arising from the reflection from the cracks and from the boundaries of the integration domain. An investigation was carried out of wave identification by means of seismograms obtained at the receiver located near the ground surface. In this study, the grid-characteristic method is employed using computational grids with triangular meshes and boundary conditions formulated at the interface between the rock and the crack, and on free surfaces in an explicit form. The proposed numerical method is extremely general and is suitable for investigations of the processes of seismic waves’ interaction with heterogeneous inclusions because it ensures the construction of the most correct computational algorithms at the boundaries of the integration domain and at the medium’s interface.     

Positivity-Preserving Discontinuous Galerkin Methods with Lax–Wendroff Time Discretizations

This work introduces a single-stage, single-step method for the compressible Euler equations that is provably positivity-preserving and can be applied on both Cartesian and unstructured meshes. This method is the first case of a single-stage, single-step method that is simultaneously high-order, positivity-preserving, and operates on unstructured meshes. Time-stepping is accomplished via the Lax–Wendroff approach, which is also sometimes called the Cauchy–Kovalevskaya procedure, where temporal derivatives in a Taylor series in time are exchanged for spatial derivatives. The Lax–Wendroff discontinuous Galerkin (LxW-DG) method developed in this work is formulated so that it looks like a forward Euler update but with a high-order time-extrapolated flux. In particular, the numerical flux used in this work is a convex combination of a low-order positivity-preserving contribution and a high-order component that can be damped to enforce positivity of the cell averages for the density and pressure for each time step. In addition to this flux limiter, a moment limiter is applied that forces positivity of the solution at finitely many quadrature points within each cell. The combination of the flux limiter and the moment limiter guarantees positivity of the cell averages from one time-step to the next. Finally, a simple shock capturing limiter that uses the same basic technology as the moment limiter is introduced in order to obtain non-oscillatory results. The resulting scheme can be extended to arbitrary order without increasing the size of the effective stencil. We present numerical results in one and two space dimensions that demonstrate the robustness of the proposed scheme.     

Numerical modeling of wave processes in layered media in the Arctic region

The aim of this work is the numerical simulation of wave propagation in media with linear-elastic and acoustic layers as exemplified by the seismic prospecting problems in the Arctic region and the explosive impact on an iceberg. The complete system of equations describing the state of a linearly elastic body and the system of equations describing the acoustic field are solved. The grid-characteristic method is used to provide the contact and boundary conditions, including the contact condition between acoustic and linear-elastic layers, to be correctly described.     

Array-based,parallel hierarchical mesh refinement algorithms for unstructured meshes

In this paper, we describe an array-based hierarchical mesh refinement capability through uniform refinement of unstructured meshes for efficient solution of PDE’s using finite element methods and multigrid solvers. A multi-degree, multi-dimensional and multi-level framework is designed to generate the nested hierarchies from an initial coarse mesh that can be used for a variety of purposes such as in multigrid solvers/preconditioners, to do solution convergence and verification studies and to improve overall parallel efficiency by decreasing I/O bandwidth requirements (by loading smaller meshes and in-memory refinement). We also describe a high-order boundary reconstruction capability that can be used to project the new points after refinement using high-order approximations instead of linear projection in order to minimize and provide more control on geometrical errors introduced by curved boundaries.The capability is developed under the parallel unstructured mesh framework “Mesh Oriented dAtaBase” (MOAB Tautges et al. (2004)). We describe the underlying data structures and algorithms to generate such hierarchies in parallel and present numerical results for computational efficiency and effect on mesh quality. We also present results to demonstrate the applicability of the developed capability to study convergence properties of different point projection schemes for various mesh hierarchies and to a multigrid finite-element solver for elliptic problems.     

Parallel unstructured tetrahedral mesh adaptation: algorithms,implementation and scalability

The use of unstructured adaptive tetrahedral meshes in the solution of transient flows poses a challenge for parallel computing due to the irregular and frequently changing nature of the data and its distribution. A parallel mesh adaptation algorithm, PTETRAD, for unstructured tetrahedral meshes (based on the serial code TETRAD) is described and analysed. The portable implementation of the parallel code in C with MPI is described and discussed. The scalability of the code is considered, analysed and illustrated by numerical experiments using a shock wave diffraction problem. Copyright © 1999 John Wiley & Sons, Ltd.     

An Algorithm for Three-Dimensional Mesh Generation for Arbitrary Regions with Cracks

An algorithm for generating unstructured tetrahedral meshes of arbitrarily shaped three-dimensional regions is described. The algorithm works for regions without cracks, as well as for regions with one or multiple cracks. The algorithm incorporates aspects of well known meshing procedures, but includes some original steps. It uses an advancing front technique, along with an octree to develop local guidelines for the size of generated elements. The advancing front technique is based on a standard procedure found in the literature, with two additional steps to ensure valid volume mesh generation for virtually any domain. The first additional step is related to the generation of elements only considering the topology of the current front, and the second additional step is a back-tracking procedure with face deletion, to ensure that a mesh can be generated even when problems happen during the advance of the front. To improve mesh quality (as far as element shape is concerned), an a posteriori local mesh improvement procedure is used. The performance of the algorithm is evaluated by application to a number of realistically complex, cracked geometries.     

Analysis and shape optimization in incompressible flows with unstructured grids

The aim of this study is to develop and validate numerical methods that perform shape optimization in incompressible flows using unstructured meshes. The three-dimensional Euler equations for compressible flow are modified using the idea of artificial compressibility and discretized on unstructured tetrahedral grids to provide estimates of pressure distributions for aerodynamic configurations. Convergence acceleration techniques like multigrid and residual averaging are used along with parallel computing platforms to enable these simulations to be performed in a few minutes. This computational frame-work is used to analyze sail geometries. The adjoint equations corresponding to the “incompressible” field equations are derived along with the functional form of gradients. The evaluation of the gradients is reduced to an integral around the boundary to circumvent hurdles posed by adjoint-based gradient evaluations on unstructured meshes. The reduced gradient evaluations provide major computational savings for unstructured grids and its accuracy and use for canonical and industrial problems is a major contribution of this study. The design process is driven by a steepest-descent algorithm with a fixed step-size. The feasibility of the design process is demonstrated for three inverse design problems, two canonical problems and one industrial problem.     

Uniformly high-order schemes on arbitrary unstructured meshes for advection–diffusion equations

The paper presents a linear high-order method for advection–diffusion conservation laws on three-dimensional mixed-element unstructured meshes. The key ingredient of the method is a reconstruction procedure in local computational coordinates. Numerical results illustrate the convergence rates for the linear equation and a non-linear hyperbolic system with diffusion terms for various types of meshes.     

基于边长的三维形状插值

形状插值在计算机图形学和几何处理中是一个极其重要而基础的问题,在计算机动画等领域有 着广泛应用。注意到在平面三角网格和三维四面体网格插值问题中,对边长平方插值等价于对回拉度量进行插 值,因此具有等距扭曲和共形扭曲同时有界的良好性质。通过将其推广至曲面三角网格,提出了一种完全基于 边长的曲面三角网格插值算法。给定边长,在重建网格阶段,使用牛顿法对边长误差能量进行优化。并且给出 了其海森矩阵的解析正定化形式,从而避免了高代价的特征值分解步骤。注意到四面体网格的边长平方插值结 果具有极低曲率,意味着只需少许修改即可将其压平从而嵌入三维空间。因此提出先将曲面三角网格四面体化, 再从四面体网格的插值结果提取表面。然后将这表面作为初始化用于边长误差能量的牛顿迭代,从而使得收敛 结果更加接近全局最优。在一系列三角网格上进行了实验,结果说明了本文方法比之前方法的边长误差更小, 且得到的结果还是有界扭曲的。     

Higher-order Finite Elements for Hybrid Meshes Using New Nodal Pyramidal Elements

We provide a comprehensive study of arbitrarily high-order finite elements defined on pyramids. We propose a new family of high-order nodal pyramidal finite element which can be used in hybrid meshes which include hexahedra, tetrahedra, wedges and pyramids. Finite elements matrices can be evaluated through approximate integration, and we show that the order of convergence of the method is conserved. Numerical results demonstrate the efficiency of hybrid meshes compared to pure tetrahedral meshes or hexahedral meshes obtained by splitting tetrahedra into hexahedra.     

支持裂纹扩展的三维有限元网格生成算法

描述了任意形状三维区域的非结构四面体网格生成算法,该算法对不含裂纹的区域、含单裂纹或多裂纹的区域都适用。算法首先使用八叉树来确定网格单元大小,然后采用前沿推进技术来生成网格。在前沿推进过程中,采用基于几何形状和基于拓扑结构的两个步骤来保证前沿向前移动过程中发生问题时仍能进行正确执行,并且使用了一种局部网格优化方法来提高网格划分的质量。最后,将算法运用到带有裂纹的复杂实体模型,实验结果表明该算法具有较强的适用性和较高的性能。     

A high efficiency parallel unstructured solver dedicated to internal combustion engine simulation

IFP-C3D, a hexahedral unstructured parallel solver dedicated to multiphysics calculation is being developed at IFP to compute the compressible combustion in internal engines. IFP-C3D uses an unstructured formalism, the finite volume method on staggered grid, time splitting, SIMPLE loop, subcycled advection, turbulent and Lagrangian spray and a liquid film model. Original algorithms and models such as the conditional temporal interpolation methodology for moving grids, the remapping algorithm for transferring quantities on different meshes during the computation enable IFP-C3D to deal with complex moving geometries with large volume deformation induced by all moving geometrical parts (intake/exhaust valve, piston). Large super-scalar machines up to 1000 processors are being widely used and IFP-C3D has been optimized for running on these Cluster machines. IFP-C3D is parallelized using the Message Passing Interface (MPI) library to distribute a calculation over a large number of processors. Moreover, IFP-C3D uses an optimized linear algebraic library to solve linear matrix systems and the METIS partitionner library to distribute the computational load equally for all meshes used during the calculation and in particular during the remap stage when new meshes are loaded. Numerical results and performance timing are presented to demonstrate the computational efficiency of the code.     

Adaptive space–time finite element methods for the wave equation on unbounded domains

Comprehensive adaptive procedures with efficient solution algorithms for the time-discontinuous Galerkin space–time finite element method (DGFEM) including high-order accurate nonreflecting boundary conditions (NRBC) for unbounded wave problems are developed. Sparse multi-level iterative schemes based on the Gauss–Seidel method are developed to solve the resulting fully-discrete system equations for the interior hyperbolic equations coupled with the first-order temporal equations associated with auxiliary functions in the NRBC. Due to the local nature of wave propagation, the iterative strategy requires only a few iterations per time step to resolve the solution to high accuracy. Further cost savings are obtained by diagonalizing the mass and boundary damping matrices. In this case the algebraic structure decouples the diagonal block matrices giving rise to an explicit multi-corrector method. An h-adaptive space–time strategy is employed based on the Zienkiewicz–Zhu spatial error estimate using the superconvergent patch recovery (SPR) technique, together with a temporal error estimate arising from the discontinuous jump between time steps of both the interior field solutions and auxiliary boundary functions. For accurate data transfer between meshes, a new enhanced interpolation (EI) method is developed and compared to standard interpolation and projection. Numerical studies of transient radiation and scattering demonstrate the accuracy, reliability and efficiency gained from the adaptive strategy.     

Space-time Generalized Riemann Problem Solvers of Order k for Linear Advection with Unrestricted Time Step

This work concerns high-order approximations of the linear advection equation in very long time. A GRP-type scheme of arbitrary high-order in space and time with no restriction on the time step is developed. In the usual GRP solvers, we consider a polynomial approximation of the solution in space in each cell at the initial time. Here, we add a second polynomial approximation of the solution in time in each interface. Thanks to this double approximation, the resulting scheme is compact. It is proved to be of order k+1 in space and time, where k is the degree of the polynomials. Thanks to the compactness of the scheme, a two-dimensional extension is detailed on unstructured meshes made of triangles. Several numerical test-cases and comparison with existing methods illustrate the excellent behaviour of the scheme.