Method for fast extraction of dfsu file point data

Abstract

The application provides a method for quickly extracting dfsu file point data, relates to the technical field of numerical model calculation, and solves the problems of low efficiency, long time consumption and limited batch automation of MIKE software dfsu file extraction of specific spatial point data. The scheme points include: converting the dfsu result file and the grid mesh file into the mat format readable by MATLAB, pre-storing the global grid topology relationship; inputting the sampling point coordinates and locating the grid m and the position thereof; dynamically generating the adjacent grid set {Cm} of m, calculating the physical quantity based on {Cm} through a boundary perception interpolation method, wherein if the sampling point is a boundary point, an extended grid set {Cm'} is created, that is, {Cm'}={Cm}∪{Me}, and interpolation is carried out based on the boundary point {Cm'}. Multi-time step operation and multi-point cycle extraction are supported. The method is mainly used for efficient analysis of large-scale simulation results of ocean and water conservancy engineering.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of numerical model calculation. More particularly, the present application relates to a method for quickly extracting point data of a dfsu file. BACKGROUND

[0002] MIKE21 software is a professional hydrodynamic calculation software developed by DHI. The hydrodynamic module is based on two / three-dimensional incompressible Reynolds-averaged Navier-Stokes equations, and is solved by finite volume method (FVM). One of the core features of the finite volume method is that physical quantities (such as flow velocity, water level, concentration, etc.) are stored in the center or face of the calculation unit (grid), while the grid nodes are mainly used to define the geometric shape. Although this storage method has advantages in numerical solving equations, such as good conservation, it brings significant inconvenience to users to extract result data at specific spatial positions (points) subsequently.

[0003] Although MIKE21 software itself provides a user-friendly graphical post-processing module, it is easy to view and export small-scale data for routine results, but as a commercial closed software, it has obvious shortcomings when batch processing and automatic processing of massive calculation result data, mainly in:

[0004] 1. Low data extraction efficiency: For a dfsu result file containing 100,000 grid units and a time span of one year, even on a desktop computer with high configuration, such as 16GB of memory, 12 CPU threads, and a main frequency of 3.10GHz, it takes minutes to extract data of a single spatial point over the entire time series. For research or engineering applications (such as model validation, parameter calibration, data assimilation, report generation) that need to process hundreds of sampling points, the time cost is too high, which becomes a bottleneck in the workflow.

[0005] 2. Limited automatic processing: When performing large-scale, customized data extraction and subsequent analysis (such as statistical analysis, spatial interpolation, and specific format output), users often need to export the dfsu file data to an intermediate format (such as text, CSV, or NetCDF) for secondary processing. This process is not only tedious and increases the risk of errors, but also limited by the export function and interface flexibility provided by the software.

[0006] The reasons for the above-mentioned deficiencies are as follows: 1, the closedness of commercial software: the core calculation engine (the underlying language is Fortran or C / C++) of MIKE21 and its data storage format (dfsu) are closed. Users cannot directly optimize the algorithm logic of internal data reading and processing. The post-processing module provided by the official mainly faces interactive operation, and the internal implementation may not be optimal for batch point data extraction. 2, the limitation of the underlying language processing mode: Fortran and C / C++ traditionally rely on explicit loops when processing large arrays. Although these languages have high performance, the loop traverses a large number of grid cells to locate the sampling points and perform interpolation calculation, which is difficult to fully utilize the parallel computing capability (such as multi-core, vector instruction set) under the modern CPU architecture, and becomes the main time-consuming link. 3, the limitation of the official MATLAB toolkit: DHI provides a MATLAB toolkit for reading dfsu format files. Although this provides a way to access data in the MATLAB environment, it is essentially only a wrapper for the underlying library. When reading files, the toolkit usually needs to load all the data in the target time period and the entire calculation domain (or selected area) into the memory, and then the user still needs to write loop logic to traverse the grid and perform interpolation in the memory array. Therefore, the toolkit mainly solves the problem of "being able to read", and does not fundamentally solve the core efficiency bottlenecks of "reading fast" and "batch point extraction fast". The performance improvement is not significant compared to directly using the MIKE post-processing module to export data, especially when processing large files. SUMMARY

[0007] The present application provides a method for quickly extracting point data of a dfsu file, which can significantly improve the positioning and calculation efficiency by pre-storing grid topology relationship, dynamically constructing adjacent grid set and boundary-aware interpolation. The problems of low efficiency of point data extraction of the dfsu file of the MIKE software and limitation of batch automation are solved.

[0008] To achieve these objects and other advantages and in view of its purposes, the present application provides a method for quickly extracting point data of a dfsu file, comprising the following steps:

[0009] S1, converting the result file in dfsu format generated by the MIKE software into a mat format file readable by MATLAB;

[0010] S2, converting the mesh file input by the MIKE software into a mat format file, pre-storing the global grid topology relationship, including: each grid number and its corresponding three node numbers, the coordinates of all nodes, the grid number set { Mi} associated with each node, and boundary marker information, wherein the boundary marker information includes boundary node identification and boundary edge identification;

[0011] S3, input the coordinates of the sampling point P, calculate the grid m where it is located and its position in the grid based on the pre-stored global grid topology relationship;

[0012] S4, dynamically generate the adjacent grid number set { Cm} of the grid m according to the node relationship of the grid m, and calculate the physical quantity data of the sampling point P based on the grid data in { Cm} through the boundary-aware interpolation method;

[0013] S5, if multiple time step data are extracted, perform array operation or loop execution on step S4 within the time step;

[0014] S6, if multiple sampling point data are extracted, loop execution is performed on steps S3 to S5, and the extraction result is output;

[0015] The boundary-aware interpolation method includes: a, if P is a boundary point, create an extended grid set { Cm'} = { Cm} ∪ { Me}; { Me} performs the following steps in order: 1) preferentially obtain all adjacent grids of the grid m sharing non-boundary edges, denoted as set {Me1}; 2) if {Me1} is an empty set, that is, the grid m has no non-boundary edge adjacent grid, then all internal grids within a radius R centered at P are obtained and denoted as set {Me2}, wherein R is 2 times the average grid size; 3) finally {Me} = {Me1} ∪ {Me2}; b, interpolation is performed based on the boundary point of { Cm'}.

[0016] Preferably, the grid m where the sampling point P is located and its position in the grid are calculated in step S3, and the specific process is as follows:

[0017] S301, according to the coordinates (x p , y p ) of the sampling point P, screen the grid set { Mp}, and the screening condition satisfies: min(X) ≤x p ≤ max(X) and min(Y) ≤y p ≤ max(Y); wherein X=(X1, X2, X3) and Y=(Y1, Y2, Y3) are the coordinates of the three vertices of the grid;

[0018] S302, traverse all grids in { Mp}, and determine the position relationship through the MATLAB command [in, on] = inpolygon(x p , y p , X, Y);

[0019] If in=0, it indicates that it is outside the grid;

[0020] If in=1 and on=0, it indicates that it is inside the grid;

[0021] in = 1 and on≥1 represent at the grid edge or node, the on value accumulates the shared edge / node number, wherein on takes values of 0, 1, 2, 3;

[0022] S303, determine the grid m and the position category according to the in and on values;

[0023] S304, if the sampling point P is located at the boundary edge or the boundary node, mark it as a boundary point.

[0024] Preferably, the interpolation method in step S4 includes any one of the following:

[0025] S401, adopt the nearest grid interpolation, and divide into three cases:

[0026] A, if the sampling point P is inside the grid m or at the non-shared edge / node, directly take H P = H m , wherein H P is the physical quantity data value of the sampling point P, and H m is the physical quantity data value of the center point of the grid m;

[0027] B, if the sampling point P is at the shared edge of two grids, linearly interpolate the associated node data according to the distance weight, specifically:

[0028] H P = [d(P, N2)H N1 + d(P, N1)H N2 ] / [d(P, N1) + d(P, N2)], wherein N1 and N2 are two nodes of the shared edge of the two grids, d(P, N1) is the distance from the sampling point P to the node N1, d(P, N2) is the distance from the sampling point P to the node N2, H N1 is the physical quantity data value of the node N1, and H N2 is the physical quantity data value of the node N2;

[0029] C, if the sampling point P coincides with the node N, take H P = H N , wherein H N is the physical quantity data value of the node N;

[0030] wherein the physical quantity data value H N of the node N is calculated by inversely distance weighted interpolation of the grid set { Mn}: , k is the number of grids in the grid set { Mn}, is the grid center physical quantity data value, is the distance between the node N and the grid Mn j center point, is the distance between the node N and the grid Mn iDistance of the center point

[0031] S402, using node value interpolation, the grid m three nodes where the sampling point P is respectively: m1 m2 m3 , first calculate the grid m three node data H Nm1 , H Nm2 , H Nm3 , and then through the barycentric coordinate method or inverse distance weighted method interpolation H P ;

[0032] S403, call MATLAB command griddata or scatteredInterpolant, to get H P ;

[0033] Wherein, the physical quantity includes wave height, wave direction, wave period, flow direction, flow velocity, water level and concentration.

[0034] Preferably, the barycentric coordinate method or inverse distance weighted method in step S402 is specifically:

[0035] a), barycentric coordinate method:

[0036] , , ;

[0037] Wherein, The area of the triangle surrounded by three nodes N m1 , N m2 , N m3 , The area of the triangle surrounded by sampling point P, node N m2 , node N m3 , m1 The area of the triangle surrounded by node N m1 , node P, node N m3 , The barycentric coordinate weight of sampling point P relative to grid node N m1 , The barycentric coordinate weight of sampling point P relative to grid node N m2 , The barycentric coordinate weight of sampling point P relative to grid node N m3 ;

[0038] b), inverse distance weighted method:

[0039] , i=1, 2, 3, d(P, N mi ) is the distance between sampling point P and node N mithe distance between the sampling point P and the node N mj , d(P, N mj ) is the distance between the sampling point P and the node N , w(P, N mi ) is the weight of the sampling point P to the node N

