Towards hybrid meta- modelling earthquake scenarios : improving high-fidelity numerical simulation with machine learning techniques.

GENCI Demande d’Attribution de Ressources Informatiques Description scientifique du projet Titre du projet : Towards hybrid metamodelling earthquake scenarios : improving high-fidelity numerical simulation with machine learning techniques. Num´ero du projet DARI 1: A0060410444 Responsable scientifique 2(nom, pr´enom) : LOPEZ-CABALLERO Fernando (Maˆıtre de Conf´erence/LRU) - [email protected] Laboratoire : CentraleSup´elec - Laboratoire MSSMat UMR CNRS 8579 Nombre d’heures demand´ees (Cpu mono-processeur) sur le projet : CINES BULL noeuds fins Occigen : 2.500.000 heures scalaires 1 R´esum´e The objective of this project is to predict the dynamic behaviour of strategical structures (such as nuclear power plants) and to assess the seismic response of large seismic-prone urban areas, at the occurrence of extreme ground shaking events. To this end, we have recently developed an efficient multi-tool computational platform, capable of simulating broad-band nonlinear seismic wave propagation in highly heterogeneous media, from the fault to the aboveground structures (i.e. at a regional scale). The core of the platform is represented by SEM3D, a Spectral Element (SE) Method code tailored to solve 3-D wave propagation in anisotropic, randomly heterogeneous and non-linear geo-materials. The platform efficiently handles the presence of complex surface topography, coastlines and bathymetry, 3-D geological configurations along with large fault discontinuities [Gatti˙et˙al˙2017a,21, 24,25,26]. The wave-propagation code has previously shown very good scalability in mesh featured by millions of hexahedral elements. For instance, we studied the seismic response of the Kashiwazaki-Kariwa Nuclear Power Plant (KKNPP), during the MW6.6 2007 Niigata earthquake (central west-coast of Japan) reaching a maximum frequency fmax of 7 Hz, in a computational domain of 60 km ×60 km×60 km, with a minimum grid size of 200 m ×200 m × 200 m, with minimal shear wave velocity of 700 m/s [25]. Similar earthquake 1. Uniquement en cas de demande de prolongation d’un projet existant. 2. Le responsable se charge du suivi du projet et fournit un bilan en fin d’ann´ee. 1 scenarios are currently being run for the Argostoli site (a European test site, where a vast geophysical and tectonic survey are currently carried out), within the framework of the ongoing A0040410444 call. The performances of SEM3D have been checked within this ongoing project (see activity report) : for the Argostoli test case, we managed to perform our largest physics-based high-fidelity numerical simulation of a main shock earthquake scenario, reliable up to a maximum frequency of 10 Hz, in the softest part of the domain (i.e. where shorter wave lengths must be correctly discretized). However, real case studies are featured by very soft soil deposits (typically reaching values of 200 m/s), requiring higher spatial resolution for the same maximum frequency propagated (the Argostoli test case is featured by a minimum shear-wave velocity of 650 m/s). Furthermore, the modal response of rigid structures (for instance the reactor buildings) shifts towards higher frequencies (≥10 Hz), that should be excited by broad-band synthetic input motion (for instance, for a Soil-Structure Interaction model at the site scale). Finally, the possibility of modelling the heterogeneous spatial fluctuations and uncertainty of the input parameters in a non-linear wave propagation context (due to the development of a very efficient random field generator over large domains, coupled with SEM3D), permits the modelling of complex scattering effects taking place at higher frequencies (i.e. higher spatial resolutions) and leads to explore different realizations of the same earthquake scenario. Therefore, we consider in this project an extension of the mentioned analyses (see, for instance, [25,26,27]). We have at least three tests cases in our disposal : the region surrounding the KKNPP, the Argostoli basin and the Volvi basin (another experimental site in Greece). This achievement will be exploited to improve probabilistic vulnerability assessment, using uncertainty propagation techniques (construction of meta-model simulations and fragility curves). As a matter of fact, high-fidelity earthquake simulations can be exploited to either construct brand-new site-specific databases including a great number of synthetic ground motions, either to integrate existing recording databases, to be exploited for several purposes. In other words, the results of the set of cumbersome simulations foreseen in this project will be a suitable basis for surrogate modelling (performed externally, in a posterior phase, by applying, for instance, machine learning techniques such as ANN2BB [29]), so as to reduce the computation time associated keeping an accurate prediction. However, in doing so, the experimental design (ED) is constructed by SEM3D analyses. 2 2 Pr´esentation g´en´erale For a long time, geophysicists have concentrated on wave propagation in simplified geological media (e.g. considering the Earth’s crust as subhorizontally layered visco-elastic half-space) neglecting crucial issues, such as : — the heterogeneities in the Earth’s crust (at regional and site scale) — the non-linearity of shallow soil deposits — the surface topography — solid-fluid interaction (coastlines, bathymetry etc) Although the mentioned numerical models were capable of reproducing reasonably well the P, S or Rayleigh waves arrival times and the main reflections within the Earth’s crust, they were generally limited to the very lowfrequency part of the radiated spectrum. Moreover, very rare are the examples of a continuous fault-to-structure interaction studies. The most reliable and efficient method was proposed by Bielak et al. [8], called the Domain Reduction Method. It consists into a two-step analysis, declined as follows : (i) regional scale wave-propagation in simplified geological medium (typically the deep visco-elastic stratified Earth’s crust), not including the non-linear site-effects ; (ii) a second run of the analysis on a smaller domain (eventually adding the structure) delimited by an artificial boundary at which equivalent inertial forces and free-field velocities, obtained in the precedent step, are applied. Quinay et al. [13] applied this method to study the structural response of the Kashiwazaki-Kariwa Nuclear Power Plant (KKNPP) in a fault-to-structure framework, up to 1 Hz. Those aspects clashes with the need of broad-band realistic synthetics to input into realistic structural models : the transient dynamics of aboveground structures (eventually with soil-structure interaction at the site scale) has been traditionally studied separately, by employing selected spectrumcompatible recordings. In any case, the seismic design of critical structures requires refined analyses providing synthetic ground motion time-histories reliable at higher frequencies and reproducing complex phenomena such as the coda in the seismic signals [2], the directivity and the incoherence of the wave motion near-source [24], the non-linear site effects [20]. The ever increasing availability of computer power paves the way to deterministic exploration of the scenarios space, compatibly with the degree of knowledge of focal issues, such as : 1. the 3-D geological structure of the Earth’s crust 2. the geomechanical properties of soil deposits and deep bedrock 3 3. the surface morphology (topography, bathymetry, coastlines) 4. the characteristics of the active faults (tectonic context, slip patches, focal mechanisms) The level of uncertainty related to this piece of information, convolved together with the numerical dispersion, determines its predictive fidelity of such physics-based analyses. Researchers are trying to address those issues, mainly through full-waveform inversion [16,22], noise correlations [Book˙Fouque˙et˙al˙2008˙Wave˙Propagation˙Random˙Layered], probabilistic models of the fluctuations [6,7], complex rheological models. In this project, instead, we try to employ the efficient multi-tool platform developed so far to reproduce numerically those observations, when considering high resolution wave propagation problem using conformal hexahedra, which are necessary for the widely used spectral finite element method [Komatitsch˙Tromp˙I˙2002,5,9]. Moreover, we aim at calibrating those models with the structural dynamics counterpart [13], for forward vulnerability assessment. Here, high resolution means meshes containing billions of elements [11,14,15]. Performing wave propagation in such meshes becomes a possibility through the unprecedented development of powerful computational resources and scalable parallel solvers. However, constructing efficiently such discretizations and running the simulations on such meshes is not an easy task. Recently some research groups have shattered the barrier of running analyses on very large meshes, exhausting all resources on very large machines all over the world. An impressive result was obtained by a Chinese research group, who run an non-linear earthquake simulation over a continental 320 km by 312 km by 40 km region, up to 18 Hz [19]. The use of octrees has made possible a high degree of scalability in mesh generation using up to 220.000 cores [11]. To tackle the mentioned obstacles, recently, we have developed a series of computational tools that help in constructing large scale numerical scenario of an earthquake event, namely : — a numerical technique that allows to generate conformal hexahedral meshes using an octree-based method similar to [12] ; — an efficient random field generator over large computational domains [28], based on a fully parallelized version of the spectral representation technique [4] ; — a parallel non-linear wave-propagation solver based on the spectral element method and coupled with the random field generator [21]. The parallel octree-based mesh generator is able to treat immersed geometries with high scalability in executions. It was tested up to 6561 processing cores. The mesh can be non-conformal, composed exclusively by hexahedral elements, or it can be made conformal using tetrahedral elements. Non4 conforming meshes have hanging nodes, (i.e. nodes located inside edges or faces), which restricts its use to numerical methods that provide support for this type of nodal interpolation. In particular, efficient parallel implementation of the SE Method (SEM) requires the use of conformal hexahedra. Owing to the SEM spectral convergence, higher resolution is obtained by prefinement (i.e. by increasing the polynomial order) rather than by increasing the number of linear elements (common practice in Finite Difference/Element schemes). The random field generator employs the sum of a series of Nϕcosines with random phases/amplitudes [4], sampled over fine regular grids [21]. The Fast Fourier Transform can be used to bring the complexity of this generation method to O(Nϕlog Nϕ). The field is then re-interpolated over the not uniformly space-distributed SE grid (featured by 5 to 10 Gauss-LobattoLegendre points per dimension). When dealing with large domains, the scalability issue is solved by generating smaller independent realizations, supported on overlapping subdomains (in a distributed memory parallel scheme), and then recomposed together into the full domain. Finally, the non-linear wave-propagation solver exploits the high resolution of the spectral element method with efficient explicit algorithm to integrate non-linear rheology representing the non-linear cyclic behaviour of soils. 3 M´ethode 3.1 M´ethode num´erique The computational platform we are intended to test in this project bares upon three numerical codes : — SEM3D : an efficient Spectral Element code tailored for wave-propagation in non-linear heterogeneous geo-materials at the regional and continental scales [21] ; — HexMesh 3: a scalable octree-based mesh generation code ; — randomField 4: a scalable random field generator over large computational domains, based on the spectral generation technique [3]. The first code is the center of our attention in this project and has already run efficiently on large scale machines, reproducing real earthquake scenarios with high accuracy. The second and third codes are used in conjunction with the spectral element code, so to mesh correctly the complex morphology of the Earth’s crust structure and its rheology with natural spatial fluctuations. 3. https://github.com/jcamata/HexMesh.git 4. https://github.com/cottereau/randomField.git 5 Non-linear wave propagation : SEM3D We develop SEM3D [18] jointly by the Institut de Physique du Globe de Paris, CEA and CentraleSup´elec, mainly for applications in geophysics, with the aim of considering realistic 3-D complex geo-structures. It is written in Fortran 95, and uses MPI. The I/O is performed using parallel HDF5 libraries. The meshes are non-structured and composed of conformal hexahedra by employing METIS 5library. The discretization in space is based on spectral elements in space and an explicit scheme in time (both for elastic and non-linear case). The spectral element method is a variant of the finite element methods which uses high-order Lagrange polynomials on GaussLobatto-Legendre (GLL) points. The corresponding quadrature ensures that the mass matrix is diagonal. Hence, inversion of the mass matrix at each time steps is not costly and very short time steps, imposed by the explicit scheme, are not an issue. The code was tested on several clusters, including Curie (by CEA) and the M´esocentre de CentraleSup´elec and ENS Paris Saclay (SGI ICE X) and GENCI Occigen (c20150447417 and the ongoing A4 allocation project A0040410444). Moreover, it has been run on the 64-nodes Bi-Opteron cluster of Service de Calcul Parall`ele de l’Institut de Physique du Globe de Paris (a model with 25 millions degrees of freedom, running over 128 cores). A previous version of the code (less efficient in terms of parallelism [23]) was tested on Jade, in a previous GENCI project (c2010046485). On these occasions, the weak scalability of the spectral element solver has been demonstrated. Finally, it should be noted that a code with similar features and structure (although developed by a completely different team) was ported on the largest clusters available (see for example [9], 15 years back). An example of weak scalability curve obtained for SEM3D is shown in Figure 1. SEM3D has been run over more than 4096 cores. Moreover, given the experience of allocation A4 (with mid-term correction), a comparison between SEM3D performances obtained on FUSION (cluster of M´esocentre Moulon) and OCCIGEN were compared between each other, to check our estimations in terms of required CPU-hours. The numerical models of the Niigata earthquake (scenario 1) and Argostoli basin (scenario 2) were exploited. Ten simulations are hereafter compared : F1-F5, run on FUSION, O1-O5, run on OCCIGEN (see the specifics in Table 1). Fig. 2shows the performances obtained, compared to the original (CP) and corrected-blind (BP) prediction (made before allocation A4). First and foremost, we noticed a very good agreement between the corrected prediction CP and F1 and BP and F4, which justifies our prediction made on FUSION. However, we needed some preliminary iterations (corresponding to simula5. http://glaros.dtc.umn.edu/gkhome/views/metis 6 Figure 1 – Weak scalability curves obtained for SEM3D, expressed as number of time rate of mesh elements processed varying the number of MPI process employed. In blue, the linear scalability curve represents the ideal optimum result (perfectly linear scalability). Green and red lines portray the worst and average performances instead. tions O1-O3) to replicate the same performances obtained on FUSION (in terms of CPU-time). Specifically, preliminary simulations O1-O3 (performed without the CINES support, before the mid-term allocation) showed a increased computation burden (expressed as tCP U /tEQK ). For instance, O2 and F4 represents the same simulation, although we obtained a CPU-time per second real time simulation tCP U /tEQ [h/s] of 1.82×103h/s for FUSION and 3.49×104h/s for OCCIGEN respectively. Iteration O4 was performed thanks to the CINES support : the same performances were obtained for FUSION and OCCIGEN on the same test case. This iteration validates SEM3D portability on OCCIGEN and allowed us to obtain mid-term extra hour allocations. In the following, we transitioned towards a higher number of DOFs by improving the numerical model of Argostoli (including the new code feature represented by the extended fault seismic source) previously run on FUSION (F5). On OCCIGEN, the same numerical model was considered (called O5), its accuracy being increased by p−refinement, i.e. by increasing two times the GLL points per edge (resulting in 10×10×10 GLL points per element), for a total ≈13.48 ·109. Expected performances were attained (quasi-linear scalability) and we successfully passed the limit of 1010 DOFs (never before achieved on FUSION), thanks to the possibility of running jobs over >720 CPU - cores (current limitation for FUSION). 7 10 810 910 10 10 11 DOFs [1] 10 1 10 2 10 3 10 4 tCP U /tEQK [h/s] F=Fusion; O=Occigen O1 O2 F1 F2 F3 F4 F5 Figure 2 – CPU-time per real earthquake simulation tCP U /tEQK , for different DOF numbers. BP and CP correspond to the blind and correctedblind predictions respectively (see midterm allocation). NEL NGLL NMP I F1 2×1065 648 F2 - 7 648 F3 - 10 648 F4 18×1065 648 F5 4.5×1065216 O1 - 5 4800 O2 - 5 648 O3 - 7 648 O4 - 5 4800 O5 4.5×10654000 CP 2×1065 216 BP 18×1065 648 Table 1 – Earthquake numerical model. NE: number of hexaedrael elements (8-nodes) ; NGLL : number of GLL points per element) ; NMP I : number of CPU cores. Octree-based Mesh Generation : HexMesh Octrees are hierarchical data structures that allows the decomposition a three-dimensional space in regular cubes, called octants. The initial octant normally surrounds all the input domain and is divided in 8 equally spaced cubes - the child octants. Each of these children could be subdivided again, in a recursive process that is usually limited by some criteria (i.e. maximum number of subdivisions, minimum octant size). An octree can be implemented by a tree structure or by a linear octree, where only octants with no children are stored. In this last representation, as also known by linear octree, each leaf is stored as a unique integer - called locational code - which identifies its position and depth level in the tree. Among the benefits of linear octrees is the implicit representation of intermediary octants (non-leaves), making it memory efficient, and better performance on sequential access. Additionally, by not making use of memory pointers, parallel implementation of the linear octree has a lower overhead in interprocess communication (in comparison to other tree structures [10]). We have implemented a set of algorithms that address the special needs of parallel octree meshing. These algorithms are : (i) Octree initialization algorithm that guarantees that all processes have at least one piece of the linear octree. (ii) Refinement algorithm that decompose 8 the octants according to some refinement criterion (i.e. maximum size or position relative to the geometry boundaries), and 2 :1 balancing algorithm that ensure that no neighbor octants differentiate in more than one level [10]. This last algorithm is important because a balanced octree results in more smooth transitions between elements, a feature that has a direct influence on the final mesh quality. It also guarantees that the resulting mesh will not have more than one hanging node for each edge or face, reducing the complexity of conforming mesh process. The resulting code has been written in plain C and makes use of MPI only. The code runs on any standard cluster. A previous version of the code (generating non-conformal hexahedra and conformal tetrahedra meshes) has been tested at Stampede Dell cluster (at TACC) and on a SGI ICE 8400 and on a Oracle cluster in Brazil. Table 2 presents a weak scalability test performed on the mesh of the Kefalonia region (Greece). The communication time (last column) attains a fairly constant percentage of the total time. The code has run over 6561 MPI cores. Cores Levels Nodes Elements Time (s) Comm(%) 81 7 85,602,744 83,102,679 12.061 21 % 243 8 769,790,232 747,937,476 32.994 20 % 729 9 6,926,153,724 6,731,438,013 105.188 25 % 6561 10 62,330,385,168 60,583,119,264 123.066 29 % Table 2 – Weak scaling of the mesh generation for the Kefalonia region in Greece. Cores : number of MPI cores ; Levels : octree refinement level ; Nodes : number of mesh nodes ; Elements : number of mesh elements ; Time : generation time ; Comm : percentage of the total elapsed time for communication operation. Random Field Generation : randomFields The heterogeneous properties of the Earth’s crust are included in our modelling strategy for large scale earthquake scenarios. In this context, by large we refer to a domain size Lmuch larger than both the correlation length ℓC(or some characteristic size of the fluctuation) and the discretization step h. It is therefore advantageous to represent those fluctuations a scalar random fields (for instance, representing the spatial distribution of the shear-modulus). Those large random fields can be effectively sampled over a coarse grid (with a step size relevant for the correlation length) and then interpolated onto the mesh of interest (provided by the GLL (GaussLobatto-Legendre) grid used in SEM3D). If the discretization step is much larger than the correlation length, the sampling becomes simple and numerically inexpensive. Indeed, for the mesh considered, the random field is 9 [19] Haohuan Fu et al. “18.9Pflopss Nonlinear Earthquake Simulation on Sunway TaihuLight: Enabling Depiction of 18-Hz and 8-meter Scenarios”. In : Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. SC ’17. Denver, Colorado : ACM, 2017, 2:1-2:12. doi :10.1145/3126908.3126910. url :http://doi.acm.org/10.1145/3126908.3126910. [20] F. Gatti, F. Lopez-Caballero, R. Paolucci et D. Clouteau. “Near-source effects and non-linear site response at Kashiwazaki-Kariwa Nuclear Power Plant, in the 2007 Chuetsu-Oki earthquake: evidence from surface and downhole records and 1D numerical simulations”. In : Bulletin of Earthquake Engineering 16.3 (2017), p. 1105-1135. doi : 10.1007/s10518-017-0255-y.url :https://doi.org/10.1007/ s10518-017-0255-y. [21] F. Gatti, L. De Carvalho Paludo, A. Svay, F. Lopez-Caballero, R. Cottereau et D. Clouteau. “Investigation of the earthquake ground motion coherence in heterogeneous non-linear soil deposits”. In : Procedia Engineering 199.Supplement C (2017). X International Conference on Structural Dynamics, EURODYN 2017, p. 2354-2359. doi :10 . 1016 / j . proeng . 2017 . 09 . 232.url :http : / / www . sciencedirect.com/science/article/pii/S1877705817336822. [22] Hiroshi Kawase, Shinichi Matsushima, Fumiaki Nagashima,Baoyintu et Kenichi Nakano. “The cause of heavy damage concentration in downtown Mashiki inferred from observed data and field survey of the 2016 Kumamoto earthquake”. In : Earth, Planets and Space 69.1 (2017), p. 3. doi :10 . 1186 / s40623 - 016 - 0591 - 1.url :https : //doi.org/10.1186/s40623-016-0591-1. [23] S. Khazaie, R Cottereau et D Clouteau. “Numerical observation of the equipartition regime in a 3D random elastic medium, and discussion of the limiting parameters”. In : Computers & Geosciences 102 (2017), p. 56-67. doi :10.1016/j.cageo.2017.02.007.url :http:// www.sciencedirect.com/science/article/pii/S0098300417301747. [24] Angkeara Svay et al. “Spatial coherency analysis of seismic ground motions from a rock site dense array implemented during the Kefalonia 2014 aftershock sequence”. In : Earthquake Engineering & Structural Dynamics 46.12 (2017). eqe.2881, p. 1895-1917. doi :10.1002/eqe. 2881.url :http://dx.doi.org/10.1002/eqe.2881. [25] F. Gatti, F. Lopez-Caballero, D. Clouteau et R. Paolucci. “On the effect of the 3-D regional geology on the seismic design of 16 critical structures: the case of the Kashiwazaki-Kariwa Nuclear Power Plant”. In : Geophysical Journal International 213.2 (2018). doi: 10.1093/gji/ggy027, p. 1073-1092. doi :10.1093/gji/ggy027.url : https://doi.org/10.1093/gji/ggy027. [26] F. Gatti, S. Touhami, F. Lopez-Caballero et D. Pitilakis. “3-D source-to-site numerical investigation on the earthquake ground motion coherency in heterogeneous soil deposits”. In : 9th European Conference on Numerical Methods in Geotechnical Engineering, 25-27 June 2018, Porto (Portugal). 2018. [27] F. Gatti et al. “Broad-band 3-D earthquake simulation at nuclear site by an all-embracing source-to-structure approach”. In : Soil Dynamics and Earthquake Engineering 115 (2018). doi: 10.1016/j.soildyn.2018.08.028, p. 263-280. doi :10.1016/j.soildyn.2018.08.028.url :http:// www.sciencedirect.com/science/article/pii/S0267726118303890. [28] L Paludo, V Bouvier et R Cottereau. “Scalable parallel scheme for sampling of Gaussian random fields over very large domains”. In : Pre-print (2018). [29] R. Paolucci, F. Gatti, M. Infantino, A. G. Ozcebe, C. Smerzini et M. Stupazzini. “Broad-band ground motions from 3D physicsbased numerical simulations using Artificial Neural Networks”. In : Bulletin of the Seismological Society of America 108.(3A) (2018), p. 12721286. doi :10.1785/0120170293. [30] Roberto Paolucci, Maria Infantino, Ilario Mazzieri, Ali G¨uney ¨ Ozcebe, Chiara Smerzini et Marco Stupazzini. “3D Physics-Based Numerical Simulations: Advantages and Current Limitations of a New Frontier to Earthquake Ground Motion Prediction. The Istanbul Case Study”. In : Recent Advances in Earthquake Engineering in Europe: 16th European Conference on Earthquake Engineering-Thessaloniki 2018. T. 46. Springer. 2018, p. 203. [31] Sara Touhami, F Gatti, Fernando Lopez-Caballero, Fabrice HOLLENDER et Edward Marc Cushing. “Argostoli site, from site test to numerical model: a holistic approach”. In : Best Practices in Physics-based Fault Rupture Models for Seismic Hazard Assessment of Nuclear Installations: issues and challenges towards full Seismic Risk Analysis. Cadarache, France, mai 2018. url :https://hal.archives-ouvertes. fr/hal-01851742. 17