Method and system for efficiently constructing Jacobian matrix in Newton iteration

Abstract

The invention discloses a method and a system for efficiently constructing a Jacobian matrix in Newton iteration. The method comprises the following steps of: constructing an (N * M) * (N * M) Jacobian matrix by adopting a CSR sparse format; pre-constructing an edge residual contribution expression, and recording a vertex variable, a material parameter and a calculation area of an edge ij by using a binary tree to form a parameterized calculation template; when multiple threads traverse grid vertexes, each thread only copies an expression template in a lightweight order once to generate a local instance, the contribution amount of an edge to a residual error and a related variable partial derivative are synchronously calculated through reverse automatic differential, and meanwhile, the symmetry of the partial derivative of an equation is utilized and is only ilt; if j, executing calculation, jgt; when i is greater than i, skipping to reduce the calculation amount, adjusting the calculation area according to the residual error and material correlation, and finally filling the matrix in parallel. According to the method, the automation advantage of reverse automatic differentiation and the high performance of near-manual coding are both considered, the cost is reduced, the expansibility of simulation software is enhanced, simulation time consumption is shortened by utilizing multi-core resources, and the method is suitable for scenes needing Newton iteration solution, such as device numerical simulation.

Description

Technical Field

[0001] This invention belongs to the field of computer-aided design technology, specifically relating to an efficient method and system for constructing Jacobian matrices in Newton iteration. Background Technology

[0002] In the field of computer-aided design for semiconductor technology, numerical simulation of semiconductor devices is a core method for analyzing the internal physical field distribution (such as electric potential and carrier concentration) and electrical performance (such as current-voltage characteristics) of devices. Its core lies in solving a coupled set of nonlinear partial differential equations consisting of the Poisson equation, the electron continuity equation, and the hole continuity equation. Among these, the finite volume method, due to its clear physical meaning and good conservation properties, is widely used to spatially discretize the above equations, resulting in a large-scale nonlinear algebraic equation system.

[0003] For this large-scale nonlinear algebraic equation system, the industry typically uses the Newton-Raphson iteration method for solution. Each iteration of the Newton-Raphson method requires calculating the residual vector of the equation system and the Jacobian matrix of the residual vector with respect to all solution variables. As a key support for Newton-Raphson iteration, the accuracy and computational efficiency of the Jacobian matrix directly determine the convergence, stability, and computational speed of the entire simulation process. Therefore, the efficient construction of the Jacobian matrix has become one of the core technical challenges in the field of semiconductor device simulation.

[0004] Currently, all existing methods for constructing the Jacobian matrix have significant shortcomings and are difficult to meet the needs of practical simulations: 1. Manual differentiation coding method: Technical personnel need to manually derive the complex analytical derivative formula of the residual with respect to the variable and code it. Although the computational efficiency is relatively high, when faced with residual equations containing complex physical models, the differentiation process is cumbersome and prone to errors. At the same time, the code is difficult to maintain in the later stage. When adding a new physical model, it needs to be re-derived and coded, and the scalability is extremely poor. 2. Numerical difference method: The derivative is approximated by applying a small perturbation to the variable. The logic is simple, but the residual needs to be calculated repeatedly to obtain the numerical difference before and after the perturbation. The computational cost is high, and the choice of perturbation amplitude directly affects the accuracy of the derivative. The accuracy control is difficult and it is easy to cause slow convergence or even failure of Newton iteration. 3. Forward Automatic Differentiation Method: By using operator overloading to propagate derivatives synchronously during the forward propagation residual calculation, although it is more efficient than the numerical difference method, the computational complexity increases linearly with the number of variables when calculating the derivative of a single residual with respect to multiple variables. In the case of multi-variable semiconductor device simulation scenarios, the efficiency still cannot meet the requirements. 4. Traditional Inverse Automatic Differentiation: As a core technology for efficiently calculating scalar gradients with respect to multiple parameters in deep learning, it theoretically meets the requirement of "differentiating a single residual with respect to multiple variables" in the Jacobian matrix. However, the traditional implementation requires dynamically constructing and traversing the computation graph (expression tree) for each calculation. In semiconductor device simulation, millions or even billions of grid cells need to be traversed. Frequent expression tree construction and destruction operations will generate huge additional time overhead, resulting in its computational efficiency often being lower than that of carefully optimized manual coding methods, making it difficult to meet the needs of efficient numerical simulation.

