A heterogeneous parallel solution method for SPH based on adaptive time-series dynamic load balancing

Abstract

This invention discloses a heterogeneous parallel solution method for SPH based on adaptive time-series dynamic load balancing, belonging to the field of computer simulation technology. The method improves the resource utilization and computational efficiency of multi-architecture parallel computing systems by introducing a dynamic load balancing mechanism based on computation time. By monitoring the actual computation time of each computing node (CPU, GPU) in real time, rather than the number of particles in each process, the method dynamically adjusts the boundary position of sub-regions, ensuring that computing units with different hardware performance can approach their optimal load state. It supports hybrid deployments of multi-core CPUs, multi-GPUs, and hybrid architectures, and the dynamic boundary adjustment is algorithmically independent of the particle attribute classification mechanism. It avoids computational bottlenecks caused by the performance lag of a single hardware component, effectively preventing the waste of overall computing power due to single-node performance bottlenecks.

Description

Technical Field

[0001] This invention belongs to the field of computer simulation technology, and specifically relates to a heterogeneous parallel solution method for SPH based on adaptive time-series dynamic load balancing. Background Technology

[0002] In the field of computational fluid dynamics (CFD) for continuum media, numerical simulation methods are generally classified into two main categories based on the point-line-surface topological configuration characteristics after the computational domain is discretized: meshed methods and meshless methods. Traditional meshed methods, such as the finite difference method (FDM), the finite element method (FEM), and the finite volume method (FVM), have been widely used in natural sciences and engineering practice, effectively handling problems related to conventional fluid flow and fluid-structure interaction. However, the node-connected topological configuration of meshed methods has inherent limitations, leading to significant bottlenecks in scenarios involving complex geometric boundaries, large deformations, multi-scale flow fields, and free surface evolution. In contrast, meshless methods, with their lack of fixed topological constraints, exhibit unique advantages in these complex scenarios.

[0003] Smoothed particle hydrodynamics (SPH), as one of the most widely used meshless methods, exhibits excellent adaptability in scenarios with complex geometric boundaries and large material deformations by discretizing the computational domain using particles. However, the computational complexity of SPH methods is generally higher than that of traditional mesh methods. As the number of simulated particles increases, the time consumption characteristics of the single-process computation mode become difficult to meet engineering requirements. Therefore, developing efficient parallel solution techniques has become crucial. Currently, SPH parallel solution schemes are mainly divided into two categories: CPU-based parallel architectures and GPU-based heterogeneous acceleration architectures.

[0004] Among them, CPU-based parallel methods include two modes: shared memory and distributed memory. Shared memory uses thread-level parallelism and achieves data interaction through shared memory. It has the advantages of low communication latency and low programming complexity, but its scalability is limited and it is difficult to support ultra-large-scale computing.

[0005] Another type is the GPU-based heterogeneous acceleration architecture. Leveraging its native advantage of massively parallel computing, GPUs exhibit computational performance far exceeding that of CPUs in solving particle mechanics governing equations. The current mainstream technology adopts a CPU-GPU heterogeneous collaboration model, offloading the serial control logic in the program to the CPU, while offloading large-scale particle computing tasks to the GPU for execution, achieving collaboration between the two through an efficient heterogeneous communication mechanism.

[0006] In distributed parallel computing systems, load balancing is a key technical challenge affecting computational efficiency. Uneven distribution of tasks across computing nodes will significantly reduce parallel computing efficiency. In particular, in particle dynamics simulations, the spatial topology of particles changes dynamically over time, necessitating the introduction of a dynamic load balancing mechanism—that is, dynamically redistributing computational tasks among different processors through real-time monitoring and strategy adjustments.

[0007] In existing technologies, dynamic load balancing schemes are mostly based on an approximate equal distribution strategy of particle numbers. By distributing the total number of particles equally (or approximately equally) according to the number of processors, the computational load of each computing node is balanced. However, this type of method has an inherent limitation: it implicitly assumes that all computing nodes have the same processing performance. In real-world scenarios, when there are differences in the hardware performance of computing nodes, even if the number of particles is evenly distributed, the actual computation time of each node will still vary due to performance differences. This phenomenon of "equal quantity but uneven time" not only leads to a decrease in overall parallel efficiency but also causes the computing resources of high-performance hardware nodes to be idle, resulting in resource waste.

