CN103530451A - Multi-grid Chebyshev parallel spectral element method with complex medium elastic wave propagation and simulation - Google Patents

Abstract

The invention discloses a multi-grid Chebyshev parallel spectrum element method. The method includes: a. Divide the calculation area into large and regular cells; b. Define the main grid and the auxiliary grid in the unit, the main grid is the Chebyshev configuration point; c. Use the main grid to divide the unit The wave field is truncated by Chebyshev, and the medium parameters and external forces in the element are truncated and expanded by using the auxiliary grid; d. Substitute the truncated expansions of the wave field, medium parameters, and external forces into the wave equation to obtain the element mass matrix and element stiffness matrix and the element external force vector; e. Based on the element mass matrix, element stiffness matrix and element external force vector of each element, the elastic wave equation is solved on the main grid.

Description

Multi-grid Chebyshev Parallel Spectral Element Method for Simulation of Elastic Wave Propagation in Complex Media

technical field

The invention relates to high-performance computing in the fields of geophysics, engineering seismology, computational acoustics, and the like.

Background technique

The simulation of elastic wave propagation in complex media plays an important role in the fields of geophysics, engineering seismology and computational acoustics. How to find a numerical simulation method that is more accurate, faster, and more suitable for parallelization has always been the goal and focus of researchers in related fields at home and abroad.

Numerical methods currently widely used in the simulation of elastic wave equations mainly include finite difference method, finite element method, pseudospectral method, and spectral element method. They have their own advantages and disadvantages and are applied in different fields.

The finite difference method approximates the derivative in the equation with a difference, thus transforming the partial differential equation into an algebraic equation for solution, and its error is determined by the number of discretized configuration points and the truncation error of the Taylor series. The finite difference method is widely used in computational geophysics and engineering earthquakes because of its simple programming and fast calculation speed. However, the low-order finite difference method requires a large number of grid points to ensure the accuracy, and the high-order finite difference method is difficult to deal with free interfaces and complex structures efficiently.

Based on the weak form of the wave equation, the finite element method divides the calculation area into a finite number of non-overlapping units. In each unit, the wave field is approximated by a low-order polynomial (such as a piecewise linear function). The finite element method is suitable for dealing with boundary and internal complex structures, and has been widely used in structural analysis and transient simulation. However, the finite element method is relatively seldom used in the field of large-scale elastic wave equation simulation. The reason is that the finite element method consumes a lot of computing resources, and the second is that the low-order finite element has dispersion phenomenon, while the usual high-order finite element has false wave.

The pseudospectral method uses infinitely differentiable global functions (such as Fourier series) to transform the problem from the physical space into the wavenumber domain for processing, and has good convergence. The pseudospectral method requires only a small number of spatial grid points to achieve high accuracy. In order to deal with boundary conditions, Chebyshev or Lagrange Legendre orthogonal polynomials are often used instead of Fourier series in space expansion. However, due to the uneven distribution of the configuration points of the orthogonal polynomial, it brings great restrictions to the choice of time step. At the same time, like the finite difference method, the pseudospectral method cannot easily deal with the problem of curved free interfaces and complex structures. In order to solve this kind of problem, some people have proposed methods such as curved coordinates and domain decomposition, but the cost is the increase of calculation amount.

The spectral element method was proposed by Patera in 1984, and it was mainly used in fluid mechanics in the early days. Seriani et al first introduced the Chebychev spectral element method into the numerical simulation of the acoustic wave equation in 1991, and researchers from various countries have done a lot of related research thereafter. It combines the finite element method and the pseudospectral method, and has both the flexibility of the finite element method for dealing with boundaries and structures and the fast convergence characteristics of the pseudospectral method. Under the premise of achieving the same accuracy, it can use more sparse element division than the traditional finite element, which reduces the calculation time and memory requirements.

The basic method of the spectral element method is to use high-order spectral expansion on each unit. The basis function represented by truncated orthogonal polynomial is selected, and interpolation is used on each unit to improve the convergence speed of the solution. Its main steps are: (1) first divide the calculated area into many sub-domains (units), each sub-domain is composed of several nodes (configuration points); (2) express the approximate solution as a truncated orthogonal Polynomial expansion; (3) Use Galerkin method to solve the variational scheme of the orthogonal problem, and obtain the global approximate solution.

The spectral element method can obtain higher precision with sparser grids and cells, but in the traditional spectral element method, there can only be a single homogeneous medium in each cell, which seriously reduces the computational efficiency in some cases. For example, when the structure of the medium is complex and the scale of change is small or even smaller than the wavelength, the units must be divided and solved according to the small scale of medium change, which will result in a huge calculation scale.

