Urban inland inundation simulation method based on CPU-GPU heterogeneous parallel computing

Abstract

The invention provides an urban inland inundation simulation method based on CPU-GPU heterogeneous parallel calculation, and aims to reduce the number of calculation grids by adopting a self-adaptive non-uniform grid, consider pipe network calculation and improve the calculation efficiency of low-fluctuation urban area flood simulation at the same time, and comprises the following steps: performing one-dimensional pipe network runoff CPU calculation, constructing a rainwater pipe network model, and constructing a rainwater pipe network model. Extracting an overflow process of an overflow node of the pipeline; dividing a non-uniform grid based on the original raster data and a gradient threshold; non-uniform grid parallel calculation is carried out; a CPU-GPU heterogeneous parallel one-dimensional and two-dimensional coupling model is constructed; and coupling water depth calculation. The method can flexibly adapt to topographic changes and reduce the number of computational grids, and solves the limitation that precision and efficiency are difficult to consider in waterlogging simulation of a uniform grid in a low-fluctuation area in combination with CPU-GPU heterogeneous parallel calculation.

Description

Technical Field

[0001] The present invention belongs to the technical field of urban waterlogging simulation in the field of hydrology and water resources, and specifically relates to an urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing. Background Art

[0002] As an important means of flood prevention and disaster reduction, urban flood simulation aims to predict and evaluate the evolution of urban floods under extreme rainfall conditions, and provide a scientific basis for urban planning, flood control measures design and emergency management.

[0003] Traditional urban flood simulation technology usually divides the surface area into regular grids and uses two-dimensional shallow water wave equations to simulate flood propagation. However, this method has limitations in accuracy and efficiency and is difficult to meet the simulation needs of multi-scale and uneven terrain in complex urban environments.

[0004] In recent years, non-uniform grid technology has attracted attention because it can automatically adjust the grid size according to the terrain, improve simulation accuracy and computational efficiency. Non-uniform grid technology uses denser grids in areas with drastic terrain changes and sparser grids in areas with gentle terrain changes, thereby reducing the amount of calculation while ensuring simulation accuracy. However, most existing non-uniform grid technologies use pre-set grid generation rules, lack the ability to adapt to terrain changes, and are difficult to cope with the rapid changes in urban terrain.

[0005] In addition, the coupling between the drainage network and surface water flow is a key factor in urban flood simulation. However, the traditional simulation method using non-uniform grid technology only considers the flow generation and convergence process of rainwater on the two-dimensional surface, but ignores the water exchange process between the drainage network and the surface.

[0006] In view of the shortcomings of the existing technology, the present invention proposes a method using CPU-GPU heterogeneous parallel computing. This method integrates non-uniform grids and GPU parallel computing while considering the calculation of the pipe network, aiming to more accurately capture the sensitive area of ​​water level calculation while reducing the number of calculation grids, which makes it possible to achieve the difficult balance between accuracy and efficiency in waterlogging simulation in low-relief areas. Summary of the invention

[0007] In view of the problems existing in the above-mentioned technologies, the present invention provides an urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing, which aims to reduce the number of computing grids by adopting an adaptive non-uniform grid, while considering the pipe network calculation, so as to improve the computational efficiency of flood simulation in low-relief urban areas.

[0008] In order to achieve the above-mentioned object of the invention, the technical solution adopted by the present invention is: A method for simulating urban waterlogging based on CPU-GPU heterogeneous parallel computing, comprising: S1: One-dimensional network runoff CPU calculation, building the rainwater network model of the target area, and extracting the overflow process at the overflow node of the pipeline; S2: Adaptive non-uniform grid division; S3: Non-uniform grid GPU parallel computing technology; S4: Construction of one-dimensional coupling model for CPU-GPU heterogeneous parallel computing; S5: Get the water depth file based on the calculation results of S4.

