Integrated circuit multi-physical field collaborative simulation method based on path tensor manifold analysis

Abstract

This invention relates to the field of circuit design technology and discloses a multi-physics co-simulation method for integrated circuits based on path tensor manifold analysis. This method reconstructs the geometric domain by analyzing the physical layout and manufacturing tolerances, and generates a spatiotemporal input tensor by fusing unsteady operating sequences and roughness features. It predicts the state tensor using a neural network with a feature-shared hidden layer and parallel independent output branches, and establishes a dynamically weighted total loss function through infinite norm statistics. It extracts the pure tangential gradient based on the velocity field spatial Jacobian matrix and tangential projection operator, and establishes microscopic interface slip boundary conditions by combining the slip coefficient calculated from roughness and temperature drift. The boundary conditions are substituted into the loss function for dynamic calibration, and the network parameters are updated through a two-stage hybrid optimization strategy of adaptive gradient descent and quasi-Newton method, outputting a multi-physics distribution tensor solution. This method can eliminate ideal boundary deviations and accurately predict local hotspot risks under dynamic power consumption impacts.

Description

Technical Field

[0001] This invention relates to the field of circuit design technology, specifically to a multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis. Background Technology

[0002] With the continuous miniaturization of integrated circuit feature sizes and the increasing integration density of packaging, the coupling effects of physical fields such as thermodynamics and electromagnetics during chip operation are becoming increasingly significant. Multiphysics co-simulation technology for integrated circuits, as a step in chip design verification, aims to predict the physical state response of the chip under real-world operating conditions by solving the partial differential control equations for multiple fields such as fluid, thermal, and electromagnetic fields using numerical calculation methods.

[0003] Current multiphysics simulations of microchannel heat dissipation architectures for integrated circuits generally rely on traditional macroscopic numerical analysis methods. These methods typically construct the simulation computational domain based on ideal, smooth geometric topological boundaries and use numerical analysis algorithms to solve the energy conservation and other related governing equations of the physical fields. In setting the boundary conditions at the fluid-structure interaction interface, mainstream approaches often employ the classical macroscopic fluid dynamics no-slip assumption, assuming that the velocity of fluid particles adhering to the solid wall is zero. Furthermore, when dealing with multiple boundary inputs, existing methods usually apply static steady-state environmental parameters and constant power consumption loads to obtain steady-state numerical solutions for the fluid velocity, pressure, and temperature fields.

[0004] However, existing multiphysics simulation methods for integrated circuits are still inadequate in handling microscopic fluid-structure interaction mechanisms. They often neglect microscopic geometric defects caused by actual manufacturing processes and frequently employ ideal no-slip assumptions in boundary condition settings. This problem stems from the fact that traditional macroscopic models fail to adequately consider the thermal-fluid coupling effect and dynamic microscopic slip mechanism at the fluid-structure interface under high-power unsteady-state conditions. The idealized approach masks the nonlinear perturbation effect of local process roughness on fluid tangential slip, leading to discrepancies between simulation results and actual physical manifolds. This discrepancy makes it difficult for existing methods to accurately resolve the dynamic thermal-fluid coupling mapping relationship at the microscopic interface, and to accurately predict the risk of local hotspots caused by etching process deviations, thus reducing the accuracy of multiphysics co-simulation in representing real-world complex conditions. Therefore, a multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis is urgently needed to address these issues. Summary of the Invention

[0005] To address the problems in related technologies, this invention provides a multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis, thereby overcoming the aforementioned technical problems in existing related technologies.

[0006] To solve the aforementioned technical problem, the present invention is achieved through the following technical solution:

[0007] This invention provides a method for multiphysics co-simulation of integrated circuits based on path tensor manifold analysis, specifically including:

[0008] Acquire integrated circuit physical layout data, material properties and fluid cooling fluid properties, introduce surface roughness tolerance matrix to construct multiphysics simulation geometric domain, and fuse dynamic working condition time series set to generate spatiotemporal input tensor;

[0009] The spatiotemporal input tensor is input into the physical information neural network, which outputs a multiphysics state tensor. A total loss function is constructed by weighting the data fitting residual tensor, the partial differential equation physical residual tensor, and the boundary condition residual tensor.

[0010] The tangential gradient driving tensor is extracted from the fluid velocity field tensor in the multiphysics state tensor, and the boundary tangential slip velocity tensor is calculated by combining the local slip coefficient tensor and the interface mask continuous weight tensor.

[0011] Substituting the boundary tangential slip velocity tensor into the boundary condition residual tensor, the parameters of the physical information neural network are optimized based on the total loss function optimization strategy, and the multiphysics distribution tensor solution is output.

[0012] As a preferred embodiment of the multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis described in this invention, the method involves acquiring integrated circuit physical layout data, material properties, and fluid cooling medium attributes; introducing a surface roughness tolerance matrix to construct the multiphysics simulation geometric domain; and fusing a dynamic operating condition time series set to generate a spatiotemporal input tensor, including:

[0013] Extract the initial three-dimensional geometric boundary of the microchannel liquid cooling structure and establish the nominal boundary coordinate set and the nominal unit normal vector field;

[0014] The surface roughness tolerance matrix is ​​analyzed to construct the interface micromorphology perturbation field. The interface micromorphology perturbation field is superimposed on the nominal boundary coordinates to generate a set of geometric domain space coordinates for multiphysics simulation.

[0015] A multi-physics dynamic testing platform is adopted. Based on the extracted minimum electromagnetic transient characteristic time constant of integrated circuit and the macroscopic equivalent thermal diffusion time constant, the time window configuration is defined, and the dynamic operating condition time series set within the specified time window is acquired synchronously.

[0016] The root mean square height and local average slope of the local surface at the set of spatial coordinates of the geometric domain of the multiphysics simulation are extracted. The root mean square height of the local surface, the local average slope and the dynamic working condition time series set are mapped to the corresponding spatial grid nodes and time steps and tensor splicing is performed to generate the spatiotemporal input tensor.

[0017] As a preferred embodiment of the multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis described in this invention, the method of constructing the interface micromorphological perturbation field by analyzing the surface roughness tolerance matrix includes:

[0018] Using the Karhunen-Loève expansion method, the integral domain is discretized by the quadrature rule and the Gaussian autocovariance function corresponding to the surface roughness tolerance matrix is ​​solved to obtain the eigenvalue sequence and the orthogonal eigenfunction set.

[0019] Define the eigenvalue energy accumulation contribution rate function, and continuously accumulate the eigenvalues ​​until the energy accumulation contribution rate reaches the preset energy retention threshold to determine the eigenvalue cutoff order;

[0020] Based on the eigenvalue sequence and orthogonal eigenfunction set, an interface micromorphological perturbation field is constructed by generating an independent random variable sequence using the Monte Carlo sampling method.

[0021] As a preferred embodiment of the integrated circuit multiphysics co-simulation method based on path tensor manifold analysis described in this invention, the spatiotemporal input tensor is input into the physical information neural network to output a multiphysics state tensor, and a total loss function is constructed by weighting the data fitting residual tensor, the partial differential equation physical residual tensor, and the boundary condition residual tensor, including:

[0022] High-dimensional features are extracted from spatiotemporal input tensors by utilizing the feature-shared hidden layer configured with continuously differentiable nonlinear activation functions in the physical information neural network.

[0023] The fluid velocity field tensor, fluid pressure field, continuous temperature field and electric potential field are output through independent tensor output branches, forming a multiphysics state tensor.

[0024] Partial derivatives of the multiphysics state tensor are obtained, and the physical residual tensor of the partial differential equation is constructed by combining the reconstructed smooth spatial domain indicator function based on sign distance. The data fitting residual tensor is constructed by combining the real state sequence. The boundary condition residual tensor is constructed based on the reconstructed system geometric domain boundary state, the gradient step scaling operator of the fluid-structure interface, and the two-sided limit decoupling operator.

[0025] The total loss function is generated by weighted summation of the data fitting residual tensor, the physical residual tensor of the partial differential equation, and the residual tensor of the boundary conditions.

[0026] As a preferred embodiment of the integrated circuit multiphysics co-simulation method based on path tensor manifold analysis described in this invention, the generation of the total loss function includes:

[0027] During the training iteration cycle, the automatic differential computation graph is called to extract the maximum absolute value of the gradient of the data fitting residual tensor, the partial differential equation physical residual tensor, and the boundary condition residual tensor relative to the feature-shared hidden layer weight matrix.

[0028] A gradient stability constant threshold and an adaptive weight hard upper limit cutoff value are introduced as stability boundary constraints. The ratio of the maximum absolute value of the gradient is calculated to determine the instantaneous target physical configuration weight value and the instantaneous target boundary condition weight value.

[0029] Based on the set initial weight constant, the instantaneous target physical configuration weight value and instantaneous target boundary condition weight value are updated smoothly over time using the exponential moving average algorithm to generate an adaptive penalty weight tensor.

[0030] The adaptive penalty weight tensor is summed as weighting coefficients to generate the total loss function.

[0031] As a preferred embodiment of the integrated circuit multiphysics co-simulation method based on path tensor manifold analysis described in this invention, the extraction of the tangential gradient driving tensor based on the fluid velocity field tensor in the multiphysics state tensor includes:

[0032] Based on the symbolic distance function of the nominal boundary and the surface projection operator, we construct the implicit surface equation of the geometric domain of the multiphysics simulation with roughness perturbation, calculate the spatial gradient manifold vector of the implicit surface equation, and generate the real local unit normal vector field.

[0033] An automatic differentiation engine based on a physical information neural network is used to obtain the Jacobian matrix of the fluid velocity field tensor with respect to spatial coordinates.

[0034] The Jacobian matrix is ​​multiplied by the real local unit normal field to generate a normal gradient distribution operator tensor.

[0035] Construct a tangential projection operator tensor consisting of the subtraction of the second-order unit tensor and the normal projection tensor;

[0036] The normal component is removed by performing an inner product operation between the tangential projection operator tensor and the normal gradient distribution operator tensor to generate the tangential gradient driving tensor.

[0037] As a preferred embodiment of the integrated circuit multiphysics co-simulation method based on path tensor manifold analysis described in this invention, the calculation of the boundary tangential slip velocity tensor by combining the local slip coefficient tensor and the interface mask continuous weight tensor includes:

