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 spectral element method. The method includes the steps that a, a calculation area is divided into large and regular units; b, main grids and auxiliary grids are defined in the units, and the main grids are Chebyshev collocation points; c, Wave fields in the units are subjected to Chebyshev truncation expansion by the utilization of the main grids, and medium parameters and external forces in the units are subjected to truncation expansion by the utilization of the auxiliary grids; d, wave field, medium parameter and external force truncation expansion formulas are substituted into a wave equation to obtain unit mass matrixes, unit stiffness matrixes and unit external force vectors; e, on the basis of the unit mass matrixes, the unit stiffness matrixes and the unit external force vectors, an elastic wave expansion is solved on the main grids.

Description

The parallel spectral element method of many grids Chebyshev of complex dielectrics elastic wave modeling

Technical field

The present invention relates to geophysics, engineering seismology, the high-performance calculation in the fields such as calculating acoustics.

Background technology

The propagation simulation of complex dielectrics Elastic Wave is in geophysics, and the fields such as engineering seismology and calculating acoustics have critical role.How to find a kind of more accurately, method for numerical simulation quicker, parallelization preferably is also state, inside and outside association area researcher's target and emphasis always.

The current numerical method of generally using in equations for elastic waves simulation, mainly contains method of finite difference, finite element method, pseudo-spectrometry and spectral element method etc.They have relative merits separately and are applied in different fields.

Method of finite difference, by the derivative difference approximation in equation, solves thereby partial differential equation is converted to algebraic equation, and its error is determined by the truncation error of discretize collocation point number and Taylor progression.Method of finite difference is because programming is simple and direct and computing velocity is comparatively fast widely used in calculating in geophysics and engineering seismology.But low order method of finite difference needs a large amount of net points to guarantee precision, and High-Order Finite-Difference Method is difficult to process efficiently free interface and labyrinth problem.

The weak form of Finite Element Method based on wave equation, is divided into zoning the unit of limited non-overlapping copies, and in each unit, wave field is similar to by low order polynomial expression (as piecewise linear function).Finite element method is applicable to processing border and inner labyrinth, aspect structure analysis and transient simulation, is obtaining a wide range of applications.But finite element method is relatively less in large-scale equations for elastic waves simulation field application, and its reason one is that finite element method need to consume a large amount of computational resources, the 2nd, there is dispersion phenomenon in low order finite element, and there is pseudo wave in common high-order limited unit.

Pseudo-spectrometry utilization infinitely can be micro-overall situation function (as Fourier Fourier progression), problem is converted into wavenumber domain from physical space and processes, there is good the holding back property of locating.Pseudo-spectrometry only needs space lattice seldom to count can reach very high precision.In order to process boundary condition, when space development, often by Chebyshev Chebyshev or Lagrangian Legendre orthogonal polynomial, replace Fourier progression.But because the collocation point distribution of orthogonal polynomial is very inhomogeneous, bring considerable restraint to the selection of time step.Meanwhile, as method of finite difference, pseudo-spectrometry can not be processed the problem of crooked free interface and labyrinth freely.For solving this class problem, have people to propose the methods such as crooked coordinate and Region Decomposition, but cost is the increase of calculated amount.

Spectral element method is proposed in 1984 by Patera, is mainly used in early days fluid mechanics.Seriani etc. were incorporated into Chebychev spectral element method ACOUSTIC WAVE EQUATION first numerical simulation Zhong， various countries researcher in 1991 has done a large amount of correlative studys thereafter.It combines finite element method and pseudo-spectrometry, has had finite element concurrently and has processed border and the dirigibility of structure and the Fast Convergent characteristic of pseudo-spectrometry.Reaching under identical precision prerequisite, can adopt the dividing elements more sparse than traditional finite element, reduced computing time and memory requirements.