[0009] Furthermore, the S1 comprises the following steps: S11: Collection of drainage network data in the target area, including drainage pipe and drainage node vector data and their field attributes, where drainage nodes include rainwater wells, rainwater grates and water outlets; S12: Divide the catchment area based on the rainwater well. Take the rainwater well as a discrete point and use the ArcGIS software to use Thiessen polygons to complete the catchment area division of the target area. S13: constructing a one-dimensional pipe network runoff model based on the rainwater pipe network data and the catchment area data divided in step S12; S14: Calculate the runoff of the rainwater pipe network on the CPU. The continuity equation and momentum equation, namely the Saint-Venant equation, are used to calculate the flow rate of the pipeline. The specific calculation formula is as follows: , , In the above formula, Represents the distance along the water flow direction, in units of ; Indicates time in units of ; Indicates the cross-sectional area of ​​water flow, in units of ; Indicates flow rate in units of ; Indicates water depth in units of ; represents the resistance slope (head loss per unit length); Represents the acceleration due to gravity in units of , S15: Extract the overflow process at the pipeline overflow node.

[0010] Further, the S2 comprises the following steps: S21: Determine the maximum refinement level based on the original DEM data. The specific calculation is as follows: , In the above formula, is the maximum refinement level, and DEM data and The number of pixels in the direction, S22: Adaptively adjust the original DEM grid resolution to , and assign the newly added grid elevation to -9999; S23: Based on the quadtree block merging algorithm, recursively merge from the finest grid to form a multi-resolution hierarchical quadtree ; S24: Resolution hierarchy quadtree based ,from Layer 1 to Layer 2 in order on each layer Define a size of , the step length is The grid window; S25: Use the defined grid window to scan step by step at the current level, and convert the slope error of the grid window with a user-specified error threshold , perform local mesh merging; S26: Based on the mesh merging result of S25, delete The redundant parent-child node relationship in the grid forms a non-uniform grid hierarchical structure quadtree ; S27: All leaf nodes in the grid are assembled into a non-overlapping grid with different resolutions, that is, a non-uniform grid.

[0011] Furthermore, S3 comprises the following steps: S31: Use Morton code and Z-order curve to establish the index relationship of non-uniform grid elements; Further, the S31 comprises the following steps: S311: Generate Morton code: for each grid element and The binary of the index is bit-interleaved to generate Morton code; S312: Create a Z-order curve: convert the Morton code corresponding to each grid element into an equivalent decimal index and connect them in order from large to small; S313: Index Mapping: Use the Z-order curve to map the indexes of all elements in the multi-resolution hierarchy into a single 1D array and store it in the GPU memory; S314: Generate a multi-resolution hierarchical structure index table and a non-uniform grid index table: by introducing Morton code, construct the connection relationship between different elements of the same refinement level, and store the information of elements in different resolution hierarchical structures and the connection relationship between the elements in the multi-resolution hierarchical structure index table and the non-uniform grid index table respectively; S315: Grid assembly: recursively traverse all leaf nodes in the non-uniform grid hierarchy quadtree from the root node and assemble them into a non-uniform grid composed of non-overlapping grid elements of different resolutions, then map the index of the elements in the non-uniform grid to a single 1D array and store it in the GPU video memory.

[0012] S32: determining the adjacent relationship of non-uniform grid elements; Further, the S32 comprises the following steps: S321: Based on the maximum refinement level (L) and the refinement level (1) where the current non-uniform grid element is located, determine the number of the finest grids contained in the grid element, that is, the number of leaf nodes corresponding to the grid element; S322: Based on multi-resolution hierarchical quadtree , obtain all row and column numbers of the finest grid contained in the current non-uniform grid element from the multi-resolution hierarchical structure index table to determine the row and column numbers of the finest grids around the current non-uniform grid element in the east, west, south and north, and obtain detailed information of these grids from the multi-resolution hierarchical structure index table according to the row and column numbers; S323: Based on multi-resolution hierarchical quadtree , recursively looking up for the ancestor nodes of these most refined grids, which are quadtrees of non-uniform grid hierarchies Nodes in S324: Delete redundant nodes to obtain all adjacent elements of the non-uniform grid.