[0040] For the barycentric coordinate method and the inverse distance weighting method, the calculation of the physical quantity data value H p of the sampling point P is as follows: H p = ×H Nm1 + ×H Nm2 + ×H Nm3 .

[0041] Preferably, before the MATLAB interpolation command is called in step S403, the following steps are performed:

[0042] 1) Create a memory pool to cache the center point coordinates and physical quantity data of the grids in the grid number set { Cm}, and the initial memory pool threshold is set to 2 GB;

[0043] 2) Add a system memory monitoring module: real-time detection of the available physical memory of the MATLAB process:

[0044] 3) Dynamic adjustment method: a, if the memory pool occupancy > the current memory pool threshold, trigger LRU release, delete the cache data that has not been accessed within the last 15 minutes, and repeat the deletion until the memory pool occupancy ≤ 90% of the current memory pool threshold; b, i) if the system available memory is less than 4 GB, suspend the memory pool function, and instead read the mat file directly, and skip the subsequent downgrade operation; ii) if the system available memory is greater than or equal to 4 GB and less than 8 GB, reduce the performance threshold N th from 5000 to 2000; compress the memory pool threshold from 2 GB to 1 GB; c) when the system available memory ≥ 8 GB for 5 minutes, restart the memory pool function, restore the memory value threshold to 2 GB, and restore the performance threshold N th to 5000.

[0045] Preferably, in step S303, when on ≥ 1, the grid set { Ms} of the shared edge / node is recorded synchronously, and the specific process is as follows:

[0046] a) When the sampling point P is located on the grid edge, according to the two node numbers N1 and N2 of the grid edge, extract the grid containing N1 or N2 from the pre-stored node-related grid set { Mi} in step S2, take the union to generate { Ms};

[0047] b) When the sampling point P coincides with the node N, directly call step S2 to prestore the node-related grid set {Mn} as {Ms};

[0048] In step S4, based on {Ms}, dynamically generate the adjacent grid set {Cm} and then perform the boundary-aware interpolation.

[0049] Preferably, a verification mechanism is added to the nearest grid interpolation in step S401: when the sampling point P is located on the grid edge or node, i.e., on≥1, the following operations are performed:

[0050] a) Synchronously use step B in step S401 and step S402 to respectively calculate the interpolation results H P1 and H P2 ;

[0051] b) Calculate the relative deviation δ = |H P1 ﹣H P2 | / max(|H P1 |, |H P2 |, ε), where ε is a zero minimum prevention quantity;

[0052] c) If δ>η, where η is a preset tolerance threshold, trigger the grid topology review:

[0053] c1) Recheck the grid m and the shared edge / node relationship determined in step S303;

[0054] c2) Review the integrity of the grid sets {M N1},{M N2} associated with the nodes N1, N2 in step S2;

[0055] d) Re-execute step B or step S402 in step S401 based on the reviewed topology structure, and output the final value with a review mark.

[0056] Preferably, in the inverse distance weighted method calculation in step S402, a distance threshold judgment mechanism is added:

[0057] a) Calculate the distances d1, d2, d3 of the sampling point P to the three nodes (N m1 , N m2 , N m3 ) of the grid m;

[0058] b) If there exists d k <β, k∈{1, 2, 3}, β is a minimum positive number related to the floating-point calculation precision, and is taken as 10 -6 -10 -10 , then it is determined that P coincides with the node N mk ;

[0059] c) Directly assign H P = H Nmk , and terminate the subsequent weighted calculation process.

[0060] d) Otherwise, calculate the inverse distance weighted value by the inverse distance weighted method of step S402.

[0061] Preferably, a sampling point pre-screening module is added before step S3, which specifically performs:

[0062] a) Based on all node coordinates stored in step S2, construct a global outer envelope rectangle R:

[0063] R = { (x, y) | x ∈ [min(X all )-Ψ, max(X all )+Ψ], y ∈ [min(Y all )-Ψ, max(Y all )+Ψ]}.

[0064] Where X all is the set of all node x coordinates, Y all is the set of all node y coordinates, and Ψ is the safety margin.

[0065] b) Create a spatial index structure: evenly divide the outer envelope rectangle R into Q×L grid cells, and record the grid number set covered by each cell.

[0066] c) For each sampling point P(x p , y p ):

[0067] c1) If (x p , y p ) R, directly mark it as an invalid point and skip the subsequent processing.

[0068] c2) If (x p , y p ) ∈ R, locate the corresponding grid cell according to the spatial index.

[0069] c3) Only load the grid set associated with this grid cell to perform grid traversal in step S3.

[0070] Preferably, when multiple sampling point data are extracted in step S6, a parallel processing procedure is performed:

[0071] S601, according to the pre-constructed spatial index structure, divide the sampling point set into blocks according to the spatial grid cells.

[0072] S602, call parallel computing units to synchronously perform grid positioning and physical quantity interpolation operations on the sampling points in each block.

[0073] S603, adopt a pipeline loading mechanism for the multi-time step data, and read next time step data asynchronously;

[0074] S604, aggregate each sub-block interpolation result and output

[0075] The present application at least includes the following beneficial effects:

[0076] First, by pre-storing the global grid topology relationship, including: each grid number and its corresponding three node numbers, the coordinates of all nodes, the grid number set {Mi} associated with each node, boundary marker information, the original dfsu file is converted into an efficient mat format, avoiding time-consuming real-time parsing. The boundary-aware interpolation dynamically expands the adjacent grid set ({Cm'}={Cm}∪{Me}), effectively solving the interpolation distortion caused by incomplete data of boundary points. Combined with multi-time step array operation and multi-point loop, batch automatic extraction is realized, and the time consumption of single-point extraction of 100,000 grids is reduced from minutes to seconds. Based on the outer envelope rectangle, the candidate grid set {Mp} is pre-screened, and the traversal range is reduced to a local area; the inpolygon function of MATLAB is used to accurately identify the position relationship (in / on value) of the point and the grid, and the boundary points are marked synchronously. Compared with global traversal, the grid search amount is reduced by more than 70%, and the boundary point identification provides key input for subsequent expansion of the interpolation set {Cm'}, avoiding redundant calculation. Three adaptive interpolation strategies are provided: nearest grid interpolation directly reads the grid center value or weighted edge node value, suitable for internal points and simple shared edge scenarios, with the highest efficiency; node value interpolation calculates the triangular point value through the barycentric method / inverse distance method, with higher accuracy; MATLAB built-in function supports complex adjacent set interpolation. Users can flexibly choose according to accuracy requirements, and the average interpolation error is controlled within 0.5%. The weight calculation formula of the barycentric coordinate method (area weight) and the inverse distance weighting method (distance reciprocal weight) is clearly defined, ensuring that both output H p The barycentric method has clear geometric meaning and is suitable for uniform grids; the inverse distance method is robust for non-uniform grids, both of which realize efficient programming through explicit formulas, with calculation time less than 1ms.

[0077] Second, the dynamic memory pool and LRU release strategy compresses the memory occupation to the local data of the adjacent grid set {Cm}. Incremental loading avoids repeated I / O; dynamic function selection (griddata / scatteredInterpolant) switches to efficient object interpolation when the number of grids exceeds the threshold, and automatically reduces N th , ensuring that an 8GB memory device can handle 50,000+ grid sets.

[0078] Synchronous generation of related grid set {Ms} when positioning the grid, all grids sharing edges / nodes, directly input interpolation module as {Cm}. Save the overhead of dynamic retrieval of adjacent grids, especially for complex boundary points, the generation speed of adjacent set is improved by 50%.

[0079] Double interpolation result cross verification (H P1 / H P2 ) and out-of-tolerance review mechanism, effectively identify topology errors such as missing grid association. Trigger topology verification when δ>η, ensure that the output value is marked with a review mark, reliability is more than 99%, suitable for high-precision scenes such as hydrological model verification. Distance threshold β avoids the risk of division by zero caused by floating-point precision. When the distance between the point and the node is <β, directly assign the node value, skip the inverse distance calculation process, eliminate numerical instability, and improve the robustness of the algorithm.

[0080] Third, the spatial index structure is divided into grid cells by the outer envelope rectangle R, which quickly excludes invalid points outside the domain. Local loading of unit related grid set compresses the grid traversal range from global to 1-2 cells, and the search efficiency is improved by 90%. Parallelization of the block and asynchronous pipeline fully utilizes multi-core CPU: sampling points are distributed to parallel threads according to spatial cells, and time step data preloading hides I / O delay. The actual time consumption of point extraction is reduced from hours to minutes, and the speedup ratio is close to the core number linear growth.

[0081] Other advantages, objects, and features of the present application will be apparent to those skilled in the art from the following description, and will be understood to be within the scope of the present application. BRIEF DESCRIPTION OF DRAWINGS

[0082] Figure 1 Flowchart of the method for quickly extracting dfsu file point data of the present application;

[0083] Figure 2 Grid relationship diagram of a specific example of the present application;

[0084] Figure 3 P-point physical quantity data obtained by the interpolation method. DETAILED DESCRIPTION

[0085] The present application will be further described in detail below, so that those skilled in the art can implement it according to the description.

[0086] It should be understood that the terms such as "have", "contain" and "include" used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0087] As Figures 1-3 shown, the present application provides a method for quickly extracting dfsu file point data, comprising the following steps:

[0088] S1, convert the result file in dfsu format generated by MIKE software into a mat format file readable by MATLAB;

[0089] S2, convert the mesh file input by MIKE software into a mat format file, prestore the global grid topology relationship, including: each grid number and its corresponding three node numbers, the coordinates of all nodes, the grid number set {Mi} associated with each node, and the boundary marker information, wherein the boundary marker information includes boundary node identification and boundary edge identification;