The basic skills of spectral element method is with higher-order spectrum, to launch on each unit.The basis function that the orthogonal polynomial of choosing to block represents utilizes collocation point interpolation on unit, to improve the speed of convergence of solution.Its key step is: (1) is first divided into many subdomains (unit) the region of calculating, and each subdomain is comprised of some nodes (collocation point); (2) orthogonal polynomial expansion that approximate solution is expressed as blocking in each subdomain; (3) by Galerkin method, solve the variation form of quadrature problem, obtain overall approximate solution.

Spectral element method can obtain degree of precision with more sparse grid and unit, but in traditional spectral element method, in each unit, can only have single uniform dielectric, the serious counting yield that reduces in some situation.Such as complicated when dielectric structure, change when yardstick is very little is even less than wavelength, must and solve by the small scale division unit of media variations, will cause so great calculating scale.

In addition, during traditional spectral element method solving wave equations, need to form overall stiffness matrix and mass matrix, while considering decay, also need to introduce damping matrix, need to expend a large amount of memory spaces, on algorithm, limit the efficiency of its parallelization.Notice in the time iteration process of wave equation, need to do repeatedly matrix and vector product, as Sp _k, S and p here _kbe respectively overall rigidity (or quality) matrix and Global Vector.Consider the sparse characteristic of overall matrix, the matrix of the overall situation and phasor multiply each other and have expended a large amount of computing times.

Above 2 have all greatly been hindered spectral element method in geophysics, and engineering seismology is calculated the large-scale practical application in the fields such as acoustics.

Summary of the invention

The object of the present invention is to provide a kind of algorithm that can overcome above-mentioned shortcoming.

For this reason, the invention provides the parallel spectral element method of a kind of many grids Chebyshev.The method 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.

Preferably, utilize Conjugate Gradient Method With Preconditioning to solve equations for elastic waves.

Further preferably, 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.

Further preferably, 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.

Preferably, 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

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

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.

Further preferably, with Conjugate Gradient Method With Preconditioning, in conjunction with the step by first technology, comprise: 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.

Preferably, medium can be contained in each inside, unit.

The embodiment of the present invention is owing to adopting multi-grid technique, dividing elements is simple, avoided the division of complicated junior unit in classic method, in unit, can adopt most suitable auxiliary grid and Interpolation-Radix-Function to describe the variation of medium, and wave field solves in more sparse main grid, reduce calculating scale, save computational resource., utilize by first technology meanwhile, can reduce memory space and calculated amount in thousandfold ground, and utilize parallel method to reach higher parallel efficiency.

Accompanying drawing explanation

Fig. 1 is the parallel spectral element method process flow diagram of many grids Chebyshev;

Fig. 2 is Conjugate Gradient Method With Preconditioning subroutine flow chart.

Embodiment

Below in conjunction with the drawings and specific embodiments, the present invention is carried out to detailed, clear, complete explanation.Obviously, described embodiment is only the present invention's part embodiment, rather than whole embodiment.Embodiment based in the present invention, those of ordinary skills, not making all other embodiment that obtain under creative work prerequisite, belong to the scope of protection of the invention.

For having single uniform dielectric in each unit of traditional spectral element method, while thering is the Simulation of elastic waves in Multi-scale model complex dielectrics, need to divide a large amount of irregular junior units, expend a large amount of memory spaces and the defect of computing time, the present invention proposes the parallel spectral element method of a kind of many grids Chebyshev, zoning is divided into larger and regular not overlapped elements, in inside, unit, with fine grid blocks, media variations is described, and adopt larger regular unit to carry out iterative computation, greatly reduce the demand to internal memory and computational resource.Utilize Chebyshev spectral element method to have and resolve the stiffness matrix of precision and the feature of mass matrix and the many grid methods parallel advantage on grid is divided, introduce by first technology simultaneously and realize high efficiency parallelization, memory space and calculated amount have further been reduced, realize the Simulation of elastic waves of high-efficiency high-accuracy, for geophysics, engineering seismology, the extensive high performance computing service in the fields such as calculating acoustics.

