Efficient differentiable fluid optimization method and system based on flow mapping and proxy gradients

Abstract

The application provides an efficient differentiable fluid optimization method and system based on flow mapping and proxy gradient, the method comprising collecting boundary conditions, current control variables and objective functions of a fluid scene to be optimized; adopting a flow mapping method to construct a forward simulation calculation graph, dividing the overall simulation time into multiple renormalization periods, each renormalization period sequentially completing three stages of velocity propulsion, flow mapping evolution and renormalization, obtaining a terminal velocity field through multi-period iterative propulsion simulation, combining the objective function to complete quantitative evaluation of the simulation process, and calculating the objective function value; relying on the forward simulation calculation graph to build a proxy gradient calculation graph, distinguishing a first adjoint path composed of a pullback operator from a second adjoint path of gradient propagation of flow mapping intermediate variables, and obtaining a proxy gradient; and iteratively updating the current control variables based on the proxy gradient until the objective function converges. The application effectively reduces the calculation cost and memory consumption of the differentiable fluid simulation optimization process.

Description

Technical Field

[0001] This invention relates to the fields of computer graphics, physical simulation, differentiable programming and fluid optimization technology, and specifically to an efficient differentiable fluid optimization method and system based on flow mapping and surrogate gradient. Background Technology

[0002] Fluid optimization technology is widely used in many technical fields such as aerodynamic shape design, robot fluid control, smoke and flame animation control, fluid parameter identification, biological fluid modeling, and fluid phenomenon data generation. These applications all require solving complex inverse fluid problems, typically including simulating the initial velocity field based on the target vortex structure, controlling the evolution of smoke morphology through keyframe parameters, and optimizing fluid lift-to-drag ratio, pressure distribution, and wake structure by adjusting object boundary control parameters.

[0003] Differentiable fluid simulation technology reconstructs various computational operators in the fluid simulation process into differentiable operators, which can provide end-to-end analytical gradients for iterative optimization algorithms. Compared with traditional gradient-free random search methods, it can provide accurate and effective parameter search directions, significantly improve the iterative convergence efficiency of high-dimensional control parameter spaces, and is the core technical means for optimizing inverse fluid problems.

[0004] However, high-fidelity fluid simulation systems include multiple high-precision computation operators such as velocity semi-Lagrange advection, pressure projection, high-order time integration, and complex boundary handling. Directly performing backpropagation on the complete forward simulation computation graph will generate a deep and computationally dense backpropagation graph. At the same time, a large amount of cache storage is required for intermediate state data throughout the simulation process. Ultimately, this leads to a significant increase in the running time of simulation optimization and overload of device memory, making it difficult to adapt to long-term, large-scale fluid optimization scenarios.

[0005] Existing methods for continuous adjoint partial differential equations can slightly reduce the storage pressure of intermediate states by optimizing the gradient calculation logic. However, the backpropagation path of this method has an inherent logical deviation from the discrete operator of the forward simulation. In complex scenarios with long-term simulations and strong vortex dynamics coupling, significant gradient shifts and distortions occur, greatly reducing the accuracy of fluid parameter optimization. Conventional flow mapping methods can effectively preserve vortex structure characteristics and improve simulation fidelity by pulling back the structure across time steps during forward simulation. However, the complete backpropagation process of this method still requires the transmission and calculation of flow mapping parameters, Jacobian matrix, and velocity gradients step by step, resulting in high redundancy in backpropagation calculations, long iteration time, and large memory consumption.

[0006] A patent search revealed an invention patent with publication number CN114021497A, which discloses a compressible turbulent fluid topology optimization method based on automatic differentiation. The method involves establishing a geometric model for the corresponding flow channel topology optimization; obtaining the basic parameters of the topology optimization object and constructing a CFD topology optimization model; constructing an objective function based on the basic parameters, forming a nonlinear programming problem; deriving the sensitivity equation based on the adjoint method; solving the fluid control equations to obtain the flow field results and the value of the objective function, and outputting the new flow field results; using automatic differentiation technology to construct the Jacobian matrix and gradient vector; constructing the adjoint equation in matrix form to obtain the adjoint multiplier; substituting the new flow field results and the adjoint multiplier into the sensitivity equation to solve for the sensitivity; using the MMA numerical optimization method, combined with the nonlinear programming problem to build a mathematical model; optimizing the flow channel, updating the design variables, obtaining the optimal solution, and outputting the optimal two-dimensional topology configuration. However, this patent only uses automatic differentiation and the adjoint method to calculate the gradient, resulting in a large computational load and high memory consumption for backpropagation, making it difficult to adapt to complex fluid optimization scenarios with long durations and high vortex dynamics.

[0007] In summary, given the problems of the existing technologies, researching an efficient differentiable fluid optimization method and system based on flow mapping and surrogate gradients has become a critical task that urgently needs to be addressed. Summary of the Invention

[0008] To address the shortcomings of existing technologies, the purpose of this invention is to provide an efficient differentiable fluid optimization method and system based on flow mapping and surrogate gradients.

[0009] An efficient differentiable fluid optimization method based on flow mapping and surrogate gradients, provided by the present invention, includes the following steps: Step S1: Collect the boundary conditions and current control variables of the fluid scenario to be optimized through physical sensors, and obtain the objective function through the user input interface; Step S2: Based on the boundary conditions and the current control variables, construct the forward simulation calculation graph using the flow mapping method, execute forward propagation covering the simulation time, obtain the endpoint velocity field, quantitatively evaluate the endpoint velocity field and the intermediate states saved during the forward propagation process based on the endpoint velocity difference, and calculate the objective function value under the current control variables. Step S3: Based on the forward simulation computation graph, build and execute the proxy gradient computation graph. With the objective function value as the optimization target for backpropagation, take the gradient of the objective function with respect to the final velocity field as the starting point of backpropagation, and propagate forward layer by layer along the proxy gradient computation graph. After the proxy gradient is constructed, the proxy gradient at the current control variable is obtained. Step S4: Based on the surrogate gradient at the current control variable, update the current control variable, and use the updated current control variable as the input for the next iteration. Repeat steps S2 to S4 until the objective function converges, and use the current control variable corresponding to the convergence as the optimization result. Step S5: Convert the optimization results into executable velocity field data and output the executable velocity field data to the fluid simulation display device to control the fluid simulation display device to display the initial velocity field that can evolve into the target vortex structure.

[0010] Preferably, in step S1, the fluid scenario to be optimized is a velocity inversion scenario, the physical sensors are a flow velocity sensor and a pressure sensor, the current control variable is the initial velocity field, and the objective function is the final velocity difference.

[0011] Preferably, in step S2, the simulation time is divided into K sequentially connected renormalization cycles. Each renormalization cycle contains N consecutive time steps. Each renormalization cycle corresponds to a cycle start velocity field and a cycle end velocity field. Within each renormalization cycle, the velocity propagation stage, the flow mapping evolution stage, and the renormalization stage are passed sequentially. The initial velocity field is used as the cycle start velocity field of the first renormalization cycle, and the cycle end velocity field of the previous renormalization cycle is used as the cycle start velocity field of the next renormalization cycle. The velocity field is updated cycle by cycle until the cycle end velocity field of the last renormalization cycle is reached, which is used as the end velocity field.

