Multi-physics field coupling numerical calculation method based on multiplexing shift Laplace preprocessing acceleration

Abstract

The invention discloses a multi-physics field coupling numerical calculation method based on multiplexing shift Laplace preprocessing acceleration. In a multi-physical field coupling calculation process induced by a frequency domain electromagnetic field, a group of electromagnetic linear systems discretized by a finite element method which changes slowly needs to be solved numerically, and generally, time and memory are consumed very much. The method can effectively reduce the number of times of inversion operation of the preprocessing matrix, and significantly reduces the demand for computing resources. Different from a traditional algorithm using approximate inverse of a preprocessing matrix, the method provided by the invention constructs accurate inverse of the preprocessing matrix, and improves the stability and convergence rate of an iteration method, thereby ensuring high efficiency and high robustness of linear system solving. According to the method, the shift Laplace preprocessing technology is introduced into the multi-physics field coupling calculation for the first time, and the electromagnetic linear system sequence in the multi-physics field coupling simulation can be efficiently calculated.

Description

Technical Field

[0001] The present invention belongs to the field of multi-physics field numerical calculation, and in particular relates to a multi-physics field coupling numerical calculation method based on multiplexed shifted Laplace preprocessing acceleration. Background Art

[0002] As power electronic circuits, radio frequency integrated circuits, and three-dimensional integrated microsystems continue to develop toward miniaturization, high integration, high frequency, and high power consumption, more and more functional units and devices are being integrated onto smaller chip areas. This leads to ever-increasing power density in circuits and microsystems, and self-heating becomes a significant challenge. During the design phase, coupled electromagnetic (EM)-thermal-stress simulation is required for electronic devices, circuits, and even entire microsystems. The finite element method (FEM) excels at handling complex structures and models with heterogeneous material properties, making it well-suited as a multiphysics numerical algorithm for electronic devices, circuits, and even entire microsystems.

[0003] However, applying the finite element method (FEM) to computational problems involving increasingly complex and large-scale real-world systems is extremely challenging. On the one hand, the linear system dimensions formed after finite element discretization can be extremely large. On the other hand, while the Poisson equation describing heat conduction and thermal stress can be efficiently solved using multigrid methods, frequency-domain electromagnetic problems are described by the time-harmonic Helmholtz equation. Specifically, the system matrix obtained from finite element discretization is a high-dimensional complex matrix with non-diagonally dominant and symmetric non-positive definite values. Furthermore, this matrix is a non-M matrix, making multigrid methods less effective. For integrated electromagnetic-thermal numerical calculations, each physical equation must be solved separately and iterated repeatedly until a steady state is reached. This means that multiple equations must be solved in a single numerical run. Therefore, the efficiency of electromagnetic problem solutions is a major bottleneck restricting the performance of multiphysics simulations. Current mainstream electromagnetic simulation software uses direct methods such as MUMPS and SuperLU, which are the default settings for robust solutions. However, due to computational complexity and memory limitations, the computational scale and efficiency of direct methods are unsatisfactory. Consequently, iterative solutions to the time-harmonic Helmholtz equation have become a research hotspot in recent years. Considering the poor numerical properties and large condition number of the system matrix of electromagnetic fields, it is crucial to develop effective preconditioning methods for iterative methods.

[0004] In recent years, various preconditioning methods for the time-harmonic Helmholtz equation have emerged in diverse fields, such as computational mathematics and computational electromagnetics. These methods are based on the traditional domain decomposition method (DDM) or multigrid algorithm, and propose corresponding preconditioning methods based on the characteristics of different applications. Among them, preconditioning techniques based on the shifted Laplace theory are intuitive and easy to implement. The key to this method is the introduction of complex shifts to the Helmholtz operator, thereby constructing a computationally inexpensive preconditioning matrix. If the shift value is too large, the shifted Laplace operator will deviate from the original equation, resulting in inefficient preconditioning. Conversely, if the shift value is too small, the preconditioning matrix is very similar to the original matrix, making the inversion cost prohibitive. Currently, research on shifted Laplace preconditioning techniques has primarily focused on the field of computational mathematics, focusing on theoretical analysis of the selection of shift values and their integration with domain decomposition or multigrid methods, with few practical application examples.