[0038] The local slip coefficient tensor is generated by multiplying the smooth interface molecular dynamics fitting benchmark value, roughness correction function and exponential nonlinear temperature correction function obtained based on the ensemble statistics of micro molecular dynamics.

[0039] The effective solid-liquid contact state coefficient is calculated based on the local average slope. The ratio of the effective solid-liquid contact state coefficient to the reference contact coefficient is substituted into the generalized sigmoid smooth step function with denominator lower bound protection. The wetting phase transition steepness and degradation phase transition threshold of the step function are dynamically mapped by the solid-liquid equilibrium contact angle and intrinsic interface adsorption energy to generate the interface mask continuous weight tensor.

[0040] Calculate the Hadamard product of the local slip coefficient tensor, the tangential gradient driving tensor, and the interface mask continuous weight tensor to generate the boundary tangential slip velocity tensor.

[0041] As a preferred embodiment of the multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis described in this invention, substituting the boundary tangential slip velocity tensor into the boundary condition residual tensor includes:

[0042] By introducing a reference velocity based on the macroscopic characteristics of the fluid, the boundary tangential slip velocity tensor is made dimensionless to generate a dimensionless target slip velocity tensor.

[0043] The target slip velocity tensor is assigned to the target slip velocity tensor interface reserved in the total loss function. The L2 norm of the difference between the fluid velocity field tensor and the target slip velocity tensor is used as the loss constraint to construct the instantiated boundary condition residual tensor.

[0044] As a preferred embodiment of the integrated circuit multiphysics co-simulation method based on path tensor manifold analysis described in this invention, the method of optimizing the parameters of the physical information neural network based on the total loss function and outputting the multiphysics distribution tensor solution includes:

[0045] Configure a cosine annealing dynamic learning rate decay function driven by the number of iteration steps;

[0046] An adaptive gradient descent algorithm based on the correction of the first-order momentum estimation tensor and the second-order momentum estimation tensor is invoked to solve the nonlinear gradient transfer of the total loss function on the network parameter tensor and update the network parameter tensor.

[0047] The algorithm monitors the relative smoothing rate of change of the total loss function and the infinite norm of the physical residual components of each partial differential equation within continuous iteration cycles. When the relative smoothing rate of change of the total loss function is lower than the quasi-Newton switching tolerance threshold and the infinite norm of the residuals is lower than the preset residual pre-convergence threshold, the adaptive gradient descent algorithm is suspended. The network parameter tensor is passed to the finite memory quasi-Newton optimization algorithm module, and the diagonal elements of the initial inverse Hessian approximation matrix are constructed using the extracted momentum estimation tensor for curvature warm-start. The second-order approximate Hessian matrix is ​​then solved for residual space fine-tuning.

[0048] As a preferred embodiment of the multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis described in this invention, the residual space fine-tuning process in which the finite-memory quasi-Newton optimization algorithm module solves for the second-order approximate Hessian matrix further includes:

[0049] Extract the term with the largest absolute value of the partial differential residual on all internal configuration nodes, construct a global absolute maximum error metric, and monitor the global absolute convergence state of the physical residual tensor components of the partial differential equation.

[0050] When the global absolute maximum error strictly decreases below the set global physical tolerance threshold, the simulation termination condition is triggered to stop the iterative reconstruction.

[0051] The multiphysics state tensors predicted by the output branches of each independent tensor of the neural network are extracted from the physical information, and then integrated after feature benchmark inverse mapping to generate the multiphysics distribution tensor solution.

[0052] The present invention has the following beneficial effects:

[0053] 1. This invention introduces a surface roughness tolerance matrix and reconstructs a random field using the Karhunen-Loève expansion method. This constructs an interface micro-morphology perturbation field characterizing the spatial correlation of micro-morphology and superimposes it onto the nominal boundary coordinates. This achieves deep integration of micro-geometric defects caused by manufacturing processes with the underlying training node set, accurately quantifying micro-geometric defects. Based on the extracted velocity field normal gradient and tangential projection operator, a fluid-structure interaction micro-interface slip boundary tensor model driven by local roughness and temperature drift is established. The tangential projection operator is applied to the normal gradient distribution operator tensor to extract the pure tangential gradient tensor, stripping away the spurious normal slip driving force caused by geometric distortion and eliminating the normal penetration velocity. This overcomes the problem of traditional macroscopic fluid dynamics' absolute no-slip assumption neglecting the fluid-structure interface thermal-fluid coupling effect and dynamic micro-slip mechanism, eliminating the deviation between simulation results and actual physical manifolds.

[0054] 2. This invention combines a test platform with high-frequency physical parameter acquisition equipment, utilizes a multi-scale time parameter set to define the time window for dynamic data acquisition, and synchronously acquires the real-time power consumption fluctuation sequence, substrate ambient temperature drift, and power supply voltage noise within the specified time window, achieving deep fusion of unsteady-state conditions at the input end. The physical information neural network captures multi-physics state tensors through feature-shared hidden layers and parallel independent tensor output branches for fluid, thermal, and electromagnetic fields. It also configures an online dynamic adjustment mechanism for adaptive penalty weight tensors based on a unified maximum norm benchmark, eliminating drastic fluctuations across magnitudes and preventing loss of numerical dominance in the master physical partial differential constraints. Combined with a two-stage hybrid optimization switching mechanism based on the gradient evolution of the loss function, it ensures numerical completeness and convergence stability during long-term iterative convergence. This allows it to maintain high-precision multi-physics prediction capabilities under actual non-ideal conditions such as dynamic power consumption shocks and ambient temperature drift, accurately predicting local hotspot risks caused by etching process deviations, and improving the accuracy of multi-physics co-simulation in representing real complex conditions.

[0055] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0056] To more clearly illustrate the technical solutions of the embodiments of the invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the invention. For those skilled in the art, the drawings can be obtained from these drawings without creative effort.

[0057] Figure 1 The flowchart of the multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis is provided by the present invention.

[0058] Figure 2 This is a schematic diagram of the process S3 provided by the present invention.

[0059] Figure 3 This is a schematic diagram of the S4 process provided by the present invention. Detailed Implementation

[0060] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0061] Example 1

[0062] Existing multiphysics simulation methods for integrated circuits often ignore microscopic geometric defects caused by actual manufacturing processes, and often adopt the ideal no-slip assumption in setting boundary conditions. They do not fully consider the thermal-fluid coupling effect and dynamic micro-slip mechanism of the fluid-solid interface under high-power unsteady-state conditions, resulting in deviations between simulation results and actual physical manifolds.

[0063] To solve the above technical problems, such as Figure 1 As shown, Embodiment 1 of the present invention provides a multi-physics co-simulation method for integrated circuits based on path tensor manifold analysis. Specifically, Embodiment 1 takes the fabrication and operation scenario of a high-power artificial intelligence chip with a microchannel liquid cooling architecture as an example: the chip's heat dissipation microchannels are manufactured using deep silicon etching technology, resulting in inherent process roughness on the solid walls, and the cooling medium is deionized water.

[0064] In the specific implementation of the above embodiment 1, firstly, the physical layout data of the integrated circuit, material properties and cooling fluid properties are acquired, and a multi-physics geometric domain is constructed by introducing a surface roughness tolerance matrix, and a spatiotemporal input tensor is generated by fusing unsteady dynamic operating parameters; a physical information neural network containing independent tensor output branches is constructed, and the multi-physics control equations are transformed into tensor differential operators and incorporated into the loss function; based on the extracted velocity field normal gradient and tangential projection operator, a microscopic interface slip boundary tensor model driven by local roughness and temperature drift is established; this boundary tensor model is substituted into the boundary condition residual term of the total loss function for dynamic calibration, driving the neural network output to approximate the real multi-field coupled manifold solution.

[0065] Furthermore, in order to better introduce the technical solution of Embodiment 1 of the present invention, the collaborative simulation method is described in detail, specifically including the following:

[0066] S1. Acquire multi-source physical and process data of integrated circuits, construct a geometric domain containing local morphological features, and integrate dynamic operating conditions to construct a spatiotemporal input tensor. This tensor is used to quantify microscopic geometric defects caused by manufacturing processes and unsteady operating conditions of the equipment, and is intrinsically embedded into the simulation origin. Specifically, this includes the following sub-steps:

[0067] S11. Analyze the physical layout data of the integrated circuit, extract the initial three-dimensional geometric boundary of the microchannel liquid cooling structure, and establish the nominal boundary coordinate set. The nominal unit normal vector field at the solid-liquid interface Material properties and fluid cooling medium properties are obtained to construct a background property matrix. Based on the manufacturing tolerances of the deep silicon etching process for semiconductors, the root mean square height of the surface is obtained. With autocorrelation length Based on manufacturing statistical characteristics, a surface roughness tolerance matrix is ​​constructed. The Karhunen-Loève expansion method is introduced as an analytical method for reconstructing the surface roughness tolerance matrix in a random field algorithm. Establish a Gaussian autocovariance function to characterize the spatial correlation of micromorphology:

[0068] ;

[0069] In the formula, Represents any two points in a spatial coordinate system and The autocovariance tensor scalar between them, with dimensions of ; Represents the root mean square height of the surface, with dimensions of ; and These represent the position vectors on the walls of the two-dimensional microchannel, both having dimensions of [missing information]. ; This represents the autocorrelation length, with dimensions of . ; This represents the 2-norm operator.

[0070] To solve the second kind of Fredholm integral equation corresponding to this Gaussian autocovariance function, the Nyström method combined with the Gauss-Legend quadrature rule is used to discretize the integration domain. Using integration nodes, we construct a standard algebraic eigenvalue solving problem and extract the eigenvalue sequence using singular value decomposition. orthogonal characteristic function set The cumulative contribution rate function of eigenvalue energy is defined to quantify the proportion of eigenvalues; its calculation formula is as follows:

[0071] ;

