Hybrid algorithm based on classical discontinuous finite element and variational quantum linear solver

Abstract

The application discloses a hybrid algorithm based on a classical discontinuous finite element and a variational quantum linear solver, relates to the technical field of quantum-classical hybrid computing crossover, and comprises the following steps: adopting a discontinuous finite element method to perform spatial dispersion on fluid control equations, using a high-order basis function to describe internal solution of a unit and exchanging unit interface information through numerical flux to form a semi-discrete form of a common differential equation group; adopting an implicit time advancing format on the semi-discrete common differential equation group to convert the physical quantity updating problem of the current time step into a solving task of a high-dimensional sparse linear equation group; applying a generalized minimum residual method on the high-dimensional sparse linear equation group on a classical computing platform, constructing a set of orthogonal basis vectors of a Krylov subspace through an Arnoldi process, and projecting the original equation group to the low-dimensional subspace to obtain a dense linear system with reduced dimension; and constructing a variational quantum solving process of an adaptive noise-containing medium-scale quantum processor based on the low-dimensional dense linear system.

Description

Technical Field

[0001] This invention relates to the field of quantum-classical hybrid computing technology, and in particular to a hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver. Background Technology

[0002] Quantum-classical hybrid computing refers to a computing paradigm that leverages the complementary strengths of classical and quantum computing and executes them collaboratively under current noisy medium-scale quantum (NISQ) hardware conditions. This technology typically involves a classical processor handling the overall algorithm flow control, data preprocessing, and post-processing, while specific subtasks suitable for quantum acceleration (such as solving, optimizing, or eigenmapping high-dimensional linear systems) are offloaded to the quantum processor for execution. Through repeated classical-quantum information interaction, it achieves computational performance that surpasses purely classical methods or solves problems that were previously difficult to handle. It serves as a crucial bridge connecting the capabilities of current quantum hardware with practical application needs.

[0003] With the increasing density of grids and the increase in polynomial order, the dimensionality of linear systems derived from implicit time progression grows dramatically, causing memory consumption and computational complexity to exceed the carrying capacity of current high-performance computing platforms. This severely restricts the efficient simulation of high Reynolds number and multi-scale flow problems. Furthermore, existing quantum linear solution patents mostly focus on solving abstract matrices, without considering the sparse structure, physical constraints, and condition number characteristics of linear systems derived from conservation laws in CFD. Directly applying these patents can easily lead to problems such as coding redundancy, gradient vanishing, or solutions without physical meaning. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver to address the problem that as mesh refinement and polynomial order increase, the dimensionality of linear systems derived from implicit time progression grows sharply, leading to memory consumption and computational complexity exceeding the capacity of current high-performance computing platforms. This severely restricts the efficient simulation of high Reynolds number and multi-scale flow problems. Furthermore, existing quantum linear solver patents mostly focus on solving abstract matrices, without considering the sparse structure, physical constraints, and condition number characteristics of linear systems derived from conservation laws in CFD. Directly applying these solutions can easily lead to problems such as coding redundancy, gradient vanishing, or solutions without physical meaning.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution: In a first aspect, the present invention provides a hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver, comprising: The fluid control equations are spatially discretized using the discontinuous finite element method. Higher-order basis functions are used to describe the solutions within the elements, and numerical flux is used to process the information exchange at the element interfaces, forming a semi-discrete set of ordinary differential equations. An implicit time-progression scheme is adopted for the semi-discrete ordinary differential equation system, which transforms the problem of updating physical quantities at the current time step into the task of solving a high-dimensional sparse linear equation system. The generalized minimum residual method is applied to the high-dimensional sparse linear equation system on a classical computing platform. A set of orthogonal basis vectors of the Krylov subspace is constructed through the Arnoldi process, and the original equation system is projected onto the low-dimensional subspace to obtain a dense linear system with reduced dimension. Based on low-dimensional dense linear systems, a variational quantum solution process adapted to noisy medium-scale quantum processors is constructed. The coefficient matrix is ​​encoded as a linear combination of Pauli operators, and parameterized quantum circuits are designed to prepare trial solution states. Based on the local loss function defined by the parameterized quantum circuit, the gradient is calculated by combining the Hadamard test with the parameter translation rule, which drives the classical optimizer to iteratively update the circuit parameters until the loss function converges, and outputs the quantum state corresponding to the optimal parameters. Quantum state tomography is performed on the quantum state to reconstruct a low-dimensional solution in classical vector form. The solution is then mapped back to the original high-dimensional space using the aforementioned orthogonal basis vectors to obtain the physical quantity correction at the current time step. The physical quantity correction is superimposed on the solution of the previous time step to complete the implicit update of the current time step, thereby achieving high-precision advancement of the flow field variables.

[0007] As a preferred embodiment of the hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver described in this invention, the following steps are taken: The fluid control equations are spatially discretized using the discontinuous finite element method, the internal solutions of the elements are described using higher-order basis functions, and the information exchange at the element interfaces is processed through numerical flux to form a semi-discrete system of ordinary differential equations. The computational region is divided into several non-overlapping computational units; Within each computational unit, a set of higher-order polynomial basis functions is introduced. These basis functions use the average value of the physical quantities at the unit center as the zero-order term, and the remaining higher-order terms are constructed from the offset of the local coordinates relative to the unit center. The conserved variables in the fluid control equations are expanded in each cell according to the higher-order basis functions and substituted into the weak-form integral expression; Perform integration by parts on the spatial derivative terms involved in the weak form, transfer the derivatives in the volume integral to the basis functions, and at the same time explicitly separate the flux integral terms on the unit boundary; At the interface shared by adjacent cells, the continuity of solutions on both sides is not forced. Instead, the numerical flux is calculated using the Roe or Lax-Friedrichs scheme, which depends only on the reconstructed values ​​of conserved variables provided by the cells on both sides of the interface. After completing the weak-form assembly of all elements, a system of semi-discrete ordinary differential equations is obtained, consisting of the time derivatives of the basis function coefficients of each element.

