A method for task decomposition and parallel solution of simulation model

Abstract

The present invention provides a method for task decomposition and parallel solution of a simulation model, which solves the two-level parallel solution problem between and within the simulation model. The method comprises: S1, extracting model differential equation description information, dynamically generating a code GPU, and establishing a solution task list; S2, the main thread traverses the task list, and dispatches a specific ordinary differential equation group solution task to the thread pool; S3, the working thread that obtains the task adopts a multi-resolution solution method to push an asynchronous parallel solution requirement to the GPU; S4, the GPU performs the solution calculation and generates a calculation task completion event; S5, the GPU task completion event is captured, and the working thread in the thread pool is awakened; S6, the working thread obtains the GPU calculation result, performs result correction, and iterates S4, S5, and S6 until the accuracy requirement is met. Through the present invention, large-scale heterogeneous ordinary differential equation groups can be processed in parallel, effectively improving the hardware resource utilization rate of the GPU during the simulation process, and achieving better load balancing.

Description

Technical Field

[0001] The present invention belongs to the application field of cross-technology of simulation and parallel computing, and in particular relates to a method for task decomposition and parallel solution oriented to a simulation model. Background Art

[0002] In many engineering fields for modeling and simulation, such as physics, chemistry, biology, etc., large-scale ordinary differential equation (ODE) initial value problem (IVP) solving is involved. Due to the use of discretization methods, under the constraint of specific solution accuracy, the local step size will be set small. In each simulation global step time period, hundreds or thousands of cumulative iterative ODE solutions are performed, consuming a lot of simulation time. Typically, for complex models such as aerospace engines, simulated combustion, and groundwater transmission, the time consumption for solving ODEs will reach more than 90% of the wall clock time, which constitutes a bottleneck for the simulation advancement speed. How to use the multi-core capabilities of modern GPUs to build an efficient model parallel solution method has become a simulation parallel computing problem that needs to be solved urgently.

[0003] Traditional parallel solution methods mainly focus on three directions:

[0004] One is to use a method similar to the transmission line model to decompose the entire simulation model into multiple components. Each component has its own solver. Through the collaborative simulation method, the operation and solution of each component are distributed to different threads to achieve the purpose of parallel solution.

[0005] The second is to extend the programming language by introducing parallel compilation instructions to provide parallel information where parallel solutions are needed. These parallel locations are mainly for loops and vector operations.

[0006] The third method is to establish a task graph based on the dependencies between codes and adopt a top-down scheduling method. Each node at the same layer in the task graph can be scheduled at the same time and dispatched to different threads.

[0007] In traditional simulation parallel solution methods, the focus is on building a multi-threaded parallel execution environment around the multi-core capabilities of the CPU, while the multi-core capabilities of the GPU are rarely considered. This is mainly due to the hardware architecture characteristics of the GPU, which can only process calculations in the form of single instruction multiple data sets (SIMD), requiring all solution processes to be completely homogeneous, such as multiple instances of the same model running simultaneously; and the various sub-models contained in complex products are usually different. For example, a car contains multiple sub-models such as engines, hydraulics, transmissions, brakes, etc. The calculation logic within the sub-models is different, and the parallel computing capabilities of the GPU cannot be brought into play, causing the GPU of the computing node to be idle during the simulation process and unable to provide computing power support, which in turn leads to an unbalanced load of computing resources within the computing node. Summary of the invention

[0008] In order to solve the above problems, the present invention provides a GPU parallel solution method for system simulation models.

[0009] The overall idea of ​​the technical solution of the present invention is to construct a solution method that combines two levels of parallelism, that is, between models and within models. First, according to the modeling description xml file of the complex product, the description information of the differential algebraic equations of all sub-models is extracted, and according to these description information, the function body and numerical solution code on the GPU side are dynamically generated and compiled, and the dynamic evaluation syntax tree of the expression on the CPU side is generated at the same time, and a list of solution tasks at the level of the equation group is established. Secondly, the main thread establishes a thread pool, traverses the list of tasks to be solved, and the main control algorithm dispatches the solution tasks of a specific ordinary differential equation group to the thread pool. The thread that obtains the solution task in the thread pool uses a multi-resolution method based on time domain decomposition to asynchronously push the specific ordinary differential equation solution calculation task to the GPU's stream processing unit; the GPU's multiple stream processors solve different sets of equations in parallel, and each stream processor solves a specific set of ordinary differential equations in parallel. When a pushed task is processed, a task completion event is generated, and a thread in the thread pool is awakened again to correct the GPU calculation result; if the corrected result does not meet the accuracy requirement, the asynchronous solution task is pushed to the GPU again, through multiple iterations until the solution accuracy requirement is met, and the result is returned to the main control algorithm to achieve the purpose of parallel and fast solution.

