Innovative multiscale numerical algorithms for highly heterogeneous media extended to seismic problems

There has been an extensive development of massively parallel computer architectures in the last decade. With single processors limited by technical issues such as heating, computers have instead been built to leverage large numbers of processors (grouped in cores) of mild speed and storage capacities. This new paradigm has led to a revision of what is expected from simulators from the viewpoint of numerical algorithms. Precision and robustness remain fundamental properties of numerical methods for extreme-scale computational science. Nevertheless, the underlying algorithms must be naturally shaped to cope with the new massively parallel architectures. In particular, failure is a certainty in long-time simulations on these new generations of parallel architectures. Consequently, numerical methods for extreme-scale computing should induce asynchronous and communication-avoiding algorithms to mitigate failures in long-term simulations better.

Elastic wave propagation in heterogeneous media is at the core of the modern techniques for seismic imagery. Taking part in the solution of inverse problems, the direct simulation of wave propagation in three-dimensional heterogeneous domain represents a big challenge for standard numerical methods. In a broad sense, very fine meshes must be adopted to approximate the high-frequency fi eld distributions as they are greatly impacted by the multiscale structures or the high contrast of the medium. (See Figure 1 for an illustration.)

Fig. 1
Figure 1. A synthetic heterogeneous medium from the HPC4E seismic test suite.


The usual idea of employing high-order polynomial interpolations on coarse meshes is not a feasible option since interfaces between different material media and faces of the partition must coincide to avoid spurious numerical modes. Consequently, numerical methods for long-time simulations require a fi ne mesh to fi t the complex geometry with a high-order approximation to minimize dispersion. The computational cost involved in such simulations is prohibitive. Therefore, wave propagation simulations for seismic imagery turn out to be the perfect prototype to assess new numerical algorithms for the forthcoming generation of architectures, the exascale computers.

Numerical methods built on the "divide-and-conquer" philosophy satisfy the architectural imperatives of the new generation of high-performance computers better than classical methods operating only on the finest scale of the discretization. Indeed, splitting the computation of extreme simulations into a set of independent problems of smaller size is a way to circumvent faults and to allow spatial and time data locality while taking full advantage (in terms of performance) of the granularity of the new generation of computer architectures.

In this context, the Multiscale Hybrid-Mixed (MHM) finite element method [1-8] appeared as an attractive "divide and conquer" option to handle heterogeneous problems. Overall, the idea relies on basis functions specially designed to upscale sub-mesh scales to an overlay coarse mesh. (See Figure 2 for an illustration on a porous media fi eld and Figure 3 for an example of a multiscale basis function.) Adapted to handle multiscale or high contrast coefficients on coarse meshes, the MHM method permits face-crossing interfaces endowed in local boundary conditions. It gives rise to a staggered algorithm, which first computes the multiscale basis functions from completely independent element-wise problems, and next obtains the degrees of freedom on faces from the global face-based formulation defi ned on a coarse mesh. As a result, the MHM method naturally embeds an upscaling procedure which is suitable to simulate wave propagation problems in highly heterogeneous media on coarse meshes.

Fig. 2
Figure 2. A coarse mesh overlying a highly heterogeneous domain [1]. Edges
are not aligned with jumping coefficients.


Fig. 3
Figure 3. A typical multiscale basis function on a triangular element [4].


To provide a detailed solution of multiscale models, the MHM method is accompanied by an a-posteriori error estimator which adds degrees of freedom only where they are really necessary. Such an estimator is closely related to the structure of the MHM method and inherits its intrinsic exibility in incorporating local space interpolations and local time stepping. As a result, it relies only on the numerical error on the faces of the coarse partition. This association produces parallel numerical algorithms shaped to handle the complex multiscale phenomena which characterize the elastostatic and elastodynamic models with rough coefficients. We illustrate in Figure 4 a typical space adaptation process. Observe that the topology of the coarse mesh stays unchanged throughout the adaptation process thus avoiding cumbersome re-meshings.

Fig. 4
Figure 4. Space adaptation process [4].


In Figure 5, we depict the initial and final distribution of the degrees of freedom associated with a coarse mesh and the resulting solutions from the action of the a-posteriori error estimator.

Fig. 5
Figure 5. An elasticity test case. The red dots represent edge partitions. The coarse mesh with the initial distribution of the degrees of freedom (top left), and its associated solution (top right). The final distribution of the degrees of freedom (bottom left) and its associated solution which is driven by the error estimator (bottom right).