[0013] S33: Implement flux exchange of non-uniform grid elements. In order to simplify the flux exchange process, the finite difference method is used to store each interface as a single continuous flow component. Based on the water depth of the grid element in the previous time step, the water depth of the grid element in the current time step is updated, and the water depth value stored in the video memory corresponding to the grid element is updated at the same time. The specific formula is as follows: , In the above formula, , Represents elements The water depth at different times, is the calculation step size, , , and are the flows at the eastern, western, southern and northern borders, Is an element area.

[0014] Further, the S4 comprises the following steps: S41: Model boundary condition file configuration. Add the overflow process at the pipeline overflow node extracted in step S1 to the bdy (i.e. node timing overflow) and bci (i.e. node information declaration) files; S42: Model parameter configuration, add bdy and bci files and parameters such as infiltration rate, initial water depth, Manning coefficient, simulation duration, time step, etc. to the config configuration file.

[0015] S43: Surface runoff calculation uses a two-dimensional shallow water equation to describe the urban surface inundation process, taking into account the overall balance of water volume in each surface unit and calculating the temporal and spatial changes of water accumulation in the surface unit. The formula is as follows: , , , In the above formula, Indicates water depth; Indicates time; Respectively represent the plane rectangular coordinate system Direction and Distance of direction; and N express Direction and Single-width flow in one direction; represents the source and sink terms in the continuity control equation, Indicates the surface water level, which is the sum of the surface elevation and the water depth; and The velocity vector in the vertical direction is Direction and Directional component; represents the Manning roughness coefficient, which is determined according to the land type; represents the acceleration due to gravity, S44: Use the GPU parallel computing technology in step S3 to solve the two-dimensional shallow water equations.

[0016] Compared with the prior art, the present invention has the following advantages and beneficial effects: 1. The non-uniform grid method used in the present invention can flexibly adapt to terrain changes and reduce the number of calculation grids. In areas with small terrain fluctuations, this method can reduce redundant calculations of a large number of identical grids, avoid the waste of computing resources, and thus improve computing efficiency; 2. By combining CPU-GPU heterogeneous parallel computing technology, the present invention couples the one-dimensional pipe network runoff model running on the CPU with the two-dimensional surface runoff model integrated with non-uniform grids and GPU parallelism. Compared with the traditional single GPU simulation solution, this heterogeneous parallel computing technology has better efficiency, thereby improving the accuracy and timeliness of urban flooding simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] The drawings described herein are used to provide a further understanding of the embodiments of the present invention, constitute a part of this application, and do not constitute a limitation of the embodiments of the present invention. In the drawings: Figure 1 It is a framework diagram of the method of the present invention; Figure 2 The multi-resolution hierarchical quadtree T1 of the present invention; Figure 3 The non-uniform grid hierarchical structure quadtree T2 of the present invention; Figure 4 It is a schematic diagram of the non-uniform grid division of the present invention; Figure 5 It is the connection relationship between different elements at the same refinement level of the present invention; Figure 6 This is a simulation result of the water depth of rainstorm waterlogging at a certain monitoring station in an embodiment of the present invention. DETAILED DESCRIPTION

[0018] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with embodiments and drawings. The exemplary embodiments of the present invention and their description are only used to explain the present invention and are not intended to limit the present invention.

[0019] Taking a certain area in a city that is prone to flooding caused by heavy rain as an example, the target area has complete drainage network data, DEM data, rainfall data, etc. The total area of ​​the area is about 31km 2 The terrain roughness is 2.57, the area with a slope of more than 10 degrees accounts for about 0.3% of the total area, the area with a slope of more than 30 degrees accounts for about 0.008% of the total area, the average slope is 0.88, and the terrain undulation is relatively small.