[0008] The information disclosed in this background section is intended only to enhance the understanding of the overall background of the invention and should not be construed as an admission or in any way implying that the information constitutes prior art known to those skilled in the art. Summary of the Invention

[0009] The purpose of this invention is to provide a heterogeneous parallel solution method for SPH based on adaptive time-series dynamic load balancing, which solves the problem of low parallel efficiency of existing dynamic load balancing technology when processor performance is inconsistent, and can ensure that the computation time of each processor process is the same, thereby improving parallel efficiency.

[0010] To achieve the above objectives, this invention provides a heterogeneous parallel solution method for SPH based on adaptive time-series dynamic load balancing, comprising the following steps:

[0011] Step 1: Discretize the computational domain of the physics problem into particles to obtain the initial particle set. The computational domain can be a regular or irregular region, and it is enclosed by a minimum cuboid. The side lengths of the cuboid in the X, Y, and Z directions are respectively... , , ;

[0012] Step 2: Based on the particle spacing and smooth radius The size of the background grid is determined, and the cuboid is filled with the background grid. The number of grids in the X, Y, and Z directions of the cuboid are respectively... , , ;

[0013] Step 3: Determine the sub-region division direction. Based on the number of grids in the X, Y, and Z directions of the cuboid, determine... , , The direction of the maximum value among the three factors is chosen as the direction for dividing the sub-region.

[0014] Step 4: Divide the computational domain into P sub-regions along the division direction determined in Step 3. The number of background grids in each sub-region satisfies |N_i - N_{i+1}| ≤ δ, where δ is a preset tolerance, and the sub-regions are numbered as follows: ;

[0015] Step 5: Assign computational tasks to each sub-region, initialize the initial physical information of particles in each sub-region, and then assign the computational tasks of the sub-regions to the corresponding processes. The computational tasks are correspondingly assigned to processes. ;

[0016] Step 6: For any computation process Determine the attributes of particles in the current sub-region and classify the particle attributes into local particles and swapped particles;

[0017] Step 7: Before each calculation begins, traverse all sub-regions, send the exchanged particles from each sub-region to the corresponding calculation processes in its neighboring sub-regions, and record the physical time consumed by each process in determining particle attributes and sending particle information. ;

[0018] Step 8: For any computation process Calculate the local maximum time step for each process based on the CFL condition. Record the physical time consumed by each process in handling the above tasks. ;

[0019] Step 9: Each process calculates the particle velocity, density, and pressure at the next moment based on the discretized fluid control equations, and records the physical time consumed by each process in handling the above tasks. ;

[0020] Step 10: Based on the particle position vector calculated in steps 6-9, determine whether the particle is still within the current sub-region. Send the information of the out-of-bounds particle to the process corresponding to the adjacent sub-region, and record the time taken by each process to process the out-of-bounds particle. ;

[0021] Step 11: Calculate the total computation time of each process, compare the total computation time difference of adjacent processes, and if the difference exceeds the first threshold, it is determined to be an unbalanced load.

[0022] Step 12: For adjacent sub-computation domains with unbalanced loads, move a row of background meshes along the partitioning direction towards the sub-computation domain corresponding to the process with longer computation time.

[0023] Step 13: Repeat steps 8 to 12 until the calculation reaches the stopping condition.

[0024] Optionally, the background mesh is a regular hexahedral mesh with a side length of... , This is the smoothing factor.

[0025] Optionally, 1.0 ≤ ≤2.0, the k value is dynamically adjusted according to the particle distribution density σ.

[0026] Optionally, in step 5, the initial physical information includes the initial time. Information on density, pressure, velocity, and position vector.

[0027] Optionally, in step 6, the particle's properties are divided into local particles and exchanged particles, specifically including:

[0028] Particle exchange: Each sub-region has an interface with its adjacent sub-regions. Particles in the first layer of the grid closest to this interface, i.e., those less than [a certain distance] from the interface, are exchanged. The particles are the exchange particles of this sub-region;

[0029] Local particles: In each sub-region, particles that do not belong to the exchanged particles are local particles.

[0030] Optionally, in step 7, the sub-region and subregions Adjacent, sub-regions The exchanged particles are sent to the process Sub-region The exchanged particles are sent to the process .