[0005] Existing technologies cannot provide sufficiently high computational performance while ensuring ease of development, code maintainability, and model scalability, thus failing to meet the actual needs of efficient and accurate simulation of semiconductor devices. An innovative Jacobian matrix construction scheme is urgently needed to resolve the above contradictions.

[0006] Therefore, those skilled in the art are dedicated to developing an efficient method and system for constructing Jacobian matrices in Newton iteration. Summary of the Invention

[0007] Therefore, in view of at least one of the above-mentioned defects or improvement needs of the prior art, the present invention aims to solve the defects of the existing methods for constructing Jacobian matrices in numerical simulation, as well as the contradiction that cannot simultaneously take into account the ease of development, code maintainability, model scalability and high computational performance, and realize a Jacobian matrix construction method and system that combines ease of development, scalability and high computational performance.

[0008] The nonlinear partial differential equations such as the Poisson equation, electron continuity equation, and hole continuity equation solved in numerical simulations take the following forms:

[0009] Discretization is performed using the finite volume method. The unknowns are defined at the vertices of the mesh elements. Elements(i) represents the mesh element to which vertex i belongs, and vertices(e) represents the vertices of mesh element e. V i Let vertex i be the control unit, and the residual of vertex i be:

[0010] As can be seen from the above formula, the residual is formed by the sum of the contributions of edge ij, and the Jacobian matrix is ​​the partial derivative of the residual with respect to the unknowns. Its construction steps include: S1: Construct the Jacobian matrix in CSR (Compressed Sparse Row) format. The Jacobian matrix is a sparse matrix of size (N×M)×(N×M), where N is the number of grid vertices and M is the total number of non - linear partial differential equations to be solved. The Jacobian matrix is divided into M×M blocks, and each block corresponds to the Jacobian matrix of a type of equation with respect to a type of variable. The number of non - zero elements in each row corresponding to a vertex in the Jacobian matrix is the number of adjacent vertices of that vertex plus 1. Allocate the row pointer array, column index array, and non - zero element array in CSR format according to the number of non - zero elements and the total number of vertices, and set the initial value of the non - zero element array to 0. S2: Traverse the M equations, and calculate the N rows of the Jacobian matrix corresponding to each equation each time. Specifically, it includes: S21: Pre - construct the edge residual contribution expression. The edge residual contribution expression is a parametric calculation template, and the process of calculating the contribution of edge ij to the residual of vertex i is recorded through a binary tree based on the variables, material parameters, and calculation region of vertices i and j of edge ij. S22: Traverse all vertices of the grid in multiple threads. Each thread processes one vertex each time to fill a row of the Jacobian matrix. For each vertex i, the thread performs the following operations: S221: Copy the edge residual contribution expression to obtain a thread - local expression instance, and each thread only copies it once. S222: Traverse all edges ij connected to vertex i, calculate the derivative contribution of edge ij to the Jacobian matrix and fill the matrix: If there is symmetry in the partial derivatives of the equation, only perform the calculation for each edge ij when i < j, and skip it when j > i to avoid double - counting the same edge. If the residual calculation is material - related, traverse the subdivision unit to which edge ij belongs and set the calculation region to the subdivision unit; if the residual calculation is not material - related, set the calculation region to all subdivision units around edge ij. During the calculation, first update the parameters of the thread - local expression instance to the actual values of edge ij, then calculate the contribution of edge ij to the residual of vertex i and the partial derivatives of the contribution with respect to all relevant variables of edge ij through reverse automatic differentiation, and finally accumulate the partial derivatives to the corresponding positions of the global Jacobian matrix. If there is symmetry in the partial derivatives of the equation, accumulate the partial derivatives to the row of the Jacobian matrix corresponding to vertex j according to the symmetric or anti - symmetric relationship. Furthermore, the step of constructing the Jacobian matrix in CSR format in step S1 specifically includes: S11: Allocate the row pointer array row_ptr. The dimension of the row pointer array is N×M + 1, and row_ptr[0]=0. S12: For each grid vertex i, where 0≤i≤N-1, obtain the number of adjacent vertices k_i of vertex i, where k_i is the total number of vertices in the grid cell to which vertex i belongs, excluding itself; the number of non-zero elements in the row of the Jacobian matrix corresponding to each vertex i is k_i+1, and for each equation m, set row_ptr[i+m×N+1]=row_ptr[i+m×N]+(k_i+1)×M, where row_ptr[N×M] is the total number of non-zero elements in the Jacobian matrix, where 0≤m≤M-1; S13: Allocate the column index array col_idx, the column index array having a dimension of row_ptr[N×M]; for each grid vertex i, obtain the adjacent vertex numbers and its own number and sort them to obtain the array col_idx_i, the array col_idx_i having a dimension of k_i+1; for each equation m, set col_idx[row_ptr[i+m×N]+j+m×(k_i+1)]=col_idx_i[j]+m×N, where 0≤j≤k_i; S14: Allocate the non-zero element array values, the non-zero element array having a dimension of row_ptr[N×M], and the initial value is set to 0; Furthermore, the variables of vertices i and j of edge ij in step S21 include electric potential, electron concentration, and hole concentration; Furthermore, in the binary tree structure of the edge residual contribution expression in step S21, the leaf nodes are input parameters, which include the variables, material parameters, and calculation regions of the vertices i and j of the edge ij; the parent nodes are arithmetic operations or elementary functions, including addition, subtraction, multiplication, division, exponentiation, dot product, cross product, modulo operation, logarithm, natural exponent, and divergence operation. The parent node performs the corresponding operation on its child nodes, and the root node is the residual contribution of the edge ij to the vertex i. Furthermore, the reverse automatic differentiation calculation process in step S222 is as follows: call the thread local expression instance after updating the parameters, traverse the binary tree from top to bottom, calculate the partial derivatives with respect to all relevant variables according to the chain rule of differentiation, and store the partial derivatives in the array der. The array der has a dimension of 2×M, the first M elements are the partial derivatives of vertex i with respect to the M equation variables of vertex i, and the last M elements are the partial derivatives of vertex i with respect to the M equation variables of vertex j. Furthermore, the process of accumulating the partial derivatives to the global Jacobian matrix in step S222 includes: Based on the vertex number i and the current equation m, locate the row corresponding to i+m×N in the row pointer array row_ptr; The first M partial derivatives of the array der are accumulated and added to the position of row_ptr[i+m×N]+i_local+m_var×(k_i+1) in the non-zero element array values; the last M partial derivatives of the array der are accumulated and added to the position of row_ptr[i+m×N]+j_local+m_var×(k_i+1) in the non-zero element array values. Where i_local is the local index number of vertex i in col_idx_i, j_local is the local index number of vertex j in col_idx_i, and m_var is the number of the equation to which the variable belongs, where 0≤m_var≤M-1; Furthermore, in step S222, if the partial derivatives of the equation exhibit symmetry, the process of accumulating the partial derivatives to the Jacobian matrix row corresponding to vertex j according to a symmetric or antisymmetric relationship includes: First, use binary search to determine the local index of vertex i in col_idx[row_ptr[j+m×N]…row_ptr[j+m×N]+k_j] and the local index of vertex j in col_idx[row_ptr[j+m×N]…row_ptr[j+m×N]+k_j]. Then, the partial derivatives are directly accumulated according to the symmetric relationship or multiplied by -1 according to the antisymmetric relationship and accumulated to the positions of the non-zero element array values[row_ptr[j+m×N]+i_local+m_var×(k_j+1)] and values[row_ptr[j+m×N]+j_local+m_var×(k_j+1)], where k_j is the number of adjacent vertices of vertex j; Furthermore, the copying of the edge residual contribution expression in step S221 is a lightweight copy, which only copies the computation graph structure of the function object and does not reconstruct the expression tree; Furthermore, the nonlinear partial differential equations mentioned in step S1 include the Poisson equation, the electron continuity equation, and the hole continuity equation. Furthermore, in step S222, when traversing the edge ij connected to vertex i, the vertex j corresponding to edge ij is obtained by traversing col_idx[row_ptr[i+m×N]…row_ptr[i+m×N]+k_i], and the vertices j=i are skipped, and only the vertices that constitute the edge are processed.