[0012] Preferably, step S2 includes the following steps: Step S2.1: Based on the initial velocity field and boundary conditions, establish the global variables required for fluid simulation. The global variables include the velocity field stored as an interlaced MAC grid, the pressure field stored as a regular grid, the fluid volume fraction, the forward flow mapping, the reverse flow mapping, the forward flow mapping Jacobian matrix, and the reverse flow mapping Jacobian matrix. Step S2.2: If it is the first renormalization cycle, the initial velocity field is used as the starting velocity field of the current renormalization cycle; otherwise, the ending velocity field of the previous renormalization cycle is used as the starting velocity field of the current renormalization cycle. At the same time, the pressure field is set to zero, the fluid volume fraction is set according to the boundary conditions, and the forward flow mapping and reverse flow mapping at the starting moment of the current renormalization cycle are initialized to the identity mapping of the grid position. The corresponding forward flow mapping Jacobian matrix and reverse flow mapping Jacobian matrix are initialized to the identity matrix. Step S2.3: In the velocity propagation phase, perform semi-Lagrangian velocity advection and pressure projection on the velocity field at the starting point of the current renormalization cycle, generate the midpoint velocity sequence step by step, and use the pressure projection to constrain the velocity field to the divergence-free velocity space and obtain the pressure field. Step S2.4: In the flow mapping evolution stage, the forward flow mapping, reverse flow mapping, forward flow mapping Jacobian matrix and reverse flow mapping Jacobian matrix at the starting moment of the current renormalization cycle are driven by the midpoint velocity sequence to record the positional correspondence and local deformation of the fluid particles from the starting moment of the current renormalization cycle to the ending moment of the cycle and from the ending moment of the cycle back to the starting moment of the cycle. Step S2.5: In the renormalization stage, the velocity field at the beginning of the current renormalization cycle is interpolated and sampled using the reverse flow mapping at the end of the current renormalization cycle, and the sampled velocity vector is stretched using the corresponding reverse flow mapping Jacobian matrix to obtain the velocity field before projection at the end of the current renormalization cycle. Step S2.6: Perform pressure projection on the velocity field before projection to obtain the velocity field at the end of the current renormalization cycle; Step S2.7: Determine whether the current renormalization cycle is the last cycle. If not, store the end velocity field of the current renormalization cycle in the buffer as the starting velocity field of the next cycle, and return to step S2.2 to continue execution. If yes, terminate the loop, use the end velocity field of the current renormalization cycle as the end velocity field of the pressure field sequence of the entire simulation, and quantitatively evaluate the end velocity field, the midpoint velocity sequence, forward flow mapping, reverse flow mapping, forward flow mapping Jacobian matrix, and reverse flow mapping Jacobian matrix saved during the simulation based on the end velocity difference, and calculate the objective function value under the current control variable.

[0013] Preferably, in step S3, based on the gradients of each computation path in the forward simulation computation graph, a gradient propagation path in the proxy gradient computation graph is constructed. The gradient propagation path is divided into a first adjoint path and a second adjoint path. The first adjoint path is a step propagation path, which forms a step connection by a pullback operator. It is used to propagate the gradient of the velocity field at the end of the period to the velocity field at the beginning of the period through the adjoint pullback operator. The second adjoint path is a path that propagates step by step through forward flow mapping, reverse flow mapping, forward flow mapping Jacobian matrix, reverse flow mapping Jacobian matrix and midpoint velocity sequence. It is used to provide auxiliary gradient correction.

[0014] Preferably, in step S3, the second adjoint path includes Eulerian-Eulerian gradient propagation, Lagrange-Lagrange gradient propagation, and Lagrange-Eulerian gradient propagation; Euler-Euler gradient propagation is used to propagate gradients between Euler velocity grids, preserving time-step propagation; Lagrange-Lagrange gradient propagation is used to propagate gradients between flow map variables and Jacobian variables. It replaces time-step propagation with long-step propagation that spans a renormalization period. Long-step propagation propagates the flow map gradient and flow map Jacobian gradient with a time step size corresponding to the period length of the renormalization period. Lagrange-Euler gradient propagation is used to couple the gradients of the flow map variables and Jacobian variables back to the Euler velocity grid. The coupling is performed at the midpoint of the renormalization period, so that the gradients of the flow map variables and Jacobian variables are propagated to the midpoint velocity gradient through the midpoint Lagrange-Euler coupling, and then propagated to the velocity field at the beginning of the renormalization period through stepwise Euler-Euler propagation.

[0015] Preferably, step S3 includes the following steps: Step S3.1: Take the last renormalization cycle as the current renormalization cycle, and at the end velocity field of the current renormalization cycle, calculate the end velocity gradient of the end velocity difference as the starting point of backpropagation. Step S3.2: Propagate the endpoint velocity gradient to the pre-projection velocity gradient through the adjoint form of the pressure projection operator, so that the reverse velocity gradient satisfies the incompressible constraint, and obtain the pre-projection velocity gradient at the end of the current renormalization cycle. Step S3.3, execute the first adjoint path: propagate the velocity gradient before projection to the velocity field at the starting point of the current renormalization period through the adjoint pullback operator to obtain the dominant gradient of the current renormalization period; Step S3.4: Propagate the dominant gradient to the backflow mapping gradient and backflow mapping Jacobian gradient at the end of the current renormalization cycle via the adjoint pullback operator, and use them as the initial gradient of the second adjoint path. Step S3.5: Under the condition of satisfying the Courant-Friedrich-Liouvi constraints, the gradient of the backflow mapping and the gradient of the backflow mapping Jacobian matrix are back-propagated from the end time of the current renormalization cycle to the start time along the Lagrange trajectory. The N time-step propagation is replaced by a single long-step propagation that spans N consecutive time steps, so as to obtain the Lagrange variable gradient after the cross-cycle propagation of the current renormalization cycle.

[0016] Step S3.6: At the midpoint of the current renormalization cycle, couple the Lagrange variable gradient back to the Eulerian velocity grid to form the midpoint velocity gradient of the current renormalization cycle. Step S3.7: Propagate the midpoint velocity gradient stepwise to the velocity field at the starting point of the current renormalization cycle in reverse order of the time steps to obtain the velocity gradient at the starting point of the current renormalization cycle, which serves as the auxiliary gradient correction for the second adjoint path. Step S3.8: Accumulate the corrections of the dominant gradient and the auxiliary gradient of the second adjoint path to obtain the surrogate gradient at the starting point of the current renormalization cycle. Step S3.9: Determine whether the current renormalization period is the first renormalization period. If yes, terminate the loop and the surrogate gradient at the start of the current renormalization period is the surrogate gradient of the initial velocity field. If no, assign the surrogate gradient at the start of the current renormalization period to the gradient of the velocity field at the end of the previous renormalization period, update the current renormalization period to the previous renormalization period, and return to step S3.2 to continue execution.

[0017] Preferably, in step S4, the Adam optimizer or gradient descent optimizer is used to update the current control variable.

[0018] Preferably, in step S1, the fluid scene to be optimized also includes a smoke control scene, the current control variable is the keyframe velocity field, the objective function is the smoke density difference, and the physical sensors include a smoke acquisition device and an image acquisition device.

[0019] This invention also provides an efficient differentiable fluid optimization system based on flow mapping and surrogate gradients, employing the aforementioned efficient differentiable fluid optimization method based on flow mapping and surrogate gradients, comprising: The data acquisition module collects the boundary conditions and current control variables of the fluid scenario to be optimized through physical sensors, and obtains the objective function through the user input interface; The forward simulation module constructs a forward simulation computation graph based on boundary conditions and current control variables using the flow mapping method, executes forward propagation covering the simulation time, obtains the final velocity field, quantitatively evaluates the final velocity field and the intermediate states saved during the forward propagation process based on the difference in the final velocity, and calculates the objective function value under the current control variables. The proxy gradient module, based on the forward simulation computation graph, builds and executes the proxy gradient computation graph. The objective function value is used as the optimization target for backpropagation. The gradient of the objective function with respect to the final velocity field is used as the starting point for backpropagation. It propagates forward layer by layer along the proxy gradient computation graph. After the proxy gradient is constructed, the proxy gradient at the current control variable is obtained. The parameter update module updates the current control variable based on the surrogate gradient at the current control variable, and uses the updated current control variable as the input for the next iteration. Steps S2 to S4 are repeated until the objective function converges, and the current control variable at the convergence point is used as the optimization result. The output module converts the optimization results into executable velocity field data and outputs the executable velocity field data to the fluid simulation display device to control the fluid simulation display device to display the initial velocity field that can evolve into the target vortex structure.