[0008] As a preferred embodiment of the hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver described in this invention, wherein: the semi-discrete ordinary differential equation system adopts an implicit time-progression scheme, transforming the physical quantity update problem of the current time step into the task of solving a high-dimensional sparse linear equation system, the specific steps of which are as follows: Select an implicit time discretization scheme, including but not limited to the first-order backward Euler scheme or the second-order implicit Runge-Kutta scheme; Integrate the semi-discrete ordinary differential equation system in the current time step, so that the time derivative term is expressed as the difference between the solution in the current time step and the solution in the previous time step divided by the time step size. The spatial residual term is linearized at the current time step, and the first derivative term with respect to the current solution is retained to form the Jacobian matrix. The coefficient matrix of the linear system is formed by multiplying the mass matrix by the inverse of the time step and then superimposing the results into the Jacobian matrix. The result of interacting the residual term of the previous time step with the mass matrix is ​​used as the right-hand side term of the linear system. This leads to the construction of a high-dimensional sparse linear system of equations with the unknown solution at the current time step as the variable.

[0009] As a preferred embodiment of the hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver described in this invention, the steps of applying the generalized minimum residual method to the high-dimensional sparse linear equation system on the classical computing platform, constructing a set of orthogonal basis vectors of the Krylov subspace through the Arnoldi process, and projecting the original equation system onto this low-dimensional subspace to obtain a dense linear system with reduced dimensionality are as follows: The initial vector for initializing the Krylov subspace is the direction after the right-hand side is normalized; Perform Arnoldi iterations sequentially: In each step, apply the coefficient matrix to the current basis vector, and then subtract its projection components on the existing basis vectors from the result vector to achieve Gram-Schmidt orthogonalization. Normalize the orthogonalized vectors and add them to the subspace as the next basis vectors; The projection coefficients generated during each orthogonalization process are recorded synchronously to form the upper Hessenberg matrix; The Arnoldi process terminates when the subspace dimension reaches the preset upper limit or the residual norm falls within the tolerance range. The original high-dimensional linear equation system is projected onto the low-dimensional dense linear system corresponding to the upper Hessenberg matrix using the generated orthogonal basis matrix.

[0010] As a preferred embodiment of the hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver described in this invention, the variable quantum solution process adapted to a noisy medium-scale quantum processor based on a low-dimensional dense linear system encodes the coefficient matrix as a linear combination of Pauli operators and designs parameterized quantum circuits to prepare trial solution states. The specific steps are as follows: The coefficient matrix of a low-dimensional dense linear system is decomposed into a linear superposition of several Pauli strings, and each Pauli string corresponds to a set of tensor product Pauli operators that act on all qubits. Based on the number of available qubits, the dimension of the solution vector is mapped to the Hilbert space of the quantum state, ensuring that the solution vector can be completely represented by the amplitude of the quantum state; A hardware-efficient simulation structure was chosen as the parameterized quantum circuit, which consists of alternating stacks of single-qubit rotating gates and controlled X-gates between adjacent qubits. Set the initial input state of the circuit to a direct product state where all qubits are in the zero state; By adjusting the angle parameters of each rotating door in the circuit, different forms of trial solution states are generated to approximate the quantum state corresponding to the true solution.

[0011] As a preferred embodiment of the hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver described in this invention, the steps are as follows: First, a local loss function is defined according to the parameterized quantum circuit. Then, the gradient is calculated using the Hadamard test combined with the parameter translation rule. This drives the classical optimizer to iteratively update the circuit parameters until the loss function converges, and the quantum state corresponding to the optimal parameters is output. The local loss function is defined as the sum of normalized squared errors between multiple Pauli expectations and the projection of the target right-hand term onto the corresponding direction; For each adjustable parameter, a fixed phase offset is added and subtracted from its current value, respectively, to form two sets of perturbation parameters; For each set of perturbation parameters, the parameterized quantum circuit is run and the Hadamard test is performed to measure the expected value of the corresponding Pauli operator; The gradient component of the parameter is approximated by the quotient of the difference between the expected values ​​under two perturbations; The gradient components of all parameters are combined into a gradient vector, which is then input into a classic ADAM-type optimizer. The optimizer updates the parameter values ​​according to the gradient direction, repeats the above measurement and update process until the loss function value is lower than the preset convergence threshold, and outputs the final parameter configuration and the corresponding quantum state.

[0012] As a preferred embodiment of the hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver described in this invention, the steps of performing quantum state tomography on the quantum state to reconstruct a low-dimensional solution in classical vector form, and mapping the solution back to the original high-dimensional space using the aforementioned orthogonal basis vectors to obtain the physical quantity correction at the current time step are as follows: Determine a complete set of Pauli bases that cover all possible combinations of m-qubit Pauli operators; For the optimal quantum state of the output, perform multiple repeated measurements under each Pauli basis and statistically analyze the frequency distribution of the measurement results for each basis. Based on the measurement frequency, the density matrix is ​​reconstructed using the maximum likelihood estimation method; Extract the real parts of the main diagonal elements of the density matrix to construct the classical low-dimensional solution vector; Multiply the low-dimensional solution vector by the orthogonal basis matrix generated in claim 4 to obtain the solution correction in the original high-dimensional space; This correction is the increase in physical quantity at the current time step relative to the previous time step.