[0031] Optionally, in step 8, the local maximum time step of each process is calculated using the following formula. :

[0032] ,

[0033] in This is an adjustable coefficient. It is usually set to 0.15. For particles in this process The speed of sound, For particles in this process smooth length, For particles in this process speed;

[0034] Arbitrary process Calculated Share with other processes and select the smallest one. As the global time step .

[0035] Optionally, in step 9, the fluid control equations include the mass conservation equation, the momentum conservation equation, and the state equation. For any calculation process... Based on the discretized fluid control equations, namely the mass conservation equation, momentum conservation equation, and state equation, specifically, the control equations are as follows:

[0036] mass conservation equation: In the formula For particles in this process and particles Speed ​​difference, For particles in this process and particles The gradient value of the smooth function;

[0037] Momentum conservation equation: superscript in the formula These represent different directions: X, Y, and Z. This is an artificial adhesive property;

[0038] Equations of state: In the formula For particles in this process The reference density is a constant.

[0039] Based on the above three equations, calculate the next moment in this process. The velocity, density, and pressure of all particles.

[0040] Optionally, after the boundary shift described in step 12, the particle distribution within the sub-region satisfies:

[0041] The boundary of the moved sub-region should maintain a buffer distance of at least 2 hours from the boundary of the adjacent sub-region.

[0042] Optionally, the first threshold is set to 10%.

[0043] Compared with existing technologies, the SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to the present invention has the following advantages or beneficial effects:

[0044] This invention significantly improves the resource utilization and computational efficiency of multi-architecture parallel computing systems by introducing a dynamic load balancing mechanism based on computation time. Compared to the traditional static load allocation strategy based on equal distribution of particles, this solution has the following technical advantages:

[0045] By monitoring the actual computation time of each computing node (CPU, GPU) in real time, rather than the number of particles within each process, the boundary positions of sub-regions are dynamically adjusted to ensure that computing units with different hardware performance can approach their optimal load state. This avoids the computational bottleneck of the entire system caused by the performance lag of a single hardware component, and effectively prevents the waste of overall computing power due to the performance bottleneck of a single node.

[0046] Breaking away from the traditional limitations of parallel computing's reliance on homogeneous hardware, it supports hybrid deployments of multi-core CPUs, multi-GPUs, and mixed architectures. The computing systems deploying software programs are more scalable, the software programs have a wider range of applications, and can fully utilize existing computing resources.

[0047] The dynamic boundary adjustment and particle attribute classification mechanism proposed in this scheme are algorithm-independent and have been verified to be compatible with mainstream particle simulation algorithms such as SPH, DEM, and MPS. Compared with traditional load balancing techniques, it is applicable to more fields. Attached Figure Description

[0048] Figure 1 A flowchart illustrating an embodiment of the SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to one of the present invention.

[0049] Figure 2 A schematic diagram of the sub-region division direction according to an embodiment of the present invention;

[0050] Figure 3 A schematic diagram of re-dividing sub-regions when the m-process is overloaded according to an embodiment of the present invention;

[0051] Figure 4 A schematic diagram illustrating the load balancing of each processor after re-dividing the regions according to an embodiment of the present invention.

[0052] Figure 5 A schematic diagram of sub-region division in a three-dimensional fluid mixing case according to an embodiment of the present invention;

[0053] Figure 6 A schematic diagram of dynamic adjustment of a sub-region in a three-dimensional fluid mixing case according to an embodiment of the present invention;

[0054] Figure 7 A schematic diagram of the three-dimensional impeller rotating sub-region division according to an embodiment of the present invention;

[0055] Figure 8A schematic diagram of the dynamic adjustment of the three-dimensional impeller rotating sub-region according to an embodiment of the present invention. Detailed Implementation

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

[0057] The terms "first," "second," "third," "fourth," etc. (if present) in the specification, claims, and accompanying drawings of this invention 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 embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein.

[0058] It should be understood that in the various embodiments of the present invention, the sequence number of each process does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention.

