A grid Xin'anjiang model parallel method based on compact index and CUDA graph

Abstract

The application discloses a grid Xin'anjiang model parallel method based on compact index and CUDA graph, constructs a downstream one-dimensional index array next and a layered compact index table based on flow direction and calculation level information, realizes compact thread mapping of effective grid of river confluence according to layers, adopts GPU memory pool one-time allocation and offset division of model calculation variable and storage area of river confluence accumulator on the GPU side, combines CPU lock page memory double buffering and GPU double buffering, realizes driving data transmission and calculation overlap in asynchronous transmission flow and calculation flow, simultaneously, through stream capture, instantiates a start sequence of a layered river confluence kernel function into an executable CUDA graph, only updates parameters and starts execution during period iteration, reduces CPU layer-by-layer start overhead, and realizes concurrent consistent accumulation of upstream boundary inflow and upstream and downstream confluence by using double-precision atomic addition.

Description

Technical Field

[0001] This invention relates to parallel computing technology for distributed hydrological models, specifically to a parallel method for the grid-based Xin'anjiang model based on compact indexes and CUDA graphs. Background Technology

[0002] Floods are characterized by their sudden onset and wide-ranging impact, often causing significant casualties and economic losses. In recent years, extreme rainfall has increased, making flood processes more complex and demanding higher levels of timeliness and accuracy in forecasting. Conducting high-resolution flood simulations and real-time forecasts can provide a scientific basis for reservoir management and flood control decision-making.

[0003] The gridded Xin'anjiang model is a distributed hydrological model that discretizes the watershed into regular grid cells. At the grid scale, it reflects the spatial heterogeneity of rainfall, evapotranspiration, and underlying surface parameters. It calculates the process time-by-time according to evapotranspiration, runoff generation, water source distribution, slope confluence, and river confluence, thus more meticulously depicting the formation and evolution of floods and improving the simulation accuracy and spatiotemporal interpretation capabilities of rainstorm floods. However, under conditions of large watersheds, high resolution, or long-term operation, the number of grid cells and the scale of state variables increase dramatically. Furthermore, river confluence is constrained by upstream and downstream dependencies, significantly increasing the overall computation time. Simultaneously, traditional Java-based models are susceptible to heap memory limitations in scenarios involving large arrays and multi-time-period state caching, posing engineering risks of excessive memory consumption or even overflow, affecting stable operational use.

[0004] At the level of general-purpose GPU parallel computing, while existing work has explored parallel acceleration under the CUDA framework, significant shortcomings remain in efficient scheduling and low-overhead execution. Previous studies have proposed GPU parallelization schemes to address the data circular dependency problem, but these methods rely on manual reconstruction of the original serial loop structure. This requires repeatedly launching the CUDA kernel function in each computation period, failing to utilize mechanisms such as computation graphs to eliminate the CPU-side scheduling and synchronization overhead caused by repeated launches. This overhead is significant in long-term iterative scenarios. Multi-GPU parallel research on hydrological and hydrodynamic processes achieves cross-GPU computation through region decomposition and overlapping sub-region communication, but its parallel granularity is coarse, mainly relying on CUDA streams or P2P communication for inter-region data exchange. It does not address kernel launch sequence solidification and thread organization optimization, and still faces efficiency losses due to frequent data allocation and transmission and step-by-step CPU scheduling during long-duration simulations. Existing GPU parallelization studies on the Xin'anjiang grid model have achieved some acceleration in independently parallel processes such as evapotranspiration, runoff generation, and slope confluence. However, for river confluence constrained by upstream and downstream dependencies, the method still relies on repeatedly scanning the entire grid layer by layer and determining the level in the kernel function. This results in an excessively high proportion of invalid threads, discontinuous memory access, and limited thread utilization and memory access efficiency. Furthermore, the lack of unified management of GPU memory allocation during time-period iterations and the failure to achieve asynchronous overlap between data transmission and computation significantly reduce the acceleration gains in long-duration simulations. In addition, none of the aforementioned existing methods effectively support concurrent and consistent injection of inflows from upstream boundaries, making it difficult to meet the actual needs of operational flood forecasting scenarios. Summary of the Invention

[0005] Purpose of the invention: The purpose of this invention is to address the shortcomings of existing technologies and provide a parallel computing method for the grid-based Xin'anjiang model based on hierarchical compact indexes and executable CUDA graphs. This method aims to solve problems such as the limited parallel efficiency caused by river confluence dependency constraints in GPU-accelerated applications of the grid-based Xin'anjiang model, the large overhead of data transfer between CPU and GPU and multiple kernel startups, and the insufficient guarantee of upstream boundary external inflow coupling and concurrent accumulation consistency.

[0006] Technical solution: The parallel method for the gridded Xin'anjiang model based on compact indexing and CUDA graphs of the present invention includes the following steps:

[0007] Step S1: Obtain the flow direction, computation level and effective mesh mask of the watershed grid cell, generate a one-dimensional index array next for the downstream grid cell, and construct a layer compact index table, so that during layer-by-layer computation, only the effective mesh of the current layer is mapped to a parallel thread.

[0008] Step S2: On the GPU side, based on the memory pool technology, allocate contiguous global memory for model parameters, state variables and other model calculation variables and river confluence accumulators at one time, and divide each storage area according to offset to complete the initialization; then, establish a dual buffer for driver data paged memory on the CPU side, and establish a corresponding dual buffer for video memory on the GPU side, and create at least one asynchronous transport stream and computation stream to realize the overlap of driver data transmission and computation.

