A Cluster-CPU-GPU Heterogeneous Parallel Method for Solving Responses of Large-Scale Anisotropic Structures

Abstract

This invention discloses a cluster-CPU-GPU heterogeneous parallel method for solving the structural response of large-scale anisotropic materials, belonging to the field of high-performance computing simulation technology. First, macroscopic and microscopic structures are established. The macroscopic model is divided into finite element meshes, and the computational domain is divided using a domain decomposition algorithm. OpenMP is used to accelerate the calculation of the stiffness matrix and displacement of macroscopic model elements. Next, the model is reconstructed at Gaussian points of the macroscopic model, and MPI technology is used for block division. Each process uses multi-threading to calculate the displacement of the microscopic model, and OpenMP is used to accelerate the calculation of the stiffness matrix of the microscopic model elements, achieving cluster MPI parallel computing. In post-processing, the multi-core advantage of the GPU is used to calculate the stress and strain of the macroscopic model for visualization, achieving cross-scale heterogeneous parallel simulation. This method, combined with relevant algorithms, can efficiently utilize resources to solve simulation problems and has broad application prospects.

Description

Technical Field

[0001] This invention belongs to the field of high-performance finite element simulation technology, specifically a cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures. Technical Background

[0002] Anisotropic materials exhibit physical properties that differ in various directions. They are mainly classified into crystals, composite materials, and fiber materials, and research on these materials requires consideration of their unique properties. Anisotropic materials currently have wide applications in fields such as medicine and aerospace.

[0003] For a long time, research on the macroscopic properties of anisotropic materials has mainly relied on experimental trial and error. However, this method is not only costly and prone to randomness, but also has a long experimental cycle, resulting in a significant waste of human and material resources. Currently, with the rapid development of computer technology, using high-performance computers to predict the properties of anisotropic materials has become a research hotspot. Compared to traditional experimental methods, numerical simulation calculations using computers not only allow for easier design of the material's microstructure, but are also simpler and more repeatable, and can more quickly yield the material's mechanical properties, facilitating further optimization and improvement. However, to obtain more accurate mechanical properties, larger-scale models are needed. But large-scale solutions often compromise the accuracy and correctness of the results. Currently, while the increase in computer clock speed is relatively small, the increase in the number of cores is significant. Fully utilizing computing resources is crucial for achieving higher accuracy and efficiency in solving these problems.

[0004] Therefore, a cluster-CPU-GPU heterogeneous parallel method based on domain decomposition is established to solve the response of large-scale anisotropic structures. This method applies the mature finite element numerical computation theory and fully utilizes the advantages of different parallel methods to accelerate computation. It can more effectively utilize computing resources to achieve rapid solutions for large-scale models. Furthermore, the cross-scale approach can better reflect the microstructure of anisotropic materials. This method can not only be applied to aerospace and other fields but also provide more accurate solutions in other areas, promoting the optimization of mechanical properties in material and structural designs. Summary of the Invention

[0005] This invention discloses a cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures, belonging to the field of finite element high-performance computing simulation technology. This method utilizes DDM to achieve domain decomposition of large-scale models and combines the parallel advantages of clusters, CPUs and GPUs to provide an efficient and high-precision solution method. By combining cross-scale analysis of the microstructure of anisotropic materials, it achieves rapid solution for large-scale anisotropic materials.

[0006] The solution provided by this invention to the above-mentioned technical problems is: a cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures, comprising the following steps:

[0007] S1. Establish macroscopic and microscopic geometric models of arbitrary scale, spatially discretize the macroscopic and microscopic geometric models to be analyzed, divide them into grids, and generate the required grid data.

[0008] S2. Set boundary conditions, constraints, and material parameters for the macroscopic model, and use CPU multithreading technology (OpenMP) to calculate the element stiffness matrix of the macroscopic model in parallel.

[0009] S3. The model is partitioned using a domain decomposition algorithm based on set values. This divides the solution mesh model into multiple computational sub-regions, ensuring that each partition carries all its parameters. The algorithm iterates through each partition based on the partition information to solve for the displacement field of the complete model. During the solution process, the open-source library STXXL is used to store data from memory on disk, reducing memory usage and enabling faster data retrieval from the disk.

[0010] S4. Based on the calculated displacement results, calculate the strain and stress at each Gaussian point of each element in the macroscopic model;

[0011] S5. Using the multi-process parallel technology (MPI) of the cluster, the units of the macro model are evenly distributed to each process, and each process calculates the displacement of the micro model at each Gaussian point of its own macro model unit.

[0012] S6. Reconstruct the model at each Gaussian point of each unit in the macroscopic model under each process, use the microscopic model as the microstructure of the macroscopic model, and use the strain at the Gaussian point to calculate the displacement boundary conditions of the microstructure at that point.

[0013] S7. In each process, CPU multithreading technology (OpenMP) is used to calculate the element stiffness matrix of each element in the microstructure in parallel, assemble the overall stiffness matrix and solve the displacement field, and calculate the strain and stress at each Gaussian point of each element in the micro model based on the displacement results.

[0014] S8. In each process, the average stress of the entire micro-model is solved in parallel using CPU multi-threading technology (OpenMP) based on the stress at each Gaussian point on the micro-model. This average stress is then used as the average stress at a single Gaussian point of a single element in the macro-model. Finally, the internal forces of the macro-model element are calculated based on the average stress at the four Gaussian points of the single element in the macro-model.

[0015] S9. Each process completes the internal force calculation of all units within the process. Process 0 collects the calculation results of the internal forces of all units in the macroscopic model to obtain the internal force vector of the macroscopic model, and calculates the displacement of the macroscopic model based on the internal force vector.

[0016] S10. Iterate the displacement judgment residuals obtained by using the displacement calculated from the displacement boundary of the macroscopic model and the internal force vector, and obtain the displacement result of the macroscopic model after the judgment converges.

[0017] S11. By using GPU to calculate the stress and strain of the macroscopic model based on the displacement results, the post-processing of the calculation is quickly realized, and a cross-scale heterogeneous parallel algorithm for anisotropic material structures is achieved.

[0018] The specific implementation method of step S1 is as follows: Establish a three-dimensional tetrahedral basic mesh element with a reasonable density based on the complexity at both the macroscopic and microscopic scales. The model saves the number of model elements, the nodes contained in each element, the number of nodes, and the three-dimensional coordinates of all nodes.

[0019] The specific implementation method of step S2 is as follows:

[0020] S21. Set material parameters, boundary conditions and constraints;

[0021] Material parameters are set for both the macroscopic and microscopic models. The macroscopic model is treated as an isotropic material, with Poisson's ratio and elastic modulus set according to requirements. The microscopic model is an anisotropic material, requiring the setting of elastic modulus and Poisson's ratio in each direction. For the macroscopic model, boundary conditions and constraints also need to be set, determining the locations of applied displacement boundaries and fixed constraints.

[0022] S22. Calculation of element stiffness matrix;

[0023] Based on the set macroscopic model material parameters, the stiffness matrix of the macroscopic model element is calculated. First, the shape function of the three-dimensional tetrahedral element is calculated. For the tetrahedral element used, its volume coordinates are defined as follows:

[0024]

[0025]

[0026] In the parametric coordinate system, the shape function is equal to the corresponding volume coordinate, and the four shape functions corresponding to the four nodes are:

[0027]

[0028] The Jacobian matrix is ​​obtained by taking the partial derivatives of the shape function with respect to each direction in the local coordinate system and then multiplying it by the nodal coordinates. The method for solving the Jacobian matrix is ​​as follows;

[0029]

[0030] The inverse of the Jacobi matrix is ​​multiplied by the partial derivative of the form function to obtain the B matrix;

