Ricard Borrell
email: ricard.borrell@bsc.es

Juan Carlos

García Cajas

Daniel Mira

Ahmed Taha

Seid Koric

Mariano Vázquez

Guillaume Houzeaux

J C Cajas

Parallel mesh partitioning based on space filling curves

Keywords: space filling curve, SFC, mesh partitioning, geometric partitioning, parallel computing

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Mesh partitioning is traditionally based on graph partitioning, which is a well-studied NP-complete problem generally addressed by means of multilevel heuristics composed of three phases: coarsening, partitioning, and un-coarsening. Different variants of them have been implemented in publicly available libraries such as Metis/ParMetis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF], Scotch/PT-Scotch [START_REF] Chevalier | Pt-scotch: A tool for efficient parallel graph ordering[END_REF] and Zoltan [START_REF] Devine | Zoltan data management services for parallel dynamic applications[END_REF]. All these libraries enable parallel partitioning, but with a limited parallel performance and a decreased quality of the parallel partition [START_REF] Karypis | A parallel algorithm for multilevel graph partitioning and sparse matrix ordering[END_REF][START_REF] Shang | Performance analysis of large scale parallel cfd computing based on codesaturne[END_REF]. Both aspects make graph-based partitioning a potential bottleneck in the simulation workflow. However, taking into account the evolution of the computing HPC systems, any potential bottleneck becomes an effective bottleneck if not addressed in time. Motivated by these circumstances, we present a fully parallel geometric partitioning alternative.

Geometric partitioning techniques obviate the topological interaction between mesh elements and perform its partition according to their spatial distribution. If we consider as unitary element the mesh cell, then its mass center can be used to determine its spatial location. A Space Filling Curve (SFC) is a continuous function used to map a multi-dimensional space into a one-dimensional space with good locality properties, i.e. it tries to preserve the proximity of elements in both spaces. The idea of geometric partitioning using SFC is to map the mesh elements into a 1D space and then easily divide the resulting segment into equally weighted sub-segments. A significant advantage of the SFC partitioning is that it can be computed very fast and it is easy to parallelize, especially when compared to graph partitioning methods. However, while the load balance of the resulting partitions can be guaranteed, the data transfer between the resulting subdomains, measured in terms of edge-cuts in the graph partitioning approach, cannot be explicitly measured and thus neither be minimized.

In this paper we describe and evaluate the performance of a parallel SFC based partitioner. We have avoided any computing or memory bottleneck that could limit the scalability of the algorithm, imposing that the solution achieved is independent (discounting rounding off errors) of the number of parallel processes used to compute it. The structure of the paper is as follows: Section 2 contains a short overview of the Hilbert SFC; in Section 3 we describe in detail the sequential and parallel versions of our SFC-based partitioner. In Section 4 the performance of the initial implementation is analyzed and some optimizations are evaluated; moreover, the partitions quality is compared with the quality of the equivalent ones computed with Metis. Finally, general conclusions are outlined section 5.

Hilbert Space Filling Curve

There are many possible definitions of SFC based on different mapping options, among them the well-known Peano [START_REF] Peano | Sur une courbe, qui remplit toute une aire plane[END_REF] and Hilbert [START_REF] Hilbert | Über die stetige Abbildung einer Linie auf ein Flächenstück[END_REF] approaches. A complete overview is given in [START_REF] Sagan | Space-Filling Curves[END_REF]. In this paper the Hilbert SFC is selected, however switching to another mapping option would be straightforward in our implementation.

The Hilbert SFC is defined by means of a geometric recursion. For the 2D case, on the p th level of the recursion a discrete function is obtained:

h : {1, ..., 2 2p } -→ {1, ..., 2 p }{1, ..., 2 p } (1)
determining the ordering of the curve through a 2 p × 2 p cartesian grid. The four initial steps of this recursion are represented in Figure 1(top). On the first step, the curve has a " " shape traversing a 2 × 2 grid. This is refined into the curve that traverses the 4 × 4 grid corresponding to the second step. This 4 × 4 grid can be decomposed into four 2 × 2 sub-grids, where the curve takes the shapes " ", " ", " " and " ", respectively. Following a specific refinement pattern for each of these basic 2 × 2 shapes we obtain the third level, and this process is continued until the desired level of refinement is achieved. Therefore, the process of generating a Hilbert curve, is based on recursively applying a specific refinement pattern to the basic 2 × 2 shapes generated at each level of the recursion. This algorithm can be implemented by cyclic lookups to two arrays: the first storing the orderings which define the basic 2 × 2 shapes, the second its refinement patterns, also referred as orientations. Details of this implementation approach can be found in [START_REF] Luitjens | Parallel space-filling curve generation through sorting[END_REF]. The same idea is followed in the 3D case, where the basic shapes are defined on 2 × 2 × 2 sub-grids. The first four steps of the 3D recursion are depicted in Figure 1 (bottom). Note that a Hilbert SFC can be defined on any rectangular or cubic domain, by using a grid with 2 p elements on each axis. A Hilbert SFC can also be generated in parallel following a multilevel approach. An initial coarse grid is defined, and the SFC-index of each of its elements, hereafter referred as coarse bins, determines the rank of the parallel process continuing the recursion within it. If the proper orientation is followed within each coarse bin, the resulting SFC generated in parallel is the same that would be obtained sequentially. This process is illustrated in Figure 2. Generating an SFC in parallel has two advantages: first it is faster and second, since there is more memory available for the overall SFC, a more refined one can be generated. Recalling that the ordering of the coarse bins is based on an SFC, this imposes a restriction to the number of parallel processes, P , generating the SFC: P , must be a power of 2 d , where d ∈ {2, 3} is the dimension of the domain.

SFC based mesh partitioning

Sequential version

Consider the problem of partitioning a mesh in M subsets using the Hilbert SFC described in Section 2. A sequential algorithm to obtain the desired partition is as follows:

1. Define a minimum bounding box enclosing the mesh. 2. Define a cartesian grid with 2 p elements in each direction, its cells are referred as bins. Note that p recursion steps will be necessary to generate the corresponding Hilbert SFC. 3. Embed the mesh elements into the grid bins. For each mesh element a weight is added to the grid bin containing it. A mesh partition can be based on the distribution of nodes (vertices) or cells, in this paper we consider cells. We use the mass center to determine the bin that contains each cell, and the weight of a cell is defined as its number of nodes. 4. Traverse the cartesian grid and label each bin with a partition number.

The SFC ordering is used to traverse the grid. The associated mapping function requires p steps for each evaluation, therefore if N = (2 d ) p the cost of this step is O(Nlog(N )). While the cartesian grid is traversed, the bins weight is accumulated on a variable and, each time the objective weight is reached, the partition number assigned to the bins is incremented and the accumulation variable is reinitialized to zero. The objective weight is initially defined as the total weight accumulated in the grid divided by M . Then it is re-evaluated each time a subset is completed, considering the remaining weight and subsets to complete, to avoid the accumulation of discrete errors. 5. Infer the mesh partition from the bins partition. This step is straightforward, the partition number of each mesh element is the one of the bin containing it.

The steps of the sequential partitioning algorithm are illustrated in Figure 3.

Parallel version

The parallel counterpart of the algorithm is described below. As initial conditions, we assume that there is a list of N c cells distributed among P parallel processes. We assume that P is a power of 2 d , i.e. P = 2 d q c , where q c are the levels of the coarse SFC recursion. The initial distribution of the mesh cells among the parallel processes could come from: i) a parallel read of the mesh file, carried out at the initialization phase of the simulation; or ii) a previous partition, when a mesh re-partitioning is considered, e.g. for load balancing purposes. The parallel SFC-based mesh partitioning algorithm has the following steps: The mesh is embedded in a cartesian grid, b) a weight is assigned to each element of the cartesian grid ("bins") depending on the number of elements contained in it, c) the cartesian grid is traversed in the SFC order, accumulating the weights of the bins and assigning subset numbers depending on the objective weights of each subset.

Define the bounding box and its cartesian grids. A minimum bounding

box is defined enclosing the mesh. This requires a collective reduction communication (MPI Allreduce), to find the maximum and minimum value for each coordinate among the mesh nodes. Then a coarse grid is defined dividing the box. Each parallel process will account for the mesh portion contained in one of the coarse bins, in particular the one with SFC-index equal to its rank. At this point it is also determined the depth of the Hilbert recursion used within the coarse bins, referred as q l . Consequently, the size of the local grids will be N l g = (2 d ) q l , and the size of the fine grid, resulting from joining the local grids, will be