[0009] Step S3: For each calculation period, the precipitation and evaporation data of the next period are asynchronously transmitted to the unoccupied drive data buffer on the GPU in the asynchronous transport stream. Simultaneously, the evapotranspiration calculation for the current period is executed in parallel in the calculation stream to update the tension water state, thereby achieving asynchronous overlap between data transmission and model calculation. The specific process is as follows:

[0010] Step S3.1: Before time step t begins calculation, confirm that the driving data for time step t has been transferred and ready in one of the GPU's dual memory buffers (e.g., Device Buffer A), and point the current driving input pointer to that buffer.

[0011] Step S3.2: While the computational stream is performing the model calculation at time step t, the precipitation field and evaporation field data at time step t+1 are asynchronously transmitted to another free buffer (e.g., Device Buffer B) in the GPU's dual memory buffer via an asynchronous transmission stream.

[0012] Step S3.3: Before entering the calculation of time step t + 1, stream synchronization is performed to ensure that the data transmission of the next time period is completed. Then, the drive input pointer is switched to point to the free buffer (Device Buffer B) that has been loaded with data of time step t + 1, thereby realizing the asynchronous overlap of data transmission and calculation.

[0013] Step S3.4: In the computational flow, evaporation calculation is performed using grid cells as parallel computing units. Each parallel thread locates the grid cell to be processed based on the global index and performs computational masking on invalid grid cells based on the effective grid mask.

[0014] Step S3.5: Read the tension water state passed from the previous time step from the state variable storage area, combine it with the rainfall field and evaporation field data of time step t, execute the three-layer evapotranspiration calculation process, and deduct the tension water corresponding to each layer according to the calculation results to update the current tension water state.

[0015] Finally, the calculated total evapotranspiration is written into the intermediate result storage area as input for the subsequent step S4;

[0016] Step S4: Perform runoff generation and water source calculation in parallel in the computational stream to calculate surface runoff, interflow and groundwater runoff, and update the free water state, surface runoff, interflow and groundwater runoff.

[0017] Step S5: In the computational flow, the slope runoff calculation is performed in parallel based on the surface runoff, interflow and groundwater runoff obtained in step S4, and the slope inflow of each grid cell into the river channel is calculated.

[0018] Step S6: Before the river confluence calculation, the obtained slope inflow is used as the local river inflow item of each grid cell for subsequent parallel river confluence calculation. The upstream boundary inflow is accumulated into the river confluence accumulator of the corresponding grid cell using atomic addition to ensure the consistency of multiple inflow points.

[0019] The specific method for upstream boundary inflow accumulation is as follows: First, obtain the spatial location and inflow process of at least one upstream boundary inflow point; then determine the target channel grid cell corresponding to the inflow point; when the grid cell corresponding to the inflow point is an invalid water system grid cell, perform a neighborhood search to obtain candidate channel grid cells, and select effective grids from the candidate grids, taking the effective candidate grid with the largest catchment area feature value as the target channel grid cell; finally, accumulate the inflow value of the current calculation period into the channel confluence accumulator of the target channel grid cell using double-precision atomic addition.

[0020] Step S7: Based on the hierarchical compact index table, slope inflow and upstream inflow in the river confluence accumulator, perform parallel calculation of river confluence layer by layer from low to high calculation level;

[0021] First, the startup sequence of the river confluence kernel function at each level is instantiated into an executable CUDA graph through a flow capture mechanism. In each computation period, the GPU time step parameters and related pointers are updated, and the resulting executable CUDA graph is started. Within each level, the Muskingum method is used to calculate the grid outflow, and the outflow is accumulated to the river confluence accumulator of the downstream grid using double-precision atomic addition, thereby reducing the scheduling overhead caused by the CPU starting layer by layer.

[0022] Step S8: Repeat steps S3 to S7 until all calculation periods are completed, and output the flow results of the watershed outlet grid cell.

[0023] Furthermore, the details of step S1 are as follows:

[0024] Step S1.1: Extract the watershed boundary based on the digital elevation model (DEM) to obtain an effective grid mask, which is used to distinguish the effective grid within the watershed from the invalid grid outside the watershed.

[0025] Expand the coordinates (row, col) of the watershed grid cell into a one-dimensional index i in row-major order, as shown in the following formula:

[0026] ;

[0027] Where row is the row coordinate of the current grid cell, col is the column coordinate of the current grid cell, and ncols is the total number of columns of the watershed grid cells;

[0028] Step S1.2: Determine the downstream grid coordinates (rowd, cold) of each grid cell based on the flow direction code, and generate a one-dimensional index array next for the downstream grid, as shown in the following formula:

[0029] ;

[0030] Where rowd is the row coordinate of the downstream grid cell and cold is the column coordinate of the downstream grid cell;

[0031] Step S1.3: When the downstream grid coordinates (rowd, cold) exceed the range of watershed grid coordinates or correspond to invalid values ​​in the effective grid mask, it is determined that there are no valid downstream grids for the grid. Valid grids that do not have valid downstream grids or whose coordinates are consistent with the preset outlet position are determined as outlet grids. Subsequently, the downstream index of the outlet grid and the downstream index of the grid cells that do not have valid downstream grids are assigned to an invalid flag or their own index i, and the one-dimensional index outlet_index of the outlet grid is recorded.

[0032] Step S1.4: Based on the downstream one-dimensional index array next, count the number of upstream incoming grid cells for each effective grid cell. Use the topological sorting algorithm to determine the computational level l of each grid cell. The level of the source grid cell with no upstream incoming cells is set to 1. Recursively calculate the level of the downstream grid cell by taking the maximum value of the levels of all its upstream grid cells plus 1. After completing the level calculation of all effective grid cells, determine the maximum level value L of each effective grid cell level value. max ;