[0020] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention identifies the dominant gradient path in a differentiable simulation system of a flow map, distinguishing between the dominant gradient of the first adjoint path and the auxiliary gradient correction term of the second adjoint path, thus abandoning the existing method of blindly differentiating the complete backpropagation graph node by node. Simultaneously, it combines long-step Lagrange-Lagrange gradient propagation, replacing the traditional multiple gradient transfer operations per time step with a single long-step gradient propagation across the renormalization period, significantly reducing the number of backpropagation iterations of flow map parameters and the Jacobian matrix, effectively reducing the overall computational load of backpropagation and improving iterative optimization efficiency.

[0021] 2. This invention constructs three types of gradient propagation coupling mechanisms: stepwise Euler-Euler, long-step Lagrange-Lagrange, and midpoint Lagrange-Euler. Relying on the physical duality of different gradient propagation modes, it fully preserves core physical gradient information such as fluid velocity advection, vortex evolution, and pressure constraints. It corrects the gradient deviation problem caused by the mismatch between gradient propagation and forward simulation operators in traditional adjoint methods. In long-term, strong vortex dynamics scenarios, it can still maintain the optimization accuracy close to that of a complete gradient calculation scheme.

[0022] 3. By reconstructing the simulation and gradient propagation logic, this invention eliminates the need to cache all intermediate state data of the complete simulation time step. It only needs to retain core parameters such as velocity field, flow mapping, and Jacobian matrix of key nodes in the renormalization cycle, which greatly reduces the amount of data in the forward simulation state cache and the backward auxiliary gradient cache. This effectively reduces the memory overhead in the differentiable fluid optimization process and improves the algorithm's adaptability to conventional hardware devices.

[0023] 4. The optimized method of this invention can be adapted to various fluid optimization tasks such as velocity inversion, smoke morphology control, and aerodynamic shape optimization. It solves the problem that existing technologies cannot adapt to long-term, high-vortex complex fluid scenarios. While ensuring optimization accuracy and computational efficiency, it greatly improves the engineering practical value of the microfluidic optimization scheme. Attached Figure Description

[0024] Other features, objects, and advantages of the present invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings: Figure 1 This is a schematic diagram of a highly efficient differentiable fluid optimization system based on flow mapping and surrogate gradient in an embodiment of the present invention; Figure 2 The following is a flowchart of an efficient differentiable fluid optimization method based on flow mapping and surrogate gradient in an embodiment of the present invention; Figure 3 This is a schematic diagram of the forward propagation calculation in an embodiment of the present invention; Figure 4 This is a schematic diagram of stream mapping transport in an embodiment of the present invention; Figure 5 This is a schematic diagram comparing the complete gradient and the proxy gradient calculation graph in an embodiment of the present invention; Figure 6 This is a schematic diagram of the proxy gradient structure of long-step L2L and step-by-step E2E in an embodiment of the present invention; Figure 7 This is a schematic diagram of the midpoint L2E coupling structure in an embodiment of the present invention; Figure 8 This is a schematic diagram of the velocity inversion results in an embodiment of the present invention; Figure 9 This is a schematic diagram of the smoke control results in an embodiment of the present invention; Figure 10 This is a schematic diagram of boundary control parameters in an embodiment of the present invention. Detailed Implementation

[0025] The present invention will now be described in detail with reference to specific embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that those skilled in the art can make several changes and improvements without departing from the concept of the present invention. These all fall within the protection scope of the present invention.

[0026] This invention discloses an efficient differentiable fluid optimization method and system based on flow mapping and surrogate gradients. The method includes: acquiring the boundary conditions and current control variables of the fluid scenario to be optimized, and obtaining a preset objective function; based on the boundary conditions and current control variables, constructing a forward simulation computation graph using the flow mapping method, dividing the overall simulation time into multiple renormalization cycles, each renormalization cycle sequentially completing three stages: velocity advancement, flow mapping evolution, and renormalization; obtaining the final velocity field through multi-cycle iterative simulation advancement, and quantitatively evaluating the simulation process in conjunction with the objective function, calculating the objective function value corresponding to the current control variable; relying on the forward... This invention constructs a proxy gradient computation graph into the simulation computation graph, distinguishing between the first adjoint path (the dominant gradient path) formed by the pullback operator and the second adjoint path (the auxiliary gradient path) of the gradient propagation of intermediate variables in the flow mapping. A proxy gradient is constructed by combining Eulerian-Eulerian gradient propagation, long-step Lagrange-Lagrange gradient propagation, and midpoint Lagrange-Eulerian gradient propagation methods to obtain the accurate proxy gradient of the current control variable. Based on the proxy gradient, the current control variable is iteratively updated until the objective function converges. The converged control variable is used as the optimization result and converted into executable velocity field data output, achieving precise control of the fluid simulation scenario. This invention can significantly reduce the computational cost and memory usage of backpropagation in differentiable fluid optimization while preserving high-fidelity gradient information and ensuring optimization accuracy. It is widely applicable to complex fluid optimization tasks with long time spans and high vortex dynamics, such as velocity inversion, smoke control, and aerodynamic shape optimization. This invention aims to solve the problems of high computational cost and difficulty in balancing gradient accuracy and computational efficiency in existing differentiable fluid optimization techniques, effectively reducing the computational cost and memory usage of the differentiable fluid simulation optimization process while ensuring gradient fidelity.

[0027] Example 1: like Figure 2 As shown, this embodiment provides an efficient differentiable fluid optimization method based on flow mapping and surrogate gradients, including the following steps: Step S1: Collect the boundary conditions and current control variables of the fluid scenario to be optimized through physical sensors, and obtain the objective function through the user input interface.

[0028] In this embodiment, the fluid scenario to be optimized is a velocity inversion scenario, the physical sensors are a flow velocity sensor and a pressure sensor, the current control variable is the initial velocity field, and the objective function is the final velocity difference.

[0029] The fluid scenarios to be optimized also include smoke control scenarios or shape optimization scenarios; current control variables also include keyframe velocity fields or boundary control quantities, which include object shape parameters, polar coordinate radius parameters, NACA airfoil parameters, angle of attack parameters, building orientation parameters, or cross-sectional spline control parameters. The objective function also includes vorticity differences, smoke density differences, object surface pressure, drag, lift-to-drag ratio, or wake vortex structure differences. These parameters are derived from flow velocity sensors, pressure sensors, smoke or image acquisition equipment, wind tunnel or water tank experimental equipment, computer-aided design models, or historical simulation databases.

[0030] This embodiment is executed by a computing device including a processor, memory, and a graphics processor, and the computing device is connected to physical sensors.

[0031] Step S2: Based on the boundary conditions and the current control variables, construct the forward simulation calculation graph using the flow mapping method, execute forward propagation covering the simulation time, obtain the endpoint velocity field, quantitatively evaluate the endpoint velocity field and the intermediate states saved during the forward propagation process based on the endpoint velocity difference, and calculate the objective function value under the current control variables.

[0032] Figure 3 This is a schematic diagram of the forward propagation computation in an embodiment of the present invention, used to illustrate the three stages of velocity advancement, flow mapping evolution and renormalization within each renormalization cycle, and how the two data paths PASS1 and PASS2 jointly generate the velocity field at the end of the cycle.

[0033] In this embodiment, the simulation time is divided into K sequentially connected renormalization cycles. Each renormalization cycle contains N consecutive time steps. Each renormalization cycle corresponds to a cycle start velocity field and a cycle end velocity field. Within each renormalization cycle, the velocity propagation stage, the flow mapping evolution stage, and the renormalization stage are passed sequentially. The initial velocity field is used as the cycle start velocity field of the first renormalization cycle, and the cycle end velocity field of the previous renormalization cycle is used as the cycle start velocity field of the next renormalization cycle. The velocity field is updated cycle by cycle until the cycle end velocity field of the last renormalization cycle is reached, which is used as the end velocity field.