[0013] As a preferred embodiment of the hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver described in this invention, the specific steps of superimposing the physical quantity correction amount onto the solution of the previous time step to complete the implicit update of the current time step and realize high-precision advancement of the flow field variables are as follows: Read the current solution state from the unit basis function coefficients stored in the previous time step; The high maintenance positives obtained in claim 7 are aligned in the same unit and basis function order; The correction amount is added to the corresponding coefficient of the previous time step one by one to obtain the updated basis function coefficients of the current time step; The distribution of physical quantities within each unit is reconstructed using the updated coefficients; After all units have been updated, proceed to the next time step loop or terminate the computation process.

[0014] In a second aspect, the present invention provides a computer device including a memory and a processor, wherein the memory stores a computer program, wherein: when the computer program is executed by the processor, it implements any step of the hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver as described in the first aspect of the present invention.

[0015] Thirdly, the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein: when the computer program is executed by a processor, it implements any step of the hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver as described in the first aspect of the present invention.

[0016] The beneficial effects of this invention are as follows: By deeply integrating the high-order discontinuous finite element method with the variational quantum linear solver, a hybrid classical-quantum algorithm architecture for implicit solutions in computational fluid dynamics is constructed. The scheme utilizes the high-dimensional sparse linear system naturally derived from the discontinuous finite element method as the explicit physical input for quantum acceleration. After dimensionality reduction through classical Krylov subspace, it adapts to the capability boundary of current medium-scale quantum hardware with noise, avoiding the resource waste and error amplification caused by the black-box processing of arbitrary matrices by general quantum linear solvers. At the same time, through the coordinated optimization of local loss function design, parameterized circuit construction, and quantum state tomography reconstruction, it ensures that while retaining the high accuracy, strong shock wave capture capability, and unconditional stability of the implicit scheme of the discontinuous finite element method, it effectively breaks through the dimensional wall limitations of classical computing in terms of memory and computing power. Attached Figure Description

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

[0018] Figure 1 This is a flowchart of a hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver.

[0019] Figure 2 A schematic diagram of a circuit designed for efficient hardware.

[0020] Figure 3 This is a schematic diagram of the initial mesh for calculating transonic viscous flow around a 2D RAE2822 system.

[0021] Figure 4 The stress diagram for the convergence of the DDG(P1)+GMRES-IVQLS algorithm.

[0022] Figure 5 The graph shows the residual convergence curve of the DDG(P1)+GMRES-IVQLS algorithm. Detailed Implementation

[0023] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0024] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0025] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0026] Reference Figures 1-5 This is one embodiment of the present invention, which provides a hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver, including: S1. The fluid control equations are spatially discretized using the discontinuous finite element method. Higher-order basis functions are used to describe the solutions within the elements, and the information exchange between the element interfaces is processed by numerical flux to form a semi-discrete set of ordinary differential equations.

[0027] S2. An implicit time-progression scheme is adopted for the semi-discrete ordinary differential equation system, which transforms the problem of updating physical quantities at the current time step into the task of solving a high-dimensional sparse linear equation system.

[0028] S3. Apply the generalized minimum residual method to the high-dimensional sparse linear equation system on a classical computing platform, construct a set of orthogonal basis vectors of the Krylov subspace through the Arnoldi process, and project the original equation system onto the low-dimensional subspace to obtain a dense linear system with reduced dimension.

[0029] S4. Based on a low-dimensional dense linear system, construct a variational quantum solution process adapted to a noisy medium-scale quantum processor, encode the coefficient matrix as a linear combination of Pauli operators, and design parameterized quantum circuits to prepare trial solution states.

[0030] S5. Define a local loss function based on the parameterized quantum circuit, calculate the gradient using the Hadamard test combined with the parameter translation rule, drive the classical optimizer to iteratively update the circuit parameters until the loss function converges, and output the quantum state corresponding to the optimal parameters.

[0031] S6. Perform quantum state tomography on the quantum state to reconstruct a low-dimensional solution in classical vector form, and use the aforementioned orthogonal basis vectors to map the solution back to the original high-dimensional space to obtain the physical quantity correction at the current time step.

[0032] S7. The physical quantity correction is superimposed on the solution of the previous time step to complete the implicit update of the current time step, thereby achieving high-precision advancement of the flow field variables.

[0033] The governing equations for compressible fluids can be expressed as:

[0034] in Operators in dimensional space for ,and A vector of conserved variables. and These represent the fluid's density and total energy, respectively. The fluid is represented by a velocity vector along the coordinate directions, where in two-dimensional space... In three-dimensional space .

[0035] For the compressible Navier-Stokes equation, its inviscid flux Defined as:

[0036] in For Kronecker notation, The pressure of a fluid can be obtained from the following formula:

[0037] in This represents the general heat flux constant, which is 1.4 for an ideal gas; viscous flux vector. Defined as:

[0038] In two-dimensional space In three-dimensional space In the above equation, the component viscous stress tensor appears. and heat flux It is given by the following formula:

[0039] in The Prandtl number is for laminar flow, taken as 0.72 for sea-level air. Molecular viscosity can be calculated using Sutherland's law as follows:

[0040] in and Let the reference kinematic viscosity coefficient and reference temperature be respectively. It is a constant of 110.5K. Indicates the temperature of the fluid.

[0041] In actual calculations, the physical quantities to be solved are four: ρ, u, v, and p. Other physical quantities of interest, such as total energy E and temperature T, can be obtained directly by solving the equations of state.

[0042] Before introducing the discontinuous finite element discretization method for the Navier-Stokes equations, we will first give the basic element notation and overall framework of the discontinuous finite element method.

[0043] First, the computational domain is meshed, and the overall mesh is denoted as a set. For any one of the grid cells ,remember for directional boundary, The unit external normal vector is denoted as (where k is the spatial dimension); then, the following piecewise smooth finite element space is introduced as the numerical solution space:

[0044] Indicates that it is defined in a p-order polynomial space on.

[0045] The finite element space is provided Next, define the unit basis functions in this space. As a function space The basis; then, the variables to be solved. In each grid cell The above can be represented by a linear combination of basis functions:

[0046] in Representing the finite element space The dimension (or the number of basis units).

[0047] In the discontinuous finite element method, the choice of basis functions is crucial. Choosing appropriate basis functions can significantly reduce computational resources for solving linear systems. Furthermore, for specific computational problems, appropriate basis functions can better preserve and distinguish the physical structure of the numerical solution. Taylor basis functions are chosen as the discrete basis functions for general problems in the finite element method, where the solution variables... The second-order finite element expansion can be expressed as

[0048] The Taylor basis functions are:

[0049] in Indicates the number of basis functions. Representation unit middle The average value, , and Representing units respectively The coordinates of the maximum and minimum values; The dimension of the space. Representation unit The center of mass.

[0050] There are two main reasons for using Taylor's basis functions: First, Taylor basis functions are easier to implement in various shaped elements (tetrahedrons, hexahedrons, etc.) in three-dimensional space, with very low dependence on mesh elements, and can be easily implemented in almost any element in any two-dimensional or three-dimensional space. Second, using Taylor basis functions, it is easier to extend to various reconstructed numerical formats and implement three-dimensional weighted intrinsically non-oscillating (WENO) limiters, etc.

[0051] To obtain the discontinuous finite element discretization of the Navier-Stokes equations, an arbitrary test function is multiplied on both sides of the equation system. And in the unit Integrating by parts, we get

[0052] Then the finite element method was used to solve the space. elements in To approximate the solution U of the original system of equations, and simultaneously replacing the test space with a finite element space (the same finite element space as the solution space), we adopt the finite element space used here. The test function Φ is selected.

[0053] All units Add the above together, and the analytical solution Replace with test function Replace with basis functions that traverse the finite element space This yields the following discrete weak form:

[0054] Non-adhesive relative flow rate It is usually constructed using an approximate Riemannian solution; for example, the Harten-Lax-van Leer-Contact (HLLC) Flux, defined as follows:

[0055] in They respectively represent the results from equation (1) Internal and external calculations obtained The trajectory. For example, at the Gauss point. It can be done To calculate, where and They represent the first One degree of freedom and Internal basis functions; similarly, from Other elements to obtain .

[0056] Another widely used approximate Riemann solver is the Local Lax-Friedrichs (LLF) Flux, which is defined as:

[0057] in It is a matrix along the interface. The eigenvalue with the largest absolute value. express The flux Jacobi matrix.

[0058] Rusanov flux can also be used, which is defined as:

[0059] in .

[0060] Viscous numerical flux It can be constructed using either the DDG method or the BR2 method; here, the DDG method is used.

[0061] The DDG method can be viewed as a multi-period internal penalty function method. Compared to internal penalty function methods, the DDG method also considers the jumps in higher-order derivatives, making it easier to directly construct viscous flux using weak solutions; viscous flux It can be written in the following form:

[0062] Viscous numerical flux represented by the DDG method based on conserved variables of the Navier-Stokes equations. and It can be represented as follows:

[0063] and:

[0064] The directional components are represented as follows:

[0065] In the above formula, some symbols need to be introduced; here, the jump operator and the averaging operator are defined as follows:

[0066] in and These represent the left and right states of the interface, respectively. Values ​​obtained from both outside and inside the interface. and As a parameter of the penalty function, it is generally taken as and Choose here and As the default value; Indicates the length of the feature facet on a uniform grid. Generally taken as the size of the cell. When applied to unstructured meshes, Take It is the best, among them Representation unit From centroid to unit The displacement vector of the centroid, Representation surface The unit outward normal vector of the center point.

[0067] Based on the above formula, the viscous term and heat flux in the Navier-Stokes equations can be calculated using the DDG method based on conserved variables and original variables; for simplicity, an example of the DG(P1) approximation in three-dimensional space is given.

[0068] Step 1: Calculate the conservative variables at the unit interface. The first derivative in the direction; for example, the first derivative in the direction of the first derivative. The first derivative in the direction is as follows:

[0069]

[0070]

[0071]

[0072]

[0073] Step 2: Determine the original variables exist Calculate the heat flux density using the first derivative in the direction. ;For example, The first derivative in the direction is as follows:

[0074]

[0075]

[0076]

[0077] Step 3: Use the results obtained in Step 2 To calculate the viscous term and heat flux at each Gaussian orthogonal point on the unit interface defined in equation (2). Therefore, viscous numerical flux ,in .

[0078] For convenience, the weakly discrete form obtained above can be written as the following system of semi-discrete ordinary differential equations:

[0079] because It's about basis functions. The linear form, i.e., the approximate form of the discontinuous finite element method:

[0080] Substituting this form into a discrete system of ordinary differential equations, we obtain the following vector form of the ordinary differential equation system:

[0081] in It is the vector of degrees of freedom of discrete conserved variables. It is the residual vector. In the above formula... It is the mass matrix of the basis functions, and its child elements are Since the discontinuous finite element space is a piecewise smooth polynomial space, the inner product of the basis functions of different elements is 0, therefore the mass matrix... It is a highly sparse diagonal block matrix; in particular, if the basis functions on each element are... If it is a set of orthogonal functions, then the mass matrix... It will become a strictly diagonal matrix.