[0090] S3, input the coordinates of the sampling point P, and calculate the grid m where it is located and its position in the grid based on the pre-stored global grid topology relationship;

[0091] S4, dynamically generate the adjacent grid number set {Cm} according to the node relationship of the grid m, and calculate the physical quantity data of the sampling point P based on the grid data in {Cm} through the boundary-aware interpolation method, the physical quantity including but not limited to wave height, wave direction, wave period, flow direction, flow velocity, water level and concentration, etc.;

[0092] S5, if multiple time step data are extracted, array operation is performed within the time step or step S4 is executed in a loop;

[0093] S6, if multiple sampling point data are extracted, steps S3 to S5 are executed in a loop, and the extraction result is output;

[0094] The boundary-aware interpolation method includes: a, if P is a boundary point, create an extended grid set {Cm'} = {Cm} ∪ {Me}; {Me} performs the following steps in order: 1) preferentially acquire all adjacent grids of grid m sharing non-boundary edges, denoted as set {Me1}; 2) if {Me1} is an empty set, i.e. grid m has no non-boundary edge adjacent grid, then append to acquire all internal grids within a radius R centered on P, denoted as set {Me2}, wherein R is 2 times the average grid size; 3) finally {Me} = {Me1} ∪ {Me2}; b, interpolate based on {Cm'} as the boundary point basis.

[0095] In the above embodiment, the dfsu result file generated by MIKE is converted into a MATLAB readable mat format. The dfs2mat function in the DHI official MATLAB toolkit can be used to perform the conversion, and the physical quantity data is stored as a three-dimensional array by time step. The mesh file is parsed into mat format by the readDfsuMesh function, pre-stored with global topological relationships: a mapping table between mesh numbers and three-node numbers, such as mesh ID: 1001 and [N201, N305, N417]; a node coordinate array, such as N201: [x=120.5, y=30.2]; a node-associated mesh set {Mi}, such as N201: [1001, 1002, 1005]; boundary identifiers, with boundary nodes marked as 1 and boundary edges marked as "open", etc. The boundary marker "open" comes directly from the inherent properties of the MIKE mesh file (.mesh) and is parsed into a string array (including "open", "closed" and empty values) by the official toolkit readDfsuMesh, which is used to identify open boundaries (connecting external water areas) and land boundaries (solid walls). During interpolation calculations, for open boundaries ("open"), topologically adjacent grids are preferentially expanded to ensure physical continuity with external water bodies; for land boundaries ("closed"), internal grids within a radius R are preferentially expanded to cope with drastic gradient changes at the boundaries.

[0096] Input sampling point coordinates P(x p ,y p ), generate a candidate grid set {Mp} that satisfies: min(X) ≤x p ≤ max(X) and min(Y) ≤ y p ≤ max(Y), where X=(X1, X2, X3) and Y=(Y1, Y2, Y3) are the coordinates of the three vertices of the mesh.

[0097] Accurate position judgment: Call MATLAB's inpolygon function to traverse {Mp} and output the (in, on) flag: If in=1 and on=0, it is an internal point; if in=1 and on≥1, it is a boundary point, and the synchronous flag isBoundary=true.

[0098] According to the node relationship of mesh m, the adjacent mesh number set {Cm} is dynamically generated, where {Cm} is generated by the following rules: a) Get the three vertex node numbers of mesh m: N m1 , N m2 , N m3 ; b) Extract the node N from the node associated grid set {Mi} pre-stored in step S2 m1 、N m2 、N m3 The grid set {MNm1}、{M Nm2}、{M Nm3}; c) Take the union of the three sets: {Cm} = {M Nm1} ∪ {M Nm2} ∪{M Nm3}; d) Remove the grid m itself from {Cm}. An example of generating an adjacent grid set is as follows: Figure 2 :For the grid m=9652, its three nodes are: N m1 =5271,N m2 =5270,N m3 =4791, get the node associated grid set from the pre-stored topology: {M 5271} = {9652, 9653, 9654, 9655, 10608, 10607}; {M 5270} ={10605, 9651, 9652, 10606, 10607}; {M 4791} = {8760, 8761, 9650, 9651, 9653, 9652}, take the union: {Cm} = {9652, 9653, ..., 10607} ∪ {10605, ...} ∪ {8760, ...}, remove mesh m (9652): {Cm} = {9653, 9654, 9655, 10608, 10607, 10605, 9651, 10606, 8760, 8761, 9650}. {Cm} only includes adjacent meshes that directly share nodes with mesh m, excluding secondary adjacent meshes (adjacent meshes of adjacent meshes). The node association set {Mi} ensures that all meshes that share vertices are obtained, and removes the own mesh to avoid repeated calculations.

[0099] The dynamic expansion rule of the {Me} set is: when there is at least one non-boundary edge adjacent grid of the grid m, i.e., {Me1} is not empty, only {Me1} is used as the expansion set, and no internal grid {Me2} within the radius R is added. Only when all the adjacent edges of the grid m are boundary edges, i.e., {Me1} is an empty set, the internal grid search within the radius R is triggered. That is, when generating the expansion set {Me}, only when there is an adjacent grid with a non-closed boundary (i.e., an open boundary or an internal edge) of the grid m, it is included in {Me1}; if all the adjacent edges of the grid m are closed boundaries (such as land solid walls), {Me1} is empty, and the internal grid {Me2} within the radius R needs to be added. This design avoids the expansion of adjacent grids for the closed boundary, while ensuring the continuity of the open boundary water area. When the sampling point P is located at a complex boundary, such as a peninsula tip or a sharp corner grid, rule 1 is preferred, and the adjacent grid sharing the non-boundary edge is used, because it is directly related to the grid topology and has clear physical meaning; only when rule 1 returns an empty set, for example, all the adjacent edges of the grid m are boundary edges, rule 2 is enabled, i.e., the internal grid within the radius R is used to supplement the data. Exemplarily, the radius R can be 200 meters. Based on the internal grid data of {Cm'}, the boundary point of the sampling point P is calculated by the boundary-aware interpolation method: {Cm'} = {Cm} ∪ {Me}.

[0100] The boundary awareness is dynamically identified by pre-stored topology (boundary marking information). The non-boundary point is interpolated using the basic adjacent set {Cm}; the boundary point is interpolated after expansion to {Cm'}, covering all types of sampling points. Therefore, for boundary-aware interpolation and batch processing, interpolation is performed based on {Cm} or {Cm'}: 1. Physical quantity calculation: for flow velocity, water level and other physical quantities, three methods can be selected: nearest grid interpolation, node barycentric method interpolation, and MATLAB scatteredInterpolant. 2. Time step processing: multi-time step data is read by circularly reading three-dimensional array slices. 3. Multi-point batch processing: the S3-S5 steps are executed for the sampling point set, and the results are output to a CSV file. The interpolation module is equipped with an internal memory pool management, and the single-point calculation time is less than 10 milliseconds. The total time for batch processing of thousands of points is about 8 seconds.

[0101] Technical effects: The mat format pre-stored topology relationship eliminates real-time analysis overhead and improves data access efficiency. The boundary awareness mechanism improves the boundary point interpolation accuracy by dynamically expanding the adjacent set. The outer envelope rectangle screening reduces the grid traversal range to a local area, speeding up the positioning process. The memory pool optimization reduces repeated I / O operations and supports large-scale time step continuous processing. The batch processing loop structure realizes automatic multi-point extraction, reducing the cost of manual operation. The multi-interpolation strategy of physical quantities adapts to different accuracy requirements, balancing the calculation resources and the reliability of the results.

[0102] In the prior art, the user needs to analyze the dfsu binary file format in real time, and each time data is extracted, the file header and grid structure need to be read repeatedly, resulting in a single point positioning time of up to several minutes. The present embodiment eliminates the real-time analysis overhead by pre-converting the mat format and storing the global topology relationship (grid-node mapping, node coordinates, node associated grid set, boundary identification). The conversion is completed once in the preprocessing stage, and the conversion of 100,000 grids takes about 3 minutes. The topology data is reused for all subsequent sampling point extraction. Compared with the traditional method, the data reading time is reduced by 90%, and the advantage is significant for multi-time step extraction.

[0103] Boundary-aware interpolation mechanism: In the traditional MIKE post-processing module, when interpolating at the boundary point, the neighborhood data is often incomplete, resulting in distorted flow rate or water level jump. The present embodiment dynamically expands the adjacent grid set: when the sampling point is identified as a boundary point, all internal grids within a radius of 200 meters are automatically merged into the interpolation base set. This operation is triggered by using the pre-stored boundary marker information (node / edge identification) and combining the grid topology relationship to retrieve the adjacent internal grids in real time. Compared with the prior art, the boundary point interpolation accuracy is improved, the physical quantity continuity is improved, and it is especially suitable for complex boundary scenarios such as shorelines and gates.

[0104] Although the existing MATLAB toolkit supports script-based extraction, the user needs to manually write the loop logic, and the global data needs to be loaded repeatedly each time. The present embodiment supports multi-point loop extraction by pre-storing topology and physical quantity arrays: 1. Perform positioning and interpolation on the sampling point set; 2. Avoid repeated I / O by using array slicing operation for multi-time step data; 3. Results are automatically aggregated to structured output. This architecture realizes automatic batch processing of thousands of sampling points, which is more than ten times more efficient than manual operation, while avoiding the risk of intermediate file transmission errors.

[0105] Compared with the prior art, the pre-stored topology mechanism reduces the single point positioning time from minutes to seconds. The boundary-aware interpolation improves the accuracy of complex terrain data and eliminates the traditional boundary distortion phenomenon. The batch processing architecture supports automation for large-scale engineering applications, reducing the need for human intervention. The memory reuse design reduces disk I / O pressure and improves the stability of long-time sequence processing.

[0106] In one specific embodiment, the specific process of calculating the grid m where the sampling point P is located in step S3 is as follows:

[0107] S301, according to the coordinates (x p , y p ) of the sampling point P, screen the grid set { Mp}, the screening condition satisfies: min(X) ≤x p ≤ max(X) and min(Y) ≤y p≤ max(Y); wherein X=(X1, X2, X3), Y=(Y1, Y2, Y3) are the coordinates of the three vertices of the grid;

[0108] S302, traverse all the grids in {Mp}, determine the position relationship by MATLAB command [in, on] = inpolygon(x p , y p , X, Y);

[0109] If in=0, it means outside the grid;

[0110] If in=1 and on=0, it means inside the grid;

[0111] If in=1 and on≥1, it means on the edge or node of the grid, and the value of on accumulates to represent the number of shared edges / nodes;

[0112] Wherein, on takes values of 0, 1, 2, 3;

[0113] S303, determine the grid m and the position category according to the values of in and on;

[0114] S304, if the sampling point P is located on the boundary edge or the boundary node, mark it as a boundary point.

[0115] In the above embodiment, the screening condition is: the extreme value of the grid vertex coordinates satisfies min(X)≤x p ≤max(X) and min(Y)≤y p ≤max(Y). Vectorized comparison instructions can be selected for parallel processing, and the screening result {Mp} is stored in the memory buffer area. This module is deployed in the grid positioning preprocessing stage.

[0116] Call the inpolygon function of MATLAB to perform position judgment:

[0117] Input parameters: sampling point coordinates (x p , y p ), grid three vertex coordinate vectors [X1, X2, X3], [Y1, Y2, Y3];

[0118] Output flags: in=0 (external point), in=1 and on=0 (internal point), in=1 and on≥1 (boundary point);

[0119] on value definition: on=1 (single edge), on=2 (two-grid shared edge intersection), on=3 (three-grid shared node). The floating point tolerance uses the default value 1×10 -12 .

[0120] Boundary point identification and location classification: when on≥1, the boundary determination is triggered: if the sampling point is located on the edge of the grid, the boundary identification (string "open" or "closed") of the edge in the pre-stored topology is queried; if it is located on the node, the boundary identification (Boolean value 1 or 0) of the node is queried.

[0121] In the above embodiment, the candidate grid set {Mp} meeting the requirements is screened out, and the search range is narrowed to a local area. When the sampling point is located on the grid edge or the node (on≥1), the boundary marking mechanism is triggered synchronously. This process relies on the complete topology data pre-stored in step S2, eliminates the calculation burden of real-time boundary analysis, and ensures the accuracy of the boundary property identification of sensitive areas such as the shoreline and the gate in the water conservancy project. The output location category (internal point / boundary point) and boundary type provide key inputs for the boundary-aware interpolation in step S4. For example, the point marked as "open boundary" will trigger the dynamic expansion of the adjacent grid set {Cm}; the "land boundary" point will introduce the internal grid supplementary data within the radius R. Combined with the location classification result (such as the on value indicating the number of shared edges), various interpolation strategies are adaptively selected to reduce redundant calculations. This embodiment significantly reduces the computational complexity of grid search and topology analysis under the premise of ensuring numerical stability through the chain processing of pre-screening-precise positioning-boundary marking, laying a foundation for efficient interpolation.

[0122] In one specific embodiment, the interpolation method in step S4 includes any one of the following:

[0123] S401, adopt the nearest grid interpolation, which is divided into three cases:

[0124] A, if the sampling point P is inside the grid m or on the non-shared edge / node, directly take H P = H m , where H P is the physical quantity data value of the sampling point P, and H m is the physical quantity data value of the center point of the grid m;

[0125] B, if the sampling point P is on the shared edge of two grids, the associated node data is linearly interpolated according to the distance weight, specifically:

[0126] H P = [d(P, N2)H N1 + d(P, N1)H N2 ] / [d(P, N1) + d(P, N2)], where N1 and N2 are the two nodes of the shared edge of the two grids, d(P, N1) is the distance from the sampling point P to the node N1, d(P, N2) is the distance from the sampling point P to the node N2, H N1 is the physical quantity data value of the node N1, and H N2 is the physical quantity data value of the node N2;

[0127] C, if the sampling point P coincides with the node N, take H P = H N , H N is the physical quantity data value of the node N;

[0128] wherein the physical quantity data value H N of the node N is calculated by inverse distance weighted interpolation of the associated grid set {Mn}: , k is the number of grids in the grid set {Mn}, is the grid center physical quantity data value, is the distance between the node N and the grid Mn j center point, is the distance between the node N and the grid Mn i center point;

[0129] S402, using node value interpolation, the grid m in which the sampling point P is located has three nodes: N m1 , N m2 , N m3 , first calculate the grid m three node data H Nm1 , H Nm2 , H Nm3 , and then interpolate H P by the barycentric coordinate method or the inverse distance weighted method;

[0130] S403, call MATLAB command griddata or scatteredInterpolant to interpolate H P from the grid center point coordinates and data in {Cm}.

[0131] In the above embodiment, the nearest grid interpolation implementation process can use double-precision floating-point units to perform distance calculation, and the distance calculation tolerance is set to 1×10 -8 . For the case that the sampling point P is located inside the grid or the non-shared edge, the grid center point physical quantity value H m is directly read and assigned to H P = H m . When P is located on the shared edge of two grids, the Euclidean distances d1 and d2 of P to the edge nodes N1 and N2 are calculated, and linear interpolation is performed. The node physical quantity value H N is calculated by inverse distance weighted calculation: taking the previous node N as the center, searching the associated grid set {Mn}, and calculating the distance d i from each grid center point to N.

[0132] The node value interpolation implementation process: the barycentric coordinate method area calculation uses the vector cross product formula, and the area tolerance is fixed to 1×10 -10 . The total area S totaland sub-triangle area S1 (△PN m2 N m3 ), S2 (△N m1 PN m3 ), weight ω1=S1 / S total , ω2=S2 / S total , ω3=1-ω1-ω2. Inverse distance weighting method sets distance threshold β=1×10 -8 , returns node value when min(d k )<β; otherwise, calculates ω k =(1 / d k ²) / Σ(1 / d i ²). Two methods finally output according to H p =Σ(ω k ×H Nmk ).

[0133] MATLAB function call control: create dynamic memory pool management adjacent grid set {Cm} data, memory threshold set to 2GB. When the grid data is not cached, incrementally load from mat file; release data according to LRU policy when it is over limit. According to the base n of {Cm}, select the interpolation function: call griddata function when n≤N th , construct scatteredInterpolant object when n>N th . Performance threshold N th default 5000, automatically down to 2000 when system available memory is less than 8GB.

[0134] Technical effects: three interpolation methods form a complementary mechanism, nearest grid method improves the efficiency of regular point calculation, node interpolation method guarantees the accuracy of boundary points, and built-in function method adapts to large-scale adjacent set scenarios. Dynamic memory management reduces data repeated loading and optimizes memory resource utilization. Threshold adaptive mechanism balances the calculation stability under different hardware environments. Float tolerance setting eliminates numerical calculation error and ensures the accuracy of geometric relationship judgment. Special computing unit accelerates core algorithm execution to meet real-time processing needs. Hierarchical processing strategy covers from simple to complex application scenarios to enhance the practicability of the method.

[0135] In one specific embodiment, the barycentric coordinate method or inverse distance weighting method in step S402 is as follows:

[0136] a) Barycentric coordinate method:

[0137] , , ;

[0138] wherein, is the three nodes N m1 , Nm2 , N m3 , the triangle area enclosed by the sampling point P, the node N m2 , the node N m3 , the triangle area enclosed by the node N m1 , the node P, the node N m3 ;

[0139] b) inverse distance weighting method:

[0140] wherein i = 1, 2, 3;

[0141] For the barycentric coordinate method and the inverse distance weighting method, the physical quantity data value H p of the sampling point P is calculated in the following manner: H p = ×H Nm1 + ×H Nm2 + ×H Nm3 .

[0142] In the above embodiments, the barycentric coordinate method weight calculation can use the vector cross product formula to calculate the triangle area, and the floating point tolerance for area calculation is set to 1 x 10 -10 . The specific process is as follows: input the grid three-node coordinates N m1 (x1, y1), N m2 (x2, y2), N m3 (x3, y3) and the sampling point P(x p , y p ). Calculate the total area S total = 0.5 x | (x2-x1) (y3-y1) - (x3-x1) (y2-y1) |; Calculate the sub-triangle area S1 = 0.5 x | (x2-x p ) (y3-y p ) - (x3-x p ) (y2-y p ) | (corresponding to △P N m2 N m3 ); S2 = 0.5 x | (x p -x1) (y3-y1) - (x3-x1) (y p -y1) | (corresponding to △N m1 P N m3 ). The weight calculation formula is ω1 = S1 / S total ,

[0143] ω2 = S2 / S total , ω3 = 1-ω1-ω2.

[0144] Inverse distance weighting method weight calculation: Set distance threshold β = 1 x 10 -8 , exponential parameter fixed as 2, used for most water flow, wave simulation scenarios. Calculate the Euclidean distance from the sampling point to the three nodes: , k = 1, 2, 3. If there is d k <β, terminate the calculation and return the node value H Nmk ; Otherwise, calculate the weight denominator D = Σ (1 / d k ²), normalize the weight ω k = (1 / d k ²) / D.

[0145] Weighted summation unified output: The two methods share the weighted summation box, Hp=ω1×H Nm1 +ω2×H Nm2 +ω3×H Nm3 .

[0146] Technical effects: Explicit area formula guarantees the geometric accuracy of the barycentric method, suitable for structured regular grid scenarios. Distance threshold mechanism improves the numerical stability of the inverse distance method, avoiding near node singular value problems. Unified output framework simplifies hardware implementation, reduces switching overhead of the two methods. Fusion multiplication design improves computing throughput, shortens the single interpolation period to 5 clock cycles. Numerical verification enhances the reliability of the results, preventing physical quantity logic errors. Pipeline structure optimizes resource utilization, supporting real-time continuous data processing. Preloading mechanism hides data reading delay, improves large-scale point set processing efficiency.