[0011] This invention also discloses an efficient system for constructing Jacobian matrices in Newton iterations, used to implement the above-mentioned efficient method for constructing Jacobian matrices in Newton iterations. The system includes core functional modules corresponding to the steps of the above method, and the efficient construction of Jacobian matrices is completed through the collaboration of the modules.

[0012] Generally speaking, compared with the prior art, the above technical solutions of the inventive concept can achieve the following beneficial effects: 1. Improve development efficiency and model scalability: By means of the reverse automatic differentiation mechanism to automatically calculate derivatives, there is no need for technicians to manually deduce complex residual-variable analytical derivative formulas, greatly reducing the difficulty of code development and the later maintenance cost; when introducing a new physical model, only corresponding items need to be added to the edge residual contribution expression template, without redesigning the derivative calculation logic, and the system can automatically adapt, significantly enhancing the expansion ability of the simulation software.

[0013] 2. Ensure high computing performance: Through the mechanisms of "pre-building the edge residual contribution expression template" and "thread-local lightweight copy", the core performance bottleneck of frequently dynamically constructing and destroying expression trees in the traversal of a large number of grid cells in traditional reverse automatic differentiation is avoided. The thread-local instance can efficiently complete parameter updates and derivative calculations, and the overall efficiency is close to or even exceeds that of carefully optimized manual coding in complex model scenarios; at the same time, the symmetry of the partial derivatives of the equation (only calculating the edges where i < j and skipping the edges where j > i) is used to further reduce the amount of calculation and improve the construction speed.

