An explicit finite element parallel computing method for seismic wave simulation

Abstract

The application discloses an explicit finite element parallel computing method for seismic wave simulation, comprising the following steps: constructing a three-dimensional computing domain, generating an unstructured finite element grid based on real terrain and velocity structure data; performing domain decomposition on the computing domain, and dividing the grid into multiple sub-regions; distributing the tasks of the sub-regions to parallel computing nodes, discretely solving a three-dimensional elastic wave equation containing a damping term by using explicit time integration and a lumped mass matrix, synchronously calculating time sequences of wave fields of the sub-regions; integrating the results to obtain seismic wave field distribution, and outputting ground motion parameters of a target region. By constructing an integrated explicit finite element parallel computing framework of a seismic source-medium-terrain, efficient, stable and high-precision simulation of seismic wave propagation under complex geological conditions is realized, and the reliability of strong ground motion prediction is significantly improved.

Description

Technical Field

[0001] This invention relates to the field of earthquake disaster simulation technology, and in particular to an explicit finite element parallel computing method for seismic wave simulation. Background Technology

[0002] Currently, earthquake motion prediction mainly relies on empirical attenuation relationships derived from statistical regression of historical earthquakes, or simplified wave propagation models based on the assumption of homogeneous or layered media. While these methods are widely used in engineering practice, they fail to accurately reflect local topographic undulations, the lateral non-uniformity of underground velocity structures, and the spatial distribution of sedimentary basins. This makes it difficult to accurately characterize the complex spatial distribution patterns and amplitude variations of near-field strong earthquakes, especially in areas with abrupt topographic changes and complex geological structures, where significant deviations often exist between predicted results and actual measurements. In numerical simulation, the finite difference method is widely used due to its relatively simple implementation. However, this method typically relies on structured meshes, often facing challenges in balancing numerical stability and accuracy when dealing with irregular surface morphologies, undulating interfaces, and complex media boundaries. Furthermore, it has poor adaptability to complex geometries. While the spectral element method achieves a better balance between accuracy and efficiency, it has stringent requirements for mesh quality (such as element geometric regularity). Mesh generation is difficult in highly heterogeneous geological models containing numerous discontinuities, resulting in high computational costs and limiting its application in practical large-scale three-dimensional problems.

[0003] In regions with intense tectonic activity, significant topographic differences, and sparse strong-motion observation stations, such as the southern edge of the Tibetan Plateau, the distribution of earthquake hazards is the result of a complex coupling of multiple physical mechanisms, including the directional effect of the hypocenter rupture process, the focusing or shielding effect of topography on seismic waves, and the amplification and resonance effect of sedimentary basins on seismic waves of specific frequencies. Existing prediction models, whether empirical or simplified numerical, struggle to achieve high-fidelity simulations of the entire process from hypocenter rupture to wave propagation and site response under such multi-factor, strongly coupled physical scenarios. Therefore, current technologies cannot provide reliable quantitative data with sufficient spatial resolution for refined seismic hazard analysis in these regions, determination of seismic motion input at major engineering sites, and seismic damage risk assessment of urban building complexes. Summary of the Invention

[0004] The purpose of this invention is to provide an explicit finite element parallel computing method for seismic wave simulation. By constructing an explicit finite element parallel computing framework that integrates the source, medium, and terrain, it achieves efficient, stable, and high-precision numerical simulation of the seismic wave propagation process under complex geological conditions, effectively improving the accuracy and reliability of the prediction of the spatial distribution of strong ground motion.

[0005] To address the aforementioned technical problems, a first aspect of this invention provides an explicit finite element parallel computation method for seismic wave simulation, comprising the following steps: A three-dimensional computational domain for simulating seismic wave propagation is constructed. The three-dimensional computational domain is a physical space covering the target epicenter area. Based on the real terrain data and three-dimensional velocity structure data of the target epicenter area, an unstructured finite element mesh that fits the geological structure is generated. The three-dimensional computational domain is decomposed into multiple sub-regions by dividing the unstructured finite element mesh. The computational task of each sub-region is assigned to multiple parallel computing threads. Each parallel computing thread uses the unstructured finite element mesh, an explicit time integration scheme and a lumped mass matrix to discretize and solve the three-dimensional elastic wave equation containing damping terms, and synchronously calculates the discretized wave field time series of each sub-region. By integrating the synchronous calculation results of all the parallel computing nodes, the seismic wavefield distribution within the three-dimensional computing domain is obtained, and the ground motion parameters of at least one target region within the three-dimensional computing domain are output.

[0006] Furthermore, the generation of an unstructured finite element mesh that conforms to the geological structure based on the real terrain data and three-dimensional velocity structure data of the target epicenter area includes: Acquire high-resolution terrain data and a three-dimensional velocity model covering the target epicenter region; Based on the high-resolution terrain data and the three-dimensional velocity model, a set of spatially discretized grid cells that conforms to the surface topography and the geometry of the underground velocity interface is constructed to form the unstructured finite element mesh.

[0007] Furthermore, the construction of a spatially discretized grid unit set that conforms to the surface topography and the geometry of the subsurface velocity interface, forming the unstructured finite element grid, includes: By integrating the high-resolution terrain data with the three-dimensional velocity model, geometric constraints are established that simultaneously constrain the surface morphology and the interface between the internal media. Based on the geometric constraints, an initial mesh composed of polyhedral elements is generated in the three-dimensional computational domain using an unstructured mesh generation algorithm. The surfaces and edges of the polyhedral elements conform to the geometry of the surface topography and the underground velocity interface. The initial mesh is locally refined, gradually becoming sparser outwards from the epicenter, to obtain an unstructured finite element mesh that conforms to the geological structure.

[0008] Furthermore, the discretization solution of the three-dimensional elastic wave equation containing damping terms using an explicit time integration scheme and a lumped mass matrix to obtain the discretized wave field time series within the three-dimensional computational domain includes: Based on the unstructured finite element mesh, the three-dimensional elastic wave equation containing damping terms is spatially discretized, and the three-dimensional elastic wave equation is transformed into a dynamic system defined on the mesh element nodes; An explicit time-integral scheme with alternating updates of velocity and displacement variables is used to recursively solve the dynamic system in the time domain. In the time-recursive solution process, the lumped mass matrix is ​​used to replace the coordination mass matrix, the mass matrix of the dynamic system is diagonalized, and the solution of the dynamic system is simplified to an explicit node-by-node scalar operation. The discretized wave field time series in the three-dimensional computational domain is generated through the spatial discretization, time recursion, and diagonalization of the mass matrix solution process.

[0009] Further, the step of spatially discretizing the three-dimensional elastic wave equation based on the unstructured finite element mesh, transforming the three-dimensional elastic wave equation into a dynamic system defined on the mesh element nodes, includes: The three-dimensional computational domain is discretized into a set of interconnected polyhedral elements using the unstructured finite element mesh, wherein the vertices of each polyhedral element constitute the mesh element node. Based on the geometric and topological connectivity of the unstructured finite element mesh, shape functions defined on the polyhedral elements are selected to construct a weak form expression of the three-dimensional elastic wave equation at the element level. The three-dimensional elastic wave equation is expressed in a weak form at the element level. Based on the node connection relationship of the unstructured finite element mesh, a global dynamic equation system is assembled with the displacement and velocity at the nodes of the mesh elements as unknown variables, thus completing the transformation of the three-dimensional elastic wave equation into the dynamic system.

[0010] Furthermore, the explicit time integration scheme employing alternating updates of velocity and displacement variables to perform a recursive solution of the dynamic system in the time domain includes: Discretize the simulation time into a series of equally spaced time steps; Within each time step, the grid cell node acceleration of the current time step is calculated using the explicit time integration scheme based on the displacement field and velocity field of the dynamic system in the previous time step. Based on the grid cell node acceleration at the current time step, the grid cell node velocity field at the current time step is updated using the explicit time integration format; Based on the grid cell node velocity field at the current time step, the grid cell node displacement field at the current time step is updated using the explicit time integration format. Following the time step sequence, the steps of calculating acceleration, updating the velocity field, and updating the displacement field are repeated to complete the recursive solution for the entire simulation period, forming the discretized wave field time series in which the displacement and velocity of the grid cell nodes change with time.

[0011] Furthermore, the motion state sequence of several grid cell nodes within the three-dimensional computational domain during the entire simulation period is calculated through the spatial discretization, time recursion, and mass matrix diagonalization process, forming the discretized wavefield time series, including: The displacement and velocity fields of all grid cell nodes in the three-dimensional computational domain at each time step are collected through the spatial discretization, time recursion and mass matrix diagonalization solution process. The displacement field and velocity field at each time step are structured according to the grid cell node number and the time step order to obtain a motion state data set with time as the dimension and spatial nodes as the index. The motion state data set is subjected to consistency verification and format standardization processing, and the output is the discretized wave field time series of the three-dimensional computational domain throughout the entire simulation period.

[0012] Furthermore, the domain decomposition of the three-dimensional computational domain, dividing the unstructured finite element mesh into multiple sub-regions, includes: Based on the element distribution and connection relationship of the unstructured finite element mesh, the three-dimensional computational domain is divided into several spatially continuous and computationally balanced sub-regions. The step of allocating the computational tasks of each sub-region to multiple parallel computation threads includes: Based on the division results of the sub-regions, the calculation tasks of the unstructured finite element mesh set and the associated discretized wave field time series corresponding to each sub-region are dynamically allocated to a corresponding parallel computing thread. The computational tasks assigned to several parallel computing threads are configured to independently perform the recursive solution of the discretized wavefield time series within the sub-regions assigned to them, based on the same explicit time integration format and lumped mass matrix, and to synchronize wavefield variables at the shared boundaries of adjacent sub-regions, thereby achieving collaborative parallel computation of the entire three-dimensional computational domain.

