Spectral element tearing and splicing region decomposition method for embedding multi-pole perfect matching layer technology

Abstract

The invention discloses a spectral element tearing and splicing region decomposition method embedded with a multi-pole perfect matching layer technology, and the efficiency and precision of seismic wave forward modeling simulation are improved through the fusion of the multi-pole perfect matching layer technology and a region decomposition strategy. The time domain convolution of the multi-pole perfect matching layer is eliminated through a recursion auxiliary equation, the boundary truncation calculation complexity is reduced, glancing waves and long-distance propagation waves are efficiently absorbed, and false reflection is restrained. According to the regional decomposition strategy, a global three-dimensional problem is converted into a two-dimensional interface problem through implicit interface constraint and explicit angular point continuity execution, meanwhile, the explicit calculation advantage of a spectral element method on an angular mass matrix is kept, and rapid solution is achieved without matrix inversion. A load balancing mechanism dynamically allocates tasks of a computational domain and a perfect matching layer, and parallel computing efficiency is guaranteed. The system supports high-order interpolation and large-scale grid division of a complex geologic model, appropriately increases memory consumption, improves engineering practicability, and provides an efficient and reliable simulation tool for deep ground resource exploration.

Description

Technical Field

[0001] The present invention belongs to the technical field of computational solid mechanics, and in particular relates to a spectral element tearing and splicing region decomposition method embedded with a multipole perfectly matched layer technology. Background Art

[0002] Geophysical survey methods, as a key technical means for geological structure detection and underground resource extraction, are among the most fundamental, strategic, and critical industries in modern science and technology and have been extensively studied. Seismic exploration methods offer inherent advantages for complex geological structures. However, traditional differential equation numerical methods often face computational inefficiencies when performing forward numerical simulations of open-domain three-dimensional seismic waves. Furthermore, computational resource limitations necessitate the use of specific artificial boundaries to truncate the computational domain, specifically eliminating a large number of non-physical false reflections. The spectral element method, due to its significant advantages in explicit computation, has become one of the most popular algorithms for seismic wave forward modeling. However, when the problem scale reaches a certain level, parallel computing remains an important acceleration strategy that must be considered. Unfortunately, commonly used numerical acceleration methods require restructuring the computational format, significantly increasing implementation difficulty because interpolation at the interface between subdomains can affect the diagonal properties of the matrix. The addition of truncation boundaries will consume additional computing resources, further exacerbating the difficulty of solving the problem. Taking the perfectly matched layer with the best absorption performance as an example, it not only requires the allocation of grids like the computational domain, but the multi-pole perfectly matched layer will also cause complex convolution operations in the time domain, which makes the computational time consumption span between the perfectly matched layer area and the computational domain huge. Even using the regional decomposition strategy for parallel computing, it is difficult to balance the load between different processes. A method to solve the above problems is urgently needed. Summary of the Invention

[0003] In order to solve the above technical problems, the present invention proposes a spectrum element tearing and splicing region decomposition method embedded with multipole perfectly matched layer technology to solve the problems existing in the above-mentioned prior art.

[0004] In a first aspect, to achieve the above-mentioned objectives, the present invention provides a method for spectral element tearing and splicing region decomposition embedded in a multipole perfectly matched layer technology, comprising the following steps:

[0005] Construct the geological model required for simulation based on the actual geological characteristics, add the perfectly matched layer area in the normal direction of the truncation boundary and perform mesh division to generate a grid file containing grid information;

[0006] Set up the simulation configuration file to define the number of subdomains, material properties, boundary types, time steps, and source and receiver locations;

[0007] Distribute the process load between the computational domain and the perfectly matched layer area according to the configuration file to balance the workload of each process;

[0008] Construct the mapping relationship between all subdomains and the global system, divide the subdomain nodes into internal nodes, interface nodes and corner nodes, and establish point, surface and unit mapping between subdomain entities;

[0009] Based on the weak form of the second-order elastic wave equation, the mass matrix and stiffness matrix are assembled for the subdomains in parallel, and only the mass matrix is ​​assembled for the perfectly matched layer subdomain, and its stiffness matrix is ​​loaded into the right-hand side in vector form through a recursive auxiliary differential equation.

[0010] Initialize displacement, velocity, acceleration field variables and load vectors, update field variables through time iteration, and modify displacement and stress in real time by combining perfectly matched layer recursive auxiliary differential equations;

[0011] Solve the global interface problem, assign the interface dual variables and corner point primitive variables to the subdomains, and complete the wavefield evolution calculation.

[0012] Optionally, a three-dimensional geological model is generated based on the actual geological conditions, and the multipolar perfectly matched layer region is expanded in the normal direction of the model boundary;

[0013] The computational domain and the perfectly matched layer region are divided into multiple subdomains. Hexahedral meshing is performed on each subdomain to generate a mesh file containing the subdomain material properties, boundary types, node coordinates, and element topology relationships.

[0014] Optionally, the process of setting the simulation configuration file includes:

[0015] Define the grid file storage path, number of subdomains, elastic parameters and density of each subdomain;

[0016] Specify the global simulation time step, total duration, spatial location of the earthquake source and waveform form, and mark the receiving point location to record the field variable waveform;

[0017] Mark the boundary type between the computational domain and the perfectly matched layer region and generate a configuration file to control the simulation process.

[0018] Optionally, the process of constructing a mapping relationship between the subdomain and the global system includes:

[0019] Define a signed Boolean matrix for each subdomain, map the subdomain interface dual variables to the global dual variables, and ensure that the variables on both sides of the same interface have opposite signs;

[0020] Define an unsigned Boolean matrix for each subdomain, map the subdomain corner points to the global corner points, and eliminate the non-uniqueness of the solutions of shared nodes;

[0021] Establish mapping relationships between entities within the subdomain, including the correspondence between local and global numbers of nodes, faces, and hexahedral elements.

[0022] Optionally, the process of assembling matrices for subdomains in parallel includes:

[0023] For the computational domain subdomain, the mass matrix, stiffness matrix and load matrix are assembled in parallel based on the Gauss-Lobatto-Legendre interpolation points;

[0024] For the perfectly matched layer subdomain, only the mass matrix is ​​assembled, and its stiffness matrix is ​​loaded to the right-hand side by a vector generated by a recursive auxiliary differential equation. The recursive equation is constructed based on a high-order complex stretching function, and the modified stress variable participates in the load update in vector form.

[0025] Optionally, the process of updating the field variables through time iteration and combining the perfectly matched layer recursive auxiliary differential equation to correct the displacement and stress in real time includes:

[0026] Based on the Newmark-beta prediction formula, update the intermediate variables of the subdomain displacement and velocity at the current moment, and update the right-hand side of the load;

[0027] For the perfectly matched layer subdomain, the displacement variables and historical auxiliary variables are input into the correction function, and the corrected displacement, stress and load vectors at the current moment are generated through the recursive auxiliary differential equation. The right-hand side terms are updated and the auxiliary variables are saved.