The new multiscale domain decomposition algorithms have been applied to the 2D and 3D elastodynamic models in the scope of the HPC4E project. It was assumed that the coefficients may be highly heterogeneous or present high contrast ratio in accordance to the HPC4E seismic test suite. Preliminary results were performed and validated the concept of these algorithms in terms of mathematical (i.e. accuracy) and computational (i.e. scalability) terms. Speci cally, the three-dimensional elasticity model is adopted to study the computational performance of the multiscale algorithms in a parallel architecture. This is a fundamental step towards the use of the MHM methodology to the elastodynamic model. In Figure 6, the MHM method is shown to scale nicely up to thousands of cores, and we anticipate that a similar efficiency will be attained for the elastodynamic model. In future works, we intend to demonstrate that the high heterogeneity in the fields provided by the HPC4E seismic test suite is accurately captured by the adoption of re fined local meshes, which is the case when the MHM methods reach their highest computational efficiency (see Figure 6).

Fig. 6
Figure 6. Weak scaling measurements.


Regarding the elastodynamic model, we first adopt a numerical test with analytical solution assuming homogeneous coefficients. The results highlight that the new domain decomposition algorithm produces optimal and high-order super-convergent displacement, velocity and stress, and conserve the total energy. They are depicted in Figures 7-8. A second test case for the elastodynamic model deals with the heterogeneous coefficients problem. We solve the model on a three-layered domain using a coarse mesh non-aligned with the interfaces (see Figure 9). The results show that the MHM method is, indeed, precise on coarse meshes through upscaling and can handle high contrast interfaces on non-aligned meshes. These features are presented in Figure 10.

Fig. 7
Figure 7. An elastodynamic test case with homogeneous coefficients and analytical solution. The MHM method achieves superconvergence for the displacement, velocity and stress variables (left). The discrete energy is conserved (right). The linear interpolation is used on faces.


Fig. 8
Figure 8. An elastodynamic test case with homogeneous coefficients. The absolute value
of the displacement in three different time steps.


Fig. 9
Figure 9. The non-homogeneous domain and the coarse mesh. The domaincontains three different layers.


Fig. 10
Figure 10. An elastodynamic test case with non-homogeneous coefficients.
Four snapshots of the solution using the MHM method.


We conclude that the MHM method, which is particularly adapted to parallel computing environments, emerges as a realistic option to handle three-dimensional wave propagation problems in highly heterogeneous media with the upshot to be used in the seismic imagery field.


[1] C. Harder, D. Paredes and F. Valentin A family of multiscale hybrid-mixed finite element methods for the Darcy equation with rough coefficients. J. Comput. Phys., Vol. 245, pp. 107-130, 2013.
[2] R. Araya, C. Harder, D. Paredes and F. Valentin Multiscale hybrid-mixed method. SIAM J. Numer. Anal., Vol. 51, No. 6, pp. 3505-3531, 2013.
[3] D. Paredes, F. Valentin, and H. M. Versieux. On the robustness of multiscale hybrid-mixed methods. Math. Comp., Vol. 86, No. 304, pp. 525-548, 2017.
[4] C. Harder, D. Paredes and F. Valentin. On a multiscale hybrid-mixed method for advective-reactive dominated problems with heterogenous coefficients. SIAM Multiscale Model. and Simul., Vol. 13, No. 2, pp. 491518, 2015.
[5] C. Harder, A. L. Madureira and F. Valentin. A hybrid-mixed method for elasticity. ESAIM: Math. Model. Num. Anal., Vol. 50, No. 2, pp. 311336, 2016.
[6] C. Harder and F. Valentin. Foundations of the MHM method. G. R. Barrenechea, F. Brezzi, A. Cangiani, E. H. Georgoulis Eds. Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, Springer, Lecture Notes in Computational Science and Engineering, 2016
[7] W. Pereira and F. Valentin. A Locking-Free MHM Method for Elasticity. Proceeding Series of the Brazilian Society of Computational and Applied Mathematics, 2016
[8] R. Araya, C. Harder, A. Poza and F. Valentin. Multiscale Hybrid-Mixed Method for the Stokes and Brinkman Equations - The Method. Comp. Meth. Appl. Mech. Eng., Vol. 324, pp. 29-53, 2017


Antonio Tadeu Gomes, Weslley Pereira, Frédéric Valentin - National Laboratory for Scientific Computing (LNCC) - Petrópolis - Brazil

Diego Paredes - Pontificia Universidad Católica (PUCV) - Valparaíso - Chile


This article also appeared in the LinkedIn GroupJoin us!

Other LinkedIn articles:

Link Wind energy in HPC4E
Link Investigation of flame structure of biomass-derived gaseous fuels
Link Finding the "not so easy" oil
Link How supercomputing can help improving the energy sector
Link Effects of fuel composition on biogas combustion in premixed laminar flames
Link New generation subsurface imaging gets a boost from HPC
Link Novel Hybridizable Discontinous Galerkin method paves the way to tackle realistic 3D problems in seismic imaging
Link Improving short-range wind intensity prediction based on multimodel meteorological ensemble forecasts and Genetic Programming
Link Biogas Utilisation for Sustainable Power Generation
Link Preparing the oil and gas industry for the Exaflop era
Link Dynamical and statistical high resolution downscaling approaches for the surface wind