[0013] Furthermore, the computational task allocated to several parallel computing threads is configured to independently perform the recursive solution of the discretized wavefield time series within its assigned sub-region based on the same explicit time integration format and lumped mass matrix, and to synchronize wavefield variables at the shared boundary of adjacent sub-regions, including: Configure each parallel computing thread with the same numerical solver kernel containing the explicit time integration format and lumped mass matrix, and load the unstructured finite element mesh set and initial field information corresponding to the sub-region to which it is assigned. Each parallel computing thread independently performs the recursive solution calculation of the discretized wavefield time series within its allocated sub-region, according to the explicit time integration format and lumped mass matrix. After each parallel computing thread completes the recursive solution for each time step, the shared boundary between its sub-region and adjacent sub-region is identified, and the current time step wave field variable of the grid cell node on the shared boundary is extracted. The wave field variables of the grid cell nodes on the shared boundary are exchanged between the parallel computing threads corresponding to the adjacent sub-regions through the parallel communication interface. Based on the boundary wave field variables obtained through exchange, the wave field states of several parallel computing threads at the shared boundary are updated for their own sub-regions, thereby achieving cross-sub-region wave field connection and global computation synchronization, and preparing for the recursive solution of the next time step.

[0014] Furthermore, the step of dynamically allocating the computation task of the unstructured finite element mesh set corresponding to each sub-region and the associated discretized wave field time series to a corresponding parallel computing thread includes: Obtain the computational load of the unstructured finite element mesh set corresponding to each sub-region. The computational load is related to the size of the mesh set and the recursive solution complexity of the discretized wave field time series. Based on the real-time processing capability and communication topology of each parallel computing thread, with the goal of load balancing and minimizing communication overhead, each sub-region is dynamically mapped to a corresponding parallel computing thread. The geometric data, physical parameters, and initial conditions for the recursive solution of the unstructured finite element mesh set corresponding to each sub-region are sent to the mapped parallel computing thread to complete the dynamic allocation of the computing task.

[0015] Accordingly, a second aspect of the present invention provides an electronic device, comprising: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described explicit finite element parallel computation method for seismic wave simulation.

[0016] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described explicit finite element parallel computing method for seismic wave simulation.

[0017] The above-described technical solutions of the embodiments of the present invention have the following beneficial technical effects: 1. By placing a finite fault source model, realistic 3D topographic data, and heterogeneous velocity structures within a unified explicit finite element framework for fully coupled simulation, an integrated modeling of the physical processes from source rupture to wave propagation and ultimately to surface response is achieved. This overcomes the limitations of traditional methods where the source, propagation medium, and topographic effects are isolated or oversimplified. It can more realistically reproduce the propagation path, energy distribution, and interference patterns of seismic waves in complex geological environments, thereby significantly improving the physical rationality and accuracy of spatial distribution predictions of strong ground motion, especially in areas with dramatic topographic relief and significant lateral changes in the medium. 2. The use of unstructured finite element mesh generation and local refinement techniques can closely match the complex geometric features of the real surface morphology and subsurface velocity interface. Compared with traditional numerical methods that rely on regular meshes, by directionally refining the mesh in steep terrain areas and areas with strong velocity gradients, the analytical capability of key features of seismic wave propagation is effectively improved without significantly increasing the overall computational scale. This ensures the numerical stability and spatial accuracy of the simulation results and provides a reliable technical means for characterizing fine phenomena such as local terrain magnification and basin edge effects. 3. By combining domain decomposition strategies with parallel computing architecture, a large-scale 3D computational domain is efficiently decomposed into multiple load-balanced sub-regions, and synchronous computation and boundary coordination are performed on each parallel node based on a unified explicit solver. This fully utilizes high-performance computing resources to complete ultra-large-scale 3D seismic wave simulation tasks containing hundreds of millions of grid cells within an acceptable timeframe, thus breaking through the computational efficiency bottleneck and achieving a leap from methodological innovation to engineering application. This provides feasible technical support for conducting high-resolution, broadband seismic hazard analysis and seismic damage scenario construction. Attached Figure Description

[0018] Figure 1 This is a flowchart of an explicit finite element parallel computing method for seismic wave simulation provided in an embodiment of the present invention; Figure 2 This is a diagram of the underground medium velocity structure provided in an embodiment of the present invention; Figure 3 This is a comparison chart of the PGV distribution simulated by the algorithm provided in this embodiment of the invention and the actual station observations; Figure 4 This is a comparison diagram of the simulated waveform and the actual observed waveform provided by the algorithm in this embodiment of the invention; Figure 5 This is a snapshot of the surface wavefield simulated by the algorithm provided in this embodiment of the invention, showing the propagation of seismic waves in a sedimentary basin. Figure 1 ; Figure 6 This is a snapshot of the surface wavefield simulated by the algorithm provided in this embodiment of the invention, showing the propagation of seismic waves in a sedimentary basin. Figure 2 . Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments and the accompanying drawings. It should be understood that these descriptions are merely exemplary and not intended to limit the scope of the invention. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0020] The calculation method of this invention will be specifically explained below using a simulated scenario of the 2025 Tingri earthquake as an example. Please refer to... Figure 1 The first aspect of this invention provides an explicit finite element parallel computation method for seismic wave simulation, comprising the following steps: Step S100: Construct a three-dimensional computational domain for seismic wave propagation simulation. The three-dimensional computational domain is a physical space covering the target epicenter area. Based on the real terrain data and three-dimensional velocity structure data of the target epicenter area, generate an unstructured finite element mesh that fits the geological structure.

[0021] First, the target epicenter area (such as Tingri County and its surrounding areas) was identified, with the epicenter as its core. Based on this, a sufficiently large three-dimensional computational domain was constructed, for example, a physical space 220 km long (east-west), 220 km wide (north-south), and extending approximately 55 km below the Moho discontinuity. Then, high-resolution SRTM30 topographic data of this region was integrated with a three-dimensional crustal velocity model (such as FWEA23) constructed based on deep exploration data. The topographic relief and major subsurface velocity interfaces (such as the sedimentary basement, crystalline basement, and Moho discontinuity) were used as geometric constraints to drive an unstructured mesh generation algorithm. This generated a set of tetrahedral or hexahedral elements with surfaces conforming to the real terrain, internal elements adapting to velocity structure changes, and locally refined in areas of abrupt topographic changes and velocity gradient zones. This formed the final unstructured finite element mesh used for computation.

[0022] Step S200: Perform domain decomposition on the three-dimensional computational domain, dividing the unstructured finite element mesh into multiple sub-regions.

[0023] After constructing a high-precision unstructured finite element mesh, this massive mesh containing hundreds of millions of elements needs to be decomposed into domains to enable parallel computing. First, based on the mesh's element distribution and topological connectivity, a graph partitioning algorithm (such as METIS or ParMETIS) is used to automatically divide the entire 3D computational domain into multiple spatially continuous sub-regions with tightly connected internal elements. The partitioning process aims to minimize the shared boundary area between sub-regions while ensuring that each sub-region contains approximately the same number of elements, thereby achieving a balanced computational load and minimizing communication overhead. For example, in a heliocentric earthquake simulation, based on the number of nodes in the high-performance computing cluster (e.g., 1024 computing cores), a 220 km × 220 km × 55 km 3D computational domain is divided into 1024 sub-regions. The computational load of each sub-region remains within a preset balanced range, and communication overhead is reduced by decreasing the number of shared boundaries between sub-regions to minimize the exchange of complex physical quantities at the boundaries. The partitioning results generate element lists, node mapping relationships, and adjacency information for each sub-region, providing a data foundation for subsequent task allocation.

[0024] Step S300: The computational task of each sub-region is assigned to multiple parallel computing threads. Each parallel computing thread uses the unstructured finite element mesh, an explicit time integration scheme and a lumped mass matrix to discretize and solve the three-dimensional elastic wave equation containing damping terms, and synchronously calculates the discretized wave field time series of each sub-region.

[0025] Based on the sub-region division results generated in step S200, the mesh element set, material parameters, initial conditions, and boundary information corresponding to each sub-region are dynamically allocated to an independent parallel computing thread (such as a CPU core or GPU computing card). Each parallel computing unit loads the same explicit finite element solver kernel, which uses the central difference method as the explicit time integration format and diagonalizes the global mass matrix using a lumped mass matrix, thereby simplifying the acceleration solution to node-by-node scalar operations and significantly improving computational efficiency. During the solution process, each node independently performs time-step iterations for its responsible sub-region: within each time step, it first calculates the nodal internal forces (considering damping terms) based on the current displacement and velocity, and then, combined with the source equivalent nodal forces, explicitly calculates the nodal acceleration using the lumped mass matrix; then, it updates the nodal velocity and displacement. After completing the calculation for each time step, the node exchanges wavefield variables on the shared boundary with nodes in adjacent sub-regions through parallel communication interfaces such as MPI, ensuring the continuity and physical consistency of the wavefield at the sub-region interface. Taking heliocentric earthquake simulation as an example, each node uses a time step of approximately 0.005 seconds (determined by the CFL condition) and progresses synchronously over thousands of time steps. This ultimately yields the displacement and velocity time histories of all grid nodes within each sub-region over the entire duration of the seismic motion (e.g., 60 seconds), resulting in a discretized wavefield time series. This parallel solution strategy enables single-machine computation tasks that would normally take months to complete to be finished in hours, providing technical support for large-scale, high-resolution seismic wave simulation.

[0026] Step S400: Integrate the synchronous calculation results of all parallel computing threads to obtain the seismic wavefield distribution in the three-dimensional computing domain, and output the ground motion parameters of at least one target sub-region in the three-dimensional computing domain.