[0010] The present invention proposes a method for decomposing and solving simulation model tasks in parallel, comprising:

[0011] S1. Extract description information of ordinary differential equations, dynamically generate code for GPU, generate abstract syntax tree for CPU, and establish a task list for solving equations at the level of the equations;

[0012] S2. Establish a thread pool, traverse the solution task list, and dispatch the solution tasks of the ordinary differential equation group to the thread pool;

[0013] S3, the worker thread of the task is obtained from the thread pool, and the computing task of solving the ordinary differential equation group is asynchronously pushed to the stream processor of the GPU using a multi-resolution method based on time domain decomposition;

[0014] S4, GPU's multiple stream processors solve different sets of equations in parallel, and each stream processor solves a specific set of ordinary differential equations in parallel. When a stream processor completes the solution, a task completion event is generated;

[0015] S5, the task completion event of the GPU is captured by the signal processing thread, waking up the working thread in the thread pool;

[0016] S6. The awakened working thread obtains the GPU calculation result and corrects the result. If the corrected result meets the accuracy requirement, the result is returned; otherwise, the asynchronous solution task is pushed to the GPU again, and steps S4, S5, and S6 are repeated.

[0017] Furthermore, in the step S1, the equation group description information is extracted from the modeling description xml file of the complex product in units of sub-models.

[0018] Furthermore, in step S1, the step of dynamically generating code for the GPU includes:

[0019] Parse the name, equation group expression, and equation group coefficient parts of the input description information to generate the corresponding function body code; establish a function pointer lookup table based on the generated function body; merge the numerical solution GPU algorithm code with the generated code, call the GPU runtime compilation instruction, and generate GPU binary code.

[0020] Furthermore, in step S1, the step of generating an abstract syntax tree for the CPU includes:

[0021] Parse the name, equation group expression, and equation group coefficient parts of the input description information, add the state variables and variables to the symbol table, convert the logical operation on the right side of the equal sign into an expression instance, and register the symbols in the symbol table to the corresponding expression instance. Then, convert the expression instance into an abstract syntax tree (AST) and cache it in memory so that the corresponding equation group calculation process can be repeated. Each equation group corresponds to an AST syntax tree, and finally obtains the syntax tree list ExprList.

[0022] Furthermore, in step S1, the step of establishing a task list for solving a system of equations includes:

[0023] Parse the initial value of the state variable, integration time, global integration step, local integration step, and concurrent number of parallel solutions of the input description information to obtain the solution parameter list SolveParamList. Traverse ExprList and SolveParamList to establish the solution task list TaskList. Each task in TaskList consists of an ordered pair Task:={expr,solveParam}.

[0024] Furthermore, in step S2, the main thread first establishes a thread pool containing a certain number of threads, and then traverses the task list to solve each task Task i Dispatched to the thread pool. Preferably, the number of threads is 2 to 3 times the total number of CPU hardware cores.

[0025] Furthermore, in step S3, a worker thread Thread of a task acquisitionk Awakened to Task i Based on the information in the time domain, a multi-resolution method based on time domain decomposition is used to call the coarse-grained numerical solution algorithm to obtain the ordinary differential equation ODE i The 0th estimate of and assign it to Based on this estimated value As a parameter, push asynchronous parallel solution tasks to the GPU i .

[0026] Furthermore, in step S4, the GPU solves different sets of equations in parallel among multiple stream processors according to its internal execution logic, and solves a specific set of ordinary differential equations in parallel within each stream processor; at a certain moment, a stream processor SP is scheduled to is the initial value, and the fine-grained integration method is used to execute the internal parallel solution task; when the task is completed, a task completion event Event is generated i .