[0072] In the formula, Indicates truncation up to the th The cumulative energy contribution rate at time step 1, dimensionless; and These represent the order of the descending order of the first... Rank and number eigenvalues ​​of order 1, all of which are units of 1. ; This indicates the truncation order of the characteristic term, and is dimensionless. This represents the total number of discrete nodes and is dimensionless. Under extreme conditions of extreme mesh refinement, if the total number of discrete nodes... As the value approaches infinity, the denominator of the formula converges to the analytic trace of a continuous Gaussian autocovariance function, ensuring the validity of the continuous physical limit of discrete energy accounting. The eigenvalues ​​are continuously accumulated until... The preset energy retention threshold (e.g., 0.95) is reached to determine the final truncation order of the feature terms. If the local autocorrelation length The extremely small value causes the eigenvalue decay to stagnate, reaching the set absolute cutoff upper limit. If the energy retention threshold is still not met, the remaining unresolved high-frequency energy is aggregated into a white noise residual field and superimposed onto the reconstructed surface to ensure strict conservation of local geometric statistical energy. A spatially continuous interface micromorphological perturbation field is constructed based on the acquired feature components.

[0073] ;

[0074] In the formula, Representing spatial coordinates The height of the microscopic morphology perturbation at the interface, with dimensions of ; This represents a sequence of independent random variables that follow a standard normal distribution and is dimensionless. Indicates the corresponding first orthogonal characteristic function of order 1, dimensionless For sequences of independent random variables The Monte Carlo sampling method is invoked, and a global random seed bound to the physical system's timestamp is configured to ensure the reproducibility of microscopic geometric samples. The number of Monte Carlo samplings is set to [number missing]. In each sampling cycle, a set of dimensions is generated. Substituting a standard normally distributed random vector into the expansion formula generates the single-instance perturbation height. This interface micro-topography perturbation field is then superimposed onto the nominal boundary coordinates to generate the set of multiphysics simulation geometric domain coordinates after a single sampling reconstruction. In the formula, This represents the set of coordinates in the reconstructed multiphysics simulation geometric domain, with dimensions of ; Denotes the set of nominal boundary coordinates, with dimensions of ; This represents the nominal unit normal vector field of the solid-liquid interface, which is dimensionless. (To be continued...) After the sub-microscopic physical reconstruction and multi-physics collaborative solution are completed, an ensemble averaging operation is performed on the output tensor set to extract stable predictive statistical moments. A non-uniform tensor product mesh with boundary layer encryption is used to spatially discretize the geometric domain space after single-instantiation reconstruction, generating a low-level training node set containing spatial topological constraints.

[0075] S12. Construct a test platform including high-frequency physical parameter acquisition equipment to obtain the dynamic operation sequence of the integrated circuit under non-ideal conditions, and extract morphological features from the surface roughness tolerance matrix based on local window integration. The specific implementation steps are as follows:

[0076] S121. Construct a multi-physics dynamic test platform including a high-frequency power analyzer, a substrate-embedded RTD array, and a high-bandwidth oscilloscope. Before the driving integrated circuit enters high-power unsteady-state operation, extract the intrinsic time scale of the physical system to determine the preset sampling frequency. With a specified time window The value boundaries are established. A time constant incorporating electromagnetic transient characteristics is established. With thermal diffusion time constant A multi-scale time parameter set.

[0077] Extract the parasitic parameter matrix of the integrated circuit power supply network and high-frequency switch, and calculate the local time constant of each node. Configure the transient energy contribution rate threshold. The set of effective response mode nodes, where the proportion of transient energy to total electromagnetic wave energy is greater than a certain threshold, is selected. The minimum local time constant corresponding to this set of nodes is taken as the minimum electromagnetic transient characteristic time constant. The calculation formula is:

[0078] ;

[0079] In the formula, The minimum electromagnetic transient characteristic time constant of the integrated circuit power supply network and high-frequency switch is represented by the dimensionless time constant. ; Indicates a local node or parasitic element branch index in the power supply network; dimensionless. Indicates the first The transient energy contribution rate of each node is dimensionless. This represents the preset transient energy contribution rate threshold, which is dimensionless. and They represent the first The equivalent series resistance and capacitance of each decoupling capacitor branch, with dimensions respectively... and ; and They represent the first The inductance and equivalent resistance of each parasitic inductance branch, with dimensions respectively and When the system is in an ideal parasitic parameter-free operating condition, i.e. or hour, tending to It characterizes the transient response as infinitely fast, ensuring the physical and logical integrity under extreme operating conditions.

[0080] The non-uniform material layer distribution in the geometric domain of multiphysics simulation is analyzed, and the macroscopic equivalent thermal diffusion time constant is extracted based on the series thermal resistance equivalent and volume weighting principle. The calculation formula is:

[0081] ;

[0082] In the formula, The macroscopic equivalent thermal diffusion time constant of the chip substrate and liquid cooling structure is represented by the dimensionless value. ; This represents the equivalent thermal resistance of a non-uniform material along the dominant heat transfer direction, with dimensions of . ; This represents the equivalent total heat capacity of the system, with dimensions of . ; This represents the equivalent discrete layer number of a non-homogeneous material and is dimensionless. Indicates the first The thickness of the layer, with dimensions of ; Indicates the first The thermal conductivity of the layer, with dimensions of ; Indicates the first The thermal conductivity cross-sectional area of ​​the layer, with dimensions of ; Indicates the first The density of the layer, with dimensions of ; Indicates the first The specific heat capacity of the layer, with dimensions of ; Indicates the first The volume of the layer, with dimensions of When the thermal conductivity of the material layer At that time, equivalent thermal resistance The macroscopic equivalent thermal diffusion time constant approaches This conforms to the physical extreme convergence logic of a superconducting heat body instantly reaching thermal equilibrium. The time window configuration for dynamic data acquisition is defined by the following mathematical formula:

[0083] , ;

[0084] In the formula, This indicates the preset sampling frequency for multi-physics synchronous operation, with dimensions of . ; This represents the sampling anti-aliasing safety factor, which is dimensionless and its parameter value is determined according to the Nyquist sampling theorem. To fully capture high-frequency voltage noise; Represents the maximum deadline for a specified time window, with units of . ; Represents the thermodynamic relaxation multiplier, dimensionless, taking parameter values. To ensure that the time window can cover more than The thermal transient evolution and relaxation process. If the integrated circuit under test evolves into a purely steady-state DC power supply with no dynamic power consumption fluctuations, the minimum electromagnetic transient characteristic time constant is achieved. sampling frequency Converging to a minimum steady-state monitoring frequency; if the system is in a physical state of infinite heat capacity or extremely high thermal resistance, causing thermodynamic evolution to stagnate, i.e., the thermal diffusion time constant... The maximum deadline for the specified time window This extends to positive infinity to cover the slow thermal drift process, ensuring dimensional consistency and mathematical model completeness across physical scales. A multiphysics dynamic testing platform is used at the aforementioned predetermined sampling frequency. Synchronously acquire specified time window Real-time power fluctuation sequence within substrate ambient temperature drift and power supply voltage noise This constitutes a dynamic operating condition time series set.

[0085] S122. Establish a spatial slip calculation window at the fluid-structure interface. The interface micromorphological perturbation field generated in step S11 is analyzed. The microchannel solid wall in coordinate system is extracted by spatial integration. The local roughness morphology parameters at the location. Calculate the root mean square height of the surface. With average slope :

[0086] ;

[0087] ;

[0088] In the formula, Representing spatial coordinates The local root-mean-square height of the surface at a given location, with dimensions of . ; Representing spatial coordinates The local average slope at a given location is dimensionless. Represented by coordinates A local two-dimensional geometric integration window centered on the locating element; The physical area of ​​the local two-dimensional geometric integration window, with dimensions of ; The dimension of the microscopic topographic perturbation at the integral element is . ; and Let each represent a spatial integral infinitesimal element in an orthogonal two-dimensional coordinate system, both having dimensions. Under extreme conditions, if the semiconductor processing technology tends towards ideal absolute smoothness, the microscopic morphology perturbation in all regions within the integration window will be highly variable. Then the integral results of the above two equations will automatically converge to It conforms to the physical extreme value convergence logic of an ideal smooth surface without roughness, ensuring dimensional consistency and the completeness of the mathematical model.

[0089] S123. Discretize the spatial coordinates of the reconstructed multiphysics simulation geometric domain into a mesh, map the dynamic working condition time series set and the extracted local roughness morphological parameters to the corresponding spatial mesh nodes and time steps, and perform tensor concatenation along the spatial coordinate dimension and the time dimension to generate a high-dimensional spatiotemporal input tensor. This is used as the bottom training node set of the physical information neural network to achieve deep integration of the micro-geometric defects of the manufacturing process and the actual unsteady working conditions of the equipment at the input end.

[0090] S2, Receive the spatiotemporal input tensor output in step S1 This is mapped to a multiphysics output manifold, a physical information neural network for multiphysics coupling solution is constructed, and a dynamically adaptive weighted total loss function is defined. The specific sub-steps are as follows:

[0091] S21. Constructing a physical information neural network architecture with multi-branch tensor output: This involves converting the discretized spatiotemporal input tensor... The input is fed into a multi-layer fully connected feedforward neural network. The topology of this network consists of a shared hidden layer in the preceding stage and parallel, independent output branches in subsequent stages. The shared hidden layer captures the low-level nonlinear coupling mapping relationships between multiple physics fields through feedforward mapping, and its evolution equation is expressed as:

[0092] ;

[0093] In the formula, Indicates the first The layer features share the output feature tensor of the hidden layer, which is dimensionless; The input is a dimensionless spatiotemporal input tensor with input dimension . ; This represents the total number of hidden layers sharing features, and is set to a fixed integer of 6. Indicates the first Layer features share the weight matrix of the hidden layer, which is dimensionless. Its dimension is ,when Its dimension is ; Indicates the first The layer features share the bias vector of the hidden layer, which is dimensionless and has a dimension of 256; This represents a continuously differentiable nonlinear activation function whose first and higher-order derivatives are continuous over the entire domain. Its computational expression is: In specific conditions such as extreme truncation or the absence of higher-order feature evolution, if the input tensor degenerates to a constant, the outputs of each hidden layer smoothly transition through the zero-point self-balancing property of the Swish function, ensuring the numerical stability of the network's forward propagation.

[0094] The last layer output of the feature-sharing hidden layer It serves as the common input for subsequent parallel independent output branches. The parallel independent output branches include a fluid field output branch, a thermal field output branch, and an electromagnetic field output branch. These three branches do not share parameters to eliminate direct interference from gradients across physical fields. Each independent output branch is internally composed of an independent hidden layer and a linear output layer connected in series. Its general mathematical expression is:

[0095] ;

[0096] ;

[0097] In the formula, These are placeholders representing the fluid field, thermal field, and electromagnetic field branches, respectively. This represents the number of independent hidden layers within each independent output branch, and is set to a fixed integer of 2. Indicates the first branch within the corresponding branch The output feature tensor of the layer is dimensionless. ; Indicates the first branch The independent weight matrix of the layer is dimensionless. Its dimension is ,when Its dimension is ; Indicates the first branch The independent bias vector of the layer is dimensionless and has a dimension of 128. and These represent the weight matrix and bias vector of the linear output layer at the end of each branch, respectively, and are dimensionless. This represents the physical field prediction tensor output by each branch after decoupling.

[0098] The linear output layer parameter dimensions and mapping relationships for each independent output branch are defined as follows: For the fluid field output branch, its output weight matrix... Dimensions Output bias vector The dimension is 4, and the output is a dimensionless state tensor. Includes fluid velocity vector field The three spatial components and the fluid pressure field For the thermal field output branch, its output weight matrix... Dimensions Output bias vector The dimensionless state tensor output is 1. For the global temperature field of solid and liquid For the electromagnetic field output branch, its output weight matrix... Dimensions Output bias vector The dimensionless state tensor output is 1. electric potential field The outputs of each branch are ultimately combined to form a complete dimensionless multiphysics state tensor. .

[0099] This embodiment establishes a skeleton network capacity capable of fully representing flow, heat transfer, and electromagnetic coupling information by finely limiting the depth and width of the feature-sharing hidden layer. This avoids underfitting due to an insufficiently small network size or declining generalization performance due to an excessively large network size. The activation function explicitly adopts the Swish function, which is unbounded above, bounded below, smooth, and non-monotonic. It is also highly differentiable across the entire domain, satisfying the stringent continuity requirements for calculating high-order partial derivatives of the Navier-Stokes equations and energy conservation equations during automatic differentiation in physical information neural networks. This eliminates the numerical breaks that occur during differentiation of conventional piecewise linear functions. Parallel and independent tensor output branches are designed in subsequent stages, with independent hidden layer structures configured within each branch, achieving complete decoupling of the mapping from shared features to specific physical fields. This multi-head independent branch architecture ensures that the specific gradient flow of the residuals of the governing equations of different physical fields mainly acts on their respective dedicated branch parameters during backpropagation. This avoids ill-conditioned gradient cancellation and vicious competition in the multi-task joint learning process of multi-field partial differential equations, and improves the convergence speed and solution accuracy of the network when handling cross-scale physical field co-simulation.

[0100] S22. Extract the multiphysics partial differential control equations within the microchannel liquid cooling and chip heating system, transform them into dimensionless tensor differential operators, and establish a physical residual tensor set. Specifically, use the automatic differentiation technique of neural networks to calculate the state tensor. With respect to the partial derivatives of the spatiotemporal coordinate components, construct the following physical control equation residuals:

[0101] Electromagnetic field Joule heating effect residual tensor and its heat generation rate source item :

[0102] , ;

[0103] In the formula, Let the spatial gradient operator be denoted by . ; The conductivity of a semiconductor or metal interconnect is expressed in units of 1. ; Represents the electric potential field, with dimensions of ; This represents the residual of the charge conservation equation, with dimensions of ; The volumetric heat generation rate tensor of space is denoted by , and its dimensions are . .

[0104] Fluid dynamic mass conservation residual tensor With momentum conservation residual tensor The mathematical expression is as follows:

[0105] , ;

[0106] In the formula, The velocity vector field tensor of the cooling working fluid has dimensions of ; The fluid pressure field tensor has dimensions of . ; Represents fluid density, with dimensions of ; Represents the dynamic viscosity of a fluid, with dimensions of ; Dimensions are ; Dimensions are .

[0107] Global energy conservation residual tensor The expression is:

[0108] ;

[0109] In the formula, Represents the continuous temperature field tensor, with dimensions of ; This represents the local density of a material, with dimensions of . ; This represents constant pressure heat capacity, with dimensions of _____. ; Represents the local thermal conductivity, with dimensions of ; Dimensions are In the solid-state computational domain, physical properties are assigned solid-state parameters and the velocity vector field tensor is forcibly set. At this point, the momentum conservation residual tensor loses its physical meaning and no longer participates in the solid domain calculation, while the convective heat transfer term... Automatically degenerate into The energy conservation residual formula smoothly transitions to a solid-state pure heat conduction partial differential equation; under the extreme specific condition of power disconnection, i.e., the electric potential field... And when there is no external heat source input, the volumetric heat generation rate tensor Automatic convergence to The system evolves towards an isothermal equilibrium state, ensuring the physical and logical integrity of the tensor operator under various preset and extreme application conditions.

[0110] In this embodiment, the local thermal conductivity Based on the material property parameters obtained in step S11, a basic thermal conductivity constitutive model is constructed. The spatial mapping of the fluid-structure interaction computational domain is performed based on the geometric spatial coordinates reconstructed in step S11. The internal node-level real-time nonlinear calculation is then performed using the continuous temperature field tensor predicted by the independent output branches of the neural network in step S21 as the independent variable.

[0111] ;

[0112] In the formula, Representing spatial coordinates The local thermal conductivity at a given temperature, with dimensions of ; and These represent the fluid cooling medium and solid wall material obtained in step S11 at the standard reference temperature, respectively. The fundamental thermal conductivity, with dimensions of ; and These represent the intrinsic thermal conductivity temperature coefficients of the above materials, with dimensions of... ; This represents the continuous temperature field tensor predicted by the thermal field output branch of the physical information neural network in step S21, with dimensions of... ; This represents the spatial domain Boolean indicator function determined by the spatial boundary of the multiphysics simulation geometry reconstructed in step S11. It is dimensionless and occurs when the spatial coordinates... When located in the internal region of a solid, the element takes the value of When located within the fluid region, the value is [value]. .

[0113] To achieve continuous differentiability of the fluid-structure interaction computational domain and ensure smooth gradient propagation in the automatic differential computation graph within the physical information neural network framework, the spatial domain indicator function is reconstructed using the signed distance function. Define the shortest orthogonal distance from any spatial coordinate in the multiphysics simulation geometry to the fluid-structure interaction interface as the sign distance, and construct a smooth step function based on the sign distance to replace the discrete Boolean assignment:

[0114] ;

[0115] In the formula, Representing spatial coordinates The smooth spatial domain indicator function at the location is dimensionless; Representing spatial coordinates The signed distance function at point , with dimensions . The convention is to take positive values ​​inside the solid region, negative values ​​inside the fluid region, and zero values ​​at the fluid-structure interaction interface; The domain transition steepness coefficient, with dimensions of This is used to control the physical width of the fluid-solid boundary transition zone. When the spatial coordinates are deep within the solid region, The exponential term tends to decay , Smooth convergence to extreme values When the spatial coordinates are deep within the fluid region, The exponent term tends to infinity. Smooth convergence to extreme values When at a fluid-structure interaction interface, , Values This achieves a smooth transition between the two phases. A continuous, globally differentiable sign-distance smoothing step function is introduced to replace the discrete step function, eliminating the Dirac delta function singularity generated when the neural network performs automatic differentiation of higher-order partial derivatives across phase boundaries. The nonlinear spatial gradient of the physical residual tensor of the partial differential equations with respect to the network input coordinate features is continuously and completely backpropagated to the feature-sharing hidden layer nodes via a smoothing spatial domain indicator function, maintaining the numerical completeness and convergence stability of the computational graph gradient flow in the multi-domain coupled solution system.

[0116] S23. Incorporate tensor differential operators as hard physical constraints into the loss function of the neural network, construct a dimensionless total loss function containing a general boundary tensor interface, and extract the statistical characteristics of the backpropagation gradient to achieve online dynamic weight adjustment. Specifically, this includes the following sub-steps:

[0117] S231. Define the dimensionless total loss function. Fitting residual tensors to data Physical residual tensor of partial differential equations With boundary condition residual tensor The dynamic weighted sum:

[0118] ;

[0119] Among them, the data fitting residual tensor Based on the real-state sequence collected by the physical testing platform in step S1, the specific mathematical formula is as follows:

[0120] ;

[0121] In the formula, This represents the total number of discrete measurement data nodes actually acquired, and is dimensionless. Representing the neural network in spatiotemporal coordinates The dimensionless predicted multi-field state tensor is output and scaled by the feature reference value. It represents the dimensionless actual physical state tensor measured at the same spacetime coordinates.

[0122] Among them, the physical residual tensor of partial differential equations The differential operators formed by the dimensionless processing of the control equations of each physical field extracted in step S22 are aggregated, and the specific mathematical formula is as follows:

[0123] ;

[0124] In the formula, Represents the total number of unlabeled configuration nodes randomly sampled within the physical computing domain; it is dimensionless. , , , Representing nodes respectively The dimensionless residual operators for charge, mass, momentum, and energy conservation are given. Under extreme conditions where there is no external heat source input or the fluid velocity field evolves to a steady state, the relevant time-domain partial derivatives and the heat source term automatically cancel each other out. This corresponds to the residual tensor contribution converging to a minimum, maintaining the numerical completeness of the multi-field operator aggregation process.

[0125] Among them, the boundary condition residual tensor To characterize the deviation of the output boundary state from the external target constraint, a general tensor interface is configured in the loss function. The specific mathematical formula is as follows:

[0126] ;

[0127] In the formula, Represents the total number of configuration nodes sampled on the geometric domain boundary of the reconstructed system, dimensionless; This represents the dimensionless fluid prediction velocity field tensor. This represents the dimensionless target slip velocity tensor interface at the fluid-structure interaction interface. This tensor interface serves as a placeholder, and its specific input data is provided by the dynamic output of subsequent generation steps. The fluid boundary operator tensor is dimensionless. This represents the set of dimensionless boundary operator tensors, including the fluid inlet / outlet boundaries, wall heat flux boundaries, and constant voltage boundaries. Under extreme test conditions where there are no external physical constraints (i.e., the process is in a state of completely free evolution), the target boundary constraint tensor interfaces all degenerate into zero tensors, and the boundary residuals are forced to zero, ensuring the mathematical completeness of the loss function structure.