[0027] After parallel computation is completed, the master node aggregates the sub-regional wavefield data generated by all computation nodes and reconstructs the seismic wavefield distribution of the complete three-dimensional computational domain. For specific target sub-regions of interest, such as the locations of townships like Changsuo and Quluo in Dingri County, or engineering sites distributed along important transportation lines, key ground motion parameters such as peak ground acceleration (PGA), peak ground velocity (PGV), and response spectrum are extracted from their surface grid nodes. These high spatial resolution simulation results can be directly used to draw seismic intensity distribution maps, assess the ground motion amplification effect at the site, and can also be used as input conditions to further predict seismic damage distribution at the township or city scale by combining them with building structural vulnerability curves.

[0028] Through the complete process described above, from data preparation, source modeling, parallel solution to result extraction, this method achieves high-fidelity numerical simulation of seismic wave propagation processes involving complex source characteristics, realistic topography, and non-uniform crustal structures within a unified framework. Its technical advantages lie in significantly improving the physical basis and spatial resolution for predicting strong ground motion in complex tectonic regions lacking sufficient strong earthquake observation records (such as the southern edge of the Tibetan Plateau). It overcomes the limitations of traditional empirical methods in characterizing topography and basin effects, and provides a quantifiable and verifiable high-efficiency technical means for seismic hazard analysis, seismic fortification of major engineering projects, and regional earthquake damage emergency assessment.

[0029] Further, in step S100, based on the real terrain data and three-dimensional velocity structure data of the target epicenter area, an unstructured finite element mesh that fits the geological structure is generated, including: Step S110: Obtain high-resolution terrain data and a three-dimensional velocity model covering the target epicenter area.

[0030] In the seismic wave simulation of this region, the first step is to acquire high-precision surface and subsurface data covering a target area with a radius of approximately 100 kilometers centered on the epicenter. Specifically, surface topographic data can be obtained using an SRTM digital elevation model with a spatial resolution of 30 meters, or even higher-precision airborne LiDAR data, to accurately reflect the dramatic topographic relief, deep valleys, and ridge morphology of the Himalayan fringe. Simultaneously, existing three-dimensional crustal velocity models for the region need to be constructed or referenced, such as wave velocity structure models constructed by integrating geophysical exploration results from deep seismic sounding, receiver functions, and background noise imaging. This model should clearly reflect the geometry and spatial variations of transverse and longitudinal wave velocities at major velocity interfaces, including the sedimentary cover, upper crustal crystalline rocks, middle and lower crust, and the Moho discontinuity. These high-resolution topographic and velocity data collectively constitute the most fundamental inputs for subsequent physical modeling and numerical discretization.

[0031] Step S120: Based on high-resolution terrain data and a three-dimensional velocity model, construct a set of spatially discretized grid cells that conform to the geometric shape of the surface terrain undulations and the underground velocity interface to form an unstructured finite element mesh.

[0032] Using the topographic surfaces and key isovelocity surfaces in the 3D velocity structure (such as the Vs=3.0 km / s interface and the Moho surface) obtained above as geometric constraints, an unstructured mesh generation algorithm is driven. This process first triangulates the topographic surface and each velocity interface into a surface mesh to ensure geometric consistency of the interface mesh. Then, based on these surface meshes as boundaries, algorithms such as constrained Delaunay tetrahedralization or the front-stepping method are used to generate a set of volumetric mesh elements that fill the entire 3D computational domain. During mesh generation, a local mesh refinement strategy is implemented, particularly for ridges and valleys with drastic surface slope changes, and areas with significant lateral variations in subsurface velocity or the presence of low-velocity sedimentary basins. This strategy gradually thins the mesh outwards from the epicenter. For example, in sedimentary areas such as the Dingri Basin, the near-surface mesh size is significantly reduced to accurately analyze the amplification and filtering effects of low-velocity layers on seismic waves. The final unstructured finite element mesh has nodes and edges of its elements (mainly hexahedrons) that strictly conform to the real terrain surface and underground velocity interface. The element size varies smoothly in space according to the geological complexity and the accuracy requirements of the wave field simulation, thus forming a high-quality spatial discretization model that can faithfully represent the geometric features of complex geological structures and meet the stability conditions of subsequent explicit finite element calculations.

[0033] By integrating high-precision observation and detection data and utilizing unstructured grid generation technology, high-fidelity spatial discretization of complex geological structures in the real world (including dramatic topographic relief, sedimentary basins, velocity discontinuities, etc.) was achieved. This fundamentally overcomes the shortcomings of traditional regular grids or simplified layered models in accurately describing actual geological morphology, ensuring the high realism of the medium model on which seismic wave propagation simulations depend. This provides a crucial prerequisite for obtaining reliable and spatially detailed simulation results of strong ground motion.

[0034] Furthermore, in step S120, the construction of a spatially discretized grid unit set that conforms to the surface topography and the geometry of the subsurface velocity interface forms an unstructured finite element mesh, including: Step S121: Integrate high-resolution terrain data with a three-dimensional velocity model to establish geometric constraints that simultaneously constrain the interface between the surface morphology and the internal medium.

[0035] To achieve accurate geometric model construction, the acquired SRTM30 topographic data needs to be spatially registered and fused with a regional 3D crustal velocity model (such as FWEA23). First, the topographic elevation data is converted into a continuous digital elevation surface (DEP). Simultaneously, key velocity isosurfaces or geological interfaces are extracted from the velocity model, such as the Vs=2.5 km / s interface representing the sedimentary basement, the upper and middle crustal boundary, and the Moho discontinuity. These surfaces, defined by different data sources, collectively constitute a multi-layered constraint framework in 3D space. Geometric processing ensures the consistency of these surfaces in spatial location (e.g., the surface topography is the highest interface), and a unified set of geometric constraint files or data structures that can be read by subsequent mesh generation algorithms is established, thus laying the geometric foundation for generating a computational mesh that closely conforms to the actual geological structure.

[0036] Step S122: Based on geometric constraints, an initial mesh composed of polyhedral elements is generated in the three-dimensional computational domain using an unstructured mesh generation algorithm. The surfaces and edges of the polyhedral elements conform to the geometry of the surface topography and the underground velocity interface.

[0037] Using the geometric constraints established above, which include surface topographic surfaces and multiple subsurface velocity interfaces, as input, an unstructured mesh generator built into software such as TetGen, Gmsh, or specialized geological modeling software is invoked. Within a predefined 3D computational domain bounding box, using these constraint surfaces as impenetrable boundaries, an initial volume mesh primarily composed of hexahedral elements is automatically generated. During this process, it is mandatory that the vertices, edges, and faces of each generated volume mesh element (polyhedron) strictly lie on or conform to these constraint surfaces, achieving what is known as "conformal conformation." Surface mesh nodes are precisely distributed along the topographic curves, and mesh nodes on both sides of the velocity interfaces are also precisely located on these interfaces, thus completely preserving the geometry of the real geological structure in the discretized mesh model and ensuring that physical properties (such as wave velocity) can change abruptly in the correct spatial location.

[0038] Step S123: Perform local refinement on the initial mesh, gradually sparsening outwards from the epicenter, to obtain an unstructured finite element mesh that fits the geological structure.

[0039] The initial mesh generated in step S122 is analyzed to identify key areas requiring higher resolution simulation. This method employs a local mesh refinement strategy that gradually thins outwards from the epicenter, ensuring the highest mesh resolution in the messeismic zone near the epicenter. This allows for the resolution of higher frequency seismic wave components and yields wider-bandwidth calculation results. Specifically, slope and curvature are calculated based on high-resolution topographic data, and steep, undulating ridgelines and canyons are marked as refinement zones. Simultaneously, the three-dimensional velocity model is analyzed, and areas with large lateral variations in shear wave velocity (such as the edge of the Dingri Basin and near major fault zones) are also marked as refinement zones. Subsequently, a mesh optimization algorithm is initiated, gradually transitioning from the epicenter to the surrounding areas. Within these marked areas, the distribution density of mesh nodes is significantly increased by iteratively subdividing existing polyhedral elements (e.g., splitting a large tetrahedron into multiple smaller tetrahedrons). The refinement process follows certain transition rules to ensure smooth mesh size changes between refined and unrefined areas, avoiding overly sharp jumps in element size to maintain the stability of numerical calculations. After optimization, the final unstructured finite element mesh has the highest resolution in the mesothermic zone near the epicenter and gradually becomes sparser towards the surrounding areas. This allows for a more detailed characterization of the topographic focusing effect and the propagation details of waves in complex velocity structures, while also taking into account computational efficiency and broadband simulation requirements.

[0040] Precise geometric constraints were established through data fusion, and an initial discrete model faithful to the geological structure was generated using conformal fitting technology. Furthermore, optimization based on terrain and physical properties intelligently concentrated computational resources on key areas affecting the accuracy of the wavefield simulation. This entire process ensures that the generated unstructured finite element mesh not only geometrically accurately reproduces the complex morphology of the surface and subsurface, but also achieves a balance between efficiency and accuracy in computational resource allocation, providing a crucial high-quality spatial discretization foundation for subsequent high-precision, high-efficiency seismic wave numerical simulations.

[0041] Further, in step S300, the three-dimensional elastic wave equation containing damping terms is discretized and solved using an explicit time integration scheme and a lumped mass matrix to obtain a discretized wavefield time series in the three-dimensional computational domain, including: Step S310: Based on the unstructured finite element mesh, the three-dimensional elastic wave equation containing damping terms is spatially discretized, and the three-dimensional elastic wave equation is transformed into a dynamic system defined on the mesh element nodes.