[0082] Right end item Let represent the discrete residual vector of the Navier-Stokes equations, assuming here that the number of equations is . The number of equations is actually the same as the number of physical equations; for example, the number of two-dimensional equations is... The number of three-dimensional equations is The residual vector can be set. Since the discontinuous finite element method is a piecewise smooth polynomial space, the inner product of the basis functions of different elements is... Therefore, the mass matrix of each element is independent of its neighboring elements, allowing for the analysis of the mass matrix form of each element independently. Simultaneously, for the residual vector in the vector equation... Or a conserved quantity For each component, the mass matrix form is also the same, and the mass matrix can be further subdivided into the form of each scalar equation. The semi-discrete form can be rewritten as follows:

[0083] in Represents cell index, The index of the equation. Represents the degree of freedom index. These represent the number of grid cells, the number of discrete equations (related to the dimension of the equations), and the number of degrees of freedom (same as the number of basis functions, related to the approximation accuracy); the sub-elements of the mass matrix are set as... For each Unit and each The mass submatrices of the equations have the same form.

[0084] It should be noted that the discontinuous finite element method used here is the discontinuous Galerkin method, and the test function form is the same as the basis function form. The above equation can be rewritten as follows:

[0085] here The vector of degrees of freedom for each unit and each equation. The residual vector corresponding to each equation for each unit. The mass submatrix for each unit and each equation; the mass submatrix is... A square matrix of order, where each element of the submatrix is ​​. .

[0086] It can be observed that the mass submatrix is ​​only related to the type and number of approximation basis functions, and the mass submatrix form is the same for each equation in each element; the same result is achieved for each equation in each grid, because the element values ​​are obtained through integration and are grid-dependent; without considering boundaries and virtual grid points, the equations mass matrix in The actual mass matrix of all meshes can be represented in the following form:

[0087] The matrix in the above formula is a diagonal block matrix, where each diagonal block is a unit. mass matrix , represented as:

[0088] Here It is also a diagonal block matrix, where the matrix on each diagonal is... It is necessary to further explain here that the submatrix elements The value of each equation is exactly the same, therefore it can be set as ; It can be rewritten as:

[0089] matrix order is , The order of the matrix is , The order is Furthermore, the matrix is ​​a large, highly sparse diagonal block matrix; given a fixed basis function, the matrix... It is invertible, and we can assume that the matrix is ​​a dense matrix. Of course, the sparsity of this submatrix block is related to the choice of basis functions. In actual calculations, it is necessary to calculate the product of a matrix and a vector, or the product of the inverse of a matrix and a vector. The reverse for:

[0090] in for

[0091] Large sparse matrix-vector multiplication can be divided into The dense submatrix is ​​the vector product of a dense matrix or its inverse. and At this point, the form of the matrix is ​​determined by the choice of basis functions and its dimension.

[0092] Submatrix elements The value of is determined by the form of the basis functions and the integral over each cell grid; in fact, Taylor basis functions of any order satisfy . This ensures the independence of the equations whose degrees of freedom are represented by the integral average, guaranteeing conservation regardless of how the other degrees of freedom are obtained through any other reconstruction method; at this point, the matrix... The form is:

[0093] For any ordinary basis function matrix It should be a dense matrix, for the mass submatrix composed of Taylor basis functions. The lower right submatrix block of a matrix is ​​usually a dense matrix. If this submatrix block is defined as... At this time, the matrix The inverse can also be expressed as:

[0094] In practical calculations, the mass matrix or its inverse is usually multiplied by something like a vector of degrees of freedom or a residual vector, assuming it is set as follows: ,in , The matrices introduced here are all square matrices. The order of the square matrix is ​​related to the dimension of the equation and the accuracy of the solution, as shown in the table below: ; To obtain the non-zero elements of the quality matrix, the element values ​​need to be calculated. Gaussian numerical integration was used to calculate the integral value. Because a polynomial space approximation was employed, a certain level of precision in Gaussian numerical integration could accurately recover the integral value. Depending on the required precision, for example, if a certain level of precision is required... To approximate the order with high precision, we need to use the highest-order polynomial of the basis functions. Therefore, at least the order needs to be recovered by Gaussian integral; A polynomial of order 1, that is, one that needs to be higher than 1. Gaussian integrals with algebraic precision; utilizing the integral mapping from the parameter domain to the physical domain, the integral is calculated as follows:

[0095] Next, we will present an implicit method for time discretization of this semi-discrete form, first transforming the spatially discretized form into a system of ordinary differential equations:

[0096] in Represents the mass matrix, For the solution vector, This is the residual vector.

[0097] The implicit time discretization scheme, namely the first-order backward Euler method, is used to solve the ordinary differential equations, accelerating the solution of the discrete system. The first-order Euler implicit time integral is used to linearize the spatial discrete equations in the time direction, specifically expressed as follows:

[0098] in It is a Jacobi matrix. It is the time step. It is the first in the direction of time. Step and The solution difference between time steps; the solution of this linear equation system at each time step is obtained using the GMRES method with the LU-SGS preprocessor.

[0099] Because DDG has a simple form for handling viscous fluxes, it can more easily obtain the exact Jacobian matrix of the numerical scheme of the Navier-Stokes equations. Using the exact Jacobian matrix in the implicit scheme can accelerate the convergence of the GMRES+LU-SGS algorithm, and in practice, it can approach the second-order convergence rate of Newton's method.

[0100] At the beginning of the time progression, strong transient phenomena will occur in the flow field due to the propagation of boundary information into the interior region.

[0101] Therefore, using smaller CFL coefficients (i.e. smaller time steps) in the initial stage allows the convection of the coefficient matrix in equation (3) to dominate, thus avoiding robustness issues.

[0102] To accelerate convergence, a local time step is used for each grid cell for steady flow:

[0103] In the above formula It is a unit The characteristic length, This is the maximum wave speed. As time progresses, the CFL coefficient increases using the following strategy:

[0104]

[0105]

[0106]

[0107] in , These are the upper and lower limits of the CFL growth factor; here we take... .