[0027] Furthermore, in step S5, the signal processing thread detects the task completion event Event i After that, wake up the worker thread Thread in the thread pool j .

[0028] Furthermore, in step S6, the awakened working thread Thread j Retrieving calculation results from GPU by The first element in is the initial value, and the coarse-grained numerical solution algorithm is called to obtain the ODE i The first estimate of And make the following corrections to the final result:

[0029]

[0030] if If the accuracy requirement is met, then return

[0031] Otherwise, the estimated value As parameters, push new asynchronous parallel solving tasks of ordinary differential equations to GPU;

[0032] Jump to step S4 and continue iterating until the accuracy requirement is met and the result is output.

[0033] Compared with the existing methods, the GPU-based heterogeneous ordinary differential equation solution method provided by the present invention has the following beneficial effects:

[0034] (1) The present invention provides a two-level parallel solution capability. The solution process of the ordinary differential equation system is mapped to different stream processors of the GPU, so that the solution calculation of the heterogeneous ordinary differential equation system can be parallelized between the stream processors; at the same time, the lightweight thread concurrent scheduling capability inside the GPU stream processor is utilized to complete the parallel solution of a single ordinary differential equation system in the form of a single instruction multiple data set (SIMD) in the stream processor. This not only solves the problem of parallel solution of heterogeneous ordinary differential equation systems, but also can obtain a better acceleration ratio in the solution of a single ordinary differential equation.

[0035] (2) The present invention has a wider application field. It is applicable to models containing large-scale heterogeneous ordinary differential equations and all computational fields involving the solution of ordinary differential equations. It is not limited to specific models and can speed up the solution speed, especially for those that require implicit numerical solution algorithms, which can effectively shorten the simulation time.

[0036] (3) The present invention can improve the utilization rate of hardware resources in the simulation process. By using the method provided by the present invention, the model simulation solution tasks can be flexibly divided as needed, and some computing tasks originally belonging to the CPU can be transferred to the GPU, so that the CPU and GPU loads of the computing node are relatively balanced, avoiding GPU idleness and waste. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] In order to more clearly illustrate the embodiments of the present disclosure or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present disclosure. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0038] Figure 1 A schematic diagram of a parallel solution process according to an embodiment of the present invention;

[0039] Figure 2 A schematic diagram for describing an ordinary differential equation according to an embodiment of the present invention;

[0040] Figure 3 A schematic diagram of dynamic generation of GPU code according to an embodiment of the present invention;

[0041] Figure 4 This is a schematic diagram of an example of ODE solving task parameters according to an embodiment of the present invention;

[0042] Figure 5 A schematic diagram of ODE solution scheduling according to an embodiment of the present invention;

[0043] Figure 6 This is a timing diagram of an ODE calculation process according to an embodiment of the present invention. DETAILED DESCRIPTION

[0044] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0045] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

[0046] The present invention provides a simulation model task decomposition and parallel solution method, such as Figure 1 As shown, the method comprises the following steps:

[0047] S1, extract the description information of the ordinary differential equation group of the simulation model, dynamically generate code for the GPU, generate an abstract syntax tree for the CPU, and establish a task list for solving the equation group level;

[0048] S2. Establish a thread pool, traverse the task list, and dispatch the task of solving the ordinary differential equation system to the thread pool;

[0049] S3, the worker thread of the task is obtained from the thread pool, and the ordinary differential equation solving calculation task is asynchronously pushed to the stream processing unit of the GPU using a multi-resolution method based on time domain decomposition;

[0050] S4, GPU's multiple stream processors solve different sets of equations in parallel, and each stream processor solves ordinary differential equations in parallel. When a stream processor completes the solution, a task completion event is generated;

[0051] S5, the task completion event of the GPU is captured by the signal processing thread, waking up the working thread in the thread pool;

[0052] S6. The awakened working thread obtains the GPU calculation result and corrects the result. If the corrected result meets the accuracy requirement, the result is returned; otherwise, the asynchronous solution task is pushed to the GPU again, and steps S4, S5, and S6 are repeated.