[0034] Specifically, step S2 includes the following sub-steps: Step S2.1: Based on the acquired initial velocity field and boundary conditions, establish the global variables required for fluid simulation in the memory of the computing device. The global variables include the velocity field stored as an interlaced MAC mesh. Pressure fields stored as regular grids Fluid volume fraction Forward Flow Mapping Reverse flow mapping Forward flow mapping Jacobian matrix and the Jacobian matrix of the reverse flow mapping (The variables are distinguished by superscripts at different times).

[0035] Step S2.2: If it is the first renormalization cycle, then the initial velocity field is taken as the starting velocity field of the current renormalization cycle. Otherwise, the velocity field at the end of the previous renormalization cycle is taken as the velocity field at the beginning of the current renormalization cycle. At the same time, the pressure field is set to zero, and the fluid volume fraction is set according to the boundary conditions. And map the forward flow at the start of the current renormalization cycle. and reverse flow mapping All are initialized to the identity mapping of the grid positions, and the corresponding forward flow is mapped to the Jacobian matrix. and the Jacobian matrix of the reverse flow mapping All are initialized to the identity matrix.

[0036] Step S2.3: In the velocity propagation phase, perform semi-Lagrangian velocity advection and pressure projection on the velocity field at the starting point of the current renormalization cycle, and generate the midpoint velocity sequence from time step 0 to N-1 within the cycle step by step. Pressure projection is used to constrain the velocity field to a divergence-free velocity space and obtain the pressure field. .

[0037] Step S2.4, in the flow mapping evolution stage, utilize the midpoint velocity sequence The forward flow mapping that drives the current renormalization cycle initiation time Reverse flow mapping Forward flow mapping Jacobian matrix and the Jacobian matrix of the reverse flow mapping Evolution is used to record the positional correspondence and local deformation of fluid particles from the beginning to the end of the cycle and from the end of the cycle back to the beginning of the cycle.

[0038] Step S2.5: In the renormalization stage, the velocity field at the beginning of the current renormalization cycle is interpolated and sampled using the reverse flow mapping at the end of the current renormalization cycle. The sampled velocity vector is stretched using the corresponding reverse flow mapping Jacobian matrix to obtain the velocity field before projection at the end of the cycle.

[0039] Figure 4 This is a schematic diagram of flow mapping transport in an embodiment of the present invention, used to illustrate how forward flow mapping, reverse flow mapping and their Jacobian describe the positional correspondence and local deformation relationship of fluid particles at different times.

[0040] Step S2.6: Perform pressure projection on the velocity field before projection to obtain the velocity field at the end of the current renormalization period that satisfies the incompressible constraint. .

[0041] The velocity field at the end of the period is expressed by formula (1): (1) In the formula, represents the interpolation matrix at the reverse flow mapping position, and ⊙ represents the block product performed on each discrete velocity vector. This represents the velocity field at the starting point of the current renormalization cycle. This represents the velocity field at the end of the current renormalization cycle. The Jacobian matrix represents the reverse flow mapping of the current renormalization cycle starting point.

[0042] Step S2.7: Determine whether the current renormalization cycle is the last cycle. If not, then set the velocity field at the end of the current renormalization cycle. Store the velocity field in the buffer as the starting velocity field for the next cycle, and return to step S2.2 to continue execution; if so, terminate the loop, take the end velocity field of the current renormalization cycle as the end velocity field of the entire simulation, and quantitatively evaluate the end velocity field and the midpoint velocity sequence, pressure field sequence, forward flow mapping, reverse flow mapping, forward flow mapping Jacobian matrix and reverse flow mapping Jacobian matrix saved during the simulation according to the objective function, and calculate the objective function value under the current control variables.

[0043] Step S3: Based on the forward simulation computation graph, build and execute the surrogate gradient computation graph. With the objective function value as the optimization target for backpropagation, take the gradient of the objective function with respect to the final velocity field as the starting point of backpropagation, and propagate forward layer by layer along the surrogate gradient computation graph. After the surrogate gradient is constructed, the surrogate gradient at the current control variable is obtained.

[0044] Based on the gradients of each computational path in the forward simulation computation graph, a gradient propagation path (also known as an adjoint path) is constructed in the surrogate gradient computation graph. This surrogate gradient propagation path is divided into a first adjoint path and a second adjoint path. The first adjoint path is a step propagation path, formed by a pullback operator, used to propagate the gradient of the velocity field at the end of the period to the velocity field at the beginning of the period via the adjoint pullback operator. The second adjoint path is a path that propagates stepwise through forward flow mapping, backward flow mapping, the forward flow mapping Jacobian matrix, the backward flow mapping Jacobian matrix, and the midpoint velocity sequence, used to provide auxiliary gradient correction.

[0045] Figure 5This is a schematic diagram comparing the complete gradient and surrogate gradient computation graphs in an embodiment of the present invention. It shows that backpropagation is divided into a dominant gradient path PASS1 and an auxiliary gradient path PASS2, and the surrogate gradient reduces the computational load by retaining the dominant path and simplifying the high-overhead auxiliary path.

[0046] Furthermore, the second adjoint path includes Eulerian-Eulerian E2E gradient propagation, Lagrange-Lagrange L2L gradient propagation, and Lagrange-Eulerian L2E gradient propagation. Eulerian-Eulerian gradient propagation is used to propagate gradients between Eulerian velocity grids, preserving time-step propagation. Lagrange-Lagrange gradient propagation is used to propagate gradients between flow map variables and Jacobian variables, replacing time-step propagation with long-step propagation spanning one renormalization period. The long-step propagation propagates the flow map gradient and flow map Jacobian gradient with a time step length corresponding to the period length of the renormalization period. Lagrange-Eulerian gradient propagation is used to couple the gradients of the flow map variables and Jacobian variables back to the Eulerian velocity grid. Coupling is performed at the midpoint of the renormalization period, allowing the gradients of the flow map variables and Jacobian variables to propagate to the midpoint velocity gradient via midpoint Lagrange-Eulerian coupling, and then to the periodic starting velocity field of the renormalization period via stepwise Eulerian-Eulerian propagation.

[0047] Based on the physical duality of different adjoint operators in the first and second adjoint paths, long-step propagation across the renormalization period is used to replace the time-step propagation in the second adjoint path. That is, a surrogate gradient computation graph is constructed based on the physical duality of three types of adjoint operators: Euler-Euler propagation, Lagrange-Lagrange propagation, and Lagrange-Euler propagation. The surrogate gradient computation graph retains the dominant gradient path and simplifies the high-overhead auxiliary gradient path. The gradient is propagated along the surrogate gradient computation graph to obtain the surrogate gradient of the current control variable.

[0048] Step S3 includes the following sub-steps: Step S3.1: Taking the last renormalization period as the current renormalization period, calculate the final velocity gradient of the objective function J at the velocity field at the end of the current renormalization period, denoted as... This serves as the starting point for backpropagation.

[0049] Step S3.2, calculate the final velocity gradient of the current renormalization cycle. The velocity gradient propagates to the pre-projection velocity gradient via the adjoint form of the pressure projection operator. This allows the reverse velocity gradient to satisfy the incompressible constraint, resulting in the velocity gradient before projection at the end of the current renormalization cycle.

[0050] Step S3.3, execute the first accompanying path PASS1: project the velocity gradient before the end of the current renormalization cycle. The dominant gradient of the current renormalization period is obtained by directly propagating the velocity field to the starting point of the current renormalization period via the accompanying pullback operator. The expression is as follows: (2) In the formula, J represents the objective function. This represents the dominant gradient of the velocity field at the start of the current renormalization period, generated by the pullback operator S. represents the interpolation matrix at the reverse flow mapping position at the end of the current renormalization cycle, and ⊙ represents the block product performed on each discrete velocity vector.