[0042] Using the aforementioned unstructured finite element mesh, the three-dimensional elastic wave equation describing seismic wave propagation, including damping terms, is spatially discretized. First, within each tetrahedral or hexahedral element, appropriate linear or quadratic functions are selected as interpolation basis functions for the displacement field. The partial differential equations of the continuous medium, including inertial forces, damping forces, and elastic restoring forces, are transformed into weak expressions in integral form that hold true on each element, involving the calculation of the element stiffness matrix, element damping matrix, and element mass matrix. The damping term is introduced by adding an integral term related to the velocity field in the weak form, typically using the Rayleigh damping assumption, expressing the damping matrix as a linear combination of the mass matrix and stiffness matrix to characterize the dissipation characteristics of seismic wave energy in the medium. Subsequently, based on the mesh topology, the stiffness, damping, and mass matrices contributed by all elements are assembled into global stiffness, global damping, and global mass matrices according to node numbers, and the source dynamic parameters are introduced into the system in the form of equivalent nodal forces. Ultimately, the continuous wave equation with damping terms is transformed into a large-scale second-order ordinary differential equation system with the displacements (and velocities) of all grid element nodes as unknowns—a discretized dynamic system. Its standard form is [M]{a} + [C]{v} + [K]{u} = {F}, where [C] is the global damping matrix and {v} is the nodal velocity vector. This system has dimensions ranging from tens of millions to hundreds of millions, precisely describing the coupling relationships and energy dissipation mechanisms of motion on discrete nodes.

[0043] Step S320: The dynamic system is solved recursively in the time domain using an explicit time integration scheme that alternately updates velocity and displacement variables.

[0044] For the aforementioned large-scale dynamic system, an explicit central difference method is employed for recursive solution in the time domain. The core of this method is to alternately and explicitly update the velocity and displacement fields of the nodes within each time step. Specifically, given the node displacements and velocities at the current time step, the node accelerations for the next time step are directly calculated using the equations of motion. Subsequently, using this acceleration and the current velocity, the velocity for the next time step is updated through time step integration. Finally, using the updated velocity and the current displacement, the displacement for the next time step is explicitly updated. The entire process does not require solving a global system of simultaneous equations, only simple vector operations. The time step is rigorously determined based on the minimum grid cell size and the highest wave velocity of the medium, using numerical stability conditions (such as the Courant-Friedrich-Lévy condition) to ensure numerical stability in seismic time history simulations lasting tens of seconds in complex, non-homogeneous media.

[0045] Step S330: In the time-recursive solution process, the lumped mass matrix is ​​used to replace the coordinated mass matrix, the mass matrix of the dynamic system is diagonalized, and the solution of the dynamic system is simplified to a node-by-node scalar operation in an explicit format.

[0046] By employing a concentrated (or clustered) mass matrix technique, the mass of each element is distributed and concentrated across its nodes according to specific rules (such as row summation), thus transforming the global mass matrix into a diagonal matrix. This crucial operation diagonalizes the mass terms multiplied by nodal accelerations in the dynamic system's equations of motion. Therefore, in each time step of step S320, the acceleration calculation for each node no longer depends on solving the system of equations, but simplifies to obtaining it independently by dividing the resultant force on that node (from the stiffness contribution and source force of surrounding elements) by its corresponding concentrated mass (a scalar). This transforms the solution of large-scale systems into fully parallel, node-by-node scalar operations, significantly reducing the computational complexity of a single step and laying the algorithmic foundation for subsequent efficient parallel computation.

[0047] Step S340: Through spatial discretization, time recursion and diagonalization of the mass matrix, a discretized wave field time series in the three-dimensional computational domain is generated.

[0048] By continuously executing the explicit recursive and simplified solution process described in steps S320 and S330, starting from the initial static state, displacement, velocity, and acceleration values ​​at all grid cell nodes in the entire three-dimensional computational domain can be obtained at each discrete time step. These fields spatially cover tens of kilometers underground to the surface, and temporally cover the entire process from rupture occurrence to seismic subsidence. Arranging and storing the node motion states of all time steps in chronological order constitutes a complete discretized wavefield time series. This series data is the original high-fidelity simulation result, from which the seismic motion time history at any location on the surface can be further extracted, or the spatial distribution field reflecting the seismic motion intensity, such as peak ground acceleration field, peak velocity field, and spectral acceleration fields of different periods, can be obtained through post-processing calculations, providing direct input for disaster analysis.

[0049] Through the systematic implementation of steps S310 to S340, the complex physical problem of seismic wave propagation is transformed into a large-scale numerical computation task that can be efficiently executed on a high-performance computer. This method ensures the geometric and physical property accuracy of the model through finite element spatial discretization on an unstructured mesh, avoids time-consuming implicit iterative solutions through explicit time integration, and minimizes computational complexity through lumped mass matrix techniques. The combination of these three aspects makes it possible to perform full-waveform simulations of large-scale three-dimensional models containing hundreds of millions of degrees of freedom, heterogeneous media, and complex terrain. This enables high-resolution revelation of the complex spatial distribution characteristics of near-fault strong ground motion, providing a powerful numerical experimental tool for understanding key scientific issues such as topographic effects and basin response, as well as for setting ground motion inputs in seismic engineering.

[0050] Further, in step S310, the three-dimensional elastic wave equation is spatially discretized based on an unstructured finite element mesh, transforming the three-dimensional elastic wave equation into a dynamic system defined on the mesh element nodes, including: Step S311: Discretize the three-dimensional computational domain into a set of interconnected polyhedral elements using an unstructured finite element mesh, where the vertices of each polyhedral element constitute mesh element nodes.

[0051] Using a pre-generated, geologically conforming unstructured finite element mesh, the entire 3D computational domain has been divided into a seamless collection of massive polyhedral elements (primarily tetrahedral elements in this model). Each tetrahedral element is uniquely identified by its four vertices, which are called "mesh element nodes." These nodes are precisely situated on the real terrain surface or subsurface velocity interfaces, ensuring that the volume and shape of each element formed by these nodes accurately reflect the geometric characteristics of its local geological space. For example, in the Dingri Basin, the elements are small and regularly shaped; in steep mountainous areas, the elements may exhibit a larger aspect ratio to adapt to rapid terrain changes. This interconnected element-node network provides a complete geometric framework for mapping continuous physical space into a discrete, numerically computable data structure.

[0052] Step S312: Based on the geometric and topological connection relationship of the unstructured finite element mesh, select the shape function defined on the polyhedral element to construct a weak form expression of the three-dimensional elastic wave equation at the element level.

[0053] Within each polyhedral element (such as a tetrahedron) established in step S311, an approximate expression of the physical field (displacement) needs to be established. To this end, a set of polynomial functions defined on each element is selected as "shape functions." Typically, for linear hexahedral elements, the shape functions are linear functions of the coordinates, ensuring that the displacement field within the element is continuous and linearly changing. Based on the minimum potential energy principle of elastic dynamics or the Galerkin weighted residual method, the physical laws (constitutive relations) and equilibrium equations describing the stress-strain relationship are integrally weighted and averaged within the element domain. This process transforms the three-dimensional elastic wave equation, originally acting at every point in the continuous medium in the form of partial differential equations, into a weak integral expression that holds true on each discrete element, with nodal displacements as unknowns. This expression naturally contains the geometric information of the element (through the Jacobian matrix) and the material physical properties (such as density and Lamé constant), ultimately forming the mathematical basis for calculating the element stiffness matrix and mass matrix.

[0054] Step S313: The three-dimensional elastic wave equation is expressed in a weak form at the element level. Based on the node connection relationship of the unstructured finite element mesh, it is assembled into a global dynamic equation system with the displacement and velocity at the mesh element nodes as unknown variables, thus completing the transformation of the three-dimensional elastic wave equation into a dynamic system.

[0055] Based on the topological connectivity of the unstructured mesh, thousands of weakly formal expressions at the element level are "assembled" into a unified global system. Specifically, this involves traversing all elements and accumulating each coefficient in their calculated element stiffness and mass matrices, according to their corresponding local node numbers, into the corresponding rows and columns of the global stiffness and mass matrices (determined by the global node number). Simultaneously, the spatiotemporal slip distribution provided by the finite fault source model, equivalent to time-history force vectors acting on relevant fault nodes, is also assembled into the global load vector. Finally, a large-scale, sparse global dynamic equation set of the form [M]{a} + [K]{u} = {F} is generated, where [M] is the global mass matrix, [K] is the global stiffness matrix, {u} and {a} are the displacement and acceleration vectors of all mesh element nodes, respectively, and {F} is the load vector. Thus, the seismic wave propagation problem controlled by continuous partial differential equations is completely transformed into an algebraic dynamic system about discrete node motion, suitable for computer solution.

[0056] Through the systematic execution of steps S311 to S313, a rigorous and faithful mathematical transformation from continuous geological media to a discrete algebraic system was achieved. First, a discrete spatial architecture consistent with the geometry of complex geological structures was established using an unstructured mesh. Then, shape functions and variational principles from the finite element method were used to construct local mathematical models conforming to physical laws on each element. Finally, these models were integrated into a global dynamic system through a mesh topology-based assembly process. This process mathematically rigorously ensures that the discrete model approximates the original continuous physical problem, laying a precise mathematical foundation for subsequent large-scale numerical solutions using efficient and stable explicit time integration algorithms. It serves as a crucial bridge connecting the complexity of real-world geology with high-performance numerical simulation.

[0057] Furthermore, step S320 employs an explicit time-integral scheme that alternately updates velocity and displacement variables to perform a recursive solution of the dynamic system in the time domain, including: Step S321: Discretize the simulation time into a series of equally spaced time steps.

