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 field distributions as they are greatly impacted by the multiscale structures or the high contrast of the medium. (See Figure 1 for an illustration.)

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 fine mesh to fit 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 field 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 defined 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**.

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.

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.

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**. Specically, 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 refined local meshes, which is the case when the MHM methods reach their highest computational efficiency (see Figure 6).

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.

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**.

**References**

[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 Group. **Join us!**

**Other LinkedIn articles:**