[0033] Step S1.5: For each level l = 1, ..., L max Filter the one-dimensional indexes of grid cells that meet the requirements of level 1 and are valid grids, form a compact index array IndexList[l], and record its length count_l, thereby constructing the layer compact index table Layer Compact Index.

[0034] Furthermore, the specific process of step S2 is as follows:

[0035] Step S2.1: On the GPU side, a contiguous global memory space is allocated at once using memory pooling technology to store model parameters, state variables and other model calculation variables, as well as the river confluence accumulator.

[0036] Step S2.2: Set the byte alignment size (Align) for video memory allocation, and calculate the alignment space size corresponding to each model's computational variables in sequence. The alignment calculation formula is:

[0037] ;

[0038] in, This represents the actual number of bytes required before variable alignment. Indicates the introduction After calculating the compensation margin, the integer part is rounded down, which means the integer number of standard aligned data blocks that the variable needs to occupy. Multiplying this by Align will give the final number of bytes allocated to satisfy the GPU's merged memory access characteristics.

[0039] Subsequently, the interval is divided by accumulation, with 0 as the initial offset. The offset of the current variable in the video memory pool is equal to the sum of the offset of the previous variable and the size of its alignment space.

[0040] Finally, the base address pointer of the contiguous global video memory is added to the offset of each variable, and the independent video memory pointer of each variable is directly divided through address addressing, and the initialization is completed;

[0041] Step S2.3: Establish a dual buffer for driver data on the CPU side using a page-locked memory allocation method. The dual buffer includes Host Buffer A and Host Buffer B.

[0042] Step S2.4: Establish a dual memory buffer on the GPU corresponding to the obtained Host Buffer A and Host Buffer B. The dual memory buffer includes Device Buffer A and Device Buffer B, and the driving data includes at least rainfall field and evaporation field data.

[0043] Step S2.5: Create at least one asynchronous transport stream and at least one computation stream.

[0044] Furthermore, the specific process of parallel calculation and updating of the free water state for runoff generation and water source distribution in step S4 is as follows:

[0045] Step S4.1: In the computational flow, runoff generation calculation is performed using grid cells as parallel computing units, and the soil moisture content status is updated;

[0046] Step S4.2: Divide the runoff results into surface runoff, interflow and groundwater runoff, and update the free water content status;

[0047] Step S4.3: Write the water source division results into the intermediate result storage area of ​​the model calculation variables as input for the subsequent step S5.

[0048] Furthermore, the specific process of parallel calculation of slope confluence in step S5 is as follows:

[0049] Step S5.1: Perform slope runoff calculation in the computation flow using grid cells as parallel computing units, and calculate the slope inflow based on the surface runoff, interflow and groundwater runoff obtained in step S4;

[0050] Step S5.2: Update the slope confluence status;

[0051] Step S5.3: Write the slope inflow into the intermediate result storage area of ​​the model calculation variables, as the input for the subsequent step S7 river confluence.

[0052] Furthermore, the specific process of performing parallel computation of river confluence based on the hierarchical compact index table and launching the executable CUDA graph in step S7 is as follows:

[0053] Step S7.1: Enable the stream capture mechanism during the model initialization phase, according to the hierarchy l = 1, ..., L max The startup sequence of the river confluence kernel function is recorded sequentially, and for each level, the corresponding compact index array IndexList[l] and its length count_l are read, and thread mapping is established according to the index array;

[0054] Step S7.2: The thread locates the current grid cell based on IndexList[l], reads the slope inflow of the grid cell and the upstream inflow in the channel confluence accumulator, calculates the channel outflow of the current grid cell using the Muskingum method, and accumulates the channel outflow into the channel confluence accumulator of the downstream grid cell using double-precision atomic addition based on the downstream one-dimensional index array.

[0055] Step S7.3: Complete levels 1 to L max After the kernel function initiates sequence capture, it ends flow capture and instantiates the hierarchical river confluence computation graph into an executable CUDA graph;

[0056] Step S7.4: During each computation period iteration, update the time step parameters and driver input pointers corresponding to the executable CUDA graph, and start the updated executable CUDA graph in the computation flow to adapt to the double buffering switching of driver data and avoid repeatedly building the execution graph.

[0057] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0058] (1) This invention uses a hierarchical compact index table to directly address and calculate the river grid layer by layer, which effectively avoids the ineffective scheduling caused by full grid scanning and threshold-based iterative advancement. By centrally managing data through a unified memory pool and cooperating with page-locked double-buffered asynchronous transmission, data preparation and calculation overlap, reducing the data transmission wait between the CPU and GPU.

[0059] (2) By fixing the kernel function startup sequence of the hierarchical river confluence to CUDA graph execution and combining upstream inflow injection and atomic accumulation consistency control, the entire process of river confluence is stably executed on the GPU, avoiding the synchronization and scheduling overhead caused by the CPU participating in hierarchical confluence with too few grids.

[0060] (3) Under the premise of maintaining the physical mechanism and computational chain integrity of the grid Xin'anjiang model, this invention is easy to connect directly with upstream boundary inflow data, realizes efficient parallel computing for large-scale grid and long time series simulation, can meet the engineering requirements of real-time flood forecast for computational timeliness and result consistency, and has good engineering application value and promotion prospects. Attached Figure Description

[0061] Figure 1 This is a schematic diagram of the overall process of the method of the present invention;

[0062] Figure 2 This is a schematic diagram illustrating the construction of the downstream one-dimensional index array and hierarchical compact index table in the embodiment;

