A method and system for dynamically updating a preconditioned GMRES site seismic simulation

Abstract

The application relates to the field of earthquake engineering and high-performance computing, and discloses a dynamic updating pre-processing GMRES site seismic simulation method and system, which comprises the following steps: a three-dimensional complex site model is established, boundary grid division and effectiveness checking are implemented; a free field is calculated, a scattering field is expressed based on single-layer potential theory, and a boundary integral equation is established; a boundary element topology structure is established, non-overlapping and overlapping subdomain division is carried out, and local boundary integral coefficient matrixes of each overlapping subdomain are generated in parallel to perform LU decomposition; the GMRES iteration solving of the matrix free restart is cooperatively executed by CPU and GPU, the matrix-vector multiplication, local preprocessing and Arnoldi orthogonalization are completed; convergence is judged according to the real relative residual, and convergence stagnation is judged according to the restart cycle residual ratio; the scattering field is calculated, is superimposed with the free field, and three-dimensional complex site seismic response is obtained. The application can reduce storage and calculation overhead, improve parallel solving efficiency and convergence stability.

Description

Technical Field

[0001] This invention relates to the fields of earthquake engineering and high-performance computing technology, specifically to a method and system for dynamically updating and preprocessing GMRES site earthquake simulation. Background Technology

[0002] Non-uniform, irregular, and complex three-dimensional sites such as mountains and valleys can induce seismic wave scattering, focusing, and local amplification effects. Existing research indicates that domain discretization methods (such as the finite element method and the finite difference method) have inherent limitations in handling such problems, including large computational scale and complex boundary treatment. The traditional indirect boundary element method (IBEM), as a boundary discretization method, only requires discretization of the site boundaries and automatically satisfies the infinite domain radiation condition, thus possessing unique advantages in simulating seismic motion in any complex three-dimensional site. However, the coefficient matrix generated by IBEM is dense and asymmetric, and the storage and computational costs increase exponentially with the number of elements, making it difficult to balance accuracy, efficiency, and memory usage in solving seismic wave fields in large-scale complex three-dimensional sites. To alleviate these bottlenecks, current research attempts to introduce GPU parallel technology to accelerate the main time-consuming processes such as Green's function calculation and large-scale matrix-vector multiplication. However, if the strategy of assembling and storing complete dense matrices is still used, it fails to fundamentally overcome the dual limitations of memory and solution efficiency. On the other hand, while the Generalized Minimal Residual (GMRES) method avoids the LU decomposition of the complete coefficient matrix, its convergence performance is highly dependent on frequency and grid size. When only a local region consistently generates large residuals, the original local solution domain struggles to effectively characterize the coupling effect between the high-residual subdomain and its surrounding subdomains, thus weakening the local preprocessor's ability to reduce the residuals of the high-residual subdomain and leading to a slower iterative convergence speed. Furthermore, global preprocessor reconstruction introduces a large amount of unnecessary local matrix calculations and LU decomposition overhead, increasing the overall computational cost. Therefore, this invention provides a method and system for dynamically updating the preprocessed GMRES site seismic simulation, solving the above technical problems. Summary of the Invention

[0003] This invention provides a dynamic update preprocessing GMRES site seismic simulation method and system, which solves the problems of large generation and storage overhead of global dense coefficient matrix, large computational load of GMRES iteration and easy convergence stagnation in seismic response simulation of arbitrarily large-scale three-dimensional complex sites. It can reduce storage and computation overhead and improve parallel solution efficiency and convergence stability.

[0004] Therefore, the present invention provides the following technical solution: A method for dynamically updating preprocessed GMRES site seismic simulation, the method comprising: Step 1: Establish a three-dimensional numerical model of the complex site on the CPU and perform boundary mesh generation and validity check. Step 2: Calculate the free field on the GPU, represent the scattering field based on the single-layer potential theory, and establish the boundary integral equation; Step 3: On the CPU side, establish the boundary unit topology and divide the non-overlapping and overlapping subdomains, and on the GPU side, generate the local boundary integral coefficient matrix of each overlapping subdomain in parallel and perform LU decomposition. Step 4: The generalized minimum residual method GMRES is used to solve the matrix free restart problem in a collaborative manner between the CPU and GPU. Matrix-vector multiplication, local preprocessing, and Arnoldi orthogonalization are performed in parallel on the GPU. Step 5: Determine convergence based on the true relative residuals, and determine convergence stall based on the residual ratio of the restart cycle. When convergence stall occurs, perform dynamic updates of local preprocessed data for the high residual subdomain. Step 6: After convergence, calculate the scattered field and superimpose it with the free field to obtain the three-dimensional seismic response of the complex site.

[0005] Optionally, in step 1, firstly, the three-dimensional complex site numerical model is constructed on the CPU, and the medium parameters and incident seismic wave information corresponding to the three-dimensional complex site numerical model are obtained. The medium parameters include density, shear wave velocity and Poisson's ratio. Subsequently, the CPU performs boundary mesh generation on the three-dimensional complex site numerical model and performs a validity check on the generated boundary mesh. The CPU then transmits the medium parameters, incident seismic wave information, and boundary mesh data corresponding to the three-dimensional complex site numerical model to the GPU. The boundary mesh data includes node coordinates, center coordinates, normal vectors, area, and node connection information. When performing validity checks on the partitioned boundary mesh, the highest computation frequency of the current computational condition is used as the basis. f max and the minimum wave velocity in the numerical model of the three-dimensional complex site c min Determine the shortest wavelength l min : , The maximum allowable feature size of the boundary unit is 1 / 10 of the shortest wavelength. h allow : , No. i Each boundary cell has a characteristic size based on its maximum side length. h i and check h i Is it less than h allow ;like hi > h allow If this condition is met, then the boundary unit is locally encrypted until the condition is satisfied. h i ≤ h allow After completing the dimensional checks of all boundary elements, proceed to step 2.

[0006] Optionally, in step 2, the GPU receives the parameter information of the three-dimensional complex site numerical model transmitted by the CPU and calculates the free field displacement and free field stress. Then, based on the single-layer potential theory in the Indirect Boundary Element Method (IBEM), the scattered field is expressed as the integral form of the product of the virtual wave source density and the Green's function on the boundary of the three-dimensional complex site numerical model. The boundary integral equation is established according to the boundary conditions, and after numerical discretization, a complex linear equation system with the virtual wave source density as the unknown is formed. The scattered field displacement and scattered field stress are characterized by the boundary integrals of the virtual wave source density and the displacement Green's function and the stress Green's function, respectively. When constructing the complex linear equation system based on the boundary conditions, the load term on the right-hand side of the complex linear equation system is determined by the negative values ​​of the free field displacement and the free field stress. When the k-th field element does not coincide with the l-th source element (i.e., k ≠ l), the GPU uses Gaussian integration to calculate the scattering field displacement and stress Green's function; when they coincide (i.e., k = l), analytical solution is used to form the self-acting boundary integral coefficients.

[0007] Optionally, in step 3, the CPU determines the adjacency relationship between boundary units based on the node connection information of the boundary units, constructs the boundary unit topology relationship, and establishes a global topology structure. Based on the global topology structure, the overall boundary units are divided into several non-overlapping subdomains, and then each non-overlapping subdomain is expanded outward based on the boundary unit topology relationship to form overlapping subdomains. The GPU generates a local coefficient matrix for each overlapping subdomain and performs LU decomposition on the local coefficient matrix of each overlapping subdomain in parallel to obtain the initial local preprocessor corresponding to each overlapping subdomain. When the GPU generates a local coefficient matrix for each overlapping subdomain, the first... s The non-overlapping subfields are denoted as C s non-overlapping subfields C s Extend outward along the topological connection relationship of the boundary unit Layers, the resulting overlapping subdomains are denoted as Oh s ( ) When there is no need to specify the exact number of layers, it is simply referred to as Ohs GPU side Oh s ( ) The boundary within the region is used as the integration domain to generate the corresponding local coefficient matrix.

[0008] Optionally, step 4 includes: Step 41: After entering the k-th GMRES restart cycle, the CPU is responsible for restart cycle management and determining the k-th restart cycle. k The virtual source density iteration value at the start of each GMRES restart cycle, and the preprocessor used in that cycle is determined and invoked. M k On the GPU side, a matrix-free approach is used to calculate the true residual at the start of the cycle based on the current virtual source density iteration value. r k,start And generate initial Krylov basis vectors. v k,1 ; Step 42: The CPU and GPU collaborate to execute the preprocessing GMRES iteration within the current cycle to solve for the virtual source density. On the GPU side, the saved local LU decomposition results and row permutation information are first called to perform preprocessing operations. Then, within the current restart cycle, the product of the global coefficient matrix and the current Krylov basis vector is calculated in a matrix-free manner, and subsequent Krylov basis vectors are generated successively using the Arnoldi orthogonalization process. Finally, the CPU receives the calculation results returned by the GPU, updates the Hessenberg matrix, solves the small-scale least squares problem, and then updates the virtual source density iteration value.

