Spectral Element Discretization Research Articles

Implementation of the spectral element discretization of the nonstationary velocity–vorticity–pressure Stokes formulation

This paper deals with the resolution of the time-dependent velocity–vorticity–pressure formulation of Stokes problem in a two- and three-dimensional domain. The discretization is based on the spectral element method with respect to the space variable and Euler’s implicit scheme with respect to the time variable. We implement the asymmetric linear system to solve the full spectral element discrete problem. Detailed numerical experiments lead to confirm the usefulness of the discretization.

The algorithmic resolution of spectral-element discretization for the time-dependent Stokes problem

We consider two algorithms for the resolution of the time-dependent Stokes problem with nonstandard boundary conditions by the domain-decomposition spectral-element method. The first algorithm (Elimination method) is based on the Uzawa method by decoupling the linear system, while the second algorithm (Global inversion) is based on the overall resolution of the system by the GMRES method. A detailed implementation is proposed and some numerical tests are carried out in two and three dimensions and where the domain is multiply connected.

A spectral element discretization for quasi‐static magnetohydrodynamic flows

AbstractThe classical staggered spectral element method (SEM) is revisited and extended to quasi‐static magnetohydrodynamic (MHD) flows. In this realm, which is valid in the limit of vanishing magnetic Reynolds number, the evaluation of the Lorentz force in the momentum equation requires the electric current density, governed by Ohm's law and a charge conservation condition derived from Ampère's law, to be determined. Once discretized with the SEM, this translates into solving one additional problem for the electric potential involving the so‐called consistent Poisson operator. The method is well suited for fully three‐dimensional flows in complex geometries. Changes in resolution requirements aside, consideration of the electromagnetic quantities is estimated to increase the computational cost associated with MHD by about 40% relative to hydrodynamics. The accuracy and the capabilities of the scheme is demonstrated on a set of common flows from the MHD literature. Exponential convergence with polynomial order is confirmed for the electric current density.

The posteriori analysis of the spectral element discretization of the wave equation

This work focuses on discretizing a second-order linear wave equation using the implicit Euler scheme for time discretization and the spectral element method for spatial discretization. We prove that optimal adaptivity can be achieved by combining adaptive time steps with adaptive spectral mesh. We introduce two sets of error indicators for time and space, respectively, and derive optimal estimates.

A fully scalable homogenization method to upscale 3-D elastic media

SUMMARY Modelling seismic wavefields in complex 3-D elastic media is the key in many fields of Earth Science: seismology, seismic imaging, seismic hazard assessment and earthquake source mechanism reconstruction. This modelling operation can incur significant computational cost, and its accuracy depends on the ability to take into account the scales of the subsurface heterogeneities varying. The theory of homogenization describes how the small-scale heterogeneities interact with the seismic waves and allows to upscale elastic media consistently with the wave equation. In this study, an efficient and scalable numerical homogenization tool is developed, relying on the similarity between the equations describing the propagation of elastic waves and the homogenization process. By exploiting the optimized implementation of an elastic modelling kernel based on a spectral-element discretization and domain decomposition, a fully scalable homogenization process, working directly on the spectral-element mesh, is presented. Numerical experiments on the entire SEAM II foothill model and a 3-D version of the Marmousi II model illustrate the efficiency and flexibility of this approach. A reduction of two orders of magnitude in terms of absolute computational cost is observed on the elastic wave modelling of the entire SEAM II model at a controlled accuracy.

Neko: A modern, portable, and scalable framework for high-fidelity computational fluid dynamics

Computational fluid dynamics (CFD), in particular applied to turbulent flows, is a research area with great engineering and fundamental physical interest. However, already at moderately high Reynolds numbers the computational cost becomes prohibitive as the range of active spatial and temporal scales is quickly widening. Specifically scale-resolving simulations, including large-eddy simulation (LES) and direct numerical simulations (DNS), thus need to rely on modern efficient numerical methods and corresponding software implementations. Recent trends and advancements, including more diverse and heterogeneous hardware in High-Performance Computing (HPC), are challenging software developers in their pursuit for good performance and numerical stability. The well-known maxim “software outlives hardware” may no longer necessarily hold true, and developers are today forced to re-factor their codebases to leverage these powerful new systems. In this paper, we present Neko, a new portable framework for high-order spectral element discretization, targeting turbulent flows in moderately complex geometries. Neko is fully available as open software. Unlike prior works, Neko adopts a modern object-oriented approach in Fortran 2008, allowing multi-tier abstractions of the solver stack and facilitating hardware backends ranging from general-purpose processors (CPUs) down to exotic vector processors and FPGAs. We show that Neko’s performance and accuracy are comparable to NekRS, and thus on-par with Nek5000’s successor on modern CPU machines. Furthermore, we develop a performance model, which we use to discuss challenges and opportunities for high-order solvers on emerging hardware.

