# Efficient Parallel Generation of Random Field of Mechanical Properties for Geophysical Application

Luciano de Carvalho Paludo, Victor Bouvier, Régis Cottereau, Didier Clouteau

► To cite this version:

Luciano de Carvalho Paludo, Victor Bouvier, Régis Cottereau, Didier Clouteau. Efficient Parallel Generation of Random Field of Mechanical Properties for Geophysical Application. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01515074

HAL Id: hal-01515074

https://hal.science/hal-01515074

Submitted on 27 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Efficient Parallel Generation of Random Field of Mechanical Properties for Geophysical Application

L. Paludo<fnref n="1" />, V. Bouvier<fnref n="1" />,
R. Cottereau<fnref n="1" />, D. Clouteau<fnref n="1" />

**Abstract** — We study the generation of random fields applied to problems where the domain is much larger than the characteristic distance of the properties fluctuation. Spectral Representation is used and a variation for isotropic fields is proposed. To improve scalability, a method to divide the domain in independent overlapping subsets is developed. The effectiveness of these modifications will be shown on tests in a 3D seismic waves propagation code running over hundreds of cores. Comparison with other schemes for large scale generation of random fields will also be performed.

**Keywords** — Random Fields, Gaussian Process, Spectral Representation, High Performance Computing.

## 1 Introduction

Modelling earth crust heterogeneities is necessary when we are interested in studying the seismic coda or higher frequency signals [1]. Its description requires an enormous amount of data. Fortunately when asymptotic regimes are considered, such as homogenization [2] or weak scattering regime [3] the solution of the mechanical problem depends only on some statistics of media properties. In this context a stochastic description of the media is adapted and the set of values for a given parameter can be represented by a random field realization.

An important requirement is that the sample generation cost (CPU time and memory) should remain small compared to the simulation time. We will consider in this work domains of dimension $L$ much larger than both the correlation length $\ell_c$ (or some characteristic size over which the fluctuations of the random field are significant) and the discretization step $h$. If $h > \ell_c$, the random field is a Gaussian white noise. We therefore restrict our attention to the case where $h < \ell_c \ll L$. In this context building a relevant field can become computationally expensive and choosing an efficient algorithm is essential.

We chose to simulate our fields using the classical algorithm by [4]. This method, based on spectral representation, can ensure the $C^\infty$ regularity of the field with little communication between processors, making it appealing to parallel implementation. The complexity of this method is $O(N^2)$, where $N$ is the size of the mesh. We developed a variation for isotropic materials, discretizing the wave number space partly uniformly and partly randomly. This approach decreases the complexity of the algorithm to $O(N^{1+1/d})$, where $d$ is the number of dimensions of space. Although spectral methods allow an easy discretization in space we are still summing contributions for every point over the whole domain. Statistically one can expect that, when the distance is larger than a given length, neglecting this residual contribution would affect in a minor way the generated field. This idea makes possible to divide the domain in several independent subsets. To handle the transition between the subsets an overlap is imposed. A method to make a statistically correct superposition is shown.

<footnote n="1">¹ Ecole Centrale Paris, Laboratoire MSSMat, UMR 8579 CNRS, grande voie des vignes, 92295 Chatenay-Malabry, France,
[luciano.de-carvalhol,regis.cottereau,didier.clouteau]@ecp.fr
[victor.bouvier]@student.ecp.fr</footnote>
## 2 Wave propagation in heterogeneous media

The equation describing elastic wave propagation in elastic (non-dissipative) media can be expressed as:

$$ \rho(\mathbf{x}) \frac{\partial^2 v}{\partial t^2}(\mathbf{x}, t) - \nabla_{\mathbf{x}} \{C(\mathbf{x}) : \nabla_{\mathbf{x}} \otimes v(\mathbf{x}, t)\} = 0, \quad (\mathbf{x}, t) \in \Omega \times \mathbb{R} \quad (1) $$

where $\rho(\mathbf{x})$ is the medium density, $C(\mathbf{x})$ is the fourth-order elastic tensor, and $v(\mathbf{x}, t)$ is the displacement field. To take into account media heterogeneity we compute $C = \{C(\mathbf{x}) : \mathbf{x} \in \Omega\}$ replacing the properties by a random field carrying inside its stochastic description.