In addition, when the traditional spectral element method solves the wave equation, it needs to form a global stiffness matrix and a mass matrix. When considering attenuation, it also needs to introduce a damping matrix, which consumes a lot of storage space and limits the efficiency of its parallelization in the algorithm. Note that in the time iteration process of the wave equation, it is necessary to do multiple matrix and vector products, such as Sp _k , where S and p _k are the global stiffness (or mass) matrix and global vector, respectively. Considering the sparse nature of the global matrix, the global matrix and phasor multiplication consumes a lot of computing time.

The above two points have greatly hindered the large-scale practical application of spectral element method in geophysics, engineering seismology, computational acoustics and other fields.

Contents of the invention

The object of the present invention is to provide an algorithm capable of overcoming the above-mentioned disadvantages.

Therefore, the present invention provides a multi-grid Chebyshev parallel spectrum element method. The method includes:

a. Divide the computing area into large and regular units;

b. Define the main grid and auxiliary grid in the unit, the main grid is the Chebyshev configuration point;

c. Use the main grid to perform Chebyshev truncation and expansion of the wave field in the unit, and use the auxiliary grid to perform truncation and expansion of the medium parameters and external forces in the unit;

d. Substitute the truncated expansion of the wave field, medium parameters and external force into the wave equation to obtain the element mass matrix, element stiffness matrix and element external force vector;

e. Based on the element mass matrix, element stiffness matrix and element external force vector of each element, the elastic wave equation is solved on the main grid.

Preferably, the elastic wave equation is solved using a preconditioned conjugate gradient method.

Further preferably, the elastic wave equation includes the product of the global matrix and the global vector; the step of using the preconditioned conjugate gradient method to solve the elastic wave equation includes: calculating the element matrix, according to the corresponding relationship between the global node and the element node from the global vector Take out the unit vector; independently calculate the unit matrix-vector product on the unit, and then superimpose the resulting unit vector according to the corresponding relationship between the unit node and the global node to form a global result vector.

Further preferably, all units are evenly distributed to P CPU computing nodes, and the step of using the preconditioned conjugate gradient method to solve the elastic wave equation includes: the process of each computing node independently reads in the medium parameters of the area, and the computing unit Matrix; by superimposing and transferring the vector values on the common nodes of adjacent units between sub-regions, the part of the global vector allocated on the calculation node is obtained; calculating the unit matrix and vector product inside the calculation node to obtain the unit result vector; superposition and transfer Neighboring cells share the value of the result vector on the node, and obtain the part of the global result vector allocated on the node.

并取其对角线元素，作为单元预条件矩阵；各进程之间传递相邻单元共有节点对应的预条件矩阵元素值并叠加，用于更新单元预条件矩阵，从而形成全局预条件矩阵在各单元上的部分；各进程独立地初始化波场，开始时间迭代；各进程独立计算向量

各进程之间传递相邻单元共有节点对应的b_n元素值并叠加，用于更新b_n，独立地存储于各CPU节点的内存中；用预条件共轭梯度法结合逐元技术，迭代求解方程组

得到波场增量δu_n；各进程独立地用上一步求解得到的波场增量δu_n更新波场及其导数；判断时间迭代是否结束，若未结束则返回计算向量b_n的步骤继续迭代。Preferably, based on the element mass matrix, element stiffness matrix and element external force vector of each element, the step of solving the elastic wave equation on the main grid includes: independently calculating the matrix according to the element mass matrix and element stiffness matrix in each process

And take its diagonal elements as the unit preconditioning matrix; each process transmits and superimposes the element values of the preconditioning matrix corresponding to the common nodes of adjacent units to update the unit preconditioning matrix, thereby forming a global preconditioning matrix in each The part on the unit; each process independently initializes the wave field and starts time iteration; each process independently calculates the vector

Each process transmits and superimposes b _n element values corresponding to the common nodes of adjacent units, which are used to update b _n and are stored independently in the memory of each CPU node; use preconditioned conjugate gradient method combined with element-by-element technology to iteratively solve equation set

Obtain the wave field increment δu _n ; each process independently uses the wave field increment δu _n obtained from the previous step to update the wave field and its derivative; judge whether the time iteration is over, if not, return to the step of calculating the vector b _n to continue the iteration .

Further preferably, the step of using the preconditioned conjugate gradient method combined with the element-by-element technique includes: each process independently initializes the solution and calculates the residual vector; each process transmits the residual vector values corresponding to the common nodes of adjacent units and superimposes to obtain an updated residual vector ; Each process independently calculates the matrix-vector product, and uses a parallel algorithm to calculate the vector inner product; each process transmits and superimposes the vector values on the common nodes of adjacent units to obtain an updated matrix-vector product; updates the solution and the residual vector.

Preferably, each unit may contain multiple media inside.