Spectral element discretization of the time-dependent Stokes problem with nonstandard boundary conditions

This work deals with the spectral element discretization of the time-dependent Stokes problem in two- and three-dimensional domains. The boundary condition is defined on the normal component of the velocity and the tangential components of the vorticity. The discretization related to the time variable is processed by a Backward Euler method. We prove through a detailed numerical analysis the well-posedness of the full discrete problem.

The error spectral element discretization analysis of the non stationary velocity vorticity pressure Stokes formulation

The error spectral element discretization analysis of the non stationary velocity vorticity pressure Stokes formulation

Three-Dimensional Forward Modeling of Transient Electromagnetic Method Considering Induced Polarization Effect Based on Spectral Element Method

The transient electromagnetic method (TEM) is widely used in the exploration of mineral, petroleum, and geothermal resources due to its sensitivity to low-resistivity bodies, limited site constraints, and strong resistance to interference. In practical applications, the TEM often uses a long wire source instead of an idealized horizontal electric dipole (HED) source as the excitation source. This is due to the complex external conditions and the relatively large distance between the receiving zone and the transmitter source. Compared to the HED, the long wire source can provide a larger excitation current, generating stronger signals to meet the requirements of a higher signal-to-noise ratio or deeper exploration. It also produces longer-duration signals, thereby providing better resolution. Additionally, for the interpretation of TEM data, three-dimensional forward modeling plays a crucial role. However, the mature traditional TEM forward method is based on a simple, sometimes inappropriate model, as it is well established that the induced polarization (IP) effect is widely present in the deep earth, especially in oil and gas reservoirs. The presence of the IP effect results in negative responses in field data that do not conform to the traditional theoretical decay law of TEM, which can significantly impact data processing and inversion results. To address this issue, a TEM forward modeling method considering the IP effect based on the spectral element method (SEM) has been developed in this study. Firstly, starting from the Helmholtz equation satisfied by the time domain electric field, we introduce the Debye model with polarization information into the forward modeling by utilizing the differential form of Ohm’s law. As a result, we derive the boundary value problem for the time domain electric field that considers the induced polarization effect. Using Gauss–Lobatto–Legendre (GLL) polynomials as the basis functions, the SEM is employed to discretize the governing equations at each time step and obtain spectral element discretization equations. Then, temporal discretization equations are derived using the second-order backward Euler formula, and the linear system of equations is solved using the Pardiso direct solver. Finally, the electromagnetic responses at any time channel are obtained via SEM interpolation and numerical integration, thereby achieving three-dimensional TEM forward modeling considering the IP effect. The results indicate that this method can effectively reflect the spatial distribution of polarizable subsurface media. It provides valuable references for studying the polarization parameters of subsurface media and performing a three-dimensional inversion of TEM data considering the induced polarization effect.

Symplectic time-domain finite element method (STD-FEM) for room acoustic modeling

A new, extended symplectic time-domain finite element method (STD-FEM) is proposed for room acoustic modeling. As a mathematical model for sound propagation and reflection, the classical time-domain wave equations are used, which can be extended with air absorption. Frequency-dependent locally reactive boundary conditions are also introduced to the model. Spatially, a third-order tensor product type spectral element discretization is applied, which allows us to use explicit time steps. For this purpose, a partitioned Runge-Kutta method is employed, which is an extension of a third-order symplectic time-discretization. The method is validated on a number of examples and the performance and stability benefits are presented. We present the implementation of the method in the soundy.ai application.

Geofluid Object Workbench (GeoFLOW) for Atmospheric Dynamics in the Approach to Exascale: Spectral Element Formulation and CPU Performance