[0028] The corrected right-hand side term is loaded into the global interface problem. The interface dual variables and corner point primal variables are solved and then distributed to the subdomains. The acceleration and velocity variables are recovered by combining the Newmark-beta correction formula.

[0029] Optionally, the process of solving the global interface problem, assigning the interface dual variables and corner primitive variables to the subdomains, and completing the wavefield evolution calculation includes:

[0030] Associate displacement, velocity, and acceleration field variables with their spatial coordinates, and write the field variables to disk at a preset time;

[0031] Perform spatial interpolation on the field variables at the receiving point position, extract the displacement waveform data and generate the simulation waveform;

[0032] Calculate the L2 norm error between the simulation waveform and the reference solution to verify the correctness of the simulation results.

[0033] In a second aspect, the present invention further provides a spectral element tearing and splicing region decomposition system embedded in a multipole perfectly matched layer technology, for implementing a spectral element tearing and splicing region decomposition method embedded in a multipole perfectly matched layer technology, the system comprising:

[0034] The geological modeling module is used to construct the 3D geological model required for simulation based on the real geological characteristics, expand the multipolar perfectly matched layer region in the normal direction of the truncation boundary, and generate a grid file containing subdomain material properties, boundary type and grid information;

[0035] The simulation configuration module is used to define the number of subdomains, material parameters, global time step, total duration, source and receiver locations, and generate configuration files to control the simulation process;

[0036] The load distribution module is used to dynamically distribute process tasks to balance the load based on the computational domain and perfectly matched layer area marked in the configuration file, the number of computing units and the field quantity to correct the time consumption;

[0037] The mapping relationship building module is used to divide the subdomain nodes into internal nodes, interface nodes and corner nodes, and establish the Boolean matrix mapping relationship between the subdomain entities and the global system;

[0038] Matrix assembly module, used to assemble the mass matrix and stiffness matrix for subdomains based on the weak form of the elastic wave equation. For the perfectly matched layer subdomain, only the mass matrix is ​​assembled, and its stiffness matrix is ​​loaded to the right-hand side through the vector generated by the recursive auxiliary differential equation;

[0039] The field variable iteration module is used to initialize the displacement, velocity, acceleration field variables and load vectors, update the field variables through time iteration, and correct the displacement and stress in real time by combining the multipole perfectly matched layer recursive auxiliary differential equation;

[0040] The global solution module is used to construct and solve the global interface problem, assign interface dual variables and corner point primitive variables to subdomains, and complete the wavefield evolution calculation.

[0041] In a third aspect, the present invention further provides a computer terminal device, comprising:

[0042] one or more processors;

[0043] a memory, coupled to the processor, for storing one or more programs;

[0044] When the one or more programs are executed by the one or more processors, the one or more processors implement a spectral element tearing and splicing region decomposition method embedded with a multipole perfectly matched layer technology.

[0045] In a fourth aspect, the present invention further provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements a spectral element tearing and splicing region decomposition method embedded with a multipole perfectly matched layer technology.

[0046] Compared with the prior art, the present invention has the following advantages and technical effects:

[0047] The present invention provides a method for spectral element tearing and splicing regional decomposition embedded with multi-pole perfect matching layer technology. The present invention significantly improves the efficiency and accuracy of seismic wave forward simulation through the deep integration of multi-pole perfect matching layer technology and regional decomposition strategy. The multi-pole perfect matching layer eliminates time domain convolution operations through recursive auxiliary equations, greatly reduces the computational complexity of boundary truncation, achieves efficient absorption of grazing waves and long-distance propagation waves, and effectively suppresses false reflections; the regional decomposition strategy converts the global three-dimensional problem into a two-dimensional interface problem through implicit interface constraints and explicit corner point continuity execution, retains the explicit calculation advantage of the diagonal mass matrix of the spectral element method, and can be quickly solved without matrix inversion. At the same time, the load balancing mechanism dynamically allocates the computational domain and perfect matching layer tasks to ensure parallel computing efficiency; the system supports high-order interpolation and large-scale grid division of complex geological models, significantly improving engineering practicality while moderately increasing memory consumption, and providing an efficient and reliable simulation tool for deep-earth resource exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] The accompanying drawings, which constitute part of the present invention, are provided to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are provided to explain the present invention and do not constitute an undue limitation of the present invention. In the accompanying drawings:

[0049] Figure 1 This is a flow chart of a regional decomposition SETI-DP method embedded in a multipole perfectly matched layer for seismic wave forward modeling according to an embodiment of the present invention;

[0050] Figure 2 This is a geometric diagram of a method for simulating the effect of infinitely propagating seismic waves in Example 1 of the first embodiment of the present invention;

[0051] Figure 3 Schematic diagram of the field distribution of the x-direction displacement at t=0.035s in Case 1 of Example 1 of the present invention;

[0052] Figure 4 Schematic diagram of the comparison results between the received waveform of the displacement field in the x-direction at the receiving point and the analytical solution in Case 1 of Example 1 of the present invention;

[0053] Figure 5 This is a geometric diagram for simulating the absorption effect of a multipole perfectly matched layer on grazing waves and its ability to cut off long-distance propagation of seismic waves in Case 2 of Example 1 of the present invention;

[0054] Figure 6 Schematic diagram of the field distribution of the z-direction displacement at time t=0.400s in Case 2 of Example 1 of the present invention;

[0055] Figure 7Schematic diagram of the field distribution of the z-direction displacement at time t=0.800s in Case 2 of Example 1 of the present invention;

[0056] Figure 8 Schematic diagram of the field distribution of the z-direction displacement at time t=1.300s in Case 2 of Example 1 of the present invention;

[0057] Figure 9 Schematic diagram of the field distribution of the z-direction displacement at time t=1.700s in Case 2 of Example 1 of the present invention;

[0058] Figure 10 Schematic diagram showing the comparison results between the received waveform of the displacement field in the z direction at the receiving point in Case 2 in Example 1 of the present invention and the traditional method;

[0059] Figure 11 This is a geometric diagram of Case 3 in Example 1 of the present invention;

[0060] Figure 12 Schematic diagram of the field distribution of the z-direction displacement at time t=1.600s in Case 3 of Example 1 of the present invention;

[0061] Figure 13 Schematic diagram of the field distribution of the z-direction displacement at time t=1.800s in Case 3 of Example 1 of the present invention;

[0062] Figure 14 Schematic diagram of the comparison results between the received waveform of the displacement field in the z direction at the receiving point of Case 3 in Example 1 of the present invention and the traditional method. DETAILED DESCRIPTION

[0063] It should be noted that, in the absence of conflict, the embodiments and features of the embodiments of the present invention can be combined with each other. The present invention will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0064] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0065] First, the technical terms involved in the following embodiments are explained.

[0066] The present invention relates to an efficient numerical simulation method for large-scale complex open-domain seismic wave forward modeling scenarios, in particular to a dual-primal spectral element tearing and interconnecting (SETI-DP) method that combines the high precision of a spectral method with the regional decomposition characteristics of a finite element tearing and merging method and embeds a multi-pole perfectly matched layer (MPML) technology.

