A method for generating adaptive grid and mapping data based on tree data structure

Abstract

The application discloses a kind of FDTD adaptive grid generation and data mapping method based on tree data structure, through the adaptive equal-ratio refinement mechanism of tree structure, the "on-demand distribution" of spatial resolution is realized, while guaranteeing the high-fidelity modeling of fine structure, the number of grids is maximally reduced.In addition, by introducing the Morton coding technology, the topological linearization rearrangement is performed on the generated hierarchical non-uniform grid, ensuring the high aggregation of physically adjacent grids in physical memory address, eliminating the complex pointer communication overhead, and ensuring that the algorithm shows the ultimate operation flux and excellent parallel speedup ratio in multi-scale complex system simulation.

Description

Technical Field

[0001] This invention relates to the fields of computational electromagnetics and high-performance numerical computing, specifically to an FDTD adaptive mesh generation and data mapping method based on a tree-structured data structure. Background Technology

[0002] With the widespread application of Computational Electromagnetics (CEM) in engineering, the geometric complexity and electrical dimensions of the models to be simulated are increasing, exhibiting significant "multi-scale" and "irregular" characteristics. Among numerous numerical solution algorithms, the Finite Difference Time Domain (FDTD) method has become the mainstream tool for analyzing transient electromagnetic fields due to its ability to directly solve the curl equation, support parallel computing, and obtain broadband responses in a single simulation. However, traditional FDTD algorithms typically rely on a globally uniform Cartesian grid to discretize the space. Limited by the Nyquist sampling theorem and the requirements for geometric modeling accuracy, the grid step size of the global computational domain is forced to depend on the smallest geometric feature (such as tiny gaps or thin-layer structures) or the shortest wavelength in the model. For large-scale multi-scale computational problems containing local fine structures, this limitation leads to excessive densification of many non-critical regions, resulting in a waste of computational memory and time resources, making it difficult to meet the needs of engineering for rapid solutions to large-scale complex problems.

[0003] To address the challenges of cross-scale simulations (where the model simultaneously contains large and fine structures), numerous scholars have researched mesh generation techniques for FDTD (Focused Fiber-Codd Transformation). In 2017, the Chinese Journal of Electronics (45(12):2944-2948) published a novel sub-mesh technique based on the FDTD algorithm. This technique, based on a spatial filtering algorithm, improves the stability of the FDTD algorithm through spatial frequency domain filtering, resulting in a highly stable FDTD algorithm. However, when performing data interaction across scales or non-uniform regions, it is easy to generate jump accesses across memory regions and additional conditional branch operations, which can affect the spatial locality of memory data to some extent, leading to a decrease in cache hit rate under modern high-performance computing architectures.

[0004] Chinese invention patent application (publication number: CN115758492A) discloses an automatic generation method for FDTD conformal meshes based on 3D CAD. By dividing the CAD model into mesh lines, and then calculating the geometric and material information of each boundary cell through longitudinal ray scanning and transverse surface scanning, this method significantly improves the generation rate and reduces computational complexity. However, while conformal meshes can partially alleviate step approximation errors when dealing with models containing fine, sharp structures, they still strongly rely on fine meshes to accurately capture field singularities at points of curvature abrupt change; if the mesh is not fine enough, non-physical oscillations in the local field will compromise the overall simulation accuracy. Summary of the Invention

[0005] To address the issues of traditional FDTD algorithms, where uniform meshes lead to significant waste of computational memory and time resources due to the limitation of the global minimum mesh size when processing multi-scale models containing fine structures, and the problems arising from the significant disruption of data memory storage continuity caused by the introduction of unstructured meshes in existing local submesh or non-uniform mesh techniques, rendering traditional array-index-based direct indexing methods ineffective, this invention provides an FDTD adaptive mesh generation and data mapping method based on a tree-structured data structure. Through an adaptive proportional refinement mechanism using a tree structure, spatial resolution is allocated "on demand," minimizing the number of meshes while ensuring high-fidelity modeling of fine structures. Furthermore, by introducing Morton coding technology, topological linearization rearrangement is performed on the generated hierarchical non-uniform meshes, ensuring a high degree of clustering of physically adjacent meshes in physical memory addresses, eliminating cumbersome pointer communication overhead, and ensuring that the algorithm exhibits extreme computational throughput and excellent parallel speedup in multi-scale complex system simulations.

[0006] To solve the above problems, the present invention adopts the following technical solution: An FDTD adaptive grid generation and data mapping method based on a tree-structured data structure includes the following steps: Step 1: Import a two-dimensional electromagnetic target model with second-order differential properties. Divide the surface contour of the electromagnetic target into several closed geometric boundary segments and calculate the local differential geometric properties at each point of the geometric boundary segments. Based on the local differential geometric properties, construct a geometric feature index tree covering the entire computational domain to enable fast querying of the geometric curvature properties of any spatial region within the entire computational domain. Step 2: Construct a tree-like initial grid covering the entire computational domain; call the geometric feature index tree from Step 1 to perform topological intersection detection on each grid node, and perform tree split evaluation based on the detection results. Perform tree split on grid nodes that meet the split conditions; generate a non-uniform grid by recursively executing the topological intersection detection and tree split evaluation process, and extract terminal grid nodes that no longer have child nodes as the final leaf nodes participating in the computation. Step 3: Rearrange all the leaf nodes that ultimately participate in the calculation based on Morton coding topology linearization, and perform pointerless FDTD electromagnetic field component time step iteration based on the rearranged one-dimensional linear storage structure. During the iteration phase, direct offset access is performed using the mapping relationship.