[0147] A specific example is given, which extracts all variables of a sampling point at all time steps from a wave.dfsu file. The existing dfsu file stores three variables Hsg, Dir and T01 of 47505 grids, 26321 nodes and 8760 time steps. The calculation uses a desktop computer, and its main performance indicators are 12 threads, 3.10 GHz main frequency, 16 GB memory, and more than 1 TB of available storage space on the hard disk.

[0148] Step 1, convert the wave.dfsu file calculated by MIKE software into wave.mat format file which is easy to operate in matlab. The converted file includes three 47505*8784 two-dimensional arrays of Hsg, Dir and T01, the array row represents the grid number, the list represents the time step, Hsg is the wave height, Dir is the direction, and T01 is the time. In order to reduce the storage space, all variables are adjusted by an integer scaling factor (Hsg scaling factor 50, Dir takes 100, T01 takes 1), and then stored in 16-character positive integer format uint16. Finally, the size of the mat file is 1.30GB, compared with the original dfsu file of 4.67GB, the space is saved by 72%.

[0149] Step 2, convert the mesh format file input by MIKE software calculation into mesh mat format file which is easy to operate in matlab. The converted file includes NtoE, tn, xn, yn, xe, ye, which respectively represent the grid connected by each node, the node of each grid, the x coordinate of each node, the y coordinate of each node, the x coordinate of three nodes of each grid, and the y coordinate of three nodes of each grid.

[0150] Step 3, input the given sampling point coordinates, judge its grid M and position in the grid. In this example, the sampling point P coordinates are x P =114.782193, y P =22.313652, after preliminary screening, there are two grids that meet the conditions, Mp={9652, 9653}, the schematic diagram is shown in Figure 2 . The preliminary screening condition is: the maximum value of the grid node horizontal coordinate max(X)>x P , the minimum value min(X)<x P , the maximum value of the vertical coordinate max(Y)>y P , and the minimum value min(Y)<y P . Then use the matlab command (in, on)= inpolygon(x P , y P , X, Y) to judge the relationship between point P and grid 9652, 9653: when grid number m=9652, in=1, on=0; when grid number m=9653, in=0, on=0. It shows that the sampling point P is located in the grid m=9652, and not on the edge line and node.

[0151] Step 4, find the grid set { Cm} which shares the node with m, and based on all the data on the grid in the grid set, interpolate the data of the sampling point. In this embodiment, three interpolation methods are used, which are nearest grid interpolation, node interpolation and matlab built-in scattered point interpolation command interpolation.

[0152] The nearest grid interpolation is used, and the sampling point directly uses the value of grid m. The Hsg, Dir, and T01 of point P are the Hsg, Dir, and T01 values ​​of grid m.

[0153] Use the node data interpolation of the grid. First find the three nodes of the grid m, which are N m1 、N m2 、N m3 , and then find the common grid { M Nm1}、{M Nm2} and { M Nm3}. Nmi} data is interpolated to N mi (i = 1, 2, 3), then use N m1 、N m2 、N m3 The physical quantity data of point P is obtained by interpolation.

[0154] In this example, find three nodes N in the grid m=9652 m1 、N m2 、N m3 They are 5271, 5270 and 4791 respectively. The grid sets corresponding to the three nodes are { M 5271}= {9652, 9653, 9654, 9655, 10608, 10607}, {M 5270}={10605, 9651, 9652, 10606, 10607} and { M 4791} ={8760, 8761, 9650, 9651, 9653,9652}, see the schematic diagram Figure 2 .

[0155] {M Nmi} data is interpolated to N mi (i = 1, 2, 3), first calculate each grid node N mi Influence weight coefficient ω i (i = a positive integer from 1 to k), ω i To ensure that the sum of all weights is 1, the inverse of the distance from node N to the center of each grid is used as the weight.

[0156] In this example, node N m1 and the meshes that share this node { M Nm1} as an example, {M Nm1 The weights of the six grids in} are 0.167, 0.186, 0.157, 0.163, 0.151, and 0.176 respectively.

[0157] Node N m1 The data such as wave height Hsg can be expressed as:

[0158] Hsg Nm1 = 0.167×Hsg 9652 + 0.186×Hsg 9653 + 0.157×Hsg 9654 + 0.163×Hsg 9655 +0.151×Hsg 10608 + 0.176×Hsg 10607 ;

[0159] The data of N m1 , N m2 , N m3 are obtained by using the barycentric coordinate method to obtain the data of P points, as shown in FIG. 2. This embodiment discloses two interpolation methods, which are the barycentric coordinate interpolation method and the inverse distance weighted interpolation method, and both can be expressed as Hsg, for example: Figure 3

[0160] Hsg P = ω1×Hsg Nm1 +ω2×Hsg Nm2 +ω3×Hsg Nm3 . The difference between the two interpolation methods lies in the different values of the weight coefficient ω.

[0161] The physical quantity data of the P points are obtained by using the matlab built-in command griddata or scatteredInterpolant, and the command is as follows, for example, the variable Hsg:

[0162] HsgP = griddata(X,Y, Hsg, method);

[0163] or F = scatteredInterpolant(X, Y, Hsg, method) ; HsgP = F(x, y)。

[0164] Wherein X, Y and H respectively represent the grid set { M Nm1}, { M Nm2} and {M Nm3 ​The center coordinates and Hsg data of each grid in the set { Cm} are created in the memory pool, and the initial memory pool threshold is set to 2 GB.

[0165] In one embodiment, the following steps are performed before calling the MATLAB interpolation command in step S403:

[0166] 1) Create a memory pool to cache the center point coordinates and physical quantity data of the grids in the set { Cm}, and set the initial memory pool threshold to 2 GB.

[0167] 2) Add a system memory monitoring module to detect the available physical memory of the MATLAB process in real time.

[0168] 3) Dynamic adjustment method: a. If the memory pool usage is greater than the current memory pool threshold, trigger LRU release, delete the cache data that has not been accessed within the last 15 minutes, and repeat the deletion until the memory pool usage is less than or equal to 90% of the current memory pool threshold; b. i) If the system available memory is less than 4 GB, suspend the memory pool function and read the mat file directly, and skip the subsequent downgrade operation; ii) If the system available memory is greater than or equal to 4 GB and less than 8 GB, reduce the performance threshold N th from 5000 to 2000; compress the memory pool threshold from 2 GB to 1 GB; c. When the system available memory is greater than or equal to 8 GB for 5 minutes, restart the memory pool function, restore the memory value threshold to 2 GB, and restore the performance threshold N th to 5000.

[0169] In the above embodiment, in MATLAB, containers.Map is used to create a key-value storage structure, and the memory pool data is managed using the key-value storage structure, where the key is the grid number and the value is a structure containing the grid center point coordinates and physical quantity data. The initial memory pool threshold is set to 2 GB. Memory usage is calculated by accumulating the total number of bytes of all cached data, and the size of a single data is estimated to be 20 bytes. When the total usage exceeds 2 GB due to the addition of new data, the system triggers the LRU release mechanism. The memory pool is deployed in the memory management module of the computing server and shares the memory space with the MATLAB process. The hash table uses open addressing method to solve the conflict, and the cache entries are stored in a continuous memory block in order of grid number. When running on a 64-bit Windows or Linux system, the memory pool realizes physical memory mapping through the memmapfile function of MATLAB.

[0170] The working process is that when the system loads the grid data, firstly, it is inquired whether the numbered cache exists in the memory pool. If it exists, the data is directly read; if it does not exist, the physical quantity data is read from the MAT file, and the coordinates and physical quantity are packaged into a structure and stored in the memory pool. Each data access updates the last access timestamp. After adding new data, the total memory pool occupancy value is calculated in real time, and when it is detected that it exceeds 2 GB, the LRU release process is started. The initial threshold of 2 GB is set according to the typical workstation configuration, and can be adapted to 50% of the safe margin of a 16 GB memory device. The calculation of a single data of 20 bytes includes: 8 bytes of longitude, 8 bytes of latitude, and 4 bytes of physical quantity value. The timestamp accuracy is millisecond level, and the system clock GetTickCount is used to realize it.

[0171] The system memory monitoring module is added, which detects the available physical memory in real time through the operating system API, and the sampling frequency is 1 time per second. The monitoring module calls the GlobalMemoryStatusEx function of the Windows system or the sysinfo function of the Linux system to obtain the available memory value accurate to MB. The module is integrated into the daemon thread of the MATLAB main process, and non-blocking query is adopted to avoid calculation delay. When it is detected that the system available memory is lower than 8192 MB, two operations are synchronously performed: the interpolation performance threshold N th It is down-regulated from 5000 to 2000, and the memory pool threshold is compressed from 2 GB to 1 GB. When the available memory is lower than 4 GB, the memory pool function is immediately suspended, and the mode of directly reading the MAT file is switched.

[0172] Function test: On a test machine equipped with 16 GB of memory, simulate the memory pressure scene: 1. Run the MATLAB process with 1.5 GB of memory occupation; 2. Start the memory loading tool to consume 12 GB of memory; 3. Record the system response: when the available memory decreases to 7.8 GB, N th It is adjusted from 5000 to 2000 within 200 ms; when the available memory decreases to 3.9 GB, the memory pool is closed within 500 ms.

[0173] When the memory pool occupancy exceeds the current threshold, the LRU (Least Recently Used) release mechanism is triggered. The system traverses all cache entries, and deletes the data whose last access time is earlier than the current time by 900 seconds. The release operation is executed in a loop until the memory pool occupancy decreases to 90% of the threshold. For example, it is released to 1843 MB when the threshold is 2 GB, and to 922 MB when the threshold is 1 GB. During the release process, two-level index optimization is adopted: firstly, an ordered linked list of last access time is established to quickly locate the timeout data; and then a hash table is used to realize O(1) complexity deletion. For the scene where the system memory is lower than 8192 MB, the performance threshold N thDown to 2000, forced to choose griddata interpolation function instead of scatteredInterpolant. When the available memory is less than 4096MB, the memory pool is completely closed, and all data requests directly access the hard disk MAT file.