[0067] Example 1

[0068] like Figure 1 As shown, this embodiment provides a spectrum element tearing and splicing region decomposition method embedded with multipole perfectly matched layer technology, including:

[0069] Construct the geological model required for simulation based on the actual geology, add the perfectly matched layer region in the normal direction of the truncation boundary and perform mesh division to generate a mesh file containing mesh information;

[0070] Set up the simulation configuration file to define the number of subdomains, material properties, boundary types, time steps, and source and receiver locations;

[0071] Distribute the process load between the computational domain and the perfectly matched layer area according to the configuration file to balance the workload of each process;

[0072] Construct the mapping relationship between all subdomains and the global system, divide the subdomain nodes into internal nodes, interface nodes and corner nodes, and establish point, surface and unit mapping between subdomain entities;

[0073] Based on the weak form of the second-order elastic wave equation, the mass matrix and stiffness matrix are assembled for the subdomains in parallel, and only the mass matrix is ​​assembled for the perfectly matched layer subdomain, and its stiffness matrix is ​​loaded into the right-hand side in vector form through a recursive auxiliary differential equation.

[0074] Initialize displacement, velocity, acceleration field variables and load vectors, update field variables through time iteration, and modify displacement and stress in real time by combining perfectly matched layer recursive auxiliary differential equations;

[0075] Solve the global interface problem, assign the interface dual variables and corner point primitive variables to the subdomains, and complete the wavefield evolution calculation.

[0076] Specifically, the present invention addresses the technical difficulties existing in seismic wave forward modeling and proposes a spectral element tearing and splicing (SETI-DP) regional decomposition method embedded with multipole perfectly matched layer technology.

[0077] The present invention provides a technical solution to solve the existing technical difficulties. The technical solution is as follows: using the concept of multi-pole perfectly matched layers and the recursive relationship of high-order complex stretching functions, a first-order recursive auxiliary differential equation applicable to any high-order perfectly matched layer is constructed. Specific auxiliary variables are added to iteratively correct the field variables in the perfectly matched layer at each moment to achieve a truncation effect. Different stiffness matrix calculation strategies are adopted in the computational domain and the perfectly matched layer region to achieve the purpose of directly separating the computational domain and the perfectly matched layer region. Lagrange multipliers are introduced at the interfaces of all subdomains, including the perfectly matched layer subdomain, and their continuity at the interfaces is implicitly enforced to achieve the purpose of field quantity continuity. On this basis, in order to eliminate the non-uniqueness of solutions at nodes shared by multiple subdomains, the concept of corner points is introduced to explicitly enforce the continuity of field quantities. Through these measures, the original three-dimensional problem is perfectly equivalent to a two-dimensional problem involving only interfaces. Finally, a unified parallel computing framework with direct separation of the computational domain and the perfectly matched layer and load balancing is generated.

[0078] As an implementation method in this embodiment, a three-dimensional geological model is generated based on real geological conditions, and the multipolar perfectly matched layer region is expanded in the normal direction of the model boundary;

[0079] The computational domain and the perfectly matched layer region are divided into multiple subdomains. Hexahedral meshing is performed on each subdomain to generate a mesh file containing the subdomain material properties, boundary types, node coordinates, and element topology relationships.

[0080] Specifically, in step one, a geological model required for simulation is constructed based on the actual geological conditions, and a perfectly matched layer region is added in the normal direction of the model boundary that needs to be truncated. Finally, meshing is performed, and a mesh file containing various mesh information is generated and exported.

[0081] More specifically, in step 1, the geological model required for simulation is constructed according to the actual geological conditions, and the perfect matching layer area is added in the normal direction of the model boundary that needs to be truncated. Finally, the grid is divided and the grid file containing various grid information is generated and exported. The specific geological model is as follows Figure 2 、 Figure 5 、 Figure 11 As shown, the generated mesh model can be seen Figure 3 、 Figure 6 、 Figure 12 .

[0082] As an implementation method of this embodiment, the process of setting the simulation configuration file includes:

[0083] Define the grid file storage path, number of subdomains, elastic parameters and density of each subdomain;

[0084] Specify the global simulation time step, total duration, spatial location of the earthquake source and waveform form, and mark the receiving point location to record the field variable waveform;

[0085] Mark the boundary type between the computational domain and the perfectly matched layer region and generate a configuration file to control the simulation process.

[0086] Specifically, in step 2, a simulation configuration file is set, including the grid storage path, the number of subdomains, subdomain materials, subdomain boundary types, global simulation time step, global simulation total duration, the locations of source and receiving points, and the form of the source.

[0087] More specifically, in step 2, you set up a simulation configuration file, including the grid storage path, the number of subdomains, subdomain materials, subdomain boundary types, global simulation time step, total global simulation duration, the locations of source and sink points, and the source type. This file controls the entire simulation environment and process.

[0088] As an implementation method of this embodiment, the following is further included: Step 3, running the configuration file, and reasonably distributing the load of each process based on the calculation area and perfect matching layer area marked in the configuration file, combined with the number of computing units and the time consumption of the perfect matching layer for field quantity correction, so that the task amount allocated to each process is at a similar level.

[0089] Specifically, in step 3, the configuration file is run. Based on the calculation area and the perfect matching layer area marked in the configuration file, combined with the number of calculation units and the time consumption of the perfect matching layer for field correction, the load of each process is reasonably distributed so that the amount of tasks assigned to each process is equal. The specific division results can be found in Figure 3 、 Figure 6 、 Figure 12 It is observed that Tables 1, 2, and 3 also give the maximum and minimum perfectly matched layer subdomains contained in the regional decomposition SETI-DP method using embedded multipolar perfectly matched layers.

[0090] As an implementation method of this embodiment, the process of establishing a mapping relationship between a subdomain and a global system includes:

[0091] Define a signed Boolean matrix for each subdomain, map the subdomain interface dual variables to the global dual variables, and ensure that the variables on both sides of the same interface have opposite signs;

[0092] Define an unsigned Boolean matrix for each subdomain, map the subdomain corner points to the global corner points, and eliminate the non-uniqueness of the solutions of shared nodes;

[0093] Establish mapping relationships between entities within the subdomain, including the correspondence between local and global numbers of nodes, faces, and hexahedral elements.

[0094] Specifically, in step four, the mapping relationship between all subdomains and the global problem is constructed, and the nodes of all subdomains are divided into three categories corresponding to internal nodes, interface nodes and corner points respectively; as well as the mapping relationship between all entities within all subdomains, including the mapping relationship between points, surfaces and units.

[0095] More specifically, step four involves constructing a mapping relationship between all subdomains and the global problem. This involves classifying all subdomain nodes into three categories: internal nodes, interface nodes, and corner nodes. This also involves mapping relationships between all entities within each subdomain, including between points, surfaces, and cells. The mapping relationship between subdomains and the global problem is key to building the SETI-DP algorithm framework. Because this step requires correlation calculations on the global grid, it can only be implemented serially. The specific data required is:

[0096] In order to construct the global problem, we introduce the signed Boolean matrix [B i Sl ], and an unsigned Boolean matrix [B i Cl ], where [B i Sl ] is used to establish the mapping relationship between the dual unknowns contained in the subdomain i in the l direction and the global dual unknowns; [B i Cl ] is used to establish the mapping relationship between the corner points contained in the subdomain i in the l direction and the global corner points. li , C li , S l , C l represents the number of dual variables of subdomain i in the l direction, the number of corner points of subdomain i in the l direction, the number of dual variables of the global l direction, and the number of corner points of the global l direction. i Sl ] has a dimension of S li ×S l , the symbols contained only need to satisfy the opposite sides of the same interface. i Cl ] is C li ×C l Specifically in three dimensions, [B i S ]=diag{[B i Sx ],[B i Sy ],[B i Sz ]},[B i C ]=diag{[B i Cx],[B i Cy ],[B i Cz ]}. And according to the definition of Lagrange multiplier, we can also find the following relationship:

[0097]

[0098] As an implementation method of this embodiment, the process of assembling matrices for sub-domains in parallel includes:

[0099] For the computational domain subdomain, the mass matrix, stiffness matrix and load matrix are assembled in parallel based on the Gauss-Lobatto-Legendre interpolation points;

[0100] For the perfectly matched layer subdomain, only the mass matrix is ​​assembled, and its stiffness matrix is ​​loaded to the right-hand side by a vector generated by a recursive auxiliary differential equation. The recursive equation is constructed based on a high-order complex stretching function, and the modified stress variable participates in the load update in vector form.

[0101] Specifically, in step five, based on the weak form of the second-order elastic wave equation, the system matrices are assembled for all subdomains in parallel, including the necessary components required for simulation, such as the mass matrix, stiffness matrix, and load matrix. If it is a perfectly matched layer area, only the mass matrix needs to be assembled.

[0102] More specifically, step five is to assemble the system matrix for all subdomains based on the weak form of the second-order elastic wave equation, including the mass matrix, stiffness matrix, and load matrix, which are necessary components for simulation. If it is a perfectly matched layer region, only the mass matrix needs to be assembled, because the stiffness matrix of the perfectly matched layer region needs to be obtained through recursive auxiliary differential equations and presented in the form of vectors. The control equations of the present invention in the computational domain are based on linear elastic theory, including momentum conservation equations, constitutive relations, and infinitesimal strain-displacement relations:

[0103]

[0104] σ ij =c ijkl ε ij (8)

[0105]

[0106] Subscripts i, j, k, l = 1, 2, 3 are coordinate axis indices, x i Represents the spatial coordinate, u i is the displacement vector u, f i represents the load vector, c ijkl is the fourth-order elasticity tensor, ε ij and σ ijare the strain and tension tensors, respectively.

[0107] The momentum formula (1) is tested using the test function, and the following is obtained by partial integration within the calculation area:

[0108]

[0109] Among them, <> Ω and <> Г denote the volume integral and area integral of the computational domain, respectively. The form of the term (n.σ) in the area integral depends on the type of boundary surface. In SEM, the computational domain is divided into hexahedral elements, each of which is mapped to a standard reference element (ξ, η, ζ)∈[-1,1]×[-1,1]×[-1,1] and a Gauss-Lobatto-Legendre (GLL) interpolation point. Then, the variables u and σ in Equation (10) can be expanded in the reference element as follows:

[0110]

[0111] Where u i (n) and σ ij (n) are the coefficients of the displacement vector and stress vector components respectively. The subscripts i, j = 1, 2, 3 are the indices of the local coordinates. N and N s are the total number of node basis functions in the computational domain Ω and on the boundary Г, respectively. n represents the GLL basis function, and its one-dimensional p-order form can be expressed as:

[0112]

[0113] Among them, L p (ξ) is the pth order Legendre polynomial, and L' p (ξ) is its derivative. j , j = 0, 1..., p} in the reference element ξ∈[-1, 1] are called GLL nodes, which are p-order polynomials (1-ξ j 2) L' p (ξ j )=0.

[0114] Through formulas (11)-(13), the matrix equation can be expressed as:

[0115]

[0116] The specific expressions of all the above matrices and vectors are as follows

[0117]

[0118] Among them, C kl It is a 3×3 matrix whose elements come from the fourth-order elasticity tensor. The subscripts k, l = 1, 2, and 3 represent the indices of the variable components respectively.

[0119] As an implementation method of this embodiment, the following further includes: Step 6, calculating and setting relevant perfect matching layer characteristic parameters based on the distance from the integration point to the computational domain in combination with the priori parameters, and then allocating auxiliary variables to all perfect matching layer subdomains to store the time-related auxiliary variables {U N} 0 ,{σ N} 0 .

[0120] Specifically, in step 6, the relevant PML characteristic parameters are calculated and set according to the distance from the integration point to the computational domain in combination with the prior parameters, and then auxiliary variables are allocated to all PML subdomains to store the time-dependent auxiliary variables {U N} 0 ,{σ N} 0 .

[0121] The complex stretching function can be transferred directly from the partial derivatives to the physics variables to avoid complicated convolution operations, as shown below:

[0122]

[0123] The change of this variable is not strictly correct because the complex stretching function varies with space, but this change does not substantially affect the performance of the perfect matching layer. So the frequency domain expression of the perfect matching layer region can be obtained:

[0124]

[0125] Assuming that the n-order perfectly matched layer region has the corrected displacement and stress variables, they are expressed as and Where i, j, k represent the directions of the original field variables, l, j are the directions of the corresponding variable stretching, so the governing equations of the perfectly matched layer region can be obtained from equations (19) and (20) by inverse Fourier transform in the time domain:

[0126]

[0127] For a perfectly matched layer of order n, the strain tensor should be solved using the displacement variable after stretching, with the following relationship:

[0128]

[0129] Next, the stretching process of the displacement variable and strain variable in the n-order perfectly matched layer region will be derived. For the displacement variable u in the perfectly matched layer domain, k Without losing generality, let’s take the stretching function along the l direction as an example:

[0130]

[0131] Among them, S (N) l is the N-order stretching function in the l direction, is the modified displacement variable of order N in the l direction. There is the following recursive relationship,

[0132]

[0133] in

[0134]

[0135] The above equation shows that the jth correction variable Can be combined with the j-1th variable u j-1 kl Without losing generality, all first-order stretching functions in formula (24) adopt CFS complex stretching functions. The j-th CFS complex stretching function can be written as:

[0136]

[0137] in

[0138]

[0139] Substitute formula (27) into formula (25) and introduce the auxiliary memory variable F j-1 , we get the following recursive equation:

[0140]

[0141] in

[0142]

[0143] Since the parameter α j l ,β j l ,d j l , has nothing to do with frequency, so formula (33) can be converted in the time domain to obtain:

[0144]

[0145] Substituting formula (34) into formula (27) yields the final recursive auxiliary differential equation:

[0146]

[0147] Similarly, we can easily get σ N ij(j) The recursive auxiliary differential equation is as follows:

[0148]