[0058] Based on the numerical stability requirements of the explicit time integration algorithm, the total duration of the entire ground motion simulation is first determined, for example, simulating the wavefield evolution within 50 seconds after a Dinghe earthquake. Then, according to the Courant-Friedrich-Lévy (CFL) stability condition, a fixed time step Δt that satisfies the global stability requirement is calculated and selected. This step is jointly determined by the minimum size of all elements in the entire three-dimensional computational domain (typically in the indented region) and the highest wave velocity (typically in the high-velocity bedrock region). In the Dinghe region model, considering the complex near-surface low-velocity layer and detailed topography, the calculated stability step is typically on the order of a few ten-thousandths of a second. Finally, the entire simulation duration is discretized into tens of thousands of continuous and equally spaced tiny computational steps, providing a time frame for subsequent stepwise recursion.

[0059] Step S322: Within each time step, the acceleration of the grid cell nodes at the current time step is calculated using an explicit time integration scheme based on the displacement field and velocity field of the dynamic system at the previous time step.

[0060] At the beginning of each time step (e.g., step n), the system knows the displacement and velocity fields of all grid element nodes at the end of the previous time step (step n-1). First, the nodal accelerations for the current time step need to be calculated. The core of this calculation is solving for the inertial terms in the dynamic system equations. Since a lumped mass matrix has been applied, and the global mass matrix is ​​diagonal, the acceleration of each node can be solved independently and explicitly. Specifically, for each node, the nodal force acting on it is calculated. This nodal force consists of two parts: first, the elastic restoring force generated by the deformation of surrounding elements (obtained by multiplying the global stiffness matrix by the nodal displacement vector); and second, the equivalent force allocated to the node at that moment by the finite fault source model. Then, dividing the nodal force by its corresponding lumped mass (i.e., the mass scalar allocated to the node) directly yields the acceleration component of the node at the current time step. This process is performed in parallel for all nodes throughout the computational domain.

[0061] Step S323: Based on the grid cell node acceleration at the current time step, update the grid cell node velocity field at the current time step using an explicit time integration format.

[0062] After obtaining the accelerations of all nodes at the current time step (step n), the velocity field is updated by time integration of the velocities using an explicit central difference scheme. The specific update formula is: the velocity vector at the current time step equals the velocity vector at the previous time step, plus the product of the acceleration vector at the current time step and the time step size Δt. This operation is completely vectorized; the velocity update of each node depends only on its own velocity at the previous moment and the currently calculated acceleration, without needing to solve simultaneously with other nodes. For example, in simulating the propagation of a heliocentric earthquake rupture, nodes located near the fault will receive significant acceleration at the moment the rupture arrives, thus their velocities are rapidly updated through this step, reflecting the arrival of seismic wave energy.

[0063] Step S324: Based on the grid cell node velocity field of the current time step, update the grid cell node displacement field of the current time step using an explicit time integration format.

[0064] After obtaining the updated velocity field at the current time step (step n), the displacement is again integrated over time using the explicit central difference scheme. The update formula is: the displacement vector at the current time step equals the displacement vector at the previous time step, plus the product of the current (updated) velocity vector and the time step size Δt. Thus, the system completes the transition from one state (displacement and velocity at step n-1) to the next new state (displacement and velocity at step n). This "displacement-velocity-acceleration" update sequence is a characteristic of the explicit central difference method, alternately using the most recently calculated field quantity to update another field quantity, ensuring the explicitness and recursiveness of the calculation.

[0065] Step S325: Repeat the above steps of calculating acceleration, updating velocity field and updating displacement field according to the time step order to complete the recursive solution for the entire simulation period and form a discretized wave field time series in which the displacement and velocity of the grid cell nodes change with time.

[0066] Starting with the set initial conditions (usually zero displacement and zero velocity, or an initial rupture disturbance applied at the fault), the loop consisting of steps S322 to S324 is repeated. Each time the loop is completed, physical time is advanced by Δt. In heliocentric earthquake simulations, this process needs to be repeated hundreds of thousands or even millions of times to cover the entire process from rupture triggering, wave propagation, surface shaking, to gradual subsidence. At each step of the loop, the system generates and stores (or caches) the displacement and velocity values ​​of all nodes at that moment. After the loop covers the entire preset simulation period, a complete dataset is obtained: the displacement and velocity history of each grid cell node at every discrete time point. This four-dimensional dataset (three-dimensional space + one-dimensional time) is the final "discrete wavefield time series," which completely records the evolution of seismic waves over time in the entire three-dimensional computational domain.

[0067] Through systematic recursion in steps S321 to S325, the solution to the complex dynamic system is transformed into a series of highly localized, parallelizable basic arithmetic operations. This method ensures numerical stability through strict time step control, simplifies acceleration calculations to scalar division using diagonalization of the lumped mass matrix, and achieves efficient time advancement without iteratively solving large-scale linear equation systems through an alternating explicit update strategy of "acceleration first, then velocity, then displacement." This entire mechanism makes it possible to perform long-term transient dynamic simulations of ultra-large-scale systems containing tens to hundreds of millions of degrees of freedom, thereby efficiently generating high-resolution, full-waveform strong ground motion time history data for earthquake engineering analysis.

[0068] Furthermore, in step S340, through spatial discretization, time recursion, and diagonalization of the mass matrix, the motion state sequence of several grid cell nodes in the three-dimensional computational domain over the entire simulation period is calculated, forming a discretized wavefield time series, including: Step S341: Collect the displacement and velocity fields of all grid cell nodes in the three-dimensional computational domain at each time step, obtained through spatial discretization, time recursion, and diagonalization of the mass matrix.

[0069] After parallel computing completes the recursive solution for all time steps, the sub-region wavefield data distributed across multiple computing nodes needs to be collected and integrated. The master node initiates a global collection operation through a high-speed interconnect network, sequentially reading the displacement and velocity vector data of all grid cell nodes within its responsible sub-region from each participating parallel node, in time step order. Considering that the diurnal simulation may involve hundreds of millions of nodes and hundreds of thousands of time steps, this is a data-intensive communication process. The collection process must ensure the order and integrity of data blocks and may be accompanied by lightweight real-time verification, such as checking whether the size of each time step data packet meets expectations, to detect silent data corruption that may be caused by memory errors or communication failures, laying the foundation for the subsequent construction of a reliable wavefield database.

[0070] Step S342: The displacement field and velocity field at each time step are structured according to the grid cell node number and time step order to obtain a set of motion state data with time as the dimension and spatial nodes as the index.

[0071] After obtaining the raw data, it needs to be converted into a standardized format that facilitates storage, retrieval, and post-processing. The massive dataset is organized according to two core dimensions: spatial and temporal. Spatially, an index is built based on the globally unique identifier of each grid cell node to ensure that data from nodes with the same name from different sub-regions are correctly merged. Temporally, data frames are arranged strictly according to the chronological order of the time steps. Finally, a hierarchical dataset is constructed, with a time series as its top-level structure. Each time point corresponds to a "snapshot" containing the displacement and velocity values ​​of all nodes in the entire computational domain. In heliocentric seismic applications, this dataset is typically written into a self-describing, compressible, and efficient scientific data format such as HDF5 or NetCDF. The dataset also stores the grid's geometric information (node ​​coordinates, cell connectivity) as metadata, enabling the complete reproduction and understanding of the wavefield sequence and its spatial carrier.

[0072] Step S343: Perform consistency verification and format standardization on the motion state data set, and output the discretized wave field time series of the three-dimensional computational domain throughout the entire simulation period.

[0073] After generating the structured dataset, final quality control and interface standardization are required. Consistency checks include, but are not limited to: checking whether the displacement and velocity fields spatially satisfy basic continuity (in non-fault regions), verifying that energy conservation is within acceptable numerical error ranges, and confirming that the time series does not exhibit non-physical oscillations or abrupt changes at critical moments (such as peak ground motion). Format standardization converts the data into specific standard formats according to the needs of the target application. For example, for use in earthquake engineering software, it may be necessary to extract the time histories of surface nodes and output them in SAC, SEED, or specific text formats; for visualization, it may be necessary to generate file sequences in VTK or XDMF formats. Ultimately, the output is a validated, uniformly formatted "discrete wavefield time series" product, which can be directly used for seismic motion parameter extraction, wavefield visualization, spectral analysis, and as input for other simulation programs (such as structural response analysis).

[0074] Through the systematic implementation of steps S341 to S343, the raw numerical results generated by large-scale parallel computing are transformed into reliable, structured, and standardized data products that can directly serve scientific analysis and engineering applications. This method ensures the integrity of simulation results through efficient data aggregation, achieves efficient management and rapid access to massive wavefield data through structured organization in the spatiotemporal dimensions, and guarantees the scientific reliability and broad compatibility of the data for downstream applications through rigorous consistency checks and format standardization. This complete data processing chain enables the accurate and convenient extraction and utilization of the rich physical information contained in high-fidelity numerical simulations, ultimately transforming the complex computational process into a quantitative knowledge base that can support earthquake hazard assessment, mechanism research, and seismic design in engineering.

[0075] Furthermore, in step S200, the three-dimensional computational domain is decomposed into multiple sub-regions, including: Step S210: Based on the element distribution and connection relationship of the unstructured finite element mesh, the three-dimensional computational domain is divided into several spatially continuous and computationally balanced sub-regions.