[0128] To address the difficulty in fitting neural networks to gradient abrupt changes at the interface caused by the discontinuous abrupt changes in the thermal conductivity of the fluid working medium and the solid wall, a two-sided limit partitioning decoupling operator and a gradient step scaling operator are introduced into the interface residual construction. This is applicable to any configured node at the interface. Introducing the local normal vector along the true local unit vector Infinitely small physical distance Constructing virtual sampling points on the fluid side Virtual sampling points on the solid side The thermal conductivity and spatial temperature gradient at the two independent limiting positions are obtained using an automatic differentiation mechanism. The residual penalty weight in the high-property-property-contrast region is adaptively amplified by the interface gradient step scaling operator, forcing the network parameters to accurately capture the sharp inflection of the thermal conductivity gradient when crossing heterogeneous phase boundaries. The specific construction formula for the heat flux density continuity residual operator is as follows:

[0129] ;

[0130] ;

[0131] In the formula, and Let these represent the real-time local thermal conductivity at the virtual sampling points on the fluid and solid sides, respectively, with dimensions of... ; and These represent the spatial gradient tensor values ​​of the continuous temperature field, calculated by automatic differentiation on the fluid and solid sides, respectively, with dimensions of... ; This represents the true local unit normal vector of the fluid-structure interface, which is dimensionless. Indicates the first The gradient step scaling operator values ​​at each interface node are dimensionless. This represents the sensitivity coefficient for discontinuity of the interface gradient, and is dimensionless. Under the specific working condition of an ideal homogeneous system with identical fluid-structure interaction properties on both sides, the property ratio term automatically decays to 0, the interface gradient step scaling operator smoothly degenerates to a conventional constant of 1, and the gradient residuals on both sides maintain a continuous operator transition, ensuring the dimensional consistency and physical logic completeness of the formula under various complex boundary conditions and multiple field abrupt changes.

[0132] This embodiment supplements the boundary residual system with an explicit constraint on the continuity of interfacial heat flux density, breaking the smoothing and passivation limitations of traditional single continuous loss functions when dealing with thermally conductive and variable property interfaces. By introducing a virtual point-pair bilateral limit sampling operator, the directional derivatives of the fluid and solid domains at the phase boundary are explicitly decoupled in the automatic differential calculation graph, allowing the network to retain reasonable physical sharp cusp characteristics at the zero-measure interface. Combined with a gradient step adaptive scaling operator based on the logarithmic mapping of property contrast, the penalty extremum can be automatically increased in the high-frequency thermal flux impact region across media, preventing the risk of local underfitting of interface boundary conditions due to global smoothing optimization, and ensuring the analytical accuracy of the multi-physics manifold across phase boundaries.

[0133] S232. Starting from the construction of the underlying numerical computation graph, the definition of the initial state, and the stability control at extreme points, an automatic differential gradient transfer and adaptive penalty weight tensor online dynamic adjustment mechanism is constructed. In the forward propagation computation stage, the discretized spatiotemporal input tensor undergoes feature extraction via a feature-shared hidden layer parameter matrix, and is then fed into each independent tensor output branch to calculate the multiphysics state tensor. Based on this, the automatic differential engine constructs an operation dependency tree containing all forward operator nodes and a unified automatic differential computation graph. In the backpropagation computation stage, the automatic differential engine executes a chain rule optimization along this unified automatic differential computation graph, fitting the data to the residual tensor. Physical residual tensor of partial differential equations With boundary condition residual tensor The scalar loss flow is propagated backwards, and each independent residual term is isolated and calculated relative to the feature-shared hidden layer parameter matrix. The partial derivative matrix is ​​used to generate the corresponding data gradient tensor, physical residual gradient tensor, and boundary gradient tensor. A fixed penalty weight is set for the data residual. It is always equal to 1. In each training iteration cycle... In this process, the maximum absolute value component of each residual gradient tensor across all shared weight dimensions is extracted using an automatic differential computation graph. Based on the adaptive matching criterion of introducing a stability guarantee operator and a maximum boundary threshold, the instantaneous target physical configuration weight values ​​within this iteration cycle are determined. Weight values ​​of instantaneous target boundary conditions The calculation formula is:

[0134] ;

[0135] ;

[0136] In the formula, and They represent the first The instantaneous target physical configuration weight value and instantaneous target boundary condition weight value are calculated in each training iteration cycle and are dimensionless. This represents the gradient partial derivative operator with respect to the feature-shared hidden layer parameter matrix; This indicates the extraction of the maximum absolute value component of the gradient tensor matrix across all parameter dimensions; it is dimensionless. This represents the set gradient stability constant threshold, which is dimensionless, for example, a constant value. This is used to prevent numerical singular divergence caused by denominator collapse when the system encounters local extrema or saddle points, causing the physical residual gradient to approach zero extremely. This represents the set hard upper limit cutoff value for adaptive weights, which is dimensionless, such as a constant value. This is used to forcibly lock the upper bound of weight growth under extreme degenerate physical boundary conditions, preventing a single physical constraint from experiencing a sharp increase in parameters and depriving other constraints of their learning dominance, thus ensuring the convergence stability of the entire network's parameter updates. When the absolute maximum value of the residual gradient of the data fitting approaches 0 due to perfect fitting, the instantaneous target weights converge linearly to 0, conforming to the dynamic multi-task balance convergence logic and possessing complete physical extreme value compatibility and dimensional consistency. The exponential moving average algorithm is used to perform time-step smooth recursive updates on the adaptive penalty weight tensor, and the derivation formula is as follows:

[0137] ;

[0138] ;

[0139] In the formula, and , respectively, represent the adaptive penalty weights that finally apply to the total loss function after smoothing through nonlinear time steps, and are dimensionless; Represents the exponentially decreasing rate of decay, dimensionless, with a value range configured as follows: For example, 0.95. To eliminate the boundary bias abrupt change problem caused by the lack of historical recursion in the initial training stage of the exponential moving average algorithm, the initial iteration period is explicitly defined. The weight boundary conditions at the point are given by: and This establishes an unbiased initial benchmark for the collaborative evolution of network parameters across multiple domains. Through the computational mechanism combining gradient graph tracing, dual stable boundary constraints, and unbiased initialization, the system achieves adaptive and stable adjustment of neural network parameter weights across physical scales, with rigorous mathematical convergence guarantees.

[0140] S3. Based on the physical field features extracted from the output manifold, a fluid-structure interaction microscopic interface slip boundary tensor model is established, which is linked by geometric roughness and dynamic working conditions. This model is used to break the absolute no-slip assumption of traditional macroscopic fluid mechanics, and to generate accurate boundary constraints using the morphological parameters extracted in step S1 and the velocity field output in step S2. Figure 2 As shown, the specific sub-steps are as follows:

[0141] S31. Extract the fluid velocity field tensor from the output of the fluid branch of the physical information neural network, calculate the true normal gradient by combining it with the local geometric features of the micro-rough interface, and introduce a tangential projection operator to eliminate the normal component. The specific implementation steps are as follows:

[0142] S311, The nominal boundary of the microchannel liquid cooling structure adopts an analytical signed distance function. Implicit representation is performed, and its spatial gradient corresponds to the nominal unit normal vector field. Combined with the interface micromorphology perturbation height reconstructed in step S11 Introducing the orthogonal surface projection operator By mapping the three-dimensional spatial coordinates to the nearest wall reference point, the implicit surface equations of the real fluid-structure interaction interface, including random process roughness perturbations, are constructed:

[0143] ;

[0144] In the formula, Represents the implicit surface function value of the reconstructed microchannel fluid-structure interface, which is dimensionless; The sign distance function value representing the nominal boundary, with dimensions of ; This represents a dimensionless mapping operator that projects three-dimensional spatial coordinates onto a nominal surface. This represents the height of the corresponding microscopic morphological perturbation, with dimensions of . By taking the partial derivatives of the implicit surface equation with respect to spatial coordinates, the spatial gradient manifold vector of the surface equation is obtained. :

[0145] ;

[0146] In the formula Represents the gradient operator in three-dimensional space, with dimensions of ; Denotes the tangential gradient operator along the nominal surface, with dimensions of ; The Jacobian matrix representing the projection operator is dimensionless. The true local unit normal field for spatial deflection is derived using the spatial gradient manifold vector. Its mathematical expression is:

[0147] ;

[0148] In the formula, Represents the true local unit normal vector, which is dimensionless; Represents the L2 norm operator. When the root mean square height of the surface... Under ideal smooth conditions approaching zero, the microscopic morphology disturbance is highly Always equal to 0, spatial gradient manifold vector Automatically degenerates to the nominal normal vector True local unit normal vector With the nominal normal vector Full alignment ensures the physical and logical completeness and dimensional consistency of boundary reconstruction.

[0149] S312. To address the potential velocity field discontinuity fitting problem caused by microscopic slippage at the fluid-structure interface, numerical compliance verification is established by clarifying the differentiability characteristics of the neural network within the fluid domain closure. The fluid velocity field tensor output by the neural network... Within and along the boundary of the fluid computational domain, the system is composed of nested, continuously globally differentiable Swish activation functions, whose mathematical form naturally maintains analytical differentiability at the spatial boundary. Microscopic interface slip is essentially a finite-value jump in boundary velocity at the fluid-solid interface, rather than a break in the manifold of the continuous velocity function within the fluid domain, thus satisfying the prerequisites for analytical applications of automatic differentiation. Utilizing the built-in automatic differentiation engine of the physical information neural network, the extracted fluid velocity field tensor is... Find the Jacobian matrix with respect to spatial coordinates. The Jacobian matrix With the real local unit normal vector field Perform tensor dot product operations to calculate the gradient distribution operator tensor of the velocity field along the local normal of the solid wall of the microchannel. The specific mathematical formula is as follows: In the formula, It is a second-order tensor with dimensions of . ; It is a dimensionless column vector; Let be the absolute rate of change vector of the velocity field along the normal direction, with dimensions . .