[0149] Then, eliminating stress variables from equations (22)-(23) and using the Galerkin method, we can obtain the weak form of the perfect matching layer region as follows:

[0150]

[0151] By comparing formula (41) with formula (14), it can be seen that the mass matrix of the perfectly matched layer region has the same form as that of the computational domain. However, the stiffness matrix of the perfectly matched layer region is different from that of the computational domain. After obtaining the displacement and stress variables in the previous step, the corresponding auxiliary variables need to be added for correction and finally reflected in the right-hand side of the matrix equation in the form of a vector. The specific form will be derived next.

[0152] First, the displacement variable u in the PML region becomes u after correction N , expanded in the reference element can be expressed as:

[0153]

[0154] Where L represents the number of integration nodes, and:

[0155]

[0156] U i N In order to use the modified displacement variables of formulas (35)-(37), this is different from the treatment method of separating the three directions in the computational domain, but it is very convenient to unify the two. Substituting formula (43) into formulas (21)-(23), the expression of the node stress field is obtained based on the discretized modified displacement field:

[0157]

[0158] in:

[0159]

[0160] Substituting Equation (46) into Equations (38)-(40), the nine modified stress variables acting on the nodes are as follows:

[0161]

[0162] Substituting formula (47) into the stiffness matrix term in formula (41) and performing volume integration, we can obtain the internal nodal force F corresponding to the modified stress field: PML .

[0163]

[0164] in:

[0165] Γ(ξ i )=Nσ N (ξ i ) (49)

[0166]

[0167] Finally, the matrix equation of the perfectly matched layer region is obtained as follows:

[0168]

[0169] Two different calculation strategies are used in the computational domain and the perfectly matched layer region, which ultimately leads to inconsistent order of field quantities under the two strategies. To solve this problem, the order needs to be reversed using the projection matrix before and after the field variables in the perfectly matched layer region are corrected.

[0170] As an implementation method of this embodiment, it also includes: Step seven, based on the mapping relationship between the nodes in all subdomains (including the perfectly matched layer subdomain) and the global system, constructing the system matrix of the global matrix equation for solving the global problem in subsequent time iterations, which are the system matrices of the global dual variables and the global primal variables respectively.

[0171] Specifically, in step seven, based on the mapping relationship between the nodes in all subdomains (including the perfectly matched layer subdomain) and the global system, the system matrix of the global matrix equation is constructed for solving the global problem in subsequent time iterations, which are the system matrices of the global dual variables and the global primal variables respectively.

[0172] The current computational domain and the perfectly matched layer domain are divided into N subdomains. ij (supports arbitrary shapes) Simultaneously applying Dirichlet and Newman-type transmission conditions can ensure the continuity of the field quantity:

[0173]

[0174] where n i and n jRepresent the unit normal vectors of subdomain i and subdomain j pointing to each other. Formula (49) ensures the continuity of acceleration at the interface between subdomain i and subdomain j. At this time, the continuity of displacement and velocity is also guaranteed. The Newman boundary data at the interface is actually unknown, so σ on the interface of the perfectly matched layer region is N No special treatment is required. The solution of each subdomain can be found under the constraint of normal acceleration. This type of constraint problem is called constrained variation problem in mathematics and can be solved using the Lagrange multiplier method.

[0175] Therefore, formula (14) and formula (51) can be written in a unified form and the Lagrange multiplier term is added to obtain the following matrix equation of the i-th subdomain at time n:

[0176]

[0177] At this time, the nodes of the second-order elastic wave equation are divided into three categories (internal nodes I, interface nodes S, and corner points C). Here, the corner points are the edges where any two interfaces of each subdomain intersect.

[0178] Formula (54) can be expanded to obtain the following three equations:

[0179]

[0180] Due to the diagonal advantage of SEM explicit calculation, the subdomain matrix equation is successfully decoupled into the three equations above.

[0181] Multiply both ends of formula (56) by [B i S ] and sum it across all interfaces to get the following expression:

[0182] [F SS ]{λ n}={d S} (58)

[0183] The above formula is called the global interface system, where:

[0184]

[0185] Multiply the left and right sides of formula (57) by [B i C ] T And summing over all corner points, we can get the following expression:

[0186]

[0187] The above formula is called the global corner point system, where

[0188]

[0189] Matrix [F SS ] and [F CC ] also has a fully diagonal format, so neither the global problem nor the local variable recovery requires matrix inversion. Finally, after solving Equations (58) and (61), we obtain the interface dual variables and the corner primal variables. These are substituted back into Equations (56) and (57) to solve for the interface acceleration and corner acceleration. Further solving Equation (55) yields the internal acceleration. This step will be iterated over time in Step 8.

[0190] As an implementation method of this embodiment, the process further includes: Step 8 of updating the field variables through time iteration and correcting the displacement and stress in real time in combination with the perfectly matched layer recursive auxiliary differential equation includes:

[0191] Based on the Newmark-beta prediction formula, update the intermediate variables of the subdomain displacement and velocity at the current moment, and update the right-hand side of the load;

[0192] For the perfectly matched layer subdomain, the displacement variables and historical auxiliary variables are input into the correction function, and the corrected displacement, stress and load vectors at the current moment are generated through the recursive auxiliary differential equation. The right-hand side terms are updated and the auxiliary variables are saved.

[0193] The corrected right-hand side term is loaded into the global interface problem. The interface dual variables and corner point primal variables are solved and then distributed to the subdomains. The acceleration and velocity variables are recovered by combining the Newmark-beta correction formula.

[0194] As an implementation method of this embodiment, the process of solving the global interface problem, allocating the interface dual variables and the corner point original variables to the subdomains, and completing the wavefield evolution calculation includes:

[0195] Associate displacement, velocity, and acceleration field variables with their spatial coordinates, and write the field variables to disk at a preset time;

[0196] Perform spatial interpolation on the field variables at the receiving point position, extract the displacement waveform data and generate the simulation waveform;

[0197] Calculate the L2 norm error between the simulation waveform and the reference solution to verify the correctness of the simulation results.

[0198] Specifically, in step eight, all field variables are initialized including the initial displacement variable {u i} 0 、Initial velocity intermediate variable {v i*} 0 、Initial velocity variable {v i} 0 、Initial acceleration variable {ai} 0 And the load vector {b i} 0 , and perform wave field simulation evolution calculation:

[0199] A. Initialize the right-hand side of the global problem {d r},{d s}, update the displacement variable {u i} n and the speed intermediate variable {v i*} n , and update the right-hand side items of all subdomains according to the current load {b i} n .

[0200]

[0201] v * =v n-1 +0.5Δta n-1 (2)

[0202] B. If the current subdomain is a multi-pole perfectly matched layer subdomain, the updated displacement variable {u i} n And the auxiliary variables defined in step 6 {U N} n-1 ,{σ N} n-1 The multipole perfect matching layer correction function is sent to correct the field variables, and the returned auxiliary variables {U N} n ,{σ N} n Save and load vector {F i PML} is assigned to the right-hand item.