[0063] Figure 3 This is a schematic diagram illustrating the double buffering of the locked page and the asynchronous transmission and calculation of the driving data in the embodiment;

[0064] Figure 4 This is a comparison diagram of the traditional layer-by-layer startup method and the executable CUDA graph method in the embodiment;

[0065] Figure 5 This is a comparison diagram of the flow process lines of the GPU parallel implementation and the CPU serial implementation in the embodiment;

[0066] Figure 6 This is a comparison chart of the computation time of the GPU parallel implementation and the CPU serial implementation in the example. Detailed Implementation

[0067] The technical solution of the present invention will be described in detail below, but the scope of protection of the present invention is not limited to the embodiments described.

[0068] This invention innovates collaboratively from two levels: GPU parallel scheduling mechanism and hydrological model application. At the parallel scheduling level, a hierarchical compact index table is proposed for direct layer-by-layer addressing of the river channel grid, fundamentally eliminating the invalid thread problem caused by full-grid scanning and layer-by-layer judgment. This ensures that each thread in the river confluence calculation is mapped to an effective grid, improving thread utilization and memory access continuity. Simultaneously, a flow capture mechanism is used to solidify the layer-by-layer kernel function startup sequence into an executable CUDA graph. In long-term iterations, each time period only needs to update parameters and start the execution graph, avoiding repeated scheduling and layer-by-layer synchronization overhead on the CPU side—a key technology not yet adopted in existing GPU parallel computing methods. Then, based on a one-time continuous memory allocation from the memory pool, combined with CPU-side paged memory double buffering and GPU-side double buffering, pipeline overlap between asynchronous data transmission and model computation is achieved, effectively eliminating data transmission latency.

[0069] At the level of hydrological model application, this invention achieves concurrent and consistent accumulation of multiple inflow points through double-precision atomic addition, supporting the injection of inflows from outside the upstream boundary. This overcomes the shortcomings of existing grid-based Xin'anjiang model in terms of the applicability of GPU parallel operation in operational scenarios. These innovations enable this invention to achieve significant acceleration in large-scale grid and long-term simulations while maintaining consistency between the outlet flow results and serial implementations, demonstrating good engineering applicability and promotional value for real-time flood forecasting.

[0070] like Figure 1 As shown, the parallel method for the gridded Xin'anjiang model based on compact indexing and CUDA graphs of the present invention includes the following steps:

[0071] Step S1: Obtain the flow direction, computation level and effective mesh mask of the watershed grid cell, generate a one-dimensional index array next for the downstream grid cell, and construct a layer compact index table, so that during layer-by-layer computation, only the effective mesh of the current layer is mapped to a parallel thread.

[0072] Step S2: On the GPU side, based on the memory pool technology, allocate contiguous global memory for model parameters, state variables and other model calculation variables and river confluence accumulators at one time, and divide each storage area according to offset to complete the initialization; then, establish a dual buffer for driver data paged memory on the CPU side, and establish a corresponding dual buffer for video memory on the GPU side, and create at least one asynchronous transport stream and computation stream to realize the overlap of driver data transmission and computation.

[0073] Step S3: For each calculation period, the precipitation and evaporation data for the next period are asynchronously transmitted to the unoccupied drive data buffer on the GPU in the asynchronous transport stream. Simultaneously, the evapotranspiration calculation for the current period is executed in parallel in the calculation stream to update the tension water state, thereby achieving asynchronous overlap between data transmission and model calculation; Figure 3 As shown, the specific process is as follows:

[0074] Step S3.1: Before time step t begins calculation, confirm that the driving data for time step t has been transferred and ready in one of the GPU's dual memory buffers (e.g., Device Buffer A), and point the current driving input pointer to that buffer.

[0075] Step S3.2: While the computational stream is performing the model calculation at time step t, the precipitation field and evaporation field data at time step t+1 are asynchronously transmitted to another free buffer (e.g., Device Buffer B) in the GPU's dual memory buffer via an asynchronous transmission stream.

[0076] Step S3.3: Before entering the calculation of time step t + 1, stream synchronization is performed to ensure that the data transmission of the next time period is completed. Then, the drive input pointer is switched to point to the free buffer (Device Buffer B) that has been loaded with data of time step t + 1, thereby realizing the asynchronous overlap of data transmission and calculation.

[0077] Step S3.4: In the computational flow, evaporation calculation is performed using grid cells as parallel computing units. Each parallel thread locates the grid cell to be processed based on the global index and performs computational masking on invalid grid cells based on the effective grid mask.

[0078] Step S3.5: Read the tension water state passed from the previous time step from the state variable storage area, combine it with the rainfall field and evaporation field data of time step t, execute the three-layer evapotranspiration calculation process, and deduct the tension water corresponding to each layer according to the calculation results to update the current tension water state.

[0079] Finally, the calculated total evapotranspiration is written into the intermediate result storage area as input for the subsequent step S4;

[0080] Step S4: Perform runoff generation and water source calculation in parallel in the computational stream to calculate surface runoff, interflow and groundwater runoff, and update the free water state, surface runoff, interflow and groundwater runoff.

[0081] Step S5: In the computational flow, the slope runoff calculation is performed in parallel based on the surface runoff, interflow and groundwater runoff obtained in step S4, and the slope inflow of each grid cell into the river channel is calculated.

[0082] Step S6: Before the river confluence calculation, the obtained slope inflow is used as the local river inflow item of each grid cell for subsequent parallel river confluence calculation. The upstream boundary inflow is accumulated into the river confluence accumulator of the corresponding grid cell using atomic addition to ensure the consistency of multiple inflow points.

