Light ray tracing based particle-like numerical simulation acceleration method and system

Abstract

The application discloses a particle numerical simulation acceleration method based on ray tracing, and the method is accelerated by a GPU-based ray tracing core, and the method comprises the following steps: expanding a search range of a to-be-searched particle into a neighborhood search area bounding box, wherein the neighborhood search area bounding box comprises a plurality of surrounding particles located around the to-be-searched particle; calling the ray tracing core of the GPU to emit rays from the centers of the surrounding particles respectively, and detecting whether the rays intersect with the neighborhood search area bounding box; and calling an intersection program to update a neighbor list when the rays intersect with the neighborhood search area bounding box. Therefore, the method can apply the ray tracing core to particle numerical simulation methods such as a discrete element method, a molecular dynamics method, a near-field dynamics method and a smoothed particle hydrodynamics method, can fully utilize GPU computing resources, can improve the computing efficiency of the particle numerical simulation method, and can promote the high-performance computing field of related industrial software.

Description

Technical Field

[0001] This invention relates to the field of high-performance computing technology for industrial software, and in particular to a method and system for accelerating particle numerical simulations based on ray tracing. Background Technology

[0002] Particle-based numerical methods (or particle-class numerical methods), such as discrete element methods, molecular dynamics, peridynamics, and smoothed particle hydrodynamics, are gaining popularity in research and engineering due to their meshless and highly parallelizable properties. In particle-based numerical methods, the motion of a particle depends on its interaction with surrounding particles (neighbors). Searching for the neighbors of each particle is a crucial step, called neighborhood search; it is computationally very expensive, often consuming a large portion of the simulation time. Notably, particle positions can change significantly and rapidly, especially for dynamic problems, making the neighborhood search computationally intensive. For example, contact detection and computation in discrete element methods (where neighborhood search computation dominates for spherical particles) can account for up to 80% of the runtime. Therefore, developing new acceleration algorithms for particle-class numerical methods to make them more practical for engineering-scale simulations is both important and urgent.

[0003] Currently, there has been considerable research on neighborhood search algorithms in numerical methods for particle systems. Given an N-particle system, such as... Figure 1 As shown in (a), neighborhood search can be implemented exhaustively, that is, traversing all the remaining N-1 particles to check if they are neighbors of a given particle. The exhaustive method is very simple, but its time complexity is O(N). 2 Therefore, it is rarely used in practical simulations. A popular alternative is the so-called sweep-trimming algorithm. This algorithm requires sorting the lower and upper bounds of the particle's axis-aligned bounding box on different coordinate axes. If the particle's axis-aligned bounding box overlaps with the given particle's axis-aligned bounding box on all axes, then the particle can be considered a neighbor particle of the given particle. Sweep-trimming neighborhood search is also effective even for moving particles, and its implementation can be easily found in some open-source program code, such as the discrete element open-source programs YADE and SudoDEM. In addition, neighborhood search also often employs two other popular algorithms: one is a cell-based method, such as... Figure 1 (b) shows an alternative tree-based approach, especially... Figure 1(c) shows the hierarchical bounding boxes. Cell-based neighborhood search relies primarily on a uniform grid, through which the mapping relationship between particles and grid cells can be efficiently obtained. Therefore, by traversing only the cells of a given particle and the particles in its neighboring cells, the search for neighboring particles throughout the entire domain can be significantly reduced. For hierarchical bounding box-based neighborhood search, a hierarchical bounding box tree needs to be constructed first by aligning bounding boxes with all axes of the particle, and then the neighbors of a given particle are found by traversing the corresponding branches and nodes of the hierarchical bounding box tree.

[0004] With the rapid development of general-purpose graphics processing unit (GPU) computing, modern GPUs can provide tens of thousands of physically concurrent threads for scientific computing, not just for initial image or video rendering purposes. For example, NVIDIA's CUDA (Compute Unified Device Architecture) platform provides a scalable programming model for GPU computing with hierarchical organization of threads. Such a rich set of GPU threads undoubtedly accelerates particle-based numerical simulations; for instance, one thread can be used for each particle's computation. In fact, there is already a large body of research and application of CUDA-based GPU-accelerated algorithms for accelerating particle-based simulations: discrete element method (DEM), molecular dynamics, peri-field dynamics, and smoothed particle hydrodynamics. To fully utilize the parallel computing capabilities of GPUs, neighborhood searches in these particle-based numerical methods are typically implemented using cell-based or BVH (Bounding Volume Hierarchy) algorithms, while sweep pruning algorithms are difficult to implement efficiently on GPUs.

[0005] With the rapid development of computer graphics and metaverse technology, NVIDIA has designed additional new hardware specifically for ray tracing rendering, complementing the CUDA cores. Ray tracing cores are becoming increasingly powerful to meet the growing rendering demands. While the aforementioned particle-based numerical methods (discrete element method, molecular dynamics, peridynamics, and smoothed particle hydrodynamics) offer significant speedups on GPUs compared to CPU computing systems, all applications are limited to using the CUDA cores.

[0006] Therefore, a new solution is needed. Summary of the Invention

[0007] The main objective of this invention is to provide a method and system for accelerating particle numerical simulations based on ray tracing.

[0008] To achieve the above objectives, this invention provides a method for accelerating particle-based numerical simulations based on ray tracing, using hardware acceleration based on the ray tracing core of a GPU, comprising the following steps:

[0009] The search range of the particle to be searched is expanded to a neighborhood search region bounding box, which includes multiple surrounding particles located around the particle to be searched.

[0010] The GPU's ray tracing core is invoked to emit rays from the centers of the surrounding particles, and the intersection of these rays with the bounding box of the neighboring search region is detected.