[0005] In multi-physics coupled calculations, it is necessary to solve a set of slowly varying electromagnetic field finite element matrices. Considering the characteristics of this problem and inspired by the shifted Laplace method, this paper proposes a reusable, shifted Laplace preconditioning method for coupled electromagnetic-thermal-stress calculations. A corresponding preconditioning solution framework is developed: a direct method is used to solve the first linear system in a sequence, and a shifted Laplace preconditioning matrix is constructed for subsequent linear systems. When the number of iterations exceeds a preset threshold, indicating a decline in preconditioning effectiveness, the direct method is automatically invoked to update the preconditioning matrix. The proposed method combines high robustness with high efficiency. Summary of the Invention

[0006] In the multi-physics field coupling calculation process under frequency-domain electromagnetic field excitation, it is necessary to solve a set of slowly changing electromagnetic linear systems, which is very time-consuming and memory-consuming. To address this problem, the present invention proposes a method for accelerating the numerical calculation of multi-physics field coupling based on multiplexed shifted Laplace preprocessing. This method can accelerate the numerical calculation of multi-physics field coupling such as electromagnetic-thermal-stress and reduce the demand for computing resources.

[0007] The technical solution adopted in the present invention is as follows:

[0008] The present invention proposes a reusable shifted Laplace preconditioning method and forms a solution framework for accelerating the numerical computation of multi-physics coupling based on the reusable shifted Laplace preconditioning. This solution is used for the numerical computation of electromagnetic linear system sequences in multi-physics coupling. This method can effectively reduce the number of preconditioning matrix inversion operations in the preconditioning Krylov subspace method, significantly reducing the consumption of computing resources. Unlike traditional methods that construct an approximate inverse of the preconditioning matrix, the present invention constructs its exact inverse, achieving both high robustness and high efficiency. First, a direct method is used to solve the first electromagnetic finite element matrix. The resulting LU matrix decomposition is stored and used as the shifted Laplace preconditioning matrix for subsequent systems. The subsequent systems are then iteratively solved using the Krylov subspace method. If the number of iterations of the system exceeds a predefined threshold, the direct method is applied to the next linear system, and the preconditioning matrix is updated. This continues until all linear systems in the sequence have been solved.

[0009] Furthermore, the method may specifically include the following steps:

[0010] The first step is to establish a geometric model of the multi-physics field coupling system and perform initial meshing to obtain the mesh topology, namely, shape function information, unit adjacency, and boundary conditions.

[0011] Step 2: Load the grid topology and use the Galerkin method to perform weak formalization on the wave equation of the frequency-domain electromagnetic field. Then, use the vector finite element method to discretize it and obtain the linear system matrix of the electromagnetic field finite element.

[0012] Step 3: Determine whether it is the first time to solve the electromagnetic matrix, or whether the number of Krylov subspace iterations Iter is greater than the preset value. If so, proceed to step 4 to solve the shifted Laplace preconditioning matrix; otherwise, proceed to step 5 to iteratively solve the preconditioning linear system using the Krylov subspace method.

[0013] Step 4: Solve the electromagnetic matrix using the direct method and store the LU decomposition as the shifted Laplace preconditioning matrix in subsequent solutions;

[0014] Step 5: Use the pre-processing Krylov iterative method to solve the shifted Laplace pre-processed linear system and record the number of solution iterations Iter;

[0015] Step 6: Calculate the electromagnetic field dissipation power as the heat source for the heat conduction problem. Construct a heat conduction boundary value problem, load the mesh topology, and perform a weak formalization using the Galerkin method. Then, discretize the problem using nodal basis function finite elements to obtain the linear system matrix of the heat conduction finite element, which is then solved using the multigrid method. After the thermal solution is completed, update the temperature-dependent material parameters based on the temperature distribution.

[0016] Step 7: Determine whether thermal stability has been achieved by determining whether the maximum temperature difference obtained from two adjacent electromagnetic-thermal coupling solutions is less than a preset value. If thermal stability has not been achieved, repeat the electromagnetic-thermal coupling iteration process from steps 2 to 6. If thermal stability has been achieved, terminate the electromagnetic-thermal coupling iteration process. For the electromagnetic-thermal-stress coupling numerical calculation method, proceed to step 8 to calculate the thermal stress distribution.