[0150] S313. Due to the local deflection of the normal vector of a real microscopic rough interface, the directly calculated normal gradient distribution operator tensor will contain a mixture of normal and tangential components. Construct a tensor from a second-order unit tensor... With normal projection tensor Tangential projection operator tensor formed by subtraction Its expression is In the formula It is a second-order unit tensor, dimensionless; This represents the tensor outer product operator. The inner product operation is performed between the tangential projection operator tensor and the normal gradient distribution operator tensor to remove the spurious normal slip driving force caused by geometric distortion, generating a pure tangential gradient driving tensor. The specific mathematical formula is as follows: In the formula, This represents the pure tangential gradient tensor after filtering out the normal component, with dimensions of . By introducing this tensor inner product operation, it is ensured that the boundary slip velocity tensor calculated subsequently strictly adheres to the tangential plane of the locally rough solid wall, eliminating the normal penetration velocity and satisfying the rigid physical constraint condition of impenetrable solid walls in fluid mechanics. The final extracted tangential gradient driving tensor... Directly used as the boundary tangential slip velocity tensor Input variables.

[0151] S32. Construct a spatially distributed and temperature-dependent local slip coefficient tensor. Its expression is:

[0152] ;

[0153] In the formula, Dimensions are ; This represents the baseline value for molecular dynamics fitting of a smooth interface; Represents the roughness correction function, where The root mean square height of the surface. The average slope; This is the roughness attenuation coefficient, with dimensions of... ; Indicated by substrate ambient temperature drift The input temperature correction function is dimensionless. This represents the interfacial viscosity temperature sensitivity coefficient, with dimensions of . .

[0154] In this embodiment, the molecular dynamics fitting benchmark values ​​for the smooth interface are solved independently in advance using a large-scale molecular dynamics simulation platform. A microscopic simulation box was constructed, containing a solid-walled monocrystalline silicon matrix and deionized water molecules. The Tersoff many-body potential function was used between the silicon atoms in the solid-walled matrix, the SPC / E water model was used for the deionized water molecules, and the Lennard-Jones potential function and Lorentz-Bertlot mixing rule were used for the fluid-solid interface interaction. The spatial dimensions of the microscopic simulation box were set to [value missing]. At the target standard reference temperature Next, configure the Nosé-Hoover hot tub for... The NVT canonical ensemble was relaxed to reach initial thermodynamic equilibrium; then the heat bath controller was suspended and the system was switched to NVE microcanonical ensemble operation. Obtain statistical production trajectories. Extract the force trajectories of fluid molecules at the solid-liquid interface. Based on Green-Kubo linear response theory, analyze the intrinsic interface friction coefficient by integrating the time autocorrelation function of the solid-liquid interface interaction forces.

[0155] ;

[0156] Based on the constitutive relation of the Nanoweig slip boundary and combined with hydrodynamic viscosity, the baseline values ​​for molecular dynamics fitting of smooth interfaces are calculated:

[0157] ;

[0158] In the formula, The intrinsic friction coefficient of the solid-liquid interface is expressed in terms of dimensions. ; The geometric projected area of ​​the solid-liquid interface of the microscopic simulation box is represented by the dimensionless geometrical projection area. ; Denotes the Boltzmann constant, with dimensions of ; Represents the absolute equilibrium temperature of the system, with dimensions of ; This represents the sum of instantaneous interaction forces exerted tangentially on a group of fluid molecules by a solid wall, with dimensions of . ; This represents the ensemble average operator, which is dimensionless. Represents the dynamic viscosity of a fluid, with dimensions of Time element Dimensions are Based on the above physical configuration extraction and system ensemble state determination, the baseline value for molecular dynamics fitting of the smooth interface between deionized water and the silicon-based surface in this embodiment was determined. .

[0159] For non-ideal processes and dynamic power consumption impact conditions, the stepped configuration of the solid-wall atomic layers is modified in molecular dynamics simulations to generate a series of surface root-mean-square height distributions. With average slope distribution in For the micro-geometric perturbation samples within the interval, the corresponding attenuation slip length is calculated by repeating the aforementioned NVE ensemble statistical process. The Levenberg-Marquardt nonlinear least squares method is then used to perform a global optimization fit on the micro-geometric sample dataset to extract the roughness attenuation coefficient. Its value range is determined as follows Simultaneously, considering the nonlinear evolution of the fluid-solid interface slip length with temperature, particularly the exponential response mechanism of the relaxation behavior of the hydrogen bond network in the confined water molecule layer at temperatures far from the reference temperature, an exponential nonlinear temperature correction function was constructed. In the simulation, the absolute equilibrium temperature of the system was set at... to Discrete multi-point sampling is performed within the interval to extract the corresponding interfacial viscosity attenuation rate. The interfacial viscosity temperature sensitivity coefficient is obtained by Arrhenius nonlinear fitting. Its value range is determined as follows This enables high-precision nonlinear dynamic correction of microscopic local slip coefficient tensors by temperature drift over a wide temperature range.

[0160] S33. Constructing the continuous weight tensor of the interface mask. The local roughness features are transformed into continuous features, and the expression is as follows:

[0161] ;

[0162] In the formula, It is a dimensionless scalar tensor; This represents the effective solid-liquid contact state coefficient caused by local roughness, and is dimensionless. This represents the reference contact coefficient, which is dimensionless. This represents the normalization function.

[0163] normalization function The generalized sigmoid smooth step function, after extreme value normalization and numerical stability correction, is used:

[0164] ;

[0165] Among them, the basic Sigmoid mapping function The mathematical definition formula is:

[0166] ;

[0167] In the formula, The tensor, representing the physical independent variable, is dimensionless, and its theoretical range of values ​​is strictly defined according to the surface geometric topological derivation. Within the range; This represents the set zero-valued stability constant, which is dimensionless and takes a constant value. ; This represents the coefficient of steepness of the interfacial wetting phase transition, which is dimensionless. This represents the slip degradation phase transition threshold, which is dimensionless.

[0168] To ensure that the boundary manifold accurately adapts to different fluid-solid material systems and eliminates the universality defects of fixed constant settings, the interface wetting phase transition steepness coefficient is... With the slip degradation phase transition threshold Based on the static equilibrium contact angle between the working fluid and the wall surface and the adsorption energy of the micro-interface, a dynamic correlation mapping is performed, and its calculation formula is as follows:

[0169] ;

[0170] ;

[0171] In the formula, This represents the set lower limit constant for the phase transition steepness. It is dimensionless, for example, a constant value of 2.0. It is used to construct the bottom-line nonlinear gradient in the evolution of the forward operator to prevent numerical dissipation due to excessively low steepness. and These represent the preset dimensionless proportional adjustment coefficients; It represents the static equilibrium contact angle of a specific fluid cooling medium on the surface of a specific solid wall material. It is dimensionless and the unit is radians. This represents the intrinsic interfacial adsorption energy of a single fluid molecule on a solid wall surface, extracted from equilibrium molecular dynamics simulations, with dimensions of... ; Denotes the Boltzmann constant, with dimensions of ; This represents the standard reference temperature, with dimensions of . When the material system switches to an extremely low adsorption energy or high contrast property region, the instantaneously calculated phase transition steepness coefficient changes. Approaching the lower limit and phase transition threshold The convergence to the boundary causes the difference between the fundamental constant and the benchmark term to... When extremely shrunk, the denominator passes through the maximum value operator. Introducing numerical stability constants Dynamic truncation is implemented to prevent the malignant amplification of singular values ​​and the backpropagation gradient explosion caused by an extremely small denominator; when the local fluid-solid contact region approaches the absolute roughness limit, the effective solid-liquid contact state coefficient is reduced. When, the independent variable , numerator Absolute zeroing, continuous weight tensor of interface mask Automatic convergence to scalar 0 indicates that the fluid in this region is completely viscously locked by the rough valleys, and the slip effect is completely degraded; when the local contact region tends to be ideally smooth and the denominator does not trigger truncation, the independent variable... The numerator and denominator are exactly equal. Strictly equal to scalar 1, representing the ideal slip contribution of the region with a fully released smooth interface. This composite construction scheme of extremum normalization and denominator stability protection endows the mask tensor with... The ability to perform high-order continuous nonlinear mapping within closed intervals ensures the physical and logical completeness and dimensional consistency of the fluid-structure interaction boundary tensor model under various heterogeneous phase boundary abrupt change conditions.

[0172] This embodiment introduces equilibrium contact angle and micro-adsorption energy terms, upgrading the characteristic parameters (steepness and threshold) controlling the evolution of slip step from manually specified parameters to adaptive variables driven by the underlying physical properties of the material. This allows the simulation model to exhibit universal physical correlations when facing different cooling media (such as deionized water and fluorocarbon electronic fluids) and different pipe wall surfaces (such as monocrystalline silicon and copper plating). Simultaneously, to address the ill-conditioned behavior of zero measure in the denominator that may be triggered by extreme value normalization in specific property ranges, a spatial operator protection based on maximum value limiting is constructed at the denominator end. This design mathematically cuts off the numerical divergence path when the baseline difference collapses due to parameter variations, ensuring that the mask gradient flow remains bounded and smooth during boundary condition loss backpropagation in the physical information neural network, thus improving the convergence success rate of co-simulation in extreme process deviation regions.

[0173] In this embodiment, the effective solid-liquid contact state coefficient Contact coefficient with reference Based on the microscopic geometric morphology features extracted in step S12, the following method is used: The irregular undulations of the solid wall morphology of the microchannel cause the actual solid-liquid contact area to be larger than the macroscopic projected area, thereby enhancing interfacial fluid-solid friction and weakening the slip effect. This is defined as follows: The ratio of the actual microscopic physical surface area to the geometrically projected area under local roughness is used, based on the local average slope extracted in step S12. Perform geometric differential relation analysis: Define the reference contact coefficient. The effective contact state coefficient is defined as follows: under ideal, absolutely smooth interface conditions, the solid-wall interface has no microscopic morphological disturbances or process defects, the local average slope is always zero, and the actual physical surface area is strictly equivalent to the geometrically projected area. The corresponding state coefficient is taken as an accurate physical constant. ;In the formula, Representing spatial coordinates The effective solid-liquid contact state coefficient at the point is dimensionless. This represents the local average slope obtained in step S12, which is dimensionless. This represents an absolute constant under reference conditions and is dimensionless. When actual semiconductor etching processes approach an ideal zero-defect state, micro-roughness characteristics disappear, and the local average slope... Automatic convergence to Effective contact state coefficient Automatically degenerates into extreme values , equal to the reference contact coefficient At this point, the ratio term in step S33 The calculation result converges to This makes the interface mask continuous weight tensor Approaching full weight The slip model degenerates non-destructively into the basic slip formula for an ideal smooth wall. However, in regions where localized deep silicon etching anomalies lead to extreme roughness, the local average slope... ,drive Rapidly magnified, continuous weight tensor of interface mask Automatic exponential decay to The local slip term is completely truncated numerically, and the physical state exhibits a completely blocked microfluidic state.