[0203] C. Construct the right-hand side of the global problem {d r},{d s}, forming a complete system of linear equations for the global interface problem and solving it to obtain the dual variables on all interfaces and the primal variables at all corner points.

[0204] D. According to the mapping relationship between all subdomains and the global problem, the original variables and dual variables obtained by solving the global problem are assigned to the corresponding subdomains, and the acceleration variables {a i} nand the speed variable {v i} n .

[0205]

[0206] E. Match all field variables with their coordinates and write them to disk at the set time for subsequent visualization processing. At the receiving point, interpolate the field variables at the set position and record them to generate the final simulation waveform.

[0207] F. Update all field variables to the evolution results of the previous step and return to step A to continue the waveform calculation at the next moment.

[0208] Step nine: visualize the calculation results and calculate the error between them and the reference solution to verify the correctness of the simulation results.

[0209] The above is the complete process of analyzing seismic wave forward propagation using the regional decomposition SETI-DP method embedded in a multipole perfectly matched layer. Based on this theory, it is programmed and simulated. To verify the feasibility, accuracy, and efficiency of the algorithm proposed in this paper, three numerical examples are given. The numerical results obtained by this method are compared with the analytical solution and the results of the traditional method, and the error analysis is performed. The complete implementation process of this algorithm is shown below:

[0210] (1) The geological model required for simulation is constructed according to the real geological conditions, the perfect matching layer area is added, and finally the grid is generated and the grid file containing various grid information is exported.

[0211] (2) Set the simulation configuration file, which controls the environmental conditions and simulation process of the entire simulation.

[0212] (3) According to the calculation area and perfect matching layer area marked in the configuration file, the grid file is imported into the algorithm proposed in the invention to reasonably distribute the load of each process.

[0213] (4) Construct the mapping relationship between all subdomains and the global problem respectively, and divide all nodes of the subdomains into three categories corresponding to internal nodes, interface nodes and corner points respectively; as well as the mapping relationship between all entities within all subdomains, including the mapping relationship between points, surfaces and units.

[0214] (5) Based on the weak form of the second-order elastic wave equation and the second-order governing equation of the multipole perfectly matched layer, the system matrices are assembled for all subdomains, including the mass matrix, stiffness matrix, and load matrix.

[0215] (6) Assign auxiliary variables to all perfectly matched layer subdomains and set the relevant perfectly matched layer characteristic parameters.

[0216] (7) According to the mapping relationship between the entities and the global system in all computational domains and perfectly matched layer regions, the system matrix of the global matrix equation is constructed;

[0217] (8) Initialize all field variables and perform wave field simulation evolution calculations;

[0218] (9) Visualize the calculation results.

[0219] In order to further verify the effectiveness of the present invention in actual geological models, the present invention is applied to three actual cases to verify its accuracy, efficiency and flexibility.

[0220] Case 1: Infinitely uniform medium model.

[0221] The geometry of the model can be found in Figure 2 The material parameters: P wave velocity, S wave velocity and formation density are set to 3000m / s, 1900m / s, and 2500kg / m^3 respectively, and the parameters of the bipolar perfect matching layer complex stretching function are p d |1=3,R0|1=0.0005,α0|1=πf0,p α |1=1,β0|1=1,p β |1=2, d0|1 is calculated by formula (31), v max =4000m / s, L=18m; p d |2=3, R0|2=0.0005, α0|2=0, p α |2=1,β0|2=1,p β |2=2,d0|2=d0|1 / 20.

[0222] from Figure 3 It can be seen that the outgoing wave is well absorbed and Figure 4 It can be seen that the error between the results calculated by the SETI-DP solver and the analytical solution is very small, and the calculated L-2 norm error is only 0.23%, which proves the high precision characteristics of the SETI-DP solver and the absorption performance of the multipole perfectly matched layer.

[0223] For Case 1, the comparison of computational resource consumption between the SETI-DP solver and the traditional method is shown in Table 1.

[0224] Table 1

[0225]

[0226] Table 1 shows that the SETI-DP solver embedded in a multipole perfectly matched layer is significantly more efficient than the traditional solver while retaining the accuracy of the traditional solver. Specifically, SETI-DP is 15 times more efficient than the traditional solver while only slightly increasing memory consumption, demonstrating the high performance of the SETI-DP solver.

[0227] Case 2: Slender infinitely uniform model.

[0228] The specific geometry can be found in Figure 5 In the case 1, except that the thickness of the perfectly matched layer is increased to L = 100 m, the other parameters are the same as those in case 1.

[0229] from Figure 6 、 Figure 7 、 Figure 8 、 Figure 9 It can be seen that when the long-distance seismic wave forward simulation is carried out, the multi-pole perfect matching layer shows extremely strong absorption performance, and no obvious reflection is observed during the whole process. Figure 10 Compared with the traditional method, the SETI-DP solver shows a similar accuracy, and the L-2 norm error is only calculated to be 1.1158×10 -10 .

[0230] For Case 2, the comparison of computational resource consumption between the SETI-DP solver and the traditional method is shown in Table 2.

[0231] Table 2

[0232]

[0233]

[0234] Table 2 shows that the SETI-DP solver is significantly more efficient than traditional solvers. Specifically, while also increasing memory consumption by only a small amount, the SETI-DP solver is over 15 times more efficient than traditional solvers.

[0235] Case 3: Large-scale realistic geological model.

[0236] The specific geometric information for this case is in Figure 11 The material parameters are as follows: p1 =2600m / s,v s1 =1700m / s, ρ1=2200kg / m^3, v p2 =3000m / s, v s2 =1900m / s, ρ2=2500kg / m^3, v p3 =2500m / s,v s1=1850m / s, ρ1=3700kg / m^3, the parameters of the bipolar perfect matching layer complex stretching function are p d |1=3,R0|1=0.0005,α0|1=πf0,p α |1=1,β0|1=1,p β |1=2, d0|1=439, L=300m; P d |2=3,R0|2=0.0005,α0|2=πf0 / 5,p α |2=1,β0|2=1,p β |2=2,d0|2=22.

[0237] from Figure 12 、 Figure 13 It can be seen that when the seismic wave forward simulation in the real geological model is carried out, the multi-pole perfect matching layer also shows a strong absorption performance. The seismic wave is well absorbed at the boundary and is absorbed by the Figure 14 Compared with the traditional method, the SETI-DP solver shows a similar accuracy, and the L-2 norm error is only calculated to be 1.2770×10 -10 .

[0238] For Case 3, the comparison of computational resource consumption between the SETI-DP solver and the traditional method is shown in Table 3.

[0239] Table 3

[0240]

[0241] Table 3 shows that the SETI-DP solver is significantly more efficient than traditional solvers in large-scale, realistic geological models. Specifically, SETI-DP is over 15 times more efficient than traditional solvers, again with only a small increase in memory consumption. Even with higher interpolation orders than traditional methods, SETI-DP maintains higher computational efficiency.