[0051] Step S3.4, generate the initial gradient of the second adjoint path PASS2: the dominant gradient of the current renormalization cycle. The backflow map gradient propagated to the end of the current renormalization cycle via the accompanying pullback operator is denoted as... The Jacobian gradient of the backflow at the end of the current renormalization cycle is denoted as . , which serves as the input for auxiliary gradient correction.

[0052] Step S3.5, perform long-step L2L propagation: Under the constraint of Courant-Friedrich-Liouville (CFL), the gradient of the backflow map and the gradient of the backflow map Jacobian matrix are propagated across the current renormalization period with long steps. The N time-step L2L propagation is replaced by long-step propagation with step size NΔt. The dominant term approximation formula is as follows: (3) (4) in, This represents the gradient of the Jacobian matrix of the backflow mapping at the (n+1)th time level; This represents the gradient of the backflow mapping at the (n+1)th time level; This represents the reverse flow mapping of the (n+1)th time layer; This represents the Jacobian matrix of the reverse flow mapping at the (n+1)th time level; This represents the adjoint gradient contribution corresponding to the flow map evolution operator; This represents the adjoint gradient contribution corresponding to the Jacobian evolution operator; Let represent the identity matrix, indicating that the dominant gradient in L2L propagation under the Courant-Friedrich-Liouville constraint (CFL) is approximately a propagation that remains invariant along the Lagrangian variables. Through this step, the gradient is back-propagated across the current renormalization period, yielding the Lagrangian gradient after propagation across the current renormalization period.

[0053] Figure 6 This is a schematic diagram of the proxy gradient structure of long-step L2L and step-by-step E2E in an embodiment of the present invention, illustrating the reason and structure for merging the continuous L2L backpropagation related to the flow mapping into long-step propagation across the renormalization period, while retaining the step-by-step propagation of the velocity field E2E.

[0054] Step S3.6, perform midpoint L2E coupling: at the midpoint of the current renormalization cycle, couple the Lagrange variable gradient back to the Eulerian velocity grid to form the midpoint velocity gradient of the current renormalization cycle.

[0055] Figure 7 This is a schematic diagram of the midpoint L2E coupling structure in an embodiment of the present invention, showing that the flow map gradient and the flow map Jacobian gradient are coupled back to the Eulerian velocity grid at the midpoint of the renormalization period, and then propagated to the velocity field at the beginning of the period through step-by-step E2E.

[0056] Step S3.7, perform stepwise E2E propagation: propagate the midpoint velocity gradient of the current renormalization cycle in reverse time step order to the starting point of the current renormalization cycle. The single-step E2E propagation is expressed as follows: (5) In the formula, Let B represent the Eulerian velocity field used for velocity advection before and after the pressure projection at the nth time step; B represents the adjoint gradient contribution corresponding to the velocity advection operator. Let represent the midpoint velocity field after projection obtained through the semi-Lagrange advection at the nth time step; Let represent the velocity field at the midpoint before projection obtained by the semi-Lagrange advection at the nth time step; The location of the semi-Lagrange backtracking source point is typically determined by the grid position and the current velocity field, time-stepped. Backtracking yielded the result; This represents the interpolation matrix used to sample the velocity or gradient from the Euler grid at the backtracking location; This represents the interpolation space gradient operator, used to describe the impact of backtracking position changes on velocity sampling results; This indicates the time step of a single simulation.

[0057] Through this step, the gradient is finally transmitted back to the velocity field at the starting point of the current renormalization cycle, and the velocity gradient at the starting point of the current renormalization cycle is output as an auxiliary gradient correction for the second adjoint path.

[0058] Step S3.8, gradient accumulation and cross-cycle propagation: The dominant gradient of the current renormalization cycle and the auxiliary gradient correction of the second adjoint path are accumulated to obtain the surrogate gradient at the starting point of the current renormalization cycle. Step S3.9: Determine whether the current renormalization period is the first renormalization period. If yes, terminate the loop and the surrogate gradient at the start of the current renormalization period is the surrogate gradient of the current control variable. If no, assign the surrogate gradient at the start of the current renormalization period to the gradient of the velocity field at the end of the previous renormalization period, update the current renormalization period to the previous renormalization period, and return to step S3.2 to continue execution.

[0059] Step S4: Based on the surrogate gradient at the current control variable, update the current control variable, and use the updated current control variable as the input for the next iteration. Repeat steps S2 to S4 until the objective function converges, and use the current control variable corresponding to the convergence as the optimization result.

[0060] Specifically, the Adam optimizer or gradient descent optimizer is used to update the current control variable. In each iteration, step S2 provides the objective function value of the current round to step S4. Step S4 compares the objective function value of the current round with the objective function value of the previous round. When the absolute change or relative change is less than a preset threshold, the objective function is determined to have converged.

[0061] Step S5: Convert the optimization results into executable velocity field data and output the executable velocity field data to the fluid simulation display device to control the fluid simulation display device to display the initial velocity field that can evolve into the target vortex structure.

[0062] Specifically, the application of the optimization results varies depending on the scenario. In velocity inversion scenarios, the optimization results are used to generate an initial velocity field that can evolve into the target vortex structure. In smoke control scenarios, the optimization results are used to generate velocity control data that drives the smoke density field to the target keyframe. In shape optimization scenarios, the optimization results are used to generate boundary control parameters for solids, which can be used by CNC machining, 3D modeling, wind tunnel verification, or fluid control equipment. For non-convex objective functions such as smoke control and velocity inversion, a multi-stage optimization strategy from coarse to fine can be adopted. The computing device converts the current control variables into executable data and outputs it to the aforementioned devices.

[0063] Example 2: This embodiment provides a velocity inversion application method. The control variable is the initial velocity field. The system reads the target vorticity field or target velocity field through memory, and calculates the loss between the final state and the target state after several simulation steps. The system uses the surrogate gradient from Embodiment 1 to calculate the gradient of the initial velocity field, and updates the initial velocity through a gradient-based optimizer. The converged initial velocity field can be output to a fluid simulation control device or a visualization display device to generate long-term evolution targets such as checkerboard vortices and trilobed vortex rings.

[0064] Figure 8This is a schematic diagram of the velocity inversion results in an embodiment of the present invention, showing that the optimized initial velocity field can reconstruct the checkerboard vortex target.

[0065] Example 3: Figure 9 This is a schematic diagram of the smoke control results in an embodiment of the present invention, showing the smoke morphology evolution effect driven by the keyframe velocity field.

[0066] This embodiment provides a smoke control application. The control variable is one or more keyframe velocity fields. The system performs semi-Lagrange advection on a high-resolution smoke density grid and calculates the smoke density distribution difference at the target keyframe. During backpropagation, the system uses aggregation to achieve E2E adjoint advection of the high-resolution passive scalar field, and scattering to achieve L2E coupling related to flow mapping, thereby improving GPU execution efficiency. The converged keyframe velocity fields can be output to fluid animation generation devices, display devices, or special effects production systems to drive the smoke morphology to evolve according to the target keyframe.

[0067] Example 4: Figure 10 This is a schematic diagram of boundary control parameters in an embodiment of the present invention, showing how solid boundaries such as bullet shape and wing shape are represented by low-dimensional geometric parameters and further rasterized into volume fraction fields that can be used by the fluid solver.

[0068] This embodiment provides a shape optimization application method. The control variables are solid boundary control quantities. For example, the shape of a bullet can be described by multiple cross-sectional radius parameters, an airfoil can be described by NACA parameters and angle of attack, and a building complex can be described by building orientation parameters. The system rasterizes the boundary control parameters into fluid volume fraction fields and couples them with the fluid solver through pressure projection under the Cut-Cell strategy. The objective function may include drag, lift-to-drag ratio, pressure, wake vortex index, etc. The converged boundary control quantities can be output to 3D modeling software, CNC machining equipment, wind tunnel experimental control system, or engineering design database for subsequent physical manufacturing, experimental verification, or engineering simulation.