[0011] When the ray intersects with the bounding box of the neighborhood search region, the intersection procedure is invoked to update the neighbor list.

[0012] In the particle numerical simulation acceleration method based on ray tracing provided by this invention, the step of calling the intersection procedure to update the neighbor list includes:

[0013] Compare the ray number i of the intersecting rays with the particle number j of the particle to be searched;

[0014] If particle number j is greater than ray number i, then update the neighbor list NL[i*offset+N1]=j; if particle number j is less than ray number i, then update the neighbor list NL[(i+1)*offset-N2-1]=i, where offset is the maximum number of neighbors for each particle, and N1 and N2 are the number of neighbors whose particle numbers are greater than and less than ray numbers, respectively.

[0015] In the particle numerical simulation acceleration method based on ray tracing provided by this invention, the step of calling the intersection procedure to update the neighbor list includes:

[0016] Compare the ray number i of the intersecting rays with the particle number j of the particle to be searched;

[0017] If particle number j is less than ray number i, then update the neighbor list NL[i*offset+N1]=j. If particle number j is greater than ray number i, then update the neighbor list NL[(i+1)*offset-N2-1]=i. Here, offset is the maximum number of neighbors for each particle, and N1 and N2 are the number of neighbors whose particle numbers are less than and greater than the ray number, respectively.

[0018] In the particle numerical simulation acceleration method based on ray tracing provided by this invention, before comparing the ray number i of the intersecting rays with the particle number j of the particle to be searched, the method further includes:

[0019] Calculate the distance between the center position of the surrounding particles corresponding to the intersecting rays and the center position of the particle to be searched; and

[0020] Determine whether the distance is less than a preset threshold, where the preset threshold is (1+δ)l, where δ is the relaxation size factor and l is the generalized search radius.

[0021] Furthermore, to achieve the above objectives, this invention also provides a particle-based numerical simulation acceleration system based on ray tracing, which uses GPU ray tracing cores for hardware acceleration, including:

[0022] An extension module is used to extend the search range of the particle to be searched to a neighborhood search region bounding box, wherein the neighborhood search region bounding box includes multiple surrounding particles located around the particle to be searched;

[0023] The detection module is used to invoke the GPU's ray tracing core to emit rays from the centers of the multiple surrounding particles, and to detect whether the rays intersect with the bounding box of the neighboring search region; and

[0024] An update module is used to call an intersection procedure to update the neighbor list when the ray intersects with the bounding box of the neighborhood search region.

[0025] In the ray tracing-based particle numerical simulation acceleration system provided by this invention, the update module includes:

[0026] A comparison unit is used to compare the ray number i of the intersecting rays with the particle number j of the particle to be searched;

[0027] The update unit is used to update the neighbor list NL[i*offset+N1]=j when the particle number j is greater than the ray number i, and to update the neighbor list NL[(i+1)*offset-N2-1]=i when the particle number j is less than the ray number i, where offset is the maximum number of neighbors for each particle, and N1 and N2 are the number of neighbors whose particle numbers are greater than and less than the ray number, respectively.

[0028] In the ray tracing-based particle numerical simulation acceleration system provided by this invention, the update module includes:

[0029] A comparison unit is used to compare the ray number i of the intersecting rays with the particle number j of the particle to be searched;

[0030] The update unit is used to update the neighbor list NL[i*offset+N1]=j when the particle number j is less than the ray number i, and to update the neighbor list NL[(i+1)*offset-N2-1]=i when the particle number j is greater than the ray number i, where offset is the maximum number of neighbors for each particle, and N1 and N2 are the number of neighbors whose particle numbers are less than and greater than the ray number, respectively.

[0031] In the ray-tracing-based particle numerical simulation acceleration system provided by this invention, the update module further includes:

[0032] The prediction unit is used to calculate the distance between the center position of the surrounding particles corresponding to the intersecting rays and the center position of the particle to be searched, and to determine whether the distance is less than a preset threshold, wherein the preset threshold is (1+δ)l, where δ is the relaxation size factor and l is the generalized search radius.

[0033] The present invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the steps of the ray-tracing-based particle numerical simulation acceleration method described above.

[0034] The ray tracing-based particle numerical simulation acceleration system and method provided by this invention have the following beneficial effects: The ray tracing-based particle numerical simulation acceleration method provided by this invention expands the search range of the particle to be searched to a neighborhood search region bounding box; it calls the ray tracing core of the GPU to emit rays from the centers of the multiple surrounding particles respectively, and detects whether the rays intersect with the neighborhood search region bounding box; when the rays intersect with the neighborhood search region bounding box, it calls the intersection procedure to update the neighbor list; thus, this method uses the ray tracing core to perform particle neighborhood search, and then performs this calculation process in parallel with the particle kinematics solution, thereby making full use of GPU computing resources, improving the computational efficiency of particle numerical simulation, and promoting the field of high-performance computing in related industrial software. Attached Figure Description

[0035] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort:

[0036] Figure 1 The diagram shows a neighborhood search algorithm for existing particle numerical methods; where (a) a particle is aligned with an axis-bounded box, (b) a particle is aligned with a uniform grid, and (c) a particle is aligned with a hierarchical bounding box.

[0037] Figure 2The diagram shows bounding boxes for neighborhood search problems using various particle numerical methods, where: (a) particles with uniform axis-aligned bounding boxes in the Discrete Element Method (DEM) or Molecular Dynamics Method (MDL); (b) the expanded axis-aligned bounding box (large rectangle) of the search region (small rectangle) in the DEM or MDL; (c) particle A, the influence range of particle A (circle), and its axis-aligned bounding box (rectangle) in the peridynamics or smooth hydrodynamics Method (PMD); and (d) the expanded axis-aligned bounding box (large rectangle) of the influence range of particle A in the peridynamics or smooth hydrodynamics Method (PMD).