[0007] Furthermore, the local differential geometric properties in step one include the tangent vector. Normal vector And local curvature used to reflect local geometric bending characteristics .

[0008] Furthermore, in step one, constructing a geometric feature index tree covering the entire computational domain involves storing all geometric boundary segments and their corresponding maximum local curvature values ​​into the geometric feature index tree. The implementation method is as follows: The system detects the spatial intersection of each geometric boundary segment with the bounding box of each tree node in the geometric feature index tree. If a geometric boundary segment passes through or is contained within the bounding box of a tree node, the maximum local curvature value corresponding to that segment is assigned to that tree node. If the bounding box of a tree node contains multiple geometric boundary segments, the maximum local curvature value among all intersecting geometric boundary segments is extracted and used as the unified curvature field value of that tree node, so as to realize the fast query of the geometric curvature characteristics of any spatial region in the entire computational domain.

[0009] Furthermore, the topological intersection detection specifically involves: calling the geometric feature index tree from step one to detect whether the geometric coverage of the grid node to be tested includes the medium boundary; if it does not include the medium boundary and does not meet the CFL stability condition, it needs to be forcibly split until the CFL stability condition is met; if it includes the medium boundary and the level reaches the preset maximum refinement level, the splitting condition is not met, and if the level does not reach the preset maximum refinement level, a tree-like splitting evaluation is performed.

[0010] Furthermore, the specific method for detecting whether the geometric coverage of the grid node under test includes the medium boundary by calling the geometric feature index tree in step one is as follows: the geometric coverage of the grid node under test is constructed as a grid axial bounding box; the geometric feature index tree described in step one is traversed from top to bottom, and the spatial bounding boxes of the grid axial bounding boxes and the nodes at each level in the index tree are calculated to see if they overlap: if they do not overlap, the branch of the index tree is directly removed and determined to be non-intersecting; if they overlap, the tree is recursively traversed downwards along the branch until the leaf node of the index tree containing the specific parameterized geometric boundary segment is located, and then the precise mathematical intersection calculation between the grid node under test and the specific geometric boundary segment is performed to determine whether the grid node under test includes the medium boundary.

[0011] Furthermore, the expression for the CFL stability condition is as follows:

[0012] in, For time step, and These are the spatial step sizes in the x and y directions, respectively. Let x be the speed of light in a vacuum. In the FDTD method, the x and y step sizes are the same, and the stability condition simplifies to: ,in It is the feature size of the current mesh.

[0013] Furthermore, the tree-like splitting evaluation specifically involves calculating the maximum theoretical geometric error generated by approximating the curve boundary with a straight line segment based on the mesh size and local geometric properties corresponding to the current mesh node to be tested; if the maximum theoretical geometric error is greater than a preset error threshold, the mesh node to be tested is determined to meet the splitting condition, and tree-like splitting is performed on the current mesh node, splitting the current mesh node into four segments with a size of... The child nodes, where L is the feature size of the current mesh.

[0014] Furthermore, the maximum theoretical geometric error is expressed as the maximum chord height error. It is represented as: ,in The feature size of the current grid node to be tested. It represents the maximum curvature of the curve within the area enclosed by the currently tested grid node.

[0015] Furthermore, the Morton coding-based topological linearization rearrangement in step three specifically involves: extracting the two-dimensional row and column logical coordinates of the final leaf nodes participating in the calculation, and converting them into equal-length binary sequences; interleaving and combining the binary bits of the row and column coordinates in descending order to generate Z-curve binary codes; converting the generated Z-curve binary codes into decimal values, which are the Morton codes corresponding to the physical spatial locations; performing topological linearization rearrangement on all leaf nodes participating in the calculation according to the Morton codes from smallest to largest, and storing them in a one-dimensional linear array in the rearrangement order.

[0016] Furthermore, the pointerless FDTD electromagnetic field component time-stepping iterative calculation in step three specifically involves: before the FDTD time-stepping iteration, using tree topology relationships to find the actual storage locations of the current mesh and its physically adjacent meshes, calculating the difference and storing it as an offset index in the current mesh attributes; for cross-scale coarse and fine mesh boundaries, multiple pointers are assigned to the coarse mesh to record the offset indices of adjacent fine meshes; during the iteration phase, direct offset access is performed using mapping relationships, and combined with direct offset index logic to complete the cross-scale boundary time stepping, thus achieving... The field value calculation has a time complexity; the mapping relationship is: NeighborIndex = CurrentIndex + offset; Where NeighborIndex is the memory space index of the neighboring grid, CurrentIndex is the memory space index of the current grid node, and offset is the offset required to access the memory location of the neighboring grid from the current memory location.