[0009] Optionally, in step 42, during the preprocessing operation, the first... j Krylov basis vectors in the next iteration v k,j During preprocessing, the GPU performs the following operations in sequence: First, the GPU calls the local preprocessor used in the current restart cycle, utilizing each overlapping subdomain. Oh s ( ) The corresponding local LU decomposition results are used to complete the local solution in parallel; then, only the non-overlapping subdomains in each local solution result are extracted. C s The corresponding components, and the non-overlapping subfields C s The corresponding component combination is the preprocessing vector. z k,jFinally, a matrix-free method is used to perform global coefficient matrix-vector multiplication. The product results are orthogonalized to generate new Krylov basis vectors. v k,j+1 ; During global coefficient matrix-vector multiplication, the GPU selects field boundary cells and source boundary cells block by block. Multiple GPU threads calculate the boundary integral coefficients within each block in parallel, and then perform the multiplication of that coefficient block with the current preprocessed vector. z k,j The local multiplication operation of the corresponding components is performed, and the result is accumulated into the output vector. The memory occupied is released immediately after the calculation of each block is completed. After traversing all blocks, the complete matrix-vector product result can be obtained, avoiding the need to store the complete global coefficient matrix.

[0010] Optionally, step 5 includes: Step 51, within the current GMRES restart cycle, according to the... j The estimated relative residual of the next iteration Determine whether to continue iteration; at the end of the current restart cycle, determine the virtual source density iteration value at the end of the restart cycle. X k The GPU uses a matrix-free approach to calculate the true residual. r k and the corresponding true relative residuals True relative residual For true residuals r k The size and known vector b The ratio of their sizes; according to Determine whether the preset convergence criterion is met. If it is met, proceed to step 6; otherwise, proceed to step 52. Step 52: The CPU calculates the residual ratio for the current GMRES restart cycle based on the actual residuals at the start and end times. r k ;like r k ≥ t 2, or consecutive Each restart cycle meets the requirements. t 1≤ r k < t 2, t 1 and t If 2 is the preset stagnation judgment threshold, then a severe stagnation or persistent stagnation is determined, a high residual subdomain is identified, and the local solution region of the high residual subdomain is dynamically expanded; otherwise, the current preprocessor remains unchanged. Step 53: The CPU determines whether the current number of restart cycles has reached the preset maximum number of restarts. If the convergence condition is not met and the current number of restart cycles has not reached the maximum number of restarts, it enters the next GMRES restart cycle and repeats steps 41, 42, and steps 51, 52. If the convergence condition is still not met after reaching the maximum number of restarts, a non-convergence prompt is output and the current calculation condition ends.

[0011] Optionally, in step 52, the restart cycle residual ratio r k Less than the first threshold t At time 1, the residual decrease is considered normal; it is not less than the first threshold. t 1 and less than the second threshold t At 2 o'clock, suspected convergence stagnation is recorded. If this condition is met for a preset number of consecutive GMRES restart cycles, it is determined to be persistent convergence stagnation; not less than the second threshold. t At time 2, it was determined to be a severe convergence stagnation; When severe or persistent convergence stagnation occurs, the GPU calculates the sum of squared magnitudes of the true residual components in the three directions for each boundary cell, and accumulates them according to the non-overlapping subdomains to which each boundary cell belongs, obtaining the residual amplitude of each non-overlapping subdomain. The GPU then transmits the residual amplitudes of each non-overlapping subdomain to the CPU. Upon receiving the data, the CPU sorts the non-overlapping subdomains in descending order of amplitude and accumulates the residuals of each non-overlapping subdomain sequentially to calculate the proportion of the currently selected non-overlapping subdomain residual to the total sum of the residuals of all non-overlapping subdomains. When this proportion reaches a preset accumulation ratio... c Or the number of selected non-overlapping subdomains has reached the maximum number of updates in a single session. n max When the selection is terminated, the selected non-overlapping subdomains are marked as high residual subdomains; where, n max The total number of all non-overlapping subdomains S The upper limit of the proportion of non-overlapping subdomains in a single update α jointly determined, and α and c All values ​​are within (0,1].

[0012] For any high residual subdomain C s If the current number of overlapping layers If the preset maximum number of overlap layers is not reached, the corresponding overlapping subdomain will be... Oh s ( ) Expand outward by one layer to construct Oh s( +1) If an update is triggered again after the maximum number of overlap layers has been reached, the adjacent non-overlapping subdomain with the largest residual magnitude will be selected from the non-overlapping subdomains that are adjacent to the current solution domain and have not yet been included in the domain. C t , and then C s and C t Merge to form a new base region, and expand layer by layer from the periphery of the base region according to the preset maximum number of overlap layers to reconstruct the local solution region; Expanding the local solution domain does not change the original non-overlapping subdomain partitioning and the corresponding virtual wave source density; on the GPU side, only the local boundary integral coefficient matrix is ​​regenerated for the high residual subdomain and LU decomposition is performed, while the local preprocessor parameters of the other non-overlapping subdomains are reused.

[0013] Optionally, in step 6, the GPU calculates the scattered field displacement in parallel based on the displacement Green's function and the virtual wave source density, and superimposes the scattered field displacement with the free field displacement to obtain the total displacement response; wherein, the displacement calculation of multiple field points or multiple frequency conditions is executed in parallel on the GPU.

[0014] A dynamically updated preprocessed GMRES site seismic simulation system, the system comprising: The 3D complex site model building module establishes a 3D complex site numerical model on the CPU and performs boundary mesh generation and validity checks. The module for establishing free field and boundary integral equations calculates the free field on the GPU, represents the scattering field based on single-layer potential theory, and establishes the boundary integral equations. The boundary unit topology and region partitioning module establishes the boundary unit topology and partitions non-overlapping and overlapping subdomains on the CPU side, and generates the local boundary integral coefficient matrix of each overlapping subdomain and performs LU decomposition in parallel on the GPU side. The CPU-GPU co-processing module performs iterative solution of the generalized minimum residual method GMRES with matrix free restart through CPU and GPU collaboration, and performs matrix-vector multiplication, local preprocessing and Arnoldi orthogonalization in parallel on the GPU. The module for judging the true residual and convergence determines convergence based on the true relative residual and determines convergence stall based on the residual ratio of the restart cycle. When convergence stall occurs, it performs dynamic updates of local preprocessed data on the high residual subdomain. The three-dimensional complex site seismic response calculation module calculates the scattered field after convergence and superimposes it with the free field to obtain the three-dimensional complex site seismic response.

[0015] The present invention provides a dynamically updated preprocessing method for GMRES site seismic simulation. Using IBEM as the physical solution thread, the entire solution process avoids explicit assembly and storage of the complete global dense coefficient matrix. Instead, it divides boundary cells into several sub-blocks, calculates the integral coefficients of each boundary cell in parallel on the GPU, and directly performs the product operation between the global coefficient matrix and vectors. Based on this, a CPU-GPU collaborative computing framework is constructed to execute the preprocessing-based restarted GMRES solution. When iterative convergence is detected to be persistently stalled or severely slowed, the algorithm further identifies high residual subdomains, dynamically adjusts their corresponding local solution regions, and updates the corresponding local preprocessor parameters, thereby significantly reducing storage requirements while effectively improving overall solution efficiency. Compared with existing technologies, the present invention has the following technical advantages: (1) Before solving, the boundary grid is checked for the shortest wavelength discretization requirement based on the highest calculation frequency and the minimum wave velocity. For grids that do not meet the requirements, the boundary element is locally refined until the requirements are met, thereby ensuring the effectiveness of the boundary grid and improving the robustness and engineering applicability of the method.

[0016] (2) The matrix free method is adopted to generate and calculate the boundary integral coefficients in parallel block by block on the GPU. After use, the coefficients are released without forming and storing the complete global dense coefficient matrix, so as to reduce the storage pressure and improve the computational efficiency in the process of solving large-scale three-dimensional complex sites.

[0017] (3) Through the collaborative division of labor between the CPU and GPU, the CPU completes mesh checking, topology establishment, subdomain partitioning and convergence status judgment, while the GPU completes local matrix generation, LU decomposition, matrix free calculation and large-scale vector operation, thereby giving full play to the respective advantages of the CPU in logic control and the GPU in parallel computing, and significantly improving the overall computing efficiency.

[0018] (4) The convergence state is determined based on the true relative residuals, and persistent or severe convergence stagnation is identified based on the residual ratio of the restart cycle. Only the local preprocessed data corresponding to the high residual subdomains are updated. The expansion of the local solution region is preferentially achieved by increasing the number of overlapping layers. After reaching the preset upper limit of the number of layers, adjacent non-overlapping subdomains with large residuals that can actually expand the region are introduced. This achieves accurate local update positioning and hierarchical execution of the expansion strategy, thereby effectively reducing unnecessary computational overhead and improving convergence stability and solution efficiency. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly described below. Obviously, the drawings described below are merely some embodiments of the present invention, and those skilled in the art can obtain other drawings based on these drawings without creative effort.