## 3 Stochastic field generation

The generation algorithm should meet the following specifications: (i) the time of generation should remain small compared to the simulation time and (ii) each generated sample should represent well the required statistics. We only consider here the sampling of Gaussian fields because they are the basic building block of a large number of numerical schemes. The first-order marginal density can be modified locally by combining a direct and inverse Rosenblatt transforms [5], although one has to pay attention to the influence on the correlation function of the resulting random field [6, 7].

A common approach to sample a random field $\{u(\mathbf{x}) : \mathbf{x} \in \Omega \subset \mathbb{R}^d\}$ with a given correlation function $R$ is to search it as a linear combination of independent and identically distributed random variables, where $\Omega \subset \mathbb{R}^d$ is the medium and $d$ is the number of dimensions of space. For instance, the spectral representation is a classic way to sample gaussian random field [8]:

$$ u(\mathbf{x}) = \int_{k \in \Omega} \hat{R}^{1/2}(\mathbf{k}) \exp(i\mathbf{k} \cdot \mathbf{x}) dW(\mathbf{k}) \quad (\mathbf{x} \in \Omega) \quad (2) $$

where $\{W(\mathbf{k}) : \mathbf{k} \in \Omega\}$ is a Brownian motion, $\hat{R}$ is the Fourier transform of $R$ and $\mathbf{k} \cdot \mathbf{x}$ the inner product between $\mathbf{k}$ and $\mathbf{x}$.

There are several methods in the literature to compute the stochastic integral (2) and we will introduce its main ideas. The spectral method by M. Shinozuka and G. Deodatis [8], proposes the quadrature (3):

$$ u^{\text{S.D.}}(\mathbf{x}) = \sum_{n=0}^{N} \hat{R}^{1/2}(\mathbf{k}_n) \exp(i\mathbf{k}_n \cdot \mathbf{x}) \sqrt{\Delta_n} \xi(n) \quad (\mathbf{x} \in \Omega) \quad (3) $$

where $\xi = \{\xi(n) : n \le N\}$ is a white noise, $\mathbf{k}_n \in \Omega_n$ for all $n \le N$, $(\Omega_n)_{0 \le n \le N}$ is a partition of $\Omega$ and $\Delta_n$ is the Lebesgue measure of $\Omega_n$. The expression (3) is obtained thanks to the approximation $\int_{\Omega_n} f(\mathbf{k})d\mathbf{k} \approx f(\mathbf{k}_n)$. This representation ensures the $C^\infty$ regularity on $\Omega$ of the random field $u$ and decouples $u(\mathbf{x}_1)$ and $u(\mathbf{x}_2)$ for $\mathbf{x}_1 \ne \mathbf{x}_2$. This local formulation makes parallelism straightforward. Nevertheless, some conditions must be respected when using the Fourier transform in discrete spaces. Assuming $\Omega = [0, L]^d$ where $L$ is the domain scale, we define $\Delta x = L/M$ and $\Delta k$ respectively the discretization steps in space and wave number space. To avoid the field periodicity we must compute $\Delta k \le 2\pi/L$. On the other hand, it is useless to consider wave vectors where $||\mathbf{k}|| \ge k_{\text{max}} = 2\pi/\Delta x = N\Delta k$ because the grid is too coarse to represent it. We deduce 1:

$$ \Delta k = \frac{2\pi}{L}, \quad N = M \quad (4) $$

The number of points in the wave number space, $N$, depends on the spacial parameters $L$ and $\Delta x$. As a result, when generating $u$ over a large domain or a refined mesh, computational cost grows rapidly. A well-known manner to perform the sum (3) computation faster ($O(N \log N)$) is to use the Fast Fourier Transform (F.F.T.):

$$ u^{\text{FFT}}(p\Delta x) = \sum_{n=0}^{N} \hat{R}^{1/2}(k_n) \exp(2i\pi \frac{p \times n}{N}) \sqrt{\Delta_n} \xi(n) = \mathcal{F}^{-1}[(S_n)_{0 \le n \le N}]_p \quad (5) $$