[0014] 3. Adapt to multi-core parallel computing: The method has natural multi-thread parallel characteristics. Each thread can independently process a grid vertex, and only needs to perform independent write operations on the values array of the global CSR format Jacobian matrix without locking and synchronization, which can make full use of the resources of modern multi-core processors and further shorten the construction time of the Jacobian matrix.

[0015] 4. Optimize the practicality of simulation: The Jacobian matrix is constructed in CSR sparse format, which not only saves memory occupancy but can also be directly used for subsequent Newton iteration solving, avoiding additional format conversion steps; at the same time, through the design of "adjusting the calculation area according to the correlation between the residual and the material" (when the material is relevant, taking the subdivision unit to which the edge belongs, and when it is irrelevant, taking all the units around the edge), the accuracy of the residual and derivative calculations is guaranteed, providing support for the accurate simulation of physical fields (electrical potential, carrier concentration) and electrical properties (current-voltage characteristics) in numerical simulation, and shortening the overall simulation time. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 is a schematic flowchart of an efficient construction method of the Jacobian matrix in the Newton iteration of the present invention; Figure 2 is a schematic structural diagram of the Jacobian matrix of an efficient construction method of the Jacobian matrix in the Newton iteration of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0018] This invention proposes an efficient method and system for constructing the Jacobian matrix in Newton iteration, such as... Figure 1 As shown, the method includes the following steps: S1: Construct the Jacobian matrix using CSR format. The Jacobian matrix is ​​a sparse matrix of (N×M)×(N×M), where N is the number of grid vertices and M is the total number of nonlinear partial differential equations to be solved (such as the Poisson equation, the electron continuity equation, and the hole continuity equation). The Jacobian matrix is ​​divided into M×M blocks, each block corresponding to a Jacobian matrix of a class of equations with respect to a class of variables. The number of non-zero elements in the row corresponding to each vertex in the Jacobian matrix is ​​the number of adjacent vertices of that vertex plus 1. Based on the number of non-zero elements and the total number of vertices, allocate a row pointer array, a column index array, and a non-zero element array in CSR format. The initial value of the non-zero element array is set to 0. In a more specific embodiment, such as Figure 2 As shown, assume the number of mesh vertices is N and the number of equations is M ( Figure 2 Taking equation 3 as an example, the Jacobian matrix J is a sparse matrix of (N×M)×(N×M). The partial derivative of the i-th residual of the mi-th equation with respect to the j-th variable of the mj-th equation is J[i+mi*N,j+mj*N]. The Jacobian matrix is ​​stored in CSR format, and its construction process includes: S11: Allocate a row pointer array row_ptr with dimensions N*M+1, row_ptr[0]=0; S12: For each grid vertex i (0…N-1), obtain the number of its adjacent vertices k_i. For each equation m (0…M-1), set row_ptr[i+m*N+1]=row_ptr[i+m*N]+(k_i+1)*M, then row_ptr[N*M] is the number of non-zero elements in the Jacobian matrix; S13: Assign a column index array col_idx with dimension row_ptr[N*M]. For each grid vertex i, obtain the number of its neighboring nodes, add it to its own number, and sort it in ascending order to obtain an array col_idx_i with dimension k_i+1 of the number of neighboring nodes. For each equation m, col_idx[row_ptr[i+m*N]+j+m*(k_i+1)]=col_idx_i[j]+m*N, j=0…k_i; S14: Allocate a non-zero element array values ​​with dimension row_ptr[N*M], and initialize the value to 0.

[0019] The Jacobian matrix is ​​constructed using the CSR format, which saves memory and allows it to be directly used for subsequent Newton iterations, avoiding additional transformation steps.

[0020] S2: Traverse the M equations, calculating the N rows of the Jacobian matrix corresponding to each equation in each iteration, specifically including: S21: Pre-constructed equation edge residual contribution expression: Pre-construct an expression for the residual contribution of the current equation for any edge ij connecting vertices i and j. This expression is a parameterized calculation template. Using operator overloading and function overloading mechanisms, the residual calculation formula is automatically converted into a binary tree expression tree. In the binary tree structure, the leaf nodes are input variables, including the equation variables of the edge vertices, material parameters, calculation regions, etc. Among them, there are 2M derivative variables, that is, the equation variables of the two vertices of the edge; the parent node is an arithmetic operation or elementary function, including binary operations such as addition, subtraction, multiplication and division, unary operations such as logarithms and natural exponents, and divergence operations, and performs corresponding operations on the child nodes. The root node is the residual result. The binary tree records how to calculate the contribution of an edge to the residuals of each equation (Poisson, electron, hole) at vertex i, based on variables such as the potential, electron concentration, and hole concentration at vertices i and j, as well as material parameters and the computational domain. This expression allows for the calculation of residuals and their derivatives with respect to variables once the variables are determined. S22: Multi-threaded traversal of all vertices in the mesh. Each thread calculates one vertex at a time and fills a row of the Jacobian matrix. For each vertex i (0…N-1), the thread execution process is as follows: S221: Copy the edge residual contribution expression template. A thread-local storage mechanism is used to copy the pre-built expression, resulting in a thread-local expression instance. This operation avoids multiple threads competing for the same expression resource, and because it copies the computation graph of a lightweight function object rather than rebuilding a complex expression tree, the overhead is minimal. If it has already been copied, it is skipped, ensuring that it is copied only once. The thread-local storage mechanism prevents multiple threads from competing for the same expression resource, ensuring that each thread holds only one independent expression instance.

[0021] S222: Traverse all edges ij connected to vertex i, calculate the derivative contribution of edge ij with respect to the Jacobian matrix, and fill the matrix: Assume the number of edges is \(k_i\). Traverse \(col\_idx[row\_ptr[i + m*N]\ldots row\_ptr[i + m*N]+k_i]\) to obtain the variable vertex number \(j\) for which the derivative is sought, as well as the local indices \(i\_local\) and \(j\_local\) starting from 0 for \(i\) and \(j\) in \(col\_idx\). Skip the vertex with \(j = i\), that is, only process the vertices that form edges. For edge \(ij\), if the partial derivative of this equation has symmetry, then only calculate it once for one edge, only calculate when \(i < j\), and skip it if \(j > i\); if the residual calculation is related to the material, further traverse the subdivision element to which the edge belongs, which is the element in the mesh element \(elements(i)\) to which vertex \(i\) belongs. Set the calculation area as one subdivision element, otherwise the calculation area is all the subdivision elements around this edge.

[0022] Among them, the process of calculating the derivative contribution of edge \(ij\) is as follows: S2221: Update the expression parameters: Update the parameters (such as vertex variable values, local material properties, calculation areas, etc.) in the thread-local expression instance to the actual values of the current edge \(ij\).

[0023] S2222: Perform reverse automatic differentiation to find the derivative: Call the updated expression instance and traverse the expression tree from top to bottom. According to the derivative chain rule, the partial derivatives with respect to all derivative variables can be calculated, and the results are saved in the array \(der[0\ldots 2*M - 1]\). The first \(M\) are the partial derivatives of vertex \(i\) with respect to the variables of vertex \(i\) in \(M\) equations, and the last \(M\) are the partial derivatives of vertex \(i\) with respect to the variables of vertex \(j\) in \(M\) equations. This not only calculates the contribution value of the current edge to the residual of vertex \(i\), but more importantly, through the reverse automatic differentiation mechanism, it simultaneously and efficiently calculates the partial derivatives of this contribution value with respect to all relevant variables (electric potential, electron concentration, hole concentration, etc.) on edge \(ij\).

[0024] S2223: Accumulate the edge derivatives to contribute to the global Jacobian matrix: (1) Locate the row array of the CSR format Jacobian matrix according to the vertex i number, and accumulate the edge derivatives to the corresponding positions in the row array. Specifically, update the i+m*N row of the Jacobian matrix. For the equation m_var(0…M-1) to which the variable belongs, accumulate the partial derivative of i with respect to i: values[row_ptr[i+m*N]+i_local+m_var*(k_i+1)]+=der[m_var], accumulate the partial derivative of i with respect to j: values[row_ptr[i+m*N]+j_local+m_var* (k_i+1)]+=der[m_var+M];(2)If the partial derivatives of this equation have symmetry, locate the row array of the CSR format Jacobian matrix according to the vertex j number, and add (symmetric) or subtract (antisymmetric) the edge derivatives to the row array, that is, update the Jacobian matrix j+m*N rows using the same method as above, that is, the derivatives of j with respect to j and j with respect to i. Assuming that the number of adjacent vertices of j is k_j, then i_local and j_local are the local index numbers in col_idx[row_ptr[j+m*N]…row_ptr[j+m*N]+k_j], which are determined by binary search. Update the der value according to symmetry. If it is symmetric, it remains unchanged. If it is antisymmetric, it is multiplied by -1. The formula for updating the Jacobian matrix is ​​to replace i with j and j with i in (1).

[0025] Correspondingly, the efficient Jacobian matrix construction system in Newton iteration described in this invention ensures efficient and accurate construction of the Jacobian matrix by working in coordination with the core functional modules corresponding to the above method steps.

[0026] The description in this specification is merely illustrative of the invention. Those skilled in the art can make various modifications or additions to the specific embodiments described or use similar methods to replace them, as long as they do not deviate from the content of this specification or exceed the scope defined in the claims, they should all fall within the protection scope of this invention.

Claims

1. An efficient method for constructing the Jacobian matrix in Newton iteration, characterized in that, The method includes the following steps: S1: Construct a Jacobian matrix in CSR format. The Jacobian matrix is a sparse matrix of size (N×M)×(N×M), where N is the number of grid vertices and M is the total number of non - linear partial differential equations to be solved. The Jacobian matrix is divided into M×M blocks, and each block corresponds to the Jacobian matrix of a class of equations with respect to a class of variables. The number of non - zero elements in each row corresponding to a vertex in the Jacobian matrix is the number of adjacent vertices of that vertex plus 1. Allocate the row pointer array, column index array, and non - zero element array in CSR format according to the number of non - zero elements and the total number of vertices. The initial value of the non - zero element array is set to 0. S2: Traverse the M equations, and calculate the N rows of the Jacobian matrix corresponding to each equation each time. Specifically, it includes: S21: Pre - construct an edge residual contribution expression. The edge residual contribution expression is a parameterized calculation template, and the process of calculating the contribution of edge ij to the residual of vertex i is recorded through a binary tree based on the variables, material parameters, and calculation region of vertices i and j of edge ij. S22: Traverse all vertices of the grid in multiple threads. Each thread processes one vertex each time to fill a row of the Jacobian matrix. For each vertex i, the thread performs the following operations: S221: Copy the edge residual contribution expression to obtain a thread - local expression instance, and each thread only copies it once. S222: Traverse all edges ij connected to vertex i, calculate the derivative contribution of edge ij to the Jacobian matrix and fill the matrix: If there is symmetry in the partial derivative of the equation, the calculation is only performed when i < j for each edge ij, and skipped when j > i to avoid double - counting the same edge. If the residual calculation is material - related, traverse the subdivision unit to which edge ij belongs and set the calculation region to the subdivision unit; if the residual calculation is material - independent, set the calculation region to all subdivision units around edge ij. During the calculation, first update the parameters of the thread - local expression instance to the actual values of edge ij, then calculate the contribution of edge ij to the residual of vertex i and the partial derivatives of the contribution with respect to all relevant variables of edge ij through reverse automatic differentiation, and finally accumulate the partial derivatives to the corresponding positions of the global Jacobian matrix. If there is symmetry in the partial derivative of the equation, accumulate the partial derivatives to the row of the Jacobian matrix corresponding to vertex j according to the symmetric or anti - symmetric relationship.

2. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 1, characterized in that, In step S1, the construction of the Jacobian matrix in CSR format specifically includes: S11: Allocate the row pointer array row_ptr. The dimension of the row pointer array is N×M + 1, and row_ptr[0]=0. S12: For each grid vertex i, where 0≤i≤N - 1, obtain the number of adjacent vertices k_i of vertex i. Here, k_i is the total number of vertices in the grid cell to which vertex i belongs except itself. The number of non - zero elements in each row corresponding to vertex i in the Jacobian matrix is k_i + 1. For each equation m, set row_ptr[i+m×N + 1]=row_ptr[i+m×N]+(k_i + 1)×M, and row_ptr[N×M] is the total number of non - zero elements in the Jacobian matrix, where 0≤m≤M - 1. S13: Allocate the column index array col_idx, the column index array having a dimension of row_ptr[N×M]; for each grid vertex i, obtain the adjacent vertex numbers and its own number and sort them to obtain the array col_idx_i, the array col_idx_i having a dimension of k_i+1; for each equation m, set col_idx[row_ptr[i+m×N]+j+m×(k_i+1)]=col_idx_i[j]+m×N, where 0≤j≤k_i; S14: Allocate the non-zero element array values, the non-zero element array has a dimension of row_ptr[N×M], and the initial value is set to 0.

3. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 1, characterized in that, The variables of vertices i and j of edge ij in step S21 include electric potential, electron concentration, and hole concentration.

4. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 1, characterized in that, In the binary tree structure of the edge residual contribution expression in step S21, the leaf nodes are input parameters, which include the variables of vertices i and j of edge ij, material parameters, and calculation regions; the parent nodes are arithmetic operations or elementary functions, including addition, subtraction, multiplication, division, exponentiation, dot product, cross product, modulo operation, logarithm, natural exponent, and divergence operation. The parent node performs the corresponding operation on its child nodes, and the root node is the residual contribution of edge ij to vertex i.

5. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 1, characterized in that, The reverse automatic differentiation calculation process in step S222 is as follows: call the thread-local expression instance after updating the parameters, traverse the binary tree from top to bottom, calculate the partial derivatives with respect to all relevant variables according to the chain rule of differentiation, and store the partial derivatives in the array der. The array der has a dimension of 2×M. The first M elements are the partial derivatives of vertex i with respect to the M equation variables of vertex i, and the last M elements are the partial derivatives of vertex i with respect to the M equation variables of vertex j.

6. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 2, characterized in that, The process of accumulating partial derivatives to the global Jacobian matrix in step S222 includes: Based on the vertex number i and the current equation m, locate the row corresponding to i+m×N in the row pointer array row_ptr; The first M partial derivatives of the array der are accumulated and added to the position of row_ptr[i+m×N]+i_local+m_var×(k_i+1) in the non-zero element array values; the last M partial derivatives of the array der are accumulated and added to the position of row_ptr[i+m×N]+j_local+m_var×(k_i+1) in the non-zero element array values. Where i_local is the local index number of vertex i in col_idx_i, j_local is the local index number of vertex j in col_idx_i, and m_var is the number of the equation to which the variable belongs, where 0≤m_var≤M-1.

7. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 1, characterized in that, In step S222, if the partial derivatives of the equation exhibit symmetry, the process of accumulating the partial derivatives to the Jacobian matrix row corresponding to vertex j according to a symmetric or antisymmetric relationship includes: First, use binary search to determine the local index of vertex i in col_idx[row_ptr[j+m×N]…row_ptr[j+m×N]+k_j] and the local index of vertex j in col_idx[row_ptr[j+m×N]…row_ptr[j+m×N]+k_j]. Then, the partial derivatives are directly accumulated according to the symmetric relationship or multiplied by -1 according to the antisymmetric relationship and accumulated to the positions of the non-zero element arrays values[row_ptr[j+m×N]+i_local+m_var×(k_j+1)] and values[row_ptr[j+m×N]+j_local+m_var×(k_j+1)], where k_j is the number of adjacent vertices of vertex j.

8. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 1, characterized in that, The copying of the edge residual contribution expression in step S221 is a lightweight copy, which only copies the computation graph structure of the function object and does not rebuild the expression tree.

9. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 1, characterized in that, The nonlinear partial differential equations mentioned in step S1 include the Poisson equation, the electron continuity equation, and the hole continuity equation.

10. The efficient method for constructing the Jacobian matrix in Newton iteration according to claim 2, characterized in that, In step S222, when traversing the edge ij connected to vertex i, the vertex j corresponding to edge ij is obtained by traversing col_idx[row_ptr[i+m×N]…row_ptr[i+m×N]+k_i], and the vertices j=i are skipped, and only the vertices that constitute the edge are processed.

Images