[0038] Figure 3 The main idea of ​​ray tracing is shown, where (a) is a schematic diagram of ray tracing; and (b) is the hierarchical bounding box acceleration for ray collision detection.

[0039] Figure 4 (a) NVIDIA ray tracing hardware core and (b) the general ray tracing workflow of OptiX 7;

[0040] Figure 5 The diagram shows: (a) the search range of particle A and the corresponding bounding box aligned with the extended axis; (b) the bounding box aligned with the extended axis of particle A and the center (point) of potential neighboring particles; (c) the point-bounding box intersection as a ray tracing problem; and (d) further solving the contact problem of potential contacting particles.

[0041] Figure 6 The image shown is a snapshot of a simulated collapse column of 800,000 particles.

[0042] Figure 7 The figure shows the change in average computation speed over time for a collapse column simulation using different bounding box magnification factors (800,000 particles); the speedup ratio is the ratio of the computation speed using the method of this invention to the computation speed of the traditional cell method (Quadro RTX 5000 GPU).

[0043] Figure 8 The image shows a snapshot of a simulated collapse column with 100,000 and 1.6 million particles.

[0044] Figure 9 The figure shows the evolution of the average computation speed of the collapsed column simulation as the number of particles changes; the speedup ratio is the ratio of the computation speed using the method of this invention to the computation speed of the traditional cell method (Quadro RTX 5000 GPU).

[0045] Figure 10 The following are examples: (a) the evolution of the reduction in computational speed of the collapse column simulation as the number of particles increases exponentially, and (b) the memory usage of the collapse column simulation with different particle numbers (Quadro RTX 5000 GPU).

[0046] Figure 11 The image shown is a snapshot of a simulation of 1.6 million particles stacked together.

[0047] Figure 12 The following figures illustrate the evolution of computational speed on different GPUs (1.6 million particles) for (a) particle stacking and (b) collapsing column simulations.

[0048] Figure 13 The following are simulations of particle stacking and collapsing columns on different GPUs: (a) final simulation computation speed using the inventive method and (b) evolution of the speed ratio (1.6 million particles);

[0049] Figure 14 The following are the debris flow simulation settings: (a) terrain and slip source, and (b) a locally enlarged slip source;

[0050] Figure 15 The image shown is a snapshot of the debris flow simulation.

[0051] Figure 16 The figure shows the evolution of average simulation speed and speedup ratio for debris flow on different GPUs. Detailed Implementation

[0052] To facilitate understanding of the present invention, a more complete description will be given below with reference to the accompanying drawings. Typical embodiments of the invention are shown in the drawings. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete.

[0053] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention.

[0054] The general idea of ​​this invention is to address the problem that existing GPU acceleration methods for particle numerical simulations only utilize the basic cores on the GPU (such as NVIDIA CUDA cores) and do not fully utilize the hardware resources of ray tracing GPUs (such as ray tracing cores specifically used for image rendering). This invention proposes a method for accelerating particle numerical simulations by using ray tracing cores, thereby fully utilizing the computing resources of the GPU, improving the computational efficiency of particle numerical simulations, and promoting the field of high-performance computing in related industrial software.

[0055] To better understand the above technical solutions, the following will describe the above technical solutions in detail with reference to the accompanying drawings and specific embodiments. It should be understood that the embodiments of the present invention and the specific features in the embodiments are detailed descriptions of the technical solutions of this application, rather than limitations on the technical solutions of this application. In the absence of conflict, the embodiments of the present invention and the technical features in the embodiments can be combined with each other.

[0056] This invention transforms the neighbor search problem in particle numerical simulation into a ray tracing problem, thereby leveraging the powerful capabilities of the ray tracing core.

[0057] (I) Provide a unified description of the neighborhood search problem in particle numerical methods.

[0058] Discrete element method (DEM) and molecular dynamics are discrete numerical methods, while perifield dynamics and smoothed particle hydrodynamics are continuous medium methods. In DEM, the overlap between two particles with a given size (radius) is used to determine the interaction forces between the particles, while molecular dynamics directly applies a cutoff distance to solve for the intermolecular forces. For perifield dynamics and smoothed particle hydrodynamics, the influence range is used to determine the interaction weights between a given particle and its neighboring particles. Among these particle-based methods, such as... Figure 2 As shown, particle j is considered a potential neighbor of particle i (particle A in the figure) when the following conditions are met.

[0059] d<(1+δ)l (1)

[0060] Where d is the distance between particles i and j; l is the generalized search radius, which for the discrete element method is equal to the sum of the radii of the two particles, for molecular dynamics it is the cutoff distance, and for near-field dynamics and smooth particle hydrodynamics it is equal to the influence range; δ is called the relaxation size factor, used to scale up the search radius l proportionally, thereby reducing the neighborhood search frequency. A neighborhood search is invoked once the cumulative displacement of any particle (since the last neighborhood search) exceeds its relaxation size (i.e., δl). Note: For non-spherical particles in the discrete element method, l can be the sum of the boundary sphere radii of the two particles under consideration.

[0061] (ii) Abstracting the ray tracing problem into a general collision detection problem