[0108] Finally, the specific algorithm steps for time advancement are given:

[0109] Next, the high-dimensional linear equations in step 4 are solved using the GMRES-IVQLS method. First, the high-dimensional linear equations are transformed into a series of low-dimensional linear equations using the GMRES method, and then the IVQLS method is used to efficiently solve the reduced-dimensional small-scale linear equations.

[0110] GMRES-IVQLS quantum computing method: set up , ,say

[0111] It is by and The generated Krylov subspace is usually abbreviated as or The GMRES method searches within the Krylov subspace. To satisfy Thus obtained That is, in the subspace The best approximation of the solution vector.

[0112] The GMRES method flow is as follows:

[0113] The GMRES method (Generalized Minimum Residual Method) transforms the high-dimensional linear equation system into searching for an approximate solution in the Krylov subspace, thus reducing the problem to solving a low-dimensional linear equation system in 27 steps. Then, the low-dimensional linear equations are solved using the IVQLS method.

[0114] Iterative Variational Quantum Linear Solver (IVQLS) method: For a given system of linear algebraic equations ,in Assuming (If not satisfied, then use) Element filling matrix with vector (to ensure that its dimension meets the requirements), the steps of the IVQLS method are as follows: ①The matrix It can be decomposed into a linear combination of a series of unitary matrices; Right now ,in , It is a unitary matrix; the method used here is to... Decomposed into a linear combination of Pauli strings, where the Pauli matrix is... The unitary matrix is ​​defined as follows:

[0115] Then matrix It can be broken down into: ,in , That is The quaternary representation of the coefficient, and the coefficient (in The trace of a square matrix is ​​the sum of the elements along its diagonal. ② Constructing and preparing solution vectors Variable quantum circuit (in ); Variable quantum circuits Used to prepare the solution vector, i.e. ; Specific variable quantum circuits There are many types; this study uses a hardware efficient ansatz circuit for construction, as shown in the attached figure. Figure 1 As shown.

[0116] Among them, the single-bit quantum gate is

[0117] It is a rotating gate operating on a single quantum bit, a series of A rotating gate is called a rotating layer; a two-qubit quantum gate is... Door, with The symbol represents the controlled qubit, and the other represents the control qubit. When the control qubit is... When in a state, the action on the controlled qubit Gate, when the control bit is In this state, no operation is performed, a series of The gate is called an entanglement layer; theoretically, by combining a sufficient number of rotational layers with entanglement layers, it is possible to express... Any unit vector in space; ③ Define the loss function to be optimized (Cost Function) The definition of loss functions mainly includes the global loss function. and local loss function Two kinds; and Both describe vectors with vector The degree of closeness, when and When they are close to being in the same direction, both loss functions tend to be .

[0118] The unnormalized global loss function is defined as:

[0119] in , , ,right Normalization yields:

[0120] in:

[0121] and:

[0122]

[0123] The unnormalized local loss function is defined as:

[0124] in:

[0125] then:

[0126]

[0127] right Normalization yields:

[0128]

[0129] Because the global loss function is prone to the barren plateau phenomenon (i.e., the gradients of each component quickly vanish to 0), a local loss function is used here to optimize the variable quantum line. For loss functions of the form ... For the elements, this study uses the Hadamard Test method for efficient computation, employing an additional auxiliary qubit. The Hadamard Test method is shown below: Step 1: Constructing a quantum circuit (in For the Hadamard door, Controlled (gate), which is applied to the quantum state Above:

[0130] Step 2: Measuring the auxiliary qubit The probability of this state occurring can be calculated as follows:

[0131] Therefore:

[0132] Similarly, replacing the quantum circuit with We can obtain:

[0133] and: Efficient approximation can be achieved through multiple Hadamard tests. ; ④ Optimization methods and gradient calculation; In this study, the ADAM optimization algorithm was selected to optimize the variable quantum circuit.

[0134] The first one to be optimized The parameters are The update formula for the Adam algorithm is as follows:

[0135] in, For parameters The gradient; for The speed and the initial value is 0; for The cumulative sum of squares and the initial value is 0; for The attenuation coefficient has a value range of [0,1], and is generally set to 0.9; for The attenuation coefficient has a value range of [0,1], and is generally set to 0.999; This represents the current iteration number; The learning rate; It is a very small constant, and its purpose is to prevent the denominator from being 0.

[0136] Use the parameter shift rule to accurately calculate the gradient: for ,in , , Calculations show that:

[0137] Pick , That is:

[0138] Compared to the traditional perturbation method for approximating gradient calculation, the parameter translation method has no error; ⑤ Settings Tolerance Iterative solution of linear equations ; Step 1: Settings ; Step 2: Calculation ,like Then set And call the ADAM optimization algorithm to optimize the variable quantum circuit. loss function Optimize; otherwise proceed to the next step. At this point, we need to solve the system of linear equations. When the following situations occur:

[0139] Time settings ,in: To achieve minimization:

[0140] To obtain (among them) ).

[0141] Return to step 1; Step 3: Obtain the final approximate solution vector .

[0142] In this study, regarding the quantum state in step 2... The reading adopted Quantum state tomography: Input: Tolerance From quantum state Preparation of quantum states quantum circuits Controlled quantum circuits ; Output Classical vectors ,satisfy and ; 1: Measurement under standard calculation basis sub-quantum state And record the quantum state. Number of times and order ; 2: For quantum states ,effect The quantum state was then obtained. Measurement under standard calculation basis sub-quantum state And record the quantum state. Number of times ; 3: When Time settings Otherwise set ; 4: Obtain the classical vector ,in .

[0143] The classic two-dimensional RAE2822 transonic viscous flow problem is used as a case study for testing. The initial mesh is shown in the attached figure. Figure 2 As shown; the test parameter is Mach number. Angle of attack Degree, Reynolds number Inviscid fluxes use Rusanov fluxes, viscous fluxes use DDG fluxes, and the mesh is an unstructured mesh, including... The example uses 32 CPU cores for parallel computation of four quadrilateral units. DDG(P1) and DDG(P2) take [time values ​​missing]. minutes and Minutes. The global pressure distribution of the final convergence result is shown in the attached figure. Figure 3 As shown in the attached figure, the residual convergence curve is as follows. Figure 4 As shown.

[0144] This embodiment also provides a computer device suitable for hybrid algorithms based on classical discontinuous finite element method and variable quantum linear solver, including: a memory and a processor; the memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions to implement the hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver as proposed in the above embodiment.

[0145] The computer device can be a terminal, comprising a processor, memory, communication interface, display screen, and input devices connected via a system bus. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, carrier networks, NFC (Near Field Communication), or other technologies. The display screen can be an LCD screen or an e-ink screen. The input devices can be a touch layer covering the display screen, buttons, a trackball, or a touchpad on the computer device's casing, or an external keyboard, touchpad, or mouse.

[0146] This embodiment also provides a storage medium on which a computer program is stored. When executed by a processor, the program implements the hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver as proposed in the above embodiment. The storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read Only Memory (EPROM), Programmable Red-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk.

[0147] In summary, this invention deeply integrates the high-order discontinuous finite element method with a variational quantum linear solver, constructing a hybrid classical-quantum algorithm architecture for implicit solutions in computational fluid dynamics. The scheme utilizes the high-dimensional sparse linear system naturally derived from the discontinuous finite element method as an explicit physical input for quantum acceleration. After dimensionality reduction using classical Krylov subspaces, it adapts to the capability boundaries of current medium-scale quantum hardware with noise, avoiding the resource waste and error amplification caused by the black-box processing of arbitrary matrices by general-purpose quantum linear solvers. At the same time, through the coordinated optimization of local loss function design, parameterized circuit construction, and quantum state tomography reconstruction, it ensures that while retaining the high accuracy, strong shock wave capture capability, and unconditional stability of the implicit scheme of the discontinuous finite element method, it effectively overcomes the dimensionality wall limitations of classical computing in terms of memory and computing power.

[0148] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver, characterized by: include: The fluid control equations are spatially discretized using the discontinuous finite element method. Higher-order basis functions are used to describe the solutions within the elements, and numerical flux is used to process the information exchange at the element interfaces, forming a semi-discrete set of ordinary differential equations. An implicit time-progression scheme is adopted for the semi-discrete ordinary differential equation system, which transforms the problem of updating physical quantities at the current time step into the task of solving a high-dimensional sparse linear equation system. The generalized minimum residual method is applied to the high-dimensional sparse linear equation system on a classical computing platform. A set of orthogonal basis vectors of the Krylov subspace is constructed through the Arnoldi process, and the original equation system is projected onto this low-dimensional subspace to obtain a dense linear system with reduced dimension. Based on low-dimensional dense linear systems, a variational quantum solution process adapted to noisy medium-scale quantum processors is constructed. The coefficient matrix is ​​encoded as a linear combination of Pauli operators, and parameterized quantum circuits are designed to prepare trial solution states. Based on the local loss function defined by the parameterized quantum circuit, the gradient is calculated by combining the Hadamard test with the parameter translation rule, which drives the classical optimizer to iteratively update the circuit parameters until the loss function converges, and outputs the quantum state corresponding to the optimal parameters. Quantum state tomography is performed on the quantum state to reconstruct a low-dimensional solution in classical vector form. The solution is then mapped back to the original high-dimensional space using the aforementioned orthogonal basis vectors to obtain the physical quantity correction at the current time step. The physical quantity correction is superimposed on the solution of the previous time step to complete the implicit update of the current time step, thereby achieving high-precision advancement of the flow field variables.

2. The hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver as described in claim 1, characterized in that: The fluid control equations are spatially discretized using the discontinuous finite element method, and the internal solutions of the elements are described using higher-order basis functions. Information exchange at the element interfaces is processed through numerical flux, forming a semi-discrete system of ordinary differential equations. The specific steps are as follows: The computational region is divided into several non-overlapping computational units; Within each computational unit, a set of higher-order polynomial basis functions is introduced. These basis functions use the average value of physical quantities at the unit center as the zero-order term, and the remaining higher-order terms are constructed from the offset of local coordinates relative to the unit center. The conserved variables in the fluid control equations are expanded in each cell according to the higher-order basis functions and substituted into the weak-form integral expression; Perform integration by parts on the spatial derivative terms involved in the weak form, transfer the derivatives in the volume integral to the basis functions, and at the same time explicitly separate the flux integral terms on the unit boundary; At the interface shared by adjacent cells, the continuity of solutions on both sides is not forced. Instead, the numerical flux is calculated using the Roe or Lax-Friedrichs scheme, which depends only on the reconstructed values ​​of the conserved variables provided by the cells on both sides of the interface. After completing the weak-form assembly of all elements, a system of semi-discrete ordinary differential equations is obtained, consisting of the time derivatives of the basis function coefficients of each element.

3. The hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver as described in claim 2, characterized in that: The semi-discrete ordinary differential equation system adopts an implicit time-progression scheme, transforming the physical quantity update problem at the current time step into the task of solving a high-dimensional sparse linear equation system. The specific steps are as follows: Select an implicit time discretization scheme, including but not limited to the first-order backward Euler scheme or the second-order implicit Runge-Kutta scheme; Integrate the semi-discrete ordinary differential equation system in the current time step, so that the time derivative term is expressed as the difference between the solution in the current time step and the solution in the previous time step divided by the time step size. The spatial residual term is linearized at the current time step, and the first derivative term with respect to the current solution is retained to form the Jacobian matrix. The coefficient matrix of the linear system is formed by multiplying the mass matrix by the inverse of the time step and then superimposing the results into the Jacobian matrix. The result of interacting the residual term of the previous time step with the mass matrix is ​​used as the right-hand side term of the linear system. This leads to the construction of a high-dimensional sparse linear system of equations with the unknown solution at the current time step as the variable.

4. The hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver as described in claim 3, characterized in that: The method described above applies the generalized minimum residual method to high-dimensional sparse linear equations on a classical computing platform. It constructs a set of orthogonal basis vectors for the Krylov subspace through the Arnoldi process and projects the original equations onto this low-dimensional subspace to obtain a denser linear system with reduced dimensionality. The specific steps are as follows: The initial vector for initializing the Krylov subspace is the direction after the right-hand side is normalized; Perform Arnoldi iterations sequentially: In each step, apply the coefficient matrix to the current basis vector, and then subtract its projection components on the existing basis vectors from the result vector to achieve Gram-Schmidt orthogonalization. Normalize the orthogonalized vectors and add them to the subspace as the next basis vectors; The projection coefficients generated during each orthogonalization process are recorded synchronously to form the upper Hessenberg matrix; The Arnoldi process terminates when the subspace dimension reaches the preset upper limit or the residual norm falls within the tolerance range. The original high-dimensional linear equation system is projected onto the low-dimensional dense linear system corresponding to the upper Hessenberg matrix using the generated orthogonal basis matrix.

5. The hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver as described in claim 4, characterized in that: The aforementioned variational quantum solution process, based on a low-dimensional dense linear system and adapted to a noisy medium-scale quantum processor, encodes the coefficient matrix as a linear combination of Pauli operators and designs parameterized quantum circuits to prepare trial solution states. The specific steps are as follows: The coefficient matrix of a low-dimensional dense linear system is decomposed into a linear superposition of several Pauli strings, and each Pauli string corresponds to a set of tensor product Pauli operators that act on all qubits. Based on the number of available qubits, the dimension of the solution vector is mapped to the Hilbert space of the quantum state, ensuring that the solution vector can be completely represented by the amplitude of the quantum state; A hardware-efficient simulation structure was chosen as the parameterized quantum circuit, which consists of alternating stacks of single-qubit rotating gates and controlled X-gates between adjacent qubits. Set the initial input state of the circuit to a direct product state where all qubits are in the zero state; By adjusting the angle parameters of each rotating door in the circuit, different forms of trial solution states are generated to approximate the quantum state corresponding to the true solution.

6. The hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver as described in claim 5, characterized in that: The process involves defining a local loss function based on a parameterized quantum circuit, calculating the gradient using the Hadamard test combined with the parameter translation rule, driving a classical optimizer to iteratively update the circuit parameters until the loss function converges, and outputting the quantum state corresponding to the optimal parameters. The specific steps are as follows: The local loss function is defined as the sum of normalized squared errors between multiple Pauli expectations and the projection of the target right-hand term onto the corresponding direction; For each adjustable parameter, a fixed phase offset is added and subtracted from its current value, respectively, to form two sets of perturbation parameters; For each set of perturbation parameters, the parameterized quantum circuit is run and the Hadamard test is performed to measure the expected value of the corresponding Pauli operator; The gradient component of the parameter is approximated by the quotient of the difference between the expected values ​​under two perturbations; The gradient components of all parameters are combined into a gradient vector, which is then input into a classic ADAM-type optimizer. The optimizer updates the parameter values ​​according to the gradient direction, repeats the above measurement and update process until the loss function value is lower than the preset convergence threshold, and outputs the final parameter configuration and the corresponding quantum state.

7. The hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver as described in claim 6, characterized in that: The process of performing quantum state tomography on the quantum state to reconstruct a low-dimensional solution in classical vector form, and then mapping this solution back to the original high-dimensional space using the aforementioned orthogonal basis vectors to obtain the physical quantity correction at the current time step, specifically involves the following steps: Determine a complete set of Pauli bases that cover all possible combinations of m-qubit Pauli operators; For the optimal quantum state of the output, perform multiple repeated measurements under each Pauli basis and statistically analyze the frequency distribution of the measurement results for each basis. Based on the measurement frequency, the density matrix is ​​reconstructed using the maximum likelihood estimation method; Extract the real parts of the main diagonal elements of the density matrix to construct the classical low-dimensional solution vector; Multiply the low-dimensional solution vector by the orthogonal basis matrix generated in claim 4 to obtain the solution correction in the original high-dimensional space; This correction is the increase in physical quantity at the current time step relative to the previous time step.

8. The hybrid algorithm based on classical discontinuous finite element method and variational quantum linear solver as described in claim 7, characterized in that: The process of adding the physical quantity correction to the solution of the previous time step to complete the implicit update of the current time step and achieve high-precision advancement of the flow field variables is as follows: Read the current solution state from the unit basis function coefficients stored in the previous time step; Align the high maintenance positives obtained in claim 7 with the same unit and basis function order; The correction amount is added to the corresponding coefficient of the previous time step one by one to obtain the updated basis function coefficients of the current time step; The distribution of physical quantities within each unit is reconstructed using the updated coefficients; After all units have been updated, proceed to the next time step loop or terminate the computation process.

9. A computer device, comprising a memory and a processor, wherein the memory stores a computer program, characterized in that: When the processor executes the computer program, it implements the steps of the hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver as described in any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by the processor, it implements the steps of the hybrid algorithm based on classical discontinuous finite element method and variable quantum linear solver as described in any one of claims 1 to 8.

Images