Because the embodiment of the present invention adopts the multi-grid technology, the unit division is simple, and the division of complex small units in the traditional method is avoided. The most suitable auxiliary grid and interpolation basis function can be used in the unit to describe the change of the medium, and the wave field Solving on the sparse main grid reduces the calculation scale and saves computing resources. At the same time, by using element-by-element technology, the storage space and calculation amount can be reduced by thousands of times, and the parallel method can be used to achieve higher parallel efficiency.

Description of drawings

Figure 1 is a flow chart of the multi-grid Chebyshev parallel spectral element method;

Figure 2 is a flow chart of the subroutine of the preconditioned conjugate gradient method.

Detailed ways

The present invention will be described in detail, clearly and completely below in conjunction with the accompanying drawings and specific embodiments. Apparently, the described embodiments are only some of the embodiments of the present invention, but not all of them. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without creative efforts fall within the protection scope of the present invention.

For the traditional spectral element method, there can only be a single homogeneous medium in each unit. When performing elastic wave simulation in complex media with multi-scale structures, it is necessary to divide a large number of irregular small units, which consumes a lot of storage space and calculation time. The present invention proposes a multi-grid Chebyshev parallel spectral element method, which divides the calculation area into relatively large and regular non-overlapping units, uses fine grids to describe medium changes inside the units, and uses larger regular units for iterative calculations. Greatly reduces memory and computing resource requirements. Using Chebyshev spectral element method has the characteristics of stiffness matrix and mass matrix with analytical precision and the parallel advantages of multi-grid method in grid division, and introduces element-by-element technology to achieve high-efficiency parallelization, further reducing storage space and calculation It can realize efficient and high-precision elastic wave simulation, and serve large-scale high-performance computing in geophysics, engineering seismology, computational acoustics and other fields.

(1) Multi-grid Chebyshev spectral element method element subdivision and multi-grid division technical scheme

Firstly, the computing area is divided into relatively large and regular non-overlapping units. For any unit, there is a reversible conversion function between the local coordinates of the unit and the reference unit, and the conversion between the actual physical area and the reference unit is completed through this transformation. The calculation of all units can be carried out on the reference unit uniformly. The unit can choose regular shape, such as square unit can be used for two-dimensional, and cube unit can be used for three-dimensional. On the one hand, this method avoids the division of complex units in the traditional method and saves calculation time. Simple, it can reduce the amount of calculation when solving the unit matrix.

Subsequently, in order to describe the different distributions of the wave field, medium and external force in the unit, several sets of independent grids can be defined in the reference unit. Among them, the grid used to describe and discretize the sound field is called the main grid, which uses Gauss-Lobatto-Chebyshev (GLC) configuration points, and uses the Chebyshev-Lagrange interpolation basis function to truncate and expand the wave field in the unit (for example, two The wave field in the dimensional unit can be expanded as

为Chebyshev-Lagrange插值基函数）。in, is the wave field value at each configuration point in the unit,

interpolation basis functions for Chebyshev-Lagrange).