[0031]

[0032] The D matrix form of an isotropic three-dimensional tetrahedral element is as follows;

[0033]

[0034] Where E and v are the elastic modulus and Poisson's ratio of the macroscopic model, respectively.

[0035] The expression for the element stiffness matrix is ​​obtained through the derivation of the energy principle. Since the expression for the stiffness matrix is ​​an integral operation, the element stiffness matrix is ​​calculated by introducing numerical integration. The Hammer numerical integration results are as follows.

[0036]

[0037] S23, OpenMP parallel computing;

[0038] Since the element stiffness matrix calculation for each element is independent and only occurs within a single loop with few variables, OpenMP multithreading technology is used for parallel computation of the element stiffness matrix. The elements of the macroscopic model are divided into blocks, and each thread calculates the element stiffness matrix for its assigned element.

[0039] This completes the setting of boundary conditions, constraints, and material parameters for the macroscopic model, and uses CPU multithreading technology (OpenMP) to calculate the element stiffness matrix of the macroscopic model in parallel.

[0040] The specific implementation method of step S3 is as follows:

[0041] S31, Model Domain Decomposition;

[0042] The number of regions N is determined by the size of the calculated model. The relationship between the number of regions and the partitioning method is usually that when N is less than 8, a multi-level recursive equal partitioning method is used, and when N is greater than or equal to 8, a multi-level K-way partitioning method is used.

[0043] After partitioning, the nodes and elements of each partition will be renumbered, and the node information, element information, and element information at the boundary of each partition before and after renumbering will be saved. At this point, the original finite element model is divided into N regions Ω0 to Ω0. N-1

[0044] S32, Calculation of partition stiffness matrix;

[0045] After completing the regional division of the model, the regional stiffness matrix of each region is calculated by integrating the regional stiffness matrix according to the regional situation. Based on the regional element information in S31, the element stiffness matrix of each region is extracted, and the nodal degrees of freedom are encoded according to the re-compiled nodal information of the regional region to construct the regional stiffness matrix.

[0046] S33. Iteratively solve the displacement field;

[0047] Firstly, based on the global Schwarz solution iteration in the domain decomposition algorithm, we have:

[0048]

[0049] Where ω is an arbitrary function mapping from the field Ω to R, u n This is an approximate solution to the Poisson problem. This formula is used to simultaneously solve subproblems in all subdomains of the domain algorithm, and its convergence is proven using the maximum principle. Based on this, the residuals are introduced as follows;

[0050] r n =f+Δ(u) n )

[0051] This allows us to define the Restricted Additive Schwarz Algorithm (RAS) on a continuous level.

[0052] In finite element elastic systems, for quasi-static solutions, according to Hooke's law and the principle of minimum potential energy, we have:

[0053] KU = F

[0054] Where K is the global stiffness matrix, F is the force matrix, and U is the displacement vector result. Based on the domain decomposition algorithm, the above linear system can be transformed into the following block matrix form:

[0055]

[0056] Since the RAS algorithm is defined based on the continuous level, while the finite element algorithm is based on the discrete level, it is necessary to apply the RAS algorithm to the discrete level from the perspective of linear algebra. Introducing the block Jacobi algorithm, and letting D be a diagonal matrix of K, we have:

[0057] DU n+1 =DU n +(F-KU n )

[0058] Define the residual of the discrete level as r n =F-KU n By dividing the Jacobi algorithm into blocks, the block Jacobi algorithm can be defined as follows:

[0059]

[0060] Among them, K ii Here is the overall stiffness matrix for each partition. Let R be the displacement vector of the i-th partition at the n-th iteration step. Then, we introduce R... i As N to N i Restricted operators, As N to N i The transpose operator then has At this point, the block Jacobi algorithm can be written in the following form:

[0061]

[0062] Based on the above formula and combined with the RAS algorithm applicable to any number of partitions at the discrete level, we have:

[0063]

[0064] For the block Jacobi algorithm, when expanded into a matrix with 2 partitions, we have:

[0065]

[0066] in, This represents the displacement of the first partition in the nth iteration step. The non-overlapping parts, This represents the displacement of the first partition in the nth iteration step. The overlapping part. Represents the residual of the first partition in the nth iteration step. Non-overlapping parts, This represents the overlapping portion of the residuals of the first and second partitions in the nth iteration step.

[0067] The simplified linear algebraic form of the above equation is:

[0068]

[0069] in, This is the local correction for solving the discrete level in the nth step, and we have

[0070] Solve the linear equations for each partition, and then multiply the results by half for the overlapping parts and sum them. This form applies to all partition counts. For regions with more than two partitions, the above formula can be rewritten by simply including the partition information of the overlapping nodes in the constant term matrix. This completes the displacement field solution using the RAS algorithm based on domain decomposition.

[0071] S34, Memory data storage;

[0072] For the displacement field solution process described above, a model with tens of millions of elements requires storing a large amount of information in memory during actual calculations. This includes model information such as element numbers, node coordinates, and the B-matrix and material matrix of each element. This information consumes a significant amount of memory during computation. Using the open-source library STXXL allows the memory data to be stored on disk, reducing memory usage and enabling faster access to disk data.

[0073] This open-source library uses the currently popular Parallel Disk Model (PDM). In this model, file I / O is explicitly handled by the program. It divides M bytes of data in memory into D blocks and transfers them between D hard disks to amortize latency. Typically, external memory algorithms assume they read input from external storage (hard disk) and write output to it. STXXL, however, uses pipelining to provide a new interface for external memory algorithms, allowing them to provide output as a data stream directly to the algorithm using the output, instead of writing it to memory first. Furthermore, STXXL uses synchronous I / O calls to run asynchronously in separate threads; for each file, there is a read request queue, a write request queue, and a thread. The main thread sends requests to the file queues, and the threads associated with the files execute the requests in FIFO order. This method allows data in memory to be stored on hard disk, reducing memory usage.

[0074] This completes the partitioning of the model using the domain decomposition algorithm based on set values. The solution mesh model is divided into multiple computational sub-regions, with each region carrying all its parameters. The algorithm iterates through each region based on the partition information to solve for the displacement field of the complete model. During the solution process, the open-source library STXXL is used to store data from memory on disk, reducing memory usage and enabling faster disk data retrieval.

[0075] The specific implementation method of step S4 is as follows:

[0076] S41. Solve for the stress and strain at the Gauss point;

[0077] The stress and strain at Gaussian points are calculated based on the three-dimensional tetrahedral element used. This three-dimensional tetrahedral element has four Gaussian points, and the stress and strain in each direction at each Gaussian point need to be calculated separately. The calculation method is as follows.

[0078] ε=BU

[0079] Where ε is the strain matrix at a single Gaussian point, recording the strain and shear strain in three directions, B is the B matrix at that Gaussian point, and U is the displacement matrix of each node on the element, sorted according to the element connection order. The stress matrix is ​​then calculated after obtaining the strain matrix.

[0080] σ=Dε

[0081] Where σ is the stress matrix at a single Gaussian point, recording the stress and shear stress in three directions, and D is the material matrix of the element.

[0082] This completes the calculation of strain and stress at each Gaussian point of each element in the macroscopic model based on the calculated displacement results;

[0083] The specific implementation method of step S5 is as follows: Since it is necessary to calculate the displacement of the micro-model at each Gaussian point, and more calculations are involved than a single loop, MPI multi-process technology is used for parallel processing. First, the macro-units are divided into blocks, and the stress, strain, and displacement of each Gaussian point on each block are distributed to each process using the MPI communication protocol.