N f g = (2 d ) q l q c .
2. Evaluate and redistribute of bins weight. Each parallel process evaluates the weight of the elements of its sublist accumulating it into the corresponding fine bins. Then the weights accumulated in the fine bins are redistributed such that each parallel process ends up with the ones contained in its local grid. In order to accumulate the fine bins weight while looping trough the sublist of cells, it would be required an array of size N f g , or an sparse data structure such as a hash table. Alternatively, in this paper we split this evaluation into different sub-steps. These are detailed below from the perspective of one of the parallel processes:

• Identify which coarse bins contain elements of its sublist (loop around sublist required).

• For each of these:

-Allocate a buffer of size N l g to accumulate the weights in the bins of the corresponding local grid.

-Loop around the mesh elements sublist accumulating weights on the buffer. -Store the non-zero weight bins on a new buffer (a loop around the first buffer of size N l g is required). • Collective communication (MPI Alltoall) so that each parallel process gets the size of the buffers it will receive from other processes

• Point-to-point communications (MPI Isend, MPI Irecv, MPI Wait) to re-distribute weighted bins.

• Accumulate the bins weights in the local grid. The local bins weight may have been received from other parallel processes or may have been evaluated from the same process sublist. Note that, the processes without any mesh element contained in its coarse bin will end up with an empty local grid.

3. Gather the coarse bin weights to all processes. A collective MPI Allgather communications is carried out to obtain the weights distribution across the parallel process. 4. Local grids partition. From the coarse weights distribution obtained in the previous step, each parallel process can evaluate the relative weight of its coarse bin, γ i , and the relative weight accumulated by the parallel processes with lower rank, Γ i = j<i γ j . Note Γ i ≤ 1 for all the parallel processes. Recall M are the subsets desired for the final partition, then Γ i M subsets will be completed by processes of rank lower than i. In order to generate the partition in a coherent way, the assignation of the starting subset for each parallel process has to be done carefully.

Consider the general case where Γ i M is not an integer number and define the residual subset r i as

r i = Γ i M -[Γ i M ],
where [] stands for the floor function. The starting subset for parallel process i will be [Γ i M ] + 1 and the objective weight for this subset will be scaled by the factor (1r i ), given that the parallel processes with ranks lower than i have to cover the r i fraction of it. Once the starting partition number and the objective weight for the first subset are known, the parallel process can start its local partition. SFC partitioning methods are flexible enough to generate subsets of different sizes. This situation is illustrated in Figure 4, where a mesh is partitioned by four parallel processes into five subsets. Each parallel process starts assigning subset numbers according to the coarse weights accumulated by the parallel processes of lower ranks. In the example, the first subset is filled by ranks 0 and 1, the second by ranks 1 and 2, the third by 2 and 3, and finally the fourth and fifth subsets are filled by rank 3. Note that joining the local SFC partitions generated in parallel, the result is the same that would be obtained sequentially. 5. Fine bins redistribution. The point-to-point communications (MPI Isend, MPI Irecv, MPI Wait) opposite to those made in step 2 are performed. As a result each parallel process obtains the partition number of the bins containing cells of its sublist. 6. Infer mesh partition. This step is straightforward, the partition number of each element of the sublist is the one of the fine bin containing it (received in previous step).

Note that each parallel process does not account for more than 1/P of the mesh elements of the fine grid at any step of the algorithm. Therefore, the algorithm does not present any memory or computational bottleneck. Regarding data transfers, there are three collective communications: in steps 1 (MPI Allreduce), in step 2 (MPI Alltoall) and in step 3 (MPI Allgather); and two point-to-point communications: steps 2 and 5 (MPI Isend, MPI Irecv, MPI Wait). Its influence on the parallel performance is analyzed in the following section.

Regarding the restriction of P being a power of 2 d , imposed by the Hilbert SFC, there are alternative SFC without such restriction, e.g. the Peano or Morton SFCs [START_REF] Sagan | Space-Filling Curves[END_REF]. In our case, we fix P as the largest power of 2 d minor or equal than M (the number of processes used for the simulation). An alternative option is to use a coarse grid with more bins than partitioning processes and assign more than one coarse bin to each parallel process.

The sequential and parallel executions of the algorithm produce an almost identical partition if the fine grid resulting from joining the local grids is the same than the one used for the sequential partition. However, since a larger grid can result in a better balanced partition, the parallel execution can produce more accurate results. An opposite behavior is observed with graph partitioners, where sequential and parallels results are generally different and the parallelization worsens the result [START_REF] Karypis | A parallel algorithm for multilevel graph partitioning and sparse matrix ordering[END_REF][START_REF] Shang | Performance analysis of large scale parallel cfd computing based on codesaturne[END_REF]. In case of identical grids, there can be a minor difference between the parallel and sequential SFC algorithms. This is because of the objective weight re-evaluation used to distribute the discrete errors at generating the subsets (see step 4 of the sequential algorithm). In the sequential case this re-evaluation is performed P -1 times and affects all the subsets, while in the parallel case it is performed locally within each coarse bin, the number of repetitions will depend on the portion of mesh elements partitioned by each parallel process.

Finally, repartitioning and load balancing based on SFC has a significant advantage: the partition adjustments consist in moving the cutting points along the 1D segment defined by the SFC. Therefore, changes occur in the extremes of the sub-segments obtained from the former partition and, in general, most of mesh elements can be retained in the same subdomain.

The steps of the parallel partition algorithm are outlined in Algorithm 1.

Algorithm 1 Parallel SFC-based partitioning 1. Definition of a bounding box and its cartesian grids 2. Evaluation and redistribution of fine bins weight 3. Gather coarse bin weights to all processes 4. Local grids partition based on SFC 5. Fine bins redistribution 6. Infer mesh partition

Computational tests and optimizations

In this section we present the computational tests carried out and the optimizations implemented according to the intermediate results obtained. We have considered two test cases, the first one is the simulation of a sniff flow in the respiratory system. The heterogeneous and anisotropic mesh illustrated in Figure 5 (left) is used for this case. It is mainly composed of prisms in the boundary layer and tetrahedra in the rest of the domain, while some pyramids are required on the transition between both regions. In total it has 17.7M cells. For the numerical solution of the sniff flow a large-eddy simulation (LES) formulation of the Navier Stokes equations is considered. This is coupled with a Lagrangian particle tracking solver to evaluate depositions. An illustrative image of the results is shown in Figure 5 (right). Further details about this case can be found in [START_REF] Calmet | Large-scale cfd simulations of the transitional and turbulent regime for the large human airways during rapid inhalation[END_REF][START_REF] Calmet | Subject-variability effects on micron particle deposition in human nasal cavities[END_REF]. The second case under consideration is the numerical simulation of a swirling combustor. The computational domain is also discretized with a highly anisotropic mesh composed of prisms, tetrahedra and pyramids; de-tails are shown in Figure 6 (left). The mesh is composed of 28.9M cells in total. The numerical approach for the fluid is also a LES while the turbulent combustion model is based on laminar premixed flamelets. Further details on this case can be found in [START_REF] Gövert | The effect of partial premixing and heat loss on the reacting flow field prediction of a swirl stabilized gas turbine model combustor[END_REF], an illustrative image of the flow dynamics is presented in Figure 6 (right).

The computational tests have been carried out on the Blue Waters supercomputer from the National Center for Supercomputing Applications (NCSA) at the University of Illinois (U.S.A). In particular we employed the XE nodes composed of two 16-core AMD 6276 "Interlagos" processors (nominal clock speed of at least 2.3 GHz), and connected by the Cray Gemini torus network.

The parallel SFC based partitioner developed in this paper and the physics solvers used in the test cases are both integrated in Alya [START_REF] Houzeaux | A massively parallel fractional step solver for incompressible flows[END_REF][START_REF] Houzeaux | Extension of fractional step techniques for incompressible flows: The preconditioned orthomin(1) for the pressure schur complement[END_REF],: the high performance computational mechanics code developed in the Barcelona Supercomputing Center. The physics solvable with the Alya system include: incompressible/compressible flow, solid mechanics, chemistry, particle transport, heat transfer, turbulence modeling, electrical propagation, etc. Alya is designed for massively parallel supercomputers [START_REF] Vázquez | Alya: Multiphysics engineering simulation towards exascale[END_REF], its parallelization includes both the MPI and OpenMP frameworks, as well as heterogeneous options including accelerators. Alya is one of the twelve simulation codes of the Unified European Applications Benchmark Suite (UEABS) of PRACE and thus complies with the highest standards in HPC.