The grid used to describe the medium parameters is called the auxiliary grid. According to the distribution and change of the medium, the most suitable auxiliary grid type, auxiliary interpolation function and order can be selected to accurately describe the distribution of the medium in the unit, (for example, the medium parameters in a two-dimensional unit can be expanded as

${\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\underset{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}^{\mathrm{<span\; class="notranslate">\; 3</span>}}{\mathrm{<span\; class="notranslate">\; \&phi;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; \&phi;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

为单元内辅助网格配置点上的介质参数值，φ_l(ξ)为定义于辅助网格上的插值函数）。比如，若介质在单元内连续变化，可以采用与主网格类似的GLC点作为辅助网格节点，但其阶数L不必与主网格阶数N相同。而若介质在单元内发生不连续的跳变，介质参数用阶梯函数来描述，可以采用均匀分布的配置点作为辅助网格节点，以单位方波脉冲函数作为插值基函数，阶数根据介质变化情况选择，若介质变化尺度小则可取较高阶数。一般情况下，需计算的介质分布是不均匀不规则的，因此需要利用辅助网格将单元内的介质参数作截断展开。in

is the medium parameter value on the auxiliary grid configuration point in the cell, φ _l (ξ) is the interpolation function defined on the auxiliary grid). For example, if the medium changes continuously in the unit, GLC points similar to the main grid can be used as auxiliary grid nodes, but its order L does not have to be the same as the main grid order N. If the medium has a discontinuous jump in the unit, the medium parameters are described by a step function, and uniformly distributed configuration points can be used as auxiliary grid nodes, and the unit square wave pulse function is used as the interpolation basis function, and the order changes according to the medium The case selection, if the medium change scale is small, a higher order can be taken. In general, the distribution of the medium to be calculated is uneven and irregular, so it is necessary to use the auxiliary grid to truncate and expand the medium parameters in the unit.

When it is necessary to describe the spatial distribution of the external force within the unit, a set of auxiliary grids can also be used to truncate and expand the part of the internal and external forces of the unit that varies with space.

${\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; K</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\underset{{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; K</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; f</span>}}_{{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; \&psi;</span>}}_{{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; \&psi;</span>}}_{{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&xi;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

The acoustic wave equation is transformed into a weak form according to the conventional method of the spectral element method:

$\frac{{\mathrm{<span\; class="notranslate">\; d</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; dt</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; ,</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; \&rho;c</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; )</span>}_{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; a</span>}{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; )</span>}_{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; )</span>}_{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}<span\; class="notranslate">\; ,</span>$

in ${(w,\frac{1}{{\mathrm{\&rho;c}}^2}u)}_{\mathrm{\&Omega;}}={\&Integral;}_{\mathrm{\&Omega;}}w\frac{1}{{\mathrm{\&rho;c}}^2}\mathrm{ud\&Omega;},---\left(4\right)$

(w,f) _Ω = ∫ _Ω wfdΩ

Substituting the three expansions (1), (2), and (3) into the weak form equation, the second-order linear ordinary differential equations expressed in matrix form are obtained after derivation:

Mü(t)+Ku(t)=F(t) (5)

In practice, the element mass matrix M ^e , the element stiffness matrix K ^e and the element load vector F ^e are calculated, which are used to solve the subsequent equations. The auxiliary grid will no longer be used after this step of calculation is completed.

The matrices and vectors associated with the helper grid can free up memory. The order of these element matrices and vectors is determined by the main grid order.

It only needs to iteratively solve the equations on the main grid, and the main calculation scale only depends on the order of the main grid. When the medium structure is complex and the change scale is small, the medium change can be accurately described by increasing the order of the auxiliary grid, and the main grid is still relatively sparse. Therefore, compared with the traditional method, the calculation scale can be significantly reduced, and the memory and memory in the calculation can be saved. Purpose of CPU demand.

(2) The technical scheme of element-by-element technology applied in the Chebyshev spectral element method

The spectral element method based on Chebyshev expansion is solved iteratively by the preconditioned conjugate gradient method.

The second-order linear ordinary differential equations of spectral element method can be solved iteratively by time integration and preconditioned conjugate gradient algorithm. An example of the steps of a time integration algorithm is as follows:

initialization: $\mathrm{no}=0,u_1=u_0,{\stackrel{\&Center\; Dot;}{u}}_1={\stackrel{\&CenterDot;}{u}}_0,{\stackrel{\&CenterDot;\&CenterDot;}{u}}_1=m^{-1}(f_0-{\mathrm{Ku}}_0)$

Iterate to the next moment:

n=n+1

${\mathrm{<span\; class="notranslate">\; b</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&Delta;tM</span>}{\stackrel{<span\; class="notranslate">\; \&Center\; Dot;</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; +</span>\frac{{\mathrm{<span\; class="notranslate">\; \&Delta;t</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; K</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; )</span>$

Depend on $(m+\frac{{\mathrm{\&Delta;t}}^2}{4}K){\mathrm{\&delta;u}}_{\mathrm{no}}=b_{\mathrm{no}}$ Solve to get δu _n

Update the solution, looping iteratively until complete:

u _n+1 =u _n +δu _n

${\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}}_{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; \&Delta;t</span>}}{\mathrm{<span\; class="notranslate">\; \&delta;\; u</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; \&CenterDot;</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}}_{\mathrm{<span\; class="notranslate">\; no</span>}}$

将其写作Sx=b，此方程用预条件共轭梯度法求解的一种算法步骤示例为：In the above iterative solution process, it is necessary to repeatedly calculate

Written as Sx=b, an example algorithmic step for solving this equation using the preconditioned conjugate gradient method is:

Initialization: k=0,r ₀ =b ₀ -Sx ₀ ,p ₀ =z ₀ =Q ^-1 r ₀ , where x ₀ is the initial solution, and Q is the preconditioning matrix. In order to achieve a fast convergence effect, a diagonal preconditioning matrix can be used here, that is, the diagonal elements of the matrix S are selected, and other elements are set to zero as the preconditioning matrix.

Iterate, updating the solution and residual vectors:

$\mathrm{<span\; class="notranslate">\; \&alpha;</span>}<span\; class="notranslate">\; =</span>\frac{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; Sp</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&alpha;</span>}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; \&alpha;</span>}{\mathrm{<span\; class="notranslate">\; Sp</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}$

Convergence check: If ||r _k+1 ||<ε||b||, the iteration terminates,

If not converged, update the conjugate gradient search direction:

${\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; Q</span>}}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>\frac{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}$

In the process of solving, it is necessary to calculate the product Sp _k of the global matrix and the global vector, and since the global matrix of the spectral element method is aggregated from the unit matrix according to the corresponding relationship between the global node number and the unit node number, from the point of view of the unit , the resulting global vector of Sp _k can be obtained through a series of unit-level calculations, as follows:

${\mathrm{<span\; class="notranslate">\; Sp</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>\underset{\mathrm{<span\; class="notranslate">\; e</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; S</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; e</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}<span\; class="notranslate">\; (</span>{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; S</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; p</span>}}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; e</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; v</span>}}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 6</span>}<span\; class="notranslate">\; )</span>$

和

分别为构成全局矩阵和全局向量的单元矩阵和向量，而

代表聚合。in,

and

are unit matrices and vectors constituting the global matrix and global vector, respectively, and

Represents aggregation.

In the above formulation, the global matrix-vector product is performed on local units.

Specific steps are as follows:

通过全局节点编号和单元节点编号的对应关系，从全局向量p_k扩散得到；a. First, calculate the unit matrix, and the unit vector

Through the corresponding relationship between the global node number and the unit node number, it is obtained from the global vector p _k diffusion;

b. Then, calculate the element matrix-vector product on the local element to obtain the local element result vector ${\tilde{v}}_k^e={\tilde{S}}^e{\tilde{p}}_k^e;$

再通过全局节点编号和单元节点编号的对应关系聚合成全局结果向量。这样，不需要生成全局矩阵S就得到了所需要计算的全局向量。c. Finally, the local cell result vector

Then the global result vector is aggregated through the corresponding relationship between the global node number and the unit node number. In this way, the global vector to be calculated can be obtained without generating the global matrix S.

Therefore, it is not necessary to integrate and store the global matrix in the program, but only need to calculate and store the unit matrix. The cell matrix is dense, which greatly reduces memory requirements and improves computational efficiency.

(3) Parallelization scheme based on multi-grid and element-by-element technology

Due to the use of multi-grid technology, the shape of the spectral element method unit is regular, which is easy to block processing in parallel algorithms, and it is easy to obtain load balance for each parallel sub-region. At the same time, due to the use of element-by-element technology, in the Chebyshev-based spectral element method program, most of the work can be limited to local elements, and the entire calculation can be completed through aggregation and diffusion. Therefore, the spectral element program using multi-grid and element-by-element technology has excellent parallel efficiency in structure, and has a special basis for large-scale parallel computing.

It is assumed that elastic wave propagation simulation calculations are performed on a parallel cluster computer system with P computing nodes, and MPI (message passing interface) is used as a parallel computing development environment. After unit division, the entire calculation model has N _e spectral unit units, which are equally distributed to P CPU processes, and each process independently reads the medium parameters in its corresponding sub-area, and uses the multi-grid method to independently calculate the quality of the unit matrix and element stiffness matrix and are stored in the CPU node's memory. In this way, in the space discrete stage of the spectral element method, through the task allocation of CPU nodes, the purpose of memory balance and calculation balance is achieved, and during this process, no data exchange is required between nodes, and high-efficiency parallelism can be achieved.

The equations formed by the Chebyshev spectral element method are solved by the preconditioned conjugate gradient method, where the preconditioned matrix is the matrix formed by the diagonals of the global matrix. In actual operation, there is no need to form a global matrix, as long as the matrix element values on the common nodes of adjacent units between the sub-domains responsible for each process are superimposed, and then passed to each process to update the corresponding element values of the diagonal of the unit matrix. . This can avoid a large amount of data transfer and improve parallel efficiency.

A large number of calculations in elastic wave propagation simulations exist in iterative calculations in time discretization. At each step of the time iteration, the system of equations is iteratively solved using the preconditioned conjugate gradient method. Due to the use of element-by-element technology, the global matrix and vector product can be converted into unit matrix and vector product. The main calculation steps are as follows:

a. Obtain the part of the global vector allocated on the node by superimposing and transferring the vector value on the common node of the adjacent units between the sub-regions;

b. Calculate the unit matrix and vector product inside the computing node to obtain the unit result vector, where the vector is obtained in the previous step, and the unit matrix remains unchanged in the iterative calculation;

c. Superimpose and transfer the result vector values on the common nodes of adjacent units to obtain the part of the global result vector allocated on the nodes, which is used to update the solution and residual vector, and enter the next iteration.

Taking a two-dimensional elastic wave simulation as an example, the implementation of the multi-grid Chebyshev parallel spectral element method is described in detail below. Among them, the medium parameters have discontinuous jumps on a small spatial scale. It is assumed that P CPU computing nodes are used for elastic wave propagation simulation calculations on a parallel cluster computer system, and MPI (message passing interface) is used as a parallel computing development environment.

Figure 1 is a flowchart of the multi-grid Chebyshev parallel spectral element method. As shown in Figure 1, the specific implementation process has the following steps:

(1) The calculation area is divided into N _e ×N _e non-overlapping quadrilateral units. These units can be evenly distributed to P CPU computing nodes, and each process independently reads the media parameters in its corresponding sub-area, and determines the parameters such as the order of the main grid and the auxiliary grid, and the location of the configuration point.

The main grid order N is selected according to the calculation accuracy and time requirements, and the auxiliary grid order L and configuration point positions used to describe the medium are selected according to the distribution of the medium. The auxiliary grid order K and The configuration points are selected according to the distribution of external forces. N and L are used in the calculation of the mass matrix and stiffness matrix, as well as the read-in medium parameters, and N and K are required in the calculation of the external force vector.

(2) Each process independently uses the multi-grid method to calculate the element mass matrix, element stiffness matrix and element external force vector.

The calculation formula of 2D elastic wave simulation is as follows:

${\mathrm{<span\; class="notranslate">\; m</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{cc}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; m</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}& \mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}& {\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; m</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}\end{array}\right]<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; K</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{cc}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; K</span>}}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}& {\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; K</span>}}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}\\ {\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; K</span>}}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; \&CenterDot;</span>}& {\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; K</span>}}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}\end{array}\right]<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; f</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{c}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}\\ {\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}\end{array}\right]<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