[0174] LRU time window 900 seconds with quartz crystal oscillator, you can choose XRCGB series crystal of Murata Manufacturing. Data deletion in the memory release process using the principle of TRIM instruction of solid state disk, avoid the actual data erase. When the system available memory continues to exceed 8192MB for 300 seconds, automatically restore the initial configuration: th =5000, memory pool threshold = 2GB, and reactivated cache function.

[0175] Technical effects: memory hierarchical management to ensure 8GB memory device stable operation of fifty thousand grid level calculation, LRU release strategy to maintain cache hit rate above 85%; threshold linkage mechanism balance interpolation speed and memory consumption; abnormal recovery function reduces the need for human intervention; direct file reading mode to avoid memory depletion collapse.

[0176] In one embodiment, in step S303, when on≥1, the grid set { Ms} sharing the edge / node is recorded synchronously, the specific process is:

[0177] a) When the sampling point P is located on the grid edge, according to the two node numbers N1 and N2 of the grid edge, extract the grid containing N1 or N2 from the node-related grid set {Mi} pre-stored in step S2, and take the union to generate { Ms};

[0178] b) When the sampling point P coincides with the node N, directly call the node-related grid set { Mn} pre-stored in step S2 as { Ms};

[0179] In step S4, based on {Ms}, the adjacent grid set {Cm} is dynamically generated, and then the boundary-aware interpolation is performed.

[0180] In the above embodiment, when the sampling point P is located on the grid edge, the associated retrieval is performed according to the two node numbers N1 and N2 of the edge. A pre-stored node-grid relationship mapping table can be selected, in which each node entry stores all grid number sets associated with it. From the mapping table, the grid set Set N1 containing node N1 and the grid set Set N2 containing N2 are extracted respectively. The union operation of the two sets generates the edge-related grid set {Ms}, that is, {Ms} = Set N1 ∪ Set N2The operation is deployed in the topology relationship processing unit of the grid positioning module, and is realized by hash table query to achieve millisecond-level response. In implementation, vectorized set operation is used to avoid performance loss caused by loop traversal.

[0181] Node-coincident grid set calling: when the sampling point P coincides with the grid node N, the pre-stored node-related grid set {Mn} is directly called as {Ms}. A pre-constructed node topology database can be selected, in which each node ID is associated with a list of all adjacent grid numbers. By indexing the database with the ID of node N, the pre-calculated {Mn} set is directly returned. This process is integrated into the output link of the position judgment module and is automatically triggered when on=3 is detected. The database is stored in the server memory cache area, and fixed-length integer arrays are used to store grid numbers. A single query takes less than 0.1 milliseconds.

[0182] The generated {Ms} set is directly input as the basis grid set {Cm} into the interpolation calculation module. The memory pool management mechanism can be reused, and only the physical quantity data of the grids in the {Ms} set is loaded. Before interpolation execution, duplicate grid numbers are automatically filtered to ensure the uniqueness of the set. This data transmission is deployed in the initialization stage of the interpolation pipeline, and the grid number array is directly transmitted through the memory pointer. In implementation, a pre-allocated continuous memory block is used to store the number list to avoid the delay caused by dynamic array expansion. Therefore, for step S4, when {Cm} is dynamically generated, if there is a pre-calculated {Ms} set, the {Ms} set is used as the basis of the node-related set to replace the re-retrieval from the nodes of grid m, that is, when the sampling point P is located at the edge / node, the pre-stored {Ms} set is called and input into the dynamic generation process as the node-related set, replacing the re-retrieval from the nodes of grid m.

[0183] Technical effects: The node topology pre-storage mechanism eliminates the calculation overhead of dynamically retrieving adjacent grids, improving the efficiency of boundary point processing. Directly calling the pre-stored set avoids repeated query of topology relationship, shortens the interpolation preparation time. Continuous memory storage optimizes data transmission efficiency and reduces the processing delay of large-scale grid sets. Set operation ensures the integrity of the shared edge-related grids, ensuring the interpolation accuracy of complex boundary scenarios. The pre-computed database reduces real-time calculation load and enhances the overall stability of the system.

[0184] In one specific embodiment, a verification mechanism is added to the nearest grid interpolation of step S401: when the sampling point P is located at the grid edge or node, i.e., on≥1, the following operations are performed:

[0185] a) Synchronously use steps B and S402 in step S401 to calculate interpolation results H P1 and H P2 , respectively;

[0186] b) calculating the relative deviation δ = |H P1 - H P2 | / max(|H P1 |, |H P2 |, ε), where ε is a zero minimum prevention value;

[0187] c) if δ > η, where η is a preset tolerance threshold, triggering a grid topology review:

[0188] c1) rechecking the grid m and the shared edge / node relationship determined in step S303;

[0189] c2) reviewing the integrity of the grid sets {M N1} and {M N2} associated with the nodes N1 and N2 in step S2;

[0190] d) re-executing step B in step S401 or step S402 with the reviewed topology structure, and outputting the final value with a review mark.

[0191] In the above embodiment, for the synchronous calculation of the double interpolation results, when the sampling point is located at a grid edge or node, two interpolation calculations are synchronously performed. The linear interpolation method can be selected to calculate the first result value H P1 , and specifically, the physical quantity values H N1 and H N2 of the shared edge nodes N1 and N2 are obtained, the distances d1 and d2 of the sampling point to the nodes are calculated, and the formula H P1 = (d2×H N1 + d1×H N2 ) / (d1+d2) is used to output. At the same time, the barycentric coordinate method is selected to calculate the second result value H P2 , and specifically, the physical quantity values of the three nodes of the grid where the sampling point is located are calculated, and the weighted sum is output through the area weight of the sub-triangle. The two calculation methods are deployed in the same interpolation processing unit, share the node data cache, and the execution time is not more than 5 milliseconds.

[0192] For the relative deviation determination and threshold triggering, the relative deviation δ of the double results is calculated, and the formula is: δ = |H P1 - H P2 | / max(|H P1 |, |H P2 |, ε), where the zero minimum prevention value ε is set to 1×10 -10 . The preset tolerance threshold η is fixed at 0.05. When δ > η, for example, η = 0.05, the topology review process is automatically triggered. The deviation calculation module is integrated in the result comparison unit, and double-precision floating-point operation is used to ensure accuracy. When implemented, two levels of determination are set: first, check if |H P1 - H P2 | > 1×10-6 (avoiding false triggering by minor fluctuations), and then the formal delta value is recalculated.

[0193] For grid topology verification execution, the verification process includes two levels of operations: first, recheck the topology relationship of the grid m where the sampling point is located, and confirm the shared edge / node ownership by backtracking the vertex coordinates and comparing them with the pre-stored grid table; second, verify the integrity of the grid set {M N1} and {M N2} associated with nodes N1 and N2, and check whether there are missing associated grids. When a topology error is found, the associated set is rebuilt and interpolation is re-executed (using the original step B or S402 method), and the result is marked with the verification label "verified". This process is deployed in a separate verification module, and the average processing time is 20 milliseconds. The verification data comes from the topology database converted from the original mesh file.

[0194] Technical effects: The double-interpolation cross-validation mechanism effectively identifies potential topology errors and improves the reliability of boundary point data. The relative deviation threshold determination avoids unnecessary verification operations, balancing calculation efficiency and accuracy. Two-level topology verification ensures accurate error positioning and covers common problems such as node association loss. The result labeling mechanism provides data quality identification for subsequent analysis. The floating-point minimum quantity setting eliminates zero-value interference and ensures the stability of deviation calculation. The automatic verification process reduces the need for manual intervention and improves the degree of batch processing automation.

[0195] In one specific embodiment, a distance threshold determination mechanism is added in the inverse distance weighting method calculation of step S402:

[0196] a) Calculate the distances d1, d2, and d3 of the sampling point P to the three nodes (N m1 , N m2 , N m3 ) of the grid m;

[0197] b) If d k <β, k ∈ {1, 2, 3}, β is a small positive number related to floating-point calculation precision, and the value is 10 -6 -10 -10 , then it is determined that P coincides with node N mk ;

[0198] c) Directly assign H P =H Nmk , and terminate the subsequent weighting calculation process;

[0199] d) Otherwise, calculate the inverse distance weighting value according to the original inverse distance weighting method of step S402.

[0200] In the above embodiment, the distances d1, d2, and d3 of the sampling point P to the three nodes (N m1 , N m2, N m3 ). The distance calculation can be performed by a double-precision floating-point operation unit, and the formula is: , where k = 1, 2, 3. The floating-point precision threshold β is set to 1 x 10 -8 .

[0201] When there is d k < β (k ∈ {1, 2, 3} ), it is determined that the sampling point P coincides with the node N mk . The pre-stored physical quantity value H Nmk of the node is directly read and assigned to H P . A pre-constructed node physical quantity database can be selected to index the value in real time through the node ID. The determination module is integrated at the output end of the distance calculation, and a parallel comparison circuit is used to synchronously detect the three distance values. A hardware-accelerated threshold comparator is set in implementation, and the response time is less than 0.01 milliseconds. The assignment operation skips all subsequent weighted calculation processes and directly outputs the result to the storage area.

[0202] If all d k ≥ β, the standard inverse distance weighted method is used for calculation. A two-stage pipeline structure can be selected: the first stage calculates the reciprocal of the square of the distance 1 / d k 2 , and the second stage performs the normalized weight ω k = (1 / d k 2 ) / Σ(1 / d i 2 ). Finally, H P = Σ(ω k x H Nmk ) is output through a multiplication-addition operation unit. This process is deployed in the interpolation core processor, which supports a double-precision floating-point fused multiply-add instruction. A zero-distance protection mechanism is added in implementation, which automatically interrupts the current process branch when d k < β.