Parallel performance of the partitioning process 4.1.1. Initial results

The first test carried out is the measurement of the time required for the partition of the respiratory system mesh (17.7M cells), varying the number #CPU initial version no redist.

no comm. [START_REF] Sagan | Space-Filling Curves[END_REF] 1,31 0,45 (34%) 1,27 (97%) 64 0,31 0,057 (18%) 0,30 (97%) 512 0,18 0,0078 (4%) 0,064 (36%) 4096 0,20 0,0037 (2%) 0,023 (11%) of CPU-cores from 8 up to 4096. Even when the execution with 4096 CPUs is quite unrealistic, given that a load of 4.3K cells per CPU-core is very low, it is useful to show the scalability barriers of the algorithm. Starting with 8, the number of CPUs is incremented by a factor 8 on the successive cases considered. In this way a SFC can be used to assign the coarse bins to the parallel processes. While the size of the coarse grid is multiplied by 8 in the successive cases, the local grid size is divided by 8, therefore the fine grid is kept constant. The partitions generated have 512 subsets, however, the partitioning time is independent to this number. Results are shown in column 2 of Table 1, and are depicted in Figure 8 (left), under the label "SFC initial", for its comparison with some optimizations described below. The partitioning process is very fast: a minimum time of 0.18 seconds is observed for partitioning this 17.7M elements mesh. However the scalability is very limited: with 512 CPU-cores the parallel efficiency is only 11% and 1% with 4096 CPU-cores. Despite the algorithm may be fast enough in some contexts we aim to further understand its scalability limitations, having in mind its potential usage for repartitioning and load-balancing processes.

Analysis

Following the well known Amdahl's law principle, the first step to do is to find the part of the algorithm which dominates the execution costs. First, we have skipped steps 2 and 5 of Algorithm 1, in which data redistributions are carried out. As shown in column 3 of Table 1, the rest of the algorithm takes only 34% of the partitioning time with 8 CPUs-cores, and ends up representing only 2% with 4096 CPU-cores. This percentage decrease reflects a better parallel efficiency of the part being considered, which can be calculated from Table 1: 90% with 512 CPU-cores and 24% with with 4096 CPU-cores. In conclusion, our SFC partitioning algorithms is dominated by the cost of the redistribution steps, the rest of the algorithm scales well -the scalability achieved with 4096 is good considering the low workload per CPU -and is extremely fast. Consequently, in the rest of this section we focus on the analysis and the optimization of steps 2 and 5 of Algorithm 1.

Algorithm 2

Step 2 of Algorithm 1 1. Loop around the sublist of mesh elements to count #elements to be sent to other proc. 2. For each proc. q to which messages have to be sent:

(a) Allocate buffer of size N l g (local grid of proc. q) (b) Loop around the sublist to accumulate weight on the buffer (c) Store in the "send to proc. q" buffer the bins with non-zero weight 3. MPI Alltoall to generate the communication scheme 4. Point-to-point:

(a) For each other process sending a buffer to me MPI Irecv (b) Send all the buffers evaluated in step 2, MPI Isend (c) MPI Waitall

Accumulate weight in local grid

The step 2 of Algorithms 1 is shown in detail in Algorithm 2. Step 5 is just the point-to-point communications opposite to those made in step 2.4 (step 4 of Alg. 2). The next test considered is the execution of the Algorithms 1 discarding only the MPI calls of steps 2.3, 2.4 and 5, respectively. Results are shown in column 4 of Table 1. The cost of the three communication episodes eliminated is minimal (3%) up to 64 CPU-cores, but represents 64% and 89% of the partitioning time with 512 and 4096 CPU-cores, respectively. It is thus clear that the principal scalability bottleneck are these three MPI communications that have been omitted in this last test.

Optimizations

The first optimization that we have considered is hiding communications under computations in Algorithm 2. In step 2.3 there is a MPI Alltoall communication used to get on each parallel process the size of the buffers that it will receive from the others. The size of each buffer is proportional to the number of weighted local bins to be sent to the corresponding process. In step 2.1 are counted the number of elements of the initial sublist contained in each coarse bin. This is an upper bound for the number of bins which are finally sent, since more than one element can be contained in the same local bin. Therefore, it is not possible to know the exact size of the buffers until they are generated (step 2.2). However the upper bound evaluated in step 2.1 is enough to start the communications because, according to MPI specifications, the non-blocking MPI Irecv instruction can be started using as "count" parameter the "maximum number of elements in receive buffer". In conclusion, communications can be started before the buffering process is carried out.

The corresponding optimized algorithm is shown in Algorithm 3. The MPI Alltoall communication has been advanced just after step 3.1 -the #elements in each coarse bin are exchanged instead of the send buffers size. Subsequently the MPI Irecv instructions are started (step 3.3). Finally, the non-blocking instructions MPI Isend (step 3.4.d) are interleaved with the weight evaluation and buffering (steps 3.4.a-3.4.c). On the other hand, in the step 5 of Algorithm 1, such an overlapping is not possible because only are performed a point-to-point communications without any substantial computation.

Algorithm 3 Optimized step 2 of Algorithm 1: overlapping 1. Loop around the sublist of mesh elements to count #elements to be sent to other proc. 2. MPI Alltoall to generate the communication scheme 3. Point-to-point: for each other process sending a buffer to me MPI Irecv 4. For each proc. q to which messages have to be sent:

(a) Allocate buffer of size N l g (local grid of proc. q) (b) Loop around the sublist to accumulate weight on the buffer (c) Store in the "send to proc. q" buffer bins with non-zero weight (d) Point-to-point: send the buffer MPI Isend 5. Point-to-point: MPI Waitall

Accumulate weight in local grid

The latency of the partitioning algorithm with this new version of step 2 is shown in column 3 of Table 2. With 8 and 64 CPU-cores the new algorithm is slightly slower than the initial version. As previously shown in Table 1, the communication costs represent a minimal part of the execution opt. 1 opt. 2 81,31 1,37 (104%) 1,37 (104%) 64 0,31 0,33 (106%) 0,33 (106%) 512 0,18 0,078 (43%) 0,071 (39%) 4096 0,20 0,28 (140%) 0,037 (18%) in these two cases (3%), therefore, there is not much to overlap. In fact, as the results demonstrate, it is better to perform all the communications together instead of interleaving them with calculations. In the execution with 512 CPU-cores, the benefit is noticeable, the partitioning time is reduced by 57%. Finally, in the extreme case with 4096 the new algorithm performs much worse than the initial one (40%). In this case the communications represent 89% of the partitioning time (see Table 1). Thus, like for the first two cases, the ratio between communications and computations is very unbalanced and there is not much to overlap. However, in this case the dominant part are the communications instead of computations. Again it is better to perform computations and communicators in two consecutive blocks instead of interleaving them. Results in terms of strong speedup are shown in Figure 8 (left). A second issue observed regarding the parallel performance is the load imbalance. In concordance with the previous analysis, we are focusing on the step 2 of Algorithm 1 and, in particular, on the communication operations performed in it. The MPI Alltoall operation does not produce any imbalance because has the same cost for all the parallel processes involved in it. However, the number of calls to the MPI Irecv and MPI Isend functions can be different among parallel processes and, eventually, produce an imbalance.

In particular, the number of calls to the MPI Isend function performed by one parallel process, depends on the number of mesh elements composing its initial sublist, and on how these are distributed in the domain. The more the collisions in the same coarse bin the less the number of messages to be sent. The more the collisions within bins of the local sub-grids the more are reduced the buffer sizes. In our tests, we have distributed the mesh elements into sublists of the same size, using the default ordering coming from the mesh generator.

On the other hand, the number of calls to the function MPI Irecv by each parallel process, depends on the number of mesh elements contained in its coarse bin, and on how these are distributed among the initial sublists. This is the most critical part in terms of load imbalance. The coarse bins enclosing zones where the mesh is more refined, contain more mesh elements and receive more messages from other parallel processes. In our example, the mesh of the respiratory system is very anisotropic and is specially dense on the nasal cavities and throat. Consequently the processes accounting for these zones receive much more messages than the others. This produces a significant imbalance. In fact, there are some parallel processes having an associated coarse bin without mesh elements -corresponding to an empty zone of the bounding box -and others receiving messages from more than 75% of the rest of parallel processes.