Abstract A new software framework using a well-established high-order spectral element discretization is presented for solving the compressible Navier–Stokes equations for purposes of research in atmospheric dynamics in bounded and unbounded limited-area domains, with a view toward capturing spatiotemporal intermittency that may be particularly challenging to attain using low-order schemes. A review of the discretization is provided, emphasizing properties such as the matrix product formalism and other design considerations that will facilitate its effective use on emerging exascale platforms, and a new geometry-independent, element boundary exchange method is described to maintain continuity. A variety of test problems are presented that demonstrate accuracy of the implementation primarily in wave-dominated or transitional flow regimes; conservation properties are also demonstrated. A strong scaling CPU study in a three-dimensional domain without using threading shows an average parallel efficiency of ≳99% up to 2 × 104 MPI tasks that is not affected negatively by expansion polynomial order. On-node performance is also examined and reveals that, while the primary numerical operations achieve their theoretical arithmetic intensity, the application performance is largely limited by available memory bandwidth. Significance Statement This work considers the need for computationally efficient, high-order, low dissipation numerics to fully leverage emerging exascale computing resources in an effort to examine and improve the accuracy of numerical treatments of atmospheric and weather phenomena. A new spectral element implementation is introduced that attempts to address the issues involved. Well-understood tests are presented that illustrate the known efficacy of the method in wave-dominated, quasi-laminar, and relatively strong shear flow regimes, and good conservation properties for mass and total energy are achieved. Importantly, the implementation is shown to exhibit encouraging performance characteristics.

A fully-implicit hybridized discontinuous Galerkin spectral element method for two phase flow in petroleum reservoirs

A hybridized discontinuous Galerkin method with high order spectral element discretization is presented for modeling multiphase immiscible 2D fluid flow in petroleum reservoir simulations. The time dependent nonlinear governing equations for reservoir simulation are presented, and a fully coupled implicit time stepping procedure is implemented using a high order W-Rosenbrock algorithm. For incompressible fluids considered in the present study, the governing PDE's include an elliptic PDE for the pressure field and a parabolic convective diffusion equation for the saturation (water volume fraction) in a two phase oil-water system. The PDE systems are coupled through the highly nonlinear relative permeability and capillary pressure functions relating the fluid velocities of individual species to the individual species pressures and saturations. With HDG/spectral element (HDGSEM) spatial discretization of the coupled governing equations, the system constitutes an index 1 differential algebraic equation. It is shown that the W-Rosenbrock algorithm in combination with the HDGSEM method yields an efficient and effective solution procedure which demonstrates robust, high order temporal accuracy.Numerous physical tests cases are presented demonstrating that the versatility of this approach is sufficient to address the challenges in modeling heterogeneous permeability fields and geologic features of arbitrary shape with discontinuous permeability functions. It is also capable of resolving the sharp time dependent displacement fronts in oil-water systems encountered in secondary oil recovery. It is demonstrated that the high order polynomial basis functions (up to 12th order) used in the HDGSEM algorithm have the ability to resolve complex geometry and achieve accurate results utilizing a simple Cartesian mesh for non-conforming geological features. Eliminating the need to re-grid and conform the mesh to the boundaries of different geological features greatly simplifies the workflow for petroleum engineers.

Iterative solution methods for 3D controlled-source electromagnetic forward modelling of geophysical exploration scenarios

We develop an efficient and robust iterative framework suitable for solving the linear system of equations resulting from the spectral element discretisation of the curl-curl equation of the total electric field encountered in geophysical controlled-source electromagnetic applications. We use the real-valued equivalent form of the original complex-valued system and solve this arising real-valued two-by-two block system (outer system) using the generalised conjugate residual method preconditioned with a highly efficient block-based PREconditioner for Square Blocks (PRESB). Applying this preconditioner equates to solving two smaller inner symmetric systems which are either solved using a direct solver or iterative methods, namely the generalised conjugate residual or the flexible generalised minimal residual methods preconditioned with the multigrid-based auxiliary-space preconditioner AMS. Our numerical experiments demonstrate the robustness of the outer solver with respect to spatially variable material parameters, for a wide frequency range of five orders of magnitude (0.1-10’000 Hz), with respect to the number of degrees of freedom, and for stretched structured and unstructured as well as locally refined meshes. For all the models considered, the outer solver reaches convergence in a small (typically < 20) number of iterations. Further, our numerical tests clearly show that solving the two inner systems iteratively using the indicated preconditioned iterative methods is computationally beneficial in terms of memory requirement and time spent as compared to a direct solver. On top of that, our iterative framework works for large-scale problems where direct solvers applied to the original complex-valued systems succumb due to their excessive memory consumption, thus making the iterative framework better suited for large-scale 3D problems. Comparison to a similar iterative framework based on a block-diagonal and the auxiliary-space preconditioners reveals that the PRESB preconditioner requires slightly fewer iterations to converge yielding a certain gain in time spent to obtain the solution of the two-by-two block system.