[0203] The processing priority rule when the sampling point is located at the boundary and the grid size is less than the distance threshold β. Specifically, when the sampling point P is marked as a boundary point and the size of the grid m is less than β, where β is a very small positive number related to the floating-point calculation precision, and its value interval is 10 -6 to 10 -10 , the system preferentially triggers the distance threshold determination mechanism. This mechanism requires calculating the distances d1, d2, and d3 of the sampling point P to the three nodes N m1 , N m2 , and N m3 of the grid m. If there is any distance d k less than β, the sampling point physical quantity H PEqual to the corresponding node physical quantity and terminate the subsequent process; if all distances d k are greater than or equal to β, continue to execute the boundary-aware interpolation process, dynamically create an extended grid set Cm' equal to the union set of the base adjacent grid set Cm and the supplementary grid set Me, and perform interpolation calculation based on it. This rule ensures that in the case of a small grid at the boundary, the efficient node direct assignment strategy is preferred, saving about 90% of the calculation time compared to boundary-aware interpolation, while the hierarchical processing takes into account the accuracy requirements of the boundary gradient region.

[0204] Technical effects: The floating-point precision threshold effectively avoids the risk of division by zero caused by small distances, enhancing the numerical stability of the algorithm. The node coincidence determination mechanism eliminates unnecessary weighted calculations, improving the efficiency of boundary point processing. The hardware-accelerated comparator shortens the response time of the determination, meeting the real-time calculation requirements. The two-stage pipeline design optimizes the throughput of inverse distance weighted calculation. The process interruption control reduces the energy consumption of invalid operations. The pre-stored node database ensures the immediacy of assignment operations, ensuring low-latency result output.

[0205] In one specific embodiment, a sample point pre-screening module is added before step S3, which specifically performs:

[0206] a) Based on all node coordinates stored in step S2, construct a global outer envelope rectangle R:

[0207] R = { (x, y) | x∈[min(X all )-Ψ, max(X all )+Ψ], y∈[min(Y all )-Ψ, max(Y all )+Ψ]},

[0208] Where X all is the set of all node x coordinates, Y all is the set of all node y coordinates, and Ψ is the safety margin;

[0209] b) Create a spatial index structure: evenly divide the outer envelope rectangle R into Q×L grid cells, and record the grid number set covered by each cell;

[0210] c) For each sample point P(x p , y p ):

[0211] c1) If (x p , y p ) R, directly mark it as an invalid point and skip the subsequent processing;

[0212] c2) If (x p , y p) ∈ R, according to the spatial index positioning the grid cell to which it belongs;

[0213] c3) Only load the grid set associated with the grid cell to perform the grid traversal of step S3.

[0214] In the above embodiment, the outer envelope rectangle R is constructed based on all node coordinates, and a double-precision floating-point calculation unit can be selected to perform extreme value operation: X min = min(X all ) - Ψ, X max = max(X all ) + Ψ; and Y min / Y max Similarly. This module is deployed in the preprocessing stage, and the output is stored in the spatial index database. When implemented, a parallel reduction algorithm is used to accelerate the extreme value search, and the processing time of a million nodes is less than 0.5 seconds.

[0215] The purpose of setting the safety margin is to expand the range of the outer envelope rectangle (Ψ>0) to provide a buffer space for the boundary region. Even if the sampling point coordinates slightly exceed the theoretical boundary due to floating-point errors (but within the range of Ψ), they will still be considered valid points. When constructing the global outer envelope rectangle, the value of the safety margin Ψ needs to take into account the non-uniformity of the grid size in the calculation domain. When the size of the triangular grid varies significantly, directly using the global average size to determine the value of Ψ may lead to mismatch of the local region margin. Therefore, a hierarchical safety margin strategy can be implemented: first, divide the calculation domain into sub-regions with similar grid sizes based on the grid size distribution characteristics, for example, distinguish the dense grid near the coast from the sparse grid in the open sea through clustering algorithm; then calculate the local average size D local for each sub-region independently, and set the region-specific margin according to the formula Ψ local = max(α × D local , Ψ min ), where the safety factor α is recommended to be 6.0 to cover typical geometric deformation and floating-point tolerance, and the minimum margin Ψ min is fixed at 10 meters to ensure basic fault tolerance. This strategy reduces Ψ local to 30 meters for fine grid regions, such as the near-shore area with an average size of 5 meters, to avoid excessive expansion of the candidate grid set; while ensuring that Ψ local is increased to 600 meters for sparse regions, such as the open sea area with an average size of 200 meters, to prevent valid points from being mistakenly excluded. For extremely distorted grids or sensitive areas with open boundaries, an additional 20% margin buffer can be added.

[0216] Divide the rectangle R uniformly into Q × L grid cells (Q = L = 200), with cell size Δx = (X max -X min ) / Q and Δy = (Y max -Y min) / L. A two-dimensional matrix storage unit-grid mapping relationship can be selected: traverse all grids, determine the covered unit region according to the outer envelope rectangle thereof, and add the grid ID to the record list of the covered unit. The index database is deployed in the memory cache area, and a dynamic array is used to store the unit grid set. The construction process is completed synchronously in the topology conversion stage.

[0217] For each sampling point P(x p , y p ): 1. Out-of-domain determination: if x p does not belong to [X min , X max ] or y p does not belong to [Y min , Y max ], mark "invalid point" and skip the subsequent process; 2. Unit positioning: calculate the unit index i = [(x p - X min ) / Δx], j = [(y p - Y min ) / Δy]; 3. Local loading: extract the unit (i, j) related grid set from the spatial index, and only load the set to perform grid traversal. This module runs on the sampling point input interface, and the invalid point determination is realized by a hardware comparator to achieve a nanosecond-level response.

[0218] Technical effects: The outer envelope rectangle quickly excludes out-of-domain points, avoiding invalid calculation resource consumption. The spatial index converts global search into local traversal, and the candidate set size is reduced to one thousandth in a million grid scene. The unit mapping mechanism supports O(1) complexity positioning, significantly improving real-time processing capability. Dynamic memory loading reduces the topology data access pressure, optimizing the stability of large-scale engineering applications. The safety margin setting prevents boundary misjudgment caused by floating-point errors, ensuring the reliability of screening.

[0219] In one of the embodiments, when extracting multiple sampling point data in step S6, a parallel processing flow is performed:

[0220] S601, according to the pre-constructed spatial index structure, divide the sampling point set into blocks according to the spatial grid unit;

[0221] S602, call the parallel computing unit, and synchronously perform grid positioning and physical quantity interpolation operations on the sampling points in each block;

[0222] S603, use a pipeline loading mechanism for multi-time step data, and asynchronously read the next time step data;

[0223] S604, aggregate the interpolation results of each block and output.

[0224] In the above embodiment, the sampling point set is grouped according to the grid cell to which it belongs based on a pre-constructed spatial index structure. The number of blocks is equal to the number of CPU cores, such as 16 cores divided into 16 blocks, and each block contains less than or equal to 1000 points. A load balancing algorithm can be used to distribute the point set: all spatial index cells are traversed, the number of sampling points in each cell is counted, the cells are sorted in ascending order of the number of points, and then merged into blocks until the upper limit of 1000 points is reached. The block operation is deployed in the task scheduling module, and the output is a block point set array stored in the shared memory area.

[0225] The parallel computing unit (such as MATLAB parallel pool) is called to process each block: 1. An independent thread is assigned to each block, and the grid positioning and physical quantity interpolation of the points in the block are executed synchronously; 2. The main thread monitors the status of each block and dynamically allocates idle threads; 3. Time step pipeline: while the current thread is processing time step t, the background thread asynchronously loads the physical quantity data of time step t+1 into the cache; the cache depth is fixed at 2 steps, and the preloaded data uses a double-buffer storage structure; the parallel scheduler runs in the operating system process management layer, and asynchronous I / O is executed through an independent DMA channel.

[0226] After each thread completes the block calculation: 1. The interpolation result matrix of the sampling points in the block is output, in the format [point ID, time step, physical quantity value]; 2. The main thread concatenates all block matrices in the original sampling point ID order; 3. The results are written to a CSV file, and the block memory is released at the same time. The aggregation module is integrated into the output interface, and a non-blocking disk write queue is used to avoid I / O waiting.

[0227] Technical effects: The spatial block strategy maximizes the use of multi-core parallel computing capability, shortens the processing time of thousands of points, and prevents single-thread overload through a load balancing mechanism, improving overall resource utilization. Double-buffer asynchronous loading hides data reading delay, ensuring the continuity of time series processing. Non-blocking output reduces the interference of result saving on the calculation process. Dynamic thread allocation adapts to system load changes, enhancing the stability of large-scale task processing. The memory release mechanism reduces the risk of resource occupation during long-time running.

[0228] The number of devices and processing scale described herein are used to simplify the description of the present application. Applications, modifications and variations of the present application are obvious to those skilled in the art.

[0229] Although the embodiments of the present application have been disclosed as above, they are not limited to the applications and embodiments listed in the specification, and can be fully applied to various fields suitable for the present application, and additional modifications can be easily realized by those skilled in the art, therefore the present application is not limited to specific details without departing from the general concept defined by the claims and equivalent scope.

Claims