[0053] In step S1, the input information describing the ordinary differential equations of the simulation model is extracted. The equation group description information can be extracted from the product modeling description xml file in units of sub-models, including the name, equation group expression, equation group coefficients, initial values ​​of state variables, integration time, local integration step, and number of concurrent parallel solutions. In an implementation example, Figure 2As shown, it contains two different sets of ordinary differential equations, namely the Lorentz system and the modulated wave system. Among them, the Lorentz system contains all the elements of the input information, while the modulated wave system contains only one equation, which can be regarded as a special set of equations containing only one equation; the coefficient description part of the equation set is optional and is related to the specific characteristics of the ordinary differential equation set. In this implementation case, the description information of the modulated wave system does not include the coefficient part.

[0054] In the GPU dynamic code generation sub-step of step S1, as Figure 3 As shown, Figure 2 The function descriptions corresponding to the two equations in are parsed, and the corresponding function bodies are dynamically generated, named Lorenz and Modulated_wave. After the function bodies are generated, the function pointer array function_List is constructed. When the numerical solution algorithm needs to call a function body of a certain ordinary differential equation group, the parameter sequence number is used as the index to dynamically locate the specific function body pointer in the function_List array and execute the specific numerical integration step.

[0055] Except for the function body and function pointer array GPU code which are dynamically generated, the numerical solution algorithm can be a pre-compiled code stored in a file. During the compilation phase, the dynamically generated code and the numerical solution algorithm code are merged to generate GPU binary code.

[0056] In the present invention, the numerical solution algorithm of the GPU adopts a high-resolution method, such as a 4th-order Runge-Kutta or higher-order algorithm; the algorithm can be either an explicit solution method or an implicit solution method, and can be either a variable step size or a fixed step size. In this embodiment, the solution algorithm adopts a fixed-step-size 4th-order Runge-Kutta method.

[0057] In the sub-step of generating the CPU abstract syntax tree of step S1, according to Figure 2 According to the description in , parse the equation group name, equation group expression, and equation group coefficient part, add the state variables and variables to the symbol table, convert the logical operation on the right side of the equal sign into an expression instance, and register the symbols in the symbol table to the corresponding expression instance. Then, convert the expression instance into an abstract syntax tree (AST) and cache it in memory so that the corresponding equation group calculation process can be repeated. Each equation group corresponds to an AST syntax tree, and finally obtains a syntax tree list ExprList containing 2 elements. Among them, ExprList[0] corresponds to the Lorentz system, and ExprList[1] corresponds to the modulated wave system.

[0058] Specifically, in order to support runtime evaluation of expressions for ordinary differential equations, the following four types of operators are supported in the abstract syntax tree:

[0059] Arithmetic operations: addition, subtraction, multiplication, division, modulus, and power;

[0060] General operations: absolute value, logarithm, exponent, maximum value, minimum value, integer, square root, threshold;

[0061] Trigonometric operations: sine, cosine, tangent, cotangent;

[0062] Boolean operations: logical AND, logical OR, logical NOT, logical XOR, bitwise AND, bitwise OR.

[0063] These operations are pre-compiled into binary code and can be called directly during the post-order traversal of the abstract syntax tree.

[0064] Dynamic expression evaluation can also be achieved by calling existing open source expression evaluation libraries, such as ExprTk.

[0065] In the sub-step of establishing the task list for solving the system of equations in step S1, according to Figure 2 , parse the input description information, including the initial value of the state variable, integration time, global integration step, local integration step, and the number of concurrent parallel solutions, to obtain the solution parameter list SolveParamList. Traverse ExprList and SolveParamList to establish a solution task list TaskList containing two elements. Each task in TaskList consists of an ordered pair Task:={expr,solveParam}. In this implementation example, the SolveParamList array contains two elements of the SolveParam type. The specific field value in each element is as follows: Figure 4 shown.

[0066] In step S2, the main thread first establishes a thread pool containing 2 threads, and then traverses the ordinary differential equation solving task list, and sets each task Task i Dispatched to the thread pool. The number of threads is usually 2 to 3 times the total number of cores of the current CPU hardware. All threads are initially in a waiting state. If the number of input equations is small, the number of equations actually required to be solved is used as the configuration of the thread pool.