[0084] The specific implementation method of step S6 is as follows: Reconstruct the model at each Gaussian point on each element of the macroscopic model, and calculate the displacement boundary of the microscopic model using the strain at that Gaussian point. Since the finite element method first calculates the stress and strain at the Gaussian point, and then uses interpolation to calculate the stress and strain at other locations of the element, the microscopic model is treated as a point on the macroscopic model. The new displacement and stress / strain of the microscopic model at this point are calculated, and the result is then treated as the new stress at that Gaussian point, thus achieving the purpose of cross-scale calculation.

[0085] The specific implementation method of step S7 is as follows:

[0086] S71, Calculation of boundary conditions for microscopic models;

[0087] After distributing the macroscopic model to each process according to the previous step, each process needs to perform a complete calculation of the microscopic model at each Gaussian point within that unit. At this point, each process contains all the model information of the microscopic model, as well as the strain at each Gaussian point. First, the boundary nodes and coordinates of the microscopic model are extracted, and the displacement boundary conditions of the microscopic model are obtained by multiplying the strain at the Gaussian point by the coordinates.

[0088] Δd=εd

[0089] Where Δd represents the displacement of each boundary node in the microscopic model, and d represents the coordinates of each boundary node in the microscopic model. This yields the boundary conditions of the microscopic model.

[0090] S72, Calculation of stiffness matrix of micro-model element;

[0091] The same finite element theory from S22 is used to calculate the shape functions, etc., in the stiffness matrix calculation of the microscopic model elements. Since the macroscopic model is treated as an isotropic material, while the microscopic model is treated as an anisotropic material, the calculation method for the material matrix D differs between the microscopic and macroscopic models. Therefore, the material matrix D for the microscopic model is calculated as follows:

[0092]

[0093] The 36 different parameters in the matrix are determined based on the different forms of anisotropic materials used. The stiffness matrix of each element is calculated using the same method as in S2.

[0094] S73, OpenMP parallel computing;

[0095] Accelerating the calculation of element stiffness matrix using OpenMP based on the same method in S23, and achieving simultaneous scheduling of processes and threads.

[0096] S74, Assembly of the overall stiffness matrix;

[0097] After completing the OpenMP calculation of the element stiffness matrix, the element stiffness matrix is ​​assembled into the global stiffness matrix according to the degree of freedom index of each node in each element to complete the assembly of the element stiffness matrix and obtain the global stiffness matrix K.

[0098] S75, Solving for displacement, stress, and strain;

[0099] Using the same method as in S41, the displacement field in the microscopic model and the stress and strain values ​​at each Gaussian point on the microscopic model are solved.

[0100] This completes the solution of displacement, stress, and strain of the micro-model at a single Gaussian point. It realizes the parallel calculation of the element stiffness matrix of each element on the microstructure using CPU multi-threading technology (OpenMP) in each process, assembles the overall stiffness matrix, solves the displacement field, and calculates the strain and stress at each Gaussian point of each element of the micro-model based on the displacement results.

[0101] The specific implementation method of step S8 is as follows:

[0102] S81, Calculation of mean stress in microscopic model;

[0103] Based on the stresses at all Gaussian points on a single microscopic model calculated in S75, the average stress on the microscopic model at that Gaussian point is calculated. Using the homogenization theory of simulating macroscopic mechanical properties through a composite material mesoscopic unit cell model, the average stress is calculated as follows:

[0104]

[0105] in, The average stress of the model is σ. ij Let be the stress at each Gaussian point.

[0106] S82 and OpenMP parallelism;

[0107] Since the mean stress calculations for each element in the microscopic model are independent, the S23 method is used to accelerate the calculation of the mean stress, and OpenMP is used to enable multithreading to accelerate the loop.

[0108] S83, Calculation of internal forces in microscopic model elements;

[0109] After obtaining the average stress of a single micro-model, return to the macro-model. That is, the stress at a single Gaussian point in the macro-model is equivalent to the average stress of the micro-model at that Gaussian point. Calculate the internal forces of the macro-model based on this average stress.

[0110]

[0111] Where B is the macroscopic element matrix B, S is the stress matrix, and J is the Jacobian matrix. The internal forces within a single element of the macroscopic model are then calculated.

[0112] This completes the process of using CPU multi-threading technology (OpenMP) to solve the average stress of the entire micro-model in parallel under each process, based on the stress at each Gaussian point on the micro-model, as the average stress on a single Gaussian point of a single element in the macro-model, and then calculating the internal force of the macro-model element based on the near-average stress at the four Gaussian points of the single element in the macro-model.

[0113] The specific implementation method of step S9 is as follows:

[0114] S91, Collect internal force vectors;

[0115] After each process completes the internal force calculation for all units within its process, it sends the complete internal force vectors of each unit to process 0 via the MPI communication protocol. Once the communication is complete, i.e., process 0 has finished collecting the data, the complete internal force vector of the macroscopic model is obtained.

[0116] S92. Calculate the displacement of the macroscopic model;

[0117] Obtain the internal force vector F in Then, the new displacement field U is calculated using the following formula. new .

[0118] KU = F

[0119] This completes the internal force calculation of all units within each process. Process 0 collects the calculation results of the internal forces of all units in the macroscopic model to obtain the internal force vector of the macroscopic model, and calculates the displacement of the macroscopic model based on the internal force vector.

[0120] The specific implementation method of step S10 is as follows: The external force vector F obtained by setting the boundary conditions to a larger number in S2 is... out The internal force vector F obtained from S9 in Adding them together gives F r , use F r Solving the stiffness matrix yields the new displacement U. n This converges with the previous judgment.

[0121] If convergence is calculated, F r Calculated U n This is the final displacement. If r does not meet the convergence condition, then the displacement U obtained from the external force calculation in the previous step will be... n As a new step, displacement U n+1 , with U n+1 Set these as boundary conditions for the microscopic model in S6, and repeat the above steps until convergence.

[0122] The specific implementation method of step S11 is as follows: after the model displacement is calculated in S10, the stress and strain of each element are calculated using the method in S41, but the data is copied to the GPU for calculation.

[0123] ε=BU

[0124] σ=Dε

[0125] Since this section involves only a large number of matrix operations, GPU parallelism, specifically CUDA, is used for post-processing in the finite element method. First, the displacement matrices, B matrices, and material matrices of all elements are copied from CPU memory to the GPU in one go. The number of blocks and threads are then set on the GPU. For strain calculation, each thread calculates one row of the B matrix and one column of the U matrix, achieving parallel calculation of strain results for all elements across all threads. For stress calculation, each thread calculates one row of the D matrix and one column of the strain matrix, achieving parallel calculation of stress results for all elements across all threads. This yields the stress and strain at each Gaussian point of all elements in the macroscopic model.

[0126] Since stress and strain are exact solutions obtained at the integration points, and the strain and stress at these integration points are located inside the element, linear interpolation is performed using nodal values ​​and shape functions to obtain the results at each nodal to improve the global strain and stress results. The calculation formula is as follows:

[0127]

[0128] When multiple elements share a single node, the values ​​at the node are interpolated from the integration points of several adjacent elements, resulting in multiple stress and strain values ​​at the node. In such cases, grinding is required. For example, the strain and stress at a common node are shared by four elements on the outer surface.

[0129] The calculation method is as follows:

[0130]

[0131] This completes the implementation of CPU-GPU heterogeneous parallelism based on clusters, enabling efficient simulation of anisotropic materials by effectively utilizing computing resources.