[0069] Example 5: The present invention also provides an efficient differentiable fluid optimization system based on flow mapping and surrogate gradient. The efficient differentiable fluid optimization system based on flow mapping and surrogate gradient can be implemented by executing the process steps of the efficient differentiable fluid optimization method based on flow mapping and surrogate gradient. That is, those skilled in the art can understand the efficient differentiable fluid optimization method based on flow mapping and surrogate gradient as a preferred embodiment of the efficient differentiable fluid optimization system based on flow mapping and surrogate gradient.

[0070] like Figure 1 As shown, this efficient differentiable fluid optimization system based on flow mapping and surrogate gradients includes: The data acquisition module collects the boundary conditions and current control variables of the fluid scenario to be optimized through physical sensors, and obtains the objective function through the user input interface; The forward simulation module constructs a forward simulation computation graph based on boundary conditions and current control variables using the flow mapping method, executes forward propagation covering the simulation time, obtains the final velocity field, quantitatively evaluates the final velocity field and the intermediate states saved during the forward propagation process based on the difference in the final velocity, and calculates the objective function value under the current control variables. The proxy gradient module, based on the forward simulation computation graph, builds and executes the proxy gradient computation graph. The objective function value is used as the optimization target for backpropagation. The gradient of the objective function with respect to the final velocity field is used as the starting point for backpropagation. It propagates forward layer by layer along the proxy gradient computation graph. After the proxy gradient is constructed, the proxy gradient at the current control variable is obtained. The parameter update module updates the current control variable based on the surrogate gradient at the current control variable, and uses the updated current control variable as the input for the next iteration. Steps S2 to S4 are repeated until the objective function converges, and the current control variable at the convergence point is used as the optimization result. The output module converts the optimization results into executable velocity field data and outputs the executable velocity field data to the fluid simulation display device to control the fluid simulation display device to display the initial velocity field that can evolve into the target vortex structure.

[0071] Specifically, the forward simulation module includes the following sub-modules: Module M2.1, based on the acquired initial velocity field and boundary conditions, establishes the global variables required for fluid simulation in the memory of the computing device. The global variables include the velocity field stored as an interlaced MAC mesh. Pressure fields stored as regular grids Fluid volume fraction Forward Flow Mapping Reverse flow mapping Forward flow mapping Jacobian matrix and the Jacobian matrix of the reverse flow mapping (The variables are distinguished by superscripts at different times).

[0072] Module M2.2, if it is the first renormalization cycle, then uses the initial velocity field as the starting velocity field of the current renormalization cycle. Otherwise, the velocity field at the end of the previous renormalization cycle is taken as the velocity field at the beginning of the current renormalization cycle. At the same time, the pressure field is set to zero, and the fluid volume fraction is set according to the boundary conditions. And map the forward flow at the start of the current renormalization cycle. and reverse flow mapping All are initialized to the identity mapping of the grid positions, and the corresponding forward flow is mapped to the Jacobian matrix. and the Jacobian matrix of the reverse flow mapping All are initialized to the identity matrix.

[0073] Module M2.3, during the velocity propagation phase, performs semi-Lagrangian velocity advection and pressure projection on the velocity field at the starting point of the current renormalization cycle, generating a midpoint velocity sequence step-by-step. Pressure projection is used to constrain the velocity field to a divergence-free velocity space and obtain the pressure field. .

[0074] Module M2.4 utilizes the midpoint velocity sequence during the flow mapping evolution stage. The forward flow mapping that drives the current renormalization cycle initiation time Reverse flow mapping Forward flow mapping Jacobian matrix and the Jacobian matrix of the reverse flow mapping Evolution is used to record the positional correspondence and local deformation of fluid particles from the beginning to the end of the cycle and from the end of the cycle back to the beginning of the cycle.

[0075] In module M2.5, during the renormalization phase, the velocity field at the beginning of the current renormalization cycle is interpolated and sampled using the reverse flow mapping at the end of the current renormalization cycle. The sampled velocity vector is then stretched using the corresponding reverse flow mapping Jacobian matrix to obtain the velocity field before projection at the end of the cycle.

[0076] Module M2.6 performs pressure projection on the velocity field before projection to obtain the velocity field at the end of the current renormalization period that satisfies the incompressibility constraint. .

[0077] The velocity field at the end of the period is expressed by formula (1): (1) In the formula, represents the interpolation matrix at the reverse flow mapping position, and ⊙ represents the block product performed on each discrete velocity vector. This represents the velocity field at the starting point of the current renormalization cycle. This represents the velocity field at the end of the current renormalization cycle. The Jacobian matrix represents the reverse flow mapping of the current renormalization cycle starting point.

[0078] Module M2.7 determines whether the current renormalization cycle is the last cycle. If not, it sets the velocity field at the end of the current renormalization cycle. Store the velocity field in the buffer as the starting velocity field for the next cycle, and return to module M2.2 to continue execution; otherwise, terminate the loop, take the end velocity field of the current renormalization cycle as the end velocity field of the entire simulation, and quantitatively evaluate the end velocity field and the midpoint velocity sequence, pressure field sequence, forward flow mapping, reverse flow mapping, forward flow mapping Jacobian matrix and reverse flow mapping Jacobian matrix saved during the simulation according to the objective function, and calculate the objective function value under the current control variables.

[0079] The proxy gradient module includes the following sub-modules: Module M3.1 takes the last renormalization period as the current renormalization period and calculates the final velocity gradient of the objective function J at the velocity field at the end of the current renormalization period, denoted as... This serves as the starting point for backpropagation.

[0080] Module M3.2 calculates the final velocity gradient of the current renormalization cycle. The velocity gradient propagates to the pre-projection velocity gradient via the adjoint form of the pressure projection operator. This allows the reverse velocity gradient to satisfy the incompressible constraint, resulting in the velocity gradient before projection at the end of the current renormalization cycle.

[0081] Module M3.3 executes the first accompanying path PASS1: projecting the velocity gradient at the end of the current renormalization cycle. The dominant gradient of the current renormalization period is obtained by directly propagating the velocity field to the starting point of the current renormalization period via the accompanying pullback operator. The expression is as follows: (2) In the formula, J represents the objective function. This represents the dominant gradient of the velocity field at the start of the current renormalization period, generated by the pullback operator S. represents the interpolation matrix at the reverse flow mapping position at the end of the current renormalization cycle, and ⊙ represents the block product performed on each discrete velocity vector.

[0082] Module M3.4 generates the initial gradient of the second adjoint path PASS2: the dominant gradient of the current renormalization cycle. The backflow map gradient propagated to the end of the current renormalization cycle via the accompanying pullback operator is denoted as... The Jacobian gradient of the backflow at the end of the current renormalization cycle is denoted as . , which serves as the input for auxiliary gradient correction.

[0083] Module M3.5 performs long-step L2L propagation: Under the constraint of Courant-Friedrich-Liouville (CFL), the gradient of the backflow map and the gradient of the backflow map Jacobian matrix are propagated across the current renormalization period with long steps. The N time-step L2L propagation is replaced by long-step propagation with a step size of NΔt. The dominant term approximation formula is as follows: (3) (4) in, This represents the gradient of the Jacobian matrix of the backflow mapping at the (n+1)th time level; This represents the gradient of the backflow mapping at the (n+1)th time level; This represents the reverse flow mapping of the (n+1)th time layer; This represents the Jacobian matrix of the reverse flow mapping at the (n+1)th time level; This represents the adjoint gradient contribution corresponding to the flow map evolution operator; This represents the adjoint gradient contribution corresponding to the Jacobian evolution operator; Let represent the identity matrix, indicating that the dominant gradient in L2L propagation under the Courant-Friedrich-Liouville constraint (CFL) is approximately a propagation that remains invariant along the Lagrangian variables. Through this module, the gradient is back-propagated across the current renormalization period, yielding the Lagrangian gradient after propagation across the current renormalization period.