[0062] Ray tracing is a technique used to simulate the reflection and refraction of light in the real world. It is widely used in computer graphics and games to render digital images with an added sense of realism. Figure 3(a) illustrates the main idea of ​​ray tracing. Each pixel in the rendered image corresponds to a ray originating from the camera. The final color of that pixel is calculated based on the colors of the objects struck by the corresponding ray. During the propagation path, a ray may strike multiple objects sequentially, or it may strike none, such as... Figure 3 As shown in (b). Typically, tens of thousands of ray hits (collision detection) may occur in a real-world scene. The ray tracing problem is essentially the collision detection problem between rays and scene objects. To facilitate ray path identification and efficient collision detection, hierarchical bounding boxes have been widely used in ray-traced image rendering.

[0063] (III) Hardware acceleration using NVIDIA GPU's ray tracing cores and OptiX library.

[0064] To accelerate real-time ray tracing in games, NVIDIA has introduced a new type of hardware called Ray Tracing Cores in its RTX GPUs since its initial architecture. In an RTX GPU, each stream processor has one Ray Tracing Core, and each Ray Tracing Core has two physics units: a bounding box intersection solver and a ray-triangle intersection solver. As the names suggest, the bounding box intersection solver is designed for traversing hierarchical bounding boxes, while the ray-triangle intersection solver is used to solve for ray collisions with triangles. Accompanying the Ray Tracing Cores is the newly designed OptiX 7 library, a ray tracing application programming interface developed by NVIDIA. Figure 4 (b) illustrates the general ray tracing workflow in OptiX 7, with yellow boxes indicating custom procedures. After starting the ray tracing pipeline, a ray with a custom position, length, and direction is generated and then passed to the accelerator structure for traversal, which iterates through the hierarchical bounding boxes of all objects in the scene. This accelerator traversal runs entirely on the bounding box intersection evaluator of the ray tracing core. Once a ray intersects a bounding box, a custom intersection procedure is invoked, followed by an arbitrary hit procedure. Once the ray hits the nearest geometry, the nearest hit procedure is invoked. When the ray misses a hit, a miss procedure is invoked. Notably, OptiX 7 implements state-of-the-art bounding box construction and traversal algorithms at the driver level, significantly improving ray tracing performance and simplifying ray tracing programming.

[0065] Based on this, the present invention provides a method for accelerating particle numerical simulations based on ray tracing, which uses the ray tracing core of a GPU for hardware acceleration, and includes the following steps:

[0066] S1. Expand the search range of the particle to be searched to a neighborhood search region bounding box, wherein the neighborhood search region bounding box includes multiple surrounding particles located around the particle to be searched;

[0067] Specifically, such as Figure 5 As shown in (a), firstly, the bounding box (small box) of the particle A to be searched or its influence region is expanded into the bounding box (large box) of the neighborhood search region, which contains multiple surrounding particles. Therefore, determining whether a surrounding particle is a potential neighbor of particle A is equivalent to testing whether the center Op of the surrounding particles is within the expanded bounding box, see... Figure 5 (b)

[0068] S2. Call the GPU's ray tracing core to emit rays from the centers of the multiple surrounding particles, and detect whether the rays intersect with the bounding box of the neighborhood search region;

[0069] Specifically, by expanding the search area, the problem can be transformed into a ray tracing problem: emitting an infinitesimal ray from the center of a given surrounding particle, and detecting whether this ray intersects with an expanded bounding box (large frame). See [link to relevant documentation]. Figure 5 (c). Given that the ray length is infinitesimal, the ray direction can be arbitrary. After the ray generation procedure, the intersection procedure is invoked once a given ray i intersects with the bounding box during hierarchical bounding box traversal.

[0070] S3. When the ray intersects with the bounding box of the neighborhood search region, the intersection procedure is called to update the neighbor list.

[0071] Specifically, Algorithm 1 demonstrates the pseudocode for the neighborhood search bounding box intersection procedure. Generally, once the ray tracing program reaches the bounding box intersection procedure, surrounding particles are contained within the bounding box of particle i corresponding to ray number i, such as... Figure 5 As shown in (d). However, another query is performed in line 6, referring to equation (1), for better comparison with traditional cell-based algorithms. Starting from line 7, the neighbor list NL i The construction method is the same as the cell-based algorithm, except that N1 and N2 are updated in the private memory of ray i.

[0072] Algorithm 1: Neighborhood Search Bounding Box Intersection Algorithm

[0073]

[0074]

[0075] Therefore, step S3 includes:

[0076] S31. Calculate the distance between the center position of the surrounding particles corresponding to the intersecting rays and the center position of the particle to be searched;

[0077] S32. Determine whether the distance is less than a preset threshold, wherein the preset threshold is calculated using Formula 1 above;

[0078] S33. When the distance is less than a preset threshold, compare the ray number i of the intersecting rays with the particle number j of the particle to be searched;

[0079] S34. If particle number j is greater than ray number i, then update the neighbor list NL[i*offset+N1] = j; if particle number j is less than ray number i, then update the neighbor list NL[(i+1)*offset-N2-1] = i, where offset is the maximum number of neighbors for each particle, and N1 and N2 are the number of neighbors whose particle numbers are greater than and less than the ray number, respectively. Optionally, in another embodiment of the present invention, if particle number j is less than ray number i, then update the neighbor list NL[i*offset+N1] = j; if particle number j is greater than ray number i, then update the neighbor list NL[(i+1)*offset-N2-1] = i, where offset is the maximum number of neighbors for each particle, and N1 and N2 are the number of neighbors whose particle numbers are less than and greater than the ray number, respectively.