(1) many grids Chebyshev spectral element method element subdivision and many grids partitioning technology scheme

First zoning is divided into larger and regular not overlapped elements.For there being a reversible transfer function between any unit ， unit local coordinate and reference unit, by this conversion, complete the conversion between actual physics region and reference unit.The calculating of all unit can be unified on reference unit and implement.Regular shape can be chosen in unit, as two dimension can be selected square shaped cells, three-dimensional available regular cube unit, the method is avoided the division of complex unit in classic method on the one hand, save computing time, the transfer function of actual cell and reference unit is also simpler on the other hand, while solving cell matrix, can reduce calculated amount.

Subsequently, in order to describe respectively wave field, medium and the external force different distributions in unit, can in reference unit, define several groups of grids independently mutually.Wherein, for describing with the grid of discretize sound field, be called main grid, it adopts Gauss-Lobatto-Chebyshev (GLC) collocation point, and (for example in two-dimentional unit, wave field is deployable is to utilize Chebyshev-Lagrange Interpolation-Radix-Function to block expansion to the wave field in unit

Wherein, for the wave field value on each collocation point in unit,

for Chebyshev-Lagrange Interpolation-Radix-Function).

And be called auxiliary grid for describing the grid of medium parameter.According to the distribution of medium and situation of change, can select most suitable auxiliary grid type, auxiliary interpolating function and exponent number accurately to describe the distribution of medium in unit, (for example in two-dimentional unit, medium parameter is deployable is

${\widetilde{\mathrm{\α}}}^e({\mathrm{\ξ}}_1,{\mathrm{\ξ}}_2,t)=\underset{l_1=0}{\overset{L}{\mathrm{\Σ}}}\underset{l_2=0}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\α}}_{l_1{l}_2}^3{\mathrm{\φ}}_{l_1}\left({\mathrm{\ξ}}_1\right){\mathrm{\φ}}_{l_2}\left({\mathrm{\ξ}}_2\right),---\left(2\right)$

Wherein

for the medium parameter value on auxiliary grid collocation point in unit, φ _l(ξ) for being defined in the interpolating function on auxiliary grid).Such as, if medium changes continuously in unit, can adopt with the similar GLC point of main grid as auxiliary grid node, but its exponent number L needn't be identical with main grid exponent number N.And if there is discontinuous saltus step in medium in unit, medium parameter is described with step function, can adopt equally distributed collocation point as auxiliary grid node, using unit square-wave pulse function as Interpolation-Radix-Function, exponent number is selected according to media variations situation, if the little desirable higher-order number of media variations yardstick.Generally, it is inhomogeneous irregular needing the dielectric distribution of calculating, and therefore need to utilize auxiliary grid that the medium parameter in unit is blocked to expansion.

When need to, when unit internal fine is described the space distribution of external force, also can, with one group of auxiliary grid, external force in unit similarly being blocked to expansion with the part of spatial variations.

${\tilde{f}}^e({\mathrm{\ξ}}_1,{\mathrm{\ξ}}_2,t)=\underset{k_1=0}{\overset{K}{\mathrm{\Σ}}}\underset{k_2=0}{\overset{K}{\mathrm{\Σ}}}{f}_{k_1{k}_2}^e\left(t\right){\mathrm{\ψ}}_{k_1}\left({\mathrm{\ξ}}_1\right){\mathrm{\ψ}}_{k_2}\left({\mathrm{\ξ}}_2\right).---\left(3\right)$

ACOUSTIC WAVE EQUATION is converted into weak form by the conventional method of spectral element method:

$\frac{d^2}{{\mathrm{dt}}^2}{(w,\frac{1}{{\mathrm{\ρc}}^2}u)}_{\mathrm{\Ω}}+a{(w,u)}_{\mathrm{\Ω}}={(w,f)}_{\mathrm{\Ω}},$