In order to reduce the imbalance generated in the data transmissions, we have designed a multilevel strategy. If P is the number of partitioning processes (or coarse bins), with P = 2 dp . We define a new grid, referred as vast grid, with Q = 2 dq elements, where q < p. The P partitioning processes can be divided into P/Q = 2 d(p-q) subsets with Q processes each. The ideas is to redistribute the elements of the initial sublists within each subset of Q processes, grouping the elements contained in the same bin of the vast grid. With this redistribution, the number of calls to the functions MPI Isend and MPI Irecv is bounded by P/Q, one per subset of Q processes. Moreover, this initial redistribution benefits the collision of element into bins of the coarse and local grids. This optimization can be seen as a modification to the input data of Algorithm 3. Q is chosen as a relatively low number and the redistribution is performed with the MPI Alltoallv function. The new version of the step 2 of Algorithm 1 is presented in Algorithm 4. To keep the coherence with this new version of the step 2, in the step 5 of Algorithm 1 are carried out first the point-to-point communications opposite to the MPI Irecv/Isend instructions of steps 4.4 and 4.5.d. And then the one opposite to the MPI Alltoallv of step 4.1. This multilevel approach is illustrated in Figure 7.

The result of this optimization is shown in column 4 of Table 2. We have fixed Q = 128, in the first two cases the redistribution is not carried out since the number of partitioning processes is lower than Q. With 512 CPU-cores the partitioning time is reduced by 61% with respect to the initial version Algorithm 4 Optimized step 2 of Algorithm 1: overlapping and sublists redistribution 1. if(P ≥ Q) First level sublists redistribution: MPI Allotoallv 2. Loop around the sublist of mesh elements to count #elements to be sent to other proc. 3. MPI Alltoall to generate the communication scheme 4. Point-to-point: for each other process sending a buffer to me MPI Irecv 5. For each proc. q to which messages have to be sent:

(a) Allocate buffer of size N l g (local grid of proc. q) (b) Loop around the sublist to accumulate weight on the buffer (c) Store in the "send to proc. q" buffer bins with non-zero weight (d) Point-to-point: send the buffer MPI Isend and, with 4096 by 82%. This final algorithm is capable of partitioning the 17.7M mesh of the respiratory system in less than 4 cents of second using 4096 CPU-cores. Results in terms of strong speedup are shown in Figure 8 (left). Finally, in Figure 8 (right) are compared the strong scalability of the initial and final versions of the algorithm for the second example, the mesh of the swirling combustor with 28.9M elements. The time measurements are shown in Table 3. The final version of the algorithms is more than 10× faster than the initial version using 4096 CPU-cores, in particular the mesh is partitioned in only 5 cents of second.

Partition quality

The present work is focused on the efficiency of the parallel partition algorithm. However, given that the ultimate goal of mesh partitioning is to enable the parallel execution of simulation codes, we compare the partition quality obtained with the presented SFC algorithm with the one obtained with Metis. To this end the Metis k-way algorithm with its default parameters was selected, which is very common in practice. In Figure 9 we compare the average time-step duration with different partitions for the two test cases under consideration: the respiratory system (left) and the swirling combustor (right). In each test case a different set of equations is solved, but in both of them the time-step is composed of an explicit part, in which the linear system is assembled, and an implicit part in which a linear solver is applied.

Results for both partitions and test cases are rather similar. Metis aims to generate balanced partitions minimizing also communication costs, while the SFC approach does not consider communication costs explicitly, but can be very precise on the load balancing. In general, in a system like the Blue Waters one, a load of ∼ 35K mesh elements per CPU-core would be considered optimal for a simulation with Alya. Therefore, the most representative executions of this test are the ones with 512 and 1024 CPU-cores. In the res- have been presented in detail. The final version is fully scalable and does not suffer any memory or computational bottlenecks, being able to partition a mesh of 28.9M cells in 5 cents of second using 4096 CPU-cores.

Moreover, our method maintains or even improves the quality of the partition at increasing the number of partitioning processes. We have compared it with Metis, both methods generate partitions for which the simulations performance is similar, differences obtained for different representative test cases are below 10%. Even when a much more comprehensive study is still necessary on this aspect, the SFC method seems to be competitive approach for parallel mesh partitioning with a clear potential for dynamic load balancing.

Figure 1 :

 1 Figure 1: Hilbert SFC recursive generation, 1D (top) and 3D (bottom)

Figure 2 :

 2 Figure 2: Hilbert SFC parallel generation. a) Initial coarse grid definition, b) orientation of each parallel process to continue the recursion, c) recursive generation of SFC in each parallel process.

Figure 3 :

 3 Figure3: Sequential algorithm for mesh partitioning using the Hilbert SFC. a) The mesh is embedded in a cartesian grid, b) a weight is assigned to each element of the cartesian grid ("bins") depending on the number of elements contained in it, c) the cartesian grid is traversed in the SFC order, accumulating the weights of the bins and assigning subset numbers depending on the objective weights of each subset.

Figure 4 :

 4 Figure 4: Parallel mesh partitioning using SFC.

Figure 5 :

 5 Figure 5: Geometric discretization of the respiratory system (left). Streamlines colored by velocity magnitude and iso-surface of Q-criterion pointing by arrows in the nasal cavities and throat(right).

Figure 6 :

 6 Figure 6: Geometric discretization of a swirling combustor (left). Instantaneous iso-surface of stoichiometric mixture fraction colored by pressure (right).

Figure 7 :

 7 Figure 7: Illustration of step 2 of Algorithm 1: initial version (top) optimized version with redistribution of sublist elements (bottom).

Figure 8 :

 8 Figure8: Strong speedup of the initial version of the partitioning algorithm and two optimizations: overlapping communications and computations (opt1), overlapping and initial redistribution (opt2). Left: mesh of the respiratory system (17.7M). Right: mesh of the swirling combustor (28.9M).

Figure 10 :

 10 Figure 10: SFC-partition of the respiratory system mesh (17.7M) elements.

  2020) under grant agreement 653838. J.C. Cajas acknowledges the financial support of the 'Consejo Nacional de Ciencia y Tecnología (CONACyT, México)' grant number 231588 290790. Ricard Borrell and Daniel Mira acknowledge the Juan de la Cierva postdoctoral grants with codes IJCI-2014-21034 and IJCI-2015-26686, respectively.

Table 1 :

 1 Partitioning time (sec.) with different numbers of CPU-cores.

	For the initial

Table 2 :

 2 Partitioning time (sec.) with different optimizations of the partitioning algorithm. Opt. 1: overlapping communications and computations. Opt 2: overlapping + initial redistribution. Respiratory system mesh (17.7M).

Table 3 :

 3 Partitioning time (sec.) with the initial version of the partitioning algorithm and its final optimization, for the swirling combustor mesh (28.9M).

Acknowledgments