[0242] Based on this, an embodiment of the present invention provides a method for decomposing spectral elements and splicing regions embedded in a multipole perfectly matched layer technology. The multipole perfectly matched layer adopts the concept of an equivalent perfectly matched layer. By extracting the complex stretching function that determines the characteristics of the perfectly matched layer directly from the outside of the partial derivative to the inside of the partial derivative and combining it with the field variables, the complex convolution operation faced by the multipole perfectly matched layer in the time domain is avoided. The field variables have recursive characteristics under different stretching orders. By introducing appropriate auxiliary variables, a universal first-order recursive auxiliary differential equation applicable to multipole perfectly matched layers of any order is constructed, which can be solved by a simple difference strategy. In addition, by introducing Lagrangian multiplier terms on the subdomain interface including the perfectly matched layer subdomain and implicitly enforcing its continuity on the interface, in addition, in order to eliminate the non-uniqueness of the shared Lagrangian multiplier solution at the nodes shared by multiple subdomains, the concept of corner points is introduced at these nodes. By explicitly enforcing the continuity of the original field variables, a unified two-dimensional global interface problem is constructed. This behavior will not have any impact on the system equations of the subdomain, so that the mass matrix of each subdomain still maintains its diagonal characteristics, which makes the solution of the unknown quantities in each subdomain only take a very short time, retaining the huge advantage of explicit calculation of the spectral element method. At each moment, it is only necessary to construct and solve the global problem at the current moment to obtain the solution of the dual variables and the original variables on the interface at the current moment, and then distribute them to the corresponding subdomain to recover the remaining field variables in the subdomain. Compared with traditional methods, not only does it provide a parallel computing format in the computational domain, but most importantly, the correction process of the field variables in the perfectly matched layer, which is the most time-consuming part of the entire calculation process, can also be executed in parallel. This greatly reduces the cost of boundary truncation and makes it possible to apply multi-pole perfectly matched layers to actual geological surveys.

[0243] Example 2

[0244] In this embodiment, a computer terminal device is provided, including:

[0245] one or more processors;

[0246] a memory, coupled to the processor, for storing one or more programs;

[0247] When the one or more programs are executed by the one or more processors, the one or more processors implement the methods in the above embodiments.

[0248] In this embodiment, a computer-readable storage medium is further provided, on which a computer program is stored. When the computer program is executed by a processor, the method in the above embodiment is implemented.

[0249] In this embodiment, an electronic device is further provided, including a memory and a processor. The memory stores a computer program, and the processor is configured to run the computer program to execute the method in the above embodiment.

[0250] The above program can be run in the processor, or it can be stored in the memory (or computer-readable medium), which includes permanent and non-permanent, removable and non-removable media and can be implemented by any method or technology to store information. The information can be computer-readable instructions, data structures, program modules or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technology, read-only compact disc read-only memory (CD-ROM), digital versatile disc (DVD) or other optical storage, magnetic cassettes, tape disk storage or other magnetic storage devices or any other non-transmission media that can be used to store information that can be accessed by a computing device.

[0251] These computer programs can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for executing on the computer or other programmable device to implement the process. Figure 1 a process or multiple processes and / or boxes Figure 1 The steps of the functions specified in one or more blocks can be implemented by different modules corresponding to different steps.

[0252] This embodiment provides such a device or system. The system is called a spectrum element tearing and splicing region decomposition system embedded with multipole perfectly matched layer technology, and includes:

[0253] The geological modeling module is used to construct the 3D geological model required for simulation based on the real geological characteristics, expand the multipolar perfectly matched layer region in the normal direction of the truncation boundary, and generate a grid file containing subdomain material properties, boundary type and grid information;

[0254] The simulation configuration module is used to define the number of subdomains, material parameters, global time step, total duration, source and receiver locations, and generate configuration files to control the simulation process;

[0255] The load distribution module is used to dynamically distribute process tasks to balance the load based on the computational domain and perfectly matched layer area marked in the configuration file, the number of computing units and the field quantity to correct the time consumption;

[0256] The mapping relationship building module is used to divide the subdomain nodes into internal nodes, interface nodes and corner nodes, and establish the Boolean matrix mapping relationship between the subdomain entities and the global system;

[0257] Matrix assembly module, used to assemble the mass matrix and stiffness matrix for subdomains based on the weak form of the elastic wave equation. For the perfectly matched layer subdomain, only the mass matrix is ​​assembled, and its stiffness matrix is ​​loaded to the right-hand side through the vector generated by the recursive auxiliary differential equation;

[0258] The field variable iteration module is used to initialize the displacement, velocity, acceleration field variables and load vectors, update the field variables through time iteration, and correct the displacement and stress in real time by combining the multipole perfectly matched layer recursive auxiliary differential equation;

[0259] The global solution module is used to construct and solve the global interface problem, assign interface dual variables and corner point primitive variables to subdomains, and complete the wavefield evolution calculation.

[0260] As an implementation method in this embodiment, the geological modeling module includes:

[0261] A model building unit, used to generate a three-dimensional geological model based on real geological conditions;

[0262] Perfectly matched layer extension unit, used to extend the multipolar perfectly matched layer region in the normal direction of the model boundary;

[0263] The meshing unit is used to perform hexahedral meshing on the computational domain and the perfectly matched layer region, and generate a mesh file containing node coordinates and element topology relationships.

[0264] As an implementation method of this embodiment, the simulation configuration module includes:

[0265] Parameter definition unit, used to set subdomain elastic parameters, density and boundary type;

[0266] The time-space configuration unit is used to specify the global simulation time step, total duration and source waveform form;

[0267] The receiving point marking unit is used to mark the receiving point position to record the field variable waveform data.

[0268] As an implementation method of this embodiment, the mapping relationship building module includes:

[0269] An interface mapping unit, used for establishing a mapping between subdomain interface dual variables and global dual variables through a signed Boolean matrix;

[0270] A corner point mapping unit, used for establishing a mapping between subdomain corner points and global corner points through an unsigned Boolean matrix;

[0271] Entity association unit is used to establish the correspondence between local and global numbering of nodes, faces and units within the subdomain.

[0272] As an implementation method of this embodiment, the matrix assembly module includes:

[0273] Computational domain matrix unit, used to assemble the mass matrix, stiffness matrix and load matrix of the computational domain subdomain in parallel based on Gauss-Lobatto-Legendre interpolation points;

[0274] The perfectly matched layer correction unit is used to assemble only the mass matrix of the perfectly matched layer subdomain and generate the stiffness correction vector through the recursive auxiliary differential equation to load it into the right-hand side term.

[0275] As an implementation method in this embodiment, the field variable iteration module includes:

[0276] Prediction update unit, used to update displacement and velocity intermediate variables based on Newmark-beta prediction formula;

[0277] The perfectly matched layer correction unit is used to input the displacement variables and historical auxiliary variables into the recursive equation to generate the corrected displacement, stress and load vectors;

[0278] The variable recovery unit is used to recover the acceleration and velocity variables in combination with the Newmark-beta correction formula.

[0279] As an implementation method of this embodiment, the global solution module includes:

[0280] An interface problem building unit, used to form a global two-dimensional interface linear equation system according to the mapping relationship;