[0020] like Figure 1 As shown, an embodiment of the present invention discloses a method for simulating urban waterlogging based on CPU-GPU heterogeneous parallel computing, which mainly includes the following steps: S1: One-dimensional pipe network runoff CPU calculation. Construct a rainwater pipe network model for the target area and use the one-dimensional pipe network runoff model to extract the overflow process at the pipe overflow node; S2: Input the terrain data (DEM) of the target area, and use the adaptive non-uniform grid division algorithm to automatically adjust the grid resolution of the terrain data according to the terrain undulation; S3: Non-uniform grid GPU parallel computing. Non-uniform grids include elements of non-uniform grid resolution, so each element can have one coarser neighbor or multiple finer neighbors. In order to enable non-uniform grids to perform GPU parallel computing, the index relationship between non-uniform grid elements must first be established to determine the neighbor relationship between elements in the non-uniform grid and perform flux exchange. On this basis, the kernel program is started to dispatch GPU threads for each grid element for parallel computing; S4: Construction of a one- and two-dimensional coupling model for CPU-GPU heterogeneous parallel computing. Construct a one- and two-dimensional hydrological and hydrodynamic coupling model for the target area, using the non-uniform grid divided in S2 as the calculation grid, the overflow process of the pipeline node as the two-dimensional surface water input, and the two-dimensional surface runoff model to model the urban surface runoff process. The GPU parallel computing technology in step S3 is used for flux exchange calculation.

[0021] S5: Get the water depth file based on the calculation results of S4.

[0022] Furthermore, the S1 comprises the following steps: S11: Drainage network data collection in the target area. The underground drainage network in the target area includes drainage pipes and drainage node data, where drainage nodes include rainwater wells, rainwater grates and outlets. The attribute information of drainage pipes includes pipe name, length, material, cross-sectional geometry, inlet (outlet) water node name and elevation. The attribute information of drainage pipes includes name, surface elevation, and bottom elevation; S12: Drainage area division based on rainwater wells. Rainwater wells are used as discrete points, and the target area is divided using Thiessen polygons with the help of ArcGIS software; S13: Construction of one-dimensional pipe network runoff model. Based on the drainage pipe network data and the catchment area data divided in S12, a one-dimensional pipe network runoff model is established; S14: Calculate the runoff of the rainwater pipe network on the CPU. The continuity equation and momentum equation, namely the Saint-Venant equation, are used to calculate the flow rate of the pipeline. The specific calculation formula is as follows: , , In the above formula, Represents the distance along the water flow direction, in units of ; Indicates time in units of ; Indicates the cross-sectional area of ​​water flow, in units of ; Indicates flow rate in units of ; Indicates water depth in units of ; represents the resistance slope (head loss per unit length); Represents the acceleration due to gravity in units of , S15: Extract the overflow process at the pipeline overflow node.

[0023] Further, the S2 comprises the following steps: S21: Determine the maximum refinement level based on the original DEM data. The specific calculation is as follows: , In the above formula, is the maximum refinement level, and DEM data and The number of pixels in the direction, S22: Adaptively adjust the original DEM grid resolution to , and assign the newly added grid elevation to -9999; S23: Based on the quadtree block merging algorithm, recursively merge from the finest grid to form a multi-resolution hierarchical quadtree , Figure 2 Shows a The multi-resolution hierarchy and the corresponding multi-resolution hierarchy quadtree ; S24: Resolution hierarchy quadtree based ,from Layer 1 to Layer 2 in order on each layer Define a size of , the step length is The grid window; S25: Use the defined grid window to scan step by step at the current level to determine whether there is a grid element with an elevation value of -9999 in the current window. If yes, the grid window moves forward without processing. If no, compare the slope error of the current grid window with a user-specified error threshold ( ), and perform local mesh merging.

[0024] when When updating The elevation value of the grid element of the parent and child nodes of the layer is the average elevation value of the four grids in the grid window. At the same time, the current Layer four grid elements with corresponding The parent-child node relationship of the layer grid elements is used to realize local grid merging; otherwise, no operation is performed and the slope error The calculation formula is as follows: , , , , , In the above formula, For grid elements exist The slope in the direction, For grid elements exist The slope in the direction, For grid elements The average slope, is the maximum refinement level, is the current refinement level, is the resolution of the finest grid, is the elevation value of the grid element. , , , The grid elements , , , The average slope, is the average slope of the current grid window, The slope error of the current grid window.