[0174] S34. Computational fluid boundary tangential slip velocity tensor at the fluid-structure interaction interface on the inner wall of a microchannel. Its expression is:

[0175] ;

[0176] In the formula, Dimensions are Under extreme conditions, such as a power outage causing the fluid to become completely still, resulting in a zero normal gradient operator tensor, or extremely high local roughness leading to... or When approaching zero, The method automatically converges to the zero tensor, eliminating the normal penetration velocity and ensuring the physical logic completeness and dimensional consistency of the formula. It enables the roughness matrix to directly modulate the boundary slip effect, achieving a strong physical-logical coupling between geometric reconstruction and the slip boundary model.

[0177] S4. Substitute the microscopic interface slip boundary tensor model with process roughness linkage characteristics into the neural network for dynamic calibration, and analyze the nonlinear coupling to output the tensor manifold. For example... Figure 3 As shown, the specific sub-steps are as follows:

[0178] S41. Substitute the microscopic interface sliding boundary tensor model established in step S3 into the boundary tensor interface of the total loss function, instantiate the boundary condition residuals, and use the automatic differentiation mechanism to analyze the gradient propagation of the dynamic boundary. Specifically, this includes the following sub-steps:

[0179] S411. Obtain the boundary tangential slip velocity tensor output in step S34. Introducing a system characteristic reference velocity scalar. right Perform dimensionless processing and calculate the dimensionless target slip velocity tensor. : ;In the formula, The boundary tangential slip velocity tensor has dimensions of . ; The reference velocity represents the macroscopic characteristics of the fluid, with dimensions of ; Represent the dimensionless target slip velocity tensor. (This refers to the tensor of the dimensionless target slip velocity.) The value is assigned to the dimensionless target slip velocity tensor interface reserved in the boundary condition residual term of the total loss function in step S231. This enables the mapping of data states from physical boundary models to numerical solution frameworks.

[0180] S412, Instantiate the fluid boundary operator tensor at the fluid-structure interaction interface The discrete norms of the predicted velocity field tensor and the target slip velocity tensor are used as loss constraints to construct an instantiated boundary condition residual tensor. The specific mathematical formula is as follows:

[0181] ;

[0182] In the formula, Represents the total number of configuration nodes sampled on the geometric domain boundary of the reconstructed system, dimensionless; Nodes representing the output of a neural network The dimensionless fluid prediction velocity field tensor at the location; Represents a node The dimensionless target slip velocity tensor at the location; It represents the set of other dimensionless boundary operator tensors, including fluid inlet / outlet boundaries, wall heat flux boundaries, and constant voltage boundaries. This represents the L2 norm operator. In extreme conditions, such as when the fluid is completely still, the dimensionless target slip velocity tensor... When it converges to the zero tensor, the instantiated boundary condition residual tensor It automatically transforms into a forced zero constraint of a traditional non-slip solid wall boundary, eliminating the risk of non-physical slip and ensuring the physical and logical completeness and dimensional consistency of the boundary control equations under various operating conditions.

[0183] S413. Invoke the automatic differential calculation graph of the physical information neural network to track the boundary condition residual tensor. The nonlinear gradient propagation path in the network is addressed. A computational graph mechanism is used to transform the boundary slip fluctuations, which dynamically evolve with spatial coordinates, time, local roughness morphological parameters, and substrate environmental temperature drift, into gradient partial derivatives propagating backward from the network output layer to the hidden layer. This involves using the instantiated boundary condition residual tensor... Substitute into the dimensionless total loss function The forward and backward propagation calculations are performed to ensure that the weight iteration of the neural network is strictly controlled by the real microscopic interface thermal flux coupling mechanism.

[0184] S42. An adaptive gradient descent optimization algorithm, incorporating dynamic estimations of first and second moments, is invoked to iteratively reconstruct the network weights and bias parameters of the physical information neural network. A hybrid optimization strategy is configured to drive the neural network output space to approximate the real multi-field coupled manifold solution. The specific implementation steps are as follows:

[0185] S421. Initialize the dimensionless weight parameter matrix of the neural network. A cosine annealing dynamic learning rate decay function driven by the number of iteration steps is configured. The preceding network construction stage has utilized system physical characteristic benchmarks (fluid macroscopic characteristic reference velocity). Maximum operating temperature difference With reference supply voltage A strongly coupled dimensionless process is performed on all state variables, and the output range and partial differential residual operator of each independent physical field are analytically mapped to... Order-of-magnitude space. This normalization mapping mechanism eliminates the numerical cross-level bias caused by the inherent dimensional differences between velocity, pressure, temperature, and potential, giving the subsequently set unified global tolerance threshold a clear and consistent mathematical convergence meaning. The cosine annealing decay formula is:

[0186] ;

[0187] In the formula, Indicates the first The dynamic learning rate for each iteration cycle is dimensionless. This represents the minimum learning rate lower bound set, and is dimensionless. This represents the initial basic learning rate, which is dimensionless. This represents the current optimization iteration step, and is dimensionless. This represents the maximum set number of iteration cycles, dimensionless. Under extreme iterative conditions, when optimizing the number of iterations... Approaching the maximum total number of iteration cycles At that time, cosine term Automatic convergence to Drive dynamic learning rate Smooth asymptotically to nonzero minimum This prevents parameter out-of-bounds oscillations triggered by multi-field coupling approaching the end, and ensures the asymptotic stability and dimensional consistency of the convergent trajectory.

[0188] S422. Use the adaptive gradient descent algorithm based on first-order and second-order momentum estimation to solve for the total loss function. For the weight parameter matrix The nonlinear gradient propagation is used to perform iterative updates of the network parameters. The specific mathematical formula is:

[0189] ;

[0190] In the formula, Indicates the first The network parameter tensor after the next iteration; Indicates the first The network parameter tensor for the next iteration; This represents the gradient first moment estimate tensor of the total loss function after bias correction; This represents the gradient second moment estimation tensor of the total loss function after bias correction; Represents a minimal positive scalar that prevents division by zero errors; This represents the Hadamard product operator; This represents the element-wise square root operator of the tensor. In the extreme convergence case where the physical model evolves to an absolutely steady state, the backpropagation gradient of the total loss function approaches the zero tensor, and the first-moment estimate of the tensor... As the matrix decays to zero, the parameter update step size automatically converges to zero, avoiding invalid oscillations of network parameters near the optimal physical manifold solution and ensuring the consistency of the physical convergence logic and dimensions of the formula.

[0191] S423. Establish a two-stage hybrid optimization switching mechanism based on multi-dimensional monitoring and momentum information inheritance. During the iterative process using the adaptive gradient descent algorithm, continuously monitor the continuous... The relative smooth rate of change of the total loss function within each iteration period Simultaneously calculate the infinite norm of the physical residual components of each partial differential equation. A finite-memory quasi-Newton optimization algorithm (L-BFGS) is introduced to replace the full Hessian matrix update strategy, and the historical vector storage size parameter is set. Store only the most recent The curvature information vector of the step avoids the risk of memory overflow caused by the second-order partial derivative matrix with tens of millions of parameters.

[0192] when Below the preset pseudo-Newton switching tolerance threshold And the norm of all independent residual components All are below the preset residual pre-convergence threshold. When the algorithm switches criteria, it suspends the adaptive gradient descent algorithm. Extract the suspension time. Network parameter tensors First-order momentum estimation tensor With the second-order momentum estimation tensor The initial inverse Hessian approximation matrix of the finite-memory quasi-Newton optimization algorithm is constructed using the extracted momentum tensor. The diagonal elements are then subjected to curvature warm-start to prevent sudden breaks in the historical features of previous optimizations. The reconstruction baseline formula is as follows:

[0193] ;

[0194] In the formula, The diagonal components of the initial inverse Hessian approximation matrix for the finite-memory quasi-Newton optimization algorithm are dimensionless. The dimensionless total loss function gradient second moment estimation tensor represents the suspension time. The dimensionless gradient first moment estimation tensor of the total loss function at the suspension time is represented. This represents a minimal positive constant scalar that is protected against division by zero. If the evolution of the physical model encounters a locally extreme smooth gradient that causes the first moment to approximately zero, this diagonal matrix relies on an approximate smoothing constant. Forced truncation prevents a singular surge in the initial search step size, ensuring dimensional homogeneity and curvature transition stability within the dimensionless optimization framework. The finite-memory quasi-Newton optimization algorithm performs deep fine-tuning of the multiphysics residual space based on the reconstructed initial matrix.

[0195] S43. In the second-order fine-tuning iteration stage of the finite-memory quasi-Newton optimization algorithm (L-BFGS), continuously monitor the limit absolute convergence state of the physical residual tensor components of the partial differential equation on the globally discrete configuration nodes. Discarding the residual passivation defect caused by using the global mean square error, extract the term with the maximum absolute value of the partial differential residual on all internal configuration nodes to construct a global absolute maximum error metric:

[0196]

[0197] In the formula, Represents the global absolute maximum error within the entire physical computational domain; it is dimensionless. It represents the set of independent partial differential residual operator components, including the electromagnetic field equation, the mass conservation equation, the momentum conservation equation, and the energy conservation equation. This represents the total number of discrete, unlabeled configuration nodes within the physical computing domain. Represents a node The first dimensionless mapping after the system physical characteristic reference is applied. Each physical residual operator term is strictly equal to the physical residual operator term under the extreme condition of absolute thermodynamic equilibrium without external source excitation. Global absolute maximum error Simultaneous convergence to the theoretical zero value ensures the physical completeness of the error norm measure. When the global absolute maximum error... Strictly reduce to the set dimensionless global physical tolerance threshold When the relative rate of change within a preset number of iterations (≥50 steps) falls below a preset stagnation threshold (i.e., the peak residual enters a long-tailed stationary decay period), the simulation termination condition is triggered to stop iterative reconstruction. The dimensionless state variables predicted by the independent tensor output branches of the physical information neural network are extracted, inversely mapped using feature references, and then integrated to generate a multiphysics distribution tensor solution, thus completing the simulation.