[0076] To handle large-scale computational tasks involving hundreds of millions of unstructured mesh elements, the first step is to effectively decompose the entire 3D computational domain into multiple parallelizable sub-regions. This process relies on specialized graph partitioning algorithms (such as METIS or Scotch). These algorithms take the topology of the unstructured finite element mesh as input, treating each mesh element as a vertex in the graph and shared surfaces between elements as connecting edges. The goal of partitioning is to generate several spatially continuous sub-regions with tightly connected internal elements, while rigorously balancing the total number of mesh elements in each sub-region (as a proxy for computational load) and minimizing the shared boundary area between sub-regions (to reduce subsequent parallel communication overhead). In the model for the Dinghe region, due to the dense mesh in near-surface basins and the relatively sparse mesh in mountainous areas, this algorithm can intelligently adjust the partitioning to ensure that each sub-region, despite its irregular geometry, has a roughly equivalent computational load, laying the foundation for efficient parallel computing.

[0077] Accordingly, step S300, which involves allocating the computational task of each sub-region to multiple parallel computation threads, includes: Step S301: Based on the sub-region division results, the calculation tasks of the unstructured finite element mesh set and the associated discretized wave field time series corresponding to each sub-region are dynamically allocated to a corresponding parallel computing thread.

[0078] After obtaining the sub-region partitioning scheme, task mapping and data distribution need to be performed on a high-performance computing platform composed of multiple computing nodes (such as CPU cluster nodes or GPU computing cards). A dynamic task allocator (usually a parallel job script or a dedicated runtime library) assigns each partitioned sub-region to the most suitable computing node based on the real-time available computing resources and memory capacity of each computing node, as well as the network topology of the cluster. Load balancing and communication proximity need to be considered during allocation to optimize overall performance. Subsequently, the system transmits all data corresponding to each sub-region, including the geometric information (node ​​coordinates, cell connections), physical properties (density, elastic parameters) of its internal unstructured mesh cell set, and the initial conditions and boundary information of the sub-region in the global dynamic system, through a network file system or high-speed interconnect network (such as InfiniBand) and loads it into the local memory of the corresponding parallel computing thread to complete the deployment of the computing task.

[0079] Step S302: The computation tasks allocated to several parallel computing threads are configured to independently perform the recursive solution of the discretized wave field time series within their assigned sub-regions based on the same explicit time integration format and lumped mass matrix, and to synchronize wave field variables at the shared boundary of adjacent sub-regions, thereby achieving collaborative parallel computation of the entire three-dimensional computational domain.

[0080] After loading their respective sub-region data, all computational nodes independently initialize the same numerical solver instance. This solver kernel is uniformly configured with the aforementioned explicit time integration format and lumped mass matrix algorithm. Within their assigned sub-region, each node independently performs the recursive solution calculation of the seismic wavefield. This process is algorithmically identical to single-machine serial computation, but only for local data. However, due to shared boundaries between sub-regions, the motion state of grid cell nodes on the boundaries depends on information from neighboring sub-regions. Therefore, after computation at each time step, all nodes must perform global synchronization: each node identifies the "outer boundary" of its own sub-region, extracts the current displacement and velocity values ​​of all nodes on the boundary, and then exchanges this boundary data with nodes holding neighboring sub-regions via point-to-point communication using a message passing interface (MPI). After the exchange, each node updates the state of its own boundary nodes with the received neighbor data to ensure the physical continuity of the wavefield propagation throughout the computational domain. This "computation-communication-update" cycle repeats at each time step, thereby enabling collaborative parallel computation of all sub-regions, which is logically equivalent to performing the simulation of the entire 3D model in a virtual global memory space.

[0081] Through the systematic implementation of steps S301 to S302, a large-scale 3D seismic wave simulation task, originally limited by single-machine computing resources, is decomposed into a set of concurrently executable, moderately sized sub-tasks, and their results are integrated through a sophisticated coordination mechanism. This method achieves balanced distribution of computational load through intelligent domain decomposition based on grid topology, optimizes computational resource utilization through dynamic task allocation, and ensures the physical correctness of distributed computing through a rigorously designed boundary data synchronization protocol. This parallel computing framework effectively overcomes the bottlenecks of single processors in terms of memory capacity and computing speed, making high-resolution, long-scale full-waveform dynamic simulation of models containing extremely complex geological structures and fine meshes a realistic and feasible technical approach, greatly expanding the capabilities of numerical simulation methods to solve practical engineering and scientific problems.

[0082] Furthermore, in step S302, the computational tasks allocated to several parallel computing threads are configured to independently perform recursive solutions of the discretized wavefield time series within their assigned sub-regions based on the same explicit time integration scheme and lumped mass matrix, and to synchronize wavefield variables at the shared boundaries of adjacent sub-regions, including: In step S302a, each parallel computing thread is configured with the same numerical solver kernel containing an explicit time integration format and a lumped mass matrix, and loads the unstructured finite element mesh set and initial field information corresponding to the sub-region to which it is assigned.

[0083] After domain decomposition and task allocation, each participating parallel node (e.g., a computing blade or a GPU accelerator card in the cluster) initializes an identical numerical solver software instance. This solver kernel is pre-compiled and pre-defined with the specific explicit time integration format (e.g., central difference method) and lumped mass matrix processing logic used in this method, ensuring strict consistency across all nodes at the algorithm level. Subsequently, each node reads its assigned sub-region data package via a parallel file system or network. This package contains the geometric description (node ​​coordinates, element connectivity), material physical parameters (density, Lamé constant), and initial displacement and velocity fields of the sub-region nodes obtained from global initial conditions mapping, covering all unstructured finite element mesh elements within that sub-region. Once loaded, each node logically possesses a complete, independently computable copy of the local physical model.

[0084] In step S302b, each parallel computing thread independently performs recursive solution calculations of the discretized wavefield time series within its assigned sub-region, according to the explicit time integration format and the lumped mass matrix.

[0085] Each node, within its own memory space, independently initiates an explicit time-recursive solution for the seismic wavefield based on its loaded local mesh model. This process is identical to the serial single-machine simulation algorithm: within each time step, the node calculates the internal forces of its internal nodes based on the topology of its local mesh, superimposes the source forces that may act on that sub-region, and explicitly solves for the node acceleration using the advantage of diagonalization of the lumped mass matrix, thereby updating the node velocity and displacement. In the Dingri earthquake simulation, a node responsible for a sub-region containing a fine mesh of the Dingri Basin sedimentary layers and another node responsible for a relatively sparse mesh of the surrounding bedrock mountain area, despite differences in the number of elements and material properties they process, both follow the exact same mathematical formulas and code paths for calculation, without interfering with each other, thus achieving perfect spatial partitioning and parallel execution of computational tasks.

[0086] Step S302c: After each parallel computing thread completes the recursive solution for each time step, the shared boundary between its sub-region and adjacent sub-region is identified, and the current time step wave field variable of the grid cell node on the shared boundary is extracted.

[0087] Once a node has completed all internal calculations for its subregion within the current time step, it must handle interactions with neighboring subregions before moving to the next time step. To do this, each node must identify all "shared boundaries" of its subregion based on the "subregion adjacency table" generated during global domain decomposition and the "boundary node markers" of the local mesh. These boundaries consist of a series of mesh element nodes whose motion depends on the calculation results of both the current and neighboring subregions. The node program then retrieves the displacement and velocity variables (i.e., wavefield variables) of these boundary nodes, which were just calculated at the current time step, from memory and packages them into a data buffer to be sent.

[0088] In step S302d, the wave field variables of the grid cell nodes on the shared boundary are exchanged between the parallel computing threads corresponding to adjacent sub-regions through the parallel communication interface.

[0089] After the data is packaged, each node initiates non-blocking point-to-point communication operations through the underlying parallel communication interface (such as MPI). Each node asynchronously sends data packets containing the wavefield variables of the relevant boundary nodes to all neighboring nodes that share boundaries with its sub-regions; simultaneously, it also prepares to receive corresponding data packets from these neighboring nodes. In high-performance computing environments for heliocentric seismic simulations, this communication is typically conducted through a high-speed InfiniBand network to minimize latency. The communication library manages these send and receive requests and ensures that all necessary boundary data exchanges are completed before continuing execution (through synchronization mechanisms such as MPI_Waitall). This process enables timely interoperability of physically adjacent local wavefield information held by different computing nodes across the entire network.

[0090] Step S302e: Based on the boundary wave field variables obtained by the exchange, update the wave field state of the sub-regions of several parallel computing threads at the shared boundary, realize the wave field connection across sub-regions and global computing synchronization, and prepare for the recursive solution of the next time step.

[0091] After successfully receiving boundary data packets from all neighboring nodes, each node immediately updates its boundary state. The received displacement and velocity data from the neighboring sub-regions' boundary nodes are fused with the corresponding boundary node values ​​calculated by the node itself. Typically, for nodes located on shared geometric interfaces, their motion state should be unique. Therefore, the values ​​received from neighbors can directly overwrite the local values, or conservative interpolation methods such as weighted averaging can be used to ensure the continuity and physical consistency of the wavefield at the sub-region interface. After the boundary update, the wavefield state of the sub-region at that moment becomes globally consistent and physically coordinated. Subsequently, the node program uses this updated state as a known condition, loops back to step S332, and begins the local independent calculation for the next time step. This process repeats until the simulation ends, thus achieving a synchronized, collaborative, parallel simulation of the wavefield evolution across the entire three-dimensional computational domain in a distributed memory environment.

[0092] By iteratively executing steps S302a to S302e, a highly efficient and robust large-scale distributed parallel computing framework is constructed. This framework ensures algorithm consistency by deploying a unified solver kernel across all computing nodes, maximizes the computational performance of individual nodes by allowing nodes to perform local computations completely independently, and addresses the core challenge of data consistency in distributed memory systems with minimal communication overhead through a carefully designed "computation-communication-update" protocol. This mechanism makes it possible to deploy ultra-large-scale, unstructured 3D seismic wave simulation problems onto modern high-performance computing clusters composed of tens of thousands of processor cores, thereby enabling high-resolution full-waveform simulations of extremely complex real-world geological scenes within an acceptable timeframe, significantly improving the scale and efficiency of numerical simulation methods in solving practical problems.