[0132] The beneficial effects of this invention are as follows: This invention achieves an efficient algorithm for anisotropic materials based on the advantages of different parallel methods. First, a large-scale macroscopic and microscopic structure is established and the mesh is discretized. Then, DDM technology is used to decompose the large-scale discrete model into regions. OpenMP is used to accelerate the parallel calculation of the stiffness matrix of each element. After iteratively solving the displacement of the macroscopic model, the stress and strain at each Gaussian point of the macroscopic model are calculated. MPI technology is used to evenly distribute the macroscopic model elements to each process. The macroscopic model strain is used to calculate the boundary conditions of the microscopic model, and then the displacement of the microscopic model is calculated. The average stress of the microscopic model is further calculated. The average stress of the microscopic model is used as the stress at a single Gaussian point of the macroscopic model. Then, the internal forces of the macroscopic model are calculated, and MPI technology is used to reduce all internal forces in each process. The new displacement is calculated using this internal force, and the residual is calculated with the original displacement. The iteration is performed until the iteration converges to obtain the final displacement. After completing the displacement calculation, CUDA on the GPU is used for post-processing in parallel. This paper proposes a large-scale heterogeneous parallel simulation method based on domain decomposition and combined with multi-process, multi-thread, and GPU parallelism. It uses the "divide and conquer" approach to solve the problems of insufficient computing resources and excessive computing time, thereby realizing the calculation of the response of large-scale anisotropic material structures. Attached Figure Description

[0133] Figure 1 This is a flowchart of the method for solving large-scale cluster-CPU-GPU heterogeneous parallelism in this invention;

[0134] Figure 2 Diagrams showing macroscopic and microscopic structures;

[0135] Figure 3 This is a schematic diagram of a tetrahedral unit cell;

[0136] Figure 4 This is a schematic diagram of the OpenMP process;

[0137] Figure 5A schematic diagram of K-way partitioning;

[0138] Figure 6 This is a diagram showing the partitioning results;

[0139] Figure 7 A schematic diagram of MPI collection and distribution;

[0140] Figure 8 This is a contour plot of the displacement results;

[0141] Figure 9 This is a schematic diagram of shared nodes;

[0142] Figure 10 This is a stress-strain contour plot; Detailed Implementation

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

[0144] like Figure 1 As shown, a cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures is characterized by the following steps:

[0145] S1. Establish macroscopic and microscopic geometric models of arbitrary scale, spatially discretize the macroscopic and microscopic geometric models to be analyzed, divide them into grids, and generate the required grid data.

[0146] The specific implementation steps are as follows: Establish a three-dimensional tetrahedral basic mesh element with a reasonable density based on the complexity at both the macroscopic and microscopic scales. The macroscopic and microscopic models are as follows: Figure 2 As shown. The model stores the number of model elements, the nodes contained in each element, the number of nodes, and the 3D coordinates of all nodes.

[0147] S2. Set boundary conditions, constraints, and material parameters for the macroscopic model, and use CPU multithreading technology (OpenMP) to calculate the element stiffness matrix of the macroscopic model in parallel.

[0148] S21. Set material parameters, boundary conditions and constraints;

[0149] Material parameters are set for both the macroscopic and microscopic models. The macroscopic model is treated as an isotropic material, with Poisson's ratio and elastic modulus set according to requirements. The microscopic model is an anisotropic material, requiring the setting of elastic modulus and Poisson's ratio in each direction. For the macroscopic model, boundary conditions and constraints also need to be set, determining the locations of applied displacement boundaries and fixed constraints.

[0150] S22. Calculation of element stiffness matrix;

[0151] Based on the pre-set material parameters of the macroscopic model, the stiffness matrix of the macroscopic model elements is calculated. First, the shape function of the three-dimensional tetrahedral element is calculated. For the tetrahedral element used, such as... Figure 3 As shown, its volume coordinates are defined as follows:

[0152]

[0153] In the parametric coordinate system, the shape function is equal to the corresponding volume coordinate, and the four shape functions corresponding to the four nodes are:

[0154]

[0155] The Jacobian matrix is ​​obtained by taking the partial derivatives of the shape function with respect to each direction in the local coordinate system and then multiplying it by the nodal coordinates. The method for solving the Jacobian matrix is ​​as follows;

[0156]

[0157] The inverse of the Jacobi matrix is ​​multiplied by the partial derivative of the form function to obtain the B matrix;

[0158]

[0159] The D matrix form of an isotropic three-dimensional tetrahedral element is as follows;

[0160]

[0161] Where E and v are the elastic modulus and Poisson's ratio of the macroscopic model, respectively.

[0162] The expression for the element stiffness matrix is ​​obtained through the derivation of the energy principle. Since the expression for the stiffness matrix is ​​an integral operation, the element stiffness matrix is ​​calculated by introducing numerical integration. The Hammer numerical integration results are as follows.

[0163]

[0164] S23, OpenMP parallel computing;

[0165] Since the element stiffness matrix calculation for each element is independent and occurs only within a single loop with few variables, OpenMP multithreading technology is used for parallel computation of the element stiffness matrix. The macroscopic model's elements are divided into blocks, and each thread calculates the element stiffness matrix for its assigned element, such as... Figure 4 As shown.

[0166] S3. The model is partitioned using a domain decomposition algorithm based on set values. This divides the solution mesh model into multiple computational sub-regions, ensuring that each partition carries all its parameters. The algorithm iterates through each partition based on the partition information to solve for the displacement field of the complete model. During the solution process, the open-source library STXXL is used to store data from memory on disk, reducing memory usage and enabling faster data retrieval from the disk.

[0167] S31, Model Domain Decomposition;

[0168] The number of regions N is determined by the size of the calculated model. The relationship between the number of regions and the partitioning method is typically as follows: when N is less than 8, a multi-level recursive equal partitioning method is used; when N is greater than or equal to 8, a multi-level K-way partitioning method is used. A partitioning diagram is shown below. Figure 5 As shown.

[0169] After partitioning, the nodes and elements of each partition are renumbered, and the node information, element information, and element information at the boundary of each partition are saved before and after renumbering. At this point, the original finite element model is divided into N regions Ω0 to Ω0. N-1 The partitioning results are as follows Figure 6 As shown.

[0170] S32, Calculation of partition stiffness matrix;

[0171] After completing the regional division of the model, the regional stiffness matrix of each region is calculated by integrating the regional stiffness matrix according to the regional situation. Based on the regional element information in S31, the element stiffness matrix of each region is extracted, and the nodal degrees of freedom are encoded according to the re-compiled nodal information of the regional region to construct the regional stiffness matrix.

[0172] S33. Iteratively solve the displacement field;

[0173] Firstly, based on the global Schwarz solution iteration in the domain decomposition algorithm, we have:

[0174]

[0175] Where ω is an arbitrary function mapping from the field Ω to R, u nThis is an approximate solution to the Poisson problem. This formula is used to simultaneously solve subproblems in all subdomains of the domain algorithm, and its convergence is proven using the maximum principle. Based on this, the residuals are introduced as follows;

[0176] r n =f+Δ(u) n )

[0177] This allows us to define the Restricted Additive Schwarz Algorithm (RAS) on a continuous level.

[0178] In finite element elastic systems, for quasi-static solutions, according to Hooke's law and the principle of minimum potential energy, we have:

[0179] KU = F

[0180] Where K is the global stiffness matrix, F is the force matrix, and U is the displacement vector result. Based on the domain decomposition algorithm, the above linear system can be transformed into the following block matrix form:

[0181]

[0182] Since the RAS algorithm is defined based on the continuous level, while the finite element algorithm is based on the discrete level, it is necessary to apply the RAS algorithm to the discrete level from the perspective of linear algebra. Introducing the block Jacobi algorithm, and letting D be a diagonal matrix of K, we have:

[0183] DU n+1 =DU n +(F-KU n )

[0184] Define the residual of the discrete level as r n =F-KU n By dividing the Jacobi algorithm into blocks, the block Jacobi algorithm can be defined as follows:

[0185]

[0186] Among them, K ii Here is the overall stiffness matrix for each partition. Let R be the displacement vector of the i-th partition at the n-th iteration step. Then, we introduce R... i As N to N i Restricted operators, As N to N i The transpose operator then has At this point, the block Jacobi algorithm can be written in the following form:

[0187]

[0188] Based on the above formula and combined with the RAS algorithm applicable to any number of partitions at the discrete level, we have:

[0189]

[0190] For the block Jacobi algorithm, when expanded into a matrix with 2 partitions, we have:

[0191]

[0192] in, This represents the displacement of the first partition in the nth iteration step. The non-overlapping parts, This represents the displacement of the first partition in the nth iteration step. The overlapping part. Represents the residual of the first partition in the nth iteration step. Non-overlapping parts, This represents the overlapping portion of the residuals of the first and second partitions in the nth iteration step.

[0193] The simplified linear algebraic form of the above equation is:

[0194]

[0195] in, This is the local correction for solving the discrete level in the nth step, and we have

[0196] Solve the linear equations for each partition, and then multiply the results by half for the overlapping parts and sum them. This form applies to all partition counts. For regions with more than two partitions, the above formula can be rewritten by simply including the partition information of the overlapping nodes in the constant term matrix. This completes the displacement field solution using the RAS algorithm based on domain decomposition.

[0197] S34, Memory data storage;

[0198] For the displacement field solution process described above, a model with tens of millions of elements requires storing a large amount of information in memory during actual calculations. This includes model information such as element numbers, node coordinates, and the B-matrix and material matrix of each element. This information consumes a significant amount of memory during computation. Using the open-source library STXXL allows the memory data to be stored on disk, reducing memory usage and enabling faster access to disk data.

[0199] This open-source library uses the currently popular Parallel Disk Model (PDM). In this model, file I / O is explicitly handled by the program. It divides M bytes of data in memory into D blocks and transfers them between D hard disks to amortize latency. Typically, external memory algorithms assume they read input from external storage (hard disk) and write output to it. STXXL, however, uses pipelining to provide a new interface for external memory algorithms, allowing them to provide output as a data stream directly to the algorithm using the output, instead of writing it to memory first. Furthermore, STXXL uses synchronous I / O calls to run asynchronously in separate threads; for each file, there is a read request queue, a write request queue, and a thread. The main thread sends requests to the file queues, and the threads associated with the files execute the requests in FIFO order. This method allows data in memory to be stored on hard disk, reducing memory usage.

[0200] S4. Based on the calculated displacement results, calculate the strain and stress at each Gaussian point of each element in the macroscopic model;

[0201] The specific implementation steps are as follows: calculate the stress and strain at the Gaussian points based on the three-dimensional tetrahedral element used. The three-dimensional tetrahedral element has four Gaussian points, and it is necessary to calculate the stress and strain in each direction at each Gaussian point separately. The calculation method is as follows.

[0202] ε=BU

[0203] Where ε is the strain matrix at a single Gaussian point, recording the strain and shear strain in three directions, B is the B matrix at that Gaussian point, and U is the displacement matrix of each node on the element, sorted according to the element connection order. The stress matrix is ​​then calculated after obtaining the strain matrix.

[0204] σ=Dε

[0205] Where σ is the stress matrix at a single Gaussian point, recording the stress and shear stress in three directions, and D is the material matrix of the element.

[0206] This completes the calculation of strain and stress at each Gaussian point of each element in the macroscopic model based on the calculated displacement results;

[0207] S5. Using the multi-process parallel technology (MPI) of the cluster, the units of the macro model are evenly distributed to each process, and each process calculates the displacement of the micro model at each Gaussian point of its own macro model unit.

[0208] The specific implementation steps are as follows: Since it is necessary to calculate the displacement of the microscopic model at each Gaussian point, and this involves more calculations than a single loop, MPI multi-process technology is used for parallel processing. First, the macroscopic elements are divided into blocks, and the stress, strain, and displacement of each Gaussian point on each block are distributed to each process using the MPI communication protocol, such as... Figure 7 As shown.

[0209] S6. Reconstruct the model at each Gaussian point of each unit in the macroscopic model under each process, use the microscopic model as the microstructure of the macroscopic model, and use the strain at the Gaussian point to calculate the displacement boundary conditions of the microstructure at that point.

[0210] The specific implementation steps are as follows: Reconstruct the model at each Gaussian point on each element of the macroscopic model, and use the strain at that Gaussian point to calculate the displacement boundary of the microscopic model. Since the finite element method first calculates the stress and strain at the Gaussian point, and then uses interpolation to calculate the stress and strain at other locations on the element, the microscopic model is treated as a point on the macroscopic model. The new displacement and stress / strain of the microscopic model at this point are calculated, and the result is then treated as the new stress at that Gaussian point, thus achieving the purpose of cross-scale calculation.

[0211] S7. In each process, CPU multithreading technology (OpenMP) is used to calculate the element stiffness matrix of each element in the microstructure in parallel, assemble the overall stiffness matrix and solve the displacement field, and calculate the strain and stress at each Gaussian point of each element in the micro model based on the displacement results.

[0212] S71, Calculation of boundary conditions for microscopic models;

[0213] After distributing the macroscopic model to each process according to the previous step, each process needs to perform a complete calculation of the microscopic model at each Gaussian point within that unit. At this point, each process contains all the model information of the microscopic model, as well as the strain at each Gaussian point. First, the boundary nodes and coordinates of the microscopic model are extracted, and the displacement boundary conditions of the microscopic model are obtained by multiplying the strain at the Gaussian point by the coordinates.

[0214] Δd=εd

[0215] Where Δd represents the displacement of each boundary node in the microscopic model, and d represents the coordinates of each boundary node in the microscopic model. This yields the boundary conditions of the microscopic model.

[0216] S72, Calculation of stiffness matrix of micro-model element;

[0217] The same finite element theory from S22 is used to calculate the shape functions, etc., in the stiffness matrix calculation of the microscopic model elements. Since the macroscopic model is treated as an isotropic material, while the microscopic model is treated as an anisotropic material, the calculation method for the material matrix D differs between the microscopic and macroscopic models. Therefore, the material matrix D for the microscopic model is calculated as follows:

[0218]

[0219] The 36 different parameters in the matrix are determined based on the different forms of anisotropic materials used. The stiffness matrix of each element is calculated using the same method as in S2.

[0220] S73, OpenMP parallel computing;

[0221] Accelerating the calculation of element stiffness matrix using OpenMP based on the same method in S23, and achieving simultaneous scheduling of processes and threads.

[0222] S74, Assembly of the overall stiffness matrix;

[0223] After completing the OpenMP calculation of the element stiffness matrix, the element stiffness matrix is ​​assembled into the global stiffness matrix according to the degree of freedom index of each node in each element to complete the assembly of the element stiffness matrix and obtain the global stiffness matrix K.

[0224] S75, Solving for displacement, stress, and strain;

[0225] Using the same method as in S41, the displacement field in the microscopic model and the stress and strain values ​​at each Gaussian point on the microscopic model are solved.

[0226] S8. In each process, the average stress of the entire micro-model is solved in parallel using CPU multi-threading technology (OpenMP) based on the stress at each Gaussian point on the micro-model. This average stress is then used as the average stress at a single Gaussian point of a single element in the macro-model. Finally, the internal forces of the macro-model element are calculated based on the average stress at the four Gaussian points of the single element in the macro-model.

[0227] S81, Calculation of mean stress in microscopic model;

[0228] Based on the stresses at all Gaussian points on a single microscopic model calculated in S75, the average stress on the microscopic model at that Gaussian point is calculated. Using the homogenization theory of simulating macroscopic mechanical properties through a composite material mesoscopic unit cell model, the average stress is calculated as follows:

[0229]

[0230] in, The average stress of the model is σ. ijLet be the stress at each Gaussian point.

[0231] S82 and OpenMP parallelism;

[0232] Since the mean stress calculations for each element in the microscopic model are independent, the S23 method is used to accelerate the calculation of the mean stress, and OpenMP is used to enable multithreading to accelerate the loop.

[0233] S83, Calculation of internal forces in microscopic model elements;

[0234] After obtaining the average stress of a single micro-model, return to the macro-model. That is, the stress at a single Gaussian point in the macro-model is equivalent to the average stress of the micro-model at that Gaussian point. Calculate the internal forces of the macro-model based on this average stress.

[0235]

[0236] Where B is the macroscopic element matrix B, S is the stress matrix, and J is the Jacobian matrix. The internal forces within a single element of the macroscopic model are then calculated.

[0237] S9. Each process completes the internal force calculation of all units within the process. Process 0 collects the calculation results of the internal forces of all units in the macroscopic model to obtain the internal force vector of the macroscopic model, and calculates the displacement of the macroscopic model based on the internal force vector.

[0238] S91, Collect internal force vectors;

[0239] After each process completes the internal force calculation for all units within its process, it sends the complete internal force vectors of each unit to process 0 via the MPI communication protocol. Once the communication is complete, i.e., process 0 has finished collecting the data, the complete internal force vector of the macroscopic model is obtained.

[0240] S92. Calculate the displacement of the macroscopic model;

[0241] Obtain the internal force vector F in Then, the new displacement field U is calculated using the following formula. new .

[0242] KU = F

[0243] S10. Iterate the displacement judgment residuals obtained by using the displacement calculated from the displacement boundary of the macroscopic model and the internal force vector, and obtain the displacement result of the macroscopic model after the judgment converges.

[0244] The specific implementation steps are as follows: The external force vector F obtained by enlarging the boundary conditions in S2 is... out The internal force vector F obtained from S9 in Adding them together gives F r , use F rSolving the stiffness matrix yields the new displacement U. n This is consistent with the previous step's convergence assessment. If convergence is calculated, F... r Calculated U n This is the final displacement. If r does not meet the convergence condition, then the displacement U obtained from the external force calculation in the previous step will be... n As a new step, displacement U n+1 , with U n+1 Set the boundary conditions for the microscopic model in S6, and repeat the above steps until convergence. The displacement results are as follows. Figure 8 As shown.

[0245] S11. By using GPU to calculate the stress and strain of the macroscopic model based on the displacement results, the post-processing of the calculation is quickly realized, and a cross-scale heterogeneous parallel algorithm for anisotropic material composition structures is achieved.

[0246] The specific implementation steps are as follows:

[0247] S1. Establish a discrete model of the composite material analysis object with arbitrary weaving method using the level set method;

[0248] The specific implementation steps are as follows: After calculating the model displacement in S10, the stress and strain of each element are calculated using the method in S41, but the data is copied to the GPU for calculation.

[0249] ε=BU

[0250] σ=Dε

[0251] Since this section involves only a large number of matrix operations, GPU parallelism, specifically CUDA, is used for post-processing in the finite element method. First, the displacement matrices, B matrices, and material matrices of all elements are copied from CPU memory to the GPU in one go. The number of blocks and threads are then set on the GPU. For strain calculation, each thread calculates one row of the B matrix and one column of the U matrix, achieving parallel calculation of strain results for all elements across all threads. For stress calculation, each thread calculates one row of the D matrix and one column of the strain matrix, achieving parallel calculation of stress results for all elements across all threads. This yields the stress and strain at each Gaussian point of all elements in the macroscopic model.

[0252] Since stress and strain are exact solutions obtained at the integration points, and the strain and stress at these integration points are located inside the element, linear interpolation is performed using nodal values ​​and shape functions to obtain the results at each nodal to improve the global strain and stress results. The calculation formula is as follows:

[0253]

[0254] When multiple elements share a single node, the values ​​at the node are interpolated from the integration points of several adjacent elements, resulting in multiple stress and strain values ​​at the node. Therefore, grinding is necessary. For strain and stress at a common node, if the outer surface is shared by four elements, such as... Figure 9 As shown:

[0255] The calculation method is as follows:

[0256]

[0257] This completes the implementation of heterogeneous parallelism between CPU and GPU based on a cluster, effectively utilizing computing resources to efficiently simulate anisotropic materials. The post-processing results are as follows: Figure 10 As shown.

[0258] This invention applies the DDM method and heterogeneous parallel computing strategy to the finite element calculation of anisotropic materials. Based on the advantages of different parallel methods, it makes full use of the numerical computing capabilities of computers to realize the calculation of large-scale models, and can realize the rapid calculation of the structural response of large-scale anisotropic materials.

[0259] The above description is not intended to limit the present invention in any way. Although the present invention has been disclosed through the above embodiments, it is not intended to limit the present invention. Any person skilled in the art can make changes or modifications to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures, characterized in that, Includes the following steps: S1. Establish macroscopic and microscopic geometric models of arbitrary scale, spatially discretize the macroscopic and microscopic geometric models to be analyzed, divide them into grids, and generate the required grid data. S2. Set boundary conditions, constraints, and material parameters for the macroscopic model, and use CPU multithreading technology OpenMP to calculate the element stiffness matrix of the macroscopic model in parallel. S3. Use the domain decomposition algorithm to partition the model according to the set values, divide the solution mesh model into multiple computational sub-regions, and make each partition carry all its own parameters. Based on the partition information, the displacement field of the complete model is solved by iterating through each partition. During the solution process, the open-source library STXXL is used to place the data in memory on the hard disk, thereby reducing memory usage and enabling fast reading of hard disk data. S4. Based on the calculated displacement results, calculate the strain and stress at each Gaussian point of each element in the macroscopic model; S5. Using the multi-process parallel technology MPI of the cluster, the units of the macro model are evenly distributed to each process, and each process calculates the displacement of the micro model at each Gaussian point of its own macro model unit. S6. Reconstruct the model at each Gaussian point of each unit in the macroscopic model under each process, use the microscopic model as the microstructure of the macroscopic model, and use the strain at the Gaussian point to calculate the displacement boundary conditions of the microstructure at that point. S7. In each process, CPU multi-threading technology OpenMP is used to calculate the element stiffness matrix of each element in the microstructure in parallel, assemble the overall stiffness matrix and solve the displacement field, and calculate the strain and stress at each Gaussian point of each element in the micro model based on the displacement results. S8. In each process, the average stress of the entire micro-model is solved in parallel using CPU multi-threading technology OpenMP based on the stress at each Gaussian point on the micro-model. This average stress is then used as the average stress at a single Gaussian point of a single element in the macro-model. Finally, the internal forces of the macro-model element are calculated based on the average stress at the four Gaussian points of the single element in the macro-model. S9. Each process completes the internal force calculation of all units within the process. Process 0 collects the calculation results of the internal forces of all units in the macroscopic model to obtain the internal force vector of the macroscopic model, and calculates the displacement of the macroscopic model based on the internal force vector. S10. Iterate the displacement judgment residuals obtained by using the displacement calculated from the displacement boundary of the macroscopic model and the internal force vector, and obtain the displacement result of the macroscopic model after the judgment converges. S11. By using GPU to calculate the stress and strain of the macroscopic model based on the displacement results, the post-processing of the calculation is quickly realized, and a cross-scale heterogeneous parallel algorithm for anisotropic material structures is achieved.

2. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S1 is specifically implemented as follows: establishing a three-dimensional tetrahedral basic mesh element with a reasonable density based on the complexity of macroscopic and microscopic scales; saving the number of model elements and the nodes contained in each element; saving the number of nodes and the three-dimensional coordinates of all nodes.

3. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S2 is specifically implemented as follows: S21. Set material parameters, boundary conditions and constraints; Material parameters are set for the macroscopic and microscopic models respectively. The macroscopic model is treated as an isotropic material, and Poisson's ratio and elastic modulus are set according to requirements. The microscopic model is an anisotropic material, and elastic modulus and Poisson's ratio need to be set in each direction. For the macroscopic model, boundary conditions and constraints also need to be set to determine the positions of displacement boundaries and fixed constraints. S22. Calculation of element stiffness matrix; Based on the set macroscopic model material parameters, the stiffness matrix of the macroscopic model element is calculated. First, the shape function of the three-dimensional tetrahedral element is calculated. For the tetrahedral element used, its volume coordinates are defined as follows: In the parametric coordinate system, the shape function is equal to the corresponding volume coordinate, and the four shape functions corresponding to the four nodes are: The Jacobian matrix is ​​obtained by taking the partial derivatives of the shape function with respect to each direction in the local coordinate system and then multiplying it by the nodal coordinates. The method for solving the Jacobian matrix is ​​as follows; The inverse of the Jacobi matrix is ​​multiplied by the partial derivative of the form function to obtain the B matrix; The D matrix form of an isotropic three-dimensional tetrahedral element is as follows; in, E and v are the elastic modulus and Poisson's ratio of the macroscopic model, respectively; The expression for the element stiffness matrix is ​​obtained through the derivation of the energy principle. Since the expression for the stiffness matrix is ​​an integral operation, the element stiffness matrix is ​​calculated by introducing numerical integration. The Hammer numerical integration results are as follows. S23, OpenMP parallel computing; Since the element stiffness matrix calculation for each element is independent and only occurs within a single loop with few variables involved, OpenMP multithreading technology is used for parallel calculation of the element stiffness matrix. The elements of the macroscopic model are divided into blocks, and each thread calculates the element stiffness matrix of the element assigned to it. This completes the setting of boundary conditions, constraints, and material parameters for the macroscopic model, and uses CPU multi-threading technology OpenMP to calculate the element stiffness matrix of the macroscopic model in parallel.

4. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S3 is specifically implemented as follows: S31, Model Domain Decomposition; The number of regions N is determined by the size of the calculated model. The relationship between the number of regions and the partitioning method is usually that when N is less than 8, a multi-level recursive equal partitioning method is used, and when N is greater than or equal to 8, a multi-level K-way partitioning method is used. After partitioning, the nodes and cells of each partition are renumbered, and the node information, cell information and cell information at the boundary of each partition are saved before and after renumbering. At this point, the original finite element model is divided into N regions Ω0~Ω N-1 S32, Calculation of partition stiffness matrix; After completing the regional division of the model, the partition stiffness matrix of each partition is integrated and calculated according to the partition situation; based on the partition element information in S31, the element stiffness matrix of each partition is extracted, and the partition stiffness matrix is ​​constructed by encoding the nodal degrees of freedom according to the partition re-compiled node information. S33. Iteratively solve the displacement field; Firstly, based on the global Schwarz solution iteration in the domain decomposition algorithm, we have: Where ω is an arbitrary function mapping from the field Ω to R, u n It is an approximate solution to the Poisson problem; this formula is used to solve subproblems in all subdomains of the domain algorithm simultaneously, and its convergence is proved using the maximum principle; based on this, the residual is introduced as follows; r n =f+Δ(u n ) This allows us to define the Restricted Additive Schwarz Algorithm (RAS) on a continuous level; In finite element elastic systems, for quasi-static solutions, according to Hooke's law and the principle of minimum potential energy, we have: KU = F Where K is the global stiffness matrix, F is the force matrix, and U is the displacement vector result; based on the domain decomposition algorithm, the above linear system is transformed into the following block matrix form: Since the RAS algorithm is defined based on the continuous level, while the finite element algorithm is based on the discrete level, it is necessary to apply the RAS algorithm to the discrete level from the perspective of linear algebra. By introducing the block Jacobi algorithm and letting D be the diagonal matrix of K, we have: YOU n+1 =YOU n +(F-KU n ) Define the residual of the discrete level as r n =F-KU n By dividing the Jacobi algorithm into blocks, the block Jacobi algorithm can be defined as follows: Among them, K ii Here is the overall stiffness matrix for each partition. Let R be the displacement vector of the i-th partition at the n-th iteration step; at this time, R is introduced. i As N to N i Restricted operators, As N to N i The transpose operator then has At this point, the block Jacobi algorithm can be written in the following form: Based on the above formula and combined with the RAS algorithm applicable to any number of partitions at the discrete level, we have: For the block Jacobi algorithm, when expanded into a matrix with 2 partitions, we have: in, This represents the displacement of the first partition in the nth iteration step. The non-overlapping parts, This represents the displacement of the first partition in the nth iteration step. The overlapping parts; Represents the residual of the first partition in the nth iteration step. Non-overlapping parts, This represents the overlapping portion of the residuals of the first and second partitions in the nth iteration step; The simplified linear algebraic form of the above equation is: in, This is the local correction for solving the discrete level in the nth step, and we have Solve the linear equations for each partition, and then multiply the results by one-half for the overlapping parts and add them together. This form is applicable to all partition numbers. For regions with more than 2 partitions, the above formula can be rewritten by simply including the partition information of the overlapping nodes in the constant term matrix. This completes the displacement field solution of the RAS algorithm based on domain decomposition. S34, Memory data storage; For the above displacement field solution process, for a model with tens of millions of elements, a large amount of information needs to be stored in memory in actual calculations, such as element number, node coordinates and other model information, as well as the B matrix, material matrix and other finite element model information for each element. This information will occupy a lot of memory in the calculation. Using the open source library STXXL, the memory data can be stored on the hard disk, thereby reducing memory usage and enabling fast reading of hard disk data. This open-source library uses the Parallel Disk Model (PDM). In this model, file I / O is explicitly handled by the program. The model divides M bytes of data in memory into D blocks and transfers them between D hard disks to amortize latency. Typically, external memory algorithms assume they read input from external storage (hard disk) and write output to it. STXXL, however, uses pipelining to provide a new interface for external memory algorithms, allowing them to provide output as a data stream directly to the algorithm using the output, instead of writing it to memory first. Furthermore, STXXL uses synchronous I / O calls to run asynchronously in separate threads. For each file, there is a read request queue, a write request queue, and a thread. The main thread sends requests to the file queues, and the threads associated with the files execute the requests in FIFO order. This approach allows data in memory to be stored on the hard disk, reducing memory usage. This completes the partitioning of the model using the domain decomposition algorithm based on set values, dividing the solution mesh model into multiple computational sub-regions while ensuring that each partition carries all its own parameters; iterating through each partition based on the partition information and solving the displacement field of the complete model; during the solution process, the open-source library STXXL is used to place the data in memory on the hard disk, thereby reducing memory usage and enabling fast reading of hard disk data.

5. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S4 is specifically implemented as follows: S41. Solve for the stress and strain at the Gauss point; The stress and strain at Gaussian points are calculated based on the three-dimensional tetrahedral element used. This three-dimensional tetrahedral element has four Gaussian points, and the stress and strain in each direction at each Gaussian point need to be calculated separately. The calculation method is as follows. ε=BU Where ε is the strain matrix at a single Gaussian point, recording the strain and shear strain in three directions, B is the B matrix at that Gaussian point, and U is the displacement matrix of each node on the element arranged in the order of element connection; after obtaining the strain matrix, the stress matrix is ​​calculated. σ=Dε Where σ is the stress matrix at a single Gaussian point, recording the stress and shear stress in three directions, and D matrix is ​​the material matrix of the element; This completes the calculation of strain and stress at each Gaussian point of each element in the macroscopic model based on the calculated displacement results.

6. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S5 is specifically implemented as follows: since it is necessary to calculate the displacement of the micro-model at each Gaussian point, and more computation is involved than a single loop, MPI multi-process technology is used for parallel processing; firstly, the macro-units are divided into blocks, and the stress, strain and displacement of each Gaussian point on each part of the unit are distributed to each process using the MPI communication protocol.

7. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S6 is specifically implemented as follows: reconstruct the model at each Gaussian point on each element of the macroscopic model, and use the strain at the Gaussian point to calculate the displacement boundary of the microscopic model; since the finite element method first calculates the stress and strain at the Gaussian point, and then uses the interpolation method to calculate the stress and strain at other locations of the element, the microscopic model is regarded as a point on the macroscopic model; calculate the new displacement and stress and strain of the microscopic model at this time, and then regard the result as the new stress at the Gaussian point, so as to achieve the purpose of cross-scale calculation.

8. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S7 is specifically implemented as follows: S71, Calculation of boundary conditions for microscopic models; After distributing the macroscopic model to each process according to the previous step, each process needs to perform a complete calculation of the microscopic model at each Gaussian point within that unit. At this time, each process contains all the model information of the microscopic model, as well as the strain at each Gaussian point. First, extract the boundary nodes and coordinates of the microscopic model, and use the strain at the Gaussian point multiplied by the coordinates to obtain the displacement boundary conditions of the microscopic model. Δd=εd Where Δd is the displacement of each boundary node on the micro-model, and d is the coordinate of each boundary node on the micro-model; at this point, the boundary conditions of the micro-model are obtained. S72, Calculation of stiffness matrix of micro-model element; The same finite element theory from S22 is used to calculate the shape functions, etc., in the stiffness matrix calculation of the microscopic model elements. Since the macroscopic model is treated as an isotropic material while the microscopic model is treated as an anisotropic material, the calculation method for the material matrix D differs between the microscopic and macroscopic models. Therefore, the material matrix D for the microscopic model is calculated as follows: The 36 different parameters in the matrix are determined based on the different forms of anisotropic materials used; and the stiffness matrix of each element is calculated in the same way as in S2. S73, OpenMP parallel computing; Accelerate the calculation of element stiffness matrix using OpenMP based on the same method in S23, and achieve simultaneous scheduling of processes and threads; S74, Assembly of the overall stiffness matrix; After completing the OpenMP calculation of the element stiffness matrix, the element stiffness matrix is ​​assembled into the global stiffness matrix according to the degree of freedom index of each node in each element to complete the assembly of the element stiffness matrix and obtain the global stiffness matrix K. S75, Solving for displacement, stress, and strain; Using the same method as in S41, the displacement field in the microscopic model and the stress and strain values ​​at each Gaussian point on the microscopic model are solved. This completes the solution of displacement, stress, and strain of the micro-model at a single Gaussian point. It achieves parallel computation of the element stiffness matrix of each element in the microstructure using CPU multi-threading technology OpenMP in each process, assembles the overall stiffness matrix, solves the displacement field, and calculates the strain and stress at each Gaussian point of each element in the micro-model based on the displacement results.

9. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S8 is specifically implemented as follows: S81, Calculation of mean stress in microscopic model; Based on the stresses at all Gaussian points on a single microscopic model calculated in S75, the average stress on the microscopic model at that Gaussian point is calculated. Using the homogenization theory of simulating macroscopic mechanical properties through a composite material mesoscopic unit cell model, the average stress is calculated as follows: in, The average stress of the model is σ. ij For the stress at each Gaussian point; S82 and OpenMP parallelism; Since the mean stress calculation of each element in the microscopic model is independent, the S23 method is used to accelerate the calculation of mean stress, and OpenMP is used to enable multithreading to accelerate the loop. S83, Calculation of internal forces in microscopic model elements; After obtaining the average stress of a single micro-model, return to the macro-model. That is, the stress at a single Gaussian point in the macro-model is equivalent to the average stress of the micro-model at that Gaussian point. Calculate the internal forces of the macro-model based on this average stress. Where B is the macroscopic element B matrix, S is the stress matrix, and J is the Jacobian matrix; at this point, the internal forces inside a single element in the macroscopic model are calculated. This completes the process of using CPU multi-threading technology (OpenMP) to solve the average stress of the entire micro-model in parallel under each process, based on the stress at each Gaussian point on the micro-model, as the average stress at a single Gaussian point of a single element in the macro-model. Then, based on the near-average stress at the four Gaussian points of a single element in the macro-model, the internal force of the macro-model element is calculated.

10. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S9 is specifically implemented as follows: S91, Collect internal force vectors; After each process completes the internal force calculation of all units within its process, the internal force vectors of each unit in each process are sent to process 0 via the MPI communication protocol; after the communication is completed, i.e., after process 0 has completed the collection, the complete internal force vector of the macroscopic model is obtained. S92. Calculate the displacement of the macroscopic model; Obtain the internal force vector F in Then, the new displacement field U is calculated using the following formula. new ; KU = F This completes the internal force calculation of all units within each process. Process 0 collects the calculation results of the internal forces of all units in the macroscopic model to obtain the internal force vector of the macroscopic model, and calculates the displacement of the macroscopic model based on the internal force vector.

11. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S10 is specifically implemented as follows: the external force vector F obtained by enlarging the boundary conditions in S2 is... out The internal force vector F obtained from S9 in Adding them together gives F r , use F r Solving the stiffness matrix yields the new displacement U. n This aligns with the previous step's convergence assessment. If convergence is calculated, F r Calculated U n This is the final displacement; if r does not meet the convergence condition, then the displacement U obtained from the previous step's external force calculation will be... n As a new step, displacement U n+1 , with U n+1 Set these as boundary conditions for the microscopic model in S6, and repeat the above steps until convergence.

12. The cluster-CPU-GPU heterogeneous parallel method for solving the response of large-scale anisotropic structures according to claim 1, wherein step S11 is specifically implemented as follows: after the model displacement is calculated in S10, the stress and strain of each element are calculated using the method in S41, but the data is copied to the GPU for calculation. ε=BU σ=Dε Since this section involves only a large number of matrix operations, GPU parallelism, specifically CUDA, is used for post-processing in the finite element method. First, the displacement matrices, B matrices, and material matrices of all elements are copied from CPU memory to the GPU in one go. The number of blocks and threads are set in the GPU. For strain calculation, each thread calculates one row of the B matrix and one column of the U matrix, achieving parallel calculation of strain results for all elements across all threads. For stress calculation, each thread calculates one row of the D matrix and one column of the strain matrix, achieving parallel calculation of stress results for all elements across all threads. This yields the stress and strain at each Gaussian point of all elements in the macroscopic model. Since stress and strain are exact solutions obtained at the integration points, and the strain and stress at the integration points are located inside the element, linear interpolation is performed using nodal values ​​and shape functions to obtain the results at each node in order to improve the global strain and stress results; the calculation formula is as follows: When multiple elements share a node, the values ​​at the node are interpolated from the integration points of several adjacent elements, resulting in multiple stress and strain values ​​at the node, which requires grinding. For example, the strain and stress at a common node are shared by four elements on the outer surface. The calculation method is as follows: This completes the implementation of CPU-GPU heterogeneous parallelism based on clusters, enabling efficient simulation of anisotropic materials by effectively utilizing computing resources.

Images