[0080] Furthermore, in OptiX, custom programs execute sequentially for each ray assigned to a single thread, while simultaneously generating and processing large numbers of rays in parallel. To maximize performance, adjacent rays should follow similar paths when traversing the hierarchical bounding box. An efficient approach is to reorder the particles according to their space-fill curves, such as Morton's curves and Hilbert's space-fill curves. For ease of implementation and without loss of generality, Morton's curves are used by default in this embodiment unless otherwise stated. Additionally, the BVH can be constructed via the OptiX API, but is not shown here for brevity. However, special attention should be paid to the degradation of BVH quality, which leads to worse BVH traversal, especially for dynamic simulations. Periodically rebuilding the BVH is an effective remedy to avoid this problem.

[0081] Furthermore, when using a ray tracing core for particle neighborhood search, existing methods are employed, such as allocating a GPU thread for each particle or particle pair's interaction to accelerate computation. Based on this, the aforementioned computation process can be parallelized with particle kinematics solving, thereby making better use of computational resources.

[0082] Accordingly, the present invention also provides a particle-based numerical simulation acceleration system based on ray tracing, which uses the GPU's ray tracing core and OptiX library for hardware acceleration, including:

[0083] An extension module is used to extend the search range of the particle to be searched to a neighborhood search region bounding box, wherein the neighborhood search region bounding box includes multiple surrounding particles located around the particle to be searched;

[0084] The detection module is used to invoke the GPU's ray tracing core to emit rays from the centers of the multiple surrounding particles, and to detect whether the rays intersect with the bounding box of the neighboring search region; and

[0085] An update module is used to call an intersection procedure to update the neighbor list when the ray intersects with the bounding box of the neighborhood search region.

[0086] Furthermore, the update module includes:

[0087] The prediction unit is used to calculate the distance between the center position of the surrounding particles corresponding to the intersecting rays and the center position of the particle to be searched, and to determine whether the distance is less than a preset threshold.

[0088] A comparison unit is used to compare the ray number i of the intersecting rays with the particle number j of the particle to be searched;

[0089] An update unit is used to update the neighbor list NL[i*offset+N1]=j when the particle number j is greater than the ray number i, and to update the neighbor list NL[(i+1)*offset-N2-1]=i when the particle number j is less than the ray number i, where offset is the maximum number of neighbors for each particle, and N1 and N2 are the number of neighbors whose particle numbers are greater than and less than the ray number, respectively. Optionally, the update unit can also update the neighbor list NL[i*offset+N1]=j when the particle number j is less than the ray number i, and to update the neighbor list NL[(i+1)*offset-N2-1]=i when the particle number j is greater than the ray number i, where offset is the maximum number of neighbors for each particle, and N1 and N2 are the number of neighbors whose particle numbers are less than and greater than the ray number, respectively.

[0090] This invention proposes a method for accelerating ray tracing using particle-based numerical methods (discrete element method, molecular dynamics method, peri-field dynamics method, and smoothed particle hydrodynamics method). The hardware acceleration of ray tracing using the discrete element method will be used as an example for illustration.

[0091] (I) Discrete Element GPU Acceleration Algorithm

[0092] Algorithm 2 summarizes the computation flow of a general-purpose parallel discrete element method on a GPU. For each discrete element iteration computation step, three main procedures are executed sequentially by launching the corresponding CUDA kernel: (a) updating the neighbor list; (b) calculating the contact forces of all contacts; and (c) calculating the motion of all particles. Note: In the case of spherical particles, updating the neighbor list is more time-consuming than the other two processes, but the execution frequency can be reduced by introducing a bounding box magnification factor in Equation (1) and caching potential neighbors in the neighbor list. To this end, once any particle moves a distance greater than a given threshold relative to the reference position, the switch variable updateNL is marked as true, thus triggering the neighbor list update in the next computation step.

[0093] Algorithm 2 Discrete Element GPU Accelerated Algorithm Calculation Flow

[0094]

[0095] (ii) Traditional cell neighborhood search algorithms that only use ordinary computing cores

[0096] Traditional cell-based neighborhood search algorithms mainly consist of two steps: 1) establishing a mapping relationship between particles and cell numbers; 2) neighborhood search. The algorithm is described in detail below:

[0097] (1) Mapping between particles and cells. Given a set of particle positions X Np (Np is the number of particles), iterate through all particles and calculate the index PC of the cell to which particle i belongs. i Similarly, create a cell-to-particle list (CP). Ns (Ns is the cell number) to map cells to particles. Assume a particle is mapped to only one cell, and the grid cell is greater than the maximum neighbor search distance (1+δ)l defined in equation (1). Then, perform a cardinal sort on the particle-cell pairs with respect to the cell number so that all particles given cell i can be accessed via the lookup list CP. i Indexing is performed. Table 1 lists the indexes based on... Figure 1 PC obtained from the configuration in (b) Np With CP Ns The relationship.

[0098] Table 1: List of Particle and Cellular Numbers (PC) Np With CP Ns Relationship diagram

[0099]

[0100] (2) Neighborhood search: The neighbors of a given particle i can be searched through the following three steps:

[0101] i) Locate the cell to which the particle belongs, i.e., the cell number is PC.i ;

[0102] ii) Find the 26 nearest neighbor cells (3D), where the internal particles are potential neighbors of a given particle i;

[0103] iii) Iterate through all particles belonging to the specified cell and its neighboring cells, and perform neighbor detection. Then create a neighbor list NL. i The sublist is then filled with the neighboring particle number j. For ease of calculation, if particle number j is greater than particle number i, then j is pushed into the sublist from the left; otherwise, it is pushed into the sublist from the right.