A priori and a posteriori error analysis of hp spectral element discretization for optimal control problems with elliptic equations

. A priori and a posteriori error estimates of hp spectral element method for an elliptic optimal control problem with L2-norm control constraint are investigated in this paper. We then present the optimality conditions and set up hp spectral element discretization scheme. Based on some important interpolation operators and suitable immediate variable, a priori error estimates for this problem in hp spectral element discretization are established carefully. Furthermore, using some interpolation operators, a rigorous posteriori error estimates of hp spectral element approximation are also proved for control and the coupled state approximation in L2−H1-norm and L2−L2-norm, respectively. Such estimators can be used to construct reliable adaptive spectral element methods for optimal control problems. Finally, the error analysis results are confirmed by numerical results.

Transient numerical approximation of hyperbolic diffusions and beyond

In this paper different types of hyperbolic diffusions and their corresponding transient Fokker–Planck equation are described and numerical solutions are presented. Diffusion models were developed that can fit both the marginal distribution and correlation structure and they have found a wide application in finance, turbulence and environmental time series. Hyperbolic diffusions have a complicated structure and variety of parameters and are extremely difficult to study and to model. We propose a numerical technique that solves the one-dimensional hyperbolic Fokker–Planck equation in the time dependent case. Note that this is the first study where transient hyperbolic diffusions are considered. The numerical technique is based on an adaptive reduced basis method using a spectral element discretization. It involves enrichment and projection stages where an optimal basis is found in a dynamic way using the singular value decomposition (SVD). The approach dramatically reduces the number of degrees of freedom required to solve the problem. The numerical simulations of the Fokker–Planck equation are verified with available stationary solutions.

Metrics for Intercomparison of Remapping Algorithms (MIRA) protocol applied to Earth system models

Abstract. Strongly coupled nonlinear phenomena such as those described by Earth system models (ESMs) are composed of multiple component models with independent mesh topologies and scalable numerical solvers. A common operation in ESMs is to remap or interpolate component solution fields defined on their computational mesh to another mesh with a different combinatorial structure and decomposition, e.g., from the atmosphere to the ocean, during the temporal integration of the coupled system. Several remapping schemes are currently in use or available for ESMs. However, a unified approach to compare the properties of these different schemes has not been attempted previously. We present a rigorous methodology for the evaluation and intercomparison of remapping methods through an independently implemented suite of metrics that measure the ability of a method to adhere to constraints such as grid independence, monotonicity, global conservation, and local extrema or feature preservation. A comprehensive set of numerical evaluations is conducted based on a progression of scalar fields from idealized and smooth to more general climate data with strong discontinuities and strict bounds. We examine four remapping algorithms with distinct design approaches, namely ESMF Regrid (Hill et al., 2004), TempestRemap (Ullrich and Taylor, 2015), generalized moving least squares (GMLS) (Trask and Kuberry, 2020) with post-processing filters, and WLS-ENOR (Li et al., 2020). By repeated iterative application of the high-order remapping methods to the test fields, we verify the accuracy of each scheme in terms of their observed convergence order for smooth data and determine the bounded error propagation using challenging, realistic field data on both uniform and regionally refined mesh cases. In addition to retaining high-order accuracy under idealized conditions, the methods also demonstrate robust remapping performance when dealing with non-smooth data. There is a failure to maintain monotonicity in the traditional L2-minimization approaches used in ESMF and TempestRemap, in contrast to stable recovery through nonlinear filters used in both meshless GMLS and hybrid mesh-based WLS-ENOR schemes. Local feature preservation analysis indicates that high-order methods perform better than low-order dissipative schemes for all test cases. The behavior of these remappers remains consistent when applied on regionally refined meshes, indicating mesh-invariant implementations. The MIRA intercomparison protocol proposed in this paper and the detailed comparison of the four algorithms demonstrate that the new schemes, namely GMLS and WLS-ENOR, are competitive compared to standard conservative minimization methods requiring computation of mesh intersections. The work presented in this paper provides a foundation that can be extended to include complex field definitions, realistic mesh topologies, and spectral element discretizations, thereby allowing for a more complete analysis of production-ready remapping packages.

3-D multicomponent full waveform inversion for shallow-seismic target: Ettlingen Line case study