[0083] The specific method for upstream boundary inflow accumulation is as follows: First, obtain the spatial location and inflow process of at least one upstream boundary inflow point; then determine the target channel grid cell corresponding to the inflow point; when the grid cell corresponding to the inflow point is an invalid water system grid cell, perform a neighborhood search to obtain candidate channel grid cells, and select effective grids from the candidate grids, taking the effective candidate grid with the largest catchment area feature value as the target channel grid cell; finally, accumulate the inflow value of the current calculation period into the channel confluence accumulator of the target channel grid cell using double-precision atomic addition.

[0084] Step S7: Based on the hierarchical compact index table, slope inflow and upstream inflow in the river confluence accumulator, perform parallel calculation of river confluence layer by layer from low to high calculation level;

[0085] First, the startup sequence of the river confluence kernel function at each level is instantiated into an executable CUDA graph through a flow capture mechanism. In each computation period, the GPU time step parameters and related pointers are updated, and the resulting executable CUDA graph is started. Within each level, the Muskingum method is used to calculate the grid outflow, and the outflow is accumulated to the river confluence accumulator of the downstream grid using double-precision atomic addition, thereby reducing the scheduling overhead caused by the CPU starting layer by layer.

[0086] Step S8: Repeat steps S3 to S7 until all calculation periods are completed, and output the flow results of the watershed outlet grid cell.

[0087] like Figure 2 As shown, the details of step S1 are as follows:

[0088] Step S1.1: Extract the watershed boundary based on the digital elevation model (DEM) to obtain an effective grid mask, which is used to distinguish the effective grid within the watershed from the invalid grid outside the watershed.

[0089] Expand the coordinates (row, col) of the watershed grid cell into a one-dimensional index i in row-major order, as shown in the following formula:

[0090] ;

[0091] Where row is the row coordinate of the current grid cell, col is the column coordinate of the current grid cell, and ncols is the total number of columns of the watershed grid cells;

[0092] Step S1.2: Determine the downstream grid coordinates (rowd, cold) of each grid cell based on the flow direction code, and generate a one-dimensional index array next for the downstream grid, as shown in the following formula:

[0093] ;

[0094] Where rowd is the row coordinate of the downstream grid cell and cold is the column coordinate of the downstream grid cell;

[0095] Step S1.3: When the downstream grid coordinates (rowd, cold) exceed the range of watershed grid coordinates or correspond to invalid values ​​in the effective grid mask, it is determined that there are no valid downstream grids for the grid. Valid grids that do not have valid downstream grids or whose coordinates are consistent with the preset outlet position are determined as outlet grids. Subsequently, the downstream index of the outlet grid and the downstream index of the grid cells that do not have valid downstream grids are assigned to an invalid flag or their own index i, and the one-dimensional index outlet_index of the outlet grid is recorded.

[0096] Step S1.4: Based on the downstream one-dimensional index array next, count the number of upstream incoming grid cells for each effective grid cell. Use the topological sorting algorithm to determine the computational level l of each grid cell. The level of the source grid cell with no upstream incoming cells is set to 1. Recursively calculate the level of the downstream grid cell by taking the maximum value of the levels of all its upstream grid cells plus 1. After completing the level calculation of all effective grid cells, determine the maximum level value L of each effective grid cell level value. max ;

[0097] Step S1.5: For each level l = 1, ..., L max Filter the one-dimensional indexes of grid cells that meet the requirements of level 1 and are valid grids, form a compact index array IndexList[l], and record its length count_l, thereby constructing the layer compact index table Layer Compact Index.

[0098] The specific process of step S2 in this embodiment is as follows:

[0099] Step S2.1: On the GPU side, a contiguous global memory space is allocated at once using memory pooling technology to store model parameters, state variables and other model calculation variables, as well as the river confluence accumulator.

[0100] Step S2.2: Set the byte alignment size (Align) for video memory allocation, and calculate the alignment space size corresponding to each model's computational variables in sequence. The alignment calculation formula is:

[0101] ;

[0102] in, This represents the actual number of bytes required before variable alignment. Indicates the introduction After calculating the compensation margin, the integer part is rounded down, which means the integer number of standard aligned data blocks that the variable needs to occupy. Multiplying this by Align will give the final number of bytes allocated to satisfy the GPU's merged memory access characteristics.

[0103] Subsequently, the interval is divided by accumulation, with 0 as the initial offset. The offset of the current variable in the video memory pool is equal to the sum of the offset of the previous variable and the size of its alignment space.

[0104] Finally, the base address pointer of the contiguous global video memory is added to the offset of each variable, and the independent video memory pointer of each variable is directly divided through address addressing, and the initialization is completed;

[0105] Step S2.3: Establish a dual buffer for driver data on the CPU side using a page-locked memory allocation method. The dual buffer includes Host Buffer A and Host Buffer B.

[0106] Step S2.4: Establish a dual memory buffer on the GPU corresponding to the obtained Host Buffer A and Host Buffer B. The dual memory buffer includes Device Buffer A and Device Buffer B, and the driving data includes at least rainfall field and evaporation field data.

[0107] Step S2.5: Create at least one asynchronous transport stream and at least one computation stream.

[0108] The specific process of parallel calculation of runoff generation and water source distribution and updating of free water state in step S4 of this embodiment is as follows:

[0109] Step S4.1: In the computational flow, runoff generation calculation is performed using grid cells as parallel computing units, and the soil moisture content status is updated;

[0110] Step S4.2: Divide the runoff results into surface runoff, interflow and groundwater runoff, and update the free water content status;