1. A method for quickly extracting point data from a dfsu file, characterized in that: The following steps are involved: S1. Convert the result file in dfsu format generated by MIKE software into a mat format file that can be read by MATLAB; S2. Convert the mesh file input by the MIKE software into a mat format file, pre-store the global mesh topology, including: each mesh number and its corresponding three node numbers, the coordinates of all nodes, the mesh number set {Mi} associated with each node, and boundary marking information, wherein the boundary marking information includes boundary node identifiers and boundary edge identifiers; S3. Input the coordinates of the sampling point P and calculate the grid m where it is located and its position in the grid based on the pre-stored global grid topology relationship; S4. Dynamically generate the adjacent grid number set {Cm} based on the node relationship of grid m, and calculate the physical quantity data of sampling point P through the boundary-aware interpolation method based on the grid data in {Cm}; S5. If extracting data of multiple time steps, perform array operations within the time step or execute step S4 in a loop; S6. If multiple sampling point data are to be extracted, steps S3 to S5 are executed in a loop and the extraction results are output; The boundary-aware interpolation method includes the following steps: a. If P is a boundary point, then create an extended grid set {Cm'} = {Cm} ∪ {Me}; {Me} acquisition step: first obtain all adjacent grids of grid m that share non-boundary edges, recorded as set {Me1}; if {Me1} is an empty set, then additionally obtain all internal grids within a radius R with P as the center, recorded as set {Me2}; finally, {Me} = {Me1} ∪ {Me2}; b. Interpolate based on {Cm'} as the boundary point; The specific process of calculating the grid m where the sampling point P is located and its position in the grid in step S3 is as follows: S301, according to the coordinates (x p ,y p ) Filter the grid set {Mp}, and the filtering condition satisfies: min(X)≤x p ≤max(X) and min(Y)≤y p ≤max(Y); where X = (X1, X2, X3) and Y = (Y1, Y2, Y3) are the coordinates of the three vertices of the mesh; S302, traverse all grids in {Mp}, and use the MATLAB command [in, on] = inpolygon (x p ,y p , X, Y) to determine the position relationship: If in=0, it means outside the grid; If in=1 and on=0, it means inside the grid; If in=1 and on≥1, it means on the grid edge or node, and the cumulative on value represents the number of shared edges / nodes; Among them, the value of on is 0, 1, 2, 3; S303, determine the grid m and location category based on the in and on values; S304: If the sampling point P is located at a boundary edge or a boundary node, it is marked as a boundary point.

2. The method for quickly extracting dfsu file point data according to claim 1, characterized in that: The interpolation method in step S4 includes any of the following: S401, using the nearest grid interpolation, divided into three cases: A. If the sampling point P is inside the grid m or on a non-shared edge / node, directly take H P =H m , where H P is the physical quantity data value of the sampling point P, H m is the physical quantity data value of the center point of grid m; B. If the sampling point P is on the shared edge of two grids, the associated node data is linearly interpolated according to the distance weight, specifically: H P =[d(P, N2)H N1 +d(P,N1)H N2 ] / [d(P, N1)+d(P, N2)], where N1 and N2 are two nodes on the shared edge of the two grids, d(P, N1) is the distance from the sampling point P to the node N1, d(P, N2) is the distance from the sampling point P to the node N2, and H N1 is the physical quantity data value of node N1, H N2 is the physical quantity data value of node N2; C. If the sampling point P coincides with the node N, take H P =H N , H N is the physical quantity data value of node N; Among them, the physical quantity data value H of node N N Compute by inverse distance weighted interpolation the associated grid set {Mn}: k is the number of grids in the grid set {Mn}, is the physical quantity data value at the center of the grid, d(N, Mn j ) is the node N and the grid Mn j Distance from the center point, d(N, Mn i ) is the node N and the grid Mn i distance from the center point; S402, using node numerical interpolation, the three nodes of the grid m where the sampling point P is located are: N m1 、N m2 、N m3 , first calculate the three-node data H of the grid m Nm1 , H Nm2 , H Nm3 , and then interpolate H by using the barycentric coordinate method or the inverse distance weighted method P ; S403, call MATLAB command griddata or scatteredInterpolant, with the coordinates of the grid center point in {Cm} and the data interpolation value H P ; Among them, physical quantities include wave height, wave direction, wave period, flow direction, flow velocity, water level and concentration.

3. The method for quickly extracting dfsu file point data according to claim 2, characterized in that: The barycentric coordinate method or the inverse distance weighted method in step S402 is specifically: a) Barycentric coordinate method: ω3=1-ω1-ω2; Among them, S(ΔN m1 N m2 N m3 ) are three nodes N m1 、N m2 、N m3 The area of ​​the triangle, S(ΔPN m2 N m3 ) is the sampling point P, node N m2 Node N m3 The area of ​​the triangle, S(ΔN m1 PN m3 ) is node N m1 , node P, node N m3 The area of ​​the triangle is ω1, which is the area of ​​the sampling point P relative to the grid node N. m1 The barycentric coordinate weight of the sampling point P is the weight of the grid node N. m2 The barycentric coordinate weight of the sampling point P is the weight of the grid node N. m3 The weight of the center of gravity coordinates; b) Inverse distance weighted method: Where i = 1, 2, 3, d(P, N mi ) is the sampling point P to node N mi The distance, d(P, N mj ) is the sampling point P to node N mj The distance, ω i is node N mi The weight of the sampling point P; For the barycentric coordinate method and the inverse distance weighted method, the physical quantity data value H of the sampling point P p The calculation method is: H p =ω1×H Nm1 +ω2×H Nm2 +ω3×H Nm3 .

4. The method for quickly extracting dfsu file point data according to claim 2, characterized in that: In step S403, before calling the MATLAB interpolation command, perform the following steps: 1) Create a memory pool to cache the center point coordinates and physical quantity data of the grids in the grid number set {Cm}. The initial memory pool threshold is set to 2GB. 2) Add a system memory monitoring module to detect the available physical memory of the MATLAB process in real time: 3) Dynamic adjustment method: a) If the memory pool occupancy is greater than the current memory pool threshold, trigger LRU release and delete cache data that has not been accessed within 15 minutes before the current time. Repeat the deletion until the memory pool occupancy is ≤ 90% of the current memory pool threshold. b) If the system available memory is less than 4GB, suspend the memory pool function and directly read the mat file instead, and skip the subsequent downgrade operation. ii) If the system available memory is greater than or equal to 4GB and less than 8GB, increase the performance threshold N th Reduced from 5000 to 2000; compressed the memory pool threshold from 2GB to 1GB; c. When the system available memory is ≥8GB for 5 minutes, restart the memory pool function, restore the memory value threshold to 2GB, and restore the performance threshold N. th is 5000.

5. The method for quickly extracting dfsu file point data according to claim 1, characterized in that: In step S303, when on≥1, the grid set {Ms} of shared edges / nodes is synchronously recorded. The specific process is as follows: a) When the sampling point P is located at the edge of the grid, according to the two node numbers N1 and N2 of the grid edge, the grid containing N1 or N2 is extracted from the node-associated grid set {Mi} pre-stored in step S2, and the union is taken to generate {Ms}; b) When the sampling point P coincides with the node N, the node-associated grid set {Mn} pre-stored in step S2 is directly called as {Ms}; In step S4, based on {Ms}, a set of adjacent meshes {Cm} is dynamically generated, and then boundary-aware interpolation is performed.

6. The method for quickly extracting dfsu file point data according to claim 2, characterized in that: A verification mechanism is added to the nearest grid interpolation in step S401: when the sampling point P is located at a grid edge or node, that is, on≥1, the following operations are performed: a) Synchronously use step B in step S401 and step S402 to calculate the interpolation result H P1 and H P2 ; b) Calculate the relative deviation δ = |H P1 ﹣H P2 | / max(|H P1 |,|H P2 |, ε), ε is the anti-zero minimum; c) If δ>η, where η is the preset tolerance threshold, a grid topology review is triggered: c1) rechecking the grid m and the shared edge / node relationship determined in step S303; c2) Review the mesh set {M associated with nodes N1 and N2 in step S2 N1 },{M N2 }completeness; d) Re-execute step B or step S402 in step S401 with the verified topological structure, and output the final value with the verification mark.

7. The method for quickly extracting dfsu file point data according to claim 3, characterized in that: In the inverse distance weighted method calculation of step S402, a distance threshold determination mechanism is added: a) Calculate the number of nodes from the sampling point P to the grid m (N m1 , N m2 , N m3 ) distances d1, d2, d3; b) If there is d k <β, k∈{1, 2, 3}, β is a minimum positive number related to floating-point calculation accuracy, and its value is 10 -6 -10 -10 , then determine whether P and node N mk coincide; c) Directly assign H P =H Nmk , and terminate the subsequent weighted calculation process; d) Otherwise, the inverse distance weighted value is calculated according to the inverse distance weighted method of the original step S402.

8. The method for quickly extracting dfsu file point data according to claim 1, characterized in that: Add a sampling point pre-screening module before step S3, specifically perform the following: a) Based on the coordinates of all nodes stored in step S2, construct the global outer envelope rectangle R: R={(x,y)|x∈[min(X all )-Ψ,max(X all )+Ψ],y∈[min(Y all )-Ψ,max(Y all )+Ψ]}, Among them, X all is the set of x coordinates of all nodes, Y all is the set of y coordinates of all nodes, Ψ is the safety margin; b) Create a spatial index structure: divide the outer envelope rectangle R evenly into Q×L grid cells and record the grid number set covered by each cell; c) For each sampling point P(x p ,y p ): c1) If Directly mark it as an invalid point and skip subsequent processing; c2) If (x p ,y p )∈R, locate the grid cell according to the spatial index; c3) Only the grid set associated with the grid unit is loaded to perform the grid traversal in step S3.

9. The method for quickly extracting dfsu file point data according to claim 8, characterized in that: In step S6, when extracting data from multiple sampling points, a parallel processing flow is executed: S601, dividing the sampling point set into blocks according to spatial grid units according to the pre-built spatial index structure; S602: Calling a parallel computing unit to synchronously perform grid positioning and physical quantity interpolation operations on sampling points in each block; S603, using a pipeline loading mechanism for multi-time step data, asynchronously reading the next time step data; S604: Aggregate and output the interpolation results of each block.

Images