[0084] Module M3.6 performs midpoint L2E coupling: at the midpoint of the current renormalization cycle, the Lagrange variable gradient is coupled back to the Eulerian velocity grid to form the midpoint velocity gradient of the current renormalization cycle.

[0085] Module M3.7 performs stepwise E2E propagation: the midpoint velocity gradient of the current renormalization cycle is propagated stepwise to the starting point of the current renormalization cycle in reverse time step order. The single-step E2E propagation is expressed as follows: (5) In the formula, Let B represent the Eulerian velocity field used for velocity advection before and after the pressure projection at the nth time step; B represents the adjoint gradient contribution corresponding to the velocity advection operator. Let represent the midpoint velocity field after projection obtained through the semi-Lagrange advection at the nth time step; Let represent the velocity field at the midpoint before projection obtained by the semi-Lagrange advection at the nth time step; The location of the semi-Lagrange backtracking source point is typically determined by the grid position and the current velocity field, time-stepped. Backtracking yielded the result; This represents the interpolation matrix used to sample the velocity or gradient from the Euler grid at the backtracking location; This represents the interpolation space gradient operator, used to describe the impact of backtracking position changes on velocity sampling results; This indicates the time step of a single simulation.

[0086] Through this module, the gradient is finally transmitted back to the velocity field at the starting point of the current renormalization cycle, and the velocity gradient at the starting point of the current renormalization cycle is output as an auxiliary gradient correction for the second adjoint path.

[0087] Module M3.8, Gradient Accumulation and Cross-Period Transmission: Accumulates the dominant gradient of the current renormalization cycle and the auxiliary gradient correction of the second adjoint path to obtain the surrogate gradient at the starting point of the current renormalization cycle. Module M3.9 determines whether the current renormalization period is the first renormalization period. If so, the loop terminates, and the surrogate gradient at the start of the current renormalization period becomes the surrogate gradient of the current control variable. If not, the surrogate gradient at the start of the current renormalization period is assigned to the gradient of the velocity field at the end of the previous renormalization period, and the current renormalization period is updated to the previous renormalization period. Then, the program returns to module M3.2 to continue execution.

[0088] It should be noted that the division of the various modules in the above system is merely a logical functional division. In actual implementation, they can be fully or partially integrated into a single physical entity, or they can be physically separated. Furthermore, these modules can be implemented entirely in software through processing element calls; they can be implemented entirely in hardware; or some modules can be implemented by processing element calls to software, while others are implemented in hardware.

[0089] For example, these modules can be one or more integrated circuits configured to implement the above methods, such as one or more ASICs, one or more DSPs, or one or more FPGAs. As another example, when a module is implemented through processing element scheduler code, the processing element can be a general-purpose processor (e.g., a CPU), a GPU, or other processor capable of calling program code for parallel or conventional computation. Furthermore, these modules can be integrated together to form a System-on-a-Chip (SoC).

[0090] In summary, this invention reduces the computational cost and memory usage of backpropagation while maintaining near-perfect gradient accuracy, making it suitable for long-duration, high-vortex dynamics fluid optimization tasks such as velocity inversion, smoke control, and shape optimization. Therefore, this invention effectively overcomes the various shortcomings of existing technologies and has high industrial applicability.

[0091] Those skilled in the art will understand that, besides implementing the system and its various devices, modules, and units provided by this invention in the form of purely computer-readable program code, the same functions can be achieved entirely through logical programming of the method steps, making the system and its various devices, modules, and units of this invention function in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers. Therefore, the system and its various devices, modules, and units provided by this invention can be considered as a hardware component, and the devices, modules, and units included therein for implementing various functions can also be considered as structures within the hardware component; alternatively, the devices, modules, and units for implementing various functions can be considered as both software modules implementing the method and structures within the hardware component.

[0092] Specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art can make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. Unless otherwise specified, the embodiments and features described in this application can be arbitrarily combined with each other.

Claims

1. An efficient differentiable fluid optimization method based on flow mapping and surrogate gradient, characterized in that, Includes the following steps: Step S1: Collect the boundary conditions and current control variables of the fluid scenario to be optimized through physical sensors, and obtain the objective function through the user input interface; Step S2: Based on the boundary conditions and the current control variables, construct a forward simulation computation graph using the flow mapping method, execute forward propagation covering the simulation time, obtain the endpoint velocity field, quantitatively evaluate the endpoint velocity field and the intermediate states saved during the forward propagation process based on the endpoint velocity difference, and calculate the objective function value under the current control variables. Step S3: Based on the forward simulation calculation graph, construct and execute the proxy gradient calculation graph. With the objective function value as the optimization target for backpropagation, take the gradient of the objective function with respect to the endpoint velocity field as the starting point of backpropagation, and propagate forward layer by layer along the proxy gradient calculation graph. After the proxy gradient is constructed, the proxy gradient at the current control variable is obtained. Step S4: Based on the surrogate gradient at the current control variable, update the current control variable, and repeat steps S2 to S4 until the objective function converges. The current control variable at the time of convergence is taken as the optimization result. Step S5: Convert the optimization result into executable velocity field data, and output the executable velocity field data to the fluid simulation display device to control the fluid simulation display device to display the initial velocity field that can evolve into the target vortex structure.

2. The efficient differentiable fluid optimization method based on flow mapping and surrogate gradient according to claim 1, characterized in that, In step S1, the fluid scenario to be optimized is a velocity inversion scenario, the physical sensors are a flow velocity sensor and a pressure sensor, the current control variable is the initial velocity field, and the objective function is the final velocity difference.

3. The efficient differentiable fluid optimization method based on flow mapping and surrogate gradient according to claim 1, characterized in that, In step S2, the simulation time is divided into K sequentially connected renormalization cycles. Each renormalization cycle contains N consecutive time steps. Each renormalization cycle corresponds to a cycle start velocity field and a cycle end velocity field. Within each renormalization cycle, a velocity propagation stage, a flow mapping evolution stage, and a renormalization stage are passed sequentially. The initial velocity field is used as the cycle start velocity field of the first renormalization cycle, and the cycle end velocity field of the previous renormalization cycle is used as the cycle start velocity field of the next renormalization cycle. The velocity field is updated cycle by cycle until the cycle end velocity field of the last renormalization cycle is reached, which is used as the end velocity field.