[0017] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention introduces a quadtree data structure as the logical basis for spatial partitioning, constructing an adaptive mesh framework with highly adaptive resolution. The hierarchical depth of tree nodes is dynamically adjusted based on local geometric features, achieving decoupling between the mesh discrete scale and the minimum feature size of the geometric model from the underlying data logic. For multi-scale electromagnetic models containing large areas of uniform background medium and local fine structures (such as dielectric waveguide structures), this "on-demand allocation" of mesh density breaks the requirements of CFL conditions and modeling accuracy on the global mesh step size in the traditional FDTD algorithm. While ensuring high-fidelity modeling of local fine structures, it significantly reduces the total number of meshes, thereby significantly saving memory usage and computation time.

[0018] 2. When solving Maxwell's curl equations using the FDTD algorithm, the time extrapolation of the electric or magnetic field components at any point in space is highly dependent on the spatial differential coupling of the electromagnetic field components in their physical neighborhood. This invention introduces Morton coding technology to perform topological linearization on non-uniform meshes, ensuring that physically adjacent mesh cells in the two-dimensional electromagnetic target model to be simulated remain highly continuous in the computer's physical memory. This optimization significantly improves the cache hit rate of CPU / GPU during numerical computation.

[0019] 3. By solidifying the topological adjacency relationships between grids, complex tree traversal and hash lookup are transformed into direct memory access based on offsets, thus achieving... Extremely fast addressing with minimal time complexity. This data structure-level optimization eliminates complex conditional branches and pointer jumps, greatly reducing communication and synchronization overhead, and ensuring that the algorithm has excellent linear speedup and scalability in ultra-large-scale parallel computing environments. Attached Figure Description

[0020] Figure 1 This is an overall flowchart of the FDTD adaptive grid generation and data mapping method based on tree data structure of the present invention.

[0021] Figure 2 This is a schematic diagram illustrating the generation of Morton codes by interleaving the row and column values ​​of the grid coordinates according to the present invention.

[0022] Figure 3 This is a schematic diagram of the tree structure and grid coordinate mapping of the present invention.

[0023] Figure 4 This is a schematic diagram of the adjacency offset index of the tree structure cross-scale grid of the present invention.

[0024] Figure 5 This is a diagram illustrating the geometric model of a Y-type dielectric waveguide.

[0025] Figure 6(a) is an adaptive mesh model generated after spatial partitioning of the Y-type dielectric waveguide using the method of the present invention.

[0026] Figure 6(b) is a magnified view of the local mesh of the bifurcation region of the Y-type dielectric waveguide in Figure 6(a).

[0027] Figure 6(c) is a magnified view of the local mesh at the dielectric boundary of the Y-type dielectric waveguide in Figure 6(a).

[0028] Figure 7 This is a diagram illustrating the geometric model of a photonic crystal waveguide.

[0029] Figure 8(a) is an adaptive mesh model generated after spatial partitioning of the photonic crystal waveguide using the method of the present invention.

[0030] Figure 8(b) is a partial enlarged view of the periodic dielectric pillar structure of the photonic crystal waveguide in Figure 8(a); Figure 8(c) is a magnified view of the defect channel in the photonic crystal waveguide in Figure 8(a); Figure 8(d) is a magnified view of the boundary grid of the photonic crystal waveguide in Figure 8(a). Detailed Implementation

[0031] The following detailed description of the FDTD adaptive grid generation and data mapping method based on tree data structure, in conjunction with the accompanying drawings and specific embodiments, provides a comprehensive overview of the present invention.

[0032] like Figure 1As shown, the present invention provides an FDTD adaptive grid generation and data mapping method based on a tree-structured data structure, comprising the following steps: Step 1: Import the two-dimensional geometric model of the target to be simulated, which has second-order differential properties, and define the boundary contour of the target to be simulated as a set of closed piecewise parametric curves: ,in It is a vector function of the parametric curve, which gives the corresponding parameters. At that time, the two-dimensional coordinate position of that point on the curve Calculate the tangent vector at any point on the boundary. Normal vector And the local curvature used to reflect the degree of geometric curvature at that location. Among them, the tangent vector Normal vector and local curvature The calculation formulas are as follows:

[0033]

[0034] in, For curve functions with respect to parameters The first-order partial derivative, Let be the unit normal vector of the plane containing the two-dimensional curve, usually taken as (-1, 0, 1). For curve functions with respect to parameters The second-order partial derivative, Let be the magnitude of the vector.

[0035] A feature index tree for the geometric model is constructed to detect the spatial intersection of each geometric boundary segment with the bounding box of each tree node. If a geometric boundary segment passes through or is contained within the bounding box of a tree node, the local maximum curvature value corresponding to the segment is assigned to the tree node. If the bounding box of a tree node contains multiple geometric boundary segments, the maximum curvature value among all intersecting segments is extracted as the unified curvature field value of the tree node, which is used to determine whether the current position needs to be refined during the refinement process.

[0036] Step 2: Perform spatial adaptive partitioning based on a quadtree data structure. Construct an initial quadtree-shaped mesh covering the entire computational domain. The initial size of the mesh nodes is preset according to the boundary range of the computational domain. The specific preset rule is as follows: calculate the global physical bounding box of the two-dimensional electromagnetic target model to be simulated, and extend it outward with the thickness of the additional absorbing boundary conditions and the width of the vacuum buffer layer; take the maximum length and width of the total computational domain after extension as the initial reference size.