The research leading to these results has received funding from the European Union Horizon 2020 Programme (2014-2020) and the Brazilian Ministry of Science, Technology and Innovation through Rede Nacional de Pesquisa (RNP) under the HPC4E Project (grant agreement No. 689772). This work is also part of the PRAC "Simulations of Aircraft Engines" supported by the National Science Foundation. It has also been financially supported by the PRACE preparatory access projects funded in part by the EU Horizon 2020 research and innovation programme (2014-

piratory system the SFC partition outperforms the Metis partition by 10% and 3%, respectively; contrarily, for the combustor the Metis partition outperforms SFC partition by 4% and 9%, respectively. The graph obtained from the discretization of the irregular respiratory system geometry makes difficult the balancing for Metis. When the load per CPU is very low communication costs become critical and graph partitions are a better option. Note that the SFC partitions have been generated in parallel, setting P (the number of partitioning processes) as the highest power of 8 lower then M (the number of subdomains generated), in other words, using the largest possible subset of the parallel processes engaged on the simulation for the mesh partitioning. On the other hand, Metis partitions have been generated sequentially but some studies show a degradation of the partition quality with its parallel version [START_REF] Shang | Performance analysis of large scale parallel cfd computing based on codesaturne[END_REF].

Finally, Figure 10 illustrates the partition of the respiratory system mesh. The partition number assigned to each non-zero weighted bin of the fine grid is plotted, this represents the output of step 4 of Algorithm 1.

Concluding Remarks

In this paper, a new SFC-based parallel mesh partitioner has been developed and implemented. The analysis and optimization of the algorithm


Hideo Aochi
email: h.aochi@brgm.fr

BRGM André Burnol

Mechanism of the M L 4.0 25th April 2016 earthquake in southwest of France in the vicinity of the Lacq gas field

The source mechanism of the ML4.0 25 th April 2016 Lacq earthquake (Aquitaine Basin, South-West France) is analyzed from the available public data and discussed with respect to the geometry of the nearby Lacq gas field. It is one of the biggest earthquakes in the area in the past few decades of gas extraction and the biggest after the end of gas exploitation in 2013. The routinely obtained location shows its hypocenter position inside the gas reservoir. We first analyze its focal mechanism through regional broad-band seismograms recorded in a radius of about 50 km epicentral distances and obtain EW running normal faulting above the reservoir. While the solution is stable using regional data only, we observe a large discrepancy between the recorded data on nearby station URDF and the forward modeling up to 1Hz. We then look for the best epicenter position through performing wave propagation simulations and constraining the potential source area by the Peak Ground Velocity (PGV).

The resulting epicentral position is a few to several km away to the north or south direction with respect to station URDF such that the simulated particle motions are consistent with the observation.

The initial motion of the seismograms shows that the epicenter position in the north from URDF is preferable, indicating the north-east of the Lacq reservoir. This study is an application of full waveform simulations and characterization of near-field ground motion in terms of an engineering factor such as PGV. The finally obtained solution gives a moment magnitude of Mw3.9 and the best focal depth of 4 km, which corresponds to the crust above the reservoir rather than its interior. This position is consistent with the tendency of Coulomb stress change due to a compaction at 5 km depth in the crust.

Therefore, this earthquake can be interpreted as a relaxation of the shallow crust due to a deeper gas reservoir compaction so that the occurrence of similar events cannot be excluded in the near future.

It would be necessary to continue monitoring such local induced seismicity in order to better understanding the reservoir/overburden behavior and better assessing the local seismic hazard even after the end of gas exploitation.

Introduction

There were 11 earthquakes of magnitude larger than 3.8 in metropolitan France in 2016 according to the available catalogue EMSC (European-Mediterranean Seismological Centre). Among them, the ML4.0 25 th April 2016 earthquake is suspected to be induced in gas exploitation field of Lacq, southwestern France, close to the Pyreneans mountains (Figure 1). Gas extraction in Lacq started in the late 1950's and ended in 2013. During the exploitation, the induced seismicity and crustal deformation has been intensively studied in the 1990's (e.g. [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Segall | Induced stresses due to fluid extraction from axisymmetric reservoirs[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF], and more recently summarized in [START_REF] Bardainne | Etude de la sismicité de Lacq et analyse des formes d'ondes par décomposition en chirplets[END_REF] and [START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. Geophysical measurement of crustal deformation (subsidence) was explained by mechanical model of reservoir compaction (e.g. [START_REF] Segall | Induced stresses due to fluid extraction from axisymmetric reservoirs[END_REF]. More than 2000 earthquakes have been recorded from 1974 to 1997 (ML < 4.2) and the spatial patterns of the seismicity were correlated to the reservoir geometry and preexisting fault structure, principally orienting WNW-ESE direction [START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. The ML4.0 2016 earthquake is quite large comparing to the past seismicity in the area (Figure 1B), although the national catalogue summarized by [START_REF] Cara | SI-Hex: a new catalogue of instrumental seismicity for metropolitan France[END_REF] may have removed recognized induced earthquakes from it. The seismicity shows that there are less large-magnitude events with time, while the increase in small events since mid-1990s probably reflects the improvement in the event detectability according to the network extension. This 2016 event was large enough to be perceived by local habitants and should be a scientific indicator in order to examine the current states of the crust and the reservoir of this old gas extraction field.

The efforts of the seismological observation have focused on the Pyrenean mountain range, about 30 km south from Lacq (Figure 1A), where French national broad-band and acceleration networks are deployed and the data are available from ORFEUS/EIDA (European Integrated Data Achieves). The ML4.0 2016 earthquake is localized in the national catalogue of Rénass (Réseau National de Surveillance Sismique) at (0.59°W, 43.44°N) and at 5km depth. One broad-band station URDF (Urdes, 0.5931°W, 43.4383°N) is located very closely to this epicenter only by a few hundred meters. Figure 2 shows examples of the available seismograms from the data archives. Note that all the stations other than URDF are at distances farther than 50 km. This earthquake is important to study, because the routinely obtained earthquake location corresponds to the inside of the reservoir of gas extraction field. In order to discuss possible relation to the reservoir and understand the underlying mechanism of this earthquake, we need to assess its source parameters. Therefore, we carry out the seismological analysis using available data. First, we will perform moment tensor inversion using very low frequencies from the regional data. The existence of the closest station URDF allows us to investigate the source parameters more precisely. Secondly, we will run ground motion simulations based on the obtained source mechanism. We will show that spatial distribution of Peak Ground Velocity (PGV) and particle motion patterns help constraining the possible earthquake location. This is the main purpose of this study; how to extract the useful information from one near-field ground motion data. Finally, for the obtained source parameters, we discuss the mechanical consistency through a calculation of Coulomb stress change due to a reservoir compaction.

Regional moment tensor inversion

Analysis procedure

The earthquake's location and magnitude are automatically obtained in the catalogue by RéNaSS.

For our first purpose of obtaining the focal mechanism and a reliable focal depth, we carry out fullwaveform moment-tensor inversions using the regional data.

We fix the epicenter position obtained in the Rénass catalogue, namely at (0.59°W, 43.44°N). We let the focal depth vary from 2 to 8 km, compared to the 5 km depth provided in the Rénass catalogue.

We use three stations at about 50 km distance from the epicenter and apply a band-pass filter between 16-32 seconds. This choice of source-receiver distances and frequency ranges are consistent with the operational moment-tensor inversion in other regions (e.g. [START_REF] Dreger | Determination of source parameters at regional distances with single station or sparse network data[END_REF][START_REF] Kubo | NIED seismic moment tensor catalogue for regional earthquakes around Japan: quality test and application[END_REF] and guarantees stable convergence. The Green's functions are calculated for each component of moment tensor in a regional 1D structure model [START_REF] Souriau | A tomographic study of the lithosphere beneath the Pyrenees from local and teleseismic data[END_REF][START_REF] Douglas | The importance of crustal structure in explaining the observed uncertainties in ground motion estimation[END_REF], as presented in Figure 3. For simplicity, we use the angular fault parametrization (strike, dip and rake angles) instead of independent six components of the moment tensor. The model parameters to be obtained are moment magnitude (Mw), origin time (t0), fault strike, fault dip and fault rake for a given focal depth.

We apply a Genetic Algorithm [START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF][START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF] (GAUL library) to find the best combination of the model parameters. GA is suitable for multi-dimensional global search problems where the search space potentially contains multiple local minimum. GA shows a rapid convergence only after a few iterations of generation among several independent model parameters in earthquake source analyses (e.g. [START_REF] Ulrich | Rapidness and robustness of finite source inversion from elliptical patches method using continuous GPS and acceleration data: 2011 Mw9.0 Tohoku earthquake[END_REF]. As we use a very low frequency band, we extend the times series by padding with zero so as to use the time window from 20 seconds before the origin time where syn and obs are the synthetic and observed signals of three components, respectively. The model parameters are searched such that misfit function is minimum.

Result -Source mechanism

The adopted GA algorithm shows a good convergence. Figure 4 shows the obtained focal mechanism for each focal depth from 2 km to 8 km. Comparing to the routinely obtained focal depth of 5 km, we obtain the minimum residual at 3 km depth. The residual is very stable from 2.5 km to 4 km depths. This relatively large depth range can be explained by a typical kilometric rupture length of a magnitude 4 earthquake. The result clearly shows that the event should be shallower than the one obtained by the routine location by one or two kilometers. The focal mechanism (strike, dip, rake)

shows a EW-running normal fault, being (78.1°, 45.2°, -100.2°) for 3 km depth location. On the other hand, the estimated moment magnitude (Mw 3.8) is close to the routinely determined local magnitude of 4.0.

Figures 5 and6 show the convergence of the GA inversion and waveform comparison for the minimum residual (focal depth of 3 km). It is found that the solution converges very quickly after about 10 generations. The GA sometimes attempts mutation process to maintain diversity (in mutation, a solution may change completely from the previous solution), but a better solution remains very stable around the final values. This reassures the reliability of the employed inversion process in this study.

One of the important purposes of this study is scrutinizing the ground motion in the near field where the populations felt this earthquake as reported in local media. Figure 7 compares the observations and synthetics at URDF in the same framework of the inversions, namely comparative to Figure 6. We show them for different focal depths. It is remarked that the obtained solution does not fit the observation at all, neither in its amplitude nor in phase. This discrepancy is because of the inadequate wave propagation in the near field probably due to the incorrect epicenter location, while the focal mechanism inversion is not bothered by this fact as the very low frequency band of the seismograms at distance (about 50 km) is only analyzed.

Near-field ground motion simulations and source parameters

Simulation Method

For our second purpose of assessing the ground motion observed in the very near field (URDF), we carry out 3D numerical simulations of seismic wave propagation. We apply a 3D finite difference method (Aochi et al., 2013a,b). One can easily calculate the ground motion at one station using semianalytical methods such as discrete wavenumber method [START_REF] Bouchon | A simple method to calculate Green's functions for elastic layered media[END_REF]. However, as described later, we aim to explorer the ground motion pattern with respect to the source position. For calculating the ground motions at many receiver positions, the finite difference method is more convenient than discrete wave number method, as it is sufficient to run a single simulation once.

We limit our model dimension to a 15 km (EW) x 15 km (NS) x 12.5 km (depth) volume. The location of URDF is less than 1 km distance from the epicenter of Rénass. We discretize the volume with a ∆s = 50 m spacing to have sufficient grids between source and receiver. The time step of finite difference simulation is taken as ∆t = 0.003 s, to satisfy the condition of stability and resolution for our interest.

This discretization in finite difference scheme is fine enough to analyze the frequency range lower than 1 Hz. The code is equipped with a Perfectly Matching Layer absorbing condition for minimizing artificial reflection from the outer boundaries [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]. Earthquake source is introduced as a point source with a source time function of a cubic B-spline function of 0.5-second duration (e.g. the same form of function given in [START_REF] Douglas | The importance of crustal structure in explaining the observed uncertainties in ground motion estimation[END_REF]. This function is suitable for simulating an earthquake of about magnitude 4 to characterize the low-frequency content of the signals (Aochi et al., 2013a). From engineering point of view of quantitative seismic hazard assessment, we calculate Peak Ground Velocity (PGV) from the calculated ground motion at each receiver, every 200 m around the hypocenter and every 500 m elsewhere after a low-pass filter of 1 Hz. The use of synthetic PGV is more reliable than using any empirical attenuation relation (e.g. Battaglia and Aki, 2003) or Ground Motion Prediction Equations (GMPEs), which provides a measure of amplitude of signals (acceleration, velocity, spectra and so on) in function of magnitude and source-receiver distance without taking into account of complex seismic wave radiation in different azimuths.

Result -source location

Let us check the near-field ground motion at URDF in more detail. We ran a forward modeling supposing a detailed 1D structure model after [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF] and [START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF].

The gas reservoir is in a 3D complex structure (low velocity zones at about 5 km and 7-8 km depths), which is not available digitally. However a local 1D structure model should be enough as epicentral distance to URDF is an order of 1 km and focal depth should be shallower than low velocity zones of the reservoir. The assumed 1D model is shown in Figure 8. The comparison of the ground motion for the best focal mechanism solution, namely Mw3.8 at 3 km depth, is shown in Figure 9. The duration of the given source time function is supposed 0.5 seconds for a typical magnitude 4 event. The waveforms are filtered between 0.05 and 1 Hz so as to verify the onset of the ground motion. Indeed, the simulated ground motions do not fit very well the observation. The phases are shifted and the relative amplitudes of component are different. We also plot the particle motions of the first five seconds in panel (C). The simulation shows a large movement to the WSW-ENE direction corresponding to the azimuth toward the source position, while the observation indicates the N-S direction. This test indicates that the direction where the wave comes from is different between the observation and the simulation and the source position may be incorrect.

Next we explore the ground motion patterns around the epicenter. We ran a simulation around the epicenter for many receivers (named internally from #0001 to #2041). Figure 10 shows the Peak Ground Velocity (the maximum value of ground velocity) for each component after a filter up to 1 Hz. This indicates a very specific ground motion pattern due to the wave radiation from the source position. In order to estimate the epicenter position, we select the area only for PGV(NS)/PGV(EW) > 3 and 0.043 cm/s < PGV(NS) < 0.387 cm/s (a factor of 3 with respect to the observed value), noting that from Figure 9, PGV(NS)/PGV(EW) = 4.64 and PGV(NS) = 0.129 cm/s for the observation. This constrains the potential area of the epicenter location, which is probably to the North or South direction by a distance of 3 upto 10 km with respect to the station URDF. We calculate the misfit using the same definition as Equation ( 1), and this is shown in Figure 11. For the local minimum positions of the misfit (#1249 in the North, and #1196 in the South), we compare the waveforms in Figure 12. The horizontal particle motions are plotted for the first five seconds of the ground motions filtered between 0.02 and 1 Hz. Both positions show the very similar particle motions dominating the North-South movements. Although the frequency band of simulations is still limited, we compare the beginning part of the NS ground motion without filtering to check the first motion direction (Figure 12C). The simulations for both receiver locations show an opposite sign of the ground motion. The observations are noisy and contain higher frequencies, but the displacement waveforms integrated once from the velocity show that the observation is rather consistent with the receiver #1196, namely ground moves first to negative (southern) direction.

Revisit of moment tensor inversion (final solution)

Our preferred location of this earthquake is different from the one routinely obtained by a few kilometers. Although the regional moment tensor inversion at very low frequencies at distances of about 50 km (Chapter 2) should not be influenced significantly due to a small difference in epincenter, it is worth re-analyzing the solution for the new location. We carry out the same procedure as Chapter 2 and Figure 13 illustrates the solutions found for different focal depths. The best solution is found for

Mw3.9 at 4 km depth (c.f. Mw3.8 at 3 km in Figure 4). The focal mechanism (strike, dip, rake) = (66.3°, 32.0°, -114.6°) is also consistent with the previous solution (78.1°, 45.2°, -100.2°) (Chapter 2). This revisit justifies our analysis procedure and the obtained solution.

Sensibility test

Here we test the sensibility of the source description on the near-field ground motions.

Magnitude 4 event briefly has a fault length L of 1 km, according to the empirical relation 𝐿𝐿 = 10 -2.01+0.50 𝑀𝑀 [START_REF] Wells | New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF]. As there are still a few or several kilometers from source to receiver, a point source assumption is adopted in the previous calculation. However, any detail in source description might affect the ground motions, and then the spatial pattern of PGV. Here, we carry out ground motion simulations adopting stochastically simulated finite earthquake scenarios among the database of Ide and Aochi (2005). Figure 14 shows two tested earthquake scenarios. Ide and Aochi (2005) consider a heterogeneous fault model, stochastically generated through Fractal feature, and let rupture to propagate spontaneously according to the stress and the friction. The initiated rupture spontaneously progresses and stops. Among more than 15000 scenarios ranging from magnitudes M1.3 to M6, the selected two scenarios have moment magnitudes of 3.78 (Model A) and 3.83 (Model B), respectively. Two earthquakes have different rupture directivities (ruptured area with respect to the epicenter) and the duration of rupture is about 0.5 seconds.

Figure 15 shows the simulated PGV map for the two models. We assume the final focal mechanism, namely (66.3°, 32.0°, -114.6°), by a difference of about 15° in each angle comparing to the previous point source simulation (Figures 10 and 11). The spatial pattern changes particularly around the epicenter in the EW component due to the extension of the fault plane. However the NS component is still dominant in the ground motions. When we apply the same restriction on the PGV map as Figure 10, namely PGV(NS)/PGV(EW) > 3 and 0.043 cm/s < PGV(NS) < 0.387 cm/s, we get almost the same area only by a slight difference of a few grids (< 2 km). This test confirms the utilization of point source approximation in our approach and also shows the stability regardless of the uncertainty of source mechanism within about 15°. The point of interest is not located above the fault plane, but is at a horizontal distance of more than 5 km. Thus the detail of the rupture process may not influence significantly on the ground motions, as the distance is several times larger than the fault dimension.

Discussion

We have studied a moderate earthquake precisely using regional networks. It is always better if local network was employed and made available; however, we have demonstrated that useful information can be extracted even from the limited data. In particular, the near-field ground motion recorded very closely to the epicenter allows us to calibrate the epicenter position through the PGV spatial distribution and particle motion. For the area of our study, there have been some other moderate earthquakes, such as the 2nd September 2013 (ML 4.0) and the 20th February 2017 (ML3.5)

after Rénass report (see also Figure 1B). We have applied the same procedure of the analysis as in the previous chapter on these two earthquakes, using the same ATE, TERF and PYLO stations, though TERF did not record the 2013 event. The obtained focal mechanisms are summarized in Table 1. The moment magnitudes are estimated as Mw3.25 and 3.21, respectively, much smaller than the 2016 event. For the Mw3.25 2013/09/02 event, we show the detailed analysis in Appendix, as this earthquake could be important in terms of reported local magnitude (ML4.0) as large as the 2016 event. Near-field ground motion simulation upto 1 Hz confirms the small magnitude of Mw3.25. The localization is preferably estimated in the south-west by a few kilometers of station URDF, and this may correspond to the top of the reservoir.

Seismic network is costly. Another alternative is to collect the perception of habitants. Hough (2014) analyzed the felt data, corrected in USGS "Did you feel it?" system, for the recent moderate earthquakes (Mw 4.0 -5.6) occurring in the central and eastern US and obtained a regression equation to correlate the felt perception and the ground motions. For smaller magnitude events, for example, during a local seismicity crisis in a post-mining site in the South-Eastern France, [START_REF] Dominique | Microseismicity in post mining context -The coal basin of Gardanne (Fuveau, Bouches du Rhone, Provence[END_REF] calibrated the seismicity localization using felt reports to cope for the absence of instruments. His implication was confirmed later by using a local monitoring network. In our study, we used a factor of PGV (Peak Ground Velocity) for the estimate of epicenter position. PGV scale is often related to the felt perception. Bommer et al. (2006) proposed a notion such as "just perceptible" (0.1 cm/s), "clearly perceptible" (0.65 cm/s) and "disturbing" (1.3 cm/s). The waveforms recorded at URDF peaked at 0.6 cm/s (Figure 2), exceeding the 0.1 cm/s threshold even after band-pass filtering between 0.05 and 1

Hz (Figure 9). This level is enough to call the perception of habitants. Furthermore, the simulation produces PGVs at around 0.4 cm/s in the frequency range between 0.05 and 1 Hz (Figure 10),

suggesting that PGV at full frequency ranges might have been large enough to reach locally the "disturbing" level of Bommer et al. (2006). Such perception is important in regions of moderate seismicity and should be quantified as a future task.

According to the available geological map, the reservoir size of Lacq is about 10 km in EW and 5 km in NS direction [START_REF] Segall | Induced stresses due to fluid extraction from axisymmetric reservoirs[END_REF][START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. Our epicenter location is rather around the northern bound of the reservoir (Figure 16). We are then interested in the stress field, which may have led to this earthquake. For this purpose, we carry out a simple elastic calculation without accounting for any complex behavior from the reservoir. Let us assume a vertical compaction of 1 cm along an elliptical horizontal interface (10 km x 6 km) at 5 km depth and then calculate the displacement and stress field at 4 km depth. This setting can be calculated in a homogeneous, semi-infinite elastic medium [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF][START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF]. Figure 17 shows the calculation at 4 km depth, for which our solution has the minimum residual. Stress field is calculated for both possible solutions of the focal mechanism, (strike, dip, rake) = (66.3°, 32.0°, -114.6°) and (274.7°, 61.2°, -75.4°). The Coulomb stress change ∆CCF is defined as

Δ𝐶𝐶𝐶𝐶𝐶𝐶 = Δ𝜏𝜏 -𝜇𝜇 × Δ𝜎𝜎 𝑛𝑛
where Δ𝜏𝜏 and Δ𝜎𝜎 𝑛𝑛 (positive is taken for compression) are the shear and normal stress changes, respectively (e.g. [START_REF] King | Static stress changes and the triggering of earthquakes[END_REF]. Although the vertical displacement is simply symmetric (Figure 17A), the stress field has a variation from negative to positive values (Figure 17B-D). The amplitudes of the source parameters (dimension and compaction) are indicative, thus so is the obtained stress level. We do not aim to discuss the absolute change of stress, but this simple calculation allows discussing the spatial pattern of ∆CCF due to the reservoir compaction and the earthquake localization we obtained.

Both areas identified in the previous chapter corresponds roughly to the southern edge of the reservoir

(5 km in the east and 3 km in the south from the elliptical reservoir) or in the northern part (4.5 km in the east and 9 km in the north). Stress change reduces very quickly with distance, such that the former position is influenced much more due to the reservoir deformation. Both positions could receive a positive change of ∆CFF. As both are located in the south and the north, respectively, the favorable fault solution is different; the south-dipping solution for the southern location and the north-dipping one for the north. This is consistent with the seismo-tectonic interpretation on the cross sections of [START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF].

Conclusion

The source mechanism of the ML4.0 25 th April 2016 Lacq, South-West France, earthquake is analyzed. The routinely obtained position shows its hypocenter position within the old gas reservoir.

We first analyze its focal mechanism through regional seismograms recorded in a 50 km radius around the epicenter and obtain EW running normal faulting. The solution is stable enough from the regional data. However a forward modeling of ground motion to the nearby station URDF does not explain the observation. We then look for the best epicenter position by performing wave propagation simulations.

We constrain the potential source area using spatial distribution of Peak Ground Velocity (PGV). This indicates that the epicenter position is preferable a few to several km away to the north or south direction with respect to station URDF. The particle motions become also consistent with the observation. From the first movement of the ground motion, the epicenter position in the north from URDF is preferable, indicating the north-eastern edge of the Lacq reservoir. This is more consistent with the induced seismicity migration to the north and the northeast of the Lacq gas field by [START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. It is shown here that a single nearby station data allows estimating better the source characteristics according to the prior information (e.g. routinely obtained earthquake location) if using the full waveforms, and an engineering factor such as PGV also facilitates the calibration. The final solution shows a moment magnitude of Mw3.9 and the focal depth of 4 km corresponding in the crust above the reservoir. This position is consistent with the Coulomb stress change pattern due to a compaction at 5 km depth in the curst. Therefore, this earthquake can be interpreted as a relaxation of the shallow crust due to a deeper reservoir compaction. The insight from our seismological and mechanical analyses will be useful for improving local and regional seismic hazard analysis. This paper shows the ability of the seismological study using regional network, but for further studies, any local monitoring network would be required even after the definitive closure of the gas field. catalogue [START_REF] Cara | SI-Hex: a new catalogue of instrumental seismicity for metropolitan France[END_REF]. (B) Seismicity in function of time for an area of longitude between 0.8°W and 0.3°W and of latitude between 43.25°N and 43.75°N from the same catalogue (magnitude in Mw).

Figure captions

The known biggest earthquake of ML4.2 1978/09/18 event is plotted by cross [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF]. The obtained source parameter is summarized in Table 1(a).

(

  4h44:10 UTC, by Rénass) for duration of 60 seconds. The misfit function (residual) is defined in

Figure 1 :

 1 Figure 1: (A) Map around the 2016/04/25 ML4.0 Lacq earthquake in the Pyrenees region in the South-Western France. The epicenter position is marked by star. The stations available from ORFEUS/EIDA are shown by triangles. The seismicity is plotted for the period 1970-2009 from SI-Hex / BCSF-LDG

  Three recent felt earthquakes of 2013/09/02, 2016/04/25 and 2017/02/20 are plotted (magnitude in ML). The period after 2010 is shaded since the SI-Hex/BCSF-LDG catalogue ends by 2009.

Figure 2 :

 2 Figure 2: East-west components of ground motions in velocity. The signals are aligned at 4:43:40 UT, and normalized by its maximum amplitude (in cm/s). Station locations are shown in Figure 1.

Figure 3 :

 3 Figure 3: 1D profile of seismic wave velocities used in this moment tensor inversions of this study.

Figure 4 :

 4 Figure 4: The obtained solution of focal mechanism for each focal depth. Residual is normalized by the minimum residual (focal depth = 3 km). The number shows moment magnitude obtained.

Figure 5 :

 5 Figure 5: GA convergence for the focal depth of z = 3 km. Population number starts from 100 at the first generation and increases up to 186 at the 21 st generation. Top panel shows the convergence of residual during the whole seeds. Bottom panels show the convergence of five model parameters with the generation. All the seeds are plotted by grey small points, and the best solutions of each generation are emphasized with large black symbols.

Figure 6 :

 6 Figure 6: Waveform fitting between the observation and synthetics at three stations for the focal depth of z = 3 km. Synthetics represent the final solution after the 21 st generation and the thin grey lines show the intermediate solutions (the best solution at each generation, shown in Figure5). The seismograms in velocity are filtered between 16-32 s. The origin time is then 30 + t0 (s). The residuals are calculated between 10 and 70 seconds.

Figure 7 :

 7 Figure 7: Comparison of waveforms at URDF for depths of z = 2, 3 and 8 km. This station is not used in the inversions. The calculation and illustrations are in the same way as Figure 6.

Figure 8 :

 8 Figure 8: 1D structure model used for near-field ground motion simulations.

Figure 9 :

 9 Figure 9: Comparison of ground motion at URDF between the observation and the simulation for the focal mechanism obtained at 3 km depth. (A) Positions of the source (Rénass) and receiver URDF. (B) Three component ground motions filtered between 0.05 and 1 Hz. (C) Horizontal particle motions for the time 0 to 5 seconds.

Figure 10 :

 10 Figure10: PGV (Peakg Ground Velocity) of each component of the ground motions simulated for the focal mechanism obtained at 3 km depth. The PGV is calculated after a filter upto 1 Hz. In the below panels, the area is restricted for PGV(NS)/PGV(EW) > 3 and 0.043 cm/s < PGV(NS) < 0.387 cm/s. The source is put at the origin of the figure. Small triangles show the synthetic receiver positions.

Figure 11 :

 11 Figure 11: Spatial distribution of misfits between the observation and synthetics for the area selected in Figure 10. The three components and the NS component only are shown. The misfit is calculated during the 10 seconds after a filter between 0.02 to 1 Hz (see text).

Figure 12 :

 12 Figure 12: The comparison of ground motions at receivers #1249 (6 km away to the North) and at #1196 (9 km away to the South). (A) The relative position of both receivers with respect to hypocenter. (B) Comparison of the waveforms in time series during 10 seconds (upper panel) and in particle motions during the first 5 seconds (bottom). A bandpass filtere between 0.05 and 1 Hz is applied. The waveforms are shifted on illustration. (C) Comparison of the waveforms in velocity and displacement for the NS component without filter.

Figure 13 :

 13 Figure 13: The finally obtained focal mechanism solution for each focal depth. Residual is normalized by the minimum residual (focal depth = 4 km). The number shows moment magnitude obtained.

Figure 14 :

 14 Figure14: Earthquake scenarios A and B on a finite fault, simulated byIde and Aochi (2005). The rupture time and final slip are shown for the two models. The slip direction is parallel to the horizontal axis. Moment magnitudes are 3.78 and 3.83, respectively. The star represents the initial point of earthquake rupture, and fixed to the hypocenter for the ground motion simulations.

Figure 15 :

 15 Figure 15: The PGV maps of three components calculated for the two finite source models A and B, shown in Figure 14. On the right, the area restricted by PGV(NS)/PGV(EW) > 3 and 0.043 cm/s < PGV(NS) < 0.387 cm/s is shown for the dominant NS component. See also caption of Figure 13.

Figure 16 :

 16 Figure 16: Obtained Localisation of the 2016 Lacq earthquake (ellipsoid), comparing to the routinely obtained solution (red circle) and station URDF (red rectangle). The dimension of the ellipsoid briefly corresponds to the area around the smallest residuals observed in Figures 4 and 15, namely +-1.5kmin EW, +-3.0 km in NS and +-1.5 km in depth. Map is adopted from Google Map (http://maps.google.fr). The mesh shows the top of the reservoir.

Figure 17 :

 17 Figure 17: Displacement and stress field at 4 km depth due to a vertical compaction of 1 cm, projected as broken elliptical area (10 km x 6 km), at 5 km depth. Stress is calculated for the both possible solution of the focal mechanisms, (strike, dip, rake) = (66.3°, 32.0°, -114.6°) and (274.7°, 61.2°, -75.42°), respectively, as fault plane 1 and 2. The two circles indicate briefly the two indicated areas of the earthquake. (A) Vertical displacement in cm, (B) shear stress in MPa to the rake direction, (C) normal stress in MPa and (D) ∆CFF in MPa to the rake direction.

Figure A1 :

 A1 Figure A1: Waveform (ground velocity in cm/s) fitting between the observation and synthetics at three stations for the 2013/09/02 earthquake with focal depth of z = 4 km. Synthetic (red curves) represent the final solution after the 21 st generation and the thin grey lines show the intermediate solutions. The seismograms in velocity are filtered between 16-32 s. The time 0 is taken at 12:36:35.

Figure A2 :

 A2 Figure A2: PGV (Peak Ground Velocity) of each component of the ground motions simulated for the focal mechanism obtained at 4 km depth. The PGV is calculated after a filter upto 1 Hz. The relative position of station URDF with respect to the supposed source position is shown. Small triangles show the synthetic receiver positions.

Figure A3 :

 A3 Figure A3: (a) Misfit between the observation and synthetics for the 10 seconds. The area is limited by PGV(NS)/PGV(EW) > 5 and 0.016 cm/s < PGV(NS) < 0.143 cm/s. Initially supposed URDF station position is shown by green. Two local minimums are shown in red and green. (b) Comparison of the waveforms filtered between 0.02 and 1 Hz. (c) Particle motion of the first 5 seconds.

  Figure 1 :

  Figure 7 :

  Figure 8 :

  Figure 9 :

Figure A3 :

  Figure A3 :

Table 1 :

 1 Obtained source parameters for (a) the ML4.0 2nd September 2013 and (b) the ML3.5 20th February 2017 earthquakes. Misfit is normalized for the minimum value of each earthquake. The

	reference time for the event is supposed at 12:36:35 and 17:03:41, respectively.
	(a)
	Focal depth z (km)	Mw	t0 (s)	Strike (°)	Dip (°)	Rake (°)	Misfit
	3.0	3.26	0.04	131.5	80.5	-34.7	1.31
	4.0	3.25	0.33	263.56	35.4	-90.3	1
	5.0	3.22	0.47	123.4	61.7	-51.5	1.11
	6.0	3.13	0.65	43.4	74.3	-176.6	1.35
	8.0	3.14	7.04	108.0	20.1	-66.9	2.29
	(b)
	Focal depth z (km)	Mw	t0 (s)	Strike (°)	Dip (°)	Rake (°)	Misfit
	3.0	3.30	0.70	105.7	81.0	67.6	1.03
	5.0	3.23	1.04	7.8	34.9	-169.6	1.00
	6.0	3.21	0.94	7.2	38.2	-170.9	1

Acknowledgments

This study has been funded by internal research & development fund from BRGM. The calculation of Green function was carried out at French national supercomputing center (GENCI/CINES) under the grants c2016-046700 (2016), A0010406700 (2017) and A0030406700 (2017-18). We thank two anonymous reviewers for their comments, which improved our manuscript.

Data

The seismological data (RESIF: Réseau sismologique et géodésique français; http://www.resif.fr/ ) used in this study are available from ORFEUS-EIDA data portal (http://www.orfeus-eu.org/data/eida/index.html) together with earthquake information provided by

Appendix

This appendix shows the detailed analysis for the 2013/09/02 ML4.0 earthquake mentioned in the discussion section. We apply the same procedure as in the text for the 2016/04/25 ML4.0 earthquake. Figure A1 shows the waveform fit for the focal mechanism inversion, supposing 4.0 km focal depth. There are only two stations ATE and PYLO available at epicenteral distance of about 50 km. The station URDF (a few kilometers from the supposed epicenter) is not used in the inversion.

Figure A2 shows the PGV map from the nearfield ground motion using the obtained focal mechanism. From the observation, the PGV in EW, NS and UD is 0.008, 0.048, and 0.019 cm/s, respectively. Figure A3 represents the grid search of better positions of source-receiver. The NS movement is dominant and this information limits the probable area. The PGV is much smaller than the 2016/04/25 earthquake and briefly consistent with our simulation of Mw3.25 source. The position #00142, namely the epicenter located in the south-west by a few kilometers, represents the minimum misfit and consistent with the polarity of the particle motions.