Wherein ${(w,\frac{1}{{\mathrm{\ρc}}^2}u)}_{\mathrm{\Ω}}={\∫}_{\mathrm{\Ω}}w\frac{1}{{\mathrm{\ρc}}^2}\mathrm{ud\Ω},---\left(4\right)$

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

By in three expansions (1), (2), (3) substitution weak form equation, through deriving, obtain the Second Order Linear Ordinary Differential Equation group with matrix representation:

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

In reality, calculate element mass matrix M ^e, element stiffness matrix K ^ewith unit loads vector F ^e, for follow-up solving equations.This step will no longer be used auxiliary grid after having calculated.

The matrix relevant to auxiliary grid and vector can releasing memories.These cell matrixs and vectorial exponent number are determined by main grid exponent number.

Only need be in main grid iterative system of equations, main calculating scale also only depends on the exponent number of main grid.When dielectric structure is complicated, change yardstick hour, can media variations accurately be described by increasing auxiliary grid exponent number, therefore it is comparatively sparse that main grid still keeps, and compares can significantly reduce calculating scale with classic method, reach save calculate in the object of internal memory and cpu demand.

(2) technical scheme of applying in Chebyshev spectral element method by first technology

The spectral element method Conjugate Gradient Method With Preconditioning iterative launching based on Chebyshev.

Second Order Linear Ordinary Differential Equation group pot life integration and the preconditioned conjugate gradient method iterative of spectral element method.A kind of step example of time integral algorithm is as follows:

Initialization: $n=0,u_1=u_0,{\stackrel{\·}{u}}_1={\stackrel{\·}{u}}_0,{\stackrel{\·\·}{u}}_1=M^{-1}(f_0-{\mathrm{Ku}}_0)$

A downward moment iteration:

n=n+1

$b_n=\mathrm{\ΔtM}{\stackrel{\·}{u}}_n+\frac{{\mathrm{\Δt}}^2}{2}(f_n-K{u}_n)$

By $(M+\frac{{\mathrm{\Δt}}^2}{4}K){\mathrm{\δu}}_n=b_n$ Solve and obtain δ u _n

More new explanation, loop iteration is until complete:

u _n+1=u _n+δu _n

${\stackrel{\·}{u}}_{n+1}=\frac{2}{\mathrm{\Δt}}{\mathrm{\δu}}_n-{\stackrel{\·}{u}}_n$

In above iterative process, need repeatedly calculate

write Sx=b, a kind of algorithm steps that this equation solves with Conjugate Gradient Method With Preconditioning is exemplified as:

Initialization: k=0, r ₀=b ₀-Sx ₀, p ₀=z ₀=Q ^-1r ₀, x wherein ₀for initial solution, Q is pre-conditional matrix.In order to reach the effect of Fast Convergent, can adopt the pre-conditional matrix in diagonal angle herein, choose the diagonal entry of matrix S, other elements get zero, as pre-conditional matrix.

Iteration, more new explanation and remaining vector:

$\mathrm{\α}=\frac{(r_k,z_k)}{(p_k,{\mathrm{Sp}}_k)},x_{k+1}=x_k+\mathrm{\α}{p}_k,r_{k+1}=r_k-\mathrm{\α}{\mathrm{Sp}}_k$

Convergence checks: if || r<sub TranNum="161">k+1</sub>||<ε || b|| iteration stops,

If do not restrain, upgrade the conjugate gradient direction of search:

$z_{k+1}=Q^{-1}{r}_{k+1},p_{k+1}=z_{k+1}+\frac{(r_{k+1},z_{k+1})}{(r_k,z_k)}{p}_k,k=k+1$

In solution procedure, need to calculate the product Sp of overall matrix and Global Vector _k, and because the overall matrix of spectral element method is that cell matrix is polymerized by global node numbering and the corresponding relation of cell node numbering, therefore from the viewpoint of unit, Sp _kresult Global Vector can be by the calculating of a series of unit level, as shown in the formula:

${\mathrm{Sp}}_k=\left(\underset{e}{\mathrm{\&Sigma;}}{\tilde{S}}^e\right)p_k=\underset{e}{\mathrm{\&Sigma;}}\left({\tilde{S}}^e{\tilde{p}}_k^e\right)=\underset{e}{\mathrm{\&Sigma;}}{\tilde{v}}_k^e=v_k---\left(6\right)$

Wherein,

with

be respectively the cell matrix and the vector that form overall matrix and Global Vector, and

represent polymerization.

In above-mentioned formula, the matrix-vector product of the overall situation completes in local unit.

Concrete steps are as follows:

A. first, computing unit matrix, and unit is vectorial

by the corresponding relation that global node is numbered and cell node is numbered, from Global Vector p _kdiffusion obtains;

B. then, in local unit, computing unit matrix-vector is long-pending, obtains local unit result vector ${\tilde{v}}_k^e={\tilde{S}}^e{\tilde{p}}_k^e;$

C. last, local unit result vector

by global node, number with the corresponding relation of cell node numbering and aggregate into global outcome vector again.Like this, do not need to generate the Global Vector that overall matrix S has just obtained required calculating.

Therefore, in program, there is no need integrated and store overall matrix, and only need calculate and memory cell matrix.And cell matrix is dense, this has just greatly reduced memory requirements and has improved counting yield.

(3) based on many grids with by the Parallelization Scheme of first technology

Owing to adopting multi-grid technique, spectral element method cell configuration rule, the piecemeal that is easy to parallel algorithm is processed, and each parallel subregion is easy to obtain load balance.Meanwhile, owing to having adopted by first technology, in the spectral element method program of launching based on Chebyshev, the work of the overwhelming majority can be confined to implement in local unit, and can complete whole calculating by polymerization and diffusion.Therefore, adopt many grids and the spectral element program by first technology structure, just to there is good parallel efficiency, possess the special basis that large-scale parallel calculates.

Suppose and carry out elastic wave modeling calculating having in the parallel cluster computer system of P computing node, and adopt MPI(message passing interface) as parallel computation development environment.Carry out after dividing elements, whole computation model has N _eindividual spectral element unit, all gives P cpu process by them, and each process is independently read in the medium parameter in its corresponding subregion, utilizes mass matrix and the element stiffness matrix of multi-grid method independence computing unit, and is stored in the storer of this cpu node.So, in the spatial spreading stage of spectral element method, by the task of cpu node, distribute, reached the object that internal memory is balanced and calculating is balanced, and in this process, between node, do not needed swap data, can reach high efficiency parallel.

The system of equations that Chebyshev spectral element method forms will solve with Conjugate Gradient Method With Preconditioning, and wherein pre-conditional matrix is the matrix that overall diagonal of a matrix forms.In practical operation, without forming overall matrix, as long as by the matrix element value stack on the total node of adjacent cells between the responsible subdomain of each process, then pass to each process, for the corresponding element value of updating block diagonal of a matrix.Can avoid so a large amount of data transmission, improve parallel efficiency.

In elastic wave modeling, a large amount of calculating is present in the iterative computation in time discrete.In each step of time iteration, need by Conjugate Gradient Method With Preconditioning iterative system of equations.Owing to having adopted by first technology, overall matrix wherein and vector product can be converted into cell matrix and vector product, its main calculation procedure is as follows:

A. by superposeing and transmitting the vector value on the total node of adjacent cells between subregion, obtain the part that Global Vector distributes on node;

B. in computing node internal compute unit matrix and vector product, obtain unit result vector, wherein vector is previous step gained, and cell matrix remains unchanged in iterative computation;

C. superpose and transmit the result vector value on the total node of adjacent cells, obtaining the part that global outcome vector distributes on node, for more new explanation and remaining vectorial, entering next iteration.