4. The efficient differentiable fluid optimization method based on flow mapping and surrogate gradient according to claim 3, characterized in that, Step S2 includes the following steps: Step S2.1: Based on the initial velocity field and the boundary conditions, establish the global variables required for fluid simulation. The global variables include the velocity field stored as an interlaced MAC grid, the pressure field stored as a regular grid, the fluid volume fraction, the forward flow mapping, the reverse flow mapping, the forward flow mapping Jacobian matrix, and the reverse flow mapping Jacobian matrix. Step S2.2: If it is the first renormalization cycle, the initial velocity field is taken as the starting velocity field of the current renormalization cycle; otherwise, the ending velocity field of the previous renormalization cycle is taken as the starting velocity field of the current renormalization cycle. At the same time, the pressure field is set to zero, the fluid volume fraction is set according to the boundary conditions, and the forward flow mapping and reverse flow mapping at the starting time of the current renormalization cycle are initialized to the identity mapping of the grid position. The corresponding forward flow mapping Jacobian matrix and reverse flow mapping Jacobian matrix are initialized to the identity matrix. Step S2.3: In the velocity propagation phase, a semi-Lagrange velocity advection and pressure projection are performed on the velocity field at the starting point of the current renormalization cycle, and a midpoint velocity sequence is generated step by step. The pressure projection is used to constrain the velocity field to a divergence-free velocity space and obtain the pressure field. Step S2.4: In the flow mapping evolution stage, the forward flow mapping, reverse flow mapping, forward flow mapping Jacobian matrix and reverse flow mapping Jacobian matrix at the starting moment of the current renormalization cycle are driven by the midpoint velocity sequence to record the positional correspondence and local deformation of the fluid particles from the starting moment of the current renormalization cycle to the ending moment of the cycle and from the ending moment of the cycle back to the starting moment of the cycle. Step S2.5: In the renormalization stage, the velocity field at the beginning of the current renormalization cycle is interpolated and sampled using the reverse flow mapping at the end of the current renormalization cycle, and the sampled velocity vector is stretched using the corresponding reverse flow mapping Jacobian matrix to obtain the velocity field before projection at the end of the current renormalization cycle. Step S2.6: Perform pressure projection on the velocity field before projection to obtain the velocity field at the end of the current renormalization cycle; Step S2.7: Determine whether the current renormalization cycle is the last cycle. If not, store the end velocity field of the current renormalization cycle into the buffer as the starting velocity field of the next cycle, and return to step S2.2 to continue execution. If so, the loop is terminated, and the velocity field at the end of the current renormalization cycle is taken as the end velocity field of the entire simulation. Based on the difference in the end velocity, the end velocity field and the midpoint velocity sequence, pressure field sequence, forward flow mapping, reverse flow mapping, forward flow mapping Jacobian matrix and reverse flow mapping Jacobian matrix saved during the simulation are quantitatively evaluated, and the objective function value under the current control variables is calculated.

5. The efficient differentiable fluid optimization method based on flow mapping and surrogate gradient according to claim 1, characterized in that, In step S3, based on the gradients of each computation path in the forward simulation computation graph, a gradient propagation path in the proxy gradient computation graph is constructed. The gradient propagation path is divided into a first adjoint path and a second adjoint path. The first adjoint path is a step propagation path, which forms a step connection by a pullback operator. It is used to propagate the gradient of the velocity field at the end of the period to the velocity field at the beginning of the period via the adjoint pullback operator. The second adjoint path is a path that propagates step by step through forward flow mapping, reverse flow mapping, forward flow mapping Jacobian matrix, reverse flow mapping Jacobian matrix and midpoint velocity sequence. It is used to provide auxiliary gradient correction.

6. The efficient differentiable fluid optimization method based on flow mapping and surrogate gradient according to claim 5, characterized in that, In step S3, the second accompanying path includes Eulerian-Eulerian gradient propagation, Lagrange-Lagrange gradient propagation, and Lagrange-Eulerian gradient propagation; The Euler-Euler gradient propagation is used to propagate gradients between Euler velocity grids, preserving time-step propagation; The Lagrange-Lagrange gradient propagation is used to propagate gradients between the flow map variables and the Jacobian variables, replacing time-step propagation with long-step propagation spanning one of the renormalization periods. The long-step propagation propagates the flow map gradient and the flow map Jacobian gradient with a time step size corresponding to the period length of the renormalization period. The Lagrange-Euler gradient propagation is used to couple the gradients of the flow mapping variable and the Jacobian variable back to the Euler velocity grid. The coupling is performed at the midpoint of the renormalization period, so that the gradients of the flow mapping variable and the Jacobian variable are propagated to the midpoint velocity gradient via the midpoint Lagrange-Euler coupling, and then propagated to the velocity field at the beginning of the renormalization period via the stepwise Euler-Euler coupling.

7. The efficient differentiable fluid optimization method based on flow mapping and surrogate gradient according to claim 6, characterized in that, Step S3 includes the following steps: Step S3.1: Take the last renormalization cycle as the current renormalization cycle, and at the end velocity field of the current renormalization cycle, calculate the end velocity gradient of the end velocity difference as the starting point of back propagation. Step S3.2: Propagate the endpoint velocity gradient to the pre-projection velocity gradient through the adjoint form of the pressure projection operator, so that the reverse velocity gradient satisfies the incompressible constraint, and obtain the pre-projection velocity gradient at the end of the current renormalization cycle. Step S3.3, execute the first adjoint path: propagate the pre-projection velocity gradient to the periodic starting velocity field of the current renormalization period through the adjoint pullback operator to obtain the dominant gradient of the current renormalization period; Step S3.4: Propagate the dominant gradient to the backflow mapping gradient and backflow mapping Jacobian gradient at the end of the current renormalization cycle via the adjoint pullback operator, and use them as the initial gradient of the second adjoint path. Step S3.5: Under the condition of satisfying the Courant-Friedrich-Liouvi constraints, the gradient of the reverse flow mapping and the gradient of the reverse flow mapping Jacobian matrix are back-propagated from the end time of the current renormalization cycle to the start time along the Lagrange trajectory. The N time-step propagation is replaced by a single long-step propagation that spans N consecutive time steps to obtain the Lagrange variable gradient after the cross-cycle propagation of the current renormalization cycle. Step S3.6: At the midpoint of the current renormalization cycle, couple the Lagrange variable gradient back to the Eulerian velocity grid to form the midpoint velocity gradient of the current renormalization cycle. Step S3.7: Propagate the midpoint velocity gradient to the velocity field at the starting point of the current renormalization cycle in reverse order of time steps to obtain the velocity gradient at the starting point of the current renormalization cycle, which serves as the auxiliary gradient correction for the second adjoint path. Step S3.8: The dominant gradient and the auxiliary gradient correction of the second adjoint path are accumulated to obtain the surrogate gradient at the starting point of the current renormalization cycle. Step S3.9: Determine whether the current renormalization period is the first renormalization period. If yes, terminate the loop, and the surrogate gradient at the starting point of the current renormalization period is the surrogate gradient of the initial velocity field. If no, assign the surrogate gradient at the starting point of the current renormalization period to the gradient of the velocity field at the end point of the previous renormalization period, update the current renormalization period to the previous renormalization period, and return to step S3.2 to continue execution.

8. The efficient differentiable fluid optimization method based on flow mapping and surrogate gradient according to claim 5, characterized in that, In step S4, the current control variable is updated using the Adam optimizer or gradient descent optimizer.

9. The efficient differentiable fluid optimization method based on flow mapping and surrogate gradient according to claim 1, characterized in that, In step S1, the fluid scenario to be optimized also includes a smoke control scenario, the current control variable is the keyframe velocity field, the objective function is the smoke density difference, and the physical sensors include a smoke acquisition device and an image acquisition device.

10. A highly efficient differentiable fluid optimization system based on flow mapping and surrogate gradients, employing the highly efficient differentiable fluid optimization method based on flow mapping and surrogate gradients as described in any one of claims 1-9, characterized in that, include: The data acquisition module collects the boundary conditions and current control variables of the fluid scenario to be optimized through physical sensors, and obtains the objective function through the user input interface; The forward simulation module constructs a forward simulation computation graph based on the boundary conditions and the current control variables using the flow mapping method, performs forward propagation covering the simulation time, obtains the endpoint velocity field, quantitatively evaluates the endpoint velocity field and the intermediate states saved during the forward propagation process based on the endpoint velocity difference, and calculates the objective function value under the current control variables. The proxy gradient module, based on the forward simulation computation graph, builds and executes a proxy gradient computation graph. With the objective function value as the optimization target for backpropagation, the gradient of the objective function with respect to the final velocity field is used as the starting point for backpropagation. It propagates forward layer by layer along the proxy gradient computation graph. After the proxy gradient is constructed, the proxy gradient at the current control variable is obtained. The parameter update module updates the current control variable based on the surrogate gradient at the current control variable, uses the updated current control variable as the input for the next iteration, and repeats steps S2 to S4 until the objective function converges. The current control variable at the time of convergence is then used as the optimization result. The output module converts the optimization results into executable velocity field data and outputs the executable velocity field data to the fluid simulation display device to control the fluid simulation display device to display the initial velocity field that can evolve into the target vortex structure.

Images