[0059] It should be understood that in this invention, "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion, for example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0060] It should be understood that in this invention, "multiple" refers to two or more. "And / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, "and / or B" can represent: A existing alone, A and B existing simultaneously, and B existing alone. The character " / " generally indicates that the preceding and following related objects are in an "or" relationship. "Contains A, B, and C", "Contains A, B, and C" means that all three A, B, and C are contained; "Contains A, B, or C" means that one of A, B, and C is contained; "Contains A, B, and / or C" means that any one, two, or three of A, B, and C are contained.

[0061] Depending on the context, "if" as used here can be interpreted as "when," "when," "in response to determination," or "in response to detection."

[0062] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0063] First embodiment:

[0064] like Figure 1 As shown, the SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to a preferred embodiment of the present invention includes the following steps:

[0065] Step 1: Discretize the computational domain of the physics problem into particles to obtain the initial particle set. The computational domain can be a regular or irregular region, and it is enclosed by a minimum cuboid with side lengths of 1 / 2 in the X, Y, and Z directions. , , ;

[0066] Step 2: Based on the particle spacing and smooth radius The size of the background grid is determined, and the cuboid is filled with the background grid. The number of grids in the X, Y, and Z directions of the cuboid are respectively... , , That is, the number of grid cells in a certain direction is equal to the corresponding side length of the cuboid divided by the side length of the grid cells, rounded up to the nearest integer. Specifically, the background grid is a regular hexahedral grid with a side length of... , Smoothing factor, 1.0 ≤ ≤2.0, the k value is dynamically adjusted according to the particle distribution density σ. It can be set to 1.5.

[0067] Step 3: Determine the sub-region division direction. Based on the number of grids in the X, Y, and Z directions of the cuboid, determine... , , The size of the three factors is considered, and the direction of the maximum value among them is chosen as the direction for dividing the sub-region. Figure 2 As shown;

[0068] Step 4: Divide the computational domain into P sub-regions along the division direction determined in Step 3. The number of background grids in each of the P sub-regions is approximately equal, and the number of background grids in each sub-region satisfies |N_i - N_{i+1}| ≤ δ, where δ is a preset tolerance. The sub-regions are numbered as follows: ;

[0069] Step 5: Assign computational tasks to each sub-region and initialize the initial physical information of the particles in each sub-region, that is, assign an initial time to the particles in each sub-region. The system collects physical information such as density, pressure, velocity, and position vector, and then assigns the computational tasks for sub-regions to the corresponding processes. Specifically, it divides the sub-regions... The computational tasks are correspondingly assigned to processes. ;

[0070] Step 6: For any computation process The properties of particles in the current sub-region are determined and classified into local particles and exchange particles. It should be noted that exchange particles are particles that are close to the first layer of the grid near the interface, i.e., particles that are less than s away from the interface. Local particles are particles that are not exchange particles in each sub-region.

[0071] Step 7: Before each calculation begins, traverse all sub-regions and send the exchange particles from each sub-region to the calculation processes corresponding to its neighboring sub-regions; specifically, sub-regions and subregions Adjacent, sub-regions The exchanged particles are sent to the process Sub-region The exchanged particles are sent to the process Record the physical time consumed by each process in determining particle properties and sending particle information. .

[0072] Step 8: For any computation process Calculate the local maximum time step for each process based on the CFL condition. Record the physical time consumed by each process in handling the above tasks. ;

[0073] Step 9: Each process calculates the particle velocity, density, and pressure at the next moment based on the discretized fluid control equations, and records the physical time consumed by each process in handling the above tasks. ;

[0074] Step 10: Calculate the total number of particles in this process based on steps 6-9. The position vector at time is used to determine whether the particle is still located in the current sub-region. Within the domain, for particles that do not meet this condition, the information of each out-of-bounds particle is sent to the process corresponding to the adjacent sub-computation domain, and the time taken by each process to process the out-of-bounds particle is recorded. ;

[0075] Step 11: Calculate the total computation time of each process and compare it with the total time difference of adjacent processes. If the difference exceeds a first threshold, it is determined to be a load imbalance. Specifically, the computation time of each process recorded in steps 7, 8, 9, and 10 is used as the basis for the calculation. , , , The sum of these two values ​​represents the total calculation time. By comparing the computation time difference between two adjacent processes, if the computation time difference between the two adjacent processes does not exceed 10% (adjustable), the load of the two processes is determined to be balanced. If the computation time difference between the two adjacent processes exceeds this threshold, the computation load of the two processes is determined to be unbalanced.