The two-dimension elastic wave simulation of take is below example, illustrates the embodiment of the parallel spectral element method of many grids Chebyshev.There is discontinuous saltus step in medium parameter wherein in the small scale of space, supposes in parallel cluster computer system and carry out elastic wave modeling calculating with P CPU computing node, and adopt MPI(message passing interface) as parallel computation development environment.

Fig. 1 is the parallel spectral element method process flow diagram of many grids Chebyshev.As shown in Figure 1, specific implementation process has following steps:

(1) zoning is divided into N _e* N _ethe quadrilateral units of individual non-overlapping copies.These unit all can be given to P CPU computing node, each process is independently read in the medium parameter in its corresponding subregion, and the parameter such as exponent number of definite main grid, auxiliary grid, collocation point position.

Main grid exponent number N is selected according to the precision of calculating and time demand etc., for describing auxiliary grid exponent number L and the collocation point position of medium, according to dielectric distribution situation, selects, and for describing auxiliary grid exponent number K and the collocation point of external force, according to external force distribution situation, selects.In the calculating of mass matrix, stiffness matrix, need to use N and L, and the medium parameter reading in, in the calculating of outer force vector, need for N and K.

(2) each process is independently with force vector outside many grid methods computing unit mass matrix, element stiffness matrix and unit.

The computing formula of two-dimension elastic wave simulation is as follows:

$M^e=\left[\begin{array}{cc}{\stackrel{\&OverBar;}{M}}^e& 0\\ 0& {\stackrel{\&OverBar;}{M}}^e\end{array}\right],K^e=\left[\begin{array}{cc}{\stackrel{\&OverBar;}{K}}_1^e& {\stackrel{\&OverBar;}{K}}_2^e\\ {\stackrel{\&OverBar;}{K}}_2^{e\&CenterDot;}& {\stackrel{\&OverBar;}{K}}_3^e\end{array}\right],F^e=\left[\begin{array}{c}{\stackrel{\&OverBar;}{F}}_1^e\\ {\stackrel{\&OverBar;}{F}}_2^e\end{array}\right],---\left(7\right)$

Wherein

${\stackrel{\&OverBar;}{M}}_{i_1{i}_2{j}_1{j}_2}^e=\frac{{\mathrm{\&Delta;}}_1^e{\mathrm{\&Delta;}}_2^e}{4}\underset{l_1,l_2=0}{\overset{L}{\mathrm{\&Sigma;}}}{\widetilde{\mathrm{\&rho;}}}_{l_1{l}_2}^e{B}_{i_1{j}_1{l}_1}{B}_{i_2{j}_2{l}_2}$

${\left[K_1^e\right]}_{i_1{i}_2{j}_1{j}_2}=\underset{l_1,l_2=0}{\overset{L}{\mathrm{\&Sigma;}}}[({\mathrm{\&lambda;}}_{l_1{l}_2}+2{\mathrm{\&mu;}}_{l_1{l}_2})\frac{{\mathrm{\&Delta;}}_2^e}{{\mathrm{\&Delta;}}_1^e}{\stackrel{\stackrel{\&OverBar;}{\&OverBar;}}{B}}_{i_1{j}_1{l}_1}{B}_{i_2{j}_2{l}_2}+{\mathrm{\&mu;}}_{l_1{l}_2}\frac{{\mathrm{\&Delta;}}_1^e}{{\mathrm{\&Delta;}}_2^e}{B}_{i_1{j}_1{l}_1}{\stackrel{\stackrel{\&OverBar;}{\&OverBar;}}{B}}_{i_2{j}_2{l}_2}]$