[0198] The preferred embodiments of the invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention.

Claims

1. A multiphysics co-simulation method for integrated circuits based on path tensor manifold analysis, characterized in that, include: Acquire integrated circuit physical layout data, material properties and fluid cooling fluid properties, introduce surface roughness tolerance matrix to construct multiphysics simulation geometric domain, and fuse dynamic working condition time series set to generate spatiotemporal input tensor; The spatiotemporal input tensor is input into the physical information neural network to output a multiphysics state tensor, and a total loss function is constructed by weighting the data fitting residual tensor, the partial differential equation physical residual tensor, and the boundary condition residual tensor. Based on the fluid velocity field tensor in the multiphysics state tensor, the tangential gradient driving tensor is extracted, and the boundary tangential slip velocity tensor is calculated by combining the local slip coefficient tensor and the interface mask continuous weight tensor. Substitute the boundary tangential slip velocity tensor into the boundary condition residual tensor, optimize the parameters of the physical information neural network based on the total loss function optimization strategy, and output the multiphysics distribution tensor solution.

2. The integrated circuit multiphysics co-simulation method based on path tensor manifold analysis according to claim 1, characterized in that, Acquire integrated circuit physical layout data, material properties, and fluid cooling medium attributes; introduce a surface roughness tolerance matrix to construct a multiphysics simulation geometric domain; and fuse dynamic operating condition time series sets to generate a spatiotemporal input tensor, including: Extract the initial three-dimensional geometric boundary of the microchannel liquid cooling structure and establish the nominal boundary coordinate set and the nominal unit normal vector field; The surface roughness tolerance matrix is ​​analyzed to construct the interface micromorphology perturbation field. The interface micromorphology perturbation field is superimposed on the nominal boundary coordinates to generate a set of geometric domain space coordinates for multiphysics simulation. A multi-physics dynamic testing platform is adopted. Based on the extracted minimum electromagnetic transient characteristic time constant of integrated circuit and the macroscopic equivalent thermal diffusion time constant, the time window configuration is defined, and the dynamic operating condition time series set within the specified time window is acquired synchronously. The root mean square height and local average slope of the local surface at the set of spatial coordinates of the geometric domain of the multiphysics simulation are extracted. The root mean square height of the local surface, the local average slope and the dynamic working condition time series set are mapped to the corresponding spatial grid nodes and time steps and tensor splicing is performed to generate the spatiotemporal input tensor.

3. The method for multiphysics co-simulation of integrated circuits based on path tensor manifold analysis according to claim 2, characterized in that, The surface roughness tolerance matrix is ​​analyzed to construct the interface micromorphological perturbation field, including: Using the Karhunen-Loève expansion method, the integral domain is discretized by the quadrature rule and the Gaussian autocovariance function corresponding to the surface roughness tolerance matrix is ​​solved to obtain the eigenvalue sequence and the orthogonal eigenfunction set. Define the eigenvalue energy accumulation contribution rate function, and continuously accumulate the eigenvalues ​​until the energy accumulation contribution rate reaches the preset energy retention threshold to determine the eigenvalue cutoff order; Based on the eigenvalue sequence and orthogonal eigenfunction set, an interface micromorphological perturbation field is constructed by generating an independent random variable sequence using the Monte Carlo sampling method.

4. The method for multiphysics co-simulation of integrated circuits based on path tensor manifold analysis according to claim 1, characterized in that, The spatiotemporal input tensor is input into the physical information neural network, which outputs a multiphysics state tensor. A total loss function is constructed, weighted by the data fitting residual tensor, the partial differential equation physical residual tensor, and the boundary condition residual tensor, including: High-dimensional features are extracted from spatiotemporal input tensors by utilizing the feature-shared hidden layer configured with continuously differentiable nonlinear activation functions in the physical information neural network. The fluid velocity field tensor, fluid pressure field, continuous temperature field and electric potential field are output through independent tensor output branches, forming a multiphysics state tensor. Partial derivatives of the multiphysics state tensor are obtained, and the physical residual tensor of the partial differential equation is constructed by combining the reconstructed smooth spatial domain indicator function based on sign distance. The data fitting residual tensor is constructed by combining the real state sequence. The boundary condition residual tensor is constructed based on the reconstructed system geometric domain boundary state, the gradient step scaling operator of the fluid-structure interface, and the two-sided limit decoupling operator. The total loss function is generated by weighted summation of the data fitting residual tensor, the physical residual tensor of the partial differential equation, and the residual tensor of the boundary conditions.

5. The integrated circuit multiphysics co-simulation method based on path tensor manifold analysis according to claim 4, characterized in that, Generate the total loss function, including: During the training iteration cycle, the automatic differential computation graph is called to extract the maximum absolute value of the gradient of the data fitting residual tensor, the partial differential equation physical residual tensor, and the boundary condition residual tensor relative to the feature-shared hidden layer weight matrix. A gradient stability constant threshold and an adaptive weight hard upper limit cutoff value are introduced as stability boundary constraints. The ratio of the maximum absolute value of the gradient is calculated to determine the instantaneous target physical configuration weight value and the instantaneous target boundary condition weight value. Based on the set initial weight constant, the instantaneous target physical configuration weight value and instantaneous target boundary condition weight value are updated smoothly over time using the exponential moving average algorithm to generate an adaptive penalty weight tensor. The adaptive penalty weight tensor is summed as weighting coefficients to generate the total loss function.

6. The method for multiphysics co-simulation of integrated circuits based on path tensor manifold analysis according to claim 1, characterized in that, Extracting the tangential gradient driving tensor from the fluid velocity field tensor in the multiphysics state tensor, including: Based on the symbolic distance function of the nominal boundary and the surface projection operator, we construct the implicit surface equation of the geometric domain of the multiphysics simulation with roughness perturbation, calculate the spatial gradient manifold vector of the implicit surface equation, and generate the real local unit normal vector field. An automatic differentiation engine based on a physical information neural network is used to obtain the Jacobian matrix of the fluid velocity field tensor with respect to spatial coordinates. The Jacobian matrix is ​​multiplied by the real local unit normal field to generate a normal gradient distribution operator tensor. Construct a tangential projection operator tensor consisting of the subtraction of the second-order unit tensor and the normal projection tensor; The normal component is removed by performing an inner product operation between the tangential projection operator tensor and the normal gradient distribution operator tensor to generate the tangential gradient driving tensor.

7. The method for multiphysics co-simulation of integrated circuits based on path tensor manifold analysis according to claim 1, characterized in that, The boundary tangential slip velocity tensor is calculated by combining the local slip coefficient tensor and the interface mask continuous weight tensor, including: The local slip coefficient tensor is generated by multiplying the smooth interface molecular dynamics fitting benchmark value, roughness correction function and exponential nonlinear temperature correction function obtained based on the ensemble statistics of micro molecular dynamics. The effective solid-liquid contact state coefficient is calculated based on the local average slope. The ratio of the effective solid-liquid contact state coefficient to the reference contact coefficient is substituted into the generalized sigmoid smooth step function with denominator lower bound protection. The wetting phase transition steepness and degradation phase transition threshold of the step function are dynamically mapped by the solid-liquid equilibrium contact angle and intrinsic interface adsorption energy to generate the interface mask continuous weight tensor. Calculate the Hadamard product of the local slip coefficient tensor, the tangential gradient driving tensor, and the interface mask continuous weight tensor to generate the boundary tangential slip velocity tensor.

8. The method for multiphysics co-simulation of integrated circuits based on path tensor manifold analysis according to claim 1, characterized in that, Substituting the boundary tangential slip velocity tensor into the boundary condition residual tensor includes: By introducing a reference velocity based on the macroscopic characteristics of the fluid, the boundary tangential slip velocity tensor is made dimensionless to generate a dimensionless target slip velocity tensor. The target slip velocity tensor is assigned to the target slip velocity tensor interface reserved in the total loss function. The L2 norm of the difference between the fluid velocity field tensor and the target slip velocity tensor is used as the loss constraint to construct the instantiated boundary condition residual tensor.

9. The method for multiphysics co-simulation of integrated circuits based on path tensor manifold analysis according to claim 1, characterized in that, The parameters of the physical information neural network are optimized based on the total loss function, and the output multiphysics distribution tensor solution is provided, including: Configure a cosine annealing dynamic learning rate decay function driven by the number of iteration steps; An adaptive gradient descent algorithm based on the correction of the first-order momentum estimation tensor and the second-order momentum estimation tensor is invoked to solve the nonlinear gradient transfer of the total loss function on the network parameter tensor and update the network parameter tensor. The algorithm monitors the relative smoothing rate of change of the total loss function and the infinite norm of the physical residual components of each partial differential equation within continuous iteration cycles. When the relative smoothing rate of change of the total loss function is lower than the quasi-Newton switching tolerance threshold and the infinite norm of the residuals is lower than the preset residual pre-convergence threshold, the adaptive gradient descent algorithm is suspended. The network parameter tensor is passed to the finite memory quasi-Newton optimization algorithm module, and the diagonal elements of the initial inverse Hessian approximation matrix are constructed using the extracted momentum estimation tensor for curvature warm-start. The second-order approximate Hessian matrix is ​​then solved for residual space fine-tuning.

10. The method for multiphysics co-simulation of integrated circuits based on path tensor manifold analysis according to claim 9, characterized in that, The residual space fine-tuning process in solving the second-order approximate Hessian matrix using the finite-memory quasi-Newton optimization algorithm module also includes: Extract the term with the largest absolute value of the partial differential residual on all internal configuration nodes, construct a global absolute maximum error metric, and monitor the global absolute convergence state of the physical residual tensor components of the partial differential equation. When the global absolute maximum error strictly decreases below the set global physical tolerance threshold, the simulation termination condition is triggered to stop the iterative reconstruction. The multiphysics state tensors predicted by the output branches of each independent tensor of the neural network are extracted from the physical information, and then integrated after feature benchmark inverse mapping to generate the multiphysics distribution tensor solution.

Images