[0025] S26: Based on the mesh merging result of S25, delete The redundant parent-child node relationship in the grid generates a new quadtree, called a non-uniform grid hierarchical quadtree. ,like Figure 3 As shown; S27: Quadtree from non-uniform grid hierarchy The root node recursively traverses all leaf nodes and assembles all leaf nodes into a non-overlapping grid of different resolutions, that is, a non-uniform grid, such as Figure 4 As shown; Furthermore, S3 comprises the following steps: S31: Establishing index relationships between non-uniform grid elements. In order to achieve efficient representation and operation on complex grids, it is crucial to establish index relationships between elements. One of the commonly used methods is to apply Morton codes and Z-order curves. This bit interleaving technology can cleverly map multi-dimensional space coordinates to one-dimensional space.

[0026] Further, the S31 comprises the following steps: S311: Generate Morton code: Bit-interleave the binary of the X and Y indexes of each grid element to generate a Morton code; S312: Create a Z-order curve: convert the Morton code corresponding to each grid element into an equivalent decimal index and connect them in order from large to small; S313: Index Mapping: Use the Z-order curve to map the indexes of all elements in the multi-resolution hierarchy into a single 1D array and store it in the GPU memory. S314: Generate a multi-resolution hierarchical structure index table and a non-uniform grid index table: By introducing Morton code, the connection relationship between different elements at the same refinement level is constructed, and the information of elements in different resolution hierarchical structures and the connection relationship between elements are stored in the multi-resolution hierarchical structure index table and the non-uniform grid index table respectively. S315: Grid assembly: recursively traverse all leaf nodes in the non-uniform grid hierarchy quadtree from the root node, assemble all leaf nodes into a non-uniform grid composed of non-overlapping grid elements of different resolutions, and map the index of the elements in the non-uniform grid to a single 1D array and store it in the GPU memory, such as Figure 5 As shown; S32: Determine the neighbor relationship of non-uniform grid elements. Considering the unique data structure of non-uniform grids, each element may have one coarser neighbor or multiple finer neighbors around it. In order to achieve flux exchange between grids, it is particularly important to clarify the neighbor relationship between elements. In non-uniform grids, the neighbor relationship of elements is usually not intuitively visible, so a series of steps are required to accurately determine the neighbor relationship between elements.

[0027] Further, the S32 comprises the following steps: S321: Based on the maximum refinement level ( L ) and the refinement level of the current non-uniform mesh element ( l ), determine the number of the finest grids contained in the grid element, that is, the number of leaf nodes corresponding to the grid element; S322: Based on the multi-resolution hierarchical quadtree T1, all row and column numbers of the finest grid contained in the current non-uniform grid element are obtained from the multi-resolution hierarchical index table to determine the row and column numbers of the finest grids around the east, west, south and north of the current non-uniform grid element, and detailed information of these grids is obtained from the multi-resolution hierarchical index table according to the row and column numbers; S323: Based on multi-resolution hierarchical quadtree , recursively search for the ancestor nodes of these finest grids, which must be non-uniform grid hierarchy quadtrees Nodes in .

[0028] S324: Delete redundant nodes to obtain all adjacent elements of the non-uniform grid.

[0029] S33: Implement flux exchange of non-uniform grid elements. In order to simplify the flux exchange process, the finite difference method is used to store each interface as a single continuous flow component. Based on the water depth of the grid element in the previous time step, the water depth of the grid element in the current time step is updated, and the water depth value stored in the video memory corresponding to the grid element is updated at the same time. The specific formula is as follows: , In the above formula, , Represents elements The water depth at different times, is the calculation step size, , , and are the flows at the eastern, western, southern and northern borders, Is an element area.