[0281] A variable allocation unit is used to allocate the obtained interface dual variables and corner point original variables to the corresponding subdomains;

[0282] The result output unit is used to associate the field variables with the spatial coordinates and write them to the disk, interpolate at the receiving point to generate the simulation waveform, and calculate the L2 norm error to verify the result.

[0283] The system or device is used to implement the functions of the method in the above-mentioned embodiment. Each module in the system or device corresponds to each step in the method, which has been explained in the method and will not be repeated here.

[0284] Through the above implementation, the problems of spectral element tearing and splicing region decomposition in the multipole perfectly matched layer technology embedded in the related art are solved, thereby ensuring that the defects existing in the existing technology are solved.

[0285] The above are merely preferred embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be included in the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. A spectrum element tearing and splicing region decomposition method embedded in multipole perfectly matched layer technology, characterized by: The following steps are involved: Construct the geological model required for simulation based on the actual geological characteristics, add the perfectly matched layer area in the normal direction of the truncation boundary and perform mesh division to generate a grid file containing grid information; Set up the simulation configuration file to define the number of subdomains, material properties, boundary types, time steps, and source and receiver locations; Distribute the process load between the computational domain and the perfectly matched layer area according to the configuration file to balance the workload of each process; Construct the mapping relationship between all subdomains and the global system, divide the subdomain nodes into internal nodes, interface nodes and corner nodes, and establish point, surface and unit mapping between subdomain entities; Based on the weak form of the second-order elastic wave equation, the mass matrix and stiffness matrix are assembled for the subdomains in parallel, and only the mass matrix is ​​assembled for the perfectly matched layer subdomain, and its stiffness matrix is ​​loaded into the right-hand side in vector form through a recursive auxiliary differential equation. Initialize displacement, velocity, acceleration field variables and load vectors, update field variables through time iteration, and modify displacement and stress in real time by combining perfectly matched layer recursive auxiliary differential equations; Solve the global interface problem, assign the interface dual variables and corner point primitive variables to the subdomains, and complete the wavefield evolution calculation.

2. The method according to claim 1, characterized in that Generate a 3D geological model based on real geological conditions and expand the multipolar perfectly matched layer region in the normal direction of the model boundary; The computational domain and the perfectly matched layer region are divided into multiple subdomains. Hexahedral meshing is performed on each subdomain to generate a mesh file containing the subdomain material properties, boundary types, node coordinates, and element topology relationships.

3. The method according to claim 1, characterized in that The process of setting up the simulation profile includes: Define the grid file storage path, number of subdomains, elastic parameters and density of each subdomain; Specify the global simulation time step, total duration, spatial location of the earthquake source and waveform form, and mark the receiving point location to record the field variable waveform; Mark the boundary type between the computational domain and the perfectly matched layer region and generate a configuration file to control the simulation process.

4. The method according to claim 1, wherein The process of constructing the mapping relationship between the subdomain and the global system includes: Define a signed Boolean matrix for each subdomain, map the subdomain interface dual variables to the global dual variables, and ensure that the variables on both sides of the same interface have opposite signs; Define an unsigned Boolean matrix for each subdomain, map the subdomain corner points to the global corner points, and eliminate the non-uniqueness of the solutions of shared nodes; Establish mapping relationships between entities within the subdomain, including the correspondence between local and global numbers of nodes, faces, and hexahedral elements.

5. The method according to claim 1, wherein The process of assembling matrices for subdomains in parallel includes: For the computational domain subdomain, the mass matrix, stiffness matrix and load matrix are assembled in parallel based on the Gauss-Lobatto-Legendre interpolation points; For the perfectly matched layer subdomain, only the mass matrix is ​​assembled, and its stiffness matrix is ​​loaded to the right-hand side by a vector generated by a recursive auxiliary differential equation. The recursive equation is constructed based on a high-order complex stretching function, and the modified stress variable participates in the load update in vector form.

6. The method according to claim 1, characterized in that The process of updating field variables through time iteration and combining the perfectly matched layer recursive auxiliary differential equation to correct displacement and stress in real time includes: Based on the Newmark-beta prediction formula, update the intermediate variables of the subdomain displacement and velocity at the current moment, and update the right-hand side of the load; For the perfectly matched layer subdomain, the displacement variables and historical auxiliary variables are input into the correction function, and the corrected displacement, stress and load vectors at the current moment are generated through the recursive auxiliary differential equation. The right-hand side terms are updated and the auxiliary variables are saved. The corrected right-hand side term is loaded into the global interface problem. The interface dual variables and corner point primal variables are solved and then distributed to the subdomains. The acceleration and velocity variables are recovered by combining the Newmark-beta correction formula.

7. The method according to claim 1, characterized in that The process of solving the global interface problem, assigning interface dual variables and corner point primitive variables to subdomains, and completing the wavefield evolution calculation includes: Associate displacement, velocity, and acceleration field variables with their spatial coordinates, and write the field variables to disk at a preset time; Perform spatial interpolation on the field variables at the receiving point position, extract the displacement waveform data and generate the simulation waveform; Calculate the L2 norm error between the simulation waveform and the reference solution to verify the correctness of the simulation results.

8. A spectrum element tearing and splicing region decomposition system embedded with multipole perfectly matched layer technology, characterized by: The system comprises: The geological modeling module is used to construct the 3D geological model required for simulation based on the real geological characteristics, expand the multipolar perfectly matched layer region in the normal direction of the truncation boundary, and generate a grid file containing subdomain material properties, boundary type and grid information; The simulation configuration module is used to define the number of subdomains, material parameters, global time step, total duration, source and receiver locations, and generate configuration files to control the simulation process; The load distribution module is used to dynamically distribute process tasks to balance the load based on the computational domain and perfectly matched layer area marked in the configuration file, the number of computing units and the field quantity to correct the time consumption; The mapping relationship building module is used to divide the subdomain nodes into internal nodes, interface nodes and corner nodes, and establish the Boolean matrix mapping relationship between the subdomain entities and the global system; Matrix assembly module, used to assemble the mass matrix and stiffness matrix for subdomains based on the weak form of the elastic wave equation. For the perfectly matched layer subdomain, only the mass matrix is ​​assembled, and its stiffness matrix is ​​loaded to the right-hand side through the vector generated by the recursive auxiliary differential equation; The field variable iteration module is used to initialize the displacement, velocity, acceleration field variables and load vectors, update the field variables through time iteration, and correct the displacement and stress in real time in combination with the multipole perfectly matched layer correction function; The global solution module is used to construct and solve the global interface problem, assign interface dual variables and corner point primitive variables to subdomains, and complete the wavefield evolution calculation.

9. A computer terminal device, characterized in that: include: one or more processors; a memory, coupled to the processor, for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors implement the spectral element tearing and splicing region decomposition method embedded in the multipole perfectly matched layer technology according to any one of claims 1 to 6.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method for spectral element tearing and splicing region decomposition of the embedded multipole perfectly matched layer technology according to any one of claims 1 to 6 is implemented.

Images