[0067] In step S3, if Figure 5 As shown, a thread that gets the task k ,k∈{0,1} is awakened with Task i Based on the information in ,i∈{0,1+, the multi-resolution method based on time domain decomposition is used to call the coarse-grained numerical solution algorithm to obtain the ordinary differential equation system ODE i The 0th estimate of and assign it to Based on this estimated value As a parameter, push asynchronous parallel solution tasks to the GPU.

[0068] In this implementation example, if the current thread is processing the task of solving the Lorentz system, it corresponds to Task0. First, according to the specific values ​​of the fields in SolveParamList[0], we can get the starting integration time t0=0.0, the global integration step ΔT0=0.2, and the local integration step δT=0.01 of Task0. Select the coarse-grained numerical solution algorithm, the second-order Runge-Kutta method. With δT as the step size, perform the cumulative integration from time t0 to time t0+ΔT0 to obtain a set of initial value estimation vectors with a length of 20 by As the initial value, push asynchronous computing tasks to the GPU.

[0069] When pushing computing tasks to the GPU, you need to specify the number of lightweight threads for parallel solving and the GPU integration step. In this implementation example, N=64 in SolveParamList[0]. However, since 64 cannot divide 20, it is necessary to make a correction and take N=20×(64÷20)=60, and the corresponding δt=δT÷(64÷20). That is, the GPU integration step is 1 / 3 of the CPU integration step. When actually setting ΔT, δT and N, the following joint conditions should be met as much as possible:

[0070]

[0071] That is, N is divisible by both (ΔT÷δT) and 32. The first condition maintains the integer multiple relationship between the coarse-grained and fine-grained integration steps to avoid interpolation or extrapolation when the result is corrected (see step S6 for details); the second condition makes the GPU's concurrent scheduling unit (wrap) run at full capacity as much as possible to reduce idle threads. If these two conditions cannot be met at the same time, condition one should be met first, and N should be kept as close to an integer multiple of 32 as possible. In order to speed up the convergence speed and shield the GPU hardware scheduling time, N should be as large as possible.

[0072] In step S4, Figure 5 , 6 As shown in FIG. 1 , multiple stream processors of the GPU solve different sets of equations in parallel, and each stream processor solves a specific set of ordinary differential equations in parallel. The GPU schedules a stream processor SP at a certain moment according to its internal execution logic, and uses a fine-grained algorithm (the 4th-order Runge-Kutta method, see step S1 for details) to vector is the initial value, the number of concurrent connections is N i ,i∈{0,1}, with an integration step of δt i ,i∈{0,1}, in δT i,i∈{0,1} time period, execute the parallel solution task of a specific set of equations; when the task is completed, a task completion event Event is generated i ,i∈{0,1}.

[0073] In step S5, Figure 5 , 6 As shown, the signal processing thread detects the task completion event Event i ,i∈{0,1}, wake up the thread Thread in the thread pool j ,j∈{0,1}.

[0074] In step S6, Figure 5 , 6 As shown, the awakened thread Thread j ,j∈{0,1} retrieves the calculation result vector from the GPU by The first element in is the initial value, and the coarse-grained numerical solution algorithm is called again to obtain the ordinary differential equation system ODE i , the first estimated value of i∈{0,1} And make the following corrections to the final result:

[0075]

[0076] Determine whether the corrected result meets the accuracy requirements:

[0077]

[0078] Where ∈ is the error limit.

[0079] if If the accuracy requirement is met, the ODE system is returned. i ,i∈{0,1} in the current global integration step ΔT i , the result vector of i∈{0,1}

[0080] Otherwise, the estimated value As a parameter, push a new asynchronous parallel solution task to the GPU, the current thread Thread j ,j∈{0,1} enters the waiting state. The algorithm operation logic goes to step S4 and continues to iterate until the accuracy requirement is met and the result is output.

[0081] For the same ODE system i ,i∈{0,1}, after the algorithm runs to step S4, all vectors with superscripts in steps S4 to S6 have their corresponding superscripts incremented during calculation or solution.

[0082] Those skilled in the art will appreciate that the foregoing embodiments are merely intended to illustrate the technical solutions of the present invention rather than to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art will appreciate that the technical solutions described in the foregoing embodiments may be modified or some or all of the technical features may be replaced by equivalents. Such modifications or replacements do not deviate the essence of the corresponding technical solutions from the scope defined by the claims of the present invention.

Claims

1. A method for task decomposition and parallel solution of a simulation model, characterized in that: Construct a solution method that combines two levels of parallelism between simulation models and within simulation models, including: S1, extract the description information of the ordinary differential equation group of the simulation model, dynamically generate code for the GPU, generate an abstract syntax tree for the CPU, and establish a task list for solving the equation group level; S2, establish a thread pool, traverse the solution task list, and dispatch ordinary differential equation solution tasks to the thread pool; S3, obtain the working thread of the equation solving task from the thread pool, and use the multi-resolution method based on time domain decomposition to push the asynchronous parallel solution requirements to the GPU; S4, multiple stream processors of the GPU solve different sets of ordinary differential equations in parallel, and each stream processor solves the set of ordinary differential equations in parallel. When the stream processor completes the solution, a task completion event is generated; S5, the task completion event of the GPU is captured by the signal processing thread, waking up the working thread in the thread pool; S6. The awakened working thread obtains the GPU calculation result and corrects the result. If the corrected result meets the accuracy requirement, the result is returned; otherwise, the asynchronous solution task is pushed to the GPU again, and steps S4, S5, and S6 are repeated. In step S3, a worker thread Thread that obtains the task k Awakened to Task i Based on the information in, the coarse-grained numerical solution algorithm is called to obtain the ordinary differential equation ODE i The 0th estimate of and assign it to Based on this estimated value As a parameter, push asynchronous parallel solution tasks to the GPU i ; In step S4, at a certain moment, a stream processor SP is scheduled to is the initial value, and the fine-grained integration method is used to execute the internal parallel solution task; when the task is completed, a task completion event Event is generated i ; In step S5, the signal processing thread detects the task completion event Event i After that, wake up the worker thread Thread in the thread pool j ; In step S6, the awakened worker thread Thread j Retrieving calculation results from GPU by The first element in is the initial value, and the coarse-grained numerical solution algorithm is called to obtain the ODE i The first estimate of And make the following corrections to the final result: if If the accuracy requirement is met, then return Otherwise, the estimated value As parameters, push new asynchronous parallel solving tasks of ordinary differential equations to GPU; Jump to step S4 and continue iterating until the accuracy requirement is met and the result is output.

2. The method according to claim 1, characterized in that In step S1, the description information of the ordinary differential equation group of the simulation model is extracted in units of sub-models according to the product modeling description xml file, including the name, equation group, equation group coefficients, initial values ​​of state variables, integration time, local integration step, and number of concurrent parallel solutions.

3. The generation method according to claim 1, characterized in that: In step S1, the step of dynamically generating code for the GPU includes: Parsing the name, equation group expression, and equation group coefficient part of the description information of the ordinary differential equation group of the simulation model to generate corresponding function body code; According to the generated function body, a function pointer lookup table is established; the numerical solution GPU algorithm code is merged with the generated code, and the GPU runtime compilation instruction is called to generate GPU binary code.

4. The generation method according to claim 1, characterized in that: In step S1, the step of generating an abstract syntax tree for the CPU includes: Parse the name, equation group expression, and equation group coefficient part of the description information of the ordinary differential equation group of the simulation model, add the state variables and variables to the symbol table, convert the logical operation on the right side of the equal sign into an expression instance, and register the symbols in the symbol table to the corresponding expression instance; Then, the expression instance is converted into an abstract syntax tree (AST) and cached in the memory so that the corresponding equation group calculation process can be repeated; each equation group corresponds to an AST syntax tree, and finally a syntax tree list ExprList is obtained.

5. The generation method according to claim 1, characterized in that: In step S1, the steps of establishing a task list for solving a system of equations include: Parse the initial value of the state variable, integration time, global integration step, local integration step, and concurrent number of parallel solutions of the input description information to obtain the solution parameter list SolveParamList; Traverse ExprList and SolveParamList to establish a solution task list TaskList, where each task in TaskList consists of an ordered pair Task:={expr,solveParam}.

6. The generation method according to claim 1, characterized in that: In step S2, the main thread first establishes a thread pool containing a certain number of threads, and then traverses the ordinary differential equation solving task list, and assigns each task Task i Dispatched to the thread pool.

7. The generation method according to claim 6, characterized in that: The number of threads is 2 to 3 times the total number of CPU hardware cores.

Images