An underwater explosion load calculation method based on discontinuous Galerkin method and MPI

Abstract

The application relates to a large-scale underwater explosion problem, in particular to a calculation method of underwater explosion load based on an LDGM and MPI. The calculation method of underwater explosion load based on the LDGM and MPI comprises the following steps: performing non-structural grid division on a water area geometric model, calculating data in the non-structural grid after spatial discretization and time discretization according to a LDGM second-order wave equation, performing regional decomposition by using METIS during the calculation process, dividing the non-structural network into multiple subdomains and distributing the subdomains to different processors, exchanging numerical flux on the subdomain boundaries between the processors, performing parallel processing, until time domain loop calculation is completed, and obtaining a pressure field of a flow field evolving with time. The LDGM is used to solve the flow field pressure load, can more accurately capture the strong gradient of a shock wave and local details of a cavitation area, avoids numerical dissipation and oscillation problems, and improves calculation efficiency through a parallel processing scheme.

Description

Technical Field

[0001] This application relates to the problem of large-scale underwater explosions, specifically a method for calculating underwater explosion loads based on the discontinuous Galerkin method and MPI. Background Technology

[0002] Underwater explosions are a crucial research area in underwater engineering, marine military, and marine resource development. The process involves various complex physical phenomena, including the propagation of shock waves in water, the generation, expansion, and collapse of bubbles, and the interaction between the shock wave and surrounding structures. High-precision numerical simulations of underwater explosions face challenges such as drastic changes in physical properties (e.g., density fields) and complex boundary conditions, placing extremely high demands on the accuracy and stability of numerical methods. These challenges are specifically manifested in strong discontinuities (shock waves), large deformations (bubbles, splash), multiphase coupling (water, gas, explosives, structures), and cavitation. Despite years of research, accurate prediction and protection against underwater explosions remain extremely challenging.

[0003] Traditional simulation methods, primarily including the finite element method, finite difference method, and finite volume method, struggle to accurately capture the discontinuities of pressure in the flow field when solving underwater explosion dynamics problems. These methods may also exhibit numerical oscillations, leading to unstable solutions. Secondly, when handling large deformations and moving interfaces using Lagrangian or Eulerian meshes, these methods suffer from issues such as fuzzy interface capture and approximate interface geometry descriptions, resulting in insufficient accuracy or even computational termination. Finally, underwater explosion simulations typically span large spatiotemporal scales, often requiring highly refined meshes to accurately capture the pressure wave propagation characteristics of the fluid, posing significant challenges to computational resources and efficiency.

[0004] The Locally Discontinuous Galerkin Method (LDGM) is a high-precision numerical method capable of efficiently handling fluid dynamics problems with inclusion, strong discontinuity, and strong nonlinearity, particularly suitable for simulating shock wave propagation. Meanwhile, MPI (Message Passing Interface) parallel computing is a standard specification for a message passing interface. It defines a series of functions, constants, and behaviors for inter-process communication on distributed memory systems. Its core idea is that multiple processes (each with its own independent address space) collaborate to complete a large computational task by sending and receiving messages, primarily used in high-performance computing. Combining the LDGM and MPI parallel computing offers a possibility for solving underwater explosion problems. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the prior art and provide an underwater explosion load calculation method based on the discontinuous Galerkin method and MPI. This method discretizes the fluid medium of the far-field underwater explosion into an unstructured mesh, uses LDGM to solve the flow field pressure load, and can more accurately capture the strong gradient of the shock wave and the local details of the cavitation region by minimizing the local weighted residual. It avoids the numerical dissipation and oscillation problems that may occur in traditional methods, and improves the computational efficiency through a parallel processing scheme.

[0006] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: a method for calculating underwater explosive loads based on the discontinuous Galerkin method and MPI, comprising the following steps:

[0007] S1. Mesh Generation: Establish the geometric model of the underwater environment in the underwater explosion problem, and perform unstructured mesh generation. The computational domain is in a Cartesian coordinate system. O-xyz It was established in [the context], and discretized into [the following]. N Non-overlapping unstructured tetrahedral units Kj { j =1, 2, ..., N This yields a mesh file containing mesh composition information and node coordinate information.

[0008] S2. Subdomain division: Based on the mesh file in step S1, obtain the adjacency relationship of the mesh cells and generate a text file for the METIS software package to divide the subdomain.

[0009] S3. Cell Information Initialization: Based on the set number of subdomains, the gpmetis program is used to perform multi-level partitioning to obtain a region partitioning file containing cell affiliation attributes; the grid file generated in step S1 is read in parallel, and the nodes and cells of the subdomains are locally renumbered and duplicate nodes are deleted in each CPU core to ensure the uniqueness of shared nodes in local storage, so as to support continuous access to subdomain data.

[0010] S4.MPI Information Initialization: Establish the adjacent topological relationships of grid cells, cell surfaces, and grid nodes. Based on the topological structure of cells and nodes in adjacent subdomains, establish and store the relationships between every two adjacent subdomains, and store them continuously in the data structure established for non-blocking communication.

[0011] S5. Order Reduction Process: Introducing Auxiliary Variables The discrete second-order wave equation of LDGM is expressed in first-order form, as shown in equations (1) and (2).

[0012]

[0013] in, p Indicates hydrodynamic acoustic pressure, superscript The time derivative is given by the fluid velocity as... , and These are the bulk modulus and density of the fluid, respectively.

[0014] S6. Solving for approximate solutions Finding approximate solutions in piecewise polynomial space Such that for all test functions In the piecewise polynomial space, the relationships of equations (3) and (6) are satisfied respectively.

[0015]

[0016] in , and These are the normal vectors n in x -、 y -and z - Components in the direction, , , and These are numerical fluxes, which are discretized tetrahedral units. A single-valued function on the boundary.

[0017] S7. Introduce numerical flux to communicate unit information:

[0018] For the cell boundaries in the internal computational domain, the numerical flux calculation method is shown in equations (7) and (8).

[0019]

[0020] In the formula, Stability coefficient Auxiliary vector coefficients , It is the unit normal vector pointing out of the element plane. , , yes exist x- , y- , z- The components in the three directions, of which the components in equations (7) and (8) are... , , and The calculation formulas are shown in equations (9) to (12), respectively.

[0021]

[0022] In the formula, the superscript " "and" "" indicates a surface shared by two adjacent tetrahedral elements. s The variable above, where the superscript " ” Tetrahedron K variables, " ” express K The variables of adjacent tetrahedral elements.

[0023] For the element boundaries on the outer surface of the computational domain, the numerical flux is calculated as shown in equations (13)-(14).

[0024] (13)

[0025] (14)

[0026] In the formula, the fluid sound velocity is , and These are the bulk modulus and density of the fluid, respectively. p Represents hydrodynamic acoustic pressure, stability coefficient g N The vector gradient represents the normal gradient of the underwater explosion incident load on the flow field boundary, and c1 and a1 are impedance coefficients defined on the non-reflective boundary. Represents the acceleration of the fluid node.

[0027] S8. Parallel Processing: The data in each unstructured grid in step S1 is computed according to different computational unit types to ensure that each block of data is concentrated in the cache as much as possible, allowing the CPU core to continuously execute the same type of computation; data communication uses MPI_ISend and MPI_IRecv, and MPI_Waitall is used to ensure that all data communication has been completed.

[0028] S9. Time Discretization: A fourth-order Runge-Kutta method is used to update the variable values ​​in the time domain, outputting the required physical quantities in parallel. Each CPU core generates corresponding mesh information and flow field physical quantities, which are then post-processed. The size of the time step is determined according to... CFL The (Courant-Friedrichs-Levy) condition is used to determine this, as shown in equation (15).

[0029]

[0030] In the formula CFL = 0.01, It is the smallest feature length among all discretized units.

[0031] S10. Parallel solution of the wave equation: The data in the unstructured mesh of the computational domain after spatial and temporal discretization are calculated according to the LDGM formula in step S6. During the calculation, METIS in steps S2 and S3 is used to perform domain decomposition, dividing the unstructured network into multiple subdomains and assigning them to different processors. The processors exchange numerical fluxes on the boundaries of the subdomains through MPI in steps S4 and S8 to perform parallel processing until the time-domain cyclic calculation is completed, and the pressure field of the flow field evolving with time is obtained.

[0032] Preferably, the volume fraction in formula (3) in step S6 is calculated using the eight-point Gaussian integral formula. This calculation requires the unstructured tetrahedral elements in the physical coordinate system to be integrated. K ( Linearly mapped to isoparametric elements in the reference coordinate system ( The isoparametric elements have vertices at coordinates (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1) in the reference coordinate system.

[0033] The coordinate mapping relationship between the physical coordinate system and the reference coordinate system is given by formula (16):

[0034]

[0035] In the formula k It is the first in the computational domain Ω k Units, , , and It corresponds to the first k The coordinates of the four vertices of each unit in the physical coordinate system.

[0036] Preferably, the basis functions in formula (3) about x , y and z The partial derivatives are calculated using formula (17):

[0037]

[0038] Preferably, the flow field is not disturbed by the shock wave loading at the initial moment, and the initial conditions are as follows: .

[0039] Preferred acoustic dynamic pressure Linear artificial viscous pressure was added. The calculation formula is (18).

[0040]

[0041] In the formula, This is the damping coefficient, unless otherwise specified. Volumetric strain rate .

[0042] Preferably, in the MPI-based three-dimensional LDGM parallel program, the cell located in a certain subdomain is called a physical cell, and the adjacent cells of the boundary cell of the subdomain are called virtual cells. When calculating the numerical flux of the boundary cell of the subdomain, it is necessary to pass the information of the virtual cell to the CPU core where its adjacent cells are located. The information transfer between different CPU cores is accomplished by means of MPI data communication.

[0043] Preferably, the MPI communication mode includes aggregated communication and point-to-point communication, with the time step determined in step S9 for aggregated communication. Summing numerical errors and calculating the minimum feature length. L min It uses functions provided by MPI. MPI_Bcast , MPI_Allgather and MPI_Allreduc accomplish.

[0044] Preferably, in a 3D LDGM MPI parallel program, there are three types of memory allocation: single memory allocation, contiguous multiple memory allocations, and non-contiguous multiple memory allocations.

[0045] The single memory allocation works by accurately calculating the size of each array and allocating memory space that matches the size of the array.

[0046] Multiple consecutive memory allocations with the same data type and data structure, and with similar memory access times, can be pre-allocated into a memory pool, and small blocks of memory can be allocated consecutively from this pool when needed. During this process, some memory addresses can be reused, thus effectively making full use of memory.

[0047] Non-contiguous memory allocation can be done using two methods: static allocation and dynamic allocation. Static allocation is performed at compile time, and memory is allocated on the stack, which is suitable for small, fixed-size, and general-purpose arrays. Dynamic allocation allocates memory on the heap, and the memory it occupies is released in a timely manner when some parts of the array are no longer needed, which is suitable for large-scale arrays.

[0048] The beneficial effects of this invention are as follows:

[0049] 1. High precision advantage, overcoming the difficulties of accurately capturing the discontinuity and numerical oscillation of flow field pressure. This invention uses high-order polynomials to approximate the solution within the element, inherently possessing high precision and low numerical dissipation. This means it can capture the shock wave front generated by underwater explosions with extremely high clarity, without over-smoothing the wave front as in traditional methods, thus calculating the peak pressure and attenuation law of the shock wave more accurately. Furthermore, traditional methods often require the introduction of large artificial viscosity to stabilize the shock wave, which "contaminates" the entire flow field, especially in the back-wave region and cavitation region, leading to physical distortion. The required artificial viscosity is minimal, better preserving the physical authenticity of the flow field.

[0050] 2. Superior computing power and efficiency overcome a series of key challenges of traditional computing models, such as large computational load and low efficiency. LDGM, based on unstructured meshes, can easily discretize complex geometries. Combined with MPI parallel computing, it can decompose the computational task of tens of millions or even hundreds of millions of mesh cells onto tens of thousands of processor cores, making it possible to simulate the impact of a full-scale underwater explosion on a real ship. Furthermore, the computation of this invention is highly localized, with communication between cells strictly limited to directly adjacent cells. This gives it near-ideal parallel scalability on distributed memory systems (supercomputers). The computational scale can increase almost linearly with the number of processor cores, greatly reducing computation time.

[0051] 3. Enhance the realism of physical simulations. Underwater explosion shockwaves create a massive negative pressure zone, leading to cavitation. This invention provides high-precision simulation of the pressure field, offering more reliable input pressure data for cavitation cutoff models or other more complex cavitation models, thus enabling more accurate prediction of the generation, extent, and collapse of cavitation zones. Attached Figure Description

[0052] Figure 1 This is a schematic diagram of the isoparametric transformation of a tetrahedral element from the physical coordinate system to the reference coordinate system.

[0053] Figure 2 A flowchart of large-scale MPI parallel computation for 3D LDGM;

[0054] Figure 3 This is a schematic diagram of physical and virtual units within a subdomain.

[0055] Figure 4 This is a schematic diagram of the data storage structure for non-blocking communication.

[0056] Figure 5 A schematic diagram of the calculation model for spherical wave load propagation testing;

[0057] Figure 6 A comparison of pressure changes at measurement points calculated by 3D AFEM and parallel 3D LDGM;

[0058] Figure 7 A comparison of parallel speedup ratios for the spherical wave propagation problem under different numbers of CPU cores;

[0059] Figure 8 This is a comparison of the parallel efficiency of the spherical wave propagation problem under different numbers of CPU cores. Detailed Implementation

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

[0061] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this application described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0062] The present invention will now be further described with reference to the accompanying drawings. The following embodiments are only used to illustrate the structure of the present invention more clearly.

[0063] In this invention, the spatial discretization of LDGM includes three parts: complex geometric network decomposition, local higher-order approximation of cells, and introduction of numerical flux to communicate cell information.

[0064] The computational domain is in the rectangular coordinate system. O-xyz It was established in [the context], and discretized into [the following]. N Non-overlapping unstructured tetrahedral units Kj { j =1, 2, ..., N To use the LDGM discrete second-order wave equation, auxiliary variables are introduced. The wave equation is reformulated in first-order form.

[0065]

[0066] In the formula, p represents hydrodynamic acoustic pressure, and the superscript... The time derivative is given by the fluid velocity as... , and These are the bulk modulus and density of the fluid, respectively.

[0067] LDGM finds approximate solutions in a piecewise polynomial space. Such that for all test functions In piecewise polynomial space, the following relations hold:

[0068]

[0069] in , and These are the normal vectors n in x -、 y -and z - Components in the direction. , , and These are numerical fluxes, which are discretized tetrahedral units. A single-valued function on the boundary.

[0070] The volume fraction in formula (3) is calculated using the eight-point Gaussian integral formula, which requires the tetrahedral elements in the physical coordinate system to be integrated. K ( Linearly mapped to isoparametric elements in the reference coordinate system ( ),like Figure 1 As shown. The coordinates of the vertices of the isoparametric elements in the reference coordinate system are (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1), respectively. The coordinate mapping relationship between the physical coordinate system and the reference coordinate system is given by the following formula (16):

[0071]

[0072] In the formula k It is the first in the computational domain Ω k Units, , , and It corresponds to the first k The coordinates of the four vertices of each element in the physical coordinate system. The basis functions in equation (3) about x , y and z The partial derivatives can be calculated using the following formula:

[0073]

[0074] The introduction of numerical flux enables communication between adjacent cells, ensuring that physical quantities flowing out of one cell can completely flow into adjacent cells. For cell boundaries within the internal computational domain, the numerical flux is calculated as follows, a method that guarantees numerical stability.

[0075]

[0076] In the formula Stability coefficient Auxiliary vector coefficients , It is the unit normal vector pointing out of the element plane. , , yes exist x- , y- , z- The components in three directions. Among them, in equations (7) and (8) , , and The calculation formulas are given by equations (9) to (12), respectively.

[0077]

[0078] The superscript " "and" "" indicates a surface shared by two adjacent tetrahedral elements. s Variables on, "Represents a tetrahedron" K variables, "express K The variables of adjacent tetrahedral elements.

[0079] For the element boundaries on the outer surface of the computational domain, their numerical flux depends on the physical boundary conditions.

[0080]

[0081] In the formula, the fluid sound velocity is , and These are the bulk modulus and density of the fluid, respectively. p Represents hydrodynamic acoustic pressure, stability coefficient g N The vector gradient represents the normal gradient of the underwater explosion incident load on the flow field boundary, and c1 and a1 are impedance coefficients defined on the non-reflective boundary. Represents the acceleration of the fluid node.

[0082] For the time discretization of LDGM, a fourth-order Runge-Kutta method is used to update the variable values ​​in the time domain. The size of the time step can be determined according to... CFL The (Courant-Friedrichs-Levy) condition is used to determine this.

[0083]

[0084] In the formula CFL = 0.01, It is the smallest feature length among all discretized elements. For tetrahedral elements, , V Indicates the volume of a unit cell. CFL = 0.01 is a very small value, designed to ensure the absolute stability of the algorithm when solving the problem of strong discontinuous shock waves.

[0085] In this invention, the flow field is initially unaffected by the shock wave loading, and the initial conditions are as follows: This initial condition provides a clean, interference-free starting point, with the flow field initially calm and uniform. After the simulation begins, the shock wave, introduced through boundary conditions and the source phase, clearly demonstrates its generation, propagation, and the interaction of the flow fields. The time-domain calculation loop is as follows: initialization (setting the initial flow field) -- determining the time step -- performing time integration -- applying physical boundary conditions to the updated flow field, repeating the steps until the preset simulation time is reached.

[0086] The CFL (Courant-Friedrichs-Levy) conditions applied in this section can automatically determine the maximum stable time step. The simulation recalculates the CFL conditions at each step and dynamically adjusts them. This ensures that the calculations are performed efficiently and stably.

[0087] To address the strong pressure discontinuities caused by shock waves and ensure the stability of the LDGM, acoustic dynamic pressure... Linear artificial viscous pressure was added. ,

[0088]

[0089] In the formula It is the damping coefficient, volumetric strain rate. In this invention, unless otherwise specified, .

[0090] The core feature of LDGM is that it allows for discontinuities in solutions between any two adjacent elements. The local polynomial solution of each element is computed relatively independently, depending only on its neighboring elements, and coupled through numerical flux. Therefore, LDGM exhibits high parallelizability and is suitable for high-performance computing. In this invention, the MPI parallel processing flow is as follows: Figure 2 As shown, the complete scheme consists of three parts: cell information initialization, MPI information initialization, and parallel solution of 3D LDGM: (1) In the cell information initialization module, the cells of the entire computational domain are almost uniformly divided into different subdomains. The cells and nodes of each subdomain are renumbered to facilitate indexing and appropriate memory allocation. (2) In the MPI information initialization module, the relationship between each two adjacent subdomains is established and stored according to the topology of cells and nodes in adjacent subdomains, and stored continuously in the data structure established for non-blocking communication. (3) In the parallel solution module of 3D LDGM, a parallel 3D LDGM model is mainly used to solve the second-order wave equation satisfied by the underwater explosion shock wave, and then the calculation results are output in parallel.

[0091] This invention employs the gpmetis program from the METIS software package to handle the region partitioning of large-scale unstructured tetrahedral elements. In the MPI-based three-dimensional LDGM parallel program, an element located within a subdomain is called a physical element, and the adjacent elements of the boundary element of that subdomain (located within other subdomains) are called virtual elements. Figure 3 A schematic diagram of physical and virtual cells is provided, with virtual cells marked in green. For LDGM, the computation of each physical cell is influenced by information from its neighboring cells, which is coupled through the numerical flux of the cell boundaries. Specifically, when calculating the numerical flux of a subdomain boundary cell, the information of the virtual cell needs to be passed to the CPU core containing its neighboring cells.

[0092] Information transfer between different CPU cores is accomplished through MPI data communication. Data communication is a crucial part of the design of large-scale MPI parallel algorithms, directly affecting the performance and efficiency of the entire program. MPI communication modes mainly include two types: aggregate communication and point-to-point communication. In this invention, aggregate communication is only used in necessary situations, such as determining the time step in equation (9). The numerical errors are summed and the minimum feature length is calculated. L min Operations such as those performed using functions provided by MPI, such as... MPI_Bcast , MPI_Allgather and MPI_Allreduc To achieve this, point-to-point communication employs [a specific method / mechanism]. To prevent deadlocks during computation and communication, point-to-point communication is implemented using [a specific method / mechanism]. MPI_ISend and MPI_IRecvTo achieve non-blocking communication, before the numerical flux calculation begins, a method is used... MPI_Waitall This is to ensure that all data communication has been completed.

[0093] Since sequential memory access can reduce latency, this invention designs as follows: Figure 4 The access-friendly data structure shown ensures efficient concurrency of multiple non-blocking communications.

[0094] Simulating large-scale underwater explosion problems using a 3D LDGM model discretizes the target computational domain into a large number of mesh elements, introducing numerous degrees of freedom and resulting in a massive computational workload. Furthermore, the computer requires significant memory to store related variable information. When numerous processes and threads operate concurrently, the time consumed by memory accesses impacts the overall computational efficiency of the MPI program. Typically, without any memory management strategies and optimization techniques, a parallel 3D LDGM model for underwater explosion problems can only handle a maximum of a few million elements. To scale up computation and further improve parallel efficiency, effective memory management strategies need to be designed within the MPI parallel framework. For MPI parallel programs, the key optimizations in memory management are latency and memory utilization. These two factors significantly influence the efficiency and scalability of MPI programs, especially when dealing with large-scale computational problems. To achieve efficient memory management, reducing latency and limiting memory utilization must be considered in MPI parallel programs.

[0095] In the MPI parallel program of 3D LDGM, non-contiguous data access is mainly caused by two aspects: (1) unoptimized data access order; (2) poor data locality. The data access order largely depends on the numbering of cells and nodes, but after the region is partitioned by the METIS algorithm, the cells and nodes allocated to each CPU core are not arranged in order. Therefore, a renumbering strategy for nodes and cells needs to be designed in each CPU core to support continuous data access. It should be noted that duplicate nodes need to be removed before renumbering. In addition, the MPI program of this invention is written in FORTRAN language and adopts column-major order storage. Therefore, in order to ensure that the accessed data is contiguous in memory, the array is traversed column by column. Data locality is divided into temporal locality and spatial locality. Temporal locality means that recently accessed data may be accessed again soon. Since caching systems usually retain recently accessed data, temporal locality is mainly enhanced by cache optimization, that is, reusing data already loaded into the cache as much as possible. This invention improves cache hit rate and enhances temporal data locality by increasing loop intervals and reducing iterations where appropriate to execute multiple instructions within a single loop frame. This approach is crucial when processing large-scale multidimensional arrays, especially in the transfer of cell information involved in data communication. Spatial locality means that when a CPU core accesses data at a certain location, nearby data is also likely to be accessed. In the MPI parallel program of 3D LDGM, the computational workload of cells within a subdomain, cells at the subdomain boundary, and cells of periodicity is not equivalent. Dividing large-scale computational tasks into blocks according to different computational cell types ensures that each block of data is concentrated in the cache as much as possible, allowing CPU cores to execute the same type of computational tasks consecutively, thereby further enhancing spatial data locality. Overall, enhanced data locality increases the probability of consecutive data access and reduces latency, thereby further improving the parallel efficiency of MPI programs. Furthermore, the supercomputing platform used in this invention is based on a Non-Unified Memory Access (NUMA) architecture, where each computing node has its own independent CPU core and memory, and the memory is evenly distributed among the CPU cores, further reducing hardware-level memory access latency.

[0096] In the MPI parallel program for 3D LDGM, there are three types of memory allocation: single memory allocation, contiguous memory allocation, and non-contiguous memory allocation. For large-scale computational problems, any improper handling of arrays can lead to poor memory utilization. The sharp increase in memory allocation can affect the normal execution of programs. Therefore, this invention designs optimization strategies for each type of memory allocation.

[0097] (1) For a single memory allocation, it is necessary to accurately calculate the size of each array and allocate memory space that matches the array size to avoid unnecessary memory allocation. At the same time, determining the size of the array to be used before the loop calculation begins can save the time consumed by subsequent data indexing operations.

[0098] (2) For multiple consecutive memory allocations, they usually have the same data type and data structure, and the memory access time is relatively close. Therefore, a memory pool can be pre-allocated, and small blocks of memory can be allocated consecutively from the memory pool when needed. In this process, some memory addresses can be reused, thereby effectively making full use of memory. This method also avoids frequent memory allocation, thereby improving the overall performance of the program.

[0099] (3) For multiple non-contiguous memory allocations, this invention considers both static and dynamic allocation methods. Static allocation is performed at compile time, and memory is allocated on the stack. This method can achieve faster access speed and is suitable for small, fixed-size, and general-purpose arrays. For large-scale arrays, this invention uses dynamic allocation, allocating memory on the heap, and the memory occupied by some arrays can be released in a timely manner when some arrays are no longer needed.

[0100] In summary, this invention employs the three strategies described above—optimizing memory allocation size, reusing allocated memory addresses, and flexible dynamic allocation—to effectively mitigate excessive increases in memory usage. Therefore, the MPI-parallel 3D LDGM program of this invention can handle simulations of several large-scale underwater explosion problems.

[0101] The invention includes the following steps:

[0102] S1. Mesh Generation: Establish the geometric model of the water area in the underwater explosion problem, perform unstructured mesh generation, and obtain a mesh file containing mesh composition information and node coordinate information.

[0103] S2. Subdomain division: Based on the mesh file in step S1, obtain the adjacency relationship of the mesh cells and generate a text file for the METIS software package to divide the subdomain.

[0104] S3. Cell Information Initialization: Based on the set number of subdomains, the gpmetis program is used to perform multi-level partitioning to obtain a region partitioning file containing cell affiliation attributes; the mesh file generated in step S1 is read in parallel, and the nodes and cells are renumbered and duplicate nodes are deleted in each CPU core to support continuous access to subdomain data.

[0105] S4.MPI Information Initialization: Establish the adjacent topological relationships of grid cells, cell surfaces, and grid nodes. Based on the topological structure of cells and nodes in adjacent subdomains, establish and store the relationships between every two adjacent subdomains, and store them continuously in the data structure established for non-blocking communication.

[0106] S5. Order Reduction Process: Introducing Auxiliary Variables The discrete second-order wave equation of LDGM is expressed in first-order form.

[0107] S6. Solving for approximate solutions Finding approximate solutions in piecewise polynomial space .

[0108] S7. Introduce numerical flux to communicate unit information: Specify the calculation methods for numerical flux at unit boundaries within the internal computational domain and on the outer surface of the computational domain.

[0109] S8. Parallel Processing: The data in each unstructured grid in step S1 is computed according to different computational unit types to ensure that each block of data is concentrated in the cache as much as possible, allowing the CPU core to continuously execute the same type of computation; data communication uses MPI_ISend and MPI_IRecv, and MPI_Waitall is used to ensure that all data communication has been completed.

[0110] S9. Time Discretization: A fourth-order Runge-Kutta method is used to update the variable values ​​in the time domain, outputting the required physical quantities in parallel. Each CPU core generates corresponding mesh information and flow field physical quantities, which are then post-processed. The size of the time step is determined according to... CFL The (Courant-Friedrichs-Levy) condition is used to determine this.

[0111] S10. Parallel solution of the wave equation: The data in the unstructured mesh of the computational domain after spatial and temporal discretization are calculated according to the LDGM formula in step S6. During the calculation, METIS in steps S2 and S3 is used to perform domain decomposition, dividing the unstructured network into multiple subdomains and assigning them to different processors. The processors exchange numerical fluxes on the boundaries of the subdomains through MPI in steps S4 and S8 to perform parallel processing until the time-domain cyclic calculation is completed, and the pressure field of the flow field evolving with time is obtained.

[0112] Example 1: To verify the accuracy of the parallel three-dimensional LDGM model in solving underwater explosion shock wave problems, this invention simulates the propagation of spherical wave loads within a watershed. The watershed size is set to... The upper surface of the watershed is set as a free surface with an atmospheric pressure of 0.101 MPa. An incident load boundary condition is applied to the lower surface, while other surfaces are set as non-reflective boundary conditions. The locations of the burst point and the source point are (0,0,0) and (0,0, -0.2 m), respectively. Throughout the simulation, the incident spherical wave loading pressure amplitude at the source point is 0.1 MPa. The bulk modulus of the fluid is... Kf and density 2.14 GPa and 1000 kg / m³ respectively 3 A schematic diagram of the computational model is shown below. Figure 5 As shown.

[0113] Based on the same simulation setup, this invention uses the three-dimensional AFEM provided in ABAQUS software and the parallel three-dimensional LDGM proposed in this invention for simulation. Figure 6 The comparison of dynamic pressure histories obtained by 3D AFEM and parallel 3D LDGM at measurement points (1m, 0, 1m) is presented. It can be seen that the pressure results obtained by the parallel 3D LDGM model are highly consistent with those calculated by 3D AFEM in ABAQUS. Therefore, it can be proven that the current parallel 3D LDGM model is sufficiently accurate to simulate the propagation of underwater explosion shock wave loads within the watershed.

[0114] To further test the parallel performance of the parallel 3D LDGM model, this invention simulates the problem using different numbers of CPU cores, ranging from 60 to 7680, while keeping the total computational load constant (strong scaling test). To expand the test scope, the computational domain was refined to 40 million tetrahedral elements. Figure 7 and Figure 8 The parallel speedup and parallel efficiency are shown for different numbers of CPU cores. It can be seen that when the number of CPU cores is less than 2000, the speedup of the parallel program exceeds the ideal value, and the parallel efficiency is greater than 100%, achieving superlinear speedup. When the number of CPU cores reaches 7680, the parallel efficiency remains at 98.3%, indicating that the parallel 3D LDGM model has high scalability. The good parallel performance and scalability demonstrate the great potential of this model in simulating large-scale 3D underwater explosion impact problems.

[0115] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the technical principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for calculating underwater explosive loads based on the discontinuous Galerkin method and MPI, characterized in that: Includes the following steps: S1. Mesh Generation: Establish the geometric model of the underwater environment in the underwater explosion problem. The computational domain is in a Cartesian coordinate system. O-xyz It was established in [the context], and discretized into [the following]. N Non-overlapping unstructured tetrahedral units K j { j =1, 2, ..., N This yields a mesh file containing mesh composition information and node coordinate information. S2. Subdomain division: Based on the mesh file in step S1, obtain the adjacency relationship of the mesh cells and generate a text file for the METIS software package to divide the subdomain; S3. Cell Information Initialization: Based on the set number of subdomains, the gpmetis program is used to perform multi-level partitioning to obtain a region partitioning file containing cell affiliation attributes; the mesh file generated in step S1 is read in parallel, and the nodes and cells of the subdomains are locally renumbered and duplicate nodes are deleted in each CPU core to ensure the uniqueness of shared nodes in local storage, so as to support continuous access to subdomain data. S4.MPI Information Initialization: Establish the adjacent topological relationships of grid cells, cell surfaces, and grid nodes. Based on the topological structure of cells and nodes in adjacent subdomains, establish and store the relationships between every two adjacent subdomains, and store them continuously in the data structure established for non-blocking communication. S5. Order Reduction Process: Introducing Auxiliary Variables The discrete second-order wave equation of LDGM is expressed in first-order form, as shown in equations (1) and (2). in, p Indicates hydrodynamic acoustic pressure, superscript The time derivative is given by the fluid velocity as... , and These are the bulk modulus and density of the fluid, respectively. S6. Solving for an approximate solution Finding approximate solutions in piecewise polynomial space This makes it possible for all basis functions In the piecewise polynomial space, the relationships of equations (3) and (6) are satisfied respectively. in, , and These are the normal vectors n in x -、 y -and z - Components in the direction, , , and These are numerical fluxes, which are discretized tetrahedral units. A single-valued function on the boundary, the fluid sound velocity is , and These are the bulk modulus and density of the fluid, respectively. S7. Introducing numerical flux: For the cell boundaries in the internal computational domain, the numerical flux calculation method is shown in equations (7) and (8). In the formula, Stability coefficient Auxiliary vector coefficients , It is the unit normal vector pointing out of the element plane. , , yes exist x- , y- , z- The components in the three directions, of which the components in equations (7) and (8) are... , , and The calculation formulas are given by equations (9) to (12), respectively. In the formula, the superscript " "and" "" represents a variable on surface s shared by two adjacent tetrahedral elements, where the superscript " ” Tetrahedron K variables, " ” express K The variables of adjacent tetrahedral elements; For the element boundaries on the outer surface of the computational domain, the numerical flux is calculated as shown in equations (13)-(14). （13） （14）； In the formula, the fluid sound velocity is , and These are the bulk modulus and density of the fluid, respectively. p Represents hydrodynamic acoustic pressure, stability coefficient g N The vector gradient represents the normal gradient of the underwater explosion incident load on the flow field boundary, and c1 and a1 are impedance coefficients defined on the non-reflective boundary. The acceleration representing the fluid node; S8. Parallel Processing: The data in each unstructured grid in step S1 is computed according to different computational unit types, ensuring that each block of data is concentrated in the cache, allowing the CPU core to continuously execute the same type of computation; data communication uses MPI_ISend and MPI_IRecv, and MPI_Waitall is used to ensure that all data communication has been completed; S9. Time Discretization: A fourth-order Runge-Kutta method is used to update the variable values ​​in the time domain, and the required physical quantities are output in parallel. Each CPU core generates corresponding mesh information and flow field physical quantities, which are then post-processed. The size of the time step is determined according to... CFL The conditions are determined, as shown in equation (15). In the formula CFL = 0.01, It is the smallest characteristic length among all discretized units, and the fluid sound velocity is... , and These are the bulk modulus and density of the fluid, respectively. S10. The data in the unstructured mesh in the computational domain after spatial and temporal discretization are calculated according to the LDGM formula in step S6. During the calculation, METIS in steps S2 and S3 is used to perform domain decomposition, dividing the unstructured mesh into multiple subdomains and assigning them to different processors. The processors exchange numerical fluxes on the subdomain boundaries through MPI in steps S4 and S8 to perform parallel processing until the time-domain cyclic calculation is completed, and the pressure field of the flow field evolving with time is obtained.

2. The underwater explosive load calculation method based on the discontinuous Galerkin method and MPI as described in claim 1, characterized in that: The volume fraction in formula (3) in step S6 is calculated using the eight-point Gaussian integral formula. This calculation requires the non-structured tetrahedral elements in the physical coordinate system to be converted into their components. K ( Linearly mapped to isoparametric elements in the reference coordinate system ( The coordinates of the vertices of the isoparametric elements in the reference coordinate system are (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1), respectively. The coordinate mapping relationship between the physical coordinate system and the reference coordinate system is given by formula (16): In the formula k It is the first in the computational domain Ω k Units, , , and It corresponds to the first k The coordinates of the four vertices of each unit in the physical coordinate system.

3. The underwater explosion load calculation method based on the discontinuous Galerkin method and MPI as described in claim 2, characterized in that: The basis functions in formula (3) about x , y and z The partial derivatives are calculated using formula (17): 。 4. The underwater explosive load calculation method based on discontinuous Galerkin method and MPI as described in claim 1, characterized in that: The flow field is initially unaffected by the shock wave loading, and the initial conditions are as follows: .

5. The underwater explosive load calculation method based on discontinuous Galerkin method and MPI as described in claim 1, characterized in that: Linear artificial viscous pressure was added to the acoustic dynamic pressure. The calculation formula is (18). In the formula, It is the damping coefficient. Volumetric strain rate The speed of sound in the fluid is , and These are the bulk modulus and density of the fluid, respectively.

6. The underwater explosive load calculation method based on the discontinuous Galerkin method and MPI as described in claim 1, characterized in that: In the MPI-based 3D LDGM parallel program, the cell located in a subdomain is called a physical cell, and the adjacent cells of the boundary cell of the subdomain are called virtual cells. When calculating the numerical flux of the boundary cell of the subdomain, the information of the virtual cell needs to be transmitted to the CPU core where its adjacent cells are located. The information transmission between different CPU cores is accomplished by MPI data communication.

7. The underwater explosive load calculation method based on the discontinuous Galerkin method and MPI as described in claim 1, characterized in that: MPI communication modes include aggregated communication and point-to-point communication. Aggregated communication is used to determine the time step in step S9. Summing numerical errors and calculating the minimum feature length. L min Each through the functions provided by MPI MPI_Bcast , MPI_Allgather and MPI_Allreduce accomplish.

8. The underwater explosive load calculation method based on the discontinuous Galerkin method and MPI as described in claim 1, characterized in that: In the MPI parallel program of 3D LDGM, there are three types of memory allocation: single memory allocation, contiguous multiple memory allocations, and non-contiguous multiple memory allocations. A single memory allocation works by accurately calculating the size of each array and allocating memory space to it that matches the array size. Multiple consecutive memory allocations have the same data type and data structure, and the memory access time is close. A memory pool is pre-allocated, and small blocks of memory are continuously allocated from the memory pool when needed. During the process, some memory addresses are reused, thus making full use of memory. Non-contiguous memory allocation can be done using two methods: static allocation and dynamic allocation. Static allocation is performed at compile time, and memory is allocated on the stack, which is suitable for small, fixed-size, and general-purpose arrays. Dynamic allocation allocates memory on the heap, and the memory it occupies is released in a timely manner when some parts of the array are no longer needed, which is suitable for large-scale arrays.

Images