SUMMARY In this work, we investigate methodological development and application of viscoelastic full waveform inversion to a multicomponent source and receiver data set at the near-surface scale. The target is located in Germany and lies in an ancient war trench from the War of the Spanish Succession that was buried at approximately 1 m depth. We present the pre-processing of the data, including a matching filter correction to compensate for different source and receiver coupling conditions during the acquisition, as well as a dedicated multistep workflow for the reconstruction of both P-wave and S-wave velocities. Our implementation is based on viscoelastic modelling using a spectral element discretization to accurately account for the wave propagation complexity in this shallow region. We illustrate the inversion stability by starting from different initial models, either based on dispersion curve analysis or homogeneous models consistent with first arrivals. We recover similar results in both cases. We also illustrate the importance of accounting for the attenuation by comparing elastic and viscoelastic results. The 3-D results make it possible to recover and locate the trench precisely in terms of interpretation. They also exhibit another trench structure, in a direction forming an angle at 45° with the direction of the targeted trench. This new structure had been previously interpreted as an artifact in former 2-D inversion results. The archaeological interpretation of this new structure is still a matter of discussion.

A spectral element approach to computing normal modes

SUMMARYWe introduce a new approach to the computation of gravito-elastic free oscillations or normal modes of spherically symmetric bodies based on a spectral element discretization of the radial ordinary differential equations. Our method avoids numerical instabilities often encountered in the classical method of radial integration and root finding of the characteristic function. To this end, the code is built around a sparse matrix formulation of the eigenvalue problem taking advantage of state-of-the-art parallel iterative solvers. We apply the method to toroidal, spheroidal and radial modes and we demonstrate its versatility in the presence of attenuation, fluid layers and gravity (including the purely elastic case, the Cowling approximation, and full gravity). We demonstrate higher-order convergence and verify the software by computing seismograms and comparing these to existing numerical solutions. Finally, to emphasize the general applicability of our code, we show spectra and eigenfunctions of Earth, Mars and Jupiter’s icy moon Europa and discuss the different types of modes that emerge.

Isogeometric collocation discretizations for acoustic wave problems

The acoustic wave problem is here discretized by collocation isogeometric analysis (IGA) in space and Newmark schemes of first and second order in time. A numerical study in the plane on both Cartesian and NURBS domains investigates the convergence rate of the proposed collocation Newmark-IGA method and its dependence on the main isogeometric parameters, the mesh size h, the spline polynomial degree p, the spline regularity k, and on the time step size Δt. In addition, a Ricker wavelet propagation test is accurately reproduced. The stability thresholds in time of the proposed method depend linearly on h and inversely on p. Therefore, the proposed collocation Newmark-IGA method retains the good convergence and stability properties of standard Galerkin IGA and spectral element discretizations of acoustic problems, as well as the high computational efficiency of collocation methods due to matrix sparsity and fast function evaluation.

An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. Part II: Subcell finite volume shock capturing

The second paper of this series presents two robust entropy stable shock-capturing methods for discontinuous Galerkin spectral element (DGSEM) discretizations of the compressible magneto-hydrodynamics (MHD) equations. Specifically, we use the resistive GLM-MHD equations, which include a divergence cleaning mechanism that is based on a generalized Lagrange multiplier (GLM). For the continuous entropy analysis to hold, and due to the divergence-free constraint on the magnetic field, the GLM-MHD system requires the use of non-conservative terms, which need special treatment.Hennemann et al. (2020) [25] recently presented an entropy stable shock-capturing strategy for DGSEM discretizations of the Euler equations that blends the DGSEM scheme with a subcell first-order finite volume (FV) method. Our first contribution is the extension of the method of Hennemann et al. to systems with non-conservative terms, such as the GLM-MHD equations. In our approach, the advective and non-conservative terms of the equations are discretized with a hybrid FV/DGSEM scheme, whereas the visco-resistive terms are discretized only with the high-order DGSEM method. We prove that the extended method is semi-discretely entropy stable on three-dimensional unstructured curvilinear meshes. Our second contribution is the derivation and analysis of a second entropy stable shock-capturing method that provides enhanced resolution by using a subcell reconstruction procedure that is carefully built to ensure entropy stability.We provide a numerical verification of the properties of the hybrid FV/DGSEM schemes on curvilinear meshes and show their robustness and accuracy with common benchmark cases, such as the Orszag-Tang vortex and the GEM (Geospace Environmental Modeling) reconnection challenge. Finally, we simulate a space physics application: the interaction of Jupiter's magnetic field with the plasma torus generated by the moon Io.