[0076] Step 12: For adjacent sub-computation domains with unbalanced loads, move a row of background meshes along the partitioning direction towards the sub-computation domain corresponding to the process with longer computation time, such as... Figure 3 and Figure 4 As shown;

[0077] Step 13: Repeat steps 8 to 12 until the calculation reaches the stopping condition.

[0078] In this embodiment of the invention, in step 5, the initial physical information includes the initial time. Information on density, pressure, velocity, and position vector.

[0079] In this embodiment of the invention, in step 8, the local maximum time step of each process is calculated using the following formula. :

[0080] ,

[0081] in This is an adjustable coefficient. It is usually set to 0.15. For particles in this process The speed of sound, For particles in this process smooth length, For particles in this process Speed; will allow any process Calculated Share with other processes and select the smallest one. As the global time step .

[0082] In this embodiment of the invention, in step 9, the fluid control equations include the mass conservation equation, the momentum conservation equation, and the state equation. For any calculation process... Based on the discretized fluid control equations, namely the mass conservation equation, momentum conservation equation, and state equation, specifically, the control equations are as follows:

[0083] mass conservation equation: In the formula For particles in this process and particles Speed ​​difference, For particles in this process and particles The gradient value of the smooth function;

[0084] Momentum conservation equation: superscript in the formula These represent different directions: X, Y, and Z. This is an artificial adhesive property;

[0085] Equations of state: In the formula For particles in this process The reference density is a constant.

[0086] Based on the above three equations, calculate the next moment in this process. The velocity, density, and pressure of all particles.

[0087] In this embodiment of the invention, after the boundary movement described in step 12, the particle distribution within the sub-region satisfies the following condition: the boundary of the moved sub-region maintains a buffer distance of at least 2h (smooth length) with the boundary of the adjacent sub-region. It should be noted that if the distance between the boundaries of adjacent sub-regions is less than 2h, the particles in the moved sub-region may not be able to fully obtain the physical quantity contributions of neighboring particles due to the truncation of the computational domain, leading to calculation errors. The buffer region can prevent discontinuities in particle distribution caused by abrupt boundary changes and reduce numerical oscillations caused by insufficient interpolation or attribute transfer.

[0088] The embodiments of the present invention can count the total computation time of each process, and then determine whether to dynamically adjust the sub-computation domain boundary based on whether the difference in the total computation time of a single step between two adjacent processes is greater than a preset threshold, thereby achieving dynamic balancing of the computational load of the processes.

[0089] Second embodiment:

[0090] Another embodiment of the SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing of the present invention includes the following implementation steps:

[0091] Step S1, simulation model establishment: First, establish an SPH simulation model for the physical scene to be simulated, including information such as fluid region, boundary conditions, and initial particle state; by introducing a dynamic boundary condition update mechanism, it can adapt to large deformation scenarios.

[0092] Step S2, Computational Domain Discretization: The physical domain is discretized using particle-based methods based on the particle spacing Δx and the smoothing radius h. Preferably, the computational domain is enclosed within a minimum rectangular box, and the side lengths (Lx, Ly, Lz) in the X, Y, and Z directions are determined. The background mesh side length dcell = αh is set, where α is typically 1.5, serving as the minimum unit for subsequent domain partitioning and load balancing adjustments. By employing a minimum rectangular box enclosing strategy, invalid computational regions are reduced, saving memory compared to the traditional bounding box method.

[0093] Step S3, sub-region partitioning direction selection: Compare Lx, Ly, and Lz, select the direction where the maximum value is located as the partitioning axis to maximize the efficiency of inter-process communication; partition along the direction of the maximum size to reduce the frequency of communication between adjacent sub-regions.

[0094] Step S4, Sub-computation Domain Partitioning: Based on the number of solver processes Nproc, the background grid is equally divided into Nproc sub-domains along the main partitioning direction. Each sub-domain has approximately the same number of background grids along the main direction, and each sub-domain is numbered i (i=1,2,...,Nproc). This equal partitioning by Nproc achieves coarse-grained allocation of computational resources, evenly distributes the number of grids along the main direction, and reduces initial load bias.

[0095] Step S5, Process Allocation and Particle Information Initialization: Assign a unique process to each sub-computation domain and initialize the physical information of all particles in the domain, including density ρ, pressure P, velocity v, position r, etc.