[0020] Figure 1 A flowchart illustrating a method for dynamically updating and preprocessing GMRES site seismic simulation according to a specific embodiment of the present invention; Figure 2 This is a schematic diagram of a three-dimensional sedimentary basin seismic response calculation model in an embodiment of the present invention; Figure 3 This is a schematic diagram of the physical solution process for the boundary integral equation of seismic response in complex sites based on IBEM in an embodiment of the present invention. Figure 4 This is a schematic diagram of the process of freely restarting GMRES using the CPU-GPU collaborative execution matrix in an embodiment of the present invention; Figure 5 This is a schematic diagram illustrating the construction of non-overlapping and overlapping subdomains in an embodiment of the present invention; Figure 6 This is a schematic diagram of adaptive update of high residual subdomains in an embodiment of the present invention; Figure 7 This is a schematic diagram of the structure of a dynamically updated preprocessing GMRES site seismic simulation system provided in a specific embodiment of the present invention. Detailed Implementation

[0021] The specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0022] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0023] Figure 1 A flowchart of a method for dynamically updating and preprocessing GMRES site seismic simulation according to a specific embodiment of the present invention is provided, the method comprising: Step S1: Establish a three-dimensional numerical model of the complex site and perform boundary mesh generation and validity check.

[0024] First, the CPU is responsible for constructing a three-dimensional numerical model of the complex site, acquiring the corresponding medium parameters and incident seismic wave information. These medium parameters include elastic parameters such as density, shear wave velocity, and Poisson's ratio. This embodiment uses an irregular three-dimensional sedimentary basin as an example. Figure 2 As shown, the computational domain of the model includes an outer soil body and an irregular three-dimensional layered sedimentary basin located inside it, which is composed of multiple sedimentary strata.

[0025] Subsequently, the CPU performs boundary mesh generation on the three-dimensional complex site numerical model, performs validity checks on the generated mesh, and transmits relevant parameters and boundary mesh data to the GPU. The boundary mesh data includes node coordinates, center coordinates, normal vectors, area, and node connection information.

[0026] Based on the highest calculation frequency under the current calculation conditions f max and the minimum wave speed in the model c min Determine the shortest wavelength l min :

[0027] The maximum allowable feature size of the boundary unit is 1 / 10 of the shortest wavelength. h allow :

[0028] The validity check refers to, the first i Each boundary cell has a characteristic size based on its maximum side length. h i And check h i ≤ h allow .like h i > h allow If this condition is met, then the boundary unit is locally encrypted until the condition is satisfied. h i ≤ h allow After completing the dimensional checks of all boundary elements, proceed to step S2.

[0029] Step S2: Calculate the free field, represent the scattering field based on single-layer potential theory, and establish the boundary integral equation. Based on the indirect boundary element method (IBEM), virtual wave sources are arranged on the discrete boundary to construct the scattering wave field, and a discrete boundary integral equation with the virtual wave source density as the unknown is established according to the boundary conditions.

[0030] like Figure 3As shown, the physical solution process of the boundary integral equation for the seismic response of a complex site based on IBEM includes: first, solving the free-field seismic response of the complex site; then, applying a virtual wave source to the boundary of the complex site and constructing expressions for the displacement and stress of the scattered field based on the Green's function; establishing the boundary integral equation according to the boundary conditions and performing numerical discretization on the boundary integral equation to form a system of complex linear equations with the virtual wave source density as the unknown; solving the system of complex linear equations to obtain the virtual wave source density, and solving the scattered field seismic response of the complex site based on the virtual wave source density; finally, superimposing the free field and the scattered field to synthesize the total wave field of the complex site.

[0031] The GPU receives model parameter information transmitted from the CPU and calculates free field displacement and stress accordingly. Based on the single-layer potential theory, the scattered field is represented as an integral of the product of the virtual wave source density and the corresponding Green's function on the model boundary. Boundary integral equations are established according to the boundary conditions, and after numerical discretization, a system of complex linear equations with the virtual wave source density as the unknown quantity can be formed.

[0032] The scattered field displacement and stress are characterized by boundary integrals of the virtual source density and the displacement Green's function and stress Green's function, respectively. A system of linear equations is constructed based on the boundary conditions, and the load term on the right-hand side is determined by the negative values ​​of the free field displacement and stress.