[0037] The minimum mesh size limit is determined based on the simulation accuracy requirements: the highest cutoff frequency of the injected electromagnetic excitation source and the maximum relative permittivity in the entire computational domain are extracted, and the shortest wavelength of the electromagnetic wave in the corresponding medium is calculated; according to the requirements of satisfying the numerical dispersion tolerance (such as taking 1 / 10 to 1 / 20 of the shortest wavelength), the minimum spatial step size to ensure simulation accuracy is calculated, and it is used as the minimum size lower limit of the quadtree mesh in the adaptive splitting process.

[0038] Starting from the root node (Level=0) that covers the entire node, execute the recursive split (RecursiveSplit(Node)) function; the specific steps are as follows: 1. Topological Intersection Detection: The geometric feature index tree from step one is invoked to query whether the geometric coverage of the current mesh node contains a medium boundary (i.e., the physical contour line where electromagnetic parameters change abruptly). Specifically, the geometric coverage of the current mesh node is constructed as a mesh axial bounding box (AABB). The geometric feature index tree from step one is traversed from top to bottom to calculate whether the mesh axial bounding box overlaps with the spatial bounding boxes of nodes at each level in the index tree. If there is no overlap, the index tree branch is directly removed, indicating a non-intersection. If there is an overlap, the tree branch is recursively traversed downwards until a leaf node containing a specific parameterized geometric boundary segment is located. Then, a precise mathematical intersection calculation is performed between the mesh node to be tested and the specific geometric boundary segment to determine whether the mesh node contains a medium boundary. If no medium boundary is contained and the CFL stability condition is met, the splitting stops. If the CFL stability condition is not met, a forced splitting is required until the CFL stability condition is met. The expression for the CFL stability condition is:

[0039] in, For time step, and Let x and y be the spatial step sizes in the x and y directions, respectively. In the FDTD method, the x and y step sizes are the same, so the CFL stability condition simplifies to:

[0040] in It is the feature size of the current mesh.

[0041] If the current grid node's geometric coverage includes a medium boundary and the refinement level reaches the preset maximum refinement level, then the splitting condition is not met, and splitting stops; if the refinement level does not reach the preset maximum refinement level, then a tree-like splitting evaluation is performed.

[0042] 2. Tree-like splitting evaluation: Extract the maximum curvature value of the curved boundary within the region enclosed by the grid. Substitute the values ​​into the following formula to calculate the maximum chord height error of the curve approximating its straight line segment. ; , in This represents the feature size of the current mesh node. If the maximum chord height error of the curve region contained in the current mesh is greater than a preset threshold, and the current refinement level has not reached the preset maximum refinement level, it indicates that the current mesh is too large and cannot fit the curvature at the required accuracy, requiring splitting. If the maximum chord height error of the curve region contained in the current mesh is less than or equal to the preset threshold, it indicates that the mesh meets the computational requirements and does not need to be split.

[0043] 3. Split Execution: Divide the node to be split into four equal parts of size 1. The child nodes are then recursively split, and the RecursiveSplit(Node) function is called again. As the mesh size is halved, the chord height error will... The speed converges rapidly.

[0044] 4. Solidify the final mesh: Solidify the mesh results after splitting to meet the CFL stability requirements and chord height error judgment to generate a non-uniform mesh; define the terminal mesh nodes that no longer have child nodes in the non-uniform mesh topology as the leaf nodes that finally participate in the calculation, as the smallest spatial discrete unit that actually carries electromagnetic field components and participates in FDTD time step iteration, and extract all the leaf nodes that finally participate in the calculation in the entire computational domain.

[0045] Step 3: Organize and map the generated leaf node meshes involved in the computation. This data organization and mapping is implemented using topological linearization rearrangement based on Morton coding. The specific rules are as follows: traverse and extract all the final leaf nodes involved in the computation generated in Step 2. Due to the inconsistent adaptive mesh size, first project the center coordinates of all leaf nodes involved in the computation onto a global refined mesh to obtain their corresponding two-dimensional logical integer row and column indices. .Will and Convert to a 32-bit binary number. binary bits and The binary bits are interleaved and combined bit by bit from most significant bit to least significant bit to generate the binary code for the Z-curve. If... , The generated 64-bit Z-curve binary code The generated Z-curve binary code is converted into a decimal value, which is the Morton code corresponding to the physical spatial location. All leaf nodes involved in the calculation are then topologically linearized and rearranged in ascending order of this code, and stored in a one-dimensional linear storage structure.

[0046] Subsequently, the relative memory offsets of adjacent grids are calculated and stored. The Morton codes and corresponding spatial hierarchy depth information of each leaf node in the one-dimensional linear storage structure are extracted. A spatial topology model of the two-dimensional adaptive non-uniform grid is reconstructed through decoding mapping to visualize the grid division results and verify the accuracy of capturing the geometric features of the two-dimensional electromagnetic target model to be simulated. Then, based on the one-dimensional linear storage structure, the memory offsets between adjacent leaf nodes are calculated using the Morton-coded mapping relationship, and pointerless FDTD electromagnetic field component time-stepping iterative deduction is performed. During the iterative deduction process or after meeting the preset iteration termination condition, the electromagnetic field time-domain data at the set observation nodes are extracted and post-processed for analysis to obtain the electromagnetic characteristic simulation results of the target model.