where $S_n = \hat{R}^{1/2}(k_n) \sqrt{\Delta_n} \xi(n)$ and $p \in \{0, \dots, M\}$. Unlike the previous method, the FFT needs all the $\mathbf{x}' \in \Omega$ to compute $u(\mathbf{x})$, making parallelism difficult to implement. The $L^2$ norm of the error committed

<footnote n="1">¹The majoration of $N$ could be improved when we consider a cut-wave vector $\mathbf{k}_c$ ($\hat{R}(\mathbf{k}_c) \approx 0$) : $N = E[\frac{k_c \times L}{2\pi}] + 1$.</footnote>
on the correlation function is dominated by $M(L) \times N^{-1/d}$ where $M(L)$ increases with $L$.

Another classic way to compute (2) is to consider it as the expectancy of a random variable $\exp(i\mathbf{K} \cdot \mathbf{x})$ : $\int_{k \in \Omega} \hat{R}(k) \exp(i\mathbf{k} \cdot \mathbf{x})d\mathbf{k} = \mathbb{E}[\exp(i\mathbf{K} \cdot \mathbf{x})]$ where $K$ follows the probability density $\hat{R}(k)$. It is called the Randomization Method [9, 10] :

$$u^R(\mathbf{x}) = \frac{1}{\sqrt{N_r}} \sum_{n=0}^{N_r} \xi(n) \exp(i\mathbf{k}_n \cdot \mathbf{x}) \quad (\mathbf{x} \in \Omega) \quad (6)$$

where $(\mathbf{k}_n)_{n \le N_r}$ is a set of $N_r$ realisations of $\mathbf{K}$. The Randomization Method does not introduce aliasing or periodicity; there is no condition on $N_r$ which involves $\Delta x$ or $M$. It can be demonstrated that $\mathbb{L}^2$ norm of the error committed on the correlation function, thanks to the Central Limit Theorem, is dominated by $\nu(\Omega) \times N_r^{-1/2}$ where $\nu(\Omega) = \|\sqrt{1-R^2}\|_{\mathbb{L}^2(\Omega^2)}$. The fact that $\nu$ depends strongly on $\Omega$ is the main limitation to use this method for sampling over large domain.

## 4 Isotropic fields in spherical coordinates

When considering isotropic fields one can reduce the complexity of Spectral Representation method from $O(N^2)$ to $O(N^{1+1/d})$ using spherical coordinates to describe the vector $\mathbf{k}$. We choose randomly the two angles $\theta_n$ and $\phi_n$ that define the direction of $\mathbf{k}_n$ and his norm $r_n$ deterministically. The deterministic radius $(r_n)_{n \le N_r}$ assures that we explore all the spectrum and the random direction reduces the integral from a volume to a line with no further drawback. $(\mathbf{k}_n)_{n \le N_r}$ is defined as :

$$\mathbf{k}_n = r_n \{\cos(\theta_n) \sin(\phi_n), \sin(\theta_n) \sin(\phi_n), \cos(\phi_n)\}^T \quad (r_n \in \mathbb{R}^+) \quad (7)$$

where $\{\theta_n : n \le N_r\}$ and $\{\phi_n : n \le N_r\}$ are respectively white noises in $[0, 2\pi]$ and $[0, \pi]$.

$$u^{\text{I.S.}}(\mathbf{x}) = \sum_{n=0}^{N_r} \sqrt{\hat{R}(r_n) r_n \sin(\phi_n) \Delta_n} \exp(i\mathbf{k}_n \cdot \mathbf{x}) \xi(n) \quad (\mathbf{x} \in \Omega) \quad (8)$$

To ensure the fact that we explore all the spectrum, we are used to choose $N_r \approx N^{1/d}$ which is a good approximation of the number of point on the medium radius. The $\mathbb{L}^2$ norm of the correlation function error is dominated by $M^{\text{I.S.}}(L) \times N^{-1/d}$ where $M^{\text{I.S.}}(L)$ increases with $L$.

## 5 Localization of the sampling

As the domain becomes larger the computational cost of generating a sample grows rapidly and without threshold. To bound this overflowing instead of performing the whole domain at once, it could be interesting to sample over several smaller subdomains. The issue here is how to ensure regularity between the fields generated on different subdomains. We address to this problem making a transition overlapping area between subdomains.