[0096] Step S6, Distinguishing between local particles and exchanged particles: In each subdomain, particles in the layer less than h from the neighborhood boundary are defined as exchanged particles, and the rest are local particles. Exchanged particles need to periodically synchronize data with adjacent processes to ensure the integrity of particle forces; by pre-generating a list of exchanged particles when calculating the attributes of local particles, communication delays can be hidden.

[0097] Step S7, Synchronization of exchanged particles: At the beginning of each calculation step, all sub-computation domains are traversed, and the exchanged particle data of this domain is sent to the corresponding process of the adjacent domain through MPI or other parallel communication protocols, and the time T1 of the exchange operation is recorded.

[0098] Step S8, Global Time Step Calculation: Each process calculates the global time step based on the particle properties of its local domain (sound speed c, smoothness length h, maximum velocity). Calculate the local maximum time step The formula used is:

[0099] ,

[0100] in The value is typically 0.15. The global time step Δt is taken for all processes. The minimum value is recorded, and the time taken for the time step calculation is T2.

[0101] Step S9, Parallel solution of the governing equations: Each process independently performs the SPH governing equations (mass conservation, momentum conservation, state equations) calculations on the particles in its domain, updating the velocity, density, pressure, and other properties at the next moment. Record the time T3 spent solving the governing equations.

[0102] Step S10, Particle Out-of-Bounds Detection and Migration: For each updated particle, check if it has exceeded its bounds. If the particle's position exceeds the subdomain of this process, the complete data of that particle needs to be transferred to the corresponding neighboring process. Record the particle migration operation time T4.

[0103] Step S11, Load Balancing Statistics: Calculate the time consumed in the above steps for each process, and accumulate it as the total single-step time Tsum = T1 + T2 + T3 + T4. If the difference in Tsum between adjacent processes is greater than 10% (the threshold is adjustable, usually between 5% and 20%), then the load is considered unbalanced and adjustments are needed.

[0104] Step S12, Subdomain Boundary Adaptive Adjustment: For processes with lighter loads, appropriately increase the subdomain length (in units of one layer of background mesh); for processes with heavier loads, correspondingly reduce the domain range to achieve a re-distribution of particles and computational load.

[0105] Step S13, Iterative loop: Repeat steps S7 to S12 until the simulation termination condition is met (such as the total number of steps, simulation time, etc. reaching a specific threshold).

[0106] Parameter range:

[0107] Particle spacing Δx: generally depends on the physical scene and the required resolution;

[0108] Smooth radius h: usually taken as 1.2~2.0 times Δx, 1.5Δx is recommended;

[0109] Background grid side length dcell: 1.5h recommended;

[0110] Value: Generally 0.1~0.3, 0.15 is recommended;

[0111] Load balancing threshold: can be adjusted according to the actual hardware environment, 10% is recommended;

[0112] Simulation termination condition: can be the total number of steps (e.g., 10,000 steps) or the total simulation time (e.g., 10 seconds).

[0113] This invention's embodiments are based on a dynamic load criterion of total time consumed per process, rather than simply the number of particles. This significantly improves parallel efficiency in heterogeneous hardware computing systems, fully utilizes the computing power of both CPU and GPU, solves the problem of high hardware consistency requirements in traditional load balancing algorithms, and greatly expands the applicability of SPH parallel simulation.

[0114] Third embodiment:

[0115] This embodiment performs parallel simulation of a large-scale three-dimensional fluid mixing process. A high-resolution fluid model with 1 million particles is used, with a spatial discretization accuracy of Δx = 0.02m, a smooth kernel radius of h = 0.03m, and a background mesh size of dcell = 0.045m. The condition number is set to 0.15 to ensure computational stability. To fully utilize heterogeneous computing resources, the simulation is executed in parallel on 8 processes (Nproc=8) consisting of 7 CPU processes and 1 GPU process.

[0116] Initially, the particles are distributed within the container. Based on the aforementioned load balancing scheme, the system automatically divides the computational domain into eight dynamic subdomains corresponding to each processor, such as... Figure 5 As shown. During the simulation, the subdomain boundaries are dynamically adjusted in real time according to changes in particle distribution to ensure a balanced computational load across all regions, such as... Figure 6 As shown. The complete simulation takes 10 seconds (approximately 5000 time steps need to be calculated).