[0111] Step S4.3: Write the water source division results into the intermediate result storage area of ​​the model calculation variables as input for the subsequent step S5.

[0112] The specific process of parallel calculation of slope confluence in step S5 of this embodiment is as follows:

[0113] Step S5.1: Perform slope runoff calculation in the computation flow using grid cells as parallel computing units, and calculate the slope inflow based on the surface runoff, interflow and groundwater runoff obtained in step S4;

[0114] Step S5.2: Update the slope confluence status;

[0115] Step S5.3: Write the slope inflow into the intermediate result storage area of ​​the model calculation variables, as the input for the subsequent step S7 river confluence.

[0116] The specific process of performing parallel computation of river confluence based on the hierarchical compact index table and starting the executable CUDA graph in step S7 of this embodiment is as follows:

[0117] Step S7.1: Enable the stream capture mechanism during the model initialization phase, according to the hierarchy l = 1, ..., L max The startup sequence of the river confluence kernel function is recorded sequentially, and for each level, the corresponding compact index array IndexList[l] and its length count_l are read, and thread mapping is established according to the index array;

[0118] Step S7.2: The thread locates the current grid cell based on IndexList[l], reads the slope inflow of the grid cell and the upstream inflow in the channel confluence accumulator, calculates the channel outflow of the current grid cell using the Muskingum method, and accumulates the channel outflow into the channel confluence accumulator of the downstream grid cell using double-precision atomic addition based on the downstream one-dimensional index array.

[0119] Step S7.3: Complete levels 1 to L max After the kernel function initiates sequence capture, it ends flow capture and instantiates the hierarchical river confluence computation graph into an executable CUDA graph;

[0120] Step S7.4: During each computation period iteration, update the time step parameters and driver input pointers corresponding to the executable CUDA graph, and start the updated executable CUDA graph in the computation flow to adapt to the double buffering switching of driver data and avoid repeatedly building the execution graph.

[0121] To verify the feasibility and beneficial effects of the technical solution of this invention, this embodiment takes the Tunxi River Basin as an example, selecting the flood process from 08:00 on June 18, 2024 to 00:00 on July 7, 2024 as the test object. The specific method is as follows:

[0122] Step S1: Obtain the flow direction, computation level, and effective mesh mask of the watershed grid cell, generate a one-dimensional index array of downstream grid cells, and construct a hierarchical compact index table so that during layer-by-layer computation, only the effective mesh of the current layer is mapped to a parallel thread.

[0123] Step S2: On the GPU side, based on the memory pool technology, allocate contiguous global memory for model parameters, state variables and other model calculation variables and river confluence accumulators at one time, and divide each storage area according to offset to complete the initialization; on the CPU side, establish a driver data paged memory double buffer, and on the GPU side, establish a corresponding video memory double buffer, and at the same time create at least one asynchronous transport stream and calculation stream to realize the overlap of driver data transmission and calculation;

[0124] Figure 2 (a) in the diagram is a schematic diagram of the mesh flow topology in this embodiment. Figure 2 (b) in the diagram is a schematic diagram of the one-dimensional index and calculation hierarchy in this implementation. Figure 2 (c) in the example is a schematic diagram of the hierarchical compact index table and thread mapping;

[0125] Step S3: For each calculation period, the rainfall field and evaporation field data of the next period are asynchronously transmitted to the unoccupied drive data buffer on the GPU in the asynchronous transport stream. At the same time, the evapotranspiration calculation of the current period is executed in parallel in the calculation stream to update the tension water state, thereby realizing the asynchronous overlap of data transmission and model calculation.

[0126] Step S4: Perform runoff generation and water source calculations in parallel within the computational stream to calculate surface runoff, interflow, and groundwater runoff, and update the free water state;

[0127] Step S5: Perform slope runoff calculations in parallel within the computational stream to calculate the slope inflow of each grid cell into the river channel;

[0128] Step S6: Before calculating the river confluence, the upstream boundary inflow is accumulated into the river confluence accumulator of the corresponding grid cell using atomic addition to ensure the consistency of multiple inflow points concurrently.

[0129] Step S7: Based on the hierarchical compact index table, perform parallel computation of river confluence layer by layer from low to high computation level; first, instantiate the startup sequence of the river confluence kernel function of each level into an executable CUDA graph through the flow capture mechanism; in each computation period, update the GPU time step parameters and related pointers, and start the executable CUDA graph; within each level, use the Muskingum method to calculate the grid outflow, and use double-precision atomic addition to accumulate the outflow to the river confluence accumulator of the downstream grid, thereby reducing the scheduling overhead caused by the CPU starting layer by layer;

[0130] Step S8: Repeat steps S3 to S7 until all calculation periods are completed, and output the flow results of the watershed outlet grid cell.

[0131] This embodiment uses the same input data and parameter settings to compare calculations on CPU serial implementation (Java) and GPU parallel implementation (this invention), respectively. Figure 4 For CPU serial implementation, under the default heap memory configuration, there is a problem of insufficient memory causing computation interruption. The computation can only be completed after the Java Virtual Machine heap memory limit is adjusted to 44GB.

[0132] In contrast, the GPU parallel implementation of this invention employs a memory organization method with dual buffering of memory pools and locked pages, which can stably complete full-time calculations under the same computational examples and output the flow results of the watershed outlet grid cells; its outlet flow process line is consistent with the CPU serial implementation results (e.g. Figure 5 ), while the average computation time is significantly reduced (e.g. Figure 6 ), the acceleration ratio is about 20 times.