[0104] Note: In practice, different particles may have different numbers of neighbors, thus requiring different amounts of memory to store the neighbor list. To avoid the overhead caused by repeated memory allocation, the neighbor list NL can be pre-allocated. i It has a sufficiently large memory segment to cover all possible neighbors. Furthermore, in most particle-based numerical methods, particle pair computations only need to be performed once for each pair of particles. Since each pair of particles is each other's neighbor, we only need to check neighbors whose particle numbers are less than a given particle number; and because two additional lists (i.e., sublist N) are used... jgi Np List of sub-lists N jli Np These are used to record the number of neighbors with particle number j that is greater than or less than particle number i, so that the target neighbors of a given particle i can be easily accessed.

[0105] (III) Test System Hardware and Software Environment

[0106] Software environment: Linux / Ubuntu 18.04 operating system, CUDA 11.2, OptiX 7.2, GCC 7.5

[0107] The graphics card hardware parameters are shown in the table below:

[0108]

[0109] (IV) Basic Ray Tracing Configuration

[0110] The entry function for ray tracing using the OptiX library is configured as follows:

[0111]

[0112] For each particle, the above ray tracing call will emit a beam of length 1x10. -16The ray is used to identify particles and rays. The ray number is initialized to the particle number. Furthermore, recently hit and any hit are disabled to improve performance by setting the OptiX macros `OPTIX_RAY_FLAG_DISABLE_CLOSESTHIT` and `OPTIX_RAY_FLAG_DISABLE_ANYHIT`. In OptiX, each ray can share data between different programs by setting a private memory variable, which is not allowed in CUDA internal functions (e.g., shared memory). Similar to the cell-based neighborhood search method, two auxiliary lists N... jgi Np and N jli Np It is introduced and passed to the ray tracing call entry function as a private memory variable for particle i (or ray i).

[0113] (V) Test Examples

[0114] The test cases include three typical problems at different scales: particle accumulation, collapsing columns, and debris flow. All material property settings for the simulations are listed in Table 2, with gravitational acceleration and time step set to 10 m / s². 2 And 0.001s. The simulation is divided into three different tests, and the relevant settings parameters are shown in Table 3.

[0115] Table 2: Material parameter values ​​in the simulation