[0030] Further, the S4 comprises the following steps: S41: Model boundary condition file configuration. Add the overflow process at the pipeline overflow node extracted in step S1 to the bdy (i.e. node timing overflow) and bci (i.e. node information declaration) files; S42: Model parameter configuration. Add bdy and bci files and parameters such as infiltration rate, initial water depth, Manning coefficient, simulation duration, time step, etc. to the config configuration file; S43: Surface runoff calculation uses a two-dimensional shallow water equation to describe the urban surface inundation process, taking into account the overall balance of water volume in each surface unit and calculating the temporal and spatial changes of water accumulation in the surface unit. The formula is as follows: , , , In the above formula, Indicates water depth; Indicates time; Respectively represent the plane rectangular coordinate system Direction and Distance of direction; and N express Direction and Single-width flow in one direction; represents the source and sink terms in the continuity control equation, Indicates the surface water level, which is the sum of the surface elevation and the water depth; and The velocity vector in the vertical direction is Direction and Directional component; represents the Manning roughness coefficient, which is determined according to the land type; represents the acceleration due to gravity; S44: Use the GPU parallel computing technology in step S3 to solve the two-dimensional shallow water equations.

[0031] Figure 6 This is an example of the simulation result of the water depth of a rainstorm at a monitoring station in this embodiment. It can be seen from the figure that the results obtained by the simulation method using CPU-GPU parallel computing are roughly consistent with the results obtained by the traditional simulation method using only CPU computing.

[0032] The above embodiments are only used to help understand the method and core idea of ​​the present invention; at the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. The urban flooding simulation method based on CPU-GPU heterogeneous parallel computing is characterized by: The method comprises the following steps: S1: One-dimensional network runoff CPU calculation, building the rainwater network model of the target area, and extracting the overflow process at the overflow node of the pipeline; S2: Adaptive non-uniform grid division; S3: Non-uniform grid GPU parallel computing; S4: Construction of one-dimensional coupling model for CPU-GPU heterogeneous parallel computing; S5: Get the water depth file based on the calculation results of S4.

2. The urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing according to claim 1 is characterized in that: The S1 comprises the following steps: S11: Collection of drainage network data in the target area, including drainage pipe and drainage node vector data and their field attributes, where drainage nodes include rainwater wells, rainwater grates and water outlets; S12: Divide the catchment area based on the rainwater well. Take the rainwater well as a discrete point and use the ArcGIS software to use Thiessen polygons to complete the catchment area division of the target area. S13: constructing a one-dimensional pipe network runoff model based on the rainwater pipe network data and the catchment area data divided in step S12; S14: Calculate the runoff of the rainwater pipe network on the CPU. The continuity equation and momentum equation, namely the Saint-Venant equation, are used to calculate the flow rate of the pipeline. The specific calculation formula is as follows: ， ， In the above formula, Represents the distance along the water flow direction, in units of ; Indicates time in units of ; Indicates the cross-sectional area of ​​water flow, in units of ; Indicates flow rate in units of ; Indicates water depth in units of ; represents the resistance slope; Represents the acceleration due to gravity in units of , S15: Extract the overflow process at the pipeline overflow node.

3. The urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing according to claim 1 is characterized in that: The S2 comprises the following steps: S21: Determine the maximum refinement level based on the original DEM data. The specific calculation is as follows: ， In the above formula, is the maximum refinement level, and DEM data and The number of pixels in the direction, S22: Adaptively adjust the original DEM grid resolution to , and assign the newly added grid elevation to -9999; S23: Based on the quadtree block merging algorithm, recursively merge from the finest grid to form a multi-resolution hierarchical quadtree ; S24: Resolution hierarchy quadtree based ,from Layer 1 to Layer 2 in order on each layer Define a size of , the step length is The grid window; S25: Use the defined grid window to scan step by step at the current level, and convert the slope error of the grid window with a user-specified error threshold , perform local mesh merging; S26: Based on the mesh merging result of S25, delete The redundant parent-child node relationship in the grid forms a non-uniform grid hierarchical structure quadtree ; S27: All leaf nodes in the grid are assembled into a non-overlapping grid with different resolutions, that is, a non-uniform grid.