Points separated by a distance larger than the correlation length are, by definition, uncorrelated. In the algorithms presented so far the mutual contribution of every point on the grid was considered. Now we subdivide the domain in smaller independent parts with a partition of unity $\psi = (\psi_i)_{i \in I}$ of $\Omega$ : $\sum_{i \in I} \psi_i(\mathbf{x}) = 1$ for all $\mathbf{x} \in \Omega$. We chose to write $u$ as :

$$u(\mathbf{x}) = \sum_{i \in I} u_{\Omega_i}(\mathbf{x}) \sqrt{\psi_i(\mathbf{x})} \quad (\mathbf{x} \in \Omega) \quad (9)$$

where $u_{\Omega_i}$ is a localized sample of $u$ over the subdomain $\Omega_i$ and for $i \ne j$, $\psi_i = 0$ over $\Omega_j$. The size of $\Omega_i$ and the overlap between $\Omega_i$ and $\Omega_j$ depends on the correlation length $\ell_c$ of the medium. An example is to decompose $\Omega = \cup_{i \in I} \Omega_i$ where each $\Omega_i$ is a sphere with $\alpha\ell_c$ for radius. The acceptable error $\varepsilon$ of the correlation function is chosen and the minimal size of intersection between the subsets can be calculated as :

$$\alpha\ell_c = M_{\psi}^{-1} \left( \frac{\varepsilon}{J \times \sup_{x \in \Omega_i} \{ \|x\|_2 R(x) \}} \right) \quad (10)$$

where $M_{\psi}(\alpha\ell_c) = \sup_{i \in I} \|\sqrt{\psi_i}\|_{\mathbb{L}^2(\Omega_i)}$. This method allows to exchange information locally, mitigating scalability issues.
# 6 Simulation

We use the generated samples to create a 3D elasticity tensor fields of a random isotropic material. For each property the statistical inputs are: first order marginal density, correlation model, correlation length, average and standard deviation. We need to communicate between the processors the seed for random number generation and the extreme coordinates (so we can find the global extremes). It ensures the regularity of the generated fields. Once each processor has this information it can perform the calculation independently. We can see one realization of a density field on Figure 1 used in a ballasted railway problem. The problem was performed in a 7,6 millions points mesh.

<figure><img src="image_1.png" /><figcaption>Figure 1: Random density field generated with isotropic spherical coordinates method (Left) and simulation snapshot of the displacement field in this media (Right)</figcaption></figure>

# References

[1] K. Aki and B. Chouet. Origin of coda waves: source, attenuation and scattering. *Journal of Geophysical Research*, 80:3322–3342, 1975.

[2] Y. Capdeville, L. Guillot, and J. J. Marigo. 1-D non periodic homogenization for the wave equation. *Geophys. J. Int.*, 181:897–910, 2010.

[3] Leonid Ryzhik, George Papanicolaou, and Joseph B. Keller. Transport equations for elastic and other waves in random media. *Wave Motion*, 24(4):327–370, December 1996.

[4] M. Shinozuka and G. Deodatis. Simulation of multi-dimensional Gaussian stochastic fields by spectral representation. *App. Mech. Rev.*, 49(1):29–53, 1996.

[5] M. Rosenblatt. Remarks on a multivariate transformation. 23(3):470–472, 1952.

[6] M. Grigoriu. Simulation of stationary non-gaussian translation processes. 124(2):121–126, 1998.

[7] B. Puig and J.-L. Akian. Non-gaussian simulation using Hermite polynomials expansion and maximum entropy principle. 19(4):293–305, 2004.

[8] M. Shinozuka and G. Deodatis. Simulation of stochastic processes by spectral representation. 44(4):191–205, 1991.

[9] P. R. Kramer, O. Kurbanmuradov, and K. Sabelfeld. Comparative analysis of multiscale gaussian random field simulation algorithms. 226:897–924, 2007.

[10] O. Kurbanmuradov, K. Sabelfeld, and P. R. Kramer. Randomized spectral and Fourier-wavelet methods for multidimensional Gaussian random vector fields. 245:218–234, 2013.