[0133] In summary, this invention explicitly defines the grid topology and computational hierarchy using a hierarchical compact index table, achieving compact thread mapping of the effective grid layer by layer. Combined with one-time continuous allocation of the memory pool and asynchronous transmission and computational overlap via page-locked double buffering, it significantly reduces memory management and data transport overhead. Furthermore, by constructing an executable CUDA graph through flow capture, it fixes the layer-by-layer initiation of layered river confluence into a single graph initiation, reducing the cost of repeated scheduling on the CPU side. This invention can stably complete full-time calculations while maintaining consistency between the outflow results and the CPU serial implementation, with a significantly reduced average execution time, demonstrating its engineering applicability and promotional value in large-scale grid, long-duration simulation, and real-time flood forecasting scenarios.

Claims

1. A parallel method for the grid-based Xin'anjiang model based on compact indexing and CUDA graphs, characterized in that, Includes the following steps: Step S1: Obtain the flow direction, computation level and effective mesh mask of the watershed grid cell, generate a one-dimensional index array next for the downstream grid cell, and construct a layer compact index table, so that during layer-by-layer computation, only the effective mesh of the current layer is mapped to a parallel thread. Step S2: On the GPU side, based on the memory pool technology, allocate contiguous global memory for the model calculation variables and the river confluence accumulator at one time, and divide each storage area according to the offset to complete the initialization; then, establish a dual buffer for driving data page memory on the CPU side, and establish a corresponding dual buffer for video memory on the GPU side, while creating at least one asynchronous transport stream and calculation stream. Step S3: For each calculation period, the precipitation and evaporation data for the next period are asynchronously transmitted to the unoccupied driver data buffer on the GPU in the asynchronous transport stream. Simultaneously, the evapotranspiration calculation for the current period is executed in parallel in the calculation stream to update the tension water state. The specific process is as follows: Step S3.1: Before time step t begins calculation, confirm that the driving data for time step t has been transferred and ready in one of the GPU's dual memory buffers, and point the current driving input pointer to that buffer. Step S3.2: While the computational stream is performing the model calculation at time step t, the precipitation field and evaporation field data at time step t+1 are asynchronously transmitted to another free buffer in the GPU's dual memory buffer via an asynchronous transmission stream. Step S3.3: Before entering the calculation of time step t + 1, stream synchronization is performed to ensure that the data transmission of the next time period is completed. Then, the drive input pointer is switched to point to the free buffer that has been loaded with the data of time step t + 1, thereby realizing the asynchronous overlap of data transmission and calculation. Step S3.4: In the computation flow, evapotranspiration calculation is performed using grid cells as parallel computing units. Each parallel thread locates the grid cell to be processed based on the global index and performs computation masking on invalid grid cells based on the effective grid mask. Step S3.5: Read the tension water state passed from the previous time step from the state variable storage area, combine it with the rainfall field and evaporation field data of time step t, execute the three-layer evapotranspiration calculation process, and deduct the tension water corresponding to each layer according to the calculation results to update the current tension water state. Finally, the calculated total evapotranspiration is written into the intermediate result storage area as input for the subsequent step S4; Step S4: Perform runoff generation and water source calculation in parallel in the computational stream to calculate surface runoff, interflow and groundwater runoff, and update the free water state, surface runoff, interflow and groundwater runoff. Step S5: In the computational flow, the slope runoff calculation is performed in parallel based on the surface runoff, interflow and groundwater runoff obtained in step S4, and the slope inflow of each grid cell into the river channel is calculated. Step S6: Before the river confluence calculation, the obtained slope inflow is used as the local river inflow item of each grid cell for subsequent parallel river confluence calculation. The upstream boundary inflow is accumulated into the river confluence accumulator of the corresponding grid cell using atomic addition to ensure the consistency of multiple inflow points. The specific method for upstream boundary inflow accumulation is as follows: First, obtain the spatial location and inflow process of at least one upstream boundary inflow point; then determine the target channel grid cell corresponding to the inflow point; when the grid cell corresponding to the inflow point is an invalid water system grid cell, perform a neighborhood search to obtain candidate channel grid cells, and select effective grids from the candidate grids, taking the effective candidate grid with the largest catchment area feature value as the target channel grid cell; finally, accumulate the inflow value of the current calculation period into the channel confluence accumulator of the target channel grid cell using double-precision atomic addition. Step S7: Based on the hierarchical compact index table, slope inflow and upstream inflow in the river confluence accumulator, perform parallel calculation of river confluence layer by layer from low to high calculation level; First, the startup sequence of the river confluence kernel function at each level is instantiated into an executable CUDA graph through a flow capture mechanism. In each computation period, the GPU time step parameters and related pointers are updated, and the resulting executable CUDA graph is started. Within each level, the Muskingum method is used to calculate the grid outflow, and the outflow is accumulated to the river confluence accumulator of the downstream grid using double-precision atomic addition, thereby reducing the scheduling overhead caused by the CPU starting layer by layer. Step S8: Repeat steps S3 to S7 until all calculation periods are completed, and output the flow results of the watershed outlet grid cell.