[0117] The results show that after the dynamic load balancing mechanism is effectively implemented, the average single-step computation time difference of each parallel process is stably controlled within 8%. Thanks to this optimized partitioning strategy for heterogeneous hardware and dynamic scenarios, the overall parallel computing efficiency is significantly improved by about 35% compared with the traditional static evenly distributed particle number method. This effectively overcomes the load imbalance problem caused by hardware differences and uneven fluid motion, providing an efficient and reliable computing framework for large-scale multiphysics coupled simulation.

[0118] Example 3

[0119] This embodiment performs efficient simulation of the three-dimensional impeller rotation process. The calculation uses a particle spacing Δx = 0.05m, a smooth kernel radius h = 0.075m, and sets... A value of 0.12 was used to ensure numerical stability. The simulation was conducted on a heterogeneous parallel computing platform, using a total of 4 parallel processes (Nproc=4), namely 3 CPU processes and 1 GPU process, to fully utilize the computational advantages of different processors, such as... Figure 7 As shown.

[0120] The core challenge of the simulation lies in handling the inherent performance differences between heterogeneous hardware (CPU and GPU). To address this, the system monitors the computational load of each process in real time. During the simulation, based on the aforementioned load balancing strategy, the amount of computational tasks allocated to each CPU and GPU process is dynamically adjusted, such as... Figure 8 As shown. This dynamic adjustment mechanism ensures that hardware units with different performance levels can work efficiently, avoiding the bottleneck effect.

[0121] The results show that, thanks to this dynamic load balancing strategy for heterogeneous platforms, computational tasks on different hardware can be completed in a timely manner throughout the simulation, achieving continuous saturation of the computational pipeline and eliminating significant inter-process waiting or synchronization bottlenecks. Compared with fixed task allocation, this method significantly improves overall resource utilization and ensures efficient simulation progress.

[0122] This invention proposes a particle algorithm parallel optimization method based on dynamic balancing of process computation time, breaking through the traditional static load allocation mode based on the number of particles. By monitoring the task completion time of each computing unit (CPU / GPU) in real time, the particle allocation strategy is dynamically adjusted to achieve approximate computation time for hardware with different performance levels. This mechanism effectively avoids the problem of low-performance hardware becoming the system bottleneck in traditional methods and supports improved resource utilization in hybrid architectures.

[0123] By eliminating reliance on homogeneous hardware, the system achieves seamless cross-platform scalability through distributed task queues and adaptive migration algorithms. It supports dynamic addition and removal of computing nodes, and its resource fragmentation rate is significantly better than traditional fixed partitioning schemes. This technology overcomes the limitations of the SPH method and can be migrated to particle simulation scenarios such as molecular dynamics (MD) and discrete element method (DEM). By abstracting 2.3-3.8 times the interaction characteristics of particles, a general load balancing framework is formed, which has been validated in mainstream software such as LAMMPS and GROMACS, resulting in improved speedup.

[0124] The foregoing description of specific exemplary embodiments of the invention is for illustrative and explanatory purposes. These descriptions are not intended to limit the invention to the precise forms disclosed, and it will be apparent that many changes and variations can be made in accordance with the foregoing teachings. The exemplary embodiments were chosen and described in order to explain the specific principles of the invention and its practical application, thereby enabling those skilled in the art to implement and utilize various different exemplary embodiments of the invention, as well as various different choices and variations. The scope of the invention is intended to be defined by the claims and their equivalents.

Claims