${\left[K_2^e\right]}_{i_1{i}_2{j}_1{j}_2}=\underset{l_1,l_2=0}{\overset{L}{\mathrm{\&Sigma;}}}({\mathrm{\&lambda;}}_{l_1{l}_2}{\stackrel{\&OverBar;}{B}}_{i_1{j}_1{l}_1}{\stackrel{\&OverBar;}{B}}_{j_2{i}_2{l}_2}+{\mathrm{\&mu;}}_{l_1{l}_2}{\stackrel{\&OverBar;}{B}}_{j_1{i}_1{l}_1}{\stackrel{\&OverBar;}{B}}_{i_2{j}_2{l}_2})---\left(8\right)$

${\left[K_3^e\right]}_{i_1{i}_2{j}_1{j}_2}=\underset{l_1,l_2=0}{\overset{L}{\mathrm{\&Sigma;}}}[({\mathrm{\&lambda;}}_{l_1{l}_2}+2{\mathrm{\&mu;}}_{l_1{l}_2})\frac{{\mathrm{\&Delta;}}_1^e}{{\mathrm{\&Delta;}}_2^e}{B}_{i_1{j}_1{l}_1}{\stackrel{\stackrel{\&OverBar;}{\&OverBar;}}{B}}_{i_2{j}_2{l}_2}+{\mathrm{\&mu;}}_{l_1{l}_2}\frac{{\mathrm{\&Delta;}}_2^e}{{\mathrm{\&Delta;}}_1^e}{\stackrel{\stackrel{\&OverBar;}{\&OverBar;}}{B}}_{i_1{j}_1{l}_1}{B}_{i_2{j}_2{l}_2}]$

${\left[{\stackrel{\&OverBar;}{F}}_1^e\right]}_{i_1{i}_2}=\frac{{\mathrm{\&Delta;}}_1^e{\mathrm{\&Delta;}}_2^e}{4}\underset{k_1,k_2=0}{\overset{K}{\mathrm{\&Sigma;}}}{A}_{i_1{k}_1}{A}_{i_2{k}_2}{\left[f_1^e\right]}_{k_1{k}_2}\left(t\right)$

And

Wherein, element mass matrix M ^e, element stiffness matrix K ^e, force vector F outside unit ^e.ρ is Media density, λ, the Lame coefficient that μ is medium, footnote l ₁l ₂be illustrated in auxiliary grid coordinate in unit. for the Chebyshev interpolation polynomial that main grid is blocked expansion, the interpolation polynomial of φ (ξ) for launching on medium auxiliary grid, the interpolation polynomial of ψ (ξ) for launching on external force auxiliary grid.

(3) independent compute matrix in each process

and get its diagonal entry, with vector mode storage, as the pre-conditional matrix in unit.Label e represents local matrix or the vector on unit.The pre-conditional matrix in diagonal angle, this is a kind of pre-conditional matrix choosing method during Conjugate Gradient Method With Preconditioning solves, and can reach the effect of Fast Convergent, and diagonal matrix can store with vector form, reduces EMS memory occupation, improves counting yield.

(4) between each process, transmit pre-conditional matrix element value the 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, be stored in independently in the internal memory of each cpu node.

(5) each process initialization wave field independently, start time iteration.

(6) 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).$

(7) 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.

(8) use Conjugate Gradient Method With Preconditioning combination 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.

(9) each process solves by previous step the wave field increment δ u obtaining independently _nupgrade wave field and derivative thereof.

(10) judge whether time iteration finishes, if do not finish, return to step (6) and continue iteration.

(11) each process is independent preserves result writing in files.

Fig. 2 is Conjugate Gradient Method With Preconditioning subroutine flow chart.As shown in Figure 2, this process comprises following sub-step:

A. each process is independent initially dissolves, and calculates remaining vector.

B. each process is transmitted remaining vector value the stack that the total node of adjacent cells is corresponding, obtains the remnants vector upgrading.

C. each process independence compute matrix vector product, uses parallel algorithm compute vector inner product.

D. each process is transmitted vector value the stack on the total node of adjacent cells, obtains the matrix-vector product upgrading.

E. more new explanation and remaining vectorial.