2. The parallel method for the gridded Xin'anjiang model based on compact indexing and CUDA graphs according to claim 1, characterized in that, The details of step S1 are as follows: Step S1.1: Extract the watershed boundary based on the digital elevation model (DEM) to obtain an effective grid mask, which is used to distinguish the effective grid within the watershed from the invalid grid outside the watershed. Expand the coordinates (row, col) of the watershed grid cell into a one-dimensional index i in row-major order, as shown in the following formula: ； Where row is the row coordinate of the current grid cell, col is the column coordinate of the current grid cell, and ncols is the total number of columns of the watershed grid cells; Step S1.2: Determine the downstream grid coordinates (rowd, cold) of each grid cell based on the flow direction encoding, and generate a one-dimensional index array next for the downstream grid, as shown in the following formula: ； Where rowd is the row coordinate of the downstream grid cell and cold is the column coordinate of the downstream grid cell; Step S1.3: When the downstream grid coordinates (rowd, cold) exceed the range of watershed grid coordinates or correspond to invalid values ​​in the effective grid mask, it is determined that there are no valid downstream grids for the grid. Valid grids that do not have valid downstream grids or whose coordinates are consistent with the preset outlet position are determined as outlet grids. Subsequently, the downstream index of the outlet grid and the downstream index of the grid cells that do not have valid downstream grids are assigned to an invalid flag or their own index i, and the one-dimensional index outlet_index of the outlet grid is recorded. Step S1.4, according to the downstream one-dimensional index array next, count the upstream converging grid number of each valid grid cell, and determine the calculation level l of each grid cell by using the topological sorting algorithm, wherein the level of the source grid cell without upstream convergence is set to 1, and the level of the downstream grid cell is recursively calculated along the flow direction, and the level of the downstream grid cell is the maximum value of the levels of all upstream grid cells plus 1, after the level calculation of all valid grid cells is completed, the maximum value in the level values of each valid grid cell is determined as the maximum level value L max ; Step S1.5: For each level l = 1, ..., L max Filter the one-dimensional indexes of grid cells that meet the requirements of level 1 and are valid grids, form a compact index array IndexList[l], and record its length count_l, thereby constructing the layer compact index table Layer Compact Index.

3. The parallel method for the gridded Xin'anjiang model based on compact indexing and CUDA graphs according to claim 1, characterized in that, The specific process of step S2 is as follows: Step S2.1: On the GPU side, a contiguous global memory space is allocated at once using memory pooling technology to store model parameters, state variables and other model calculation variables, as well as the river confluence accumulator. Step S2.2: Set the byte alignment size (Align) for video memory allocation, and calculate the alignment space size corresponding to each model's computational variables in sequence. The alignment calculation formula is: ； in, This represents the actual number of bytes required before variable alignment. Indicates the introduction After calculating the compensation margin, round down. Subsequently, the interval is divided by accumulation, with 0 as the initial offset. The offset of the current variable in the video memory pool is equal to the sum of the offset of the previous variable and the size of its alignment space. Finally, the base address pointer of the contiguous global video memory is added to the offset of each variable, and the independent video memory pointer of each variable is directly divided through address addressing, and the initialization is completed; Step S2.3: Establish a dual buffer for driver data on the CPU side using a page-locked memory allocation method. The dual buffer includes Host Buffer A and Host Buffer B. Step S2.4: Establish a dual memory buffer on the GPU corresponding to the obtained Host Buffer A and Host Buffer B. The dual memory buffer includes Device Buffer A and Device Buffer B, and the driving data includes at least rainfall field and evaporation field data. Step S2.5: Create at least one asynchronous transport stream and at least one computation stream.

4. The parallel method for the gridded Xin'anjiang model based on compact indexing and CUDA graphs according to claim 1, characterized in that, The specific process of parallel calculation of runoff generation and water source distribution and updating of free water state in step S4 is as follows: Step S4.1: In the computational flow, runoff generation calculation is performed using grid cells as parallel computing units, and the soil moisture content status is updated; Step S4.2: Divide the runoff results into surface runoff, interflow and groundwater runoff, and update the free water content status; Step S4.3: Write the water source division results into the intermediate result storage area of ​​the model calculation variables as input for the subsequent step S5.

5. The parallel method for the gridded Xin'anjiang model based on compact indexing and CUDA graphs according to claim 1, characterized in that, The specific process of parallel calculation of slope confluence in step S5 is as follows: Step S5.1: Perform slope runoff calculation in the computation flow using grid cells as parallel computing units, and calculate the slope inflow based on the surface runoff, interflow and groundwater runoff obtained in step S4; Step S5.2: Update the slope confluence status; Step S5.3: Write the slope inflow into the intermediate result storage area of ​​the model calculation variables, as the input for the subsequent step S7 river confluence.

6. The parallel method for the gridded Xin'anjiang model based on compact indexing and CUDA graphs according to claim 1, characterized in that, The specific process of performing parallel computation of river confluence based on the hierarchical compact index table and launching the executable CUDA graph in step S7 is as follows: Step S7.1: Enable the stream capture mechanism during the model initialization phase, according to the hierarchy l = 1, ..., L max The startup sequence of the river confluence kernel function is recorded sequentially, and for each level, the corresponding compact index array IndexList[l] and its length count_l are read, and thread mapping is established according to the index array; Step S7.2: The thread locates the current grid cell based on IndexList[l], reads the slope inflow of the grid cell and the upstream inflow in the channel confluence accumulator, calculates the channel outflow of the current grid cell using the Muskingum method, and accumulates the channel outflow into the channel confluence accumulator of the downstream grid cell using double-precision atomic addition based on the downstream one-dimensional index array. Step S7.3: Complete levels 1 to L max After the kernel function initiates sequence capture, it ends flow capture and instantiates the hierarchical river confluence computation graph into an executable CUDA graph; Step S7.4: During each computation period iteration, update the time step parameters and driver input pointers corresponding to the executable CUDA graph, and start the updated executable CUDA graph in the computation flow to adapt to the double buffering switching of driver data and avoid repeatedly building the execution graph.

Images