[0017] Step 8: Calculate the temperature rise, construct a thermal stress boundary value problem, load the mesh topology, and perform a weak formalization using the Galerkin method. Then, discretize the problem using nodal basis function finite elements to obtain the linear system matrix of the thermal stress finite element, and solve it using the multigrid method.

[0018] Step 9: Post-processing output results.

[0019] In the above method, the temperature-dependent characteristics of material parameters bring nonlinear effects to the electromagnetic-thermal coupling problem, requiring iterative solutions for the electric and thermal fields until a thermal steady state is reached. The governing equations for the frequency-domain electric and steady-state thermal fields for the i-th iteration are as follows:

[0020]

[0021] where μ r is the relative magnetic permeability, E is the electric field intensity, k0 is the wave number in vacuum, ε r is the relative dielectric constant, T i-1 is the temperature distribution in the i-1th iteration step, Z0 is the free space wave impedance, J imp is the current density, k is the thermal conductivity, and P is the heat source.

[0022] It is necessary to repeatedly solve the partial differential equations (PDEs) of the electric and thermal fields and update the material parameters that change with temperature. For thermal field calculations, the governing equation is a Poisson equation, which can be efficiently implemented using the multigrid method; but for electric field calculations, the governing equation is a time-harmonic Helmholtz equation, which is the electromagnetic linear system K obtained by finite element discretization. i It is a large complex matrix with non-diagonal dominant and non-positive definite properties. Solving this matrix will consume a lot of memory and time. Moreover, the coupled iteration needs to solve the electromagnetic linear system sequence rather than a single linear system, that is,

[0023] K i e i =b i , i≥1

[0024] Among them, e i is the solution vector, b i is the right-hand term;

[0025] The present invention proposes a reusable preprocessing method that can accelerate the solution of electromagnetic linear system sequences. In the above electromagnetic-thermal coupling iterative process, although a set of slowly varying time-harmonic Helmholtz equations needs to be solved, the coefficients in the equations are only According to the shifted Laplace theory, the time-harmonic Helmholtz equation solved in the first step can be regarded as the shifted Laplace problem of all subsequent Helmholtz equations in the sequence. Therefore, the inverse of K0 can be used as the preconditioning matrix of the electromagnetic linear system sequence.

[0026] Considering the matrix sequence K i The similarity between i=0,1,2,..., the present invention performs LU decomposition on K0 in the fourth step and As a preprocessing matrix for subsequent problems; then, construct the preprocessing problem Use Krylov subspace method to solve iteratively. i When the gap between K0 and K0 is too large, the pre-processing effect will decline, and the number of iterations Iter will exceed the preset threshold value. In this case, the third step of the present invention will automatically jump to the direct method solution step to update the pre-processing matrix.

[0027] It should be pointed out that in the pre-conditioned Krylov subspace method, the most time-consuming and memory-intensive step is the matrix decomposition inversion operation of the pre-conditioned matrix. In order to improve computational efficiency, this method constructs a reusable pre-conditioned matrix. According to the judgment in the third step, only a few LU decomposition inversion operations need to be completed during the initial electromagnetic solution or the pre-conditioned matrix update process. Because the explicit forms of the L and U decompositions of the pre-conditioned matrix can be directly called, the remaining linear systems in the sequence can be efficiently solved using the Krylov subspace method. Therefore, the K in the Krylov iteration process is s x=y can be converted to L(Ux)=y, and the solution can be completed with extremely low computational cost through forward substitution and backward substitution.

[0028] The beneficial effects of the present invention are:

[0029] This invention effectively reduces the number of preconditioning matrix inversion operations, significantly lowering computational resource requirements. Unlike conventional algorithms that use an approximate inverse of the preconditioning matrix, this invention constructs an exact inverse of the preconditioning matrix, improving the stability and convergence speed of the iterative method, thereby ensuring high efficiency and robustness in solving linear systems. This invention introduces shifted Laplace preconditioning technology into multiphysics coupling calculations for the first time, enabling efficient calculation of electromagnetic linear system sequences in multiphysics coupling simulations. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 This is the specific process for accelerating the electromagnetic-stress-coupling calculation described in the present invention;

[0031] Figure 2 is the pseudo code of the novel reusable shifted Laplace preprocessing algorithm described in the present invention;