Judge whether to reach the condition of convergence, if do not restrain, return to step c loop iteration.

The embodiment of the present invention is owing to adopting multi-grid technique, dividing elements is simple, avoided the division of complicated junior unit in classic method, in unit, can adopt most suitable auxiliary grid and Interpolation-Radix-Function to describe the variation of medium, and wave field solves in more sparse main grid, reduce calculating scale, save computational resource., utilize by first technology meanwhile, can reduce memory space and calculated amount in thousandfold ground, and utilize parallel method to reach higher parallel efficiency.

Take a length of side as L two-dimension square shape zoning Ω be example, the elastic wave minimum wavelength of wherein propagating is λ _min, in medium, have yardstick is 1/10 λ simultaneously _minnon-uniform Distribution.

By traditional spectral element method, calculate, must be divided into N by medium partition of the scale heterogeneous unit ，Jiang zoning Ω _e* N _ethe quadrilateral units of individual non-overlapping copies, the unit length of side is 1/10 λ _min, unit Chebyshev launches to get N rank, and total computing node number is (N _en+1) * (N _en+1), the size of overall matrix is S _g=(2 * (N _en+1) * (N _en+1)) ².With the overall matrix-vector product of calculated load maximum, weigh computation requirement, and only pay close attention to wherein relatively large multiplying consuming time, in traditional spectral element method, needing the multiplying total degree carrying out is T _g=(2 * (N _en+1) * (N _en+1)) ².

And calculate by the parallel spectral element method of many grids Chebyshev, can adopt larger unit, unit size can be greater than the minimum wavelength λ of propagated elastic wave conventionally _min.If get λ herein _minas the unit length of side, zoning is divided into (N _e/ 10) * (N _e/ 10) individual unit, on unit, Chebyshev launches still to get N rank, and cell matrix size is

owing to adopting by first technology, without forming overall matrix, the memory requirements of memory cell matrix is

S _e=(N _e/10)×(N _e/10)×(2×(N+1)×(N+1)) ²。On computation requirement, utilize and by first technology, overall matrix-vector is multiplied each other and is converted into (N _e/ 10) * (N _e/ 10) sub-cell matrix-vector multiplies each other, and required total multiplying number of times is T _e=(N _e/ 10) * (N _e/ 10) * (2 * (N+1) * (N+1)) ².

By matrix size, weigh the memory requirements of two kinds of modes:

$\frac{S_g}{S_e}=\frac{{(2\&times;(N_{e}N+1)\&times;(N_{e}N+1))}^2}{(N_e/10)\&times;(N_e/10)\&times;{(2(N+1)\&times;(N+1))}^2}\&ap;{100N}_e^2{\left(\frac{N}{N+1}\right)}^4---\left(10\right)$

With multiplying number of times in matrix-vector product, weigh the calculated amount of two kinds of modes:

$\frac{T_g}{T_e}=\frac{{(2\&times;(N_{e}N+1)\&times;(N_{e}N+1))}^2}{(N_e/10)\&times;(N_e/10)\&times;{(2\&times;(N+1)\&times;(N+1))}^2}\&ap;{100N}_e^2{\left(\frac{N}{N+1}\right)}^4---\left(11\right)$

If get N _e=100, N=6, that is to say, by adopting many grids and by first technology, memory requirements and computation requirement have all at least reduced by 500,000 times.

In addition, while adopting MPI concurrent development environment, determine that the principal element of parallel efficiency is the data volume of transmitting between load balance and cpu node, and the parallel spectral element method load balancing of many grids Chebyshev, need the data volume of transmission little, can reach more than 90% parallel efficiency.

Above-described embodiment; object of the present invention, technical scheme and beneficial effect are further described; institute is understood that; the foregoing is only the specific embodiment of the present invention; the protection domain being not intended to limit the present invention; within the spirit and principles in the present invention all, any modification of making, be equal to replacement, improvement etc., within all should being included in 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