4. The urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing according to claim 1 is characterized in that: S3 adaptive non-uniform grid GPU parallel computing includes the following steps: S31: Morton code and Z-order curve are used to establish the index relationship of non-uniform grid elements, and the index of the elements in the non-uniform grid is mapped to a single Array and stored in GPU memory; S32: determining the adjacent relationship of non-uniform grid elements; S33: Implement the flux exchange of non-uniform grid elements. Use the finite difference method to store each interface as a single continuous flow component. Update the water depth of the grid element in the current time step based on the water depth of the grid element in the previous time step. At the same time, update the water depth value stored in the video memory corresponding to the grid element. The specific formula is as follows: ， In the above formula, , Represents elements The water depth at different times, is the calculation step size, , , and are the flows at the eastern, western, southern and northern borders, Is an element area.

5. The urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing according to claim 4 is characterized in that: S31 uses Morton code and Z-order curve to establish the index relationship of non-uniform grid elements, including the following steps: S311: Generate Morton code: for each grid element and The binary of the index is bit-interleaved to generate Morton code; S312: Create a Z-order curve: convert the Morton code corresponding to each grid element into an equivalent decimal index, and connect them in order from large to small; S313: Index Mapping: Use the Z-order curve to map the indexes of all elements in the multi-resolution hierarchy into a single 1D array and store it in the GPU memory; S314: Generate a multi-resolution hierarchical structure index table and a non-uniform grid index table: by introducing Morton code, construct the connection relationship between different elements of the same refinement level, and store the information of elements in different resolution hierarchical structures and the connection relationship between the elements in the multi-resolution hierarchical structure index table and the non-uniform grid index table respectively; S315: Grid assembly: recursively traverse all leaf nodes in the non-uniform grid hierarchy quadtree from the root node and assemble them into a non-uniform grid composed of non-overlapping grid elements of different resolutions, then map the index of the elements in the non-uniform grid to a single 1D array and store it in the GPU video memory.

6. The urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing according to claim 4 is characterized in that: S32 determines the adjacent relationship of non-uniform grid elements, including the following steps: S321: Based on the maximum refinement level and the refinement level where the current non-uniform grid element is located, determine the number of the finest grids contained in the grid element, that is, the number of leaf nodes corresponding to the grid element; S322: Based on multi-resolution hierarchical quadtree , obtain all row and column numbers of the finest grid contained in the current non-uniform grid element from the multi-resolution hierarchical structure index table to determine the row and column numbers of the finest grids around the current non-uniform grid element in the east, west, south and north, and obtain detailed information of these grids from the multi-resolution hierarchical structure index table according to the row and column numbers; S323: Based on multi-resolution hierarchical quadtree , recursively looking up for the ancestor nodes of these most refined grids, which are quadtrees of non-uniform grid hierarchies Nodes in S324: Delete redundant nodes to obtain all adjacent elements of the non-uniform grid.

7. The urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing according to claim 1 is characterized in that: S4 includes the following steps: S41: Model boundary condition file configuration, adding the overflow process at the pipeline overflow node extracted in step S1 to the bdy, i.e., node timing overflow volume, and bci, i.e., node information declaration file; S42: Model parameter configuration, add the bdy and bci files and the infiltration rate, initial water depth, Manning coefficient, simulation duration, and time step parameters to the config configuration file. S43: Surface runoff calculation uses a two-dimensional shallow water equation to describe the urban surface inundation process, taking into account the overall balance of water volume in each surface unit and calculating the temporal and spatial changes of water accumulation in the surface unit. The formula is as follows: ， ， ， In the above formula, Indicates water depth; Indicates time; Respectively represent the plane rectangular coordinate system Direction and Distance of direction; and N express Direction and Single-width flow in one direction; represents the source and sink terms in the continuity control equation, Indicates the surface water level, which is the sum of the surface elevation and the water depth; and The velocity vector in the vertical direction is Direction and Directional component; represents the Manning roughness coefficient, which is determined according to the land type; represents the acceleration due to gravity, S44: Use the GPU parallel computing technology in step S3 to solve the two-dimensional shallow water equations.

8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing as described in any one of claims 1 to 7 above is implemented.

9. A computer-readable storage medium having computer instructions stored thereon, characterized in that: When the computer instructions are executed by the processor, the urban waterlogging simulation method based on CPU-GPU heterogeneous parallel computing as described in any one of claims 1 to 7 is implemented.

Images