[0032] Figure 3 It is a schematic diagram of X-band low-temperature co-fired ceramic (LTCC) filter;

[0033] Figure 4 This is the S-parameter simulation result of the X-band low-temperature co-fired ceramic filter under thermal steady state;

[0034] Figure 5 This is the curve of the maximum temperature of the X-band low-temperature co-fired ceramic filter changing with frequency obtained by simulation;

[0035] Figure 6 This is the curve of the maximum Fan's stress of the X-band low-temperature co-fired ceramic filter obtained by simulation as a function of frequency;

[0036] Figure 7 This is a schematic diagram of the W-band packaged filter balun;

[0037] Figure 8 The comparison of solution time between the method proposed by the present invention and the traditional direct method under different power injection conditions is shown in FIG.

[0038] Figure 9 The memory consumption comparison between the method proposed in this invention and the traditional direct method is shown in Figure 2.

[0039] Figure 10 The surface temperature distribution of the filter balun in thermal steady state at 94 GHz and 14 W power injection obtained by simulation using the multi-physics field solution method of the present invention;

[0040] Figure 11 The surface displacement distribution of the filtering balun in thermal steady state at 94 GHz and 14 W power injection was simulated using the multi-physics field solution method described in the present invention. DETAILED DESCRIPTION

[0041] The technical solution of the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0042] In order to deal with the nonlinearity in the electromagnetic-thermal coupling process, it is necessary to iteratively solve the electric field and thermal field until a thermal steady state is reached. Figure 1 This is the operational flow of the present invention combining the shifted Laplace preprocessing method to accelerate iteration.

[0043] The first step is to establish a geometric model of the multi-physics coupling system and perform initial meshing to obtain the mesh topology, namely, shape function information, cell adjacency, and boundary conditions. The governing equations for the i-th electromagnetic-thermal coupling iteration are as follows:

[0044]

[0045] The following steps 2 to 7 implement electromagnetic-thermal field coupling iterative calculation;

[0046] Step 2: Load the mesh topology and use the Galerkin method for the frequency-domain electromagnetic wave equation of the i-th electromagnetic-thermal coupling iteration to obtain the weak form:

[0047]

[0048] Then, the vector finite element is used to complete the discretization and obtain the electromagnetic field finite element linear system K i e i =b i ;

[0049] Step 3: Determine whether it is the first time to solve the electromagnetic matrix, or whether the number of Krylov subspace iterations Iter is greater than the preset value. If so, proceed to step 4 to solve the shifted Laplace preprocessing matrix; otherwise, proceed to step 5 to use the Krylov subspace method to iteratively solve the preprocessed linear system. Figure 2 As shown in , the electromagnetic linear system sequence can be solved faster. In the above electromagnetic-thermal coupling iterative process, although a set of time-harmonic Helmholtz equations need to be solved, the coefficients in the equations are only It changes slowly with increasing temperature. According to the shifted Laplace theory, the inverse of K0 is an efficient shifted Laplace preconditioning matrix for the electromagnetic linear system sequence.

[0050] Step 4: Solve the electromagnetic matrix directly and store the LU decomposition as the shifted Laplace preprocessing matrix in the subsequent solution; considering the matrix sequence K i The similarity between i=0, 1, 2, ..., the present invention preprocesses the matrix in the fourth step to complete the LU decomposition, and As a preprocessing matrix for subsequent problems. Then construct the preprocessing problem The Krylov subspace method is used to iteratively solve the calculation. i When the gap between K0 and K1 is too large, the preprocessing effect will decrease, and the number of iterations Iter will exceed the preset threshold Θ, which can be set to 20 based on experience. In this case, the present invention uses the direct method to decompose K i+1 , to update the preprocessing matrix

[0051] Step 5: Krylov iteration method to solve the preprocessing problem Record the number of electromagnetic solution iterations Iter. In the pre-processing Krylov subspace method, the most time-consuming and memory-intensive step is the matrix decomposition inversion operation of the pre-processing matrix. like Figure 2 The PGMRES function in the algorithm. In order to improve the computational efficiency, the present invention constructs a reusable pre-processing matrix. The present invention only completes the direct LU decomposition inversion operation in the first step and the pre-processing matrix update step. Because the explicit form of the L and U decomposition of the pre-processing matrix can be directly called, the remaining linear systems in the sequence can be efficiently solved using the Krylov subspace method. In this process, the Krylov iteration Convertible to The latter can be solved at a very low computational cost through forward substitution and backward substitution operations.