[0116] parameter Value <![CDATA[Particle mass density [kg / m 3 > 2000 <![CDATA[Particle contact stiffness k n , k t [N / m]]]> <![CDATA[1x10 5 ]]> <![CDATA[Wall contact stiffness k n , k t [N / m]]]> <![CDATA[1x10 7 ]]> interparticle friction coefficient 0.5 Wall friction coefficient 1.0 Viscous damping coefficient 0.5 mass damping coefficient 0.1

[0117] Table 3: Test Setup for the Three Groups

[0118]

[0119] 1. Test 1: Simulation of collapsing columns with different bounding box magnification factors

[0120] A total of 800,000 spheres with a radius of 0.1m are placed on a 71x80x141 grid, with a grid spacing of 2.4 times the radius of the spheres. Figure 6 As shown in (a). To prevent particles from stacking and forming a lattice configuration, all particles are randomly offset from the grid nodes, specifically in the x and y directions within a range of 0 to 0.1 times the particle radius. The particles are then allowed to settle under gravity until a stable accumulation is formed within 40 seconds (i.e., 40,000 steps). The final state is shown in [the diagram]. Figure 6 (b) It can be seen that the particle column collapses into a relatively straight slope with a curved tail, which is consistent with existing numerical and experimental observations. At the same time, the angle of repose of 27.4° is also consistent with the report in the literature, indicating that the developed discrete element GPU acceleration program is correct.

[0121] Five simulations were performed on a Quadro RTX 5000 GPU with bounding box magnification factors s of 0, 0.1, 0.2, 0.3, and 0.4, respectively. Larger factors indicate a lower frequency of neighborhood searches in the simulation. Note: In the extreme case where s is zero, the neighbor list needs to be updated at every time step, which is typically time-consuming and generally not used in practical simulations. To measure the simulation performance at a given time, the average computation speed was obtained by dividing the cumulative time step by the current total time elapsed. Therefore, dynamic speed can be recorded during the simulation. Figure 7 As shown in (a), the simulation speed drops rapidly at the start of the simulation, attributed to the rapid increase in coordination number (the number of contacts per particle) during column collapse. Then, as the stacking approaches static, the simulation speed converges to a relatively stable value for each simulation setting. Furthermore, both the initial and stable simulation speeds increase with increasing s. This is because a larger s helps cache more neighbors, some of which may not be in contact with a given particle, thus reducing the frequency of neighbor list updates. However, when s exceeds 0.2, the overall speed cannot be improved further. Therefore, a bounding box scaling factor s of 0.2 was chosen for the remaining simulations in this embodiment. Note: Increasing s may have the side effect of potentially increased memory allocation; we set a fixed number of neighbors of 32 for each particle. That is, the neighbor list can only cache a maximum of 32 neighbors per particle. When more neighbors need to be cached, the memory for the neighbor list must be increased.

[0122] The above analysis applies to both traditional cell-based neighborhood search algorithms and the ray tracing algorithm proposed in this invention. Furthermore, regardless of the value of s, simulations using the inventive method are faster than traditional simulations, indicating that the inventive method achieves higher computational efficiency. Quantitatively, the speed ratio of traditional and conventional methods is plotted on... Figure 7 In (b), it can be seen that, except for the case where s=0, the speedup ratio first increases rapidly and then gradually decreases. For the case where s=0, the speedup ratio increases monotonically, reaching a maximum of 1.9, which means that when we have to update the neighbor list at every time step, the inventive method can be 1.9 times faster. This means that we can achieve better performance while maximizing the use of ray tracing cores. In fact, as s increases, the use of ray tracing cores may decrease because the update frequency of the neighbor list decreases. This also shows that the inventive method is suitable for dynamic simulations, especially in the first few hundred time steps when the speedup ratio increases rapidly, at which point the particles rearrange themselves violently due to collapse.

[0123] In summary, in the current tests on the Quadra RTX 5000 GPU, the stable simulation calculation speed of 800,000 particles can reach 450 steps per second. Using the method of this invention can achieve a speed improvement of at least 1.2 to 1.4.

[0124] 2. Test 2: Simulation of Collapsed Columns with Different Number of Particles

[0125] Similar to the setup in Test 1, but this test will simulate a collapsing column using 100,000 to 1.6 million particles. For better comparison, only the y-direction depth of the particle column will be adjusted to obtain different particle counts, such as... Figure 8 As shown.

[0126] Figure 9 The results show the changes in average computation speed and speedup ratio of the conventional and inventive methods under different particle numbers during the collapsing column simulation. As the particle number increases from 100,000 to 1.6 million, the initial simulation speed drops significantly from 5,000 steps per second to 450 steps per second, while the stable simulation speed range decreases from 3,000 steps per second to 220 steps per second. In all cases, the overall performance of the inventive method is superior to the conventional method. Furthermore, with the increase in particle number, the ray tracing core can be utilized more extensively, and the speedup ratio also increases. Particularly for the 1.6 million particle case, the peak and final speedup ratios reach 1.6 and 1.3, respectively.

[0127] To compare the scalability of traditional and invented methods, a computational speed reduction weight is introduced, defined as the ratio of the simulation speed with 2N particles to that with N particles. As the name suggests, the computational speed reduction weight shows the degree to which the overall simulation speed is reduced when the number of particles doubles. A computational speed reduction weight of 0.5 indicates that the system can scale linearly without degrading computational performance. A computational speed reduction weight greater than 0.5 indicates that the hardware is not being utilized to its full potential, while a weight less than 0.5 indicates hardware deficiencies. Figure 10 As shown in (a), as expected, larger particle systems exhibit a smaller computational speed reduction ratio. In particular, in all cases, simulations based on the inventive method show a larger computational speed reduction ratio than conventional methods, indicating better scalability of the inventive method. Note: Here, hardware refers to all hardware resources in the computing system, such as CUDA cores, ray tracing cores, GPU memory size, and bus width. A detailed investigation of these resources is complex and beyond the scope of this embodiment. Figure 10 (b) only shows the memory allocation, where it is observed that memory allocation increases almost linearly with the number of chips. Furthermore, the simulation based on the inventive method uses slightly more memory than the conventional method.

[0128] 3. Test 3: Simulation under different GPU configurations

[0129] In addition to the collapsed column simulations in Tests 1 and 2, Test 3 also simulated two other problems: particle accumulation and debris flow. Building on the findings and discussions in Test 2, all simulations in Test 3 will use 1.6 million particles.

[0130] 3.1 Simulation of Particle Accumulation and Collapsed Columns

[0131] The setup for particle accumulation simulation is almost identical to that for collapsed column simulation, except that a sidewall is added to the right surface of the collapsed column. Figure 11 It shows a series of snapshots of 1.6 million particles during the deposition process.

[0132] Figure 12 The evolution of simulation computation speed for particle stacking and collapse pillar tests using 1.6 million particles on different GPUs is shown. It can be seen that the particle stacking and collapse pillar simulations initially had almost the same computation speed (before approximately 2000 time steps). This is because most particles settle during this period for both cases. Furthermore, for all GPUs, simulations based on the inventive method are faster than those based on the conventional method for the stacking and collapse problem. Specifically, the simulation computation speed is positively correlated with the overall performance of the GPU, such as... Figure 13 (a). Interestingly, the accumulation simulation is faster than the collapse column simulation, while the collapse column simulation has a larger speedup ratio than the accumulation simulation, such as Figure 13 As shown in (b). The main reason is that the collapse column simulation is more dynamic due to the large-scale rearrangement of particles, which means that the neighbor list is updated more frequently. This also verifies that maximizing the utilization of ray tracing cores can achieve higher computational efficiency, as described in Test 1. Furthermore, it is worth noting that the performance gain (in terms of speedup ratio) based on the inventive method depends on the proportion of ray tracing cores in the overall performance of the GPU, not just the number of ray tracing cores. For example, the RTX 2080 Ti has the most ray tracing cores, but its speedup ratio is the lowest.

[0133] 3.2 Debris Flow Simulation

[0134] Debris flow simulations are more complex because they involve large-scale moving particles and static, highly irregular terrain boundaries, such as... Figure 14 As shown. In this embodiment, the terrain is modeled using a triangular mesh created based on a digital elevation model (DEM). Debris particles are generated from the grid mesh of the sliding region, i.e., the terrain before sliding minus the terrain after sliding. The grid spacing is 0.25m, corresponding to approximately 1.6 million particles with a radius of 0.125m. Note: Due to the grid regularity of the DEM, adjacent triangular faces can be effectively identified using the position of a given particle; therefore, neighboring triangular faces of a particle are not cached to improve efficiency. If a particle contacts multiple triangles, the contact force is calculated by weighting the particle-triangle contact area for accuracy.

[0135] The simulation time step for debris flow is 100,000 steps, with a fixed time step of 0.001 seconds. Figure 15Several snapshots of the debris flow at 0s, 20s, 40s, and 100s are shown. The simulation captures the general flow evolution of debris flow well. The simulation computation speedup and speedup ratio of the debris flow simulation are recorded in [reference needed]. Figure 16 Similar to the collapsed column scenario, the simulation speed drops sharply in the first few hundred time steps before stabilizing. Throughout the simulation, the speedup varies between 1.1 and 1.6, depending on the specific GPU specifications. In terms of memory consumption, the simulation based on the invented method requires approximately 3138MB, slightly less than the simulation based on traditional cells (3418MB). (Comparison) Figure 10 The memory consumption results shown indicate that the simulations based on the inventive method and the traditional method consume almost the same amount of memory, although there may be some differences depending on the specific problem. Overall, the inventive method exhibits better computational performance (in terms of computational efficiency and memory consumption) than the traditional cell-based method.

[0136] (VI) Test Conclusion

[0137] The simulation efficiency based on the method of this invention is 1.1 to 1.6 times that of traditional methods, and the required GPU memory is almost the same as that of traditional cell-based simulations. For large particle systems and dynamic problems, the method of this invention can make maximum use of ray tracing cores and GPU hardware resources, resulting in a more significant performance improvement in simulation calculations.

[0138] It should be noted that the above embodiments are illustrative of the invention and not restrictive, and that those skilled in the art can devise alternative embodiments without departing from the scope of the appended claims. In the claims, any reference signs placed between parentheses should not be construed as limiting the claims. The word "comprising" does not exclude the presence of elements or steps not listed in the claims. The word "a" or "an" preceding an element does not exclude the presence of a plurality of such elements. The invention can be implemented by means of hardware comprising several different elements and by means of a suitably programmed computer. In the unit claims enumerating several means, several of these means may be embodied by the same item of hardware. The use of the words first, second, and third, etc., does not indicate any order. These words can be interpreted as names.

Claims

1. A light-ray tracing based particle-like numerical simulation acceleration method, characterized in that, Hardware acceleration based on a GPU-based ray tracing core, comprising the following steps: extending a search range of a to-be-searched particle to a neighborhood search region bounding box, the neighborhood search region bounding box comprising a plurality of surrounding particles located around the to-be-searched particle; calling a ray tracing core of a GPU to respectively emit rays from centers of the plurality of surrounding particles, and detecting whether the rays intersect the neighborhood search region bounding box; and when the rays intersect the neighborhood search region bounding box, calling an intersection program to update a neighbor list; the step of calling the intersection program to update the neighbor list comprises: comparing a ray number i of the intersected ray with a particle number j of the to-be-searched particle; if the particle number j is greater than the ray number i, updating the neighbor list NL[i*offset + N1] = j, and if the particle number j is less than the ray number i, updating the neighbor list NL[(i+1)*offset - N2-1] = i, wherein offset is a maximum number of neighbors of each particle, N1 and N2 are numbers of neighbors whose particle numbers are greater than and less than the ray number, respectively; if the particle number j is less than the ray number i, updating the neighbor list NL[i*offset + N1] = j, and if the particle number j is greater than the ray number i, updating the neighbor list NL[(i+1)*offset - N2-1] = i, wherein offset is a maximum number of neighbors of each particle, N1 and N2 are numbers of neighbors whose particle numbers are less than and greater than the ray number, respectively.

2. The light-ray tracing based particle-like numerical simulation acceleration method of claim 1, wherein, before the step of comparing the ray number i of the intersected ray with the particle number j of the to-be-searched particle, further comprising: calculating a distance between a center position of a corresponding surrounding particle of the intersected ray and a center position of the to-be-searched particle; and judging whether the distance is less than a preset threshold value, the preset threshold value being (1+δ)l, wherein δ is a relaxation size factor, and l is a generalized search radius.

3. A light ray tracing based particle-like numerical simulation acceleration system, characterized in that, Hardware acceleration based on a GPU-based ray tracing core, comprising: an extension module configured to extend a search range of a to-be-searched particle to a neighborhood search region bounding box, the neighborhood search region bounding box comprising a plurality of surrounding particles located around the to-be-searched particle; a detection module configured to call a ray tracing core of a GPU to respectively emit rays from centers of the plurality of surrounding particles, and detect whether the rays intersect the neighborhood search region bounding box; and an update module configured to, when the rays intersect the neighborhood search region bounding box, call an intersection program to update a neighbor list; the update module comprises: a comparison unit configured to compare a ray number i of the intersected ray with a particle number j of the to-be-searched particle; The updating unit is configured to update the neighbor list NL[i*offset+N1]=j when the particle number j is greater than the ray number i, and update the neighbor list NL[(i+1)*offset-N2-1]=i when the particle number j is less than the ray number i, where offset is a maximum number of neighbors of each particle, N1 and N2 are numbers of neighbor particles whose numbers are greater and less than the ray number, respectively. The updating unit is configured to update the neighbor list NL[i*offset+N1]=j when the particle number j is greater than the ray number i, and update the neighbor list NL[(i+1)*offset-N2-1]=i when the particle number j is less than the ray number i, where offset is a maximum number of neighbors of each particle, N1 and N2 are numbers of neighbor particles whose numbers are greater and less than the ray number, respectively.

4. The light-ray tracing based particle-like numerical simulation acceleration system of claim 3, wherein, The updating module further includes: The pre-judging unit is configured to calculate a distance between a center position of a corresponding surrounding particle intersecting the ray and a center position of the particle to be searched, and determine whether the distance is less than a preset threshold, the preset threshold being (1+δ)l, where δ is a relaxation size factor, and l is a generalized search radius.

5. A computer readable storage medium, characterized in that, A computer readable storage medium stores a computer program, and the computer program, when executed by a processor, implements the steps of the particle-based numerical simulation acceleration method based on ray tracing according to any one of claims 1 to 2.

6. An optical ray tracing based particle-like numerical simulation acceleration apparatus, characterized by, A computer readable storage medium stores a computer program, and the computer program, when executed by a processor, implements the steps of the particle-based numerical simulation acceleration method based on ray tracing according to any one of claims 1 to 2.

Images