[0033] When the k-th field unit is related to the k-th field unit When the source units do not coincide ( On the GPU side, Gaussian integrals are used to calculate the displacement and stress Green's function; when the two coincide ( If the boundary condition is such that the self-acting boundary integral coefficients are obtained by analytical solution, then analytical solution is used.

[0034] In one specific embodiment, based on the single-layer potential theory, in complex site boundaries... C Set virtual wave source density φ j ( y The displacement of the scattered field is expressed as:

[0035] The scattered field stress is expressed as:

[0036] In the formula, the superscript s represents the scattered field. x As the venue, y As the source point, i、j The subscripts for the direction components correspond to... x , y , z Three coordinate directions, φ j ( y ) is a point y Placej Virtual wave source density in the direction, G ij ( x,y ) is the displacement Green's function. T ij ( x,y ) is the stress Green's function. C Represents the boundary of a complex site.

[0037] Based on the boundary conditions, establish a system of boundary integral equations with the virtual wave source density as the unknown. C Taking 1 as an example, its zero-stress boundary condition can be expressed as:

[0038] In the formula, t f i ( x ) indicates field point x Place i The free field boundary stress in the direction is represented by the integral term, which represents the corresponding scattering field boundary stress.

[0039] Surface boundary C 1 Discretized into N For each element, the above integral expression can be further expressed in the following discrete form:

[0040] In the formula, y m Indicates the first m One source point, x n Indicates the first n Each venue N Indicates the number of boundary units.

[0041] It can be seen that the system of linear equations established by the boundary conditions can be simplified to the following form:

[0042] in, A This is the global boundary integral coefficient matrix corresponding to the discrete boundary integral equation; X The vector to be determined is composed of the virtual wave source densities of each boundary element in three directions; b Given a vector, composed of the negative values ​​of free-field displacement and stress. In actual matrix operations, the GPU neither constructs nor stores the complete global dense coefficient matrix. A Instead, it dynamically calculates the required boundary integral coefficients based on the block combination of field point boundary elements and source point boundary elements.

[0043] For the k The field boundary unit and the first Each source boundary element, and the corresponding 3×3 boundary integral coefficient sub-blocks are denoted as... The corresponding calculations are performed by the GPU. At that time, the Gaussian integral is used to calculate the integral of the stress Green's function over the non-self-acting element pair to obtain the non-self-acting coefficient sub-block. ;when In this case, the self-action coefficient sub-block is solved analytically. A kk .

[0044] Step S3: Establish the boundary unit topology, divide the non-overlapping subdomains and construct overlapping subdomains, and further generate the local boundary integral coefficient matrix directly for each overlapping subdomain and perform LU decomposition.

[0045] The CPU determines the adjacency of boundary cells based on their node connection information, constructs the boundary cell topology, and establishes the global topology accordingly. Based on this, the overall boundary cell is divided into several non-overlapping subdomains, and each subdomain is then expanded outwards based on the topological connection relationships to form overlapping subdomains. The GPU generates a local coefficient matrix for each overlapping subdomain and performs parallel LU decomposition on the local coefficient matrices of each subdomain to obtain the initial local preprocessor for each subdomain, ready for use in step S5.

[0046] The first s The non-overlapping subfields are denoted as C s Its topological connection relationship extends outward along the boundary unit. Layers, the resulting overlapping subdomains are denoted as Oh s ( ) When there is no need to specify the exact number of layers, it is simply referred to as Oh s The GPU side directly uses Oh s ( ) The boundary within the region is used as the integration domain to generate the corresponding local coefficient matrix.

[0047] In one embodiment, such as Figure 5 As shown, the first s The non-overlapping subfields are denoted as C s The corresponding overlapping subfield is denoted as Oh s ( ) .in, C sDetermine the scope of the local solution results. Oh s Determine the computational range of the local boundary integral coefficient matrix. Indicates the number of overlapping layers.

[0048] In step S3, the processes of establishing the boundary unit topology, dividing non-overlapping subdomains, and constructing overlapping subdomains include: Step S31: The CPU obtains the set of boundary cells. B Boundary cell node connection information, total number of boundary cells N Number of target subdomain units n c and initial number of overlapping layers 0; Based on the boundary element node connection information, the connection with the boundary element will be... b i The actual adjacent boundary cells form its first-order topological adjacency set. G i And based on each first-order topological adjacency set G i Establish boundary cell topology G .

[0049] Step S32, set the boundary elements B All boundary cells in the set are included in the unassigned boundary cell set. U From set U Select a boundary element i Starting from 0, the boundary units are... i 0 is added to the current non-overlapping subfield. C s and from the set U Remove from the middle; along the boundary cell topological connection relationship, sequentially select the non-overlapping subdomains. C s For each boundary cell that is adjacent to another boundary cell that has not yet been assigned, select that boundary cell and add it to the current non-overlapping subdomain. C s and from the set U Remove from the middle until the current non-overlapping subfield. C s The number of boundary units in the target subdomain reaches the target number of units. n c Or there are no unassigned adjacent boundary units; save the current non-overlapping subdomain. C s Repeat the above process until a set is reached. U Empty.

[0050] Step S33, for boundary elements less than nc / 2 Tail subdomain C s Determine the set of non-overlapping subdomains adjacent to it. N s And calculate the tail subfield. C s With each adjacent non-overlapping subdomain C r Number of topological connections between e sr ; the tail subdomain C s Merge into the adjacent non-overlapping subdomain with the largest number of topological connections. C t Delete the original trailing subdomain C s Then update the topological connectivity relationships between non-overlapping subdomains, the set of adjacent subdomains, the subdomain numbers, and the number of topological connections; repeat the above process until there are no boundary cells with a number less than [a certain value]. n c / 2 The tail subdomains are checked, and each non-overlapping subdomain is checked to ensure that they do not overlap and completely cover all boundary cells.

[0051] Step S34, for each non-overlapping subdomain C s , the initial overlapping subdomain Oh s (0) And the new boundary unit in layer 0 F s (0) All were determined to be C s Expand outwards sequentially according to the topological connection relationship of the boundary units. Floor 0, of which, the first New boundary unit in layer F s ( ) For the first The adjacent cells of the newly added boundary cell in layer -1 have not yet entered the overlapping subdomain. Oh s ( -1) The boundary unit will be the first New boundary unit in layer F s ( ) Add to current overlapping subfield Ohs ( -1) The overlapping subfields are obtained. Oh s ( ) ;Finish After expanding to level 0, each non-overlapping subdomain is obtained. C s Corresponding overlapping subfields Oh s Thus, the boundary unit topology is obtained. G non-overlapping subdomain set { C s} and overlapping subdomain set { Oh s}

[0052] For the s Overlapping subdomains, the GPU executes the local boundary integration coefficient matrix. A s The parallel computation and assembly of this matrix are independent of the global coefficient matrix generation, rather than being obtained by truncating the global matrix. Further integration of the local boundary coefficient matrix... A s Perform LU decomposition, and store the decomposition results and corresponding row permutation information for use in the initial preprocessor construction. After the local solution process is complete, only the non-overlapping subdomains are extracted. C s The corresponding components are then placed in the corresponding degrees of freedom positions within the global vector to obtain the preprocessed vector. The entire process does not construct a complete global preprocessing matrix.

[0053] In step S3, the local boundary integral coefficient matrix generation and LU decomposition process includes: Step S35, for the first s overlapping subdomains Oh s The CPU will overlap subdomains Oh s The boundary cell number in the data is transmitted to the GPU; the GPU reads the overlapping subfield. Oh s The node coordinates, center coordinates, normal vectors, and area of ​​the inner boundary elements are obtained, and the calculation frequency and medium parameters are also acquired.

[0054] Step S36, GPU side targets overlapping subdomains Oh s Each field boundary element within k and each source boundary unit Process them separately. When At that time, the stress Green's function was calculated using Gaussian integral on the GPU. T Integrating over the corresponding non-self-acting boundary element pair yields the first... k The field boundary unit and the first 3×3 boundary integral coefficient block between source boundary elements ;when At that time, the self-acting boundary integral coefficient block is processed in the same way as in step S2. A kk .

[0055] Step S37, overlap subdomains Oh s Boundary integral coefficient blocks corresponding to each field boundary element and source boundary element The local boundary integral coefficient matrix is ​​obtained by assembly. A s The local boundary integral coefficient matrix A s By overlapping subfields Oh s The boundary cells within the region are generated directly, without being truncated from the complete global boundary integral coefficient matrix; non-overlapping subdomains are determined. C s In overlapping subdomains Oh s The corresponding position and the corresponding range of degrees of freedom.

[0056] Step S38, the GPU integrates the coefficient matrix of the local boundary. A s Performing LU decomposition yields the desired result. P s A s = L s U s row permutation matrix P s and local LU decomposition results L s and U s ;save P s , L s , U s And the data required for local solution, without calculating the local boundary integral coefficient matrix. A s The inverse matrix.

[0057] Step S39: Based on the local LU decomposition results corresponding to each overlapping subdomain and the data required for local solution, an initial local preprocessor is formed.M 0, and the local preprocessor used in the first GMRES restart cycle is determined to be... M 1, among which, M 1 = M 0.

[0058] Step S4: After entering the k-th GMRES restart cycle, the CPU is responsible for restart cycle management and determining the k-th restart cycle. k The virtual source density iteration value at the start of each GMRES restart cycle, and the preprocessor used in that cycle is determined and invoked. M k On the GPU side, a matrix-free approach is used to calculate the true residual at the start of the current cycle based on the iterative value of the current virtual wave source density. r k,start And generate initial Krylov basis vectors. v k,1 .

[0059] In one embodiment, before entering the first GMRES restart cycle, an initial value for the initial virtual source density is given. X 0, and the GPU calculates the initial true residual using a matrix-free method:

[0060] In the formula, A Let be the global boundary integral coefficient matrix corresponding to the discrete boundary integral equation. b Let be a known vector determined by the negative values ​​of free field displacement and free field stress. X 0 represents the initial value of the given virtual source density. In actual calculations, the GPU uses a matrix-free computation method. A The product of the virtual source density iteration value does not generate and store the complete global boundary integral coefficient matrix. A .

[0061] No. k Before the start of each GMRES restart cycle, X k-1 As the virtual source density iteration value at the start of the restart cycle, and to calculate the true residual at the start of the cycle: .

[0062] Step S5: The CPU and GPU work together to execute the preprocessing GMRES iteration within the current restart cycle to solve for the virtual wave source density.

[0063] The CPU and GPU collaborate to execute the preprocessing GMRES iteration within the current cycle to solve for the virtual source density. On the GPU, preprocessing is first performed using the saved local LU decomposition results and row permutation information. Then, within the current restart cycle, the product of the global coefficient matrix and the current Krylov basis vectors is calculated using a matrix-free method, and subsequent Krylov basis vectors are generated sequentially using the Arnoldi orthogonalization process. Finally, the CPU receives the calculation results returned by the GPU, updates the Hessenberg matrix, solves a small-scale least squares problem, and thus updates the virtual source density iteration value.

[0064] When performing preprocessing operations, for the first j Krylov basis vectors in the next iteration v k,j During preprocessing, the GPU performs the following operations in sequence: First, the GPU calls the local preprocessor used in the current restart cycle, utilizing each overlapping subdomain. Oh s ( ) The corresponding local LU decomposition results are used to complete the local solution in parallel; then, only the non-overlapping subdomains in each local solution result are extracted. C s The corresponding components are combined into a preprocessed vector. z k,j Finally, a matrix-free method is used to perform global coefficient matrix-vector multiplication. The product results are orthogonalized to generate new Krylov basis vectors. v k,j+1 .

[0065] When performing global coefficient matrix-vector multiplication, the GPU selects field boundary cells and source boundary cells block by block. Multiple GPU threads calculate the boundary integral coefficients within each block in parallel, and then perform the multiplication of that coefficient block with the current preprocessed vector. z k,j The local multiplication operation of the corresponding components is performed, and the result is accumulated into the output vector. After the calculation of each block is completed, its video memory is released immediately. After traversing all blocks, the complete matrix-vector multiplication result can be obtained, thus avoiding the need to store the complete global coefficient matrix.

[0066] In one embodiment, during each GMRES restart cycle, the GPU completes local preprocessing and solves matrix-vector multiplication and Arnoldi large-scale vector operations using a matrix-free approach, while the CPU updates the Hessenberg matrix and solves small-scale least squares problems.

[0067] For the j Krylov basis vectorsv k,j The GPU calls the stored local LU decomposition results in each overlapping subdomain and performs parallel local preprocessing equations to solve the corresponding components in each subdomain, thus obtaining the local preprocessing vectors corresponding to each non-overlapping subdomain. z k,j :

[0068] in, M k This represents a local preprocessor formed by the local LU decomposition results corresponding to each overlapping subdomain within the current restart cycle. The following calculations are then performed to generate the orthogonalization vectors required to generate the next Krylov basis vectors. w k,j :

[0069] in, A This is the global boundary integral coefficient matrix corresponding to the discrete boundary integral equation. In actual computation, the GPU sequentially generates boundary integral coefficient sub-blocks corresponding to each subdomain, multiplies the current sub-block with the input vector and accumulates it to the result vector, and releases the sub-block after the operation is completed. After traversing all sub-blocks, a matrix-vector multiplication operation without having to fully store the global coefficient matrix is ​​achieved.

[0070] Specifically, such as Figure 4 As shown, enter the... k After each GMRES restart cycle, the initialization in step S4 is performed first. The initialization process includes: Step S41, proceed to the next step k After each GMRES restart cycle, initialization is performed, and the iteration sequence number is set. j =1, CPU side obtains the first k Virtual source density iteration value at the start of each GMRES restart cycle X k-1 Initial true residual of the period r k,start Convergence tolerance e tol Given vectors b Current local preprocessor M k and restart length m And fix the local preprocessor used in the current GMRES restart cycle. M k ; Calculate the initial true residual size for the calculation period:

[0071] in, β Represents the initial true residual of the periodr k,start Size, It represents the square root of the sum of the squares of the magnitudes of the vector components.

[0072] when β =0 or β / |b|2≤ e tol season X k = X k-1 And end the current GMRES restart cycle; otherwise, generate the first Krylov basis vector:

[0073] Step S42, initialize the Krylov basis vector set V Preprocessed vector set Z and Hessenberg matrix H The initial value of the minimum estimated relative residual for this period is set to positive infinity, and the optimal iteration step is initialized. and its corresponding coefficient vector Among them, the superscript " "" indicates the optimal iteration step and its corresponding coefficient vector that have been saved in this cycle.

[0074] In step S5, the CPU-GPU collaborative computing process within the current GMRES restart cycle includes: Step S51, in the current GMRES restart cycle... j In this iteration, the GPU calls the current local preprocessor. M k Utilizing each overlapping subdomain Oh s The corresponding local LU decomposition results for Krylov basis vectors v k,j In each overlapping subdomain Ω s Local solutions are performed on the components within:

[0075] in, y s For the first s The intermediate vector obtained by performing forward substitution on the overlapping subdomains q s To further execute the local solution obtained from back-substitution; the vertical lines and their subscripts indicate the extraction of the vector components corresponding to the respective subdomains. v k , j | ΩsRepresents the global Krylov basis vectors v k , j Extracting overlapping subdomains Oh s The corresponding components, Indicates the result from the local solution. q s Extracting non-overlapping subfields C s The corresponding components.

[0076] Extract only the non-overlapping subdomain C from each local solution result s Corresponding components:

[0077] Each non-overlapping subdomain C s The corresponding components form the preprocessing vector. z k,j The overall process of local solution, non-overlapping subdomain component extraction, and combination described above can be represented as follows:

[0078] Step S52: The GPU calculates the global boundary integral coefficient matrix and the preprocessing vector using a matrix-free method. z k,j The product of the two elements is then generated. On the GPU side, field boundary cells and source boundary cells are selected block by block to directly generate the corresponding boundary integral coefficient blocks. These coefficient blocks are then multiplied by the preprocessed vector. z k,j The corresponding components are multiplied and accumulated into the output vector; after completing the multiplication of the current boundary integral coefficient block, the coefficient block is immediately released, and after traversing all blocks, the vector to be orthogonalized is obtained. w k,j The matrix-vector multiplication described above using the matrix-free method is mathematically represented as:

[0079] Step S53, GPU-side vector to be orthogonalized w k,j Perform the first Arnoldi orthogonalization:

[0080] Then a second Arnoldi orthogonalization is performed:

[0081] in, i The indices of the existing Krylov basis vectors participating in the current orthogonalization, with values ​​ranging from 1 to ... jThe superscript H indicates conjugate transpose. h ij The projection coefficients obtained from the first orthogonalization are Δ. h ij This is the projection correction amount obtained from the second orthogonalization, after correction. h ij The Hessenberg matrix is ​​the first... i Line number j Column elements.

[0082] Calculate the orthogonalized vector to be orthogonalized w k,j Size:

[0083] in, h j+1 , j The vector to be orthogonalized after orthogonalization w k , j The size of the Hessenberg matrix, and also the number of... j The second diagonal element of the column.

[0084] when h j+1 , j If the value is not less than the preset minimum, generate the next Krylov basis vector:

[0085] The preset minimum value is used to determine whether the Arnoldi process should terminate and to avoid... h j+1 , j Normalize when the value is too small; v k , j+1 This is the next Krylov basis vector obtained after normalization.

[0086] Otherwise, record the Arnoldi process termination flag and do not generate the next Krylov basis vector. v k,j+1 .

[0087] Step S54: The CPU receives the calculation results returned by the GPU, updates the Hessenberg matrix, and solves the small-scale least squares problem in the current Krylov subspace to obtain the result. The coefficient vector that achieves the minimum value y (j) Further calculate the current estimated relative residuals:

[0088] in, j For the first j The iterations formed ( j +1)× j Hessenberg matrix of order 1; e 1 represents the first unit vector of the corresponding dimension; the estimated relative residual is used for iterative judgment within the current GMRES restart cycle and does not replace the true relative residual. When the current estimated relative residual is less than the minimum estimated relative residual already saved in this cycle, the current estimated relative residual and iteration step are saved. j and its corresponding coefficient vector y (j) And these are respectively used as the minimum estimated relative residual and the optimal iteration step for this period. and corresponding coefficient vector .

[0089] Step S55, when the current estimated relative residual Not greater than the convergence tolerance At that time, based on the virtual source density iteration value at the beginning of the period. X k-1 , the first to the second j The preprocessing vector and the current coefficient vector obtained in the next iteration y (j) Calculate the current virtual source density iteration value:

[0090] The GPU uses a matrix-free computation method. X k,j The corresponding true residual:

[0091] Further calculation of the true relative residual:

[0092] Step S56, when the true relative residual calculated in step S55 is obtained... e k,j Not greater than the convergence tolerance e tol At that time, X k,j The virtual source density iteration value is determined at the end of the current GMRES restart cycle. X k and output X k The minimum estimated relative residual for this period; if convergence is not determined based on the true relative residual, and if the Arnoldi process termination flag has been recorded, or the current iteration step... j Restart length reached mIf the iteration fails, the current iteration within the GMRES restart cycle ends; otherwise, the next iteration continues.

[0093] Step S57: When the iteration within the current GMRES restart cycle ends and the true relative residual convergence condition is not met, the optimal iteration step saved in this cycle is used. j and its corresponding coefficient vector y Update the virtual source density iteration value:

[0094] The minimum estimated relative residual saved in this period is used as the output value for this period. X k The relative residuals with the minimum estimate for this period.

[0095] Step S6: After the current GMRES restart cycle ends, determine the virtual source density iteration value at the end of the restart cycle. X k The GPU uses a matrix-free approach to calculate the true residual. r k and their corresponding true relative residuals and according to Determine whether the preset convergence criterion is met. If it is met, proceed to step S9; otherwise, proceed to step S7.

[0096] True relative residual The actual residual at the end of the current GMRES restart cycle r k With the known vector b The ratio; when step S6 determines that convergence has not occurred, the restart cycle residual ratio is calculated in step S7. r k The restart cycle residual ratio is the actual residual at the end of the current GMRES restart cycle. r k Compared with the actual residual at the start of the restart cycle r k,start The ratio of their sizes.

[0097] In one embodiment, let the first k The virtual source density iteration value at the end of each GMRES restart cycle is X k The GPU uses a matrix-free approach to calculate the true residual:

[0098] True relative residual e k Defined as:

[0099] GPU-based calculation of the true relative residual CPU side will Compare with the preset convergence tolerance. If If the value is not greater than the preset convergence tolerance, then the convergence condition is satisfied and the process proceeds to step S9 to perform subsequent operations; otherwise, step S7 is executed.

[0100] Step S7: Calculate the restart cycle residual ratio, determine whether persistent or severe convergence stagnation has occurred, and maintain or update the local preprocessor based on the determination result. In the event of stagnation, dynamic updates of the local preprocessed data are performed only on the high residual subdomains to avoid repeatedly performing local matrix generation and LU decomposition on all subdomains.

[0101] The CPU calculates the residual ratio for the current GMRES restart cycle based on the actual residuals at the start and end times. r k .like r k ≥ t 2, or consecutive Each restart cycle meets the requirements. t 1≤ r k < t 2 (of which) t 1. t If 2 is the preset stagnation judgment threshold, then a severe or persistent convergence stagnation is determined, a high residual subdomain is identified, and its local solution region is dynamically expanded; otherwise, the current preprocessor remains unchanged. After completing the above processing, step S8 is executed.

[0102] Restart Cycle Residual Ratio r k < t At time 1, the residual decrease is considered normal; when t 1≤ r k < t At 2 o'clock, a suspected convergence stagnation is recorded. If this condition is met for a preset number of consecutive GMRES restart cycles, it is determined to be a persistent convergence stagnation; when r k ≥ t At time 2, it was determined to be a severe convergence stagnation.

[0103] When a pause occurs, the GPU calculates the sum of the squared magnitudes of the true residual components in the three directions for each boundary cell, and accumulates them according to the non-overlapping subdomains to which each boundary cell belongs, obtaining the residual magnitude of each non-overlapping subdomain. The GPU then transmits the residual magnitudes of each subdomain to the CPU. Upon receiving the data, the CPU sorts the subdomains in descending order of magnitude and accumulates the residuals of each subdomain sequentially to calculate the proportion of the currently selected subdomain residual to the total sum of the residuals of all subdomains. When this proportion reaches a preset accumulation ratio... c Or the number of selected subdomains has reached the maximum number of updates in a single session. n max When this happens, the selection process terminates, and the selected subdomain is marked as a high residual subdomain. Among these, n max The total number of all non-overlapping subdomains S Maximum ratio of subdomains updated per transaction α jointly determined, and α and c All values ​​are within (0,1].

[0104]

[0105] For any high residual subdomain C s If the current number of overlapping layers If the preset maximum number of overlap layers is not reached, the corresponding overlapping subdomain will be... Oh s ( ) Expand outward by one layer to construct Oh s ( +1) If an update is triggered again after the maximum number of overlap layers has been reached, the adjacent non-overlapping subdomain with the largest residual magnitude will be selected from the non-overlapping subdomains that are adjacent to the current solution domain and have not yet been included in the domain. C t , and then C s and C t The regions are merged to form a new base region, and then expanded layer by layer from the periphery of the base region according to the preset maximum number of overlap layers to reconstruct the local solution region.

[0106] Expanding the local solution domain does not change the original non-overlapping subdomain partitioning and its corresponding virtual source density; on the GPU side, only the local boundary integral coefficient matrix is ​​regenerated for the high residual subdomain and LU decomposition is performed, while the local preprocessor parameters of the other subdomains are reused.

[0107] In one embodiment, the CPU calculates the restart cycle residual ratio based on the actual residuals at the beginning and end of the current restart cycle. r k :

[0108] It should be noted that the restart cycle residual ratio represents the degree of decrease in the true residual within a restart cycle and is not replaced by the Arnoldi internal residual estimate.

[0109] Preferably, the convergence tolerance is 10. - ³, First threshold t 1. Set to 0.90, the second threshold. t If the value is 0.98, and the condition of suspected convergence stagnation is met for two consecutive restart cycles, it is determined to be a persistent convergence stagnation.

[0110] If no persistent or severe convergence stagnation occurs, the current local preprocessor remains unchanged, and the next GMRES restart cycle begins. Conversely, if persistent or severe convergence stagnation occurs, the GPU extracts the true residual components r of each boundary unit in the three orthogonal directions. i,x r i,y r i,z And calculate its sum of squares. Specifically, the first... i The sum of the squared moduli of the true residual components of each boundary element in three orthogonal directions. E i Defined as:

[0111] No. s The residual size of each non-overlapping subdomain E s For each boundary unit it contains E i sum:

[0112] CPU side according to E s Sort the data in descending order and select non-overlapping subdomains sequentially. Let the total number of non-overlapping subdomains be . S The maximum percentage of subdomains that can be updated in a single instance is [missing information]. α Then the maximum number of subdomains can be updated in a single operation. n max for:

[0113] For each non-overlapping subdomain selected by the CPU, its subdomain residual is added to the sum of the residuals of the selected subdomains; when this value reaches a preset cumulative proportion as a percentage of the total residuals of all subdomains... c Or the number of selected subdomains has reached the maximum number of subdomains that can be updated in a single update. n max When the selection stops, the already selected non-overlapping subdomains are identified as high residual subdomains. In this embodiment, a preset cumulative ratio is used. c Take 0.60, maximum update ratio α Take 0.20.

[0114] like Figure 6 As shown, high residual subdomain C s When the current overlap layer has not reached the maximum overlap layer, the local solution region is expanded outward by one layer along the topological connectivity. When an update is triggered again after reaching the maximum overlap layer, the adjacent non-overlapping subdomains that can actually expand the local solution region and have the largest subdomain residual size are selected from the adjacent non-overlapping subdomains that are not yet included in the current local solution region. C t ; then, C s and C t The union of the sets is used to form a new base region. Then, the region is expanded layer by layer from the periphery of this base region according to the maximum number of overlap layers to reconstruct the local solution region. Expanding the local solution region only changes the local solution range of the current high residual subdomain, without changing the original non-overlapping subdomain division and its corresponding unknown range.

[0115] On the GPU side, only the local boundary integral coefficient matrix of the high residual subdomain is regenerated and LU decomposition is performed, while other subdomains continue to reuse the original local preprocessed data.

[0116] In step S7, the processes of determining the convergence stagnation state, identifying the high residual subdomain, and updating the local preprocessed data include: Step S71: The CPU obtains the actual residuals at the beginning and end of the current GMRES restart cycle. r k , start and r k , end Boundary unit topology G non-overlapping subdomain set { C s}, Current overlapping subdomain set { Oh s}, Current local preprocessor M k Current number of overlapping layers { s} and maximum number of overlapping layers maxThe set of selected neighboring subdomains for each non-overlapping subdomain, the number of consecutive stalls, and the first threshold. t 1. Second threshold t 2. Number of confirmations q Preset cumulative ratio c And the upper limit α of the proportion of a single update to a subdomain; calculate the residual ratio of the current restart cycle:

[0117] Set the current local preprocessor M k Copy as a local preprocessor for use in the next GMRES restart cycle M k+1 .

[0118] Step S72, when r k Less than the first threshold t When 1, the continuous stagnation count is updated to 0, and the high residual subdomain update is not triggered; when r k Not less than the first threshold t 1 and less than the second threshold t At time 2, the continuous stop count is incremented by 1. The continuous stop count has not reached the number of confirmations. q High residual subdomain updates are not triggered until the required number of confirmations is reached. q When the continuous stall count is updated to 0, a high residual subdomain update is triggered; when r k Not less than the second threshold t At time 2, the continuous stall count is updated to 0, and a high residual subdomain update is immediately triggered. If a high residual subdomain update is not triggered, the local preprocessor is maintained. M k+1 constant.

[0119] Step S73, when the high residual subdomain update is triggered, the GPU side calculates the first... i The sum of squared magnitudes of the true residual components in the three directions of each boundary element is calculated, and this sum of squared magnitudes is determined as the sum of squared magnitudes of the first boundary element. i The residual magnitude of each boundary cell:

[0120] Accumulate according to the non-overlapping subdomains to which each boundary unit belongs, to obtain the first... s Non-overlapping subdomains C s The residual amplitude:

[0121] The GPU transmits the residual magnitudes of each non-overlapping subdomain to the CPU, and the CPU then... Es Sort the non-overlapping subfields from largest to smallest.

[0122] Step S74, the CPU side calculates the total number of non-overlapping subdomains. S and the upper limit of the proportion of subdomains updated in a single instance α Calculate the maximum number of subdomains that can be updated in a single operation:

[0123] The CPU selects non-overlapping subdomains sequentially according to the sorting results and accumulates the sum of squared residual magnitudes of the selected subdomains; when the ratio of the sum of squared residual magnitudes of the selected subdomains to the total sum of residual magnitudes of all subdomains reaches a preset accumulation ratio... c Or the number of selected subdomains has reached the maximum number of subdomains that can be updated in a single update. n max When the selection stops, the previously selected non-overlapping subdomains are identified as high residual subdomains.

[0124] Step S75, for each high residual subdomain C s When its current number of overlapping layers s Less than the maximum number of overlapping layers max At that time, the current number of overlapping layers will be... s Increase by 1, and extend the corresponding local solution domain outward by one layer along the boundary element topological connection relationship to obtain the updated overlapping subdomain. Oh s .

[0125] Step S76, when the high residual subdomain C s Current number of overlapping layers s Maximum number of overlapping layers reached max At that time, determine the set of effective neighboring subdomains that are not yet included in the current local solution domain and can actually expand the local solution domain. N s eff The set of effective adjacent subdomains N s eff When empty, the local preprocessed data corresponding to the current high residual subdomain is not updated; the set of effective neighboring subdomains N s eff When not empty, select the adjacent non-overlapping subdomain with the largest residual magnitude. C t ,Will C t Add to the set of selected neighboring subdomains of the current high residual subdomain; with the current high residual subdomain C s Together with the updated set of selected adjacent subdomains, all non-overlapping subdomains serve as the starting point for region expansion, based on the maximum number of overlap layers. max Reconstructing overlapping subdomains Oh s .

[0126] Step S77, for the corresponding overlapping subdomain Oh s High residual subdomains that have been updated C s The GPU regenerates its corresponding local boundary integral coefficient matrix. A s And perform LU decomposition:

[0127] Save the updated local LU decomposition results and replace the local preprocessor with the updated local preprocessed data. M k+1 The data corresponding to the current high residual subdomain is selected; the remaining subdomains continue to be retained by the local preprocessor. M k The original local preprocessed data was copied.

[0128] Step S78: Output a local preprocessor for use in the next GMRES restart cycle. M k+1 The updated set of overlapping subdomains { Oh s}, Updated current overlap layer number { s The updated set of selected neighboring subdomains for each non-overlapping subdomain and the updated count of consecutive stalls.

[0129] Step S8: The CPU determines whether the current number of restart cycles has reached the preset maximum number of restarts. If the convergence condition is not met and the current number of restart cycles has not reached the maximum number of restarts, the CPU enters the next GMRES restart cycle and repeats steps S4 to S7. If the convergence condition is still not met after reaching the maximum number of restarts, a non-convergence message is output, and the current calculation state ends.

[0130] In step S9, the GPU calculates the scattered field based on the obtained virtual wave source density and superimposes it with the free field to output the three-dimensional complex site seismic response.

[0131] On the GPU, the scattered field displacement is calculated in parallel based on the displacement Green's function and the obtained virtual source density, and then superimposed with the free field displacement to obtain the total displacement response. Displacement calculations for multiple field points or multiple frequency conditions can also be performed in parallel on the GPU.

[0132] In one embodiment, the total displacement response is taken as an example. After the convergence condition is met, the GPU calculates the scattered field displacement in parallel according to the obtained virtual source density, divided into blocks based on the field points. and free field displacement The total displacement response is obtained by superposition:

[0133] In one specific embodiment, the boundary integral coefficients are calculated using single-precision complex numbers, and the field boundary element blocks... I p and source boundary unit blocks J q All sizes are 512; GMRES restart length m The maximum number of restarts is 20, and the convergence tolerance is 30. e tol 10 - ³; Number of target units in non-overlapping subdomains n c The initial number of overlapping layers is 128. 0 is equivalent to 1; first threshold t 1. Second threshold t 2 are 0.90 and 0.98 respectively, continuously q If two restart cycles meet the suspected stagnation condition, it is determined to be a persistent convergence stagnation, with a preset cumulative ratio. c The maximum percentage of subdomains updated in a single transaction is 0.60. α It is 0.20.

[0134] The above parameters can be adjusted based on GPU memory, number of boundary units, computing frequency, and media parameters.

[0135] This method is applicable to any complex three-dimensional terrain, such as irregular mountains, valleys, and canyons.

[0136] Accordingly, such as Figure 7 As shown, the present invention also provides a dynamically updated preprocessing GMRES site seismic simulation system, in which each module sequentially completes its corresponding calculation, and the relevant steps are executed collaboratively by the CPU and GPU. The system includes: The 3D complex site model construction module 701 is used to construct a computational model based on the boundary geometry, medium parameters and incident seismic wave information of the 3D complex site; to obtain the boundary mesh of the 3D complex site, to perform validity checks on the boundary mesh, and to transmit basic parameters and boundary mesh data to the GPU.

[0137] The module 702 for establishing free field and boundary integral equations is used to calculate the displacement and stress of the free field and the displacement and stress of the scattered field. Based on the single-layer potential theory, it constructs the boundary integral equations for characterizing the scattered field by introducing virtual wave source density. The boundary unit topology and region partitioning module 703 is used to construct the boundary unit topology relationship, divide non-overlapping subdomains, and construct overlapping subdomains corresponding to each non-overlapping subdomain based on the topology connection relationship; and directly generate the corresponding boundary integral coefficient matrix based on the local overlapping subdomain information, perform LU decomposition on it, obtain the initial local preprocessor corresponding to each overlapping subdomain, and store the obtained local preprocessing parameters.

[0138] The CPU-GPU co-processing module 704 is used to manage the GMRES restart cycle, determine the local preprocessor used in the current restart cycle, call the local LU decomposition results to perform preprocessing, complete global matrix-vector multiplication and Arnoldi orthogonalization in a matrix-free manner, update the Hessenberg matrix, solve small-scale least squares problems, and update the virtual wave source density iteration value. The true residual and convergence judgment module 705 is used to calculate the true residual and true relative residual at the end of the current restart cycle in a matrix-free manner, and to determine whether the preset convergence condition is met based on the true relative residual; and based on the true relative residual... To determine whether convergence has occurred, the residual ratio during the restart cycle is used as a reference. r k Determine whether persistent or severe convergence stagnation has occurred, and when convergence stagnation occurs, identify the high residual subdomain and expand the corresponding local solution domain; The three-dimensional complex site seismic response calculation module 706 is used to calculate the scattered field based on the obtained virtual wave source density after the convergence condition is met, and then superimpose the scattered field with the free field to obtain the three-dimensional complex site seismic response.

[0139] This invention also provides an electronic device, which includes a processor, a communication interface, a memory, and a bus. The processor, the communication interface, and the memory communicate with each other via the bus. The processor includes a CPU and a GPU, and the memory stores a computer program. When the computer program is executed by the processor, it implements the aforementioned dynamically updated preprocessing GMRES site seismic simulation method. In practical applications, the electronic device can be a computer or a server.

[0140] The present invention also provides a non-transitory computer-readable storage medium storing a computer program; when the computer program is executed by a processor including a CPU and a GPU, it implements the above-described dynamically updated preprocessing GMRES site seismic simulation method. The non-transitory computer-readable storage medium can be a USB flash drive, a portable hard drive, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disk.

[0141] It should be noted that, for the sake of simplicity, the foregoing method embodiments are all described as a series of actions. However, those skilled in the art should understand that the present invention is not limited to the described order of actions, because according to the present invention, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions and modules involved are not necessarily essential to the present invention.

[0142] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, as a computer program product. The computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer program are loaded or executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data provider to another website, computer, server, or data provider via wired or wireless means.

[0143] The embodiments of the present invention have been described in detail above. Specific implementation methods have been used to illustrate the present invention. The descriptions of the embodiments above are only for the purpose of helping to understand the method of the present invention, and are merely some, not all, embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention, and the content of this specification should not be construed as a limitation of the present invention. Therefore, any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for dynamically updating preprocessed GMRES site seismic simulation, characterized in that, The method includes: Step 1: Establish a three-dimensional numerical model of the complex site on the CPU and perform boundary mesh generation and validity check. Step 2: Calculate the free field on the GPU, represent the scattering field based on the single-layer potential theory, and establish the boundary integral equation; Step 3: On the CPU side, establish the boundary unit topology and divide the non-overlapping and overlapping subdomains, and on the GPU side, generate the local boundary integral coefficient matrix of each overlapping subdomain in parallel and perform LU decomposition. Step 4: The generalized minimum residual method GMRES is used to solve the matrix free restart problem in a collaborative manner between the CPU and GPU. Matrix-vector multiplication, local preprocessing, and Arnoldi orthogonalization are performed in parallel on the GPU. Step 5: Determine convergence based on the true relative residuals, and determine convergence stall based on the residual ratio of the restart cycle. When convergence stall occurs, perform dynamic updates of local preprocessed data for the high residual subdomain. Step 6: After convergence, calculate the scattered field and superimpose it with the free field to obtain the three-dimensional seismic response of the complex site.

2. The dynamically updated preprocessing GMRES site seismic simulation method according to claim 1, characterized in that, In step 1, firstly, the three-dimensional complex site numerical model is constructed on the CPU, and the medium parameters and incident seismic wave information corresponding to the three-dimensional complex site numerical model are obtained. The medium parameters include density, shear wave velocity and Poisson's ratio. Subsequently, the CPU performs boundary mesh generation on the three-dimensional complex site numerical model and performs a validity check on the generated boundary mesh. The CPU then transmits the medium parameters, incident seismic wave information, and boundary mesh data corresponding to the three-dimensional complex site numerical model to the GPU. The boundary mesh data includes node coordinates, center coordinates, normal vectors, area, and node connection information. When performing validity checks on the partitioned boundary mesh, the highest computation frequency of the current computational scenario is used as the basis. f max and the minimum wave velocity in the numerical model of the three-dimensional complex site c min Determine the shortest wavelength λ min : ， The maximum allowable feature size of the boundary unit is 1 / 10 of the shortest wavelength. h allow : ， No. i Each boundary cell has a characteristic size based on its maximum side length. h i and check h i Is it less than h allow ;like h i > h allow If this condition is met, then the boundary unit is locally encrypted until the condition is satisfied. h i ≤ h allow After completing the dimensional checks of all boundary elements, proceed to step 2.

3. The dynamically updated preprocessing GMRES site seismic simulation method according to claim 2, characterized in that, In step 2, the GPU receives the parameter information of the three-dimensional complex site numerical model transmitted from the CPU and calculates the free field displacement and free field stress. Then, based on the single-layer potential theory in the Indirect Boundary Element Method (IBEM), the scattered field is expressed as the integral form of the product of the virtual wave source density and the Green's function on the boundary of the three-dimensional complex site numerical model. The boundary integral equation is established according to the boundary conditions, and after numerical discretization, a complex linear equation system with the virtual wave source density as the unknown is formed. The scattered field displacement and scattered field stress are characterized by the boundary integrals of the virtual wave source density and the displacement Green's function and the stress Green's function, respectively. When constructing the complex linear equation system based on the boundary conditions, the load term on the right-hand side of the complex linear equation system is determined by the negative values ​​of the free field displacement and the free field stress. When the k-th field element and the l-th source element do not coincide, and k ≠ l, the GPU uses Gaussian integration to calculate the scattering field displacement and stress Green's function; when they coincide, and k = l, analytical solution is used to form the self-acting boundary integral coefficients.

4. The dynamically updated preprocessing GMRES site seismic simulation method according to claim 1, characterized in that, In step 3, the CPU determines the adjacency relationship between boundary units based on the node connection information of the boundary units, constructs the boundary unit topology relationship, and establishes a global topology structure. Based on the global topology structure, the overall boundary units are divided into several non-overlapping subdomains. Then, based on the boundary unit topology relationship, each non-overlapping subdomain is expanded outward to form overlapping subdomains. The GPU generates a local coefficient matrix for each overlapping subdomain and performs LU decomposition on the local coefficient matrix of each overlapping subdomain in parallel to obtain the initial local preprocessor corresponding to each overlapping subdomain. When the GPU generates a local coefficient matrix for each overlapping subdomain, the first... s The non-overlapping subfields are denoted as C s non-overlapping subfields C s Extend outward along the topological connection relationship of the boundary unit Layers, the resulting overlapping subdomains are denoted as Ω s ( ) When there is no need to specify the exact number of layers, it is simply referred to as Ω s GPU side Ω s ( ) The boundary within the region is used as the integration domain to generate the corresponding local coefficient matrix.

5. The dynamically updated preprocessing GMRES site seismic simulation method according to claim 4, characterized in that, Step 4 includes: Step 41: After entering the k-th GMRES restart cycle, the CPU is responsible for restart cycle management and determining the k-th restart cycle. k The virtual source density iteration value at the start of each GMRES restart cycle, and the preprocessor used in that cycle is determined and invoked. M k On the GPU side, a matrix-free approach is used to calculate the true residual at the start of the cycle based on the current virtual source density iteration value. r k,start And generate initial Krylov basis vectors. v k,1 ; Step 42: The CPU and GPU collaborate to execute the preprocessing GMRES iteration within the current cycle to solve for the virtual source density. On the GPU side, the saved local LU decomposition results and row permutation information are first called to perform preprocessing operations. Then, within the current restart cycle, the product of the global coefficient matrix and the current Krylov basis vector is calculated in a matrix-free manner, and subsequent Krylov basis vectors are generated successively using the Arnoldi orthogonalization process. Finally, the CPU receives the calculation results returned by the GPU, updates the Hessenberg matrix, solves a small-scale least squares problem, and then updates the virtual source density iteration value.

6. The dynamically updated preprocessing GMRES site seismic simulation method according to claim 5, characterized in that, In step 42, during the preprocessing operation, the first... j Krylov basis vectors in the next iteration v k,j During preprocessing, the GPU performs the following operations in sequence: First, the GPU calls the local preprocessor used in the current restart cycle, utilizing each overlapping subdomain. Ω s ( ) The corresponding local LU decomposition results are used to complete the local solution in parallel; then, only the non-overlapping subdomains in each local solution result are extracted. C s The corresponding components, and the non-overlapping subfields C s The corresponding component combination is the preprocessing vector. z k,j Finally, a matrix-free method is used to perform global coefficient matrix-vector multiplication. The product results are orthogonalized to generate new Krylov basis vectors. v k,j+1 ; During global coefficient matrix-vector multiplication, the GPU selects field boundary cells and source boundary cells block by block. Multiple GPU threads calculate the boundary integral coefficients within each block in parallel, and then perform the multiplication of that coefficient block with the current preprocessed vector. z k,j The local multiplication operation of the corresponding components is performed, and the result is accumulated into the output vector. The memory occupied is released immediately after the calculation of each block is completed. After traversing all blocks, the complete matrix-vector multiplication result can be obtained, avoiding the need to store the complete global coefficient matrix.

7. The dynamically updated preprocessing GMRES site seismic simulation method according to claim 5, characterized in that, Step 5 includes: Step 51, within the current GMRES restart cycle, according to the... j The estimated relative residual of the next iteration Determine whether to continue iteration; at the end of the current restart cycle, determine the virtual source density iteration value at the end of the restart cycle. X k The GPU uses a matrix-free method to calculate the true residual. r k and the corresponding true relative residuals True relative residual For true residuals r k The size and known vector b The ratio of their sizes; according to Determine whether the preset convergence criterion is met. If it is met, proceed to step 6; otherwise, proceed to step 52. Step 52: The CPU calculates the residual ratio for the current GMRES restart cycle based on the actual residuals at the start and end times. ρ k ;like ρ k ≥ τ 2, or consecutive Each restart cycle meets the requirements. τ 1≤ ρ k < τ 2, τ 1 and τ If 2 is the preset stagnation judgment threshold, then a severe stagnation or persistent stagnation is determined, a high residual subdomain is identified, and the local solution region of the high residual subdomain is dynamically expanded; otherwise, the current preprocessor remains unchanged. Step 53: The CPU determines whether the current number of restart cycles has reached the preset maximum number of restarts. If the convergence condition is not met and the current number of restart cycles has not reached the maximum number of restarts, it enters the next GMRES restart cycle and repeats steps 41, 42, and steps 51, 52. If the convergence condition is still not met after reaching the maximum number of restarts, a non-convergence prompt is output and the current calculation condition ends.

8. The dynamically updated preprocessing GMRES site seismic simulation method according to claim 7, characterized in that, In step 52, restart cycle residual ratio ρ k Less than the first threshold τ At time 1, the residual decrease is considered normal; it is not less than the first threshold. τ 1 and less than the second threshold τ At 2 o'clock, suspected convergence stagnation is recorded. If this condition is met for a preset number of consecutive GMRES restart cycles, it is determined to be persistent convergence stagnation; not less than the second threshold. τ At time 2, it was determined to be a severe convergence stagnation; When a severe or persistent convergence stall occurs, the GPU calculates the sum of squared magnitudes of the true residual components in the three directions for each boundary cell, and accumulates them according to the non-overlapping subdomains to which each boundary cell belongs, to obtain the residual magnitude of each non-overlapping subdomain. The GPU then transmits the residual magnitude of each non-overlapping subdomain to the CPU. After receiving the residual magnitude, the CPU sorts the non-overlapping subdomains in descending order according to the magnitude, and accumulates the residuals of each non-overlapping subdomain in sequence to calculate the proportion of the residual of the currently selected non-overlapping subdomain to the total residual of all non-overlapping subdomains. When this ratio reaches the preset cumulative ratio γ Or the number of non-overlapping subdomains selected has reached the maximum number of updates in a single session. n max When the selection is terminated, the selected non-overlapping subdomains are marked as high residual subdomains; where, n max The total number of all non-overlapping subdomains S The upper limit of the proportion of non-overlapping subdomains in a single update α jointly determined, and α and γ All values ​​are within (0,1]. For any high residual subdomain C s If the current number of overlapping layers If the preset maximum number of overlap layers is not reached, the corresponding overlapping subdomain will be... Ω s ( ) Expand outward by one layer to construct Ω s ( +1) ; If an update is triggered again after the maximum number of overlap layers has been reached, the neighboring non-overlapping subdomain with the largest residual magnitude will be selected from the non-overlapping subdomains adjacent to the current solution domain that have not yet been included in the current domain. C t , and then C s and C t Merge to form a new base region, and expand layer by layer from the periphery of the base region according to the preset maximum overlap layer to reconstruct the local solution region; Expanding the local solution domain does not change the original non-overlapping subdomain partitioning and the corresponding virtual wave source density; on the GPU side, only the local boundary integral coefficient matrix is ​​regenerated for the high residual subdomain and LU decomposition is performed, while the local preprocessor parameters of the other non-overlapping subdomains are reused.

9. The dynamically updated preprocessing GMRES site seismic simulation method according to claim 7, characterized in that, In step 6, the GPU calculates the scattered field displacement in parallel based on the displacement Green's function and the virtual wave source density, and superimposes the scattered field displacement with the free field displacement to obtain the total displacement response; wherein, the displacement calculation of multiple field points or multiple frequency conditions is executed in parallel on the GPU.

10. A dynamically updated preprocessed GMRES site seismic simulation system, characterized in that, The system includes: The 3D complex site model building module establishes a 3D complex site numerical model on the CPU and performs boundary mesh generation and validity checks. The module for establishing free field and boundary integral equations calculates the free field on the GPU, represents the scattering field based on single-layer potential theory, and establishes the boundary integral equations. The boundary unit topology and region partitioning module establishes the boundary unit topology and partitions non-overlapping and overlapping subdomains on the CPU side, and generates the local boundary integral coefficient matrix of each overlapping subdomain and performs LU decomposition in parallel on the GPU side. The CPU-GPU co-processing module performs iterative solution of the generalized minimum residual method GMRES with matrix free restart through CPU and GPU collaboration, and performs matrix-vector multiplication, local preprocessing and Arnoldi orthogonalization in parallel on the GPU. The module for judging the true residual and convergence determines convergence based on the true relative residual and determines convergence stall based on the residual ratio of the restart cycle. When convergence stall occurs, it performs dynamic updates of local preprocessed data on the high residual subdomain. The three-dimensional complex site seismic response calculation module calculates the scattered field after convergence and superimposes it with the free field to obtain the three-dimensional complex site seismic response.

Images