[0052] Step 6: Calculate the electromagnetic field dissipation power as the heat source for the heat conduction problem. Construct a heat conduction boundary value problem, load the mesh topology, and perform a weak formalization using the Galerkin method. Then, discretize the problem using nodal basis function finite elements to obtain the linear system matrix for the heat conduction finite element. Solve the problem using the multigrid method. After the thermal field solution is completed, update the temperature-dependent material parameters based on the temperature distribution.

[0053] Step 7: Determine whether thermal stability has been achieved by determining whether the maximum temperature difference between two consecutive electromagnetic-thermal coupling solutions is less than a preset value. If thermal stability has not been achieved, repeat the electromagnetic-thermal coupling iteration process from steps 2 to 6. If thermal stability has been achieved, terminate the electromagnetic-thermal coupling iteration process and proceed to step 8 to calculate thermal stress.

[0054] Step 8: Calculate the temperature rise, construct a thermal stress boundary value problem, load the mesh topology, perform weak formalization using the Galerkin method, and then discretize using nodal basis function finite elements to obtain the linear system matrix of the thermal stress finite element, and solve it using the multigrid method.

[0055] Step 9: Post-processing output results.

[0056] by Figure 3 The X-band filter in the example is used as an example. The multi-physics field numerical calculation method proposed in this invention is used to perform electromagnetic-thermal-stress coupling simulations at 201 evenly distributed frequency points within the frequency range of 8 to 12 GHz. The input power is 10 W. Figure 4 、 5As shown in Figures 6 and 7, the numerical calculation results of the electromagnetic, thermal and stress fields are obtained respectively, which are in good agreement with the results obtained by the simulation software COMSOL Multiphysics. At each frequency point, the electromagnetic-thermal coupling calculation is first completed. When the thermal stability is reached, the thermal stress is calculated based on the temperature rise. This process requires solving multiple electric field finite element matrices. The use of the preprocessing method proposed in the present invention can accelerate the solution process. The following table compares the solution time of the method of the present invention and COMSOL Multiphysics under the same computing resource conditions. The proposed method significantly shortens the solution time.

[0057]

[0058] On the other hand, as the input power increases, the nonlinear effect in the multi-physics coupling process becomes more significant, which is manifested as an increase in the number of electromagnetic-thermal coupling iterations. The number of electromagnetic finite element matrices that need to be solved increases, which makes multi-physics simulation under high-power electromagnetic injection face considerable challenges. The reusable shifted Laplace preprocessing algorithm proposed in this invention can effectively solve this problem. Figure 7 The filter balun in the example is used as an example. As the input power increases, the number of coupling iterations will increase. Figure 8 As shown, with the traditional direct method, the computation time increases with the number of iterations, while with the proposed method, the computation time shows little change with the number of iterations. This result is primarily due to the reusable preprocessing matrix constructed in the present invention, which reduces the computation time to a few LU decomposition operations. The remaining EM solution process is a series of well-preprocessed problems that can be solved quickly and efficiently using iterative methods.

[0059] For further comparison, the same simulation settings and computing resources are used, 1W power input at 94GHz, 8 computing nodes on a Pluto2 cluster with a total of 128 CPU cores; the performance of the traditional direct method and the method proposed in this paper for solving electromagnetic linear system sequences in multi-physics field coupling is compared. Figure 9 As shown, the maximum memory consumption of the proposed method is the same as that of the direct method. During the first electromagnetic solution, the memory consumption curves of the two methods overlap. However, in the subsequent three electromagnetic solutions, the reusable shift preprocessing algorithm retains the LU decomposition from the first solution for preprocessing, maintaining the same memory consumption and completing the three electromagnetic solutions in a very short time. In contrast, the direct method requires four repeated matrix decomposition calculations, which are time-consuming and memory-intensive. Figure 10 and Figure 11 The surface temperature and displacement distribution of the filtering balun in thermal steady state at 94 GHz and 14 W power injection obtained using the method of the present invention are respectively given.

Claims