[0093] Furthermore, step S301, which involves dynamically allocating the computation task of the unstructured finite element mesh set and the associated discretized wave field time series corresponding to each sub-region to a corresponding parallel computing thread, includes: Step S301a: Obtain the computational load of the unstructured finite element mesh set corresponding to each sub-region. The computational load is related to the size of the mesh set and the recursive solution complexity of the discretized wave field time series.

[0094] To efficiently allocate tasks, the computational load of each sub-region needs to be quantified first. In this simulation scenario, the computational load mainly depends on two core factors: first, the number and type of unstructured finite element mesh elements contained in the sub-region (e.g., tetrahedral elements have different computational costs than hexahedral elements), which directly determines the amount of element stiffness matrix computation required in each step; second, the complexity of the material in the sub-region. For example, in the Dingri Basin, the low seismic wave propagation velocity in the low-velocity sedimentary layer leads to a smaller stabilization time step, and nonlinear or complex constitutive relations may need to be handled, all of which increase the number of computation steps required per unit physical time and the computational complexity of each step. The system analyzes the mesh statistics (number of elements, number of nodes, element type distribution) and material property distribution (such as velocity structure) of each sub-region, and combines the characteristics of the selected explicit time integration algorithm to estimate the relative computational load (measured by floating-point operations or computation time) required for each sub-region to complete the entire simulation period, forming a detailed load evaluation table to provide a basis for subsequent dynamic mapping.

[0095] Step S301b: Based on the real-time processing capability and communication topology of each parallel computing thread, with the goal of load balancing and minimizing communication overhead, dynamically map each sub-region to a corresponding parallel computing thread.

[0096] After understanding the load of each sub-region and the status of available parallel computing threads in the cluster, an optimization mapping process needs to be executed. This process considers the real-time processing capabilities of each computing node (such as the number of CPU cores, clock speed, memory bandwidth, or the number of stream processors and memory size of the GPU) and its load in the current job scheduling system, while also considering the cluster's communication topology (such as whether the network connection between nodes is a fat tree, mesh, or other topology, and the latency differences caused by physical distance). The mapping algorithm (such as a graph partitioning-based load balancing algorithm) has a dual goal: First, to balance the total computing load of the sub-region set allocated to each computing node as much as possible, to avoid the "barrel effect" where some nodes finish prematurely while others are still computing; second, to map sub-regions that are physically adjacent and therefore need to frequently exchange boundary data during the simulation to computing node pairs with lower network communication overhead (e.g., different processor cores within the same physical machine, or nodes directly connected via high-speed links), thereby minimizing cross-node communication latency and bandwidth consumption and improving overall parallel efficiency.

[0097] In step S301c, the geometric data, physical parameters, and initial conditions for recursive solution of the unstructured finite element mesh set corresponding to each sub-region are sent to the mapped parallel computing thread to complete the dynamic allocation of computing tasks.

[0098] Once the optimized mapping from sub-regions to computing nodes is complete, the system initiates large-scale data distribution. For each sub-region, its corresponding complete data package, including node coordinates and element connection tables describing the mesh geometry, density and elastic parameter fields defining material properties, initial force vectors and motion states on the nodes of that sub-region obtained from the global source model and initial conditions mapping, and the connection relationships (shared boundary information) between that sub-region and its neighboring sub-regions, is efficiently transmitted to the local storage or memory of the target computing node to which it is mapped, via a parallel file system (such as Lustre or GPFS) or point-to-point communication coordinated by the master node. During data transmission, data compression and verification are typically performed to ensure transmission efficiency and accuracy. Once the target node confirms receipt of all sub-region data it is responsible for and has completed loading, the node enters a ready state, awaiting a global synchronization signal to begin parallel computation. At this point, the initial dataset for the entire large-scale 3D seismic simulation has been distributed across the nodes of the entire high-performance computing cluster, completing the dynamic allocation and deployment of computational tasks.

[0099] Through the systematic execution of steps S301a to S301c, intelligent decomposition and resource optimization of ultra-large-scale computing tasks are achieved. This method provides a quantitative basis for load balancing by accurately assessing the computing load of sub-regions based on grids and physical models; it effectively controls the crucial communication overhead in parallel computing by comprehensively considering a dynamic mapping strategy that integrates computing node performance and network topology; and it ensures that massive amounts of initial data can be accurately and quickly deployed to the distributed computing environment through an efficient and reliable data distribution mechanism. This task allocation and data deployment process lays a solid foundation for the subsequent collaborative and conflict-free execution of high-resolution seismic wave simulations by thousands or even tens of thousands of computing cores, serving as a key link in realizing the transition from complex physical problems to efficient parallel computing practices.

[0100] This invention has been fully implemented and its effectiveness verified through the case of the 2025 Tingri earthquake. The specific implementation and effects of this method can be systematically described below with reference to the accompanying drawings: This invention begins with the fusion of multi-source data and high-precision modeling, such as Figure 2 As shown, the FWEA23 three-dimensional velocity model was used. Figure 2 The method clearly characterizes the significant non-uniformity of P-wave and S-wave velocities in the Earth's crust, providing realistic physical medium parameters for wave simulation. Simultaneously, the SRTM30 topographic data defines the complex free surface morphology at the top of the computational domain. Based on this, the method generates an unstructured mesh that closely matches the real terrain and subsurface interface. This mesh is divided into approximately 400 million cells within the computational domain, thus providing a discrete spatial foundation for simulating wave field systems with up to 1.2 billion degrees of freedom.

[0101] In terms of solution strategy, the core of this method is to use an explicit velocity-displacement alternating update scheme to perform time integration on the elastic wave equation. This algorithm choice ensures numerical stability in large-scale computations. To practically handle this ultra-large-scale computational problem, the technical solution implements domain decomposition and parallel computation. For example... Figure 3 The holistic demonstration, with its detailed simulation results covering a range of 220 km × 220 km, is made possible by decomposing the entire computational domain into multiple sub-regions and distributing them to parallel computing threads for collaborative completion. Each node, based on a unified explicit algorithm kernel, independently computes its local sub-domain and synchronously shares wavefield variables at the boundaries through communication, ultimately efficiently integrating them into a global result.

[0102] like Figure 4 As shown, the waveform comparison diagrams provide conclusive evidence. For example, at stations such as XZD6007, the simulated waveforms and observed records exhibit a high degree of consistency in the arrival time, waveform morphology, and amplitude of the main seismic phases. This directly confirms the reliable accuracy of the entire chain of this method, from source setting and wave propagation simulation to surface motion calculation. Furthermore, Figure 2 It not only macroscopically displays the spatial distribution pattern of peak ground velocity, but also directly marks the specific numerical comparison between simulated and observed values ​​at the station locations (such as "Sim: 0.40 m / s; Obs: 0.42 m / s"). This quantitative agreement strongly supports the authenticity of the simulation results.

[0103] Based on reliable calculation results, this method effectively reveals the spatial complexity and physical mechanism of strong ground motion. For example... Figure 2 The color distribution in the data indicates that the strong seismic activity was not uniformly distributed, but rather exhibited a distinct strip-shaped high-value zone along the ridgeline north of the fault. This directly reflects the coupled amplification effect of the rupture directionality and topographic focusing. Simultaneously, a significant high-PGV region also appeared in the basin area south of the fault. To further analyze this phenomenon, a parameterized basin model was constructed for simulation testing. Correlation analysis shows that, as... Figure 5 and Figure 6 The study revealed that sedimentary layers (especially low-velocity layers) can trigger significant basin edge effects, resonating and amplifying seismic waves in specific frequency bands, thereby exacerbating seismic damage in the region. These findings were quantified using the simulation capabilities of this method.

[0104] In summary, this method integrates explicit finite element algorithms, unstructured meshes, and parallel computing techniques, and closely combines multi-source measured data to construct a complete technical process from physical modeling and efficient solution to result verification. The velocity model, waveform comparison, ground motion distribution, and mechanism analysis diagrams shown in the attached figures collectively constitute strong empirical evidence of the method's ability to achieve high-precision seismic wave simulation and reveal the laws governing complex ground motion fields. This technical solution provides a key tool for quantitative seismic hazard analysis and seismic damage prediction in complex tectonic zones where observations are scarce.

[0105] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor; and a memory connected to the at least one processor. The memory stores instructions executable by the at least one processor, which, when executed by the at least one processor, cause the at least one processor to perform the aforementioned explicit finite element parallel computation method for seismic wave simulation.

[0106] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described explicit finite element parallel computing method for seismic wave simulation.

[0107] The embodiments of this invention aim to protect an explicit finite element parallel computation method for seismic wave simulation, which has the following effects: 1. By placing a finite fault source model, realistic 3D topographic data, and heterogeneous velocity structures within a unified explicit finite element framework for fully coupled simulation, an integrated modeling of the physical processes from source rupture to wave propagation and ultimately to surface response is achieved. This overcomes the limitations of traditional methods where the source, propagation medium, and topographic effects are isolated or oversimplified. It can more realistically reproduce the propagation path, energy distribution, and interference patterns of seismic waves in complex geological environments, thereby significantly improving the physical rationality and accuracy of spatial distribution predictions of strong ground motion, especially in areas with dramatic topographic relief and significant lateral changes in the medium. 2. The use of unstructured finite element mesh generation and local refinement techniques can closely match the complex geometric features of the real surface morphology and subsurface velocity interface. Compared with traditional numerical methods that rely on regular meshes, by directionally refining the mesh in steep terrain areas and areas with strong velocity gradients, the analytical capability of key features of seismic wave propagation is effectively improved without significantly increasing the overall computational scale. This ensures the numerical stability and spatial accuracy of the simulation results and provides a reliable technical means for characterizing fine phenomena such as local terrain magnification and basin edge effects. 3. By combining domain decomposition strategies with parallel computing architecture, a large-scale 3D computational domain is efficiently decomposed into multiple load-balanced sub-regions, and synchronous computation and boundary coordination are performed on each parallel node based on a unified explicit solver. This fully utilizes high-performance computing resources to complete ultra-large-scale 3D seismic wave simulation tasks containing hundreds of millions of grid cells within an acceptable timeframe, thus breaking through the computational efficiency bottleneck and achieving a leap from methodological innovation to engineering application. This provides feasible technical support for conducting high-resolution, broadband seismic hazard analysis and seismic damage scenario construction.