[0047] like Figure 2 The diagram illustrates the process of generating Morton codes by interleaving the row and column values ​​of grid coordinates. It details the underlying mapping process for converting the logical space coordinates of a two-dimensional non-uniform discrete grid to a one-dimensional linear index. Specifically, for any grid within the computational domain, the integer indices of its logical row and column coordinates are first obtained, and these two decimal coordinate values ​​are converted into equal-length binary sequences. Then, following the order from most significant bit to least significant bit, the Morton codes are... coordinates and The binary data bits of the coordinates are interleaved and combined bit by bit to generate a new binary code that is twice as long. This binary code, after being converted into a decimal value, is the Morton code corresponding to that grid.

[0048] like Figure 3 The diagram illustrates the mapping between a tree structure and grid coordinates. This binary bit cross-mapping mechanism enables the dimensionality reduction transformation from two-dimensional to one-dimensional space. Its core technological advantage lies in the fact that grids that were originally adjacent in two-dimensional physical space retain a high degree of consistency and proximity in their numerical values ​​after being converted to one-dimensional Morton encoding. This maximizes the preservation of locality in physical space, thereby significantly eliminating memory jump accesses.

[0049] like Figure 4The diagram illustrates the adjacency offset index of a tree-structured cross-scale grid. This diagram visually demonstrates how this invention replaces the traditional, complex cross-scale grid search logic with integer offset operations in a one-dimensional memory array. All leaf node grids ultimately participating in the calculation are linearized using Morton coding and then uniformly stored in a one-dimensional linear array Array 1D. To achieve fast addressing, during the preprocessing stage before the FDTD time stepping officially begins, the system pre-assigns and fills in neighbor relationships for each grid. Specifically, the topological relationships of the quadtree are used to accurately find the adjacent grids in physical space, and then the index difference between these two adjacent grids in the one-dimensional linear array (Offset = n or Offset = m) is calculated and directly stored as the "offset" in the attribute field of the current grid.

[0050] Before iteration, traverse each grid in Array 1D, using a quadtree to search for its physical neighbors. Assuming the right neighbor is found, the required offset is... for: , in, This represents the absolute index of the neighboring node's memory location in the one-dimensional array. This is the absolute index of the current node's memory location in the one-dimensional array.

[0051] against Figure 4 The cross-scale boundary (i.e., the region where coarse and fine grids are adjacent) shown in the diagram is addressed by assigning multiple pointers to the coarse grid boundary, each pointing to one of the two adjacent fine grid indices. The offset of the coarse grid relative to each of the two adjacent fine grids is then calculated. After establishing this adjacency relationship, subsequent FDTD iterative calculations are completely independent of complex tree structures or pointer lists. When updating the electromagnetic field, the FDTD update equations do not require complex conditional judgments and spatiotemporal coordinate differences at the coarse-fine grid boundary. A simple integer offset read operation using the current grid's offset is sufficient to directly locate and obtain the field data of the surrounding grids. This direct memory offset-based access mechanism completely eliminates the search overhead during iteration, achieving extremely high computational efficiency consistent with uniform grids.

[0052] The present invention will be further described in detail below with reference to specific embodiments.

[0053] Example 1: Adaptive Mesh Generation of Y-type Dielectric Waveguide FDTD This embodiment takes a two-dimensional Y-type dielectric waveguide structure as the research object to illustrate the adaptive mesh partitioning based on tree data structure described in this invention.

[0054] The two-dimensional Y-type dielectric waveguide model constructed in this embodiment is defined in a rectangular coordinate system. Below, the calculation area is: , The background medium is air, and its relative permittivity is... Y-type dielectric waveguides are composed of high refractive index dielectrics with a relative permittivity of [missing information]. The input main waveguide width is bifurcation angle The main waveguide extends along the positive X-axis and splits into two symmetrical branches at the bifurcation point.

[0055] This example first constructs a first-order uniform background mesh, dividing the entire computational domain into... , That is, the background grid size is: .

[0056] Based on the quadtree data structure, spatial adaptive partitioning is performed. The partitioning criteria adopted in this embodiment strictly follow the mesh refinement requirements in step two of the above technical solution.

[0057] The maximum depth of the refinement level is set as follows: For each quadtree node, determine whether the refinement condition is met; if it is, divide the current unit into four sub-units; each sub-unit is added to the tree structure as a new child node. The above judgment is recursively performed on the child nodes. If the refinement condition is not met or the maximum level is reached, then the node is defined as a leaf node participating in the calculation.

[0058] Through the above process, a high-density fine mesh is formed in the bifurcation region of the Y-type dielectric waveguide and near the dielectric boundary, while a coarse mesh structure is maintained in the air homogeneous region.

[0059] Each quadtree node contains the following data fields: node level, spatial coordinate center. Unit size , node pointer child[4], Morton code value Z-Code. In order to improve data locality and storage efficiency, this embodiment sorts all leaf nodes involved in the calculation by Morton code. For a node in the Lth layer, its two-dimensional integer index is: ,pass Morton codes are generated using a cross-arrangement method. Each code uniquely corresponds to a node's spatial location. The tree structure is then sorted according to a preorder traversal, and numbered in ascending order. The sorted Morton codes are then mapped to the leaf nodes used in the computation after the traversal.