in

${\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; m</span>}}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{\mathrm{<span\; class="notranslate">\; 4</span>}}\underset{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; \&rho;</span>}}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}^{\mathrm{<span\; class="notranslate">\; e</span>}}{\mathrm{<span\; class="notranslate">\; B</span>}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\mathrm{<span\; class="notranslate">\; B</span>}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}$

${<span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; ]</span>}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; =</span>\underset{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; [</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&mu;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; )</span>\frac{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{\stackrel{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{<span\; class="notranslate">\; \&OverBar;</span>}}{\mathrm{<span\; class="notranslate">\; B</span>}}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\mathrm{<span\; class="notranslate">\; B</span>}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&mu;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}\frac{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{\mathrm{<span\; class="notranslate">\; B</span>}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\stackrel{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{<span\; class="notranslate">\; \&OverBar;</span>}}{\mathrm{<span\; class="notranslate">\; B</span>}}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; ]</span>$

${<span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; ]</span>}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; =</span>\underset{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; B</span>}}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; B</span>}}}_{{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&mu;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; B</span>}}}_{{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; B</span>}}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; )</span>$

${<span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; ]</span>}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; =</span>\underset{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}<span\; class="notranslate">\; [</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&mu;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; )</span>\frac{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{\mathrm{<span\; class="notranslate">\; B</span>}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\stackrel{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{<span\; class="notranslate">\; \&OverBar;</span>}}{\mathrm{<span\; class="notranslate">\; B</span>}}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&mu;</span>}}_{{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}\frac{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{\stackrel{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{<span\; class="notranslate">\; \&OverBar;</span>}}{\mathrm{<span\; class="notranslate">\; B</span>}}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\mathrm{<span\; class="notranslate">\; B</span>}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; j</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; l</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; ]</span>$

${<span\; class="notranslate">\; [</span>{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; ]</span>}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}}{\mathrm{<span\; class="notranslate">\; 4</span>}}\underset{{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; K</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; A</span>}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}{\mathrm{<span\; class="notranslate">\; A</span>}}_{{\mathrm{<span\; class="notranslate">\; i</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}{<span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; ]</span>}_{{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span>$

as well as

Among them, element mass matrix M ^e , element stiffness matrix K ^e , and element external force vector F ^e . ρ is the density of the medium, λ, μ are the Lame coefficients of the medium, and the subscripts l ₁ l ₂ represent the coordinates of the auxiliary grid in the unit. The Chebyshev interpolation polynomial expanded by truncating the main grid, φ(ξ) is the interpolation polynomial expanded on the medium auxiliary grid, and ψ(ξ) is the interpolated polynomial expanded on the external force auxiliary grid.

并取其对角线元素，以向量方式存储，作为单元预条件矩阵。标号e表示单元上的局部矩阵或向量。对角预条件矩阵，这是预条件共轭梯度法求解中的一种预条件矩阵选取方法，可以达到快速收敛的效果，并且对角矩阵可以以向量形式存储，减少内存占用，提高计算效率。(3) Calculate the matrix independently in each process

And take its diagonal elements, store them as vectors, and use them as unit preconditioning matrices. The reference e denotes a local matrix or vector on the unit. Diagonal preconditioning matrix, which is a preconditioning matrix selection method in the preconditioning conjugate gradient method, can achieve fast convergence, and the diagonal matrix can be stored in the form of vectors, reducing memory usage and improving computing efficiency.

(4) The preconditioning matrix element values corresponding to the common nodes of adjacent units are transferred between each process and superimposed to update the unit preconditioning matrix, thereby forming the part of the global preconditioning matrix on each unit, which is independently stored in each CPU in the node's memory.

(5) Each process independently initializes the wave field and starts time iteration.

(6) Each process calculates the vector independently $b_{\mathrm{no}}=\mathrm{\&Delta;tM}{\stackrel{\&CenterDot;}{u}}_{\mathrm{no}}+\frac{{\mathrm{\&Delta;t}}^2}{2}(f_{\mathrm{no}}-K{u}_{\mathrm{no}}).$

(7) The b _n element values corresponding to the shared nodes of adjacent units are transferred among the processes and superimposed for updating b _n , which are independently stored in the memory of each CPU node.

(8) Use the preconditioned conjugate gradient method combined with the element-by-element technique to iteratively solve the equation system $(m+\frac{{\mathrm{\&Delta;t}}^2}{4}K){\mathrm{\&delta;u}}_{\mathrm{no}}=b_{\mathrm{no}}.$ Get the wave field increment δu _n .

(9) Each process independently updates the wave field and its derivatives with the wave field increment δu _n obtained from the previous step.

(10) Determine whether the time iteration is over, if not, return to step (6) to continue the iteration.

(11) Each process independently writes the result to a file for storage.

Figure 2 is a flow chart of the subroutine of the preconditioned conjugate gradient method. As shown in Figure 2, this process includes the following sub-steps:

a. Each process independently initializes the solution and calculates the residual vector.

b. Each process transmits the residual vector values corresponding to the common nodes of adjacent units and superimposes them to obtain an updated residual vector.

c. Each process independently calculates the matrix-vector product, and uses a parallel algorithm to calculate the vector inner product.

d. Each process transmits and superimposes the vector values on the common nodes of adjacent units to obtain an updated matrix-vector product.

e. Update the solution and residual vectors.

Judging whether the convergence condition is reached, if not, return to step c for loop iteration.

Because the embodiment of the present invention adopts the multi-grid technology, the unit division is simple, and the division of complex small units in the traditional method is avoided. The most suitable auxiliary grid and interpolation basis function can be used in the unit to describe the change of the medium, and the wave field Solving on the sparse main grid reduces the calculation scale and saves computing resources. At the same time, by using element-by-element technology, the storage space and calculation amount can be reduced by thousands of times, and the parallel method can be used to achieve higher parallel efficiency.

Take a two-dimensional square calculation area Ω with side length L as an example, in which the minimum wavelength of the propagating elastic wave is λ _min , and the medium has a non-uniform distribution with a scale of 1/10λ _min .

To calculate with the traditional spectral element method, the unit must be divided according to the non-uniform scale of the medium, and the calculation area Ω is divided into N _e ×N _e non-overlapping quadrilateral units, the unit side length is 1/10λ _min , and the unit Chebyshev expansion takes N order, the total number of computing nodes is (N _e N+1)×(N _e N+1), and the size of the global matrix is S _g =(2×(N _e N+1)×(N _e N+1 )) ² . The global matrix-vector product with the largest calculation load is used to measure the calculation requirements, and only focus on the relatively time-consuming multiplication operation. The total number of multiplication operations required in the traditional spectral element method is T _g =(2×(N _e N+ 1)×(N _e N+1)) ² .

由于采用逐元技术，无需形成全局矩阵，则存储单元矩阵的内存需求为However, with the multi-grid Chebyshev parallel spectral element method, larger elements can be used, usually the element size can be larger than the minimum wavelength λ _min of the propagating elastic wave. Here, if λ _min is taken as the unit side length, the calculation area is divided into (N _e /10)×(N _e /10) units, and the Chebyshev expansion on the unit is still taken as N order, then the size of the unit matrix is

Due to the use of element-by-element technology, there is no need to form a global matrix, and the memory requirement of the storage cell matrix is

S _e =(N _e /10)×(N _e /10)×(2×(N+1)×(N+1)) ² . In terms of computing requirements, the global matrix-vector multiplication is transformed into (N _e /10)×(N _e /10) unit matrix-vector multiplications by element-by-element technology, and the total number of multiplication operations required is T _e =( N _e /10)×(N _e /10)×(2×(N+1)×(N+1)) ² .

Use the matrix size to measure the memory requirements of the two methods:

$\frac{{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; g</span>}}}{{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}}<span\; class="notranslate">\; =</span>\frac{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span>{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; \&ap;</span>{\mathrm{<span\; class="notranslate">\; 100</span>}\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; (</span>\frac{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span>$

The calculation amount of the two methods is measured by the number of multiplication operations in the matrix-vector product:

$\frac{{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; g</span>}}}{{\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}}<span\; class="notranslate">\; =</span>\frac{{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span>{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; \&ap;</span>{\mathrm{<span\; class="notranslate">\; 100</span>}\mathrm{<span\; class="notranslate">\; N</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; (</span>\frac{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 11</span>}<span\; class="notranslate">\; )</span>$

If N _e =100, N=6, then That is to say, by adopting multi-grid and element-by-element techniques, both memory requirements and computing requirements are reduced by at least 500,000 times.

In addition, when using the MPI parallel development environment, the main factors determining the parallel efficiency are load balance and the amount of data transferred between CPU nodes, while the multi-grid Chebyshev parallel spectrum element method load balances, the amount of data to be transferred is small, and can reach More than 90% parallel efficiency.

The specific embodiments described above have further described the purpose, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above descriptions are only specific embodiments of the present invention and are not intended to limit the scope of the present invention. Protection scope, within the spirit and principles of the present invention, any modification, equivalent replacement, improvement, etc., shall be included in the protection scope of the present invention.

Claims (7)

1. the Chebyshev of grid more than the Chebyshev spectral element method that walks abreast, its feature comprises:

A. zoning is divided into unit large and rule;

B. in unit, define main grid and auxiliary grid, main grid is Chebyshev collocation point;

C. utilize main grid that wave field in unit is blocked to expansion as Chebyshev, utilize auxiliary grid that medium parameter in unit and outer masterpiece are blocked to expansion;

D. by wave field, medium parameter and external force block expansion substitution wave equation, obtain force vector outside element mass matrix, element stiffness matrix and unit;

E. force vector outside the element mass matrix based on each unit, element stiffness matrix and unit solves equations for elastic waves in main grid.

2. many grids Chebyshev as claimed in claim 1 spectral element method that walks abreast, wherein utilizes Conjugate Gradient Method With Preconditioning to solve equations for elastic waves.

3. many grids Chebyshev as claimed in claim 2 spectral element method that walks abreast, wherein equations for elastic waves comprises the product of overall matrix and Global Vector; The described step of utilizing Conjugate Gradient Method With Preconditioning to solve equations for elastic waves comprises: computing unit matrix, and the corresponding relation retrieval unit of pressing global node and cell node from Global Vector is vectorial; Independent computing unit matrix-vector product on unit, then the corresponding relation stack of the unit vector of result being pressed to cell node and global node, form global outcome vector.

4. many grids Chebyshev as claimed in claim 2 spectral element method that walks abreast, wherein, P CPU computing node all given in all unit, the described step of utilizing Conjugate Gradient Method With Preconditioning to solve equations for elastic waves comprises: the process of each computing node is independently read in the medium parameter in this region, computing unit matrix; By superposeing and transmitting the vector value on the total node of adjacent cells between subregion, obtain the part that Global Vector distributes on computing node; In computing node internal compute unit matrix and vector product, obtain unit result vector; Result vector value on stack and the total node of transmission adjacent cells, obtains the part that global outcome vector distributes on node.

5. many grids Chebyshev as claimed in claim 1 spectral element method that walks abreast, wherein, force vector outside the element mass matrix based on each unit, element stiffness matrix and unit, the step that solves equations for elastic waves in main grid comprises:

In each process according to element mass matrix, element stiffness matrix independence compute matrix

and get its diagonal entry, as the pre-conditional matrix in unit;

Between each process, transmit pre-conditional matrix element value stack that the total node of adjacent cells is corresponding, for the pre-conditional matrix of updating block, thereby form the part of the pre-conditional matrix of the overall situation on each unit;

Each process is initialization wave field independently, start time iteration;

Each process independence compute vector $b_n=\mathrm{\&Delta;tM}{\stackrel{\&CenterDot;}{u}}_n+\frac{{\mathrm{\&Delta;t}}^2}{2}(f_n-K{u}_n);$

Between each process, transmit total b corresponding to node of adjacent cells _nelement value stack, for upgrading b _n, be stored in independently in the internal memory of each cpu node;

With Conjugate Gradient Method With Preconditioning in conjunction with by first technology, iterative system of equations $(M+\frac{{\mathrm{\&Delta;t}}^2}{4}K){\mathrm{\&delta;u}}_n=b_n.$ Obtain wave field increment δ u _n;

Each process solves by previous step the wave field increment δ u obtaining independently _nupgrade wave field and derivative thereof;

Judge whether time iteration finishes, if do not finish, return to compute vector b _nstep continue iteration.

6. many grids Chebyshev as claimed in claim 5 spectral element method that walks abreast, wherein, comprises in conjunction with the step by first technology with Conjugate Gradient Method With Preconditioning:

Each process is independent initially to be dissolved, and calculates remaining vector;

Each process is transmitted remaining vector value the stack that the total node of adjacent cells is corresponding, obtains the remnants vector upgrading;

Each process independence compute matrix vector product, by parallel algorithm compute vector inner product;

Each process is transmitted vector value the stack on the total node of adjacent cells, obtains the matrix-vector product upgrading;

More new explanation and remaining vectorial.

7. many grids Chebyshev as claimed in claim 1 spectral element method that walks abreast, wherein, medium can be contained in each inside, unit.

Images