[0108] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. An explicit finite element parallel computation method for seismic wave simulation, characterized in that, Includes the following steps: A three-dimensional computational domain for simulating seismic wave propagation is constructed. The three-dimensional computational domain is a physical space covering the target epicenter area. Based on the real terrain data and three-dimensional velocity structure data of the target epicenter area, an unstructured finite element mesh that fits the geological structure is generated. The three-dimensional computational domain is decomposed into multiple sub-regions by dividing the unstructured finite element mesh. The computational task of each sub-region is assigned to multiple parallel computing threads. Each parallel computing thread uses the unstructured finite element mesh, an explicit time integration scheme and a lumped mass matrix to discretize and solve the three-dimensional elastic wave equation containing damping terms, and synchronously calculates the discretized wave field time series of each sub-region. By integrating the synchronous calculation results of all the parallel computing threads, the seismic wavefield distribution within the three-dimensional computational domain is obtained, and the ground motion parameters of at least one target region within the three-dimensional computational domain are output.

2. The explicit finite element parallel computation method for seismic wave simulation according to claim 1, characterized in that, The process of generating an unstructured finite element mesh that conforms to the geological structure based on the real terrain data and three-dimensional velocity structure data of the target epicenter area includes: Acquire high-resolution terrain data and a three-dimensional velocity model covering the target epicenter region; Based on the high-resolution terrain data and the three-dimensional velocity model, a set of spatially discretized grid cells that conforms to the surface topography and the geometry of the underground velocity interface is constructed to form the unstructured finite element mesh.

3. The explicit finite element parallel computation method for seismic wave simulation according to claim 2, characterized in that, The construction of a spatially discretized grid unit set that conforms to the geometry of the surface topography and underground velocity interface forms the unstructured finite element mesh, including: By integrating the high-resolution terrain data with the three-dimensional velocity model, geometric constraints are established that simultaneously constrain the surface morphology and the interface between the internal media. Based on the geometric constraints, an initial mesh composed of polyhedral elements is generated in the three-dimensional computational domain using an unstructured mesh generation algorithm. The surfaces and edges of the polyhedral elements conform to the geometry of the surface topography and the underground velocity interface. The initial mesh is locally refined, gradually becoming sparser outwards from the epicenter, to obtain an unstructured finite element mesh that conforms to the geological structure.

4. The explicit finite element parallel computation method for seismic wave simulation according to claim 1, characterized in that, The method employs an explicit time integration scheme and a lumped mass matrix to discretize and solve the three-dimensional elastic wave equation containing damping terms, obtaining a discretized wavefield time series within the three-dimensional computational domain, including: Based on the unstructured finite element mesh, the three-dimensional elastic wave equation containing damping terms is spatially discretized, and the three-dimensional elastic wave equation is transformed into a dynamic system defined on the mesh element nodes; An explicit time-integral scheme with alternating updates of velocity and displacement variables is used to recursively solve the dynamic system in the time domain. In the time-recursive solution process, the lumped mass matrix is ​​used to replace the coordination mass matrix, the mass matrix of the dynamic system is diagonalized, and the solution of the dynamic system is simplified to an explicit node-by-node scalar operation. The discretized wave field time series in the three-dimensional computational domain is generated through the spatial discretization, time recursion, and diagonalization of the mass matrix solution process.

5. The explicit finite element parallel computation method for seismic wave simulation according to claim 4, characterized in that, The step of spatially discretizing the three-dimensional elastic wave equation based on the unstructured finite element mesh, transforming the three-dimensional elastic wave equation into a dynamic system defined on the mesh element nodes, includes: The three-dimensional computational domain is discretized into a set of interconnected polyhedral elements using the unstructured finite element mesh, wherein the vertices of each polyhedral element constitute the mesh element node. Based on the geometric and topological connectivity of the unstructured finite element mesh, shape functions defined on the polyhedral elements are selected to construct a weak form expression of the three-dimensional elastic wave equation at the element level. The three-dimensional elastic wave equation is expressed in a weak form at the element level. Based on the node connection relationship of the unstructured finite element mesh, a global dynamic equation system is assembled with the displacement and velocity at the nodes of the mesh elements as unknown variables, thus completing the transformation of the three-dimensional elastic wave equation into the dynamic system.

6. The explicit finite element parallel computation method for seismic wave simulation according to claim 4, characterized in that, The explicit time integration scheme employing alternating updates of velocity and displacement variables to recursively solve the dynamic system in the time domain includes: Discretize the simulation time into a series of equally spaced time steps; Within each time step, the grid cell node acceleration of the current time step is calculated using the explicit time integration scheme based on the displacement field and velocity field of the dynamic system in the previous time step. Based on the grid cell node acceleration at the current time step, the grid cell node velocity field at the current time step is updated using the explicit time integration format; Based on the grid cell node velocity field at the current time step, the grid cell node displacement field at the current time step is updated using the explicit time integration format. Following the time step sequence, the steps of calculating acceleration, updating the velocity field, and updating the displacement field are repeated to complete the recursive solution for the entire simulation period, forming the discretized wave field time series in which the displacement and velocity of the grid cell nodes change with time.

7. The explicit finite element parallel computation method for seismic wave simulation according to claim 4, characterized in that, The process of spatial discretization, time recursion, and diagonalization of the mass matrix is ​​used to calculate the motion state sequence of several grid cell nodes in the three-dimensional computational domain throughout the entire simulation period, forming the discretized wavefield time series, including: The displacement and velocity fields of all grid cell nodes in the three-dimensional computational domain at each time step are collected through the spatial discretization, time recursion and mass matrix diagonalization solution process. The displacement field and velocity field at each time step are structured according to the grid cell node number and the time step order to obtain a motion state data set with time as the dimension and spatial nodes as the index. The motion state data set is subjected to consistency verification and format standardization processing, and the output is the discretized wave field time series of the three-dimensional computational domain throughout the entire simulation period.

8. The explicit finite element parallel computation method for seismic wave simulation according to claim 1, characterized in that, The domain decomposition of the three-dimensional computational domain, dividing the unstructured finite element mesh into multiple sub-regions, includes: Based on the element distribution and connection relationship of the unstructured finite element mesh, the three-dimensional computational domain is divided into several spatially continuous and computationally balanced sub-regions. The step of allocating the computational tasks of each sub-region to multiple parallel computation threads includes: Based on the division results of the sub-regions, the calculation tasks of the unstructured finite element mesh set and the associated discretized wave field time series corresponding to each sub-region are dynamically allocated to a corresponding parallel computing thread. The computational tasks assigned to several parallel computing threads are configured to independently perform the recursive solution of the discretized wavefield time series within the sub-regions assigned to them, based on the same explicit time integration format and lumped mass matrix, and to synchronize wavefield variables at the shared boundaries of adjacent sub-regions, thereby achieving collaborative parallel computation of the entire three-dimensional computational domain.

9. The explicit finite element parallel computation method for seismic wave simulation according to claim 8, characterized in that, The computational task, which is allocated to several parallel computing threads, is configured to independently perform the recursive solution of the discretized wavefield time series within its assigned sub-region based on the same explicit time integration format and lumped mass matrix, and to synchronize wavefield variables at the shared boundary of adjacent sub-regions, including: Configure each parallel computing thread with the same numerical solver kernel containing the explicit time integration format and lumped mass matrix, and load the unstructured finite element mesh set and initial field information corresponding to the sub-region to which it is assigned. Each parallel computing thread independently performs the recursive solution calculation of the discretized wavefield time series within its allocated sub-region, according to the explicit time integration format and lumped mass matrix. After each parallel computing thread completes the recursive solution for each time step, the shared boundary between its sub-region and adjacent sub-region is identified, and the current time step wave field variable of the grid cell node on the shared boundary is extracted. The wave field variables of the grid cell nodes on the shared boundary are exchanged between the parallel computing threads corresponding to the adjacent sub-regions through the parallel communication interface. Based on the boundary wave field variables obtained through exchange, the wave field states of several parallel computing threads at the shared boundary are updated for their own sub-regions, thereby achieving cross-sub-region wave field connection and global computation synchronization, and preparing for the recursive solution of the next time step.

10. The explicit finite element parallel computation method for seismic wave simulation according to claim 8, characterized in that, The step of dynamically allocating the computation task of the unstructured finite element mesh set corresponding to each sub-region and the associated discretized wave field time series to a corresponding parallel computing thread includes: Obtain the computational load of the unstructured finite element mesh set corresponding to each sub-region. The computational load is related to the size of the mesh set and the recursive solution complexity of the discretized wave field time series. Based on the real-time processing capability and communication topology of each parallel computing thread, with the goal of load balancing and minimizing communication overhead, each sub-region is dynamically mapped to a corresponding parallel computing thread. The geometric data, physical parameters, and initial conditions for the recursive solution of the unstructured finite element mesh set corresponding to each sub-region are sent to the mapped parallel computing thread to complete the dynamic allocation of the computing task.

Images