1. A multi-physics field coupling numerical calculation method based on multiplexed shifted Laplace preprocessing acceleration, characterized in that: include: First, the direct method is used to accurately solve the first electromagnetic linear system, and the decomposed LU matrix is stored. Then, the exact inverse obtained from the above solution is used for the shifted Laplace preprocessing of the subsequent electromagnetic linear systems. The Krylov subspace method is further used to iteratively solve the subsequent electromagnetic linear systems. If the number of iterations of the system exceeds a preset threshold, the direct method is used to solve the next linear system to update the preprocessing matrix until the entire sequence of electromagnetic linear systems is solved.

2. The multi-physics field coupling numerical calculation method based on multiplexed shifted Laplace preprocessing acceleration according to claim 1 is characterized in that: The specific steps include: Step 1: Establish a geometric model of the multi-physics coupling system and perform initial meshing to obtain the mesh topology, i.e., shape function information, cell adjacency, and boundary conditions. Step 2: Load the grid topology and use the Galerkin method to perform weak formalization on the frequency-domain electromagnetic wave equation. Then, use the vector finite element method to discretize the equation and obtain the linear system matrix of the electromagnetic finite element. Step 3: Determine whether it is the first time to solve the electromagnetic matrix, or whether the number of Krylov subspace iterations Iter is greater than the preset value. If so, proceed to step 4; otherwise, proceed to step 5. Step 4: Solve the electromagnetic matrix using the direct method and store the LU decomposition as the shifted Laplace preconditioning matrix in subsequent solutions; Step 5: Pre-process the electromagnetic matrix using the Krylov iterative method and record the number of electromagnetic solution iterations Iter; Step 6: Calculate the electromagnetic field dissipation power as the heat source for the heat conduction problem, construct a heat conduction boundary value problem, and load the grid topology. Use the Galerkin method for weak formalization, then discretize using the node basis function finite element method to obtain the linear system matrix of the heat conduction finite element. Use the multigrid method to solve, and after completing the thermal field calculation, update the temperature-dependent material parameters according to the temperature distribution. Step 7: Whether the maximum temperature difference obtained from two consecutive electromagnetic-thermal coupling numerical calculations is less than the preset value is used to determine whether thermal steady state has been reached. If thermal steady state has not been reached, proceed to step 2 and start the next electromagnetic-thermal coupling numerical calculation iteration step; If the thermal steady state is reached, the electromagnetic-thermal coupling iteration process ends.

3. The multi-physics field coupling numerical calculation method based on multiplexed shifted Laplace preprocessing acceleration according to claim 2 is characterized in that: The frequency domain electric field wave equation solved in the first step is regarded as the shift Laplace problem of the subsequent multiple frequency domain electric field wave equations. Therefore, the electromagnetic finite element matrix inverse solved in the first step is is used as the preconditioning matrix and the following preconditioning problem is numerically solved using the Krylov subspace iteration method: The preprocessed matrix can be reused to reduce the number of matrix inversion operations that consume the most time and memory in the iterative solution.

4. The multi-physics field coupling numerical calculation method based on multiplexed shifted Laplace preprocessing acceleration according to claim 2, characterized in that: The Krylov subspace method is used to solve the electromagnetic solution matrix after the i-th preprocessing. When the number of iterations Iter is greater than the preset threshold, the preprocessing matrix is automatically updated in the next step 5. The multi-physics field coupling numerical calculation method based on multiplexed shifted Laplace preprocessing acceleration according to claim 2, characterized in that: The multi-physics field coupling is a multi-physics field coupling calculation process induced by the frequency domain electromagnetic field.

6. The multi-physics field coupling numerical calculation method based on multiplexed shifted Laplace preprocessing acceleration according to claim 2, characterized in that: For the electromagnetic-thermal-stress coupling numerical calculation method, after reaching thermal steady state in step 7 and ending the electromagnetic-thermal coupling iterative process, step 8 is entered to calculate the thermal stress distribution; Step 8: Calculate the temperature rise, construct the thermal stress boundary value problem, load the mesh topology, use the Galerkin method, perform weak formalization, and then discretize using nodal basis function finite elements to obtain the linear system matrix of the thermal stress finite element and solve it using the multigrid method; Step 9: Post-process the output results.

Images