[0060] Traditional FDTD adaptive meshes require maintaining cross-array index mapping relationships at the boundary between coarse and fine meshes. This embodiment eliminates cross-array index mapping in the following way: 1. All leaf nodes involved in the calculation are stored in a single linear array; 2. Adjacent nodes maintain spatial locality in the Morton-coded sequence; 3. Construct an adjacency offset table: ,Right now: , This transforms the traditional multi-array mapping relationship into integer offset operations. This method significantly reduces the time complexity of index lookup.

[0061] Figure 5 This is a schematic diagram of the geometric model of the Y-shaped dielectric waveguide, the electromagnetic target to be simulated, used in a specific embodiment of the present invention. The model is defined in a two-dimensional Cartesian coordinate system. This geometric model allows for the extraction of the local curvature of the boundaries and the abrupt changes in the dielectric material, thus providing the basic input parameters for the subsequent adaptive refinement of the tree-like mesh.

[0062] Figures 6(a), 6(b), and 6(c) illustrate the application of the adaptive grid generation method based on tree data structure described in this invention to... Figure 5 The image shows the final mesh distribution generated after spatial partitioning of the Y-shaped dielectric waveguide. Figure 6(a) shows the global mesh distribution effect in the computational domain, within the geometric boundary of the original structure and the background region (air region) where material properties are uniformly distributed. The algorithm of this invention confirms that there is no medium boundary crossing in the region through topological intersection detection, and the preset theoretical error of the chord height is also considered. Much smaller than the maximum permissible geometric error Therefore, the algorithm did not trigger hierarchical splitting and automatically retained the side length. The initial mesh structure. This processing method effectively avoids oversampling in non-critical areas, thereby significantly reducing the amount of redundant mesh and saving computational memory and time resources.

[0063] Figures 6(b) and 6(c) show magnified local mesh views of regions with fine geometry, curvature edges, and abrupt changes in material properties (such as the bifurcation apex of a Y-type dielectric waveguide), respectively. To achieve high-fidelity fitting of complex boundaries and ensure simulation accuracy, this embodiment sets a strict recursive refinement criterion: setting the maximum allowable geometric error for abrupt change regions. And will maximize the refinement level With a 7-layer configuration, the local geometric curvature was analytically extracted in the region near the waveguide bifurcation tip. (Corresponding radius of curvature is) If the side length is directly used... The initial background mesh approximates the sharp corner, and by substituting it into the chord height theoretical formula, we can obtain... The error value far exceeds the set threshold. The algorithm then triggers recursive splitting. When the mesh recursively splits to the 6th layer (i.e., the mesh size...), the algorithm continues this process. When this happens, the theoretical chord height error is recalculated. It is still greater than the threshold. Therefore, the algorithm of this invention triggers a split at layer 7, further reducing the local mesh size to... At this point, the theoretical chord height error drops to If the constraints are successfully met, the node in that region is fixed as a leaf node participating in the calculation and the splitting stops.

[0064] The number of grids at different refinement levels was statistically analyzed, and the results are shown in Table 1.

[0065] Table 1 Comparison of mesh counts for Y-type dielectric waveguide models with uniform mesh and different refinement levels.

[0066] As shown in Table 1, while ensuring the same local spatial resolution at the fine structure, the adaptive meshing mechanism of this invention can significantly eliminate redundant meshes in a uniform air background. With the adaptive refinement level... As the refinement level increases, the mesh compression effect becomes increasingly significant. When the maximum refinement level is 7, the total number of meshes generated by this invention is only 236,929, accounting for only 5.48% of the total number of traditional global uniform meshes (4,323,328). This extremely high data compression ratio intuitively demonstrates that the method of this invention can greatly save computer memory resources when processing multi-scale structures with large areas of background free space.

[0067] Example 2: Adaptive Mesh Partitioning of Photonic Crystal Waveguide FDTD This embodiment uses a two-dimensional dielectric cylindrical photonic crystal waveguide as an example to illustrate the adaptive mesh partitioning, Morton coding linearization storage, and offset index construction method based on a tree-like data structure described in this invention, focusing on demonstrating the multi-region hierarchical refinement capability. For example... Figure 7 The figure shown is a geometric model diagram of the photonic crystal waveguide of this invention: In this embodiment, a two-dimensional dielectric cylindrical photonic crystal model is established, and the computational region is defined as follows: , The background medium is air. The photonic crystal consists of periodically arranged high-refractive-index cylinders with a relative permittivity of... The crystal lattice type is a two-dimensional square lattice. Lattice constant. The radius of the medium column is A defective waveguide is formed by removing a column of dielectric pillars along the propagation direction (x-direction). The calculation employs a two-dimensional TE mode FDTD method, dividing the entire computational domain into a first-order uniform background mesh. ,

[0068] Get the background unit size: .

[0069] This example takes into account the characteristics of dielectric cylindrical photonic crystal waveguide structures and sets the following multi-region refinement strategy: 1. All units falling within the overall region of the photonic crystal lattice are refined to at least level 2 to ensure the analytical accuracy of the periodic structure.