1. A heterogeneous parallel solution method for SPH based on adaptive time-series dynamic load balancing, characterized in that, Includes the following steps: Step 1: Discretize the computational domain of the physics problem into particles to obtain the initial particle set. The computational domain can be a regular or irregular region, and it is enclosed by a minimum cuboid. The side lengths of the cuboid in the X, Y, and Z directions are respectively... , , ; Step 2: Based on the particle spacing and smooth radius The size of the background grid is determined, and the cuboid is filled with the background grid. The number of grids in the X, Y, and Z directions of the cuboid are respectively... , , ; Step 3: Determine the number of grid cells in the X, Y, and Z directions of the cuboid. , , The direction of the maximum value among the three factors is chosen as the direction for dividing the sub-region. Step 4: Divide the computational domain into P sub-regions along the division direction determined in Step 3, where the number of background grid cells in each sub-region satisfies |N_ - N_ | ≤ δ, where δ is the preset tolerance, and the sub-region is numbered as follows: ; Step 5: Assign computational tasks to each sub-region, initialize the initial physical information of particles in each sub-region, and then assign the computational tasks of the sub-regions to the corresponding processes. The computational tasks are correspondingly assigned to processes. ; Step 6: For any computation process Determine the attributes of particles in the current sub-region and classify the particle attributes into local particles and swapped particles; Step 7: Before each calculation begins, traverse all sub-regions, send the exchanged particles from each sub-region to the corresponding calculation processes in its neighboring sub-regions, and record the physical time consumed by each process in determining particle attributes and sending particle information. ; Step 8: For any computation process Calculate the local maximum time step for each process based on the CFL condition. Record the physical time consumed by each process in handling the above tasks. ; Step 9: Each process calculates the particle velocity, density, and pressure at the next moment based on the discretized fluid control equations, and records the physical time consumed by each process in handling the above tasks. ; Step 10: Based on the particle position vector calculated in steps 6-9, determine whether the particle is still within the current sub-region. Send the information of the out-of-bounds particle to the process corresponding to the adjacent sub-region, and record the time taken by each process to process the out-of-bounds particle. ; Step 11: Calculate the total computation time of each process, compare the total computation time difference of adjacent processes, and if the difference exceeds the first threshold, it is determined to be an unbalanced load. Step 12: For adjacent sub-computation domains with unbalanced loads, move a row of background meshes along the partitioning direction towards the sub-computation domain corresponding to the process with longer computation time. Step 13: Repeat steps 8 to 12 until the calculation reaches the stopping condition. The background mesh is a regular hexahedral mesh with a side length of [missing information]. , This is the smoothing factor.

2. The SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to claim 1, characterized in that, 1.0≤ ≤2.0, the k value is dynamically adjusted according to the particle distribution density σ.

3. The SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to claim 1, characterized in that, In step 5, the initial physical information includes the initial time. Information on density, pressure, velocity, and position vector.

4. The SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to claim 1, characterized in that, In step 6, the particle properties are divided into local particles and exchanged particles, specifically including: Particle exchange: Each sub-region has an interface with its adjacent sub-regions. Particles in the first layer of the grid closest to this interface, i.e., those less than [a certain distance] from the interface, are exchanged. The particles are the exchange particles of this sub-region; Local particles: In each sub-region, particles that do not belong to the exchanged particles are local particles.

5. The SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to claim 4, characterized in that, In step 7, sub-region and subregions Adjacent, sub-regions The exchanged particles are sent to the process Sub-region The exchanged particles are sent to the process .

6. The SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to claim 1, characterized in that, In step 8, the local maximum time step of each process is calculated using the following formula. : ； in This is an adjustable coefficient. Set to 0.15, For particles in this process The speed of sound, For particles in this process smooth length, For particles in this process speed; Arbitrary process Calculated Share with other processes and select the smallest one. As the global time step .

7. The SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to claim 1, characterized in that, In step 9, the fluid control equations include the mass conservation equation, the momentum conservation equation, and the state equation. For any calculation process... Based on the discretized fluid control equations, namely the mass conservation equation, momentum conservation equation, and state equation, specifically, the control equations are as follows: mass conservation equation: In the formula For particles in this process and particles Speed ​​difference, For particles in this process and particles The gradient value of the smooth function; Momentum conservation equation: superscript in the formula These represent different directions: X, Y, and Z. This is an artificial adhesive property; Equations of state: In the formula For particles in this process The reference density is a constant; Based on the above three equations, calculate the next moment in this process. The velocity, density, and pressure of all particles.

8. The SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to claim 1, characterized in that, After the boundary shift described in step 12, the particle distribution within the sub-region satisfies: The boundary of the moved sub-region should maintain a buffer distance of at least 2 hours from the boundary of the adjacent sub-region.

9. The SPH heterogeneous parallel solution method based on adaptive time-series dynamic load balancing according to claim 1, characterized in that, The first threshold is set to 10%.

Images