[0070] 2. If the element intersects with the circular boundary of the medium column, recursively refine it until the maximum refinement level is reached (level=5).

[0071] 3. If a cell contains interfaces between materials with different dielectric constants, further refinement is performed. This strategy ensures that the region of abrupt electromagnetic field change between the high-refractive-index dielectric pillar and air is finely resolved. No refinement is performed in uniform air regions far from the photonic crystal region. For each mesh, the mesh refinement operation described in step two of the technical solution is executed.

[0072] Figures 8(a) to 8(d) show the global adaptive mesh of the photonic crystal waveguide model, and local magnified views of the periodic dielectric pillars and defect channel regions and boundary regions, respectively. Since the photonic crystal in this embodiment is composed of periodically arranged high-refractive-index dielectric pillars (… Embedded in air background ( The electromagnetic field at its circular boundary exhibits strong locality and spatial gradient. To ensure high-fidelity simulation of wave propagation at the photonic bandgap and defects, this embodiment sets a maximum permissible geometric error threshold. And limit the maximum refinement level to At the circular physical boundary of the dielectric column, its geometric radius... The corresponding constant local curvature If side length is used The initial first-level coarse mesh directly approximates the circular arc boundary, with a theoretical chord height error. Far exceeding the threshold Furthermore, there is a significant jump in the dielectric constant of the material at this location, which triggers hierarchical recursive splitting in the algorithm.

[0073] When the mesh splits down to the 4th layer (mesh size) When this happens, the theoretical chord height error is recalculated. Since this value is still greater than the preset threshold, the fifth layer of refinement is triggered. After the fifth layer of refinement is completed, the boundary mesh size... This ensures that the geometric approximation error is completely converged within the tolerance, and the node at that point is then fixed as a leaf node participating in the calculation.

[0074] The method described in this embodiment is used to perform mesh generation on the photonic crystal waveguide. Table 2 shows a comparison of the number of meshes in the photonic crystal waveguide model under traditional uniform mesh and adaptive mesh at different refinement levels according to this invention. As shown in Table 2, for complex structures containing a large number of tiny periodic dielectric pillars, traditional uniform mesh is forced to refine the entire domain to meet the fitting accuracy of the boundary of a single dielectric pillar, resulting in a total of 1,048,576 meshes. However, the adaptive mesh generation strategy of this invention, due to the implementation of multi-region hierarchical refinement, achieves a higher number of meshes at the maximum refinement level. At that time, the total number of grids was significantly reduced to 176,884, and the number of grids decreased to 16.87%. This further verifies that the present invention can still maintain excellent grid reduction capability when dealing with complex electromagnetic models with high density and multi-scale fine feature superposition.

[0075] Table 2 Comparison of mesh counts in the photonic crystal waveguide model with uniform mesh and mesh counts at different refinement levels.

[0076] Comparative data analysis of two specific embodiments shows that, while maintaining the same geometric boundary fitting accuracy, the method of this invention reduces the total mesh size of the Y-type waveguide model to 5.48% of the original and the total mesh size of the photonic crystal waveguide model to 16.87%. This order-of-magnitude reduction in mesh size significantly reduces the resource consumption of multi-scale models in terms of computer memory capacity and iteration time. Existing local subgrid or non-uniform mesh techniques, while reducing the total mesh size, often lead to the destruction of physical continuity of memory, increased pointer addressing overhead, and increased complexity of boundary data interaction. This invention, by introducing Morton coding for topological linearization and adjacency offset indexing, achieves the aforementioned mesh reduction rate of 83.13% to 94.52% while mapping the unstructured mesh across scales into a one-dimensional continuous array with extremely high spatial locality. In the FDTD iterative calculation stage, this method eliminates the overhead of non-contiguous memory jump access and pointer traversal, realizing direct offset addressing based on integer addition. This optimization of the underlying mechanism effectively maintains the cache hit rate under modern processor architecture while significantly reducing the total computational load, thus providing reliable data structure support for the efficient deployment and linear expansion of this algorithm in a large-scale parallel computing environment.

Claims

1. A method for adaptive grid generation and data mapping based on tree-structured FDTD, characterized in that, Includes the following steps: Step 1: Import a two-dimensional electromagnetic target model with second-order differential properties. Divide the surface contour of the electromagnetic target into several closed geometric boundary segments and calculate the local differential geometric properties at each point of the geometric boundary segments. Based on the local differential geometric properties, construct a geometric feature index tree covering the entire computational domain to enable fast querying of the geometric curvature properties of any spatial region within the entire computational domain. Step 2: Construct a tree-like initial grid covering the entire computational domain; call the geometric feature index tree from Step 1 to perform topological intersection detection on each grid node, and perform tree split evaluation based on the detection results. Perform tree split on grid nodes that meet the split conditions; generate a non-uniform grid by recursively executing the topological intersection detection and tree split evaluation process, and extract terminal grid nodes that no longer have child nodes as the final leaf nodes participating in the computation. Step 3: Rearrange all the leaf nodes that ultimately participate in the calculation based on Morton coding topology linearization, and perform pointerless FDTD electromagnetic field component time step iteration based on the rearranged one-dimensional linear storage structure. During the iteration phase, direct offset access is performed using the mapping relationship.

2. The FDTD adaptive mesh generation and data mapping method based on tree data structure according to claim 1, characterized in that, The local differential geometric properties in step one include the tangent vector. Normal vector And local curvature used to reflect local geometric bending characteristics .

3. The FDTD adaptive mesh generation and data mapping method based on tree data structure according to claim 2, characterized in that, In step one, constructing a geometric feature index tree covering the entire computational domain involves storing all geometric boundary segments and their corresponding maximum local curvature values ​​into the geometric feature index tree. The implementation method is as follows: The system detects the spatial intersection of each geometric boundary segment with the bounding box of each tree node in the geometric feature index tree. If a geometric boundary segment passes through or is contained within the bounding box of a tree node, the maximum local curvature value corresponding to that segment is assigned to that tree node. If the bounding box of a tree node contains multiple geometric boundary segments, the maximum local curvature value among all intersecting geometric boundary segments is extracted and used as the unified curvature field value of that tree node, so as to realize the fast query of the geometric curvature characteristics of any spatial region in the entire computational domain.

4. The FDTD adaptive mesh generation and data mapping method based on tree data structure according to claim 1, wherein, The topological intersection detection specifically involves: calling the geometric feature index tree from step one to detect whether the geometric coverage of the grid node to be tested includes the medium boundary; if it does not include the medium boundary and does not meet the CFL stability condition, it needs to be forcibly split until the CFL stability condition is met; if it includes the medium boundary and the level reaches the preset maximum refinement level, the splitting condition is not met, and if the level does not reach the preset maximum refinement level, a tree-like splitting evaluation is performed.

5. The FDTD adaptive mesh generation and data mapping method based on tree data structure according to claim 4, characterized in that, The specific method for calling the geometric feature index tree in step one to detect whether the geometric coverage of the grid node to be tested includes the medium boundary is to construct the geometric coverage of the grid node to be tested as a grid axis bounding box. The geometric feature index tree described in step one is traversed from top to bottom to calculate whether the bounding box of the mesh axis overlaps with the spatial bounding boxes of each level node in the index tree. If they do not overlap, the branch of the index tree is directly removed and determined to be non-intersecting. If they overlap, the tree branch is recursively traversed downwards until the leaf node of the index tree containing the specific parameterized geometric boundary segment is located. Then, the precise mathematical intersection calculation between the mesh node to be tested and the specific geometric boundary segment is performed to determine whether the mesh node to be tested contains a medium boundary.

6. The FDTD adaptive mesh generation and data mapping method based on tree data structure according to claim 4, characterized in that, The expression for the CFL stability condition is: in, For time step, and These are the spatial step sizes in the x and y directions, respectively. Let x be the speed of light in a vacuum. In the FDTD method, the x and y step sizes are the same, and the stability condition simplifies to: ,in It is the feature size of the current mesh.

7. The FDTD adaptive mesh generation and data mapping method based on tree data structure according to any one of claims 1-6, characterized in that, The tree-like splitting evaluation specifically involves calculating the maximum theoretical geometric error generated by approximating the curve boundary with a straight line segment based on the mesh size and local geometric properties of the current mesh node to be tested. If the maximum theoretical geometric error is greater than a preset error threshold, the mesh node to be tested is determined to meet the splitting condition, and tree-like splitting is performed on the current mesh node, splitting it into four segments with a size of... The child nodes, where L is the feature size of the current mesh.

8. The FDTD adaptive mesh generation and data mapping method based on tree data structure according to claim 7, characterized in that, The maximum theoretical geometric error is the maximum chord height error. It is represented as: ,in The feature size of the current grid node to be tested. It represents the maximum curvature of the curve within the area enclosed by the currently tested grid node.

9. The tree data structure based FDTD adaptive mesh generation and data mapping method of claim 1, wherein, The Morton coding-based topological linearization rearrangement in step three specifically involves: extracting the two-dimensional row and column logical coordinates of the final leaf nodes participating in the calculation, and converting them into equal-length binary sequences; interleaving and combining the binary bits of the row and column coordinates in descending order to generate the binary code of the Z-curve; converting the generated binary code of the Z-curve into a decimal value, which is the Morton code corresponding to the physical spatial location; performing topological linearization rearrangement on all leaf nodes participating in the calculation according to the Morton code from smallest to largest, and storing them in a one-dimensional linear array in the rearrangement order.

10. The FDTD adaptive grid generation and data mapping method based on tree data structure according to claim 9, characterized in that, The pointerless FDTD electromagnetic field component time step iteration calculation in step three is as follows: before the FDTD time step iteration, the actual storage location of the current grid and its physically adjacent grids is found by using the tree topology relationship, the difference is calculated and stored in the current grid attribute as an offset index; For the cross-scale coarse-fine grid boundary, multiple pointers are assigned to the coarse grid to record the offset index of the adjacent fine grid; in the iteration stage, direct offset access is performed using the mapping relationship, and cross-scale boundary time stepping is completed combined with the direct offset index logic, realizing field value calculation of time complexity; wherein the mapping relationship is: NeighborIndex = CurrentIndex + offset; Where NeighborIndex is the memory space index of the neighboring grid, CurrentIndex is the memory space index of the current grid node, and offset is the offset required to access